CN103592681A - 一种基于信号分类的地震图像层位追踪方法 - Google Patents

Abstract

本发明公开了一种基于信号分类的地震图像层位追踪方法，在无需人工干预情况下，使用FSSCEM算法对得到的特征参数进行优化选择并聚类，并且将倾角纹理属性结合相关查找方法应用于聚类后的后续处理，实现了二维剖面上的全层位追踪，然后将半监督分类方法、倾角指向和基于相关性的层位追踪算法结合，实现了在三维地震图像中层位曲面精确追踪。本发明的方法无需大量人工干预；在降低聚类算法的计算量的同时保证了聚类的可靠性；增强了层位识别的连续性。

Description

一种基于信号分类的地震图像层位追踪方法

技术领域

本发明涉及地震图像处理领域，具体是一种基于信号分类的地震图像层位追踪方法。 

背景技术

随着经济的发展，社会对石油、天然气等资源的需求不断增大，然而要进行油气的开采，首先必须清楚地下的地质构造，因此地震解释在地质勘探过程中发挥了重要作用。地震层位追踪在地震解释工作中起到了非常重要的作用，通过地震层位追踪的方法可以找到可能含有油气的地层和油气的储存深度。 

由于地震图像中的同相轴是对勘探区域实际地质结构中的地质层位的映射，因此可以通过在地震图像中进行的同相轴追踪来实现地震层位追踪。地震波在向下传播过程中由于不同地层的组织构造不同、密度不同，会在地质界面上发生反射形成记录，最终识别和追踪该记录的过程就是层位追踪，即同相轴追踪。 

传统的层位追踪方法是通过解释人员手工完成的，处理时间长，结果不可验证。为了解决这些问题，近年来相关研究人员提出了很多关于层位自动追踪算法，使得在这方面的研究得到了迅速发展和广泛的关注。目前地震层位追踪所使用的方法主要有：相关方法、人工神经网络、遗传算法和基于有限混合高斯的全层位追踪算法等。 

（1）相关方法：通过提取待识别波形的实域段向量、自相关变换向量、最小相位分解向量和K-L变换向量等特征向量，并计算各个特征向量和确定层位的相同特征向量的相关系数，接着计算这些相关系数的均值求得一个综合相关系数作为相似性的判断标准，最后综合相关系数最大的实域段就是算法要追踪的层位。该方法法总体上可以比较清楚地分辨层位的结构，但是该方法存在的一个缺陷就是最终层位的不完整性。 

（2）人工神经网络层位追踪算法：主要将跨断层等不连续地质结构的层位追踪问题视为模式识别的问题，然后采用神经网络进行模式识别，但这种方法跨越复杂地质环境能力与训练样本状况密切相关，如果训练样本包含这种复杂地质状况，则能很好追踪，否则不能。而且训练样本获得需要大量人工干预标示。 

（3）遗传算法层位追踪算法：其具体思路是采用基于模型的方法来处理跨越断层的问题，然后将问题转化为约束最优化的问题，然后使用遗传算法来求解，这种方法存在两个问题：一是对复杂地质状况的逼近程度；二是求解精度和次优解问题。 

（4）基于有限混合高斯的全层位追踪算法，能适合复杂的地质环境，但需要人工干预，不能完全实现层位的自动提取和分离。 

发明内容

本发明所要解决的技术问题是，针对上述现有技术的不足，提供一种基于信号分类的地震图像层位追踪方法，实现二维和三维地震图像上的全层位自动精确追踪。 

为解决上述技术问题，本发明所采用的技术方案是：一种基于信号分类的地震图像层位追踪方法，该方法包括二维地震波形层位追踪方法和三维地震图像层位追踪方法： 

所述二维地震波形层位追踪方法为： 

1)寻找二维地震波形中的极值点：若二维地震波形中的某个点的值大于该点两边的点的值，且该点的值大于门限值，则该点为极大值点；若二维地震波形中的某个点的值小于该点两边的点的值，且该点的值小于门限值，则该点为极小值点； 

2)以上述各极值点为中心分别提取出一段二维地震波形，将切比雪夫零点作为拟合样本点计算二维地震波形拟合特征系数； 

3)使用SFS方法删选上述二维地震波形拟合特征系数，产生二维地震波形拟合特征系数的特征子集； 

4)采用CEM方法对二维地震波形拟合特征系数的特征子集聚类，计算散步分离性策略值，通过该策略值来评价上述特征子集和聚类效果； 

5)重复步骤4），直到遍历完所有特征子集，求得散步分离性策略值最大的特征子集，得到第i个极值点属于第j个聚类的后验概率P(c_j|X⁽ⁱ⁾)，若P(c_l|x⁽ⁱ⁾)>P(c_j|X⁽ⁱ⁾),

