地震层位自动追踪方法及装置与流程

本发明涉及地震勘探技术领域，具体涉及一种地震层位自动追踪方法及装置。

背景技术：

地震资料常常存在较多噪声，这些噪声大大干扰了自动层位追踪的结果。地震层位自动追踪方法主要根据地震波动力学和运动学特性进行对比追踪。常用的方法有基于波形特征的自动追踪和基于相关的自动追踪。基于波形特征的自动追踪只寻找搜索时窗内特征点的相似波形结构形态，不进行地震道之间的相关计算，得到的结果不太可靠。基于相关的自动追踪以种子点为中心，按照相关时窗范围，选取一段地震道，把这段地震数据与相邻道搜索时窗内的地震数据做相关分析，满足条件的点作为新的种子点。但是在地震信噪比低、同相轴不稳定的情况下，搜索不到满足条件的点，影响追踪结果。

为此，如何在信噪比较低的情况下取得较好的地震层位自动追踪结果成为当前需要解决的技术问题。

技术实现要素：

针对现有技术中的缺陷，本发明提供了一种地震层位自动追踪方法及装置，该方法能够较好的对噪音进行压制，在信噪比较低的情况下依然可以取得较好的自动追踪结果。

依据本发明一个方面，提供了一种地震层位自动追踪方法，其特征在于，包括以下步骤：

获取地震剖面数据，将所述地震剖面数据存储至二维数组a_i,j中,其中，i为line的个数，j为cdp的个数；

根据全变分(TV,Total Variation)算法对所述地震剖面数据进行去噪处理；

基于相似系数的相干体算法，将去噪处理后的地震剖面数据进行相干性计算以及倾角扫描；

在地震层位中选取起点，进行同相轴追踪。

可选地，所述根据全变分(TV,Total Variation)算法对所述地震剖面数据进行去噪处理包括：

根据以下迭代公式进行去噪处理：其中，n为迭代次数、i＝1,……,N、表示第n次迭代的结果，Δt为时间步长、边界条件满足：$a_{0,j}^n=a_{1,j}^n,a_{N,j}^n=a_{N-1,j}^n,$$a_{i,0}^n=a_{1,N}^n=a_{i,N-1}^n.$

可选地，所述基于相似系数的相干体算法，将去噪处理后的地震剖面数据进行相干性计算以及倾角扫描包括：

确定分析时窗，所述分析时窗为椭圆窗或矩形窗；

所述分析时窗内中心点的相干值为：

$C_s(t,p,q)=\frac{\underset{k=-K}{\overset{K}{\∑}}{\[\frac{1}{J}\underset{j=1}{\overset{J}{\∑}}{u}_j(t+\mathrm{k\Δt}-p{x}_j-q{y}_j)\]}^2}{\underset{k=-K}{\overset{K}{\∑}}\frac{1}{J}\underset{j=1}{\overset{J}{\∑}}{\[u_j(t+\mathrm{k\Δt}-p{x}_j-q{y}_j)\]}^2}$

其中，下标j表示落在分析时窗内的第j道，u_j表示分析时窗内第j道，x_j、y_j分别表示第j道与分析时窗内中心点t在x、y方向的距离，p、q分别表示分析时窗内中心点t在x、y方向的倾角，k为分析点；

进行倾角扫描，改变倾角p、q计算相干值，当所述相干值最大时获取相应的相干值及倾角p、q并存储。

可选地，所述在地震层位中选取起点，进行同相轴追踪包括：

根据层位设定追踪起始参考道，检索使该点处相干值最大对应的倾角p，再以下一道为参考道，循环计算至最后一道，获得同相轴结果。

依据本发明另一方面，提供了一种地震层位自动追踪装置，其特征在于，所述装置包括：

数据获取模块，适于获取地震剖面数据，将所述地震剖面数据存储至二维数组a_i,j中,其中，i为线数，j为cdp数；

去噪处理模块，适于根据全变分(TV,Total Variation)算法对所述地震剖面数据进行去噪处理；

相干计算模块，适于相似系数的相干体算法，将去噪处理后的地震剖面数据进行相干性计算以及倾角扫描；

同相轴追踪模块，适于在地震层位中选取起点，进行同相轴追踪。

