CN104268412A - 一种角道集射线层析偏移速度分析方法及装置 - Google Patents

Abstract

本发明涉及一种角道集射线层析偏移速度分析方法及装置，包括：利用速度模型对预处理后的野外采集地震数据进行偏移处理，获取偏移剖面和角度域成像道集；判断角度域成像道集中的同相轴是否水平；若同相轴不水平，更新速度模型；更新的方法包括：扫描偏移剖面的反射面，获取反射面倾角；选取偏移剖面的反射面上的点作为反射点，获取反射点的剩余曲率，以剩余曲率为约束条件，在角度域成像道集上自动拾取反射点的成像深度，获取不同反射角度的成像深度差；筛选出满足质量监控准则的反射点；利用筛选出的反射点对应地反射面倾角追踪出不同反射角度对应地射线路径，获取反射层析核函数，建立层析方程组；获取速度模型的更新量，更新初始速度模型。

Description

一种角道集射线层析偏移速度分析方法及装置

技术领域

本发明涉及地震勘探资料处理技术领域，特别涉及一种角道集射线层析偏移速度分析方法及装置。

背景技术

背景速度场的质量是偏移质量的前提条件。速度建模在工业界的应用已经从分析类方法发展到反演类方法。相对分析类方法，层析反演方法的假设更小，速度建模的精度更高。相对于其它的反演类速度建模方法，角道集中的射线层析偏移速度分析作为一种反演类的速度建模方法被广泛应用于实际生产中。

20世纪60年代，医学领域发展出了层析技术，1970年代引入地球物理学科中，并在1980年代发展完善了射线层析技术。从应用数据类型上来说，层析可以分作初至波层析、井间层析和反射层析，前两者应用透射数据反演速度模型，可称为透射层析，反射层析利用反射数据，在实施的策略上与透射层析有根本的不同。反射层析的主流方法有成像域反射层析与倾斜层析，成像域反射层析由Stork(1992)首次提出，并在1990年代得到广泛发展，也被称作层析偏移速度分析。层析偏移速度分析可以在偏移距域成像道集(ODCIG)或角度域成像道集(ADCIG)实现。由于偏移技术的限制，ODCIG层析偏移速度分析在1990年代就已经发展成熟，并集成于软件中，ADCIG层析偏移速度分析在2000年代才逐渐发展完善。

目前，现有的ADCIG层析偏移速度分析技术方案均存在如下问题：

1、在实际应用中，现有的ADCIG层析技术无法实现大规模稀疏矩阵存储；

2、在传统地ADCIG层析技术中，估计反射点真深度的误差比较大，严重影响层析偏移速度分析的精度。

需对现有技术方案进行改进，解决上述技术问题。

发明内容

为解决现有技术的问题，本发明提出一种角道集射线层析偏移速度分析方法及装置，解决了传统方法中需要估计反射点真深度的问题，并实现三维层析偏移速度分析中的大规模矩阵的稀疏存储。

为实现上述目的，本发明提供了一种角道集射线层析偏移速度分析方法，该方法包括：

利用速度模型对预处理后的地震数据进行偏移处理，获取偏移剖面和角度域成像道集；

判断角度域成像道集中的同相轴是否水平；如果同相轴水平，则表明所述速度模型满足精度；否则，更新所述速度模型；

其中，所述更新速度模型的方法包括：

扫描偏移剖面的反射面，获取反射面倾角；

选取偏移剖面的反射面上的点作为反射点，获取所述反射点处的剩余曲率，以所述剩余曲率为约束条件，在所述角度域成像道集上自动拾取所述反射点的不同反射角度对应的成像深度，利用所述成像深度获取不同反射角度的成像深度差；

筛选出满足质量监控准则的反射点，质量监控选用包括反射点的倾角、连续性、成像剖面能量、道集相似性在内的属性；；

利用筛选出的反射点对应地反射面倾角追踪出不同反射角度对应地射线路径，根据所述不同反射角度对应地射线路径获取反射层析核函数；

利用筛选出的反射点对应地成像深度差、所述反射层析核函数建立层析方程组；

求解所述层析方程组，获取所述速度模型的更新量，实现层析偏移速度分析。

优选地，所述层析方程组中的方程表达式为：

