CN105243388A - 基于动态时间规整和划分算法的波形分类方法 - Google Patents

Abstract

本发明公开一种基于动态时间规整和划分算法的波形分类方法，通过确定地震数据样本集；再根据地震数据样本集的地震相种类，确定地震数据样本集的分类数目k；从地震数据样本集中选择k个样本作为初始质心；然后基于DTW距离，将未被选为初始质心的样本数据分配到对应质心所在的簇；迭代地更新簇的质心；判断是否达到迭代次数上限，是则结束，得到最终分配的k个簇；否则根据步骤S5得到的簇的更新的质心继续执行步骤S4，进行重新分配。采用动态时间规整来对齐地震数据，减小了层位解释误差的影响，更精确地度量了两道地震数据之间的相似性；利用簇的质心来作为簇的中心，相较于传统划分算法中定义的簇的中心，更准确，更能代表一个簇的特征。

Description

基于动态时间规整和划分算法的波形分类方法

技术领域

本发明属于地震数据处理领域，具体涉及一种地震波形分类技术。

背景技术

地震波形是地震数据的基本性质，它包含了所有的定性和定量信息，如反射模式、相位、频率和振幅等信息，是地震信息的总体特征，其动态变化蕴含了丰富的内在信息，能够真实地反映地下结构的特征。波形分类法是最常采用的地震相分析方法，通过对地震信号波形进行分类，可以实现对地震相的划分。波形分类针对含有各种波形的地震数据样本集，通过合适的分类或者聚类方法，将样本划分到不同的类别，以达到区分波形样本的目的。

波形分类技术分为聚类分析和统计分类。聚类分析是把从图像或者信号中提取到的输入特征属性以“聚类”的形式分开，最后达到区分这些图像或者信号的目的。聚类分析是无监督的，不需要使用标签信息，没有样本数据的先验知识，挖掘的是样本之间的内在联系。统计分类主要是基于概率统计模型的方法，通过现实有标签的样本集得到各种特征属性和内在规律作为后续的分类依据，然后再对样本数据进行分类和识别，属于有监督的学习过程。

目前的波形分类技术主要基于无监督分类算法。聚类分析不需要先验信息，只需要输入工区的振幅数据，就能进行波形分类，对地震相进行划分。聚类分析的策略，是让同一类的成员间相似性最大，同时不同类的成员间相似性最小。通常采用的相似性度量方法是计算各个波形之间的“距离”。距离计算方法不同，各个样本之间的相似度就不同，最终的划分结果也就随之变化。因此，距离的计算方法在聚类分析中占有举足轻重的地位。

在波形分类中，相似性度量通常使用欧氏距离。欧氏距离是最易于理解的一种距离计算方法，源自欧氏空间中两点间的距离公式，也是使用最广泛的相似性度量方法。但是，欧氏距离存在着很明显的缺陷，当层位不准确且未校准层位时，在沿层位截取数据之后，同相轴没有对齐，此时如果用欧式距离计算相似性，会使原本比较相似的两道数据距离增大，即求得的相似性比实际的相似性要小很多，从而在后续的聚类过程中，可能把原本应该聚到同一个簇的两个样本，分到了两个不同的簇，从而出现错误的结果。欧氏距离对层位解释错误非常敏感，而层位解释错误在地震数据中是普遍存在的，尤其对于一些开发尚不成熟的工区。如果要使用欧氏距离来衡量波形的相似性，则需要另外的步骤来校正层位。

无监督的波形分类即为聚类，聚类分析最简单、最基本的版本是划分。基于划分的聚类算法是在给定n道地震数据集D以及要区分的地震相种类数k的条件下，把地震数据组织成k(k≤n)个分区，其中，一个分区称为一个簇，每个簇代表一种地震相。这些簇的形成旨在优化一个客观划分准则，如基于距离的相异性函数，使得根据数据集的属性，在同一个簇中的对象是“相似的”，而不同簇中的对象是“相异的”。最著名、最常用的划分方法有：k-means和k-medoids，这些名称是根据选定的簇中心的方法不同来区别的。k-means算法中，每个簇的中心都用簇中所有对象的均值来表示；而k-medoids是一种基于中心点或中心对象进行划分的算法。但现有的定义族的中心的方法，不能很好地代表一个族的特征。