可选地，所述去噪处理模块还适于：

根据以下迭代公式进行去噪处理：其中，n为迭代次数、i＝1,……,N、表示第n次迭代的结果，Δt为时间步长、边界条件满足：$a_{0,j}^n=a_{1,j}^n,a_{N,j}^n=a_{N-1,j}^n,$$a_{i,0}^n=a_{1,N}^n=a_{i,N-1}^n.$

可选地，所述相干计算模块还适于：

确定分析时窗，所述分析时窗为椭圆窗或矩形窗；

所述分析时窗内中心点的相干值为：

$C_s(t,p,q)=\frac{\underset{k=-K}{\overset{K}{\∑}}{\[\frac{1}{J}\underset{j=1}{\overset{J}{\∑}}{u}_j(t+\mathrm{k\Δt}-p{x}_j-q{y}_j)\]}^2}{\underset{k=-K}{\overset{K}{\∑}}\frac{1}{J}\underset{j=1}{\overset{J}{\∑}}{\[u_j(t+\mathrm{k\Δt}-p{x}_j-q{y}_j)\]}^2}$

其中，下标j表示落在分析时窗内的第j道，u_j表示分析时窗内第j道，x_j、y_j分别表示第j道与分析时窗内中心点t在x、y方向的距离，p、q分别表示分析时窗内中心点t在x、y方向的倾角，k为分析点；

进行倾角扫描，改变倾角p、q计算相干值，当所述相干值最大时获取相应的相干值及倾角p、q并存储。

可选地，所述同相轴追踪模块还适于：

根据层位设定追踪起始参考道，检索使该点处相干值最大对应的倾角p，再以下一道为参考道，循环计算至最后一道，获得同相轴结果。

由上述技术方案可知，本发明提供的种地震层位自动追踪方法及装置，通过TV方法对信噪比低的地震数据进行降噪，较好的实现了去噪和保护边缘的统一，同时能保护同相轴，基于相似系数的相干体算法，利用地震同相轴横向的连续性较好，在相干剖面上突出显示同相轴，从而方便进行同相轴追踪和倾角拾取。

附图说明

图1是根据本发明实施例提供的地震层位自动追踪方法的流程图；

图2是根据本发明的实施例提供的原始地震剖面数据图；

图3是根据本发明的实施例提供的直接对原始数据进行追踪获得的追踪结果图；

图4是根据本发明的实施例提供的使用全变分(TV)算法对所述地震剖面数据进行去噪处理后的地震数据图；

图5是根据本发明的实施例提供的TV方法结合相似性相干计算后追踪的结果图；

图6是根据本发明的实施例提供的地震层位自动追踪装置的结构图。

具体实施方式

下面结合附图，对发明的具体实施方式作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

图1示出了本发明实施例提供的地震层位自动追踪方法的流程图，如图1所示，本实施例的地震层位自动追踪方法包括：

101、将地震裸数据存储在二维数组中。

在本实施例中，获取地震剖面数据，将所述地震剖面数据存储至二维数组a_i,j中,其中，i为线数(line数)，j为cdp的个数。

102、根据全变分(TV,Total Variation)算法对所述地震剖面数据进行去噪处理；

在本实施例中，表示第n次迭代的结果，Δt为时间步长。用差商代替偏导数，可得各方向导数的差分如下：

${\left(a_x\right)}_{i,j}^n=\frac{1}{2}(a_{i+1,j}^n-a_{i-1,j}^n)$

${\left(a_x\right)}_{i,j}^n=\frac{1}{2}(a_{i+1,j}^n-a_{i-1,j}^n)$

${\left(a_{\mathrm{xx}}\right)}_{i,j}^n=a_{i+1,j}^n-2a_{i,j}^n+a_{i-1,j}^n$

${\left(a_{\mathrm{yy}}\right)}_{i,j}^n=a_{i+1,j}^n-2a_{i,j}^n+a_{i-1,j}^n$

${\left(a_{\mathrm{xy}}\right)}_{i,j}^n=\frac{1}{4}(a_{i+1,j+1}^n{-a}_{i-1,j+1}^n-a_{i+1,j-1}^n-a_{i-1,j-1}^n)$

散度的计算公式为：

