一种地震波形聚类方法与装置与流程

本发明属于地震勘探基础应用领域，具体涉及一种地震波形聚类方法与装置。

背景技术：

地震波形是一种综合属性参数，包含了振幅、频率、相位、几何形状等信息，用于聚类分析可以更为可靠地反映真实的地下地质情况；但是同样地震波形也包含了很多噪音和冗余信息，导致聚类结果受噪音干扰严重，分辨率较低。

无监督地震波形聚类是一种应用广泛的地震相分析方法，在地震沉积相分析和储层预测中发挥重要作用，常用的方法是将地震波形输入到自组织神经网络(self-organizingmap，som)中进行学习(learning)然后聚类。

作者matos等人提出了一种基于小波变换的无监督地震相分析方法(geophysics，2007，vol.72，no.1，p：9-21)。该方法首次将时频技术应用到了模式识别系统中，利用了wtmmla(wavelettransformmodulusmaximalineamplitudes)变换对地震信号的奇点进行了表征，再将整个wtmmla曲线作为输入数据进行som波形聚类，以提高聚类结果的解释容错性。该方法在模型实验中展示了较高的容忍解释错误的能力，但由于用wtmmla曲线替代了地震数据，大幅度减少了地震数据的原有信息，聚类结果可靠性差，在实际地震数据中的应用效果不明显。

作者saraswat和sen提出了一种基于ais的地震相分析方法(geophysics，2012，vol.77，no.4，p：45–53)，该方法利用ais(artificialimmunesystem)对地震数据进行了大幅度的缩减，去除地震数据的噪音和冗余信息，再将缩减以后的地震数据用于聚类分析，聚类结果不容易受到噪音和冗余数据的影响。该方法在实际地震数据(荷兰f3区块)聚类中，表现出较好的抗噪性，但ais缩减了地震数据50％的信息，在没有质控的情况下很难保证去除的都是噪音和冗余数据，也可能去除掉一些主体频带信息，可靠性差。

作者杜浩坤等人提出了一种基于emd的地震相分析方法(journalofappliedgeophysics，2015，no.112，p：52-61)。该方法利用emd技术将地震信号分解成有限的、频带范围不同的固有模态函数(intrinsicmodefunction,imf)，然后筛选出显示噪音较少、与原始地震信号相关性较高的分量进行重构，最后利用重构数据进行地震波形聚类。该方法具备较好的抗噪效果，但是该方法在对地震数据进行分解的过程中，分量的个数、频带范围由地震数据和算法的收敛条件决定，一般为7-10个分量，导致单个分量的频带范围太窄，与原始剖面存在较大差别，筛选和处理分量时缺乏有效的评判准则；另一方面，在对分量进行筛选时，对于噪音较多的分量都是直接丢弃，这些分量里面可能会包含一些有效信息。所以该方法在分解算法和分量处理思路上存在局限性，影响了聚类结果的可靠性。

综上所述，如何在保证可靠性的同时，提高波形聚类的抗噪性和分辨率，是近年来聚类分析一直研究的课题；将时频分析方法应用到无监督地震波形聚类中是最有潜力的研究方向之一，但目前使用的时频分析方法存在各自的局限性，仍需针对地震信号的特性进行改善。

技术实现要素：

本发明的目的是提供一种地震波形聚类方法与装置，在对地震数据进行预处理时，能够加强质量监控、优化分量筛选思路、去除噪声强点，保证聚类结果的抗噪性和可靠性，解决现有地震波形聚类方法无法同时保证聚类的可靠性和抗噪性的问题。

为解决上述技术问题，本发明提出一种地震波形聚类方法，包括以下解决方案：

方法方案一，包括如下步骤：

1)获取目标层地震数据中地震道信号的频带范围，将频带范围划分，至少包括第一设定频带范围和第二设定频带范围；其中，第一设定频带范围为地震道信号的主体频带范围，主体频带范围为振幅能量强度降至最大值的一半时的频带范围；第二设定频带范围为除去第一设定频带范围外、所述目标层地震数据在剖面上保持层状展布特征的频带范围；