发明内容

本发明所要解决的技术问题为：当存在层位解释错误时，如何准确地度量地震数据之间的相似性以及在基于划分的聚类算法中如何准确确定簇的中心；为解决所述技术问题，提出一种基于动态时间规整和划分算法的波形分类方法，采用动态时间规整(DTW,DynamicTimeWrapping)来求两道地震数据之间的距离，由于使用DTW距离来度量两道数据之间的相似性，因此在形成簇之后，采用了相应的簇中心更新算法。

本发明采用的技术方案是：基于动态时间规整和划分算法的波形分类方法，包括以下步骤：

S1：确定地震数据样本集；

S2：根据地震数据样本集的地震相种类，确定地震数据样本集的分类数目k；

S3：从地震数据样本集中选择k个样本作为初始质心；

S4：基于DTW距离，将未被选为初始质心的样本数据分配到对应质心所在的簇；

S5：迭代地更新簇的质心；

S6：判断是否达到迭代次数上限，是则结束，得到最终分配的k个簇；否则根据步骤S5得到的簇的更新的质心执行步骤S4，对未被选为初始质心的样本数据进行重新分配。

进一步地，步骤S1所述地震数据样本集为：叠前地震数据样本集或叠后地震数据样本集。

更进一步地，所述叠前地震数据样本集通过将叠前地震数据进行预处理，得到叠前地震数据样本集，具体过程为：由每一道叠前地震数据包括m个不同的偏移距，每个偏移距下有n个采样点；初始化若干值为0的列向量，且将每一列向量分成m^、个部分，每部分长度都为n^、，将每一道原始地震数据的偏移距按照从小到大进行排序，依次将同一偏移距的数据存入列向量对应序号的部分；最后得到叠前地震数据样本集；

其中，m≦m^、，n≦n^、。

更进一步地，对叠前地震数据样本集，基于DTW距离，求解最短距离时，必须满足：起点为(1,1)，终点为(m*n,m*n)；且，必经过点(j*n,j*n)；

其中，j为m个不同偏移距的序号，j＝1,2,……,(m-1)。

进一步地，步骤S2中所述地震数据样本集的地震相种类，由历史数据资料以及经验值，通过统计分析得到。

进一步地，所述步骤S3具体为：

S31：从输入的地震数据样本集中随机选择一道样本数据作为初始质心，并计初始质心个数p＝1；

S32：对于样本集中，未被选为初始质心的每一道样本数据，通过DTW距离，计算出在初始质心集合中与该道样本数据距离最小的初始质心，并记录下对应距离d(i)；

其中，i表示样本集中未被选为初始质心的样本数据的序号；

S33：将由步骤S32得到的所有d(i)相加得到SUM，并生成一个属于[0,SUM]的随机数Random；

S34：执行Random＝Random-d(i)；

S35：当Random≤0，此时的i所对应的样本数据为新的初始质心，则p＝p+1，并继续步骤S36；否则转至步骤S34；

S36：当p＝k时，得到k个初始质心，否则转至步骤S32。

进一步地，所述步骤S4具体为：根据由步骤S3选出的k个初始质心，基于DTW距离，计算得到每一个未被选为初始质心的样本数据与k个初始质心的最小距离，将样本数据分配到最小距离对应初始质心所在的簇，得到k个新的簇。

进一步地，所述步骤S5具体为：

S51：遍历簇中每一个样本数据，与该簇质心对齐，得到该簇中各样本数据与该簇质心对齐的点对；

S52：根据步骤S51中得到的点对，得到该簇各样本数据与该簇质心对其的的对齐点，并将对齐点振幅相加，取平均值得到该簇新的质心；

