CN104101902B - 地震属性聚类方法及装置 - Google Patents

Abstract

本发明公开了一种地震属性聚类方法及装置，其中地震属性聚类方法采用基于快速K‑均值的地震属性聚类，输入为预先聚类别个数及待聚类的地震属性数据，输出为每个地震属性数据所属类别标号，处理时根据三角不等式的原理省去了原始K‑均值方法每一次循环迭代都要计算每个地震属性数据到类别中心的距离用来更新每一个地震属性数据所属类别标号中的一部分不必要的计算。并且还可以包括在输入待聚类的地震属性数据之前，对每一个地震属性数据进行高斯归一化处理，以及，在进行高斯归一化处理之前，对地震属性数据中记录错误的原始时间序列信号异常数值进行剔除。本发明可以在短时间内完成大批量地震属性数据聚类分析，为勘探家提供进一步详细地质分析的基础。

Description

地震属性聚类方法及装置

技术领域

本发明涉及地球物理勘探技术领域，尤其涉及地震属性聚类方法及装置。

背景技术

在油气勘探中，只有对地下的地质情况有了充分的了解和熟悉后，才能对勘探区域的油气储藏状况做出判断。获取地质信息的一个重要手段就是分析地震数据经数学变换后得到的各种地震属性数据。地震属性数据通常是叠前或叠后地震数据，有关地震波的几何形态、运动学特征、动力学特征等参数。通过对这些参数的研究，可以得到勘探区域地下介质的结构、岩性、流体等的特征，进而推断油气的储藏信息。从这些获取的地震属性数据中经过一系列分析推断出地下地质状况的这一过程通常称为地震属性分析。其中最常用的一种方法就是聚类。所谓聚类就是根据地下介质处获取的地震属性数据之间差异的大小，将它们分为若干类别，每一类内的数据间相差较小而类别之间差别较大。通过对收集的地震属性进行聚类，可以将这些地震属性数据可以分为几大类别，进而可以对勘探区域的地质情况进行进一步的分析。例如对目标区域进行地质相带划分：根据聚类结果和测井解释结果的对照分析，来确定每个类别所对应的相带。尤其在储层预测的过程中，地震属性聚类分析是一个非常必要的步骤，起着比较重要的作用。

随着数据采集技术的不断发展和完善以及对地震属性的认识不断进步，可以得到的属性越来越多。随着勘探的区域范围以及目标地层深度不断增加，所获得的地震属性数据正在高速爆炸式地增长。如何对这些大批量的地震属性数据进行快速聚类是一个很重要的研究课题。通常这些数据的量级已经超出在一般计算机内存和处理器所能承受的范围。即使计算机的配置可以承受如此大量的数据，那么传统聚类方法运行速度也会很慢，大大超出了人们所能承受的时间范围，对使用地震属性数据进行地质分析如相带划分和储层预测都造成了很大的困难。

地震数据本身具有不同于其它的数据集的特点：首先地震属性数据通常数量会非常庞大，有时是二维类平(曲)面属性，也有时是从一个三维属性体中切出的体数据。无论二维还是三维范围的数据，一般其数量都会非常大。另一个特点地震属性数据的维数通常不会太高，即地震属性的个数会远远小于地震数据的个数。上面的举例中已经列举了三个最常用的地震信号属性了，虽然之前的文献显示目前可获取的地震属性达上百种，但通常情况下用于聚类的地震属性最多也只有几十种的量级，相对于地震数据的数量非常小。这个特性说明了对数量庞大但维数不高的地震数据进行聚类分析非常有必要。

很早之前就已经有人运用模式识别中的聚类方法进行地震属性分析了。曾经运用过的聚类方法有K-均值、自组织特征映射神经网络等方法。这些方法都已经取得了不错的效果。但是K-均值算法在进行海量地震数据聚类时的速度很慢。SOM(Self OrganizingFeature Maps，自组织特征映射神经网络)属性聚类无需事先指定聚类个数，聚类个数和聚类结果一起作为结果输出，但其需要较高的存储空间代价且给出结果速度也比较慢，通常比原始K－均值方法更慢，尤其对几百万或千万个以上的地震属性聚类是一般硬件配置的计算机无法完成的。

发明内容