其中，s_m是偏移速度的慢度，是反射层倾角，Γ^θ是反射点的反射角θ对应地射线路径，Δs是待求解的慢度更新量，z^θ是角度θ对应地成像深度，θ₁和θ₂是同一反射点对应地两个不同的反射角。

优选地，所述方法还包括：

所述不同反射角度对应地射线路径构成层析矩阵，并对所述层析矩阵进行存储；所述层析矩阵的存储结构为变长度的二维矩阵形式；所述变长度的二维矩阵的第一列每行的数据是所述层析矩阵对应行非零值的个数；所述变长度的二维矩阵的第2n列每行的数据是所述层析矩阵对应行第n个非零值的列号，所述变长度的二维矩阵的第2n+1列每行的数据是所述层析矩阵对应行第n个非零值；所述变长度的二维矩阵每行的长度等于所述层析矩阵对应行非零值个数乘以2再加1，其中，n为自然数。

优选地，所述剩余曲率的表达式为：

$Z^{\mathrm{mig}}\left(\mathrm{\θ}\right)=\frac{z^{\mathrm{mig}}\left(0\right)}{\mathrm{\γ}}\sqrt{{\mathrm{\γ}}^2+({\mathrm{\γ}}^2-1){\tan }^2\mathrm{\θ}}$

其中，z^mig(θ)是角度θ对应地成像深度，γ是剩余曲率，θ是反射点对应地反射角。

为实现上述目的，本发明还提供了一种角道集射线层析偏移速度分析装置，该装置包括：

预处理单元，利用速度模型对预处理后的地震数据进行偏移处理，获取偏移剖面和角度域成像道集；

分析单元，用于判断角度域成像道集中的同相轴是否水平；如果同相轴水平，则表明所述速度模型满足精度；否则，更新所述速度模型；

其中，所述分析单元包括：

扫描模块，用于扫描偏移剖面的反射面，获取反射面倾角；

成像深度差模块，用于选取偏移剖面的反射面上的点作为反射点，获取所述反射点处的剩余曲率，以所述剩余曲率为约束条件，在所述角度域成像道集上自动拾取所述反射点的不同反射角度对应的成像深度，利用所述成像深度获取不同反射角度的成像深度差；

筛选模块，用于筛选出满足质量监控准则反射点，质量监控可以选用包括反射点的倾角、连续性、成像剖面能量、道集相似性在内的属性；

反射层析核函数模块，用于利用筛选出的反射点对应地反射面倾角追踪出不同反射角度对应地射线路径，根据所述不同反射角度对应地射线路径获取反射层析核函数；

层析方程组模块，用于利用筛选出的反射点对应地成像深度差、所述反射层析核函数建立层析方程组；

求解模块，用于求解所述层析方程组，获取所述速度模型的更新量，实现层析偏移速度分析。

优选地，所述层析方程组模块获取的层析方程组中的方程表达式为：

其中，s_m是偏移速度的慢度，是反射层倾角，Γ^θ是反射点的反射角θ对应地射线路径，Δs是待求解的慢度更新量，z^θ是角度θ对应地成像深度，θ₁和θ₂是同一反射点对应地两个不同的反射角。

优选地，所述装置还包括：存储单元；

所述存储单元将所述不同反射角度对应地射线路径构成的层析矩阵进行存储；所述层析矩阵的存储结构为变长度的二维矩阵形式；所述变长度的二维矩阵的第一列每行的数据是所述层析矩阵对应行非零值的个数；所述变长度的二维矩阵的第2n列每行的数据是所述层析矩阵对应行第n个非零值的列号，所述变长度的二维矩阵的第2n+1列每行的数据是所述层析矩阵对应行第n个非零值；所述变长度的二维矩阵每行的长度等于所述层析矩阵对应行非零值个数乘以2再加1，其中，n为自然数。

优选地，所述成像深度差模块获取的剩余曲率的表达式为：

$Z^{\mathrm{mig}}\left(\mathrm{\θ}\right)=\frac{z^{\mathrm{mig}}\left(0\right)}{\mathrm{\γ}}\sqrt{{\mathrm{\γ}}^2+({\mathrm{\γ}}^2-1){\tan }^2\mathrm{\θ}}$

其中，z^mig(θ)是角度θ对应地成像深度，γ是剩余曲率，θ是反射点对应地反射角。