S53：判断是否达到第二迭代次数上限，是则结束，得到该簇更新后的质心；否则，以步骤S52得到的该簇新的质心，执行步骤S51至步骤S52。

本发明的有益效果：本发明的基于动态时间规整和划分算法的波形分类方法，针对存在层位解释误差的数据，采用动态时间规整来对齐地震数据，减小了层位解释误差的影响，更精确地度量了两道地震数据之间的相似性；本发明基于划分算法，利用簇的质心来作为簇的中心，相较于传统划分算法中定义的簇的中心，更准确，更能代表一个簇的特征；本发明将DTW距离和划分算法相结合，形成一个完整的波形分类流程，能更好地划分地震相。

附图说明

图1为本发明提供的序列对齐示意图。

图2为本发明提供的方案流程图。

具体实施方式

为便于本领域技术人员理解本发明的技术内容，下面结合附图对本发明内容进一步阐释。

动态时间规整(DTW,DynamicTimeWrapping)是一种衡量两个时间序列的相似度的方法，与欧氏距离不同的是，它不仅可以比较两个等长的时间序列之间的相似性，对于长度不同的时间序列也能比较相似性，同时还能消除序列之间的相位作用。

在波形分类中，根据层位加窗取得的数据是一道波形，由于存在层位解释错误、按层位截取出来的三维数据，每一道的第一个点与实际的层位剖面距离并不相同，存在一定的相位误差。为此，采用动态时间规整的方法来度量两道地震波形的相似性，从而消除这种相位误差。我们以一道数据为基准，取需要与该道计算相似性的一道作为对比道，遍历基准道中的每一个点，按照一定的规则，在对比道中寻找与基准道中的点最相似的点。这个过程，称为对齐。这样不仅可以校正层位，消除相位误差，把与实际层位剖面距离相同的点校正到同一个剖面，同时，也是一个计算相似性的过程。如图1所示，是两个波形C、Q对齐的例子，图A是需要比较相似性的两个波形，图B是将C、Q对齐的过程，图C是将两个序列对齐的结果显示。具体的DTW距离算法为现有的技术，因此在此处不作详细描述。

如图2所示为本发明方案的流程图，本发明的基于动态时间规整和划分算法的波形分类方法具体包括：

S1：确定地震数据样本集；

S2：根据地震数据样本集的地震相种类，确定地震数据样本集的分类数目k；

S3：从地震数据样本集中选择k个样本作为初始质心；

S4：基于DTW距离，将未被选为初始质心的样本数据分配到对应质心所在的簇；

S5：迭代地更新簇的质心；

S6：判断是否达到迭代次数上限，是则结束，得到最终分配的k个簇；否则根据步骤S5得到的簇的更新的质心执行步骤S4，对未被选为初始质心的样本数据进行重新分配。

步骤S1所述的地震数据样本集为：叠前地震数据样本集或叠后地震数据样本集；叠后地震数据即为样本集格式名不需进行数据预处理，直接进行后续处理。

所述叠前地震数据样本集通过将叠前地震数据进行预处理，得到叠前地震数据样本集，具体过程为：由每一道原始地震数据包括m个不同的偏移距，每个偏移距下有n个采样点；初始化若干值为0的列向量，且将每一列向量分成m^、个部分，每部分长度都为n^、，将每一道原始地震数据的偏移距按照从小到大进行排序，依次将同一偏移距的数据存入列向量对应序号的部分；最后得到地震数据样本集；

其中，m≦m^、，n≦n^、。

例如，对于第i道叠前地震数据，假设有m个不同的偏移距，每个偏移距下有n个采样点；则对第i道叠前地震数据的预处理过程包括以下步骤：

S11：初始化一个值为0，长度为m^、*n^、的列向量i^、；

S12：将列向量i^、分成m^、个部分，每部分长度都为n^、，每一部分用来存放偏移距相同的数据，即从第j^、*n^、+1到第j^、*(n^、+1)个位置存放偏移距相同的数据；

其中，j^、＝0,1,2,……,m^、-1；

