CN102819040A - 基于中心扩散加倾角属性的三维地震层位自动追踪方法 - Google Patents

Abstract

本发明公开了一种基于中心扩散加倾角属性的三维地震层位自动追踪方法，包括如下步骤：调整种子点位置；数据初始化准备；利用种子点的倾角属性确定映射点的位置，在映射点的时窗范围内寻找振幅极大值点，在每个振幅极大值点周围取一段波形数据，和种子点的波形数据进行相关系数计算，然后找出相关系数最大的映射点，并计算该映射点和种子点的振幅比；如果振幅比小于等于设定的门限值，则将该映射点确定为追踪点。本发明的积极效果是：利用中心扩散算法，有利于解决种子点提供不充分和种子点位置分布不均匀情况下，沿剖面追踪导致的层位大空洞；利用倾角属性在追踪平面的角度投影，更加精确了映射点的位置，提高了追踪的精确性效果。

Description

基于中心扩散加倾角属性的三维地震层位自动追踪方法

技术领域

本发明属于地震勘探技术领域，特别涉及一种基于中心扩散加倾角属性的三维地震层位自动追踪方法。

背景技术

地震层位追踪时地震资料解释的基础工作，准确性对地震资料后续处理和释有着很大的影响，它在结构分析，特征识别和地址评估等方面都有着重要的意义，是油气监测的基础，而油气检测是地震勘探的主要目标之一。因此，地震层位追踪是很有必要和意义的。

对地震层位的追踪是地震解释工作的基础，目前有很多种地震层位的追踪方法，大体可以分为人工手动追踪和自动追踪。

早期的层位追踪主要采用人工追踪的方法，解释人员根据地震波运动学和动力学的特征进行人工对比追踪，虽然人工追踪具有拾取准确、受干扰因素影响小的优点，但是随着油气勘探的深入，地震勘探数据量的增大，传统的人工追踪方法在大数据量面前存在效率偏低，容易受不同分析人永远的偏差和错误影响的缺点。

后来出现半自动追踪则是由解释人员选择一个种子点，定义好追踪参数，由计算机在单个剖面上追踪符合条件的点作为层位点。通过种子点搜索相邻道的“次种子点”，依次类推。虽然半自动相对于人工自动追踪来说工作效率较高，但面对大数据量地震数据时，仍然相当费时费力。

作为油气检测的基础，尤其是对于大地震数据的处理，地震层位自动追踪的工作越来越重要。因此，现在主要广泛研究地震层位的自动追踪方法。目前自动追踪的方式有两类：特征追踪和相关追踪。特征追踪在制定的倾角窗口内寻找阳电的相似结构，但在道与道之间不进行任何相关操作。基于相关的自动追踪截取地震道上种子拾取周围的一部分，把它和相邻道进行相关比较，确定层位点。显然，相关追踪虽然计算量更大，但得到的结果可靠性更高。

为了提高层位追踪的工作效率，人们先后提出了很多地震层位自动追踪的方法，如人工神经网络方法，模式识别方法，互相关方法等。其中，利用线性插值和相关的概念进行追踪是采用最广泛的方法。

层位自动追踪时需要根据一定数量的种子点进行追踪，现在大多数的层位追踪方法都采取根据剖面进行追踪，在种子点数目和位置分布合理的情况下，这样自动追踪能将层位面完整追踪出来。但是，如果在种子点数目较少或者位置分布不合理的情况下，这样的按剖面进行追踪，追踪出来的层位面很有可能不完整。因此，找到一种能在种子点不充分的情况下，仍能完整追踪出层位的方法很有必要。中心扩散算法恰好能解决种子点不充分情况下层位追踪的不完整性。

由于在自动追踪中，相关追踪获取的结果更具有可靠性，因此大多数自动追踪算法都使用相关追踪层位。在相关追踪时，大多数都将种子点直接投影到相邻道上，然后在投影点附近获取波形来和种子点进行相关。在横向波形变化不大的剖面，这样做是合理的，但是如果在横向波形变化较大的剖面，这样直接平移种子点到相邻道，可直接影响后面相关波形的选择，影响后续追踪效果。因此，为了获取更加精确的效果，如何选择相邻道的投影点变得很有意义。因此，可以考虑利用地质中的倾角概念来帮助确定投影点的位置。