即第i个极值点属于第l个聚类的后验概率最大，则第i个极值点属于第l个层位，从而将步骤1）中的极值点分为若干个聚类，所述若干个聚类组成一个类集合，类集合中的每一个聚类代表一个层位，实现二维地震波形的层位追踪；其中： 

其中，

为第i个极值点的波形拟合特征系数，

表示第i个极值点的第k个波形拟合特征系数，N为拟合的阶数；c_j表示第j个聚类；μ_l为第l个层位的N×1的波形拟合特征系数的均值向量； n1为第l个层位上的极值点数目；μ_j为第j个层位的波形拟合特征系数的均值向量；Σ_j为第j个层位的N×N的协方差矩阵，∑_j=cov((X^(j))^T,X^(j))；参数α_l表示第l个层位的分布概率，取值范围为 0≤α_l≤1，

参数α_j表示第j个层位的分布概率，取值范围为0≤α_j≤1；

表示所有层位的混合高斯分布的概率密度函数，其中j=1,2,3...,K，

表示第j个层位的高斯分布的概率密度函数；K是层位数目； 

所述三维地震图像层位追踪方法为： 

1)截取三维地震图像上的一个二维剖面，并为该二维剖面提供至少两个种子点：若提供的种子点连线是垂直于xline_time平面的，则沿着inline方向提取二维剖面；若提供的种子层位线垂直于inline_time的平面，则沿着xline方向提取二维剖面； 

2)在上述种子点指定窗口大小内上下各搜索一个极值点作为拟修正位置，然后比较这两个拟修正位置与种子点的距离，选择距离近的拟修正位置，并将所述种子点移动到该拟修正位置；如果距离相等，且待修正的种子点为第一个种子点，则选取待修正的种子点指定窗口上深度值低的极值点作为拟修正位置，并将待修正的种子点移动到该拟修正位置，若待修正的种子点不是第一个种子点，则根据前一个种子点修正的方式进行种子点修正（即：待修正的第一个种子点移动到指定窗口上的极值点上；若第二个种子点的两个拟修正位置距离第二个种子点位置相等，此时按照第一个种子点修正，如果第一个种子点修正到指定窗口上的极值点上，则第2个种子点也修正到指定窗口上的极值点上；反之，如果第一个种子点修正到指定窗口下的极值点上，则第2个种子点也修正到指定窗口下的极值点上；依次类推）：如果前一二维剖面上修正后的种子点深度值更大，则将现有二维剖面的种子点往深度值更大的拟修正位置修正，如果前一二维剖面修正后的种子点深度值更小，则修正方向相反；若在指定窗口内不能找到极值点可以作为种子点修正后的位置，则认为层位在现有二维剖面数据上是没有点的，层位在现有二维剖面上被断开；所述指定窗口大小为种子点上下各15个点； 

3)在二维剖面中所有相邻两个修正后的种子点上下各偏移15个点所围成的平行四边形组合而成的区域内寻找极值点：若二维剖面中的某个点的值大于该点所在二维剖面两边的点的值，且该点的值大于门限值，则该点为极大值点；若二维剖面中的某个点的值小于该点所在二维剖面两边的点的值，且该点的值小于门限值，则该点为极小值点； 

4)以步骤2）修正后的种子点或步骤3）得到的极值点为中心提取出一段二维剖面波形，将切比雪夫零点作为拟合样本点计算二维剖面波形拟合特征系数； 

5)采用CEM方法对二维剖面波形拟合特征系数特征子集聚类，计算散步分离性策略值，通过该策略值来评价上述二维剖面波形拟合特征系数特征子集和聚类效果； 

6)重复步骤5），直到遍历完所有特征子集，求得散步分离性策略值最大的特征子集，得到第m个极值点或种子点属于第n个聚类的后验概率Q(c_m|Y⁽ⁿ⁾)，若Q(c_l|Y^(m))>Q(c_n|Y^(m)),即第m个极值点或种子点属于第l'个聚类的后验概率最大，则第m个极值点或种子点属于第l'个层位，从而将步骤3）中的极值点和步骤2）中修正后的种子点分为若干个聚类；其中： 

其中，

为第m个极值点或种子点的波形拟合特征系数，M为拟合的阶数，

表示第m个极值点或种子点的第k个波形拟合特征系数，M是拟合阶数；c_n、c_l'分别表示第n个、第l'个聚类，μ_l'为第l'个层位的M×1的波形拟合特征系数的均值向量，

k为第l个层位上的极值点或种子点数目；μ_n为第n个层位的波形拟合特征系数的均值向量；Σ_n为第n个层位的M×M的协方差矩阵，Σ_n=cov((Y⁽ⁿ⁾)^T,Y⁽ⁿ⁾)；参数α_l'表示第l'个层位的分布概率，取值范围为