上述技术方案具有如下有益效果：本技术方案先“拾取成像深度”和“扫描反射界面的反射倾角”，然后对反射点依据获取的信息阈值再进行筛选，对符合要求的反射点进行后续的层析。避免了现有技术中人工参与，减少人的主观能动性，提高了技术方案的精度。

在本技术方案进行层析偏移速度分析的过程中，计算成像深度与真深度之差，需要估计反射点的真深度，依现有技术对反射点的真深度的估计误差通常比较大，严重影响了层析偏移速度分析的精度。本技术方案把上述绝对深度差用相对深度差代替，消除了真深度估计误差对反演结果的影响。

本技术方案涉及的层析矩阵是规模巨大的稀疏矩阵，对层析矩阵的存储结构，比常规方法节省至少1/3内存，使得本技术方案比现有技术更具备实用性。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本技术方案的工作原理示意图；

图2为本发明提出的一种角道集射线层析偏移速度分析方法流程图；

图3为本发明提出的一种角道集射线层析偏移速度分析装置框图；

图4为本发明的装置中涉及层析矩阵存储结构示意图；

图5a为本实施例初始速度模型示意图；

图5b为本实施例初始速度模型更新后的SEG盐丘速度模型示意图；

图6a为本实施例利用图5a的初始模型进行叠前深度偏移处理得到的成像剖面示意图；

图6b为本实施例利用图5b更新后的速度模型重新进行叠前深度偏移处理得到的成像剖面示意图；

图7a为本实施例基于图6a对应地第300条Inline线cdp位置分别为200、300、400的角度域道集示意图；

图7b为本实施例基于图6b对应地第300条Inline线cdp位置分别为200、300、400的角度域道集示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，为本技术方案的工作原理示意图。本技术方案的工作原理是：从初始速度模型出发，对预处理后的野外采集地震数据进行叠前深度偏移处理，获取偏移剖面和角度域成像道集，通过判断角度域成像道集中的同相轴是否水平来衡量速度模型的质量，如果角度域成像道集中同相轴是水平的，则表明速度场已比较准确，否则利用同相轴的弯曲信息来更新速度场。角度域成像道集中的层析偏移速度分析以成像剖面和角度域成像道集为基础，在成像剖面中选取界面上的点作为反射点，计算反射点处的剩余曲率，并扫描得到反射点对应的反射面倾角。根据Snell定理，结合反射面倾角追踪出不同反射角对应的射线路径，从而计算出反射层析核函数。同时利用剩余曲率作为约束，在角度域成像道集上自动拾取出该反射点不同反射角度的成像深度，从而计算不同反射角度的成像深度之差。结合层析核函数和成像深度差建立层析线性方程组，求解层析方程组得到速度模型的更新量用于更新偏移速度，实现层析偏移速度分析的一次迭代。利用更新后的偏移速度进行叠前深度偏移得到新的成像剖面和角度域成像道集，如果对成像结果和角度域成像道集不满意，重复上述过程，进行下一次的层析偏移速度分析的迭代。

基于上述工作原理，本发明提出了一种角道集射线层析偏移速度分析方法，如图2所示，该方法包括：

步骤201)：利用速度模型对预处理后的地震数据进行偏移处理，获取偏移剖面和角度域成像道集；

步骤202)：判断角度域成像道集中的同相轴是否水平；如果同相轴水平，则表明所述速度模型满足精度；否则，转至步骤203；

步骤203)：扫描偏移剖面的反射面，获取反射面倾角；

步骤204)：选取偏移剖面的反射面上的点作为反射点，获取所述反射点处的剩余曲率，以所述剩余曲率为约束条件，在所述角度域成像道集上自动拾取所述反射点的不同反射角度对应的成像深度，利用所述成像深度获取不同反射角度的成像深度差；

步骤205)：筛选出满足质量监控准则的反射点，质量监控可以选用包括反射点的倾角、连续性、成像剖面能量、道集相似性在内的属性；

步骤206)：利用筛选出的反射点对应地反射面倾角追踪出不同反射角度对应地射线路径，根据所述不同反射角度对应地射线路径获取反射层析核函数；

步骤207)：利用筛选出的反射点对应地成像深度差、所述反射层析核函数建立层析方程组；

步骤208)：求解所述层析方程组，获取所述速度模型的更新量，实现层析偏移速度分析。

优选地，所述层析方程组中的方程表达式为：