与本发明相关的现有技术包括：

因为作为油气检测的基础，地震层位的自动追踪的工作越来越重要，为了提高层位追踪的工作效果，人们提出了许多地震层位自动追踪的方法。比较典型的方法有线性插值追踪和相关模式追踪。

相关模式追踪方法，先由操作人员拾取一些控制线，然后在全区内实现插值，再以插值点作为参考点，以其上部或者下部一段小时窗内与种子点进行相关，得到的相关系数最大所对应的时移值就是要获取的层位点。

赵成喜等（“线性插值法地震剖面反射层位追踪”，西部探矿工程，2007年）使用了线性插值法自动追踪技术对地震层位进行追踪，在地震时间剖面，针对目的层位的种子点的邻区，选择小视窗和小道窗的数据体采用线性插值的数学方法进行追踪。该方法在特定的地质条件下取得了效果，其工作示意图如图1所示。但是，由于是在通过地质剖面的二维方式进行追踪，则有可能在没有足够多种子点的情况下，追踪出来的层位数据也局限在较小的区域，扩散性能不好。另外，在实际处理时，地震层位是否平行对于追踪效果有较大的影响。线性插值算法在地震层位不是水平时必须通过合理的编辑，才能得到较好的追踪曲线。

温庆庆等（“可视化地震资料解释系统的研究和开发”，西安科技大学硕士毕业论文，2008年4月）利用相关追踪算法进行地震层位追踪，先由人工选择一些控制线，然后再在全区进行插值，以插值点为参考点，在其上部或者下部一段时窗内与种子点作相关，去相关系数最大的点作为层位点。其追踪结果能清晰地看到层位构造走向。追踪模拟过程如图2所示。但是在处理过程中，对于相邻道点的选择具有随意性，如图2所示，这种追踪方式只是将种子点的位置投影到相邻道上，然后在该位置附近选取一个时窗进行波形相关计算，这在横向波形变化大时，容易产生计算误差，从而影响以后的追踪效果。

现有的基于相关追踪模型进行的三维层位追踪的算法研究一般都是根据一个剖面一个剖面进行追踪，并且依赖种子点的数量和其位置合理分布，然而有时种子点的提供并不充分，位置也不一定是按照一定规则来排布，如果采用沿剖面依种子点道进行追踪，可能导致大面积空洞的出现，这是不合理的。

其次，利用相关追踪相邻道时，映射点的位置简单的由种子点位置平移过来，在地震层位水平变化大的情况下，导致追踪的波形效果不佳。

发明内容

为了克服现有技术的上述缺点，本发明提供了一种基于中心扩散加倾角属性的三维地震层位自动追踪方法，基于种子点情况不充分的考虑，提出了采用中心扩散算法进行自动追踪的思路，很好地解决了沿剖面依种子点道追踪导致的大面积空洞的出现；同时，本发明还利用地质中的倾角属性的概念，引入倾角属性来确定相邻道映射点的位置，很好地解决了映射点位置的随意性问题。

本发明解决其技术问题所采用的技术方案是：一种基于中心扩散加倾角属性的三维地震层位自动追踪方法，包括如下步骤：

步骤一、调整种子点位置：

如果种子点位置不在局部振幅最大的位置则对种子点位置按如下方法进行调整：在种子点的时窗范围内寻找离种子点最近的极大值点，如果找到极大值点，则将种子点移动到极大值点的位置；如果未找到极大值点，则种子点位置不动；

步骤二、数据初始化准备：

首先，计算出整个三维地震数据体的梯度向量；其次，将各初始的种子点的位置放入队列queue中存放，并用一个数组array来标识每个种子点是否已经进行过相关追踪；

步骤三、中心扩散追踪：

1）计算种子点的倾角属性：

利用步骤二计算出的梯度向量，计算出数据体中每个种子点的局部梯度方差矩阵，再求出该矩阵的特征向量，将特征向量投影到追踪平面上，得到特征向量在追踪平面的投影角度；

2）确定映射点：根据第1）步确定的投影角度算出映射点的位置；

3）确定追踪点：

