基于无监督分类的二维地震数据全层位自动追踪方法

摘要

提供了一种基于无监督分类的二维地震数据全层位自动追踪方法，所述方法包括：读入二维地震数据以及控制层位数据；查找二维地震数据波形的极值点；基于查找的极值点，通过切比雪夫多项式来拟合地震波形，将拟合系数作为特征向量；对拟合的地震波形进行基于特征向量的无监督分类，从而得到二维地震数据的全层位自动追踪结果。所述方法还包括：通过相关查找来补齐间断的地震层位。所述方法还包括：通过使用层位片段融合方法来得到完整的层位线。所述方法还包括：通过相关查找来消除重叠现象。所述方法消除了人工干预机制，实现完全的二维地震数据全层位自动追踪。

说明书

技术领域

本发明属于地震勘探资料解释领域，具体地说，涉及一种二维地震数据的层位追踪方法。

背景技术

层位追踪（即，层位解释）是地质资料解释的重要环节。层位追踪就是通过地震勘探得到的地震数据来分析得到地底构造。长期以来地震层位追踪都是依靠人工提取的。但是人工提取存在较大的问题：首先，人工提取依赖于解释人员的长期以来的工作经验，主观因素影响很大；其次是人工解释存在很大的效率问题，只能对少数目地层位进行解释，无法实现对所有层的追踪，通常难以为精细的地震资料分析（如地震地层学解释）提供基础数据。随着层位追踪方法的不断研究和改进，依赖于计算机实现的层位追踪算法不断提出，追踪的效率和效果都有所改进。

现有的层位追踪方法包括：（1）P.Alberts等在2002年提出了一种基于人工神经网络层位追踪算法，主要将模式识别引入到跨断层等不连续地质结构的层位追踪，然后采用神经网络进行模式识别；（2）Reda Benbernou等在P.Alberts的工作基础上采用模糊手段进行判决，组成混合的层位自动追踪方法；（3）M.Aurnhammer等在2002年提出了一种遗传算法层位追踪算法，其具体思路是采用基于模型的方法来处理跨越断层的问题，然后将问题转化为约束最优化的问题，然后使用遗传算法来求解；（4）F.Admasu等在2004年采用了模拟退火的方法来解约束最优问题，在2006年讨论了贝叶斯方法求解约束最优化问题；（5）F.Admasu等在2006年将子波变换引入了层位追踪问题中，先将地震数据小波多尺度分解，然后再进行层位追踪，依然采用贝叶斯方法；（6）Yingwei Yu等在2011年使用定向矢量场获得层位极值信息，将目标层位建模为无向连通图，然后使用最小生成树获得目标层位，但这样思路并不能进行全层位追踪；（7）Hilde G.Borgos等在2005年引入基于有限混合高斯的全层位追踪算法，由于将层位追踪转化为分类（包括有监督分类和半监督分类）问题，所以不存在跨断层问题，能适合复杂的地质环境。

然而，现有的层位追踪方法存在以下问题：（1）P.Alberts和RedaBenbernou的神经网络方法跨越复杂地质环境能力与训练样本状况密切相关，如果训练样本包含这种复杂地质状况，则能很好追踪，否则不能，但对于多变的复杂地质环境，训练样本往往无法包含全部的状况，而且训练样本的获得需要大量人工干预标示；（2）M.Aurnhammer和F.Admasu的模型方法存在对复杂地质状况的逼近程度的问题以及求解方法的求解精度和次优解问题；（3）有监督分类主要有最大似然和贝叶斯分类，使用大量的手工标识的数据进行训练，然后用未标识的数据来测试其性能，分类器具有很高的分类精度，半监督分类相对于有监督分类需要标识的样本少了很多，但是相应的分类的精度会有所下降，但是基于有监督分类和半监督分类的分类方法需要人工进行干预，不能实现完全的自动全层位追踪。

因此，需要一种基于无监督分类的二维地震数据全层位自动追踪方法。

发明内容