其中，s_m是偏移速度的慢度，是反射层倾角，Γ^θ是反射点的反射角θ对应地射线路径，Δs是待求解的慢度更新量，z^θ是角度θ对应地成像深度，θ₁和θ₂是同一反射点对应地两个不同的反射角。

优选地，所述方法还包括：

所述不同反射角度对应地射线路径构成层析矩阵，并对所述层析矩阵进行存储；所述层析矩阵的存储结构为变长度的二维矩阵形式；所述变长度的二维矩阵的第一列每行的数据是所述层析矩阵对应行非零值的个数；所述变长度的二维矩阵的第2n列每行的数据是所述层析矩阵对应行第n个非零值的列号，所述变长度的二维矩阵的第2n+1列每行的数据是所述层析矩阵对应行第n个非零值；所述变长度的二维矩阵每行的长度等于所述层析矩阵对应行非零值个数乘以2再加1，其中，n为自然数。

优选地，所述剩余曲率的表达式为：

$Z^{\mathrm{mig}}\left(\mathrm{\θ}\right)=\frac{z^{\mathrm{mig}}\left(0\right)}{\mathrm{\γ}}\sqrt{{\mathrm{\γ}}^2+({\mathrm{\γ}}^2-1){\tan }^2\mathrm{\θ}}$

其中，z^mig(θ)是角度θ对应地成像深度，γ是剩余曲率，θ是反射点对应地反射角。

如图3所示，为本发明提出的一种角道集射线层析偏移速度分析装置框图。该装置包括：

预处理单元301，用于利用速度模型对预处理后的地震数据进行偏移处理，获取偏移剖面和角度域成像道集；

分析单元302，用于判断角度域成像道集中的同相轴是否水平；如果同相轴水平，则表明所述速度模型满足精度；否则，更新所述速度模型；

其中，所述分析单元302包括：

扫描模块3021，用于扫描偏移剖面的反射面，获取反射面倾角；

成像深度差模块3022，用于选取偏移剖面的反射面上的点作为反射点，获取所述反射点处的剩余曲率，以所述剩余曲率为约束条件，在所述角度域成像道集上自动拾取所述反射点的不同反射角度对应的成像深度，利用所述成像深度获取不同反射角度的成像深度差；

筛选模块3023，用于筛选出满足质量监控准则反射点，质量监控可以选用包括反射点的倾角、连续性、成像剖面能量、道集相似性在内的属性；

反射层析核函数模块3024，用于利用筛选出的反射点对应地反射面倾角追踪出不同反射角度对应地射线路径，根据所述不同反射角度对应地射线路径获取反射层析核函数；

层析方程组模块3025，用于利用筛选出的反射点对应地成像深度差、所述反射层析核函数建立层析方程组；

求解模块3026，用于求解所述层析方程组，获取所述速度模型的更新量，实现层析偏移速度分析。

优选地，所述层析方程组模块3025获取的层析方程组中的方程表达式为：

其中，s_m是偏移速度的慢度，是反射层倾角，Γ^θ是反射点的反射角θ对应地射线路径，Δs是待求解的慢度更新量，z^θ是角度θ对应地成像深度，θ₁和θ₂是同一反射点对应地两个不同的反射角。

优选地，所述装置还包括：存储单元；

所述存储单元将所述不同反射角度对应地射线路径构成的层析矩阵进行存储。如图4所示，为本发明的装置中涉及层析矩阵存储结构示意图。所述层析矩阵的存储结构为变长度的二维矩阵形式；所述变长度的二维矩阵的第一列每行的数据是所述层析矩阵对应行非零值的个数；所述变长度的二维矩阵的第2n列每行的数据是所述层析矩阵对应行第n个非零值的列号，所述变长度的二维矩阵的第2n+1列每行的数据是所述层析矩阵对应行第n个非零值；所述变长度的二维矩阵每行的长度等于所述层析矩阵对应行非零值个数乘以2再加1，其中，n为自然数。

优选地，所述成像深度差模块3022获取的剩余曲率的表达式为：

$Z^{\mathrm{mig}}\left(\mathrm{\θ}\right)=\frac{z^{\mathrm{mig}}\left(0\right)}{\mathrm{\γ}}\sqrt{{\mathrm{\γ}}^2+({\mathrm{\γ}}^2-1){\tan }^2\mathrm{\θ}}$

