CN104216015A - 基于隐含狄利克雷分布的三维地震信号分类方法 - Google Patents

Abstract

本发明公开了一种基于隐含狄利克雷分布的三维地震信号分类方法，包括以下步骤：S1：输入原始三维地震数据；S2：数据预处理：根据所要分析的层位，提取出所要进行波形分类分析的目的层段数据；S3：利用切比雪夫多项式拟合进行属性提取：对每一道地震数据进行切比雪夫多项式拟合，用得到的拟合系数来表示原来那一道地震道数据，经过拟合后，得到了一个三维系数体；S4：对步骤S3得到的三维系数体进行LDA分类和分类标示。本发明主要运用切比雪夫多项式拟合、期望最大算法(EM)和隐含狄利克雷分布算法(LDA)，在无人工干预的情况下，实现了三维地震信号的快速分类。

Description

基于隐含狄利克雷分布的三维地震信号分类方法

技术领域

本发明涉及地震信号分析和分类，特别涉及一种基于隐含狄利克雷分布的三维地震信号分类方法。

背景技术

在油气的地震勘探中，地震资料解释的目的是为了从地震数据中提取更多的信息进行地下构造解释以及地层和岩性特征的描述。从地震资料中获取这些信息最有效的方法之一就是地震属性特征的提取分析及波形分类技术。而基于三维地震信号的波形分类问题的难度由于地层环境的复杂度而相当大。

在地震信号的波形分类问题上，目前国内外有较多的实现方案。主要有两大类，一类是无监督分类，一类是有监督分类：

(1)无监督分类如SOM，要有很好的初始化条件；

(2)有监督分类如SVM，花费大量的内存，和需要很好的带标签样本；

以上现有技术，在实际的地震资料分析中都得到了很大程度的应用，各自也存在很多缺陷，主要的缺陷有算法复杂度太大，耗时间又需要非常大的内存，还需要非常好的初始化条件，这些缺陷影响这些方法在实际应用。

油气地震勘探面临着复杂地表和复杂地质构造的难题，在这些地区，地震波场复杂，地质构造变化剧烈，使得有效地识别弱信号，消除噪声干扰并提高地震资料信噪比成为勘探难点；在老油田，易于寻找的油气藏越来越少，取而代之的是以岩性地层为主、一般技术手段难以发现的隐蔽性特殊储层油气藏，采用常规的构造勘探方法很难取得突破，需要更好更加精细的勘探方法和技术，获取油藏地震属性信息，以进一步刻画老的油气藏和发现新的油气藏，地震信号分析方法和技术是一条重要的途径。

发明内容

本发明的目的在于克服现有技术的不足，提供一种提出了一种基于隐含狄利克雷分布的波形分类方法，主要运用切比雪夫多项式拟合、期望最大算法(EM)和隐含狄利克雷分布算法(LDA)，在无人工干预的情况下，实现三维地震信号分类的基于隐含狄利克雷分布的三维地震信号分类方法。

本发明的目的是通过以下技术方案来实现的：基于隐含狄利克雷分布的三维地震信号分类方法，包括以下步骤：

S1：输入原始三维地震数据；

S2：数据预处理：根据所要分析的层位，提取出所要进行波形分类分析的目的层段数据；

S3：利用切比雪夫多项式拟合进行属性提取：对每一道地震数据进行切比雪夫多项式拟合，用得到的拟合系数来表示原来那一道地震道数据，经过拟合后，得到了一个三维系数体；

S4：对步骤S3得到的三维系数体进行LDA分类和分类标示。

进一步地，所述的步骤S2数据预处理具体方法为：

S21：判断原始数据是否为两层位数据，若是则分析两层位间，确定两层位，并提取沿层间数据，得到目标层段数据；否则分析单层位，确定时窗及分析层位，并提取沿层时窗数据，得到目标层段数据；

S22：判断目标层段数据是否为正规的四方体，若不是则利用三次样条插值对得到的目标层段数据进行处理，使得每一道数据的维度相等，三维地震数据成为一个规整的四方体数据。

进一步地，所述的步骤S3中利用切比雪夫多项式拟合包括以下子步骤：

S31：指定切比雪夫多项式拟合阶数N；

S32：选取N+1个切比雪夫点：