S13：将偏移距按照从小到大进行排序，将第i道叠前地震数据偏移距排序为第j+1的数据存入列向量i^、的第j^、+1部分；

其中，j＝0,1,2,……,m-1；

S14：重复步骤S11至步骤S13，得到叠前地震数据样本集。

步骤S2具体为：根据查询现有的技术资料，以及相关技术人员的经验数据，通过统计分析得到原始地震数据的地震相的种类数，根据得到的地震相的种类数设置地震数据样本集的分类数目。

步骤S3具体为：为了使不同聚类中的成员之间的距离最大，选取初始质心时，采取的策略是：使它们之间的距离尽可能远。本申请选用DTW距离。

首先从输入的地震数据样本集中随机地选择一道作为初始质心；然后对于数据集中，未被选为初始质心的每一个样本数据，找出各个样本数据与初始质心距离最小的质心并记录下距离d(i)；将所有的d(i)相加得到SUM，再生成一个属于[0,SUM]的随机数Random，然后用Random＝Random-d(i)，直到Random≤0，此时的样本就是下一个初始质心，直到选定所有的k个初始质心。详细步骤如下：

S31：从输入的地震数据样本集中随机选择一道样本数据作为初始质心，并计初始质心个数p＝1；

S32：对于样本集中，未被选为初始质心的每一道样本数据，通过DTW距离，计算出在初始质心集合中与该道样本数据距离最小的初始质心，并记录下对应距离d(i)；

其中，i表示样本集中未被选为初始质心的样本数据的序号；

S33：将由步骤S32得到的所有d(i)相加得到SUM，并生成一个属于[0,SUM]的随机数Random；

S34：执行Random＝Random-d(i)；

S35：当Random≤0，此时的i所对应的样本数据为新的初始质心，则p＝p+1，并继续步骤S36；否则转至步骤S34；

S36：当p＝k时，得到k个初始质心，否则转至步骤S32。

步骤S4具体为：根据由步骤S3选出的k个初始质心，基于DTW距离，计算得到每一个未被选为初始质心的样本数据与k个初始质心的最小距离，将样本数据分配到最小距离对应初始质心所在的簇，得到k个新的簇。在采用DTW距离，计算叠前地震数据样本集的样本数据，求解最短路径需同时满足以下两个条件：

1)起点为(1,1)，终点为(m*n,m*n)；

2)必经过点(j*n,j*n)；

其中，j为m个不同偏移距的序号，j＝1,2,……,(m-1)。

步骤S5具体为：遍历簇中每一个样本数据，得到与原始质心对齐的点，将得到的各点的振幅相加，取平均值得到该簇新的质心。本发明中，选定簇的质心为簇的中心。对于一个包含若干样本的簇，更新其质心的算法如下：S51：遍历簇中每一个样本数据，与该簇质心对齐，得到该簇中各样本数据与该簇质心对齐的点对，并记录对齐的点对情况；

S52：根据步骤S51中得到的点对，得到该簇各样本数据与该簇质心对其的的对齐点，并将对齐点振幅相加，取平均值得到该簇新的质心；

通过遍历簇中每一道样本数据，取出与原始质心那一道样本数据第w个点对齐的若干点，把这些点的振幅相加，取平均值即可，得到新的质心的第w个点，从而得到该簇的更新的质心。

S53：判断是否达到第二迭代次数上限，是则结束，得到该簇更新后的质心；否则，以步骤S52得到的该簇新的质心，执行步骤S51至步骤S52。所述的第二迭代次数上限，本申请中取值为30.

本申请判断迭代停止的依据为：当各个簇中的样本在前后两次迭代中不再变化或者是达到迭代上限，则停止迭代；通过实验验证，在不断迭代的过程中，簇中的样本数据变化会越来越小，因此，本申请以迭代次数的上限作为迭代停止的判断依据。在本领域中通常将迭代停止的上限设置在10的附近，本申请取迭代上限为10。

得到了一个簇的质心后，需要根据新的质心把样本聚成新的簇，再更新质心，如此循环，直到聚成的簇不再变化或者迭代次数达到上限为止，得到最终分配的k个簇。

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