根据本发明的一方面，提供了一种基于无监督分类的二维地震数据全层位自动追踪方法，所述方法使用无监督分类来实现二维地震数据的全层位自动追踪。

根据本发明的另一方面，提供了一种基于无监督分类的二维地震数据全层位自动追踪方法，所述方法使用相关查找或层位片断融合来实现层位补齐或消除重叠。

根据本发明的一方面，提供了一种基于无监督分类的二维地震数据全层位自动追踪方法，所述方法包括：读入二维地震数据以及控制层位数据；查找二维地震数据波形的极值点；基于查找的极值点，通过切比雪夫多项式来拟合地震波形，将拟合系数作为特征向量；对拟合的地震波形进行基于特征向量的无监督分类，从而得到二维地震数据的全层位自动追踪结果。

所述极值点可包括极大值或极小值或者可查找二维地震数据波形的过零点并基于查找的过零点通过切比雪夫多项式来拟合地震波形。

对拟合的地震波形进行基于特征向量的无监督分类的步骤可包括：根据特征向量的概率密度函数建立统计模型：

$>f\left(a\right|c,u,\Sigma )=\underset{k\in ϵ}{\Pi }f\left(a_k\right|c_k,u_{c_k},{\Sigma }_{c_k})$>

$>=\underset{k\in ϵ}{\Pi }{\left(2\pi \right)}^{-n_A/2}{\left|{\Sigma }_{c_k}\right|}^{-1/2}\exp \{{-(a_k-u_{c_k})}^{\prime }{{\Sigma }^{-1}}_{c_k}(a_k-u_{c_k})/2\}---\left(3\right),$>其中，是特征向量A_k的边缘分布，c_k∈{1,2,…,n_c}表示分类的类标识，总共分为n_c类，表示高斯分布的均值，表示相应的方差；采用最大似然方法对特征向量进行聚类分析：

$>\mathrm{ML}=\underset{\theta }{\max }\left(f\left(A\right|\theta )\right)=\underset{\theta }{\max }\left(\underset{k\in ϵ}{\Pi }f\left(a_k\right|c_k,u_{c_k},{\Sigma }_{c_k})\right)---\left(4\right),$>通过对公式（4）进行求导得到：

$>\widehat{u_c}=\frac{1}{\left|{ϵ}_c\right|}\underset{k\in {ϵ}_c}{\Sigma }{A}_k---\left(5\right),$>

$>\widehat{{\Sigma }_c}=\frac{1}{\left|{ϵ}_c\right|}\underset{k\in {ϵ}_c}{\Sigma }{(A_k-\widehat{u_c})}^2---\left(6\right),$>

$>\widehat{c_k}=\arg \underset{c_k}{\max }f\left(a_k\right|c_k,u_c,{\Sigma }_c)---\left(7\right),$>其中，和分别为u，∑和c的估计值；通过EM最大期望算法不断迭代得到相应的估计值：

$>\widehat{u_c}=\frac{1}{\left|ϵ\right|}\underset{k\in ϵ}{\Sigma }f\left(c\right|A_k,\widehat{u_c},\widehat{{\Sigma }_c})A_k---\left(8\right)$>

$>\widehat{{\Sigma }_c}=\frac{1}{\left|ϵ\right|}\underset{k\in ϵ}{\Sigma }f\left(c\right|A_k,\widehat{u_c},\widehat{{\Sigma }_c}){(A_k-\widehat{u_c})}^2---\left(9\right),$>通过公式（8）和（9）得到各个聚类团的均值和协方差矩阵；将各个分类样本点映射到各自对应的聚类团中，即将各个极值点分配到不同的类型中，这些类型对应于初始得到的各个层位片段所经过的层位点，从而得到二维地震数据的全层位自动追踪结果。

所述方法还可包括：在对拟合的地震波形进行基于特征向量的无监督分类后得到的层位存在间断的情况下，通过相关查找来补齐间断的地震层位。