其中，z^mig(θ)是角度θ对应地成像深度，γ是剩余曲率，θ是反射点对应地反射角。实施例：

步骤1)：在测量工区或探区采集二维或三维地震数据；

步骤2)：对采集的地震数据进行常规预处理；

步骤3)：对预处理后地震数据进行时间域速度分析，获取较高精度的均方根速度模型；

步骤4)：将获取的均方根速度模型利用迪克斯公式转换为深度域层速度模型，作为角度域层析反演的初始模型；

初始模型也可以通过其他方法获得，作为模型测试，将真实速度模型乘以0.9获得初始模型。如图5a所示。

步骤5)：利用获得的深度域初始速度模型，对预处理后的地震数据进行叠前深度偏移处理，得到成像剖面和角度域成像道集；

当速度模型误差较大时，角度域成像道集上的同相轴会呈现出明显的上翘或下弯特征。如图6a所示，为本实施例利用图5a的初始模型进行叠前深度偏移处理得到的第300条Inline线的成像剖面示意图。如图7a所示，为本实施例基于图6a对应地第300条Inline线cdp位置分别为200、300、400的角度域道集示意图。道集同相轴呈现明显向上弯曲，表明初始速度偏低。

步骤6)：在成像剖面上扫描Inline和Crossline方向界面倾角及相似系数，存储相应结果，在本实施例中采用了多道互相关方法扫描界面倾角；

三维介质中反射点倾角的计算通过扫描InLine和CrossLine方向的二维倾角实现。把两个二维切向方向向量写成(l_x,l_y,l_z)和(c_x,c_y,c_z)，三维法向方向向量写成(t_x,t_y,t_z)，则有

$\left\{\begin{array}{c}{t}_x\&CenterDot;l_x+t_y\&CenterDot;l_y+t_z\&CenterDot;l_z=0\\ t_x\&CenterDot;c_x+t_y\&CenterDot;c_y+t_z\&CenterDot;c_z=0\\ t_x^2+t_y^2+t_z^2=1\\ t_z<0\end{array}\right.---\left(1\right)$

这样，就把求三维法向量的问题变成了求解InLine和CrossLine方向的两个2维切向方向向量的问题，使得问题的难度大大降低。下面是我们确定2维同相轴方向的方法。我们首先定义互相关函数。用表示地震数据体，为空间坐标，t为时间坐标。对于给定的数据点该点所对应的互相关函数(CCF)定义如下：

$\begin{array}{c}\mathrm{CCF}(\stackrel{\&RightArrow;}{g},{\stackrel{\&RightArrow;}{x}}_i,t)=\\ \underset{0<{\left|\right|{\stackrel{\&RightArrow;}{x}}_i-{\stackrel{\&RightArrow;}{x}}_k\left|\right|}_2<D^{t-\mathrm{\&Delta;}}}{\mathrm{\&Sigma;}}\underset{t-\mathrm{\&Delta;}}{\overset{t+\mathrm{\&Delta;}}{\&Integral;}}s({\stackrel{\&RightArrow;}{x}}_i,t^{\&prime;})s({\stackrel{\&RightArrow;}{x}}_k,t^{\&prime;}+({\stackrel{\&RightArrow;}{x}}_k-{\stackrel{\&RightArrow;}{x}}_i)\&CenterDot;\stackrel{\&RightArrow;}{g}){\mathrm{dt}}^{\&prime;}\end{array}---\left(2\right)$

上式中D为空间窗的半径，2Δ为时窗的长度，为同相轴的法线方向，同垂直的超平面方向即为同相轴的方向。可以看出，当是常数时，上式的互相关函数是一个关于的函数。然后沿着以向量为法线方向的超平面对数据过行重采样，求出重采样信号的方差COV：

$\begin{array}{c}\mathrm{CCF}(\stackrel{\&RightArrow;}{g},{\stackrel{\&RightArrow;}{x}}_i,t)=\underset{{\left|\right|{\stackrel{\&RightArrow;}{x}}_i-{\stackrel{\&RightArrow;}{x}}_k\left|\right|}_2<D}{\mathrm{\&Sigma;}}{\left|\right|s({\stackrel{\&RightArrow;}{x}}_k,t+({\stackrel{\&RightArrow;}{x}}_k-{\stackrel{\&RightArrow;}{x}}_i)\&CenterDot;\stackrel{\&RightArrow;}{g})-E\left|\right|}_2\\ E=\frac{1}{n}\underset{{\left|\right|{\stackrel{\&RightArrow;}{x}}_i-{\stackrel{\&RightArrow;}{x}}_k\left|\right|}_2<D}{\mathrm{\&Sigma;}}s({\stackrel{\&RightArrow;}{x}}_k,t+({\stackrel{\&RightArrow;}{x}}_k-{\stackrel{\&RightArrow;}{x}}_i)\&CenterDot;\stackrel{\&RightArrow;}{g})\end{array}---\left(3\right)$

其中，E为重采样信号的期望，n为空间窗内地震道的数量。利用上面求得的互相关函数和方差，我们可以将原始数据每一个采样点处的值由下面的公式确定下来：

$({\stackrel{\&RightArrow;}{x}}_i,t)\underset{g}{\underset{\&RightArrow;}{\min }}\frac{\mathrm{CCF}(\stackrel{\&RightArrow;}{g},{\stackrel{\&RightArrow;}{x}}_i,t)}{\mathrm{COV}(\stackrel{\&RightArrow;}{g},{\stackrel{\&RightArrow;}{x}}_i,t)+1}---\left(4\right)$

其中，表示在点处求得的最终的值，其值为一个向量。有了我们就可以得到空间任意点处最佳的同相轴法向向量：

$\begin{array}{c}\mathrm{grad}(\stackrel{\&RightArrow;}{x},t)=\arg \underset{g}{\underset{\&RightArrow;}{\min }}\underset{{\left|\right|{\stackrel{\&RightArrow;}{x}}_i-\stackrel{\&RightArrow;}{x}\left|\right|}_2<D}{\mathrm{\&Sigma;}}\left|\right|{\stackrel{\&RightArrow;}{g}-g\_\mathrm{get}({\stackrel{\&RightArrow;}{x}}_i,t+({\stackrel{\&RightArrow;}{x}}_i-\stackrel{\&RightArrow;}{x})\&CenterDot;\stackrel{\&RightArrow;}{g})\left|\right|}_2\end{array}---\left(5\right)$

在二维情况下，的值为同相轴的斜率，在高维线性空间中，其值为超平面的法向向量。有了InLine和CrossLine的两个方向的二维切向量，成像点处的三维法向量可利用(1)式计算。

步骤7)：利用希尔伯特变换计算成像剖面包络，记录剖面包络极大值点；