参数α_n表示第n个层位的分布概率，取值范围为0≤α_n≤1；表示所有层位的混合高斯分布的概率密度函数，其中

n=1,2,3...,T，

表示第n个层位的高斯分布的概率密度函数；T是层位数目； 

7)选择修正后种子点最多的聚类作为层位，删去其余的聚类，完成了一个二维剖面上的层位追踪； 

8)重复三维地震图像层位追踪方法的步骤1）～步骤7），直到追踪完三维地震图像上所有的二维剖面。 

所述二维地震波形层位追踪方法还包括以下步骤： 

6)计算二维地震波形的梯度向量g： 

$g=\▿u(x,y)=\left[\begin{array}{c}\frac{\∂u(x,y)}{\∂x}\\ \frac{\∂u(x,y)}{\∂y}\end{array}\right]=\left[\begin{array}{c}{g}_x\\ g_y\end{array}\right],$

其中，g_x、g_y、分别表示x、y方向的梯度向量，u(x,y)表示二维地震图像中坐标为(x,y)点所对应的像素值； 

7)利用梯度向量g计算张量矩阵T，T=g·g^T； 

8)计算上述张量矩阵T的特征值λ_i，求得二维地震波形的倾角dip；其中i=1,2； 

$\mathrm{dip}=\arctan \left(\frac{v\left(y\right)}{\left|v\left(x\right)\right|}\right),$

其中，v(x)、v(y)为张量矩阵T的特征向量； 

9)利用所述倾角dip进行二维地震波形层位补齐，消除层位间断； 

10)消除每一个聚类中不同层位上的极值点； 

11)对层位进行片段融合。 

所述三维地震图像层位追踪方法还包括以下步骤： 

9)对于每一个二维剖面，选择相邻的两个种子点作为起点和终点，分别从前往后和从后往前两个方向通过相关查找的方法进行层位补齐； 

10)在提取出修正后的种子点所在层位在所有二维剖面上的位置后，将所得的二维剖面组合形成三维地震图像中的三维层位曲面。 

二维地震波形层位补齐的具体步骤如下： 

1)设定r=1； 

2)判断r是否小于等于类集合中的聚类的总数，若是，则进入3）；否则，进入8）； 

3)判断第r类对应的层位中是否存在缺口，若是，进入4）；否则，r的值加1，返回2）； 

4)查找层位中缺口处的层位点； 

5)根据所述倾角dip的指向查找与步骤4）中层位相邻的二维地震波形上的极值点； 

6)以缺口处层位点和上述步骤5）中各极值点为中心，上下各取15个点提取一段地震波形片段，分别计算以各极值点为中心的地震波形片段与以所述缺口处层位点 为中心的地震波形片段的相关系数，将该相关系数作为各极值点与所述缺口处层位点的相关系数；所述相关系数取值范围为0～1； 

7)判断是否存在大于门限值的相关系数，若是，则选择距离所述缺口处层位点最近的极值点添加到第r类对应的层位中，并将距离所述缺口处层位点最近的极值点作为新的缺口处的层位点，返回4）；否则，将未能补齐的层位划分出来单独作为一个层位，并将聚类的总数加1，返回3）； 

8)结束。 

消除每一个聚类中不同层位上的极值点的具体步骤为： 

1)设定s=1； 

2)判断s是否小于等于类集合中的聚类的总数，若是，则进入3）；否则，进入7）； 

3)判断第s类对应的层位是否有重叠，若有，则进入4）；否则，s的值加1，返回2）； 

4)选取第s类对应的层位上第一个二维地震波形上的点S作为原点； 

5)根据所述原点搜索所述起始道的相邻道上半径为W的窗口内的所有极值点集合D；所述半径W取值范围为5～30；本发明中半径W取15。 

6)判断D中是否含有与S同类的极值点S′，若是，则将S′作为原点，返回5）；否则，将第s类中的剩余点作为新的类，并将该新的类作为类集合的最后一个类，将聚类的总数加1，返回3）； 

7)结束。 

对层位进行片段融合的具体步骤如下： 

1)设定h=1； 

2)判断h是否小于等于类集合中的聚类的总数，若是，则进入3）；否则，进入6）； 

3)判断第h类最后一个点的搜索窗口内是否存在其他类，若是，则进入4）；否则，h的值加1，返回2）；所述搜索窗口大小为第h类最后一个点上下各15个点；； 

4)合并第h类和第t类中的极值点，作为新的类r； 

5)在类集合中删除第t类，聚类的总数减1，返回3）； 

6)结束。 

对每一个二维剖面进行层位补齐的具体步骤如下： 

1)设定p=1； 

2)判断p是否小于种子点数目，若是，则进入3）；否则，进入9）； 

3)以第p个种子点和第p+1个种子点为两端点； 