通过相关查找来补齐间断的地震层位的步骤可包括：从断开处的层位点，依据倾角的指示，找到下一道上的指示点，搜索该指示点两端一定长度的矩形窗口，找到窗口内相关系数大于设定门限的极值点，并且选择距离最近的极值点添加到本层位点集中，直到补齐缺口或者搜索窗口中不存在相关系数大于门限的极值点，其中，所述相关系数表示指示点至断开处的层位点的向量与窗口内的极值点至断开处的层位点的向量的振幅相关系数。

所述方法还可包括：在无法在窗口内找到相关系数大于设定门限的极值点的情况下，通过使用层位片段融合方法来得到完整的层位线。

通过使用层位片段融合方法来得到完整的层位线的步骤可包括：依次以各个层位片段的最后一个层位点为中心确定一个矩形窗口，然后检查在该窗口内是否存在其他层位片段的起始层位点，若存在则将这两个层位片段连接形成新的层位片段，以这种方法不断的检测，直到所有的层位片段最后一个层位点的窗口内都不存在另外一个层位片段的起始层位点，最终得到完整的层位线。

所述方法还可包括：在对拟合的地震波形进行基于特征向量的无监督分类后得到的层位存在重叠现象的情况下，通过在预定窗口中查找最近的同类层位点来消除重叠现象。

通过在预定窗口中查找最近的同类层位点来消除重叠现象的步骤可包括：预设一个窗口，从本类的起始道的层位点开始，将该层位点平移到相邻道上，查找该道上同一类且距离平移的点最近的层位点划分到本层位点集中，然后继续查找下一道上的层位点，直到找完整个层位，剩下的点集作为另一个聚类团，即另外一个层位的数据点。

根据本发明的基于无监督分类的二维地震数据全层位自动追踪方法的方法消除了人工干预机制，实现完全的二维地震数据全层位自动追踪。所述方法追踪层位的精度高、速度快，有效地解决了地震解释中存在的效率低的问题，能够为层序地层学研究提供地质等时格架。所述方法能够广泛地应用于石油地震勘探、地质地矿和建筑工程地况调查。

附图说明

通过结合附图，从实施例的下面描述中，本发明这些和/或其它方面及优点将会变得清楚，并且更易于理解，其中：

图1是根据本发明的基于无监督分类的二维地震数据全层位自动追踪方法的流程图；

图2是根据本发明的地震波形极值点的示意图；

图3是根据本发明的通过切比雪夫多项式来拟合地震波形的示意图；

图4是根据本发明的基于无监督分类得到的地震层位存在间断的示意图；

图5是根据本发明的通过相关查找来补齐图4所示的间断的地震层位的示意图；

图6是根据本发明的通过层位片段融合来补齐地震层位的示意图；

图7是根据本发明的基于无监督分类得到的地震层位出现重叠现象的示意图；

图8根据本发明的通过在预定窗口中查找最近的同类层位点来消除图7所示的重叠现象的示意图；

图9是通过使用本发明的基于无监督分类的二维地震数据全层位自动追踪方法而获得的全层位追踪的示例的效果图；

图10是通过使用本发明的基于无监督分类的二维地震数据全层位自动追踪方法而获得的全层位追踪的另一示例的效果图。

具体实施方式

提供参照附图的以下描述以帮助对由权利要求及其等同物限定的本发明的实施例的全面理解。包括各种特定细节以帮助理解，但这些细节仅被视为是示例性的。因此，本领域的普通技术人员将认识到在不脱离本发明的范围和精神的情况下，可对描述于此的实施例进行各种改变和修改。此外，为了清楚和简洁，省略对公知的功能和结构的描述。

图1是根据本发明的基于无监督分类的二维地震数据全层位自动追踪方法的流程图。

参照图1，在步骤101，读入二维地震数据以及控制层位数据。其中，二维地震数据可包括地震剖面的起始CDP和终止CDP以及各个深度的振幅值，可用一个二维数组表示。控制层位数据可以是控制层在地震剖面的各个CDP号对应下的深度值。