2)将所述目标层地震数据分解，生成至少两个固有模态函数分量，包括第一固有模态分量和第二固有模态分量，其中，第一固有模态分量的频带范围为所述第一设定频带范围，第二固有模态分量的频带范围为所述第二设定频带范围；

3)对第二固有模态分量进行去噪处理，处理后与第一固有模态分量叠加重构，得到重构地震数据体，并对该重构地震数据体进行学习和聚类，得到波形地震相图。

方法方案二，在方法方案一的基础上，所述第一设定频带范围是对目标层地震数据进行频谱分析后，获取目标层的主体频带范围；所述第二设定频带范围的上限值为地震道信号的截止频率。

方法方案三、四，分别在方法方案一、二的基础上，所述第一设定频带范围是对目标层地震数据进行同步挤压小波变换、获取地震道信号的时频图后，确定的地震道信号的主体频带范围。

方法方案五，分别在方法方案一的基础上，采用同步挤压小波变换将所述目标层地震数据分解。

方法方案六，在方法方案一的基础上，所述学习为自组织神经网络学习。

为解决上述技术问题，本发明还提出一种地震波形聚类装置，包括以下解决方案：

装置方案一，包括处理器，该处理器用于执行实现以下方法的指令：

1)获取目标层地震数据中地震道信号的频带范围，将频带范围划分，至少包括第一设定频带范围和第二设定频带范围；其中，第一设定频带范围为地震道信号的主体频带范围，主体频带范围为振幅能量强度降至最大值的一半时的频带范围；第二设定频带范围为除去第一设定频带范围外、所述目标层地震数据在剖面上保持层状展布特征的频带范围；

2)将所述目标层地震数据分解，生成至少两个固有模态函数分量，包括第一固有模态分量和第二固有模态分量，其中，第一固有模态分量的频带范围为所述第一设定频带范围，第二固有模态分量的频带范围为所述第二设定频带范围；

3)对第二固有模态分量进行去噪处理，处理后与第一固有模态分量叠加重构，得到重构地震数据体，并对该重构地震数据体进行学习和聚类，得到波形地震相图。

装置方案二，在装置方案一的基础上，所述第一设定频带范围是对目标层地震数据进行频谱分析后，获取目标层的主体频带范围；所述第二设定频带范围的上限值为地震道信号的截止频率。

装置方案三、四，分别在装置方案一、二的基础上，所述第一设定频带范围是对目标层地震数据进行同步挤压小波变换、获取地震道信号的时频图后，确定的地震道信号的主体频带范围。

装置方案五，分别在装置方案一的基础上，采用同步挤压小波变换将所述目标层地震数据分解。

装置方案六，在装置方案一的基础上，所述学习为自组织神经网络学习。

本发明的有益效果是：本发明将原始地震数据分解后重构得到至少包括第一模态函数分量和第二模态函数分量，将第二摸态函数分量去噪处理后与第一模态函数分量叠加重构，对重构后的数据体进行学习和聚类，获得地震相图。本发明可以自定义设置模态函数分量的个数和频带范围，将对应主体频带范围的第一模态函数分量保留，将对应第二设定频带范围内噪音较多、但仍然表现为层状展布特征的第二模态函数分量不再直接丢弃，而是进行去噪处理后，与第一模态函数分量叠加，得到重构地震数据体用于生成地震相图，能同时保证地震波形聚类的可靠性和抗噪性，并有效去除了噪音和冗余信息。

进一步，为了得到精确的地震道的主体频带范围，本发明利用同步挤压小波变换，获取每一个地震道信号的时频图后，对每一地震道分别进行时频分析，确定每一地震道的主体频带范围和模态函数分量个数后，再逐道重构模态函数分量。但对实际地震资料进行同步挤压小波变换处理时，地震数据往往包含了几十万个地震道，为减少计算量，提高数据处理效率，所有地震道信号的主体频带范围均采用目标层段的主体频带范围，目标层段的主体频带范围是通过对目标层段的地震数据进行频谱分析后确定的。