步骤8)：综合利用步骤6和步骤7中信息，设定阈值进行反射点筛选。筛选过程包括：

1)首先假定剖面上每个点的位置都是有效反射点；

2)给定反射界面倾角最小值和最大值，滤除超出范围的点；

3)给定剖面相似系数最小值，在2)结果中滤除相似系数低于给定反射界面倾角最小值的点；

4)在步骤7得到的记录剖面包络极大点信息中，设置阈值，删除数值小于阈值的包络极大值点；设置横向连续样点数最小值，删除横向连续样点数小于设置最小值的点，在本实施例中，横向连续定义为相邻道反射点误差不超过一个样点。

5)从3)和4)的筛选结果中选出共有的点作为最终有效反射点。

步骤9)：根据步骤8中确定的反射点，拾取各个反射点不同角度的成像深度。

上述在成像道集中“拾取成像深度”利用互相关计算。成像深度的自动拾取在剩余曲率的约束下进行，约束曲线的计算公式为：

$Z^{\mathrm{mig}}\left(\mathrm{\&theta;}\right)=\frac{z^{\mathrm{mig}}\left(0\right)}{\mathrm{\&gamma;}}\sqrt{{\mathrm{\&gamma;}}^2+({\mathrm{\&gamma;}}^2-1){\tan }^2\mathrm{\&theta;}}---\left(6\right)$

其中，z^mig(θ)是成像深度，γ是剩余曲率，θ是反射角。然后针对每个要拾取剩余深度差的成像点，在道集内每一道的理论成像深度(约束曲线)附近开一个窗(窗的长度作为一个参数在参数卡中给出)，用滑动时窗内的值与上一道做相关，公式如下

$r_{\mathrm{xy}}\left(\mathrm{\&tau;}\right)=\underset{n=1}{\overset{\mathrm{trace}\_\mathrm{len}}{\mathrm{\&Sigma;}}}{x}_n{y}_{n-\mathrm{\&tau;}}---\left(7\right)$