$\▿\·\left(\frac{\▿a}{|\▿a|}\right)=\frac{a_y^2{a}_{\mathrm{xx}}-2a_x{a}_y{a}_{\mathrm{xy}}+a_x^2{a}_{\mathrm{yy}}}{{|\▿a|}^3}$

则求解方程的离散迭代格式为：

$a_{i,j}^{n+1}=a_{i,j}^n-\mathrm{\Δt\λ}\left(a_{i,j}^n{-a}_{i,j}^0\right)+\mathrm{\Δt}\[\▿\·\left(\frac{\▿a_{i,j}^n}{|\▿a_{i,j}^n|}\right)\]$其中，n为迭代次数；i＝1,……,N。边界条件满足$a_{0,j}^n=a_{1,j}^n,a_{N,j}^n=a_{N-1,j}^n,a_{i,0}^n=a_{1,N}^n=a_{i,N-1}^n.$

103、计算相干剖面，基于相似系数的相干体算法，将去噪处理后的地震剖面数据进行相干性计算以及倾角扫描；

在本实施例中，在进行相干性计算时，首先定义一个分析窗口(椭圆窗或者矩形窗)，在分析窗之内以分析点为中心包含J个地震道。定义相似系数σ(t,p,q)为分析窗内沿一定倾角p，q各道的平均能量与所有道能量的平均值的比值，我们把一个时间窗口内所有相似系数的平均值作为该窗口内中心点的相干值C_s(t,p,q)，即

$C_s(t,p,q)=\frac{\underset{k=-K}{\overset{K}{\∑}}{\[\frac{1}{J}\underset{j=1}{\overset{J}{\∑}}{u}_j(t+\mathrm{k\Δt}-p{x}_j-q{y}_j)\]}^2}{\underset{k=-K}{\overset{K}{\∑}}\frac{1}{J}\underset{j=1}{\overset{J}{\∑}}{\[u_j(t+\mathrm{k\Δt}-p{x}_j-q{y}_j)\]}^2}$

其中u_j表示分析窗内第j道，x_j，y_j分别表示该道距离中心道的距离，p，q分别表示x，y方向的倾角。基于相似系数的相干性计算利用多道的信息。相干时窗的大小取T/2～3T/2。为减小误差，在计算时给定一个阈值，小于此值的相干值将被赋0。计算相干性的过程中对倾角进行扫描，使用与界面倾角p、q相一致的倾角进行相干性计算时得到的相干性最大，可以获得界面倾角。

104、拾取层位起点，即在地震层位图像中选取起点，可以经由人工手动选取，但本发明并不限制于此。

105、层位追踪、倾角拾取，根据层位设定追踪起始参考道的深度h,检索使该点处相干性最大对应的,计算时延t并记录。以下一道为参考道，将时延t赋给h。循环计算直至最后一道。采用链式存储方式 存储倾角扫描和同相轴追踪的结果，其优点是存储量较小，同时可以根据每一条同相轴的编号检索到每一条同相轴。

106、显示自动追踪结果。

107、确定所述追踪结果是否满足用户要求，此步骤可以由用户来确定，但本发明并不限制于此。如果追踪结果不满足要求，则可以通过设置时窗和倾角扫描大小调整追踪结果。

108、将符合要求的追踪结果进行存储，以备后续使用。

下面以一个实际地震数据为例对本发明进行进一步的说明。

原始地震剖面数据如图2所示，如果直接对原始数据进行追踪则获得图3所示的追踪结果。

根据本发明的方法，使用全变分(TV,Total Variation)算法对所述地震剖面数据进行去噪处理后，获得如图4所示的地震数据图，在此基础上使用TV方法结合相似性相干计算后追踪结果图得到的追踪结果如图5所示。

与前述方法实施例相应，本发明实施例还公开了一种地震层位自动追踪装置，参照图6所示的结构图，具体可以包括：

数据获取模块，适于获取地震剖面数据，将所述地震剖面数据存储至二维数组a_i,j中,其中，i为线数(line的个数)，j为cdp数；

去噪处理模块，适于根据全变分(TV,Total Variation)算法对所述地震剖面数据进行去噪处理；