4)由两个端点同时向中间查找是否存在缺口；若存在缺口，则进入5)；否则，将p的值加1，返回2）； 

5)将缺口处层位点平移到与所述缺口处层位点所在的二维剖面相邻的二维剖面上，在相邻的二维剖面上，查找平移后的点上下各15个点范围内的极值点。 

6)以缺口处层位点和上述步骤5）中各极值点为中心，上下各取15个点提取一段地震波形片段，分别计算以各极值点为中心的地震波形片段与以所述缺口处层位点为中心的地震波形片段的相关系数，将该相关系数作为各极值点与所述缺口处层位点的相关系数；所述相关系数取值范围为0～1； 

7)判断是否存在大于门限值的相关系数，若是，选择与缺口处层位点相关系数最大的极值点，将该点作为新的缺口处层位点和新的端点，返回4）；否则，将p的值加1，返回2）； 

8)结束。 

与现有技术相比，本发明所具有的有益效果为：本发明基于信号分类的方法，实现了二维和三维地震图像上的全层位自动追踪。对于二维地震图像：使用FSSCEM算法对得到的特征参数进行优化选择并聚类，在降低聚类算法的计算量的同时保证了聚类的可靠性，无需人工干预就可以自动地将指定范围内的层位从二维地震图像中提取出来；通过倾角指向和相关查找的方法实现的层位补齐、层位去重叠和层位片段融合等后续处理方法，增强了层位识别的连续性；对于三维地震图像：在二维地震剖面图像上追踪层位，将三维地震图像层位追踪问题转化为了在二维地震剖面图像上的层位追踪问题；使用半监督聚类方法（需要种子点的方法就是半监督聚类），只需少量种子点就可以将层位提取出来，无需大量人工干预；通过相关查找的方法实现层位补齐，增强了层位识别的连续性。 

附图说明

图1为地震图像中的层位； 

图2为本发明二维地震波形全层位追踪流程图； 

图3为本发明地震道波形极值查找示意图； 

图4为本发明根据倾角的相关查找示意图； 

图5为本发明层位重叠去除流程图； 

图6为本发明层位片段融合流程图； 

图7为对二维剖面进行层位补齐的流程图； 

图8为四川某工区地震图像人工追踪效果图； 

图9为四川某工区地震图像本发明方法追踪流程图； 

图10为本发明三维地震图像层位追踪流程图； 

图11为本发明波形特征提取与半监督聚类流程图； 

图12为三维地震图像人工层位面追踪效果图； 

图13为本发明三维地震图像追踪方法追踪的三维地震图像效果图。 

具体实施方式

以下结合附图说明本发明的具体实施方式。 

1、本发明的二维地震图像中的全层位追踪方法 

地震波向下传播过程中，遇到地下界面发生反射，反射界面就是地层的分界线，即层位，如图1所示，在地震图像中对该分界线的提取就是层位追踪。本发明首先使用多项式拟合的方法重建地震波形，其多项式拟合系数作为波形的特征参数；其次，使用FSSCEM（Feature Subset Selection Component-wise EM，特征子集模型自动选择最大期望，EM（expectation-maximization）为最大期望）算法对得到的特征参数进行优化选择并聚类，以获得最优化特征参数集和相应的最优化聚类结果；最后，通过倾角指向和相关查找的方法实现的层位补齐、层位去重叠和层位片段融合等后续处理方法，实现在二维剖面的指定区域中所有层位的提取。总流程图如图2所示，具体方案如下： 

第一步、地震波形的极值查找：由于我们得到的地震图像数据都是离散的点，因此在这里我们通过比较离散值进行判断，若地震波形中某个点的值大于其两边的点的值，则这个点为极大值点，若小于两边的点的值则为极小值。 

由于在层位追踪的方法中需要遵循同相性准则，即同一层位的反射波在相邻地震记录上出现的时间是相似的，同相轴应该是平滑的，且延伸较长。因此在极值点查找中往往需要加入一个限制条件，即查找得到的极大值的振幅值要大于某个门限值，极小值的振幅值小于某个门限值，在本发明方法中我们选择0作为门限值。以一道地震数据为例最终得到的极值查找结果，如图3所示。对于查找得到的极值点我们将其存入到一个与地震数据同样大小的矩阵中，记录各道数据中的极值点位置。 

由于地震层位通常在极值点处，通过这一步，得到了地震图像中的所有层位点。 

第二步、波形特征提取：找到地震极值点后，以各极值点为中心分别提取出一段波形，本发明方法中选取波形长度为35，即极值点左右各取17个点。由于等间隔拟合容易引起龙格现象，本发明将切比雪夫零点（(i=0,1,…,n-1)）作为拟合样本点计算拟合特征系数，为之后的波形聚类作准备。在整个方法里面，切比雪夫拟合的阶数是影响很大的因素，如果阶数过小，得不到层位完整的波形，如果阶数过大，则包含太大的波形，本 发明方法选取的拟合阶数为9。 