在步骤102，查找二维地震数据波形的极值点。

具体地说，由于层位追踪通常是在二维地震数据的极大值，极小值和过零点面上进行的，因此进行层位追踪的第一步需要找到这些极大值、极小值或者过零点。二维地震数据波形的极大值和极小值被称为地震极值。如图2所示，图2示出地震波形极值点的示例。地震极值可以定义为：

$>e\left(x\right)=\{t:\frac{\mathrm{dD}(x,t)}{\mathrm{dt}}=0\}---\left(1\right).$>

其中，D＝{D(x,t)}表示地震剖面，x为CDP号或者线号，t为双程旅行时或者深度，则D(x₀,t)可表示为一条地震道。在本发明的一个实施例中，主要使用地震极大值作为层位自动追踪的基础，但应该理解，本发明并不限于使用地震极大值作为层位自动追踪的基础，在本发明的其它实施例中，还可使用地震极小值、地震过零点来作为层位自动追踪的基础。

返回参照图1，在步骤103，基于查找的极值点，通过切比雪夫多项式来拟合地震波形。

具体地说，地震道可通过切比雪夫多项式拟合重建，得到相应的拟合系数，因此二维地震数据波形的极值点周围的波形可以通过拟合得到的系数表示。通过切比雪夫多项式拟合得到的极大值两端局部地震波形可表示为：

S(z)=b₀p₀(z)+b₁p₁(z)+…+b_Np_N(z)                （2）

其中，N表示拟合阶数，$>\left(\begin{array}{c}{p}_0\left(z\right)=1\\ p_1\left(z\right)=z\\ p_{n+1}\left(z\right)=2·z·p_n\left(z\right)-p_{n-1}\left(z\right)\end{array}\right)$>为切比雪夫多项式，b_j(j=0,1,..，N)为拟合所得系数，

$>b_0=\frac{1}{N+1}\underset{k=0}{\overset{N}{\Sigma }}f\left(z_k\right)*p_0\left(z_k\right)=\frac{1}{N+1}\underset{k=0}{\overset{N}{\Sigma }}f\left(z_k\right),$>

$>b_j=\frac{2}{N+1}\underset{k=0}{\overset{N}{\Sigma }}f\left(z_k\right)*p_j\left(z_k\right)=\frac{2}{N+1}\underset{k=0}{\overset{N}{\Sigma }}f\left(z_k\right)*\cos \left(\frac{\mathrm{j\pi }(2k+1)}{2N+2}\right)j=\mathrm{1,2},...,N.$>

在地震波形拟合过程中，S(z)表示某一道上一段波形，z为深度，即等式（1）中的t，S(z)为对应的振幅值。通过切比雪夫多项式拟合的方法重建地震道，得到拟合系数，用这些拟合系数特征化波形，并将拟合系数作为特征向量。如图3所示，图3示出根据本发明的通过切比雪夫多项式来拟合地震波形的示例。

返回参照图1，在步骤104，对拟合的地震波形进行基于特征向量的无监督分类，从而得到二维地震数据的全层位自动追踪结果。

具体地说，由于相同层位上的极值点周围的波形通常非常相似，因此可以通过将特征向量分为不同的类，来得到的聚类团。聚类团中的特征向量对应的极值点表示在同一类中，每个聚类团表示一个层位片段。

对特征向量的分类是建立在特征向量的统计模型的基础上的。特征向量的联合分布被认为服从高斯混合分布。因此，根据这些独立的特征向量的概率密度函数建立统计模型：

$>f\left(a\right|c,u,\Sigma )=\underset{k\in ϵ}{\Pi }f\left(a_k\right|c_k,u_{c_k},{\Sigma }_{c_k})$>

$>=\underset{k\in ϵ}{\Pi }{\left(2\pi \right)}^{-n_A/2}{\left|{\Sigma }_{c_k}\right|}^{-1/2}\exp \{{-(a_k-u_{c_k})}^{\prime }{{\Sigma }^{-1}}_{c_k}(a_k-u_{c_k})/2\}---\left(3\right),$>