Claims (8)

1.基于动态时间规整和划分算法的波形分类方法，其特征在于，包括以下步骤：

S1：确定地震数据样本集；

S2：根据地震数据样本集的地震相种类，确定地震数据样本集的分类数目k；

S3：从地震数据样本集中选择k个样本作为初始质心；

S4：基于DTW距离，将未被选为初始质心的样本数据分配到对应质心所在的簇；

S5：迭代地更新簇的质心；

S6：判断是否达到第一迭代次数上限，是则结束，得到最终分配的k个簇；否则根据步骤S5得到的簇的更新的质心执行步骤S4，对未被选为初始质心的样本数据进行重新分配。

2.根据权利要求1所述的基于动态时间规整和划分算法的波形分类方法，其特征在于，步骤S1所述地震数据样本集为：叠前地震数据样本集或叠后地震数据样本集。

3.根据权利要求2所述的基于动态时间规整和划分算法的波形分类方法，其特征在于，所述叠前地震数据样本集通过将叠前地震数据进行预处理，得到叠前地震数据样本集，具体过程为：由每一道叠前地震数据包括m个不同的偏移距，每个偏移距下有n个采样点；初始化若干值为0的列向量，且将每一列向量分成m`个部分，每部分长度都为n`，将每一道原始地震数据的偏移距按照从小到大进行排序，依次将同一偏移距的数据存入列向量对应序号的部分；最后得到叠前地震数据样本集；

其中，m≦m`，n≦n`。

4.根据权利要求3所述的基于动态时间规整和划分算法的波形分类方法，其特征在于，对叠前地震数据样本集，基于DTW距离，求解最短距离时，必须满足：起点为(1,1)，终点为(m*n,m*n)；且，必经过点(j*n,j*n)；

其中，j为m个不同偏移距的序号，j＝1,2,……,(m-1)。

5.根据权利要求1所述的基于动态时间规整和划分算法的波形分类方法，其特征在于，步骤S2中所述地震数据样本集的地震相种类，由历史数据资料以及经验值，通过统计分析得到。

6.根据权利要求1所述的基于动态时间规整和划分算法的波形分类方法，其特征在于，所述步骤S3具体为：

S31：从输入的地震数据样本集中随机选择一道样本数据作为初始质心，并计初始质心个数p＝1；

S32：对于样本集中，未被选为初始质心的每一道样本数据，通过DTW距离，计算出在初始质心集合中与该道样本数据距离最小的初始质心，并记录下对应距离d(i)；

其中，i表示样本集中未被选为初始质心的样本数据的序号；

S33：将由步骤S32得到的所有d(i)相加得到SUM，并生成一个属于[0,SUM]的随机数Random；

S34：执行Random＝Random-d(i)；

S35：当Random≤0，此时的i所对应的样本数据为新的初始质心，则p＝p+1，并继续步骤S36；否则转至步骤S34；

S36：当p＝k时，得到k个初始质心，否则转至步骤S32。

7.根据权利要求1所述的基于动态时间规整和划分算法的波形分类方法，其特征在于，所述步骤S4具体为：根据由步骤S3选出的k个初始质心，基于DTW距离，计算得到每一个未被选为初始质心的样本数据与k个初始质心的最小距离，将样本数据分配到最小距离对应初始质心所在的簇，得到k个新的簇。

8.根据权利要求1所述的基于动态时间规整和划分算法的波形分类方法，其特征在于，所述步骤S5具体为：

S51：遍历簇中每一个样本数据，与该簇质心对齐，得到该簇中各样本数据与该簇质心对齐的点对；

S52：根据步骤S51中得到的点对，得到该簇各样本数据与该簇质心对其的的对齐点，并将对齐点振幅相加，取平均值得到该簇新的质心；

S53：判断是否达到第二迭代次数上限，是则结束，得到该簇更新后的质心；否则，以步骤S52得到的该簇新的质心，执行步骤S51至步骤S52。