找出最大相关值对应窗口移动量τ，从而确定下一道拾取的成像深度z_m＝(n+τ)*d_z(d_z为深度上的采样间隔)。对一个共成像点道集中的所有道进行循环，则完成该成像点偏移深度的自动拾取。

步骤10)：建立层析方程，传统层析偏移速度分析中一个反射射线建立的方程可写成：

其中，s_m是偏移速度的慢度，是反射层倾角，Γ^θ是反射角θ对应的射线路径，Δs是待求解的慢度更新量，z^θ是角度θ对应的成像深度，z^true是真深度。由于真深度的估计误差比较大，本专利中利用相同反射点不同反射角对应的方程(8)建立相对深度差层析方程组，此时方程变为：

其中，θ₁和θ₂是两个不同的反射角。方程(9)所示的层析方程中不再需要反射点的真深度。

步骤11)：层析方程组系数存储。层析矩阵的行数为n_ray，即反射射线的个数，列数为模型参数点的个数n_model，一般三维介质层析中，n_ray和n_model的数量级为10⁷～10⁹，所以该矩阵必须压缩存储。因为层析矩阵是稀疏矩阵，即矩阵元素大部分为零，可以仅存储矩阵中的非零值实现压缩存储，常规压缩存储方法存储非零值的行号、列号及元素值。本技术方案中设计了一种新的稀疏矩阵压缩存储方式，该存储方式可认为是一个变长度的二维数组，所述变长度的二维矩阵的第一列每行的数据是所述层析矩阵对应行非零值的个数；所述变长度的二维矩阵的第2n列每行的数据是所述层析矩阵对应行第n个非零值的列号，所述变长度的二维矩阵的第2n+1列每行的数据是所述层析矩阵对应行第n个非零值；所述变长度的二维矩阵每行的长度等于所述层析矩阵对应行非零值个数乘以2再加1，其中，n为自然数。

具体用C程序语言的结构体定义实现。从矩阵的稀疏存储结构的定义中可看出，该方法只需存储矩阵非零元素的列号和值，比传统定义方法少了非零元素的行号，一般来说可节省1/3的存储量。

步骤12)：对步骤10建立的层析方程进行求解，可以采用SIRT、LSQR等多种方程组求解方法进行求解，得到慢度更新量。

步骤13)：对步骤12得到的慢度更新量进行平滑，转换为速度，对初始速度进行更新，图5b即为更新后的SEG盐丘速度模型。

步骤14)：利用更新后的速度模型重新进行叠前深度偏移，得到新的成像剖面，如图6b所示。同时，获取角度域成像道集，如图7b所示。从图7b可以看出，速度更新后，角度域成像道集同相轴水平，满足精度要求，停止迭代。

以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种角道集射线层析偏移速度分析方法，其特征在于，该方法包括：

利用速度模型对预处理后的地震数据进行偏移处理，获取偏移剖面和角度域成像道集；

判断角度域成像道集中的同相轴是否水平；如果同相轴水平，则表明所述速度模型满足精度；否则，更新所述速度模型；

其中，所述更新速度模型的方法包括：

扫描偏移剖面的反射面，获取反射面倾角；

选取偏移剖面的反射面上的点作为反射点，获取所述反射点处的剩余曲率，以所述剩余曲率为约束条件，在所述角度域成像道集上自动拾取所述反射点的不同反射角度对应的成像深度，利用所述成像深度获取不同反射角度的成像深度差；

筛选出满足质量监控准则的反射点，质量监控选用包括反射点的倾角、连续性、成像剖面能量、道集相似性在内的属性；

利用筛选出的反射点对应地反射面倾角追踪出不同反射角度对应地射线路径，根据所述不同反射角度对应地射线路径获取反射层析核函数；

利用筛选出的反射点对应地成像深度差、所述反射层析核函数建立层析方程组；

求解所述层析方程组，获取所述速度模型的更新量，实现层析偏移速度分析。

2.如权利要求1所述的方法，其特征在于，所述层析方程组中的方程表达式为：

其中，s_m是偏移速度的慢度，是反射层倾角，Γ^θ是反射点的反射角θ对应地射线路径，Δs是待求解的慢度更新量，z^θ是角度θ对应地成像深度，θ₁和θ₂是同一反射点对应地两个不同的反射角。