在第2）步确定的映射点的时窗范围内寻找振幅极大值点，在每个振幅极大值点周围取一段波形数据，和种子点的波形数据进行相关系数计算，然后找出相关系数最大的映射点，并计算该映射点和种子点的振幅比；判断振幅比是否大于设定的门限值，如果是，则放弃该映射点，并将该映射点所在的种子道标记为已追踪；如果否，则将该映射点确定为追踪点。

与现有技术相比，本发明的积极效果是：

（1）利用中心扩散算法，有利于解决种子点提供不充分和种子点位置分布不均匀情况下，沿剖面追踪导致的层位大空洞；

（2）利用倾角属性在追踪平面的角度投影，更加精确了映射点的位置，提高了追踪的精确性效果。

附图说明

本发明将通过例子并参照附图的方式说明，其中：

图1是线性插值法剖面自动追踪工作示意图；

图2是相关追踪模拟过程示意图；

图3是本发明的基于倾角属性的三维层位自动追踪算法；

图4是中心扩散追踪的示意图；

图5是地震数据体集合示意图；

图6是种子道集合和地震数据体集合的关系示意图；

图7是倾角示意图；

图8是time-xLine平面倾角属性投影图；

图9是time-inLine平面倾角属性投影图；

图10是确定映射点的示意图；

图11是波形相似比较示意图；

图12是中心扩散追踪的具体的流程图。

具体实施方式

一种基于中心扩散加倾角属性的三维地震层位自动追踪方法，如图3所示，包括如下步骤：

步骤一、调整种子点位置：严格来说，种子点的位置应该属于数据准备阶段，但是由于种子点位置对后续追踪有重要的影响，因此调整种子点的位置很有必要。

一般认为，振幅大的地方存在层位的可能性更大，因此，如果种子点位置不在局部振幅最大的位置则应该调整。调整种子点位置时，需要在种子点附近时窗范围内（如上下10个点）寻找离种子点最近的极大值点，如果找到极大值点，则将种子点移动到极大值点的位置；如果未找到极大值点，则种子点位置不动，但是也不从该种子点扩散追踪。

步骤二、数据初始化准备：在开始层位自动追踪前，要进行数据初始化准备。

首先，由于后续中心扩散中需要引入倾角属性，因此需要计算出整个三维地震数据体的三个方向（xLine,inline,time）的梯度向量。

其次，为了后续算法的使用，需要将各初始的种子点的位置放入队列queue中存放，该队列存放没有进行过中心扩散的种子点的位置。

另外，还需要一个数组array来标识每个点是否已经进行过相关追踪，初始化为0，表示没有进行过相关追踪。如果该点进行过相关追踪，无论是否找到符合要求的层位点，则标识为找过，记为1。

步骤三、中心扩散追踪：

为了解决沿剖面追踪因为种子点提供不足和位置不合理而导致的追踪空洞出现，本发明使用了中心扩散追踪算法。

下面首先介绍扩散的原理

中心扩散算法是指以种子点为中心，向周围四个方向进行追踪。种子点位于种子道数据里，由种子点向其前、后、左、右四个方向的相邻道进行追踪。通过扩散找到的点如果符合设定的门限条件，则作为下次扩散的种子点。

在三维数据体中，中心扩散追踪的示意图如图4所示，在地震数据体中，假设种子点的位置为（xLine,inLine），则扩散的方向为前、后、左、右四个方向，即（xLine,inLine-1）,(xLine,inline+1),(xLine-1,inLine),(xLine+1,inLine)。在实际的三维地震数据体数据中，中心扩散的过程涉及到数据集在不同集合中的转移交换。

三维地震数据体由很多道数据组成，组成每道数据的是一个个离散的采样点x_i。每道数据集合用A_i表示，则A_i＝{x₁,x₂,...,x_k,x_k+1，...,x_n-1，x_n}。对于地震数据体V，由多道数据组成，即V＝{A₁,A₂,...,A_k,A_k+1，...,A_n-1,A_n}。示意图如图5所示。标记的种子道的集合V₁初始化时由拥有种子点的道数据组成V₁＝{A₂,A₅,...,A_i,...,A_k}。

V₁中的种子道个数随机，并且其位置分布不需要具有规律性。这和一般的按剖面追踪时要求种子道数据的位置具有规律性有更大的实用范围，可以适应不同的应用场景。