其中，是特征向量A_k的边缘分布，c_k∈{1,2,…,n_c}表示分类的类标识，总共分为n_c类，表示高斯分布的均值，表示相应的方差。

随后，可采用最大似然方法对特征向量进行聚类分析：

$>\mathrm{ML}=\underset{\theta }{\max }\left(f\left(A\right|\theta )\right)=\underset{\theta }{\max }\left(\underset{k\in ϵ}{\Pi }f\left(a_k\right|c_k,u_{c_k},{\Sigma }_{c_k})\right)---\left(4\right),$>

通过对公式（4）进行求导可得到：

$>\widehat{u_c}=\frac{1}{\left|{ϵ}_c\right|}\underset{k\in {ϵ}_c}{\Sigma }{A}_k---\left(5\right),$>

$>\widehat{{\Sigma }_c}=\frac{1}{\left|{ϵ}_c\right|}\underset{k\in {ϵ}_c}{\Sigma }{(A_k-\widehat{u_c})}^2---\left(6\right),$>

$>\widehat{c_k}=\arg \underset{c_k}{\max }f\left(a_k\right|c_k,u_c,{\Sigma }_c)---\left(7\right),$>

由于和分别为u，∑和c的估计值，且彼此之间是相互关联的，因此需要通过EM（最大期望）算法不断迭代得到相应的估计值：

$>\widehat{u_c}=\frac{1}{\left|ϵ\right|}\underset{k\in ϵ}{\Sigma }f\left(c\right|A_k,\widehat{u_c},\widehat{{\Sigma }_c})A_k---\left(8\right)$>

$>\widehat{{\Sigma }_c}=\frac{1}{\left|ϵ\right|}\underset{k\in ϵ}{\Sigma }f\left(c\right|A_k,\widehat{u_c},\widehat{{\Sigma }_c}){(A_k-\widehat{u_c})}^2---\left(9\right)$>

分类后可得到各个聚类团的均值和协方差矩阵等相关信息，然后将各个分类样本点映射到各自对应的聚类团中，即将各个极值点分配到不同的类型中，而这些类型则对应于初始得到的各个层位片段所经过的层位点，从而可得到二维地震数据的全层位自动追踪结果。

基于特征向量的无监督分类而得到的聚类团对应的点理论上表示一个地震层位，但是实际上可能由于地域跨度大，在一定地理位置间隔同一层位的两个点附近的波形可能相似性不高，因此在基于特征向量的无监督分类中会被分到不同的聚类团里面，这样得到的层位存在间断。如图4所示，图4示出根据本发明的基于无监督分类得到的地震层位存在间断的示例。

在这种情况下，需要根据本发明的相关查找的方式将出现间断的层位补齐。因此，根据本发明的基于无监督分类的二维地震数据全层位自动追踪方法还可包括：在步骤105，通过相关查找来补齐间断的地震层位。以下，将参照图5具体描述根据本发明的通过相关查找来补齐间断的地震层位的方法。

图5是根据本发明的通过相关查找来补齐图4所示的间断的地震层位的示意图。参照图5，对于同一层位的断开问题，从断开处的层位点E，依据倾角的指示，找到下一道上的指示点A，搜索A两端一定长度d的矩形窗口，找到窗口内相关系数大于设定门限（例如，Cov）的极值点，并且选择距离最近的极值点（例如，A2）添加到本层位点集中，直到补齐缺口或者搜索窗口中不存在相关系数大于门限的极值点，其中，所述相关系数表示指示点A至断开处的层位点E的向量与窗口内的极值点至断开处的层位点E的向量的振幅相关系数，例如，向量（A，E）与向量（A1，E）的振幅相关系数、向量（A，E）与向量（A2，E）的振幅相关系数等。