通过这一步，将极值点附近的一段波形用少数几个拟合特征表示，能减少后面聚类的计算量并且增加聚类的可靠性。 

第三步、波形特征的无监督聚类：本发明中使用FSSCEM（Feature Subset Selection Compo-nent-wise EM）算法对波形拟合特征系数聚类，将地震图像中的极值点分为若干个聚类。首先，使用SFS（sequential forward search，顺序向前搜索）算法删选波形拟合特征系数，产生特征子集，降低聚类算法的计算量；然后，采用CEM（Component-wise EM，模型自动选择最大期望）算法聚类，计算scatter separability策略值，通过该策略值来评价选择的特征子集和聚类效果；重复上述步骤，直到遍历完所有特征子集，求得策略值最大的特征子集。最后得到第i个极值点属于第j个聚类的后验概率： 

其中，

为第i个极值点的波形拟合特征系数，

表示第i个极值点的第k个波形拟合特征系数，N为拟合的阶数；c_j表示第j个聚类；μ_l为第l个层位的N×1的波形拟合特征系数的均值向量，

n1为第l个层位上的极值点数目；μ_j为第j个层位的波形拟合特征系数的均值向量；Σ_j为第j个层位的N×N的协方差矩阵，∑_j=cov((X^(j))^T,X^(j))；参数α_l表示第l个层位的分布概率，取值范围为0≤α_l≤1，

参数α_j表示第j个层位的分布概率，取值范围为0≤α_j≤1； K是层位数目。其中

通过这一步，将第三步得到的波形拟合系数分成若干个聚类，也就相当于将地震数据中的极值点成了若干个聚类，每一个聚类代表一个层位。以上三步基本实现了层位追踪，但是由于第三步中会产生误差，可能将一个层位分成多个类，即层位有中断，还需要后面的优化步骤。 

第四步、倾角计算：本发明中将倾角作为层位追踪的指向功能引入到地震图像中的层位信号分类方法研究中，处理聚类引起的误差以及追踪结果的不完整性，首先，计算地震图像的梯度向量g；然后，利用梯度向量g计算张量矩阵T；最后，计算张量矩阵的特征值λ_i,i=1,2,3，并求得倾角dip。 

通过这一步，得到了地震图像中的倾角属性，用于第五步层位补齐。 

第五步、层位补齐：聚类将地震图像中的极值点分为若干个聚类，但同一层位上的极值点可能被分到不同的类里面，这样得到的层位存在间断。在本发明中我们将相关查找方法引入到层位缺口的补齐中，具体流程如图4所示，其中相关系数门限取值范围为0到1，本发明方法中取0.3。 

第六步、重叠消除：由于特征向量的分类是基于相似性进行的，所以在分类的结果中很难避免将本属于两个不同层位的，但是波形很相似的层位点划分在同一类中。对于不同层位划分到同一类中出现重叠的问题，本发明采用的方法流程如图5所示，其中半径W曲子范围为5-30，本发明方法中取15。 

第七步、层位片段融合：在实际采集到的地震数据中，可能存在信号弱的区域，造成这些区域的层位追踪中断，即使第五步层位补齐也对这些区域无效。对于这种情况在本发明中提出了一种层位片段融合的方法，其主要流程如图6所示。 

图8是四川某工区地震图像的人工层位追踪效果，图9是本发明方法追踪效果。虽然相对人工追踪效果，本发明方法得到的层位连续性较差，但是仅存在于信号较弱的地方，对于这些区域由于不能提取出有用的极值点的特征信息，从而才导致了这些区域层位的中断。对于其他区域，本发明方法能准确高效的实现层位数据的提取，由于人工追踪的主观因素造成误差，因此较人工追踪的方法本发明方法具有更高的精确性。 

2、本发明三维地震图像中的全层位追踪方法。 

本发明将三维空间层位追踪问题转化为一系列二维平面内的层位追踪问题，最后通过组合二维剖面中的结果实现了三维图像中的层位追踪，即先在三维地震图像中每一个二维剖面中追踪层位，最后再将所有追踪得到的层位组合起来，得到三维地震图像中的层位。本发明中提出了一种基于半监督分类的三维地震数据的单一层位信号提取算法，算法流程如图9所示。 

第一步、剖面提取：本发明将三维层位追踪转换为二维剖面上的追踪，因此在截取剖面的时候需要考虑到为每一个剖面提供种子点（一般都是隔10个像素点左右一个种子点），以便在二维剖面上引导层位曲线的追踪。若提供的种子层位线是垂直于xline_time方向，则提取的剖面方向垂直于inline_time方向；类似地，若提供的种子层位线垂直于inline_time方向，则提取的剖面方向垂直于xline_time方向。 