$x_k=\cos \mathrm{\π}\frac{2K+1}{2N+2},k=\mathrm{0,1,2},...,N$

S33：计算切比雪夫点处的取值f(x_k)；

S34：计算切比雪夫系数：

$c0=\frac{1}{N+1}\underset{k=0}{\overset{N}{\mathrm{\Σ}}}f\left(x_k\right)T_0\left(x_k\right)=\frac{1}{N+1}\underset{k=0}{\overset{N}{\mathrm{\Σ}}}f\left(x_k\right)$

$\mathrm{cj}=\frac{2}{N+1}{\mathrm{\Σ}}_{k=0}^{N}f\left(x_k\right)T_j\left(x_k\right)=\frac{2}{N+1}{\mathrm{\Σ}}_{k=0}^{N}f\left(x_k\right)\cos \frac{\mathrm{j\π}(2k+1)}{2N+2}.$

进一步地，所述的步骤S4包括以下步骤：

S41：建立LDA模型：假设整个三维地震数据是由K个类模型产生的，一道地震数据是由K个类模型中的某一个类模型生成的，这些类模型服从参数为θ的多项式分布，每个类模型又与地震数据中的V个数据的一个多项分布相对应，将这个分布记为对于地震数据的任意一道数据，LDA定义了如下的生成过程：

S411：对每一道数据，从类模型分布中抽出一个类模型；

S412：从上述被抽到的类模型对对应的数据分布中抽取处一个数据；

S413：重复步骤S411和步骤S412，直至遍历地震道中的每一个数据；

S42：计算估计类模型参数θ和每一个类模型的参数具体包括以下子步骤：

S421：设计最大迭代次数和判断参数稳定的条件，用于在训练参数的时候确定训练什么时候终止或者参数是否稳定；

S422：初始化K个类模型的概率分布，因为这些模型服从参数为θ的多项式分布，所以开始先设每个类模型出现的概率相等，均为1/；

S423：初始化每一个类模型的概率分布，因为每个类模型的概率分布服从参数为的多项式分布，开始先假设对于任意一个类模型生成每一个数据的概率都是相等的，先预处理数据，使得所有数据的范围是在1到1000的某一个整数，所有刚开始每一个类模型生成每一个数据的概率为0.001，为了让模型能进行下去，初始化每个数据的生成概率为0.001加上一个随机数，加上一个随机数是为了让每个类模型有些差别，这样刚开始训练数据的时候，对于一道数据每个类模型的概率不相同，方便选择最优结果；

S424：训练模型：把所有的地震道数据输入模型中，并且记录每一道地震数据在该模型的现有的参数下生成概率。

S425：根据所有地震数据在现有模型下的概率，更新参数θ，参数，更新的准则是根据极大似然估计的原理的，更新参数使得模型生成所有地震道数据的概率最大；

S426：重复上述步骤S423～S425，直到模型的参数稳定或者迭代的次数超过我们预先设定的大小，得到最后训练的参数θ，训练出LDA的模型；

S43：分类标示：根据得到的LDA模型，把每一道数据重新再次输入模型中，假设第i道数据输入了模型，比较这一道数据在各个类模型的生成概率，若这一道数据是由第j类模型的生成概率最大，就把第i道数据的分类标签设为j，遍历所有的地震道数据，最后得出了所有地震道的分类标签，实现三维地震信号的分类。

具体地，所述的步骤S425中更新参数θ，参数的具体方法为：假设数据大小1到1000，第i道数据由第K个类模型生成的，表示第K个类模型的概率分布第i道数据为{n₁,n₂,n₃,n₄,n₅……n_m}，其中，m为每一道数据的维度，更新对于的系数每个系数加上一个小步长，然后再对做处理，使得里面的系数都在0和1之间，并且求和为1；然后再更新θ，其中θ＝{θ₁，θ₂，θ₃，……，θ_n}，其中，n为类模型个数，因为是第K个模型生成的，因此让θ_i加上一个步长，然后再对θ处理，使得所有θ₁+θ₂+θ₃+……+θ_n＝1，且任意一个θ_i大于0。

本发明的有益效果是：