初始时，V和V₁的关系如图6所示。种子道集合和三维地震数据体的关系为

黑色箭头表示V₁集合的扩充。随着中心扩散追踪的深入，由于通过中心扩散追踪出来的其他道将作为“次种子道”，即下次扩散的种子道数据，因此将被添加到V₁集合中，于是V₁的元素数逐渐增多。当然，以前已经进行过中心扩散追踪的种子道，将被标记，不再被用来追踪。

对于整个三维地震数据体的中心扩散追踪，追踪结束的标志是V₁=V，即V₁所包含的道数据集合的个数和集合V中的一样。

在具体的实现中，V₁集合用一个队列queue来表示，对已经追踪过的种子道在一个数组array中标记。相应方向追踪完毕，如果该道没有进行中心扩散，则加入队列queue，同时因为该道已经追踪过，在数组array中标记为1。

在中心扩散追踪过程中，对于中心扩散过程中的种子点，需要首先计算其倾角属性，通过得到的倾角属性，来确定向外扩散的相邻道的映射点的位置。利用倾角属性进行追踪，是首次在层位自动追踪过程中使用。

1）计算种子点的倾角属性：

本发明使用的中心扩散算法中，需要在扩散过程中，计算相应种子道的种子点的倾角属性。倾角属性的计算，是在层位追踪算法中的首次使用。利用倾角属性，可以在层位走势比较陡峭的时候，更加合理地选择映射点的位置，确保波形相关计算的效果。

根据Trygve Randen等提出的倾角属性的计算（“Three-DimensionalTexture Attributes for Seismic Data Analysis”，2000 SEG Annual Meeting,August 6-11,2000,Calgary,Alberta），在得到整个三维地震数据体中的梯度向量后，利用整体的梯度向量，在中心扩散的过程中，计算出数据体中每个种子点的局部梯度方差矩阵。每个点的倾角示意图如图7所示，θ为倾角，求出梯度方差矩阵后，再求出该矩阵的特征向量(λ₁,λ₂,λ₃)，其中λ₁表示种子点在xLine方向的坐标，λ₂表示种子点在inLine方向的坐标，λ₃表示种子点在time方向的坐标。得到特征向量后，在追踪时，根据追踪的方向，将特征向量投影到相应的平面上去。如果追踪的是time-xLine方向，则投影到time-xLine方向的平面，则其角度示意图如图8所示，如图所示，通过投影，得到特征向量在time-xLine平面的投影角度α为

$\mathrm{\α}=\arctan \left(\frac{{\mathrm{\λ}}_3}{{\mathrm{\λ}}_1}\right)$

如果追踪的是time-inLine方向，则投影到time-inLine方向的平面上，得到特征向量在time-inLine平面的投影角度β为

$\mathrm{\β}=\arctan \left(\frac{{\mathrm{\λ}}_3}{{\mathrm{\λ}}_2}\right)$

角度示意图如图9所示。

得到这两个平面的角度α,β后，从图4可知，追踪方向就主要是time-xLine和time-inLine这两个平面。

2）确定映射点：

相关波形追踪是利用类似于温庆庆等（“可视化地震资料解释系统的研究和开发”，西安科技大学硕士毕业论文，2008年4月）提出的相关追踪方法，来进行波形相关追踪。但本发明在映射点的选择上使用了倾角属性来帮助确定映射点位置，具体步骤如下所示：

映射点的选择需要利用已经计算出的倾角属性。首先，将算出的该种子点的倾角投影到要追踪的面上，如图6或图7所示。根据投影得到的角度，算出该映射点的位置。其示意图如图10所示。

因为倾角属性对于层位的总体走势有一个指向作用，利用如图6或图7的投影示意图可以计算出相邻道的映射点，而不再像温庆庆等提出的如图2所示那样直接平移种子点位置，更加精确地确定映射点。

如图10所示，h表示种子道和相邻道点的位置距离，在实际中，相邻道之间的位置距离为1，α为特征向量(λ₁,λ₂,λ₃)映射到追踪平面后，和该平面的夹角。d则是根据算出的夹角和h计算出来的映射点和种子点相距的位置点数,即d=h*tan(α)。将d算出来后，在地震数据中，这样，假设种子点的位置为a₀,则图10中所示的映射点的位置为a₀±d。