第二步、种子点修正：由于第三步使用的是半监督聚类，需要给定种子点。给定的种子点一般都是解释人员根据工作经验人工标定的，其准确性不高，往往偏离极值点。因此在追踪算法使用种子点之前需要对种子点的位置进行修正，根据规则将种子点移动到适当 的极值点上。 

根据层位追踪遵循的准则，种子点修正方案如下：在种子点指定窗口大小内上下各搜索一个极值点作为拟修正位置，然后比较这两个拟修正位置与种子点的距离，选择距离近的拟修正位置；如果距离相等，根据前一道种子点修正的方式进行种子点修正：如果前一道修正后的种子点深度值更大，则将现有道的种子点往深度值更大的拟修正位置修正，如果前一道修正后的种子点深度值更小则修正方向相反。若在窗口内不能找到极值点可以作为种子点修正后的位置，则认为层位在该道数据上是没有点的，层位在该道上被断开。 

第三步：波形特征提取及半监督聚类：流程图如图11所示，极值查找和特征值提取的方法采用二维地震图像中的层位追踪相同的方法，查找极值的局部区域划定为所有相邻两个种子点上下偏移一定距离所围成的平行四边形组合而成的区域；标识种子点，将种子点特征和查找得到的极值点特征值一起组成待聚类的样本集合；对样本集合聚类，聚类方法同二维中的聚类方法，通过公式（1）确定种子点所属的类型，并将种子点所属的类型作为目标层位所对应的类型，选择的类型为种子点分布最多的那个类。 

以上三步基本实现了层位追踪，但是由于第三步中会产生误差，追踪得到的层位可能有中断，因此还需要第四步将中断处补齐。 

第四步、层位补齐（见图7）：对于每个剖面，选择两个相邻种子点作为起点和终点，分别从前往后和从后往前两个方向通过相关查找的方法进行补齐，由于从后往前进行查找时不能使用倾角，因此在此我们通过直接平移的方法，补齐方法与二维层位补齐相同。 

第五步、组合：在提取出待识别层位在所有二维剖面上的位置后，再将所得的所有曲线组合形成三维地震图像中的三维层位曲面，获得最终结果。 

图12是三维图像人工层位面追踪效果，图13是本发明方法追踪出来的层位面，可以发现本发明方法得到的追踪结果更加精确。 

Claims (10)

1.一种基于信号分类的地震图像层位追踪方法，其特征在于，该方法包括二维地震波形层位追踪方法和三维地震图像层位追踪方法： 

所述二维地震波形层位追踪方法为： 

1)寻找二维地震波形中的极值点：若二维地震波形中的某个点的值大于该点两边的点的值，且该点的值大于门限值，则该点为极大值点；若二维地震波形中的某个点的值小于该点两边的点的值，且该点的值小于门限值，则该点为极小值点； 

2)以上述各极值点为中心分别提取出一段二维地震波形，将切比雪夫零点作为拟合样本点计算二维地震波形拟合特征系数； 

3)使用SFS方法删选上述二维地震波形拟合特征系数，产生二维地震波形拟合特征系数的特征子集； 

4)采用CEM方法对二维地震波形拟合特征系数的特征子集聚类，计算散步分离性策略值； 

5)重复步骤4），直到遍历完所有特征子集，求得散步分离性策略值最大的特征子集，得到第i个极值点属于第j个聚类的后验概率P(c_j|X⁽ⁱ⁾)，若P(c_l|x⁽ⁱ⁾)>P(c_j|X⁽ⁱ⁾),即第i个极值点属于第l个聚类的后验概率最大，则第i个极值点属于第l个层位，从而将步骤1）中的极值点分为若干个聚类，所述若干个聚类组成一个类集合，类集合中的每一个聚类代表一个层位，实现二维地震波形的层位追踪；其中： 

其中，

为第i个极值点的波形拟合特征系数， 

k=1,2,...N表示第i个极值点的第k个波形拟合特征系数，N为拟合的阶数；c_j表示第j个聚类，1≤c_j≤K；μ_l为第l个层位的N×1的波形拟合特征系数的均值向量，

n1为第l个层位上的极值点数目；μ_j为第j个层位的波形拟合特征系数的均值向量；Σ_j为第j个层位的N×N的协方差矩阵，∑_j=cov((X^(j))^T,X^(j))；参数α_l表示第l个层位的分布概率， 取值范围为0≤α_l≤1，

参数α_j表示第j个层位的分布概率，取值范围为0≤α_j≤1；

表示所有层位的混合高斯分布的概率密度函数，其中

j=1,2,3...,K，

表示第j个层位的高斯分布的概率密度函数；K是层位数目； 