1、提出了一种基于隐含狄利克雷分布的波形分类方法，主要运用切比雪夫多项式拟合、期望最大算法(EM)和隐含狄利克雷分布算法(LDA)：将地震信号做切比雪夫多项式拟合，这种拟合去地震数据冗余也无疏漏，拟合后的系数能完好的表征原始地震数据；在无需人工干预情况下，使用EM算法对得到的特征参数进行优化选择并聚类；定义了针对三维地震数据的LDA分类模型，该模型能够很好的对三维地震信号进行分类；在无人工干预的情况下，实现了三维地震信号的分类。

2、采用本发明的算法能够快速得到三维地震数据的分类，在一定程度上能和传统的分类算法起到互补作用。

附图说明

图1为EM算法流程图；

图2为本发明的LDA工作流程图；

图3为本发明的分类算法流程图；

图4为本发明的提取目的层段数据流程图；

图5为本实施例的LDA生成模型图；

图6为本实施例的训练LDA模型流程图。

具体实施方式

下面对本发明所运用的一些用语进行解释说明：

1、期望最大算法(EM)

EM描述如下：已知一个概率模型包括隐变量集Z、观测集X和参数集S，目标为了得到p(X|S)最大化的S。

最大期望算法经过两个步骤交替进行计算：

E步骤：计算期望(E)，利用对隐藏变量的现有估计值，计算其最大似然估计值；

M步骤：最大化(M)，最大化在E步上求得的最大似然值来计算参数的值。

M步上找到的参数估计值被用于下一个E步计算中，这个过程不断交替进行，流程图如图1所示，总体来说，EM的算法流程如下：

(1)初始化分布参数(用户可以根据自己的需求选择初始化的方法，可以选择随机初始化，也可以根据输入的数据来初始化)；

(2)重复步骤(3)、(4)的E、M步骤直到收敛或者到达我们迭代的最大次数；

(3)最大期望算法中的E步骤：估计未知参数的期望值，给出当前的参数估计；

(4)最大期望算法中的M步骤：重新估计分布参数，以使得数据的似然性最大，给出未知变量的期望估计。

2、隐含狄利克雷(LDA)分布

LDA的工作原理如图2所示，其中的参数原理如下：

α为类模型概率的概率分布、Dirichlet参数，θd为地震道d下的类模型概率分布，Z_d,n为第n个数据在地震道d产生的类模型，W_d,n为类模型产生的数据，β为类模型与数据概率分布，η为使β满足Dirichlet分布的参数。

下面结合附图和具体实施例进一步说明本发明的技术方案，但本发明所保护的内容不局限于以下所述。

如图3所示，基于隐含狄利克雷分布的三维地震信号分类方法，包括以下步骤：

S1：输入原始三维地震数据；

S2：数据预处理：根据所要分析的层位，提取出所要进行波形分类分析的目的层段数据；

S3：利用切比雪夫多项式拟合进行属性提取：对每一道地震数据进行切比雪夫多项式拟合，用得到的拟合系数来表示原来那一道地震道数据，经过拟合后，得到了一个三维系数体；

S4：对步骤S3得到的三维系数体进行LDA分类和分类标示。

如图4所示，本实施例所述的步骤S2数据预处理具体方法为：

S21：判断原始数据是否为两层位数据，若是则分析两层位间，确定两层位，并提取沿层间数据，得到目标层段数据；否则分析单层位，确定时窗及分析层位，并提取沿层时窗数据，得到目标层段数据；

S22：判断目标层段数据是否为正规的四方体，若不是则利用三次样条插值对得到的目标层段数据进行处理，使得每一道数据的维度相等，三维地震数据成为一个规整的四方体数据。

其中，各个名词的定义为：

层位：指在地层层序中的某一特定位置，地层的层位可以是地层单位的界线，也可以是属于某一特定时代的标志层等。

断层：地壳岩层因受力达到一定强度而发生破裂，并沿破裂面有明显相对移动的构造称为断层。

插值：利用已知点来计算未知点的过程。

目的层段：是指在三维地震信号数据中所要分析的沿单个层位的时窗数据或所要分析的两个层位间的层间数据。

属性提取的主要目的是为了优化目的层段数据中的特征，这样不仅会降低计算内存的占用率和后面分类模型设计的复杂度，提供效率，也能够在一定程度上提高分类精度，减少非相关特征对分类结果的影响，本实施例选择的属性提取算法为切比雪夫多项式拟合算法。所述的步骤S3中利用切比雪夫多项式拟合包括以下子步骤：