3）确定追踪点：

找到映射点后，在该映射点的一个设定的时窗范围（比如该点上下15个点）内，即A_i={x_k-15,x_k-14,x_k-13,...,x_k,x_k+1,...,x_k+15}内寻找该范围内的所有振幅极大值点，每个振幅极值点周围取一段波形数据，和种子点的波形数据进行相关系数计算，其示意图如图11所示。

相关系数r的计算公式如下式所指：

$r=\frac{E\left(\right[X-E\left(X\right)\left]\right[Y-E\left(Y\right)\left]\right)}{\sqrt{D\left(X\right)*D\left(Y\right)}}$

其中，X和Y分别是每个点的振幅值，E(X)和E(Y)分别是每段波形的平均振幅值，D(X)和D(Y)分别表示X和Y向量的方差。

最后找出相关系数最大的那个点。找到点后，本发明还需要比较该点和种子点的振幅比B，如下式所示，a表示目标道点的振幅，b表示种子道点的振幅。

$B=\frac{|a-b|}{|a+b|}$

由该公式可以看出，当两个点的振幅值相差越小时，B值越小。反之，则越大。判断振幅比B是否大于设定的门限值（比如0.3），如果是，则放弃该点，并将该点所在的种子道标记为已追踪；如果否，则将该点确定为追踪点。

中心扩散追踪的具体的流程图如图12所示。

Claims (4)

1.一种基于中心扩散加倾角属性的三维地震层位自动追踪方法，其特征在于：包括如下步骤：

步骤一、调整种子点位置：

如果种子点位置不在局部振幅最大的位置则对种子点位置按如下方法进行调整：在种子点的时窗范围内寻找离种子点最近的极大值点，如果找到极大值点，则将种子点移动到极大值点的位置；如果未找到极大值点，则种子点位置不动；

步骤二、数据初始化准备：

首先，计算出整个三维地震数据体的梯度向量；其次，将各初始的种子点的位置放入队列queue中存放，并用一个数组array来标识每个种子点是否已经进行过相关追踪；

步骤三、中心扩散追踪：

1）计算种子点的倾角属性：

利用步骤二计算出的梯度向量，计算出数据体中每个种子点的局部梯度方差矩阵，再求出该矩阵的特征向量，将特征向量投影到追踪平面上，得到特征向量在追踪平面的投影角度；

2）确定映射点：根据第1）步确定的投影角度算出映射点的位置；

3）确定追踪点：

在第2）步确定的映射点的时窗范围内寻找振幅极大值点，在每个振幅极大值点周围取一段波形数据，和种子点的波形数据进行相关系数计算，然后找出相关系数最大的映射点，并计算该映射点和种子点的振幅比；判断振幅比是否大于设定的门限值，如果是，则放弃该映射点，并将该映射点所在的种子道标记为已追踪；如果否，则将该映射点确定为追踪点。

2.根据权利要求1所述的基于中心扩散加倾角属性的三维地震层位自动追踪方法，其特征在于：所述用数组array来标识每个种子点是否已经进行过相关追踪的方法是：初始化为0，表示没有进行过相关追踪；如果进行过相关追踪，则标识为1，表示进行过相关追踪。

3.根据权利要求1所述的基于中心扩散加倾角属性的三维地震层位自动追踪方法，其特征在于：根据第1）步确定的投影角度计算映射点位置的方法为：

首先按如下公式计算映射点和种子点相距的位置点数d：

d=h*tan(α)

其中，h表示种子道和相邻道点的位置距离，α为特征向量在追踪平面的投影角度；则，映射点的位置为a₀±d，其中，a₀为种子点的位置。

4.根据权利要求1所述的基于中心扩散加倾角属性的三维地震层位自动追踪方法，其特征在于：所述相关系数的计算公式为：

$r=\frac{E\left(\right[X-E\left(X\right)\left]\right[Y-E\left(Y\right)\left]\right)}{\sqrt{D\left(X\right)*D\left(Y\right)}}$

其中，X和Y分别是每个映射点的振幅值，E(X)和E(Y)分别是每段波形的平均振幅值，D(X)和D(Y)分别表示X和Y向量的方差。

Images