进一步，本发明设置的第一设定频带范围为4-39hz，该频带范围对应的第一模态函数分量占据了地震信号的主体能量信息，保留且不采取任何措施。设置第二设定频带范围为40-80hz，该频带范围对应的第二模态函数分量表现为层状展布特征，去噪处理后保留其中的有用信息，将大于80hz无法保持层状展布特征的其他模态函数分量丢弃。在抑制噪音等信息的同时突出地震道信号主体能量信息的反映。

进一步，重构地震数据体进行学习前，设定了能刚好包含目标层段的信息的时窗大小，在该时窗内的波形信息不熟噪音影响，且不包含非目的层段的信息，不影响聚类结果的可靠性。还设置了大于沉积相的数目的聚类数目，充分、精确的反映了沉积相信息。

附图说明

图1是一种地震波形聚类方法流程图；

图2是川西地区目的层段频谱分析图，目的层主体频带为4-39hz；

图3是sst生成的目的层段单道时频分析图；

图4是sst分解的imf分量1示意图，频带范围为4-39hz；

图5是sst分解的imf分量2示意图，频带范围为40-80hz；

图6是sst分解的imf分量3示意图，频带范围为大于80hz；

图7是经过去噪处理以后的imf分量2示意图；

图8是经过sst处理后imf分量1叠加去噪后的imf分量2的重构地震剖面图；

图9是原始地震剖面图；

图10是本发明基于sst的无监督地震波型聚类图；

图11是基于emd的无监督地震波型聚类图。

具体实施方式

下面以川西地区雷口坡组顶部岩溶储层的地震波形聚类为例，结合附图对本发明进行进一步描述。本发明实施方法如图1所示，目的层段为雷口坡组顶部，包括步骤如下：

1、对目的层段的地震数据体进行频谱分析。由分析结果图2可知，该地区雷口坡组顶部地震数据体主体频带为4-39hz，评判标准是振幅能量强度降至最大值的一半时的频带范围；中心频率为21.5hz。

2、对工区内目标层段雷口坡组顶部的地震道信号逐道进行同步挤压小波变换(synchrosqueezingtransform，sst)，分析各道生成的时频图，重构生成自定义频带范围的imf分量。

2.1、对工区内目标层段的地震道信号逐道进行同步挤压小波变换(sst)，分析各道生成的时频图。同步挤压小波变换是2011年作者daubechies等人提出的一种结合小波变换与频率再分配算法的自适应非线性非平稳信号分解算法。

对地震道信号s(t)进行连续小波变换(continuouswavelettransform，cwt)：

式中，a为比例因子，b为时移因子，ψ^*为源小波的复合共轭，ws(a,b)系数代表了用来提取瞬时频率的集中时频图像。

据plancherel定理，方程(1)可写作：

式中，ξ代表角频率，和分别代表ψ(t)和s(t)的傅里叶变换，表示对ψ(t)的频谱取共轭，

考虑到一个单谐波信号s(t)＝acos(ωt)的傅里叶对可以表示为：

等式(2)可以变换为：

因为子波在中心频率ω0的附近提取，所以ws(a,b)将在水平线a＝ω0/ω附近提取，a＝ω0/ω表示小波的中心频率与信号的中心频率的比值。但实际上ws(a,b)总是围绕这条水平线展开，在时间尺度表达中会产生模糊投影。这种模糊主要发生在沿着时移因子b的刻度尺度上。沿着比例轴的维度b很少发现这种模糊。可以证明，当b维度中的模糊可以忽略时瞬时频率ωs(a,b)能利用下面等式计算出来：

对于任意点(a,b)，ws(a,b)≠0，将信息从时标平面映射到时频平面，并转换每个点(b,a)到(b,ωs(a,b))，这一步骤就叫做挤压同步。

因为a和b是离散值，缩放步骤δak＝ak-1-ak会遇到任意ak上。然后同步挤压变换的时频分布可以仅在频带范围[ωl-δω/2,ωl+δω/2]的中心ωl确定：

式中，δω＝ωl-ωl-1，ts(ωl,b)为在中心频率ωl处的sst时频分布。