如果通过步骤105无法在窗口内找到相关系数大于设定门限的极值点，则无法通过所述相关查找来补齐间断的地震层位，此时，仍会出现间断的层位片段。在这种情况下，可采用根据本发明的层位片段融合方法来得到完整的层位线。因此，根据本发明的基于无监督分类的二维地震数据全层位自动追踪方法还可包括：在步骤106，通过层位片段融合来补齐地震层位。以下，将参照图6具体描述根据本发明的层位片段融合方法。

图6是根据本发明的通过层位片段融合来补齐地震层位的示意图。参照图6，某一层位具有三条地震道T_x-1、T_x和T_x+1。该层位在T_x-1和T_x+1道上存在层位点A和A2，但是在T_x道上通过步骤105时在窗口内无法找到相关系数大于设定门限的极值点。在这种情况下，依次以各个层位片段的最后一个层位点为中心确定一个边长为Dw的矩形窗口，然后检查在该窗口内是否存在其他层位片段的起始层位点，若存在则将这两个层位片段连接形成新的层位片段。以这种方法不断的检测，直到所有的层位片段最后一个层位点的窗口内都不存在另外一个层位片段的起始层位点，最终得到完整的层位线。

另外，由于特征向量的无监督分类是基于相似性进行的，所以在无监督分类的结果中很难避免将属于两个不同层位但波形很相似的层位点划分在同一类中。如图7所示，图7示出根据本发明的基于无监督分类得到的地震层位出现重叠现象的示例。

在这种情况下，需要根据本发明的通过在预定窗口中查找最近的同类层位点的方式消除重叠现象。因此，根据本发明的基于无监督分类的二维地震数据全层位自动追踪方法还可包括：在步骤107，通过在预定窗口中查找最近的同类层位点来消除重叠现象。以下，将参照图8具体描述根据本发明的通过在预定窗口中查找最近的同类层位点来消除重叠现象的方法。

图8是根据本发明的通过在预定窗口中查找最近的同类层位点来消除图7所示的重叠现象的示意图。参照图8，对于不同层位划分到同一类中出现重叠的问题，可以设定一个窗口，从本类的起始道的层位点S开始，将该点平移到相邻道上S’，查找该道上同一类且距离S’最近的层位点S1划分到本层位点集中，然后继续查找下一道上的层位点（即沿S→S1方向查找下一道上的层位点），直到找完整个层位，剩下的点集作为另一个聚类团，即另外一个层位的数据点。按照这种方法依次去除所有聚类团中的重叠现象。

图9是通过使用本发明的基于无监督分类的二维地震数据全层位自动追踪方法而获得的全层位追踪的示例的效果图。图10是通过使用本发明的基于无监督分类的二维地震数据全层位自动追踪方法而获得的全层位追踪的另一示例的效果图。通与人工追踪的层位进行对比结果表明，自动追踪的层位与人工解释的层位基本一致，证明了所述方法的有效性。

本发明提供了一种基于无监督分类的二维地震数据全层位自动追踪方法。所述方法使用无监督分类来实现二维地震数据的全层位自动追踪。此外，所述方法还使用相关查找或层位片断融合来实现层位补齐或消除重叠。所述方法消除了人工干预机制，实现完全的二维地震数据全层位自动追踪。所述方法追踪层位的精度高、速度快，有效地解决了地震解释中存在的效率低的问题，能够为层序地层学研究提供地质等时格架。所述方法能够广泛地应用于石油地震勘探、地质地矿和建筑工程地况调查。

可根据计算机程序指令来执行根据本发明的上述方法。由于这些程序指令可被包括在计算机、专用处理器或可编程或专用硬件中，因此在其中执行的指令可有利于上述的功能的执行。如本领域的技术人员所理解的，计算机、处理器或可编程硬件包括可存储或接收软件或计算机代码的存储器件，所述软件或计算机代码在被计算机、处理器或硬件访问和执行时实现在本发明中描述的方法。

虽然本发明是参照其示例性的实施例被显示和描述的，但是本领域的技术人员应该理解，在不脱离由权利要求及其等同物限定的本发明的精神和范围的情况下，可以对其形式和细节进行各种改变。