S31：指定切比雪夫多项式拟合阶数N；

S32：选取N+1个切比雪夫点：

$x_k=\cos \mathrm{\π}\frac{2K+1}{2N+2},k=\mathrm{0,1,2},...,N$

S33：计算切比雪夫点处的取值f(x_k)；

S34：计算切比雪夫系数：

$c0=\frac{1}{N+1}\underset{k=0}{\overset{N}{\mathrm{\Σ}}}f\left(x_k\right)T_0\left(x_k\right)=\frac{1}{N+1}\underset{k=0}{\overset{N}{\mathrm{\Σ}}}f\left(x_k\right)$

$\mathrm{cj}=\frac{2}{N+1}{\mathrm{\Σ}}_{k=0}^{N}f\left(x_k\right)T_j\left(x_k\right)=\frac{2}{N+1}{\mathrm{\Σ}}_{k=0}^{N}f\left(x_k\right)\cos \frac{\mathrm{j\π}(2k+1)}{2N+2}$

在这里，我们按照上述步骤对-1<x<1分别进行5阶、8阶、10阶多项式拟合。切比雪夫多项式能够很好拟合原始波形，因此切比雪夫多项式拟合是一种有效的特征提取方法。虽然目前存在很多的特征提取方法，但是一个合适的特征提取方法应该是既能很好保留原始信息，又能有效降低信号冗余。符合这两个准则的特征提取方法分类结果总体来说差别不会太大，仅会存在一些细微差别。

应用到实际数据时，对每一道地震数据进行10阶的切比雪夫多项式拟合，这样就得到了十个多项式系数c1,c2……c10，用这十个系数代表原来的一道地震数据，最后得到一个以多项式系数为元素的三维体。

本实施例所述的步骤S4包括以下步骤：

S41：建立LDA模型：假设整个三维地震数据是由K个类模型产生的，一道地震数据是由K个类模型中的某一个类模型生成的，这些类模型服从参数为θ的多项式分布，每个类模型又与地震数据中的V个数据的一个多项分布相对应，将这个分布记为，聪的生成模型如图5所示，对于地震数据的任意一道数据，LDA定义了如下的生成过程：

S411：对每一道数据，从类模型分布(类模型的概率分布是按参数为θ的多项式分布)中抽出一个类模型；

S412：从上述被抽到的类模型对对应的数据分布中抽取处一个数据；

S413：重复步骤S411和步骤S412，直至遍历地震道中的每一个数据；

也就是说，对于任意一道地震数据D中的每一个数据，从该地震数据所对应的多项分布θ中抽取一个类模型Z，然后再从类模型z所对应的多项分布中抽取一个地震数据w，将这个过程重复N次，就产生了地震道数据D。这个过程就有两个参数需要推断。一个是地震道数据与类模型的分布θ；一个是每一个类模型-地震数据的分布

S42：计算估计类模型参数θ和每一个类模型的参数估计参数的算法主要是期望最大算法(EM)，得到整个地震数据的LDA生成模型，再用极大似然估计来判断每一道数据是由哪一个类模型生成的，这样就实现了的分类；如图6所示，参数估计具体包括以下子步骤：

S421：设计最大迭代次数和判断参数稳定的条件，用于在训练参数的时候确定训练什么时候终止或者参数是否稳定；

S422：初始化K个类模型的概率分布，因为这些模型服从参数为θ的多项式分布，所以开始先设每个类模型出现的概率相等，均为1/；

S423：初始化每一个类模型的概率分布，因为每个类模型的概率分布服从参数为的多项式分布，开始先假设对于任意一个类模型生成每一个数据的概率都是相等的，先预处理数据，使得所有数据的范围是在1到1000的某一个整数，所有刚开始每一个类模型生成每一个数据的概率为0.001，为了让模型能进行下去，初始化每个数据的生成概率为0.001加上一个随机数，加上一个随机数是为了让每个类模型有些差别，这样刚开始训练数据的时候，对于一道数据每个类模型的概率不相同，方便选择最优结果；