2.2、根据sst时频图分析每道地震信号的主要能量分布范围。如图3所示，为工区内目的层段某单道地震信号的时频图，优选出该地震道信号主要能量的频带范围为6-42hz。为了减少噪声和冗余信息的影响，这里对该地震道信号sst分解后生成的固有模态函数(imf)重构生成三个imf分量。imf分量1确定为主体频带范围，是地震道信号能量的主体成分；imf分量2确定为高频成分，从地震道信号主体频带到主要能量的截止频率；大于截止频率的部分作为imf分量3。作为其他实施方式，可以对该地震道信号sst分解后生成imf分量1，和imf分量2，不生成imf分量3。

上述方法是确定频带范围的精确方法，对每一道都进行时频分析，得到主体频带信息以后再逐道重构imf。但对实际地震数据进行sst处理时，地震数据往往包含了几十万个地震道。为了减少计算量，提高数据处理效率，作为其他实施方式，所有地震道信号均采用步骤1中确定的主体频带范围为4-39hz，作为整个地震数据体统一使用的主体频带范围。图4为imf分量1，频带范围为4-39hz；图5为imf分量2，频带范围为40-80hz；图6为imf分量3，频带范围为大于80hz。

等式(6)表明，信号的新时频表示仅沿频率轴同步。cwt的系数由sst重新分配，获得时频平面上的图像。瞬时频率也从新的时频表示中提取出来。imf分量sk可以通过离散同步挤压变换的进行重构：

式中，是依赖于所选小波的常数，re(·)为表示的实部；是ts(ωl,b)的离散版本tm是离散时间，tm＝t0+mδt，m＝0,1,...,n-1；n是离散信号中的样本总数；δt是采样率。

3、对步骤2获取的imf分量进行处理：保留imf分量1的所有信息；对imf分量2进行去噪处理；去除imf分量3。然后将处理后的imf分量1和imf分量2相加，获得重构地震数据体。

3.1、imf分量2是高频分量，频带范围为40-80hz，具备较多的噪音，但是依然表现为层状展布特征，说明其包含了可用信息，可用信息的评判标准是在剖面上是否能保持层状展布特征。因此对imf分量2不再直接舍弃，而是进行去噪处理后保留，图7为去噪以后的imf分量2；imf分量3，大于80hz，展示出的剖面基本由噪音成分组成，无法保持层状展布特征，因此直接将其去除；imf分量1(4-39hz)占据了地震信号的主体能量信息，保留且不采取任何措施。这样的处理旨在抑制噪音等信息的同时突出地震道信号主体能量信息的反映。

上述分量处理主要是保留主体频带部分，对噪音较多的高频成分进行去噪处理，主体频带范围和高频成分频带范围会根据实际资料有所不同。

3.2、将步骤3.1获取的去噪以后的imf分量2与imf分量1相加，得到经过sst处理后的重构地震数据体，图8和图9为重构地震剖面与原始地震剖面的对比，能看到重构地震剖面保留了原始剖面的主要特征信息，减少了噪音并突出了主体能量信息的反映。

4、对目标层段雷口坡组顶选择合适的时窗和聚类数目，将步骤3中获取的重构地震信号进行som学习，学习完成以后再进行k均值聚类，获取地震波形聚类图。

4.1、对目标层段选择合适的时窗。时窗选择直接影响了地震波形聚类的结果，过大的时窗会包含非目的层段的信息，影响聚类结果；过小的时窗包含的波形信息不全且容易受到噪音的影响。理想的情况是选取的时窗能刚好包含目标层段的信息；雷口坡组顶部岩溶储层一般发育在不整合风化壳之下90m范围内，按6000m/s的速度计算，时窗应选择雷口坡组顶部(不整合面)下延30ms内。

4.2、选取合适的聚类数目。由于地震相信息反映的是沉积相信息，所以选取的地震相数目应大于沉积相数目，一方面是为了更充分、精确的反映沉积相信息；一方面是噪音或冗余信息会占据一种或几种地震相，因此聚类数目宜多不宜少，一般为5-7类。根据工区内雷口坡组顶部的沉积相和地震相特征，将波形聚类数目定为7类。

4.3、将步骤3.2获取的经sst处理后的重构地震数据体进行som无监督学习，获得学习完成后的原型向量。