3.如权利要求1或2所述的方法，其特征在于，所述方法还包括：

所述不同反射角度对应地射线路径构成层析矩阵，并对所述层析矩阵进行存储；所述层析矩阵的存储结构为变长度的二维矩阵形式；所述变长度的二维矩阵的第一列每行的数据是所述层析矩阵对应行非零值的个数；所述变长度的二维矩阵的第2n列每行的数据是所述层析矩阵对应行第n个非零值的列号，所述变长度的二维矩阵的第2n+1列每行的数据是所述层析矩阵对应行第n个非零值；所述变长度的二维矩阵每行的长度等于所述层析矩阵对应行非零值个数乘以2再加1，其中，n为自然数。

4.如权利要求1或2所述的方法，其特征在于，所述剩余曲率的表达式为：

$z^{\mathrm{mig}}\left(\mathrm{\&theta;}\right)=\frac{z^{\mathrm{mig}}\left(0\right)}{\mathrm{\&gamma;}}\sqrt{{\mathrm{\&gamma;}}^2+({\mathrm{\&gamma;}}^2-1){\tan }^2\mathrm{\&theta;}}$

其中，z^mig(θ)是角度θ对应地成像深度，γ是剩余曲率，θ是反射点对应地反射角。

5.一种角道集射线层析偏移速度分析装置，其特征在于，该装置包括：

预处理单元，利用速度模型对预处理后的地震数据进行偏移处理，获取偏移剖面和角度域成像道集；

分析单元，用于判断角度域成像道集中的同相轴是否水平；如果同相轴水平，则表明所述速度模型满足精度；否则，更新所述速度模型；

其中，所述分析单元包括：

扫描模块，用于扫描偏移剖面的反射面，获取反射面倾角；

成像深度差模块，用于选取偏移剖面的反射面上的点作为反射点，获取所述反射点处的剩余曲率，以所述剩余曲率为约束条件，在所述角度域成像道集上自动拾取所述反射点的不同反射角度对应的成像深度，利用所述成像深度获取不同反射角度的成像深度差；

筛选模块，用于筛选出满足质量监控准则反射点，质量监控选用包括反射点的倾角、连续性、成像剖面能量、道集相似性在内的属性；

反射层析核函数模块，用于利用筛选出的反射点对应地反射面倾角追踪出不同反射角度对应地射线路径，根据所述不同反射角度对应地射线路径获取反射层析核函数；

层析方程组模块，用于利用筛选出的反射点对应地成像深度差、所述反射层析核函数建立层析方程组；

求解模块，用于求解所述层析方程组，获取所述速度模型的更新量，实现层析偏移速度分析。

6.如权利要求5所述的装置，其特征在于，所述层析方程组模块获取的层析方程组中的方程表达式为：

其中，s_m是偏移速度的慢度，是反射层倾角，Γ^θ是反射点的反射角θ对应地射线路径，Δs是待求解的慢度更新量，z^θ是角度θ对应地成像深度，θ₁和θ₂是同一反射点对应地两个不同的反射角。

7.如权利要求5或6所述的装置，其特征在于，所述装置还包括：存储单元；

所述存储单元将所述不同反射角度对应地射线路径构成的层析矩阵进行存储；所述层析矩阵的存储结构为变长度的二维矩阵形式；所述变长度的二维矩阵的第一列每行的数据是所述层析矩阵对应行非零值的个数；所述变长度的二维矩阵的第2n列每行的数据是所述层析矩阵对应行第n个非零值的列号，所述变长度的二维矩阵的第2n+1列每行的数据是所述层析矩阵对应行第n个非零值；所述变长度的二维矩阵每行的长度等于所述层析矩阵对应行非零值个数乘以2再加1，其中，n为自然数。

8.如权利要求5或6所述的装置，其特征在于，所述成像深度差模块获取的剩余曲率的表达式为：

$z^{\mathrm{mig}}\left(\mathrm{\&theta;}\right)=\frac{z^{\mathrm{mig}}\left(0\right)}{\mathrm{\&gamma;}}\sqrt{{\mathrm{\&gamma;}}^2+({\mathrm{\&gamma;}}^2-1){\tan }^2\mathrm{\&theta;}}$

其中，z^mig(θ)是角度θ对应地成像深度，γ是剩余曲率，θ是反射点对应地反射角。