S424：训练模型：把所有的地震道数据输入模型中，并且记录每一道地震数据在该模型的现有的参数下生成概率，这个步骤就是期望最大算法(EM)中的E步；

S425：根据所有地震数据在现有模型下的概率，更新参数θ，参数，更新的准则是根据极大似然估计的原理的，更新参数使得模型生成所有地震道数据的概率最大，这个步骤就是EM算法中的M步；

S426：重复上述步骤S423～S425，直到模型的参数稳定或者迭代的次数超过我们预先设定的大小，得到最后训练的参数θ，训练出LDA的模型；

S43：分类标示：根据得到的LDA模型，把每一道数据重新再次输入模型中，假设第i道数据输入了模型，比较这一道数据在各个类模型的生成概率，若这一道数据是由第j类模型的生成概率最大，就把第i道数据的分类标签设为j，遍历所有的地震道数据，最后得出了所有地震道的分类标签，实现三维地震信号的分类。

本实施例所述的步骤S425中更新参数θ，参数的具体方法为：假设数据大小1到1000，第i道数据由第K个类模型生成的，表示第K个类模型的概率分布第i道数据为{n₁,n₂,n₃,n₄,n₅……n_m}，其中，m为每一道数据的维度，更新对于的系数每个系数加上一个小步长(步长的大小由用户设计，一般为1/10000)，然后再对做处理，使得里面的系数都在0和1之间，并且求和为1；然后再更新θ，其中θ＝{θ₁，θ₂，θ₃，……，θ_n}，其中，n为类模型个数，因为是第K个模型生成的，因此让θ_i加上一个步长(步长的大小也是由用户设计)，然后再对θ处理，使得所有θ₁+θ₂+θ₃+……+θ_n＝1，且任意一个θ_i大于0。

下面简述本发明在F3叠后地震信号分类的应用

(1)数据的介绍：F3叠后地震信号是荷兰北海油田一块地区采集得到的，该数据采集时间是1987年并向研究人员公开，F3数据现已成为该研究领域的常用示例数据，其原始数据可以在网上下载到。采集该数据的目的是通过分析该数据去发现侏罗纪与白垩纪地层之间的油气，研究人员最终的确在该地区发现了油气储藏。

(2)数据范围：波形分类采用LDA分类算法对F3叠后地震信号名为MSF4D的层位(注：因为地下结构大多呈层状叠加，层位可以简单理解为沿某一地层的二维切片，大小为593*943)，取上下区间为[-8,+24]ms的总共33个采样点的窗口长度的信号，这样就相当于取出了一个拥有593*943个样本，每个样本长度为33个采样点的三维地层。进而采用本发明的算法，得到分类结果。

本发明的算法能在很多的时间内得出分类结果，并且需要的计算机内存小，该处分类结果自然连续且成块，并且把当地的详细地层信息挖掘出来，所以该算法的效果明显。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (5)

1.基于隐含狄利克雷分布的三维地震信号分类方法，其特征在于，包括以下步骤：

S1：输入原始三维地震数据；

S2：数据预处理：根据所要分析的层位，提取出所要进行波形分类分析的目的层段数据；

S3：利用切比雪夫多项式拟合进行属性提取：对每一道地震数据进行切比雪夫多项式拟合，用得到的拟合系数来表示原来那一道地震道数据，经过拟合后，得到了一个三维系数体；

S4：对步骤S3得到的三维系数体进行LDA分类和分类标示。

2.根据权利要求1所述的三维地震信号分类方法，其特征在于，所述的步骤S2数据预处理具体方法为：

S21：判断原始数据是否为两层位数据，若是则分析两层位间，确定两层位，并提取沿层间数据，得到目标层段数据；否则分析单层位，确定时窗及分析层位，并提取沿层时窗数据，得到目标层段数据；

S22：判断目标层段数据是否为正规的四方体，若不是则利用三次样条插值对得到的目标层段数据进行处理，使得每一道数据的维度相等，三维地震数据成为一个规整的四方体数据。

3.根据权利要求1所述的三维地震信号分类方法，其特征在于，所述的步骤S3中利用切比雪夫多项式拟合包括以下子步骤：

S31：指定切比雪夫多项式拟合阶数N；

S32：选取N+1个切比雪夫点：