本发明实施例提供一种地震属性聚类方法，用以在短时间内完成大批量地震属性数据聚类分析，为勘探家提供进一步详细地质分析的基础，该方法包括：

输入预先聚类别个数k及待聚类的地震属性数据{x_i,i＝1,...,N}共N个；

在N个地震属性数据中随机选取其中k个，作为k个类别中心的初始点{m_p,p＝1,...,k}；

为每一个地震属性数据到k个类别中心的距离d(x_i,m_p)估计下界l(x_i,m_p)，初始化均为0；对每一个地震属性数据计算d(x_i,m_p)，选取距离最近的类别中心代表的类别作为该地震属性数据的所属类别，并设置地震属性数据所属类别中心对每一个地震属性数据设置到所属类别中心的距离估计值变量u(x_i)＝min_p＝1,...,kd(x_i,m_p)；

循环执行如下步骤直至类别中心收敛或循环迭代至设定次数，输出每个地震属性数据所属类别标号：

对所有类别中心计算d(m_q,m_p)，q＝1,...,k，并设

确认所有地震属性数据u(x_i)≥s(c(x_i))；

对于所有同时满足条件m_p≠class(x_i)、u(x_i)＞l(x_i,m_p)及的地震属性数据x_i和类别中心m_p：如果标志r(x_i)为真，计算d(x_i,c(x_i))并置r(x_i)为假，否则赋值d(x_i,c(x_i))＝u(x_i)；如果d(x_i,c(x_i))＞l(x_i,m_r)，m_r≠c(x_i)或者计算d(x_i,m_p)；如果d(x_i,m_p)<d(x_i,c(x_i))，赋值c(x_i)＝m_p；

求取每个类别所属地震属性数据的平均值mean(m_p)；