相干计算模块，适于相似系数的相干体算法，将去噪处理后的地震剖面数据进行相干性计算以及倾角扫描；

同相轴追踪模块，适于在地震层位中选取起点，进行同相轴追踪。

本发明通过全变分(TV)方法对信噪比低的地震数据进行降噪， 提高地震剖面分辨率，较好的实现了去噪和保护边缘的统一，同时能保护同相轴。基于相似系数的相干体算法，利用地震同相轴横向的连续性较好，在相干剖面上突出显示同相轴。基于相似系数的相干算法获取构造倾角，用倾角扫描法进行同相轴追踪，从而方便进行同相轴追踪和倾角拾取。本方法分辩率较高，且对噪音的压制效果较好，在信噪比较低的情况下依然可以取得对实际地震剖面资料进行自动层位追踪较好的结果。

本发明的说明书中，说明了大量具体细节。然而，能够理解，本发明的实施例可以在没有这些具体细节的情况下实践。在一些实例中，并未详细示出公知的方法、结构和技术，以便不模糊对本说明书的理解。

类似地，应当理解，为了精简本发明公开并帮助理解各个发明方面中的一个或多个，在上面对本发明的示例性实施例的描述中，本发明的各个特征有时被一起分组到单个实施例、图、或者对其的描述中。然而，并不应将该公开的方法解释呈反映如下意图：即所要求保护的本发明要求比在每个权利要求中所明确记载的特征更多的特征。更确切地说，如下面的权利要求书所反映的那样，发明方面在于少于前面公开的单个实施例的所有特征。因此，遵循具体实施方式的权利要求书由此明确地并入该具体实施方式，其中每个权利要求本身都作为本发明的单独实施例。

本领域技术人员可以理解，可以对实施例中的设备中的模块进行自适应性地改变并且把它们设置在于该实施例不同的一个或多个设备中。可以把实施例中的模块或单元或组件组合成一个模块或单元或组件，以及此外可以把它们分成多个子模块或子单元或子组件。除了这样的特征和/或过程或者单元中的至少一些是互相排斥之处，可以采用任何组合对本说明书(包括伴随的权利要求、摘要和附图)中公开的所有特征以及如此公开的任何方法或者设备的所有过程或单元 进行组合。除非另外明确陈述，本说明书(包括伴随的权利要求、摘要和附图)中公开的每个特征可以由提供相同、等同或相似目的的替代特征来代替。

此外，本领域的技术人员能够理解，尽管在此所述的一些实施例包括其它实施例中所包括的某些特征而不是其它特征，但是不同实施例的特征的组合意味着处于本发明的范围之内并且形成不同的实施例。例如，在下面的权利要求书中，所要求保护的实施例的任意之一都可以以任意的组合方式来使用。

本发明的各个部件实施例可以以硬件实现，或者以在一个或者多个处理器上运行的软件模块实现，或者以它们的组合实现。本领域的技术人员应当理解，可以在实践中使用微处理器或者数字信号处理器(DSP)来实现根据本发明实施例的一种浏览器终端的设备中的一些或者全部部件的一些或者全部功能。本发明还可以实现为用于执行这里所描述的方法的一部分或者全部的设备或者装置程序(例如，计算机程序和计算机程序产品)。这样的实现本发明的程序可以存储在计算机可读介质上，或者可以具有一个或者多个信号的形式。这样的信号可以从因特网网站上下载得到，或者在载体信号上提供，或者以任何其他形式提供。

应该注意的是上述实施例对本发明进行说明而不是对本发明进行限制，并且本领域技术人员在不脱离所附权利要求的范围的情况下可设计出替换实施例。在权利要求中，不应将位于括号之间的任何参考符号构造成对权利要求的限制。单词“包含”不排除存在未列在权利要求中的元件或步骤。位于元件之前的单词“一”或“一个”不排除存在多个这样的元件。本发明可以借助于包括有若干不同元件的硬件以及借助于适当编程的计算机来实现。在列举了若干装置的单元权利要求中，这些装置中的若干个可以是通过同一个硬件项来具体体现。单词第一、第二、以及第三等的使用不表示任何顺序。可将这些单词解释 为名称。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围，其均应涵盖在本发明的权利要求和说明书的范围当中。