$x_k=\cos \mathrm{\&pi;}\frac{2K+1}{2N+2},k=\mathrm{0,1,2},...,N$

S33：计算切比雪夫点处的取值f(x_k)；

S34：计算切比雪夫系数：

$c0=\frac{1}{N+1}\underset{k=0}{\overset{N}{\mathrm{\&Sigma;}}}f\left(x_k\right)T_0\left(x_k\right)=\frac{1}{N+1}\underset{k=0}{\overset{N}{\mathrm{\&Sigma;}}}f\left(x_k\right)$

$\mathrm{cj}=\frac{2}{N+1}{\mathrm{\&Sigma;}}_{k=0}^{N}f\left(x_k\right)T_j\left(x_k\right)=\frac{2}{N+1}{\mathrm{\&Sigma;}}_{k=0}^{N}f\left(x_k\right)\cos \frac{\mathrm{j\&pi;}(2k+1)}{2N+2}.$

4.根据权利要求1所述的三维地震信号分类方法，其特征在于，所述的步骤S4包括以下步骤：

S41：建立LDA模型：假设整个三维地震数据是由K个类模型产生的，一道地震数据是由K个类模型中的某一个类模型生成的，这些类模型服从参数为θ的多项式分布，每个类模型又与地震数据中的V个数据的一个多项分布相对应，将这个分布记为对于地震数据的任意一道数据，LDA定义了如下的生成过程：

S411：对每一道数据，从类模型分布中抽出一个类模型；

S412：从上述被抽到的类模型对对应的数据分布中抽取处一个数据；

S413：重复步骤S411和步骤S412，直至遍历地震道中的每一个数据；

S42：计算估计类模型参数θ和每一个类模型的参数具体包括以下子步骤：

S421：设计最大迭代次数和判断参数稳定的条件，用于在训练参数的时候确定训练什么时候终止或者参数是否稳定；

S422：初始化K个类模型的概率分布，因为这些模型服从参数为θ的多项式分布，所以开始先设每个类模型出现的概率相等，均为1/k；

S423：初始化每一个类模型的概率分布，因为每个类模型的概率分布服从参数为的多项式分布，开始先假设对于任意一个类模型生成每一个数据的概率都是相等的，先预处理数据，使得所有数据的范围是在1到1000的某一个整数，所有刚开始每一个类模型生成每一个数据的概率为0.001，为了让模型能进行下去，初始化每个数据的生成概率为0.001加上一个随机数，加上一个随机数是为了让每个类模型有些差别，这样刚开始训练数据的时候，对于一道数据每个类模型的概率不相同，方便选择最优结果；

S424：训练模型：把所有的地震道数据输入模型中，并且记录每一道地震数据在该模型的现有的参数下生成概率；

S425：根据所有地震数据在现有模型下的概率，更新参数θ，参数，更新的准则是根据极大似然估计的原理的，更新参数使得模型生成所有地震道数据的概率最大；

S426：重复上述步骤S423～S425，直到模型的参数稳定或者迭代的次数超过我们预先设定的大小，得到最后训练的参数θ，训练出LDA的模型；

S43：分类标示：根据得到的LDA模型，把每一道数据重新再次输入模型中，假设第i道数据输入了模型，比较这一道数据在各个类模型的生成概率，若这一道数据是由第j类模型的生成概率最大，就把第i道数据的分类标签设为j，遍历所有的地震道数据，最后得出了所有地震道的分类标签，实现三维地震信号的分类。

5.根据权利要求4所述的三维地震信号分类方法，其特征在于，所述的步骤S425中更新参数θ，参数的具体方法为：假设数据大小1到1000，第i道数据由第K个类模型生成的，表示第K个类模型的概率分布第i道数据为{n₁,n₂,n₃,n₄,n₅……n_m}，其中，m为每一道数据的维度，更新对于的系数每个系数加上一个小步长，然后再对做处理，使得系数都在0和1之间，并且求和为1；然后再更新θ，其中θ＝{θ₁，θ₂，θ₃，……，θ_n}，其中，n为类模型个数，因为是第K个模型生成的，因此让θ_i加上一个步长，然后再对θ处理，使得所有θ₁+θ₂+θ₃+……+θ_n＝1，且任意一个θ_i大于0。