所述三维地震图像层位追踪方法为： 

1)截取三维地震图像上的一个二维剖面，并为该二维剖面提供至少两个种子点：若提供的种子点连线是垂直于xline_time平面的，则沿着inline方向提取二维剖面；若提供的种子层位线垂直于inline_time的平面，则沿着xline方向提取二维剖面； 

2)在上述种子点指定窗口大小内上下各搜索一个极值点作为拟修正位置，然后比较这两个拟修正位置与种子点的距离，选择距离近的拟修正位置，并将所述种子点移动到该拟修正位置；如果距离相等，且待修正的种子点为第一个种子点，则选取待修正的种子点指定窗口上深度值低的极值点作为拟修正位置，并将待修正的种子点移动到该拟修正位置，若待修正的种子点不是第一个种子点，则根据前一个种子点修正的方式进行种子点修正：如果前一二维剖面上修正后的种子点深度值更大，则将现有二维剖面的种子点往深度值更大的拟修正位置修正；如果前一二维剖面修正后的种子点深度值更小，则修正方向相反；若在指定窗口内不能找到极值点可以作为种子点修正后的位置，则认为层位在现有二维剖面数据上是没有点的，层位在现有二维剖面上被断开；所述指定窗口大小为种子点上下各15个点； 

3)在二维剖面中所有相邻两个修正后的种子点上下各偏移15个点所围成的平行四边形组合而成的区域内寻找极值点：若二维剖面中的某个点的值大于该点所在二维剖面两边的点的值，且该点的值大于门限值，则该点为极大值点；若二维剖面中的某个点的值小于该点所在二维剖面两边的点的值，且该点的值小于门限值，则该点为极小值点； 

4)以步骤2）修正后的种子点或步骤3）得到的极值点为中心提取出一段二维剖面波形，将切比雪夫零点作为拟合样本点计算二维剖面波形拟合特征系数； 

5)采用CEM方法对二维剖面波形拟合特征系数特征子集聚类，计算散步分离性策略值； 

6)重复步骤5），直到遍历完所有特征子集，求得散步分离性策略值最大的特征子集，得到第m个极值点或种子点属于第n个聚类的后验概率Q(c_m|Y⁽ⁿ⁾)，若Q(c_l'|Y^(m))>Q(c_n|Y^(m)),即第m个极值点或种子点属于第l'个聚类的后验概率最大，则第m个极值点或种子点属于第l'个层位，从而将步骤3）中的极值点和步骤2）中修正后的种子点分为若干个聚类；其中： 

其中，

为第m个极值点或种子点的波形拟合特征系数，M为拟合的阶数，

k=1,2,...M表示第m个极值点或种子点的第k个波形拟合特征系数，M是拟合阶数；c_n、c_l'分别表示第n个、第l'个聚类，μ_l'为第l'个层位的M×1的波形拟合特征系数的均值向量，

k为第l个层位上的极值点或种子点数目；μ_n为第n个层位的波形拟合特征系数的均值向量；Σ_n为第n个层位的M×M的协方差矩阵，Σ_n=cov((Y⁽ⁿ⁾)^T,Y⁽ⁿ⁾)；参数α_l'表示第l'个层位的分布概率，取值范围为0≤α_l'≤1，

参数α_n表示第n个层位的分布概率，取值范围为0≤α_n≤1；

表示所有层位的混合高斯分布的概率密度函数，其中

n=1,2,3...,T，

表示第n个层位的高斯分布的概率密度函数；T是层位数目； 

7)选择修正后种子点最多的聚类作为层位，删去其余的聚类，完成了一个二维剖面上的层位追踪； 

8)重复三维地震图像层位追踪方法的步骤1）～步骤7），直到追踪完三维地震图像上所有的二维剖面。 

2.根据权利要求1所述的基于信号分类的地震图像层位追踪方法，其特征在于，所述门限值为0。 

3.根据权利要求1或2所述的基于信号分类的地震图像层位追踪方法，其特征在于，所述二维地震波形层位追踪方法还包括以下步骤： 

6)计算二维地震波形的梯度向量g： 

其中，g_x、g_y、分别表示x、y方向的梯度向量，u(x,y)表示二维地震图像中坐标为(x,y)点所对应的像素值； 

7)利用梯度向量g计算张量矩阵T，T=g□g^T； 

8)计算上述张量矩阵T的特征值λ_i，求得二维地震波形的倾角dip；其中i=1,2； 

其中，v(x)、v(y)为张量矩阵T的特征向量； 

9)利用所述倾角dip进行二维地震波形层位补齐，消除层位间断； 

10)消除每一个聚类中不同层位上的极值点； 

11)对层位进行片段融合。 