将地震数据x＝[x1,x2,...xn]用n维空间向量来表示，并输入到二维网格中开始学习。这个二维网格是由原型向量(prototypevector)组成。原型向量的排列方式为矩形或六边形；其个数可以自由设定；其维度与地震信号的维度相同。

定义二维网格中原型向量的个数为y，b＝[b1,b2,...bi]，i＝1,2,…,y,并且所有输入的原型向量呈六边形排列相互连接。将具有与输入向量x最小距离的原型向量表示为bm，这一原型向量叫做最佳匹配单元(bestmatchunit)，可以利用等式8寻找bm：

||x-bm||＝min||x-bi||(8)

更新最佳匹配单元对应的原型向量及其相邻原型向量，朝着属于空间里的优胜向量移动，更新规则为：

bi(t+1)＝bi(t)+λ(t)hbi(t)[x-bi(t)](9)

λ(t)为学习参数，hbi(t)是以最佳匹配单元的邻域函数，表示为：

σ(t)是相邻向量的有效宽度，||rb-ri||²为优胜向量bm及其相邻向量bi之间的横向距离。

这些步骤将迭代进行，直到学习结束。

4.4、将步骤4.3中获取的学习完成的原型向量输入地震数据进行分类，得到无监督地震波形聚类图，即波形地震相图。通常使用k均值聚类，定义k个类簇点，目标函数满足：

式中，这里n代表输入数据点和类簇点的维度。

最终基于sst的无监督地震波形聚类结果如图10所示，图11为基于emd的无监督波形聚类结果图。比较图10和图11可以看到，利用本发明获取的地震波形聚类结果噪音更少，断层等细节刻画更为精确，同时与原始地震数据的聚类结果相符度高，可靠性好。认为本发明充分利用了sst的优势，自主定义了imf分量的个数和频带范围，并针对各分量的特点进行了处理，同时保持了良好的质量监控，地震波形聚类结果表现出较好的可靠性和抗噪性。

本发明充分利用了用同步挤压小波变换方法的优势，在将地震道信号分解为固有模态函数分量的过程中，针对不同地震数据、不同层段的频谱特征，可以自主定义分量的个数和频带范围。同时优化固有模态分量的处理，对噪音多的高频分量不再直接遗弃，而是进行去噪处理。

分解算法(同步挤压小波变换方法)和分量处理思路的改进，有助于更好的保留或突出主体能量的信息，去除噪音和冗余信息。同时，能提供有效的质量监控上，重构过程中筛选/丢弃的分量和重构地震数据体均能成图进行质量监控，检查去除的分量是否都为噪音或冗余信息，主体频带信息是否保留，大幅提高了波形聚类的可靠性。最后，将重构地震数据体转换到由原型向量组成的二维潜在空间网格中学习；将学习之后的原型向量利用k均值聚类进行聚类分组，得到聚类分布图。

本发明还提出了一种地震波形聚类装置，包括处理器，该处理器用于执行实现以下方法的指令：

1)获取目标层地震数据中地震道信号的频带范围，将频带范围划分，至少包括第一设定频带范围和第二设定频带范围；其中，第一设定频带范围为地震道信号的主体频带范围，主体频带范围为振幅能量强度降至最大值的一半时的频带范围；第二设定频带范围为除去第一设定频带范围外、所述目标层地震数据在剖面上保持层状展布特征的频带范围；

2)将所述目标层地震数据分解，生成至少两个固有模态函数分量，包括第一固有模态分量和第二固有模态分量，其中，第一固有模态分量的频带范围为所述第一设定频带范围，第二固有模态分量的频带范围为所述第二设定频带范围；

3)对第二固有模态分量进行去噪处理，处理后与第一固有模态分量叠加重构，得到重构地震数据体，并对该重构地震数据体进行学习和聚类，得到波形地震相图。

上述所指的地震波形聚类装置，实际上是基于本发明方法流程的一种计算机解决方案，即一种软件构架，可以应用到处理器中，上述装置即为与方法流程相对应的处理进程。由于对上述方法的介绍已经足够清楚完整，故不再详细进行描述。