更新每个地震属性数据到所属类别中心的距离下界：l(x_i,c(x_i))＝max{0,l(x_i,c(x_i)-d(c(x_i),mean(c(x_i)))}；

更新每个地震属性数据相关的上下界估计值：u(x_i)＝u(x_i)+d(c(x_i),mean(c(x_i)))，置r(x_i)为真；

对k个类别中心进行更新m_p＝mean(m_p)，使用每个类别所占有的地震属性数据的均值作为新的类别中心；

输入待聚类的地震属性数据之前，还包括对每一个地震属性数据进行高斯归一化处理；

所述高斯归一化处理按如下公式进行：

$y_i=\frac{x_i-\mathrm{\&mu;}}{\mathrm{\&sigma;}};$

其中，{x_i,i＝1,2,...,N}为地震属性数据，{y_i,i＝1,2,...,N}为高斯归一化后的地震属性数据，

一个实施例中，进行所述高斯归一化处理之前，还包括：对地震属性数据中记录错误的原始时间序列信号异常数值进行剔除。

一个实施例中，所述待聚类的地震属性数据包括：瞬时振幅、瞬时频率、瞬时相位、叠后瞬时振幅、反演纵波速度、反演纵横波速度比、振幅随偏移距的变化AVO截距属性、AVO梯度属性其中之一或任意组合。

本发明实施例还提供一种地震属性聚类装置，用以在短时间内完成大批量地震属性数据聚类分析，为勘探家提供进一步详细地质分析的基础，该装置包括：

输入模块，用于输入预先聚类别个数k及待聚类的地震属性数据{x_i,i＝1,...,N}共N个；

类别中心初始点选取模块，用于在N个地震属性数据中随机选取其中k个，作为k个类别中心的初始点{m_p,p＝1,...,k}；

地震属性数据至类别中心距离处理模块，用于为每一个地震属性数据到k个类别中心的距离d(x_i,m_p)估计下界l(x_i,m_p)，初始化均为0；对每一个地震属性数据计算d(x_i,m_p)，选取距离最近的类别中心代表的类别作为该地震属性数据的所属类别，并设置地震属性数据所属类别中心对每一个地震属性数据设置到所属类别中心的距离估计值变量u(x_i)＝min_p＝1,...,kd(x_i,m_p)；

循环迭代处理及输出模块，用于循环执行如下步骤直至类别中心收敛或循环迭代至设定次数，输出每个地震属性数据所属类别标号：

对所有类别中心计算d(m_q,m_p)，q＝1,...,k，并设

确认所有地震属性数据u(x_i)≥s(c(x_i))；

对于所有同时满足条件m_p≠class(x_i)、u(x_i)＞l(x_i,m_p)及的地震属性数据x_i和类别中心m_p：如果标志r(x_i)为真，计算d(x_i,c(x_i))并置r(x_i)为假，否则赋值d(x_i,c(x_i))＝u(x_i)；如果d(x_i,c(x_i))＞l(x_i,m_r)，m_r≠c(x_i)或者计算d(x_i,m_p)；如果d(x_i,m_p)<d(x_i,c(x_i))，赋值c(x_i)＝m_p；

求取每个类别所属地震属性数据的平均值mean(m_p)；

更新每个地震属性数据到所属类别中心的距离下界：l(x_i,c(x_i))＝max{0,l(x_i,c(x_i)-d(c(x_i),mean(c(x_i)))}；

更新每个地震属性数据相关的上下界估计值：u(x_i)＝u(x_i)+d(c(x_i),mean(c(x_i)))，置r(x_i)为真；

对k个类别中心进行更新m_p＝mean(m_p)，使用每个类别所占有的地震属性数据的均值作为新的类别中心；

高斯归一化处理模块，用于在输入待聚类的地震属性数据之前，对每一个地震属性数据进行高斯归一化处理；

所述高斯归一化处理模块具体用于按如下公式进行所述高斯归一化处理：

$y_i=\frac{x_i-\mathrm{\&mu;}}{\mathrm{\&sigma;}};$

其中，{x_i,i＝1,2,...,N}为地震属性数据，{y_i,i＝1,2,...,N}为高斯归一化后的地震属性数据，

一个实施例中，上述地震属性聚类装置还包括：

异常数值剔除模块，用于在进行所述高斯归一化处理之前，对地震属性数据中记录错误的原始时间序列信号异常数值进行剔除。

一个实施例中，所述输入模块具体用于输入包括瞬时振幅、瞬时频率、瞬时相位、叠后瞬时振幅、反演纵波速度、反演纵横波速度比、AVO截距属性、AVO梯度属性其中之一或任意组合的待聚类的地震属性数据。

本发明实施例可以对海量的地震属性数据聚类，只需要能够放入地震属性数据的内存即可，相对于除了地震数据本身还需存储大量网络模型参数的自组织映射神经网络，节省了大量的存储空间；并且在继承了原始K-均值方法线性计算复杂度的优点同时，利用了三角不等式节省了很多不必要的二点间的距离计算，大大提高了方法的运算速度；相比于当前流行自组织映射神经网络聚类方法可以在占用更少量硬件资源的情况下更加快速给出聚类结果。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。在附图中：

图1为本发明实施例中地震属性聚类方法中一个具体实例的实施过程示意图；

图2为本发明实施例中地震属性聚类方法聚类过程中输入输出数据的一个示例的示意图；

图3为本发明实施例中三维数据体中一条线上的聚类结果示意图；

图4为本发明实施例中地震属性聚类装置的结构示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚明白，下面结合附图对本发明实施例做进一步详细说明。在此，本发明的示意性实施例及其说明用于解释本发明，但并不作为对本发明的限定。

为了在短时间内完成大批量地震属性数据聚类分析，为勘探家提供进一步详细地质分析的基础，本发明实施例提例一种地震属性聚类方法，运用快速K-均值方法对大规模地震属性完成快速聚类，可称为基于快速K-均值的地震属性聚类方法。

本发明实施例的基于快速K-均值的地震属性聚类方法是对原始K-均值方法的改进，其主要贡献在于减少了原始K-均值方法中很多不必要的运算，使得运算速度大大提高。

原始K-均值方法如下：

输入：预先聚类别个数k；待聚类的地震属性数据{x_i,i＝1,...,N}共N个；

输出：每个地震属性数据所属的类别标号；

1、在N个地震属性数据中随机选取其中k个，作为k个类别中心的初始点{m_p,p＝1,...,k}；

2、对每一个地震属性数据，计算它到k个类别中心的距离d(x_i,m_p)，i＝1,...,N，p＝1,...,k，选取距离它最近的那个类别中心代表的类别作为这个地震属性数据的所属类别，class(x_i)＝argmin_p＝1,...,kd(x_i,m_p)；

3、对k个类别中心进行更新使用每个类别所占有的地震属性数据的均值来作为新的类别中心；

4、回到步骤2，循环直到类别中心的变化小于某个阈值或循环到一定次数为止。

原始K-均值方法每一次循环迭代都要计算每个地震属性数据到类别中心的距离d(x_i,m_r)用来更新每一个地震属性数据所属类别标号。如果地震属性数据个数很多，那么每次迭代过程中都会进行很多次这样的距离计算，如此频繁计算非常耗费时间。本发明实施例的基于快速K-均值的地震属性聚类方法正是针对这一问题作出改进，根据三角不等式的原理省去了一部分不必要的计算，从而大大提高了运行速度，因此十分适合海量地震属性数据聚类。本发明实施例的基于快速K-均值的地震属性聚类方法如下：

输入：预先聚类别个数k；待聚类的地震属性数据{x_i,i＝1,...,N}共N个；

输出：每个地震属性数据所属的类别标号；

1、在N个地震属性数据中随机选取其中k个，作为k个类别中心的初始点{m_p,p＝1,...,k}；

2、为每个d(x_i,m_p)估计下界l(x_i,m_p)，初始化均为0；对每一个地震属性数据x_i计算它到k个类别中心的距离d(x_i,m_p)，i＝1,...,N，p＝1,...,k，选取距离它最近的那个类别中心代表的类别作为这个地震属性数据的所属类别，并设置地震属性数据所属类别中心对每一个地震属性数据设置变量u(x_i)＝min_p＝1,...,kd(x_i,m_p)，即地震属性数据到它所属类别中心的距离估计值；

3、对所有类别中心计算d(m_q,m_p)，q,p＝1,...,k，并设

4、确认所有地震属性数据u(x_i)≥s(c(x_i))；

5、对于所有同时满足以下条件的地震属性数据x_i和类别中心m_p：

(1)m_p≠class(x_i)；

(2)u(x_i)＞l(x_i,m_p)；

(3)

进行如下操作：如果标志为r(x_i)为真，计算d(x_i,c(x_i))并置r(x_i)为假，否则赋值d(x_i,c(x_i))＝u(x_i)；如果d(x_i,c(x_i))＞l(x_i,m_r)，m_r≠c(x_i)或者计算d(x_i,m_p)；如果d(x_i,m_p)<d(x_i,c(x_i))，那么赋值c(x_i)＝m_p；

6、求每个类别所属地震属性数据的平均值mean(m_p)，p＝1,...,k；

7、更新每个地震属性数据到所属类别的中心距离下界：l(x_i,c(x_i))＝max{0,l(x_i,c(x_i)-d(c(x_i),mean(c(x_i)))}；

8、更新每个地震属性数据相关的上下界估计值：u(x_i)＝u(x_i)+d(c(x_i),mean(c(x_i)))，置r(x_i)为真；

9、对k个类别中心进行更新m_p＝mean(m_p)，使用每个类别所占有的地震属性数据的均值来作为新的类别中心；

10、回到步骤3，直到类别中心收敛或者循环迭代到一定次数。

本发明实施例的基于快速K-均值的地震属性聚类方法中，需要输入的是地震属性的原始数据，即从接收地震反射波记录仪器上记录下来的原始时间序列经过变换抽取得到的一系列属性值。这里可以用到八个属性值，都是瞬时属性：瞬时振幅、瞬时频率、瞬时相位、叠后瞬时振幅、反演纵波速度、反演纵横波速度比、AVO(Amplitude Versus Offset，振幅随偏移距的变化)截距属性(P属性)、AVO梯度属性(G属性)。这几个属性可以部分或者全部用于聚类分析，其中前四个属性主要能够较好地反映地质情况差异，而后四个属性对油气与非油气能够做到明显的区分功能。当然，也可以使用其它的已有或者自定义的属性。如果属性本身不是数值的可以通过某些变换将其转换为数值类型。当然属性的选择需要考虑属性是否能够反映一些地质或者油气本身的一些特有性质，这取决于油气勘探的具体目标。

输出聚类结果是每个属性值矢量所属于的类别标号。如果事先定义了待勘探区域有k个类别(k是正整数)，那么每一个属性矢量所对应的类别标号(正整数)应该是正整数集合{1,2,...,k}中的一个。每一个正整数标号对应着一种地质类别。例如标号1代表着泥岩，标号2代表着碳酸岩，标号3代表着油等等。具体每一个标号代表着什么需要聚类结果给出后由专家或其它方式给出，本发明实施例旨在通过聚类分析将地下各个地点地质属性之间的差异性和相似性刻画出来。

具体实施时，本发明实施例的基于快速K-均值的地震属性聚类方法中，还可以包括为快速聚类而设计的地震属性数据预处理。

快速K-均值算法通常要求输入数据能有相似的范围，即每个属性都有一致或相似的数值分布范围，而原始地震属性数据可能不具备这一性质。输入的原始数据即为采集的地震时间序列信号经过变换得到的若干特征属性值。因为这些特征属性值可能会分布在不同的数值范围，它们之间的差异可能会很大，为了便于统一处理它们，使得所有属性不会因为它们的数值分布范围的差异对聚类结果产生影响，因此可以对原始属性数据进行归一化处理。例如瞬时振幅这个属性变化范围会非常大，最大值和最小值相差可高达10个以上的数量级，在数据中，二个瞬时振幅最大值达到了10⁹而最小在10^-4，叠后振幅情况也类似。相比之下瞬时相位属性值仅在[-π,π]之间变化而且变化相对更为均匀。因此可以采用高斯归一化的方法对每一个属性进行处理。

具体的，高斯归一化处理可以按如下公式进行：

$y_i=\frac{x_i-\mathrm{\&mu;}}{\mathrm{\&sigma;}};$

其中，{x_i,i＝1,2,...,N}为地震属性数据，{y_i,i＝1,2,...,N}为高斯归一化后的地震属性数据，公式中输入{x_i,i＝1,2,...,N}是其中某一个属性原始数值序列，输出{y_i,i＝1,2,...,N}是这个属性经过变换后的数值序列，序列{x_i,i＝1,2,...,N}和{y_i,i＝1,2,...,N}个数是相同的。

具体实施时，在高斯归一化之前，还可以对一些因信号记录设备故障或其他原因而导致的记录错误的一些原始时间序列信号异常数值进行剔除。通常这些异常值是不在属性的合法值域范围内的，例如异常属性值会出现非数值的字符，在使用的地震属性数据中通常会出现“NaN”或“#”字符来表示非法属性值；或者明显超出属性的合法范围，在使用的瞬时相位应该位于[-π,π]，但有的数据已经在100以上这必然是非法的；此类属性还有瞬时频率应该大于0，反演纵波速度、纵横波速比大于0。因此首先可以检查每一个具体属性值是否是数值的，若是非数值将其去掉。然后如果哪个属性有确定的数值范围，对这个属性的每一个具体数值进行越界检查，若确实越界也将其去掉。去掉后使用和它位于同一时间序列上一时刻的属性值来代替。在经过以上的这些异常数值剔除、高斯归一化预处理后，每个属性的数值就会基本分布在一致范围了。

图1为本发明实施例中地震属性聚类方法中一个具体实例的实施过程示意图。参见图1，并结合前述实施例可知，本发明实施例中地震属性聚类方法较佳的实施过程是包括两个主要步骤：地震数据预处理和地震数据聚类。最初的输入是地震属性数据和需要聚类的类别数，输出为每个地震数据所属的类别。

本发明实施例中地震属性聚类方法可以将预处理后地震属性数据运用快速K-均值方法进行聚类，然后将聚类结果所有地震属性矢量的类别标号输出。通常相近位置的类别标号相同的可能性较大，因为地下的地质突变只在很少的地方发生，在大部分区域都能保证地质特性的连续相似性。图2给出了本发明实施例中地震属性聚类方法聚类过程中输入输出数据的一个示例，输入为四个二维的地震属性数据并且设定分为2个类别，显然第一个与第二个比其它两个更相似，而第三个属性矢量与第四个比其它更相似。因此将前两个归为标签为1的一类，而后两个归为标签为2的一类。

为了说明本发明实施例中基于快速K-均值属性聚类的速度优势，这里进行实验对比。将原始K-均值的效果和快速K-均值的效果进行比较。实验平台硬件配置是16核服务器72G内存，64位Windows Server 2003操作系统。使用的数据是某井区的一个正方形切块，水平x坐标(inline方向)和水平y坐标(crossline方向)上均有401个点，时间坐标(向地下的纵深方向，每2毫秒采一点)250个点。总共获取401×401×250＝40200250个点的属性。采用了四组地震属性数据：瞬时振幅、瞬时频率、瞬时相位、瞬时叠后振幅。属性组已经经过预处理。预先设定聚类个数为5类，原始K-均值和快速K-均值聚类的结果是完全一致的，因为是三维立体数据不方便显示，只在某一inline方向坐标值上抽取一个二维切片A(共有100250个采样点)显示其聚类结果，如图3所示。

为了体现快速K-均值的速度优势，给出了原始K-均值和快速K-均值二种聚类方法在取切片A上不同数量采样点的运行时间(缺省单位为秒)。可以从表一的运行时间对比中看到，当属性点个数大于1000时，快速K-均值速度至少比原始K-均值速度快了2个时间单位数量级。

表一原始K-均值聚类和快速K-均值聚类运行时间对比

	原始K-均值	快速K-均值
			A中100个点	0.0285(秒)	0.0044(秒)
A中1000个点	1.264	0.028
			A中10000个点	45.81	0.49
切片A所有100250个点	666.6	4.6
			三维数据体40200250个点	0.5(小时)	大于72(小时)

基于同一发明构思，本发明实施例中还提供了一种地震属性聚类装置，如下面的实施例所述。由于地震属性聚类装置解决问题的原理与地震属性聚类方法相似，因此地震属性聚类装置的实施可以参见地震属性聚类方法的实施，重复之处不再赘述。

图4为本发明实施例中地震属性聚类装置的结构示意图。如图4所示，本发明实施例中地震属性聚类装置可以包括：

输入模块401，用于输入预先聚类别个数k及待聚类的地震属性数据{x_i,i＝1,...,N}共N个；

类别中心初始点选取模块402，用于在N个地震属性数据中随机选取其中k个，作为k个类别中心的初始点{m_p,p＝1,...,k}；

地震属性数据至类别中心距离处理模块403，用于为每一个地震属性数据到k个类别中心的距离d(x_i,m_p)估计下界l(x_i,m_p)，初始化均为0；对每一个地震属性数据计算d(x_i,m_p)，选取距离最近的类别中心代表的类别作为该地震属性数据的所属类别，并设置地震属性数据所属类别中心对每一个地震属性数据设置到所属类别中心的距离估计值变量u(x_i)＝min_p＝1,...,kd(x_i,m_p)；

循环迭代处理及输出模块404，用于循环执行如下步骤直至类别中心收敛或循环迭代至设定次数，输出每个地震属性数据所属类别标号：

对所有类别中心计算d(m_q,m_p)，q＝1,...,k，并设

确认所有地震属性数据u(x_i)≥s(c(x_i))；

对于所有同时满足条件m_p≠class(x_i)、u(x_i)＞l(x_i,m_p)及的地震属性数据x_i和类别中心m_p：如果标志r(x_i)为真，计算d(x_i,c(x_i))并置r(x_i)为假，否则赋值d(x_i,c(x_i))＝u(x_i)；如果d(x_i,c(x_i))＞l(x_i,m_r)，m_r≠c(x_i)或者计算d(x_i,m_p)；如果d(x_i,m_p)<d(x_i,c(x_i))，赋值c(x_i)＝m_p；

求取每个类别所属地震属性数据的平均值mean(m_p)；

更新每个地震属性数据到所属类别中心的距离下界：l(x_i,c(x_i))＝max{0,l(x_i,c(x_i)-d(c(x_i),mean(c(x_i)))}；

更新每个地震属性数据相关的上下界估计值：u(x_i)＝u(x_i)+d(c(x_i),mean(c(x_i)))，置r(x_i)为真；

对k个类别中心进行更新m_p＝mean(m_p)，使用每个类别所占有的地震属性数据的均值作为新的类别中心。

具体实施时，本发明实施例中地震属性聚类装置还可以包括：

高斯归一化处理模块，用于在输入待聚类的地震属性数据之前，对每一个地震属性数据进行高斯归一化处理。

具体实施时，高斯归一化处理模块具体可以用于按如下公式进行所述高斯归一化处理：

$y_i=\frac{x_i-\mathrm{\&mu;}}{\mathrm{\&sigma;}};$

其中，{x_i,i＝1,2,...,N}为地震属性数据，{y_i,i＝1,2,...,N}为高斯归一化后的地震属性数据，

具体实施时，本发明实施例中地震属性聚类装置还可以包括：

异常数值剔除模块，用于在进行所述高斯归一化处理之前，对地震属性数据中记录错误的原始时间序列信号异常数值进行剔除。

具体实施时，输入模块401具体可以用于输入包括瞬时振幅、瞬时频率、瞬时相位、叠后瞬时振幅、反演纵波速度、反演纵横波速度比、AVO截距属性、AVO梯度属性其中之一或任意组合的待聚类的地震属性数据。

综上所述，本发明实施例是采用一种原始K-均值方法的改进版本——快速K-均值方法来对海量的地震属性数据聚类。这种快速K-均值方法只需要能够放入地震属性数据的内存即可，相对于除了地震数据本身还需存储大量网络模型参数的自组织映射神经网络，节省了大量的存储空间。另外快速K-均值在继承了原始K-均值方法线性计算复杂度的优点同时，利用了三角不等式节省了很多不必要的二点间的距离计算，大大提高了方法的运算速度。所以说，本发明实施例相比于当前流行自组织映射神经网络聚类方法可以在占用更少量硬件资源的情况下更加快速给出聚类结果。

本发明实施例拥有着很广的适用范围，可以进行数值化的地震属性都可以使用快速K-均值算法进行聚类分析。本发明实施例可以应用在可以提取出数值化地震属性的油气勘探区域的地质分析，与勘探区域本身的地质条件复杂程度无关，拥有很广阔的应用前景。

本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种地震属性聚类方法，其特征在于，包括：

输入预先聚类别个数k及待聚类的地震属性数据{x_i,i＝1,...,N}共N个；

在N个地震属性数据中随机选取其中k个，作为k个类别中心的初始点{m_p,p＝1,...,k}；

为每一个地震属性数据到k个类别中心的距离d(x_i,m_p)估计下界l(x_i,m_p)，初始化均为0；对每一个地震属性数据计算d(x_i,m_p)，选取距离最近的类别中心代表的类别作为该地震属性数据的所属类别，并设置地震属性数据所属类别中心对每一个地震属性数据设置到所属类别中心的距离估计值变量u(x_i)＝min_p＝1,...,kd(x_i,m_p)；

循环执行如下步骤直至类别中心收敛或循环迭代至设定次数，输出每个地震属性数据所属类别标号：

对所有类别中心计算d(m_q,m_p)，q＝1,...,k，并设

确认所有地震属性数据u(x_i)≥s(c(x_i))；

对于所有同时满足条件m_p≠class(x_i)、u(x_i)>l(x_i,m_p)及的地震属性数据x_i和类别中心m_p：如果标志r(x_i)为真，计算d(x_i,c(x_i))并置r(x_i)为假，否则赋值d(x_i,c(x_i))＝u(x_i)；如果d(x_i,c(x_i))>l(x_i,m_r)，m_r≠c(x_i)或者计算d(x_i,m_p)；如果d(x_i,m_p)<d(x_i,c(x_i))，赋值c(x_i)＝m_p；

求取每个类别所属地震属性数据的平均值mean(m_p)；

更新每个地震属性数据到所属类别中心的距离下界：l(x_i,c(x_i))＝max{0,l(x_i,c(x_i)-d(c(x_i),mean(c(x_i)))}；

更新每个地震属性数据相关的上下界估计值：u(x_i)＝u(x_i)+d(c(x_i),mean(c(x_i)))，置r(x_i)为真；

对k个类别中心进行更新m_p＝mean(m_p)，使用每个类别所占有的地震属性数据的均值作为新的类别中心；

输入待聚类的地震属性数据之前，还包括对每一个地震属性数据进行高斯归一化处理；

所述高斯归一化处理按如下公式进行：

$y_i=\frac{x_i-\mathrm{\&mu;}}{\mathrm{\&sigma;}};$

其中，{x_i,i＝1,2,...,N}为地震属性数据，{y_i,i＝1,2,...,N}为高斯归一化后的地震属性数据，

2.如权利要求1所述的地震属性聚类方法，其特征在于，进行所述高斯归一化处理之前，还包括：对地震属性数据中记录错误的原始时间序列信号异常数值进行剔除。

3.如权利要求1至2任一项所述的地震属性聚类方法，其特征在于，所述待聚类的地震属性数据包括：瞬时振幅、瞬时频率、瞬时相位、叠后瞬时振幅、反演纵波速度、反演纵横波速度比、振幅随偏移距的变化AVO截距属性、AVO梯度属性其中之一或任意组合。

4.一种地震属性聚类装置，其特征在于，包括：

输入模块，用于输入预先聚类别个数k及待聚类的地震属性数据{x_i,i＝1,...,N}共N个；

类别中心初始点选取模块，用于在N个地震属性数据中随机选取其中k个，作为k个类别中心的初始点{m_p,p＝1,...,k}；

地震属性数据至类别中心距离处理模块，用于为每一个地震属性数据到k个类别中心的距离d(x_i,m_p)估计下界l(x_i,m_p)，初始化均为0；对每一个地震属性数据计算d(x_i,m_p)，选取距离最近的类别中心代表的类别作为该地震属性数据的所属类别，并设置地震属性数据所属类别中心对每一个地震属性数据设置到所属类别中心的距离估计值变量u(x_i)＝min_p＝1,...,kd(x_i,m_p)；

循环迭代处理及输出模块，用于循环执行如下步骤直至类别中心收敛或循环迭代至设定次数，输出每个地震属性数据所属类别标号：

对所有类别中心计算d(m_q,m_p)，q＝1,...,k，并设

确认所有地震属性数据u(x_i)≥s(c(x_i))；

对于所有同时满足条件m_p≠class(x_i)、u(x_i)>l(x_i,m_p)及的地震属性数据x_i和类别中心m_p：如果标志r(x_i)为真，计算d(x_i,c(x_i))并置r(x_i)为假，否则赋值d(x_i,c(x_i))＝u(x_i)；如果d(x_i,c(x_i))>l(x_i,m_r)，m_r≠c(x_i)或者计算d(x_i,m_p)；如果d(x_i,m_p)<d(x_i,c(x_i))，赋值c(x_i)＝m_p；

求取每个类别所属地震属性数据的平均值mean(m_p)；

更新每个地震属性数据到所属类别中心的距离下界：l(x_i,c(x_i))＝max{0,l(x_i,c(x_i)-d(c(x_i),mean(c(x_i)))}；

更新每个地震属性数据相关的上下界估计值：u(x_i)＝u(x_i)+d(c(x_i),mean(c(x_i)))，置r(x_i)为真；

对k个类别中心进行更新m_p＝mean(m_p)，使用每个类别所占有的地震属性数据的均值作为新的类别中心；

高斯归一化处理模块，用于在输入待聚类的地震属性数据之前，对每一个地震属性数据进行高斯归一化处理；

所述高斯归一化处理模块具体用于按如下公式进行所述高斯归一化处理：

$y_i=\frac{x_i-\mathrm{\&mu;}}{\mathrm{\&sigma;}};$

其中，{x_i,i＝1,2,...,N}为地震属性数据，{y_i,i＝1,2,...,N}为高斯归一化后的地震属性数据，

5.如权利要求4所述的地震属性聚类装置，其特征在于，还包括：

异常数值剔除模块，用于在进行所述高斯归一化处理之前，对地震属性数据中记录错误的原始时间序列信号异常数值进行剔除。

6.如权利要求4至5任一项所述的地震属性聚类装置，其特征在于，所述输入模块具体用于输入包括瞬时振幅、瞬时频率、瞬时相位、叠后瞬时振幅、反演纵波速度、反演纵横波速度比、AVO截距属性、AVO梯度属性其中之一或任意组合的待聚类的地震属性数据。