4.根据权利要求1或2所述的基于信号分类的地震图像层位追踪方法，其特征在于，所述三维地震图像层位追踪方法还包括以下步骤： 

9)对于每一个二维剖面，选择相邻的两个种子点作为起点和终点，分别从前往后和从后往前两个方向通过相关查找的方法进行层位补齐； 

10)在提取出修正后的种子点所在层位在所有二维剖面上的位置后，将所得的二维剖面组合形成三维地震图像中的三维层位曲面。 

5.根据权利要求3所述的基于信号分类的地震图像层位追踪方法，其特征在于，所述步骤9）中，二维地震波形层位补齐的具体步骤如下： 

1)设定r=1； 

2)判断r是否小于等于类集合中的聚类的总数，若是，则进入3）；否则，进入8）； 

3)判断第r类对应的层位中是否存在缺口，若是，进入4）；否则，r的值加1，返回2）； 

4)查找层位中缺口处的层位点； 

5)根据所述倾角dip的指向查找与步骤4）中层位相邻的二维地震波形上的极值点； 

6)以上述步骤4）缺口处层位点和上述步骤5）中各极值点为中心，上下各取15个点提取地震波形片段，分别计算以各极值点为中心的地震波形片段与以所述缺口处层位点为中心的地震波形片段的相关系数，将该相关系数作为各极值点与所述缺口处层位点的相关系数；所述相关系数取值范围为0～1； 

7)判断是否存在大于门限值的相关系数，若是，则选择距离所述缺口处层位点最近的极值点添加到第r类对应的层位中，并将距离所述缺口处层位点最近的极值点作为新的缺口处的层位点，返回4）；否则，将未能补齐的层位划分出来单独作为一个层位，并将聚类的总数加1，返回3）； 

8)结束。 

6.根据权利要求5所述的基于信号分类的地震图像层位追踪方法，其特征在于，所述相关系数取0.3。 

7.根据权利要求3所述的基于信号分类的地震图像层位追踪方法，其特征在于，所述步骤10）中，消除每一个聚类中不同层位上的极值点的具体步骤为： 

1)设定s=1； 

2)判断s是否小于等于类集合中的聚类的总数，若是，则进入3）；否则，进入7）； 

3)判断第s类对应的层位是否有重叠，若有，则进入4）；否则，s的值加1，返回2）； 

4)选取第s类对应的层位上第一个二维地震波形上的点S作为原点； 

5)根据所述原点搜索所述起始道的相邻道上半径为W的窗口内的所有极值点集合D；所述半径W取值范围为5～30； 

6)判断D中是否含有与S同类的极值点S′，若是，则将S′作为原点，返回5）；否则，将第s类中的剩余点作为新的类，并将该新的类作为类集合的最后一个类，将聚类的总数加1，返回3）； 

7)结束。 

8.根据权利要求7所述的基于信号分类的地震图像层位追踪方法，其特征在于，所述半径W取15。 

9.根据权利要求3所述的基于信号分类的地震图像层位追踪方法，其特征在于，所述步骤11）中，对层位进行片段融合的具体步骤如下： 

1)设定h=1； 

2)判断h是否小于等于类集合中的聚类的总数，若是，则进入3）；否则，进入6）； 

3)判断第h类最后一个点的搜索窗口内是否存在其他类，若是，则进入4）；否则，h的值加1，返回2）；所述搜索窗口大小为第h类最后一个点上下各15个点； 

4)合并第h类和第t类中的极值点，作为新的类r； 

5)在类集合中删除第t类，聚类的总数减1，返回3）； 

6)结束。 

10.根据权利要求4所述的基于信号分类的地震图像层位追踪方法，其特征在于，对每一个二维剖面进行层位补齐的具体步骤如下： 

1)设定p=1； 

2)判断p是否小于种子点数目，若是，则进入3）；否则，进入9）； 

3)以第p个种子点和第p+1个种子点为两端点； 

4)由两个端点同时向中间查找是否存在缺口；若存在缺口，则进入5)；否则，将p的值加1，返回2）； 

5)将缺口处层位点平移到与所述缺口处层位点所在的二维剖面地震道相邻的二维剖面地震道上，在相邻的二维剖面地震道上，查找平移后的层位点上下各15个点范围内的极值点； 

6)以缺口处层位点和上述步骤5）中各极值点为中心，上下各取15个点提取地震波形片段，分别计算以各极值点为中心的地震波形片段与以所述缺口处层位点为中心的地震波形片段的相关系数，将该相关系数作为各极值点与所述缺口处层位点的相关系数；所述相关系数取值范围为0～1； 

7)判断是否存在大于门限值的相关系数，若是，选择与缺口处层位点相关系数最大的极值点，将该点作为新的缺口处层位点和新的端点，返回4）；否则，将p的值加1，返回2）； 

8)结束。 

Images