本说明书涉及一种频率域多数据体变时窗波形聚类方法及其装置。
背景技术：
地震道波形是地下地质体岩性、物性、流体性质、结构、构造的综合响应，任何与地震波传播有关的地质和地球物理参数都可以通过地震道波形的变化表现出来。因此，根据地震道的波形变化进行波形分类，找出地震道波形变化的总体规律，就有可能认识沉积相和岩相的变化规律。
地震波形分类技术本质上可转化为对地震波形的聚类问题。然而，传统的地震波形分类方法大都是以等时窗内地震道波形相似性来划分地震相，而陆相地层厚度变化大，如果在时间域内将地震信号重采样成等长信号，容易发生波形畸变，无法对目的层段厚度变化大的情况进行精确的地震相分析。
技术实现要素：
本说明书的目的在于提供一种频率域多数据体变时窗波形聚类方法及其装置，其可处理地层厚度变化大、岩性复杂地区的波形分类问题。
为达到上述目的，一方面，本说明书提供了一种频率域多数据体变时窗波形聚类方法，该方法包括：
确定目的层段的分析时窗，并获取多种不同的地震属性数据体；
对上述多种不同的地震属性数据体进行降维，获得融合数据体；
将所述融合数据体从时间域变换为频率域，然后对频率域信号重采样成等长信号，获得频率域数据体；
对所述频率域数据体进行聚类分析。
在上述频率域多数据体变时窗波形聚类方法中，优选地，在所述对上述多种不同的地震属性数据体进行降维，获得融合数据体步骤中，采用pca方法降维。
在上述频率域多数据体变时窗波形聚类方法中，优选地，在所述将所述融合数据体从时间域变换为频率域，然后对频率域信号重采样成等长信号，获得频率域数据体的步骤中，采用频率域信号插值法将所述融合数据体从时间域变换为频率域；所述频率域信号插值法中采用以下频率域sinc插值公式：
上述公式中，a(f)为插值函数值；an为频率域信号；n为离散点号；δf为频率域采样间隔；f为待插值点处频率值。
在上述频率域多数据体变时窗波形聚类方法中，优选地，所述频率域sinc插值公式改进后的离散sinc插值公式为：
上述公式中，a(f)为插值函数值；an为频率域信号；n为离散点号；δf为频率域采样间隔；f为待插值点处频率值。
在上述频率域多数据体变时窗波形聚类方法中，优选地，在所述对所述频率域数据体进行聚类分析步骤中，采用fcm法进行聚类分析。
另一方面，本说明书还提供了一种频率域多数据体变时窗波形聚类装置，该装置包括：
第一模块，所述第一模块用于确定目的层段的分析时窗，并获取多种不同的地震属性数据体；
第二模块，所述第二模块用于对上述多种不同的地震属性数据体进行降维，获得融合数据体；
第三模块，所述第三模块用于将所述融合数据体从时间域变换为频率域，然后对频率域信号重采样成等长信号，获得频率域数据体；
第四模块，所述第四模块用于对所述频率域数据体进行聚类分析。
在上述频率域多数据体变时窗波形聚类装置中，优选地，在所述第二模块中，采用pca方法降维。
在上述频率域多数据体变时窗波形聚类装置中，优选地，在所述第三模块中，采用频率域信号插值法将所述融合数据体从时间域变换为频率域；所述频率域信号插值法中采用以下频率域sinc插值公式：
上述公式中，a(f)为插值函数值；an为频率域信号；n为离散点号；δf为频率域采样间隔；f为待插值点处频率值。
在上述频率域多数据体变时窗波形聚类装置中，优选地，所述频率域sinc插值公式改进后的离散sinc插值公式为：
上述公式中，a(f)为插值函数值；an为频率域信号；n为离散点号；δf为频率域采样间隔；f为待插值点处频率值。
在上述频率域多数据体变时窗波形聚类装置中，优选地，在所述第四模块中，采用fcm法进行聚类分析。
本说明书提供的频率域多数据体变时窗波形聚类方法及其装置，在确定分析时窗后，先通过提取多种不同属性的数据体，以体现沉积体的不同特性；然后通过降维处理，从多个属性数据体中，提取出与研究对象相关性好的信息并将其融合成一个数据体；再将融合数据体从时间域变换到频率域，并对频率域信号重采样成等长信号，可有效克服时间域地震信号重采样导致的波形畸变；最后通过频率域数据体进行聚类分析，形成一套频率域多数据体变时窗波形聚类方案。该方案可在地层厚度变化大的情况下开展波形聚类分析，并有效结合了多种信息，可以更好地分析沉积类型。
附图说明
为了更清楚地说明本说明书实施方式或现有技术中的技术方案，下面将对实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本说明书中记载的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。在附图中：
图1为本说明书一种实施方式中频率域多数据体变时窗波形聚类方法流程示意图；
图2为本说明书一种实施方式中2x地区w1井柱状图；
图3为本说明书一种实施方式中3x地区滨4油组沉积相图；
图4为本说明书一种实施方式中x地区滨4油组时间域波形聚类图；
图5为本说明书一种实施方式中过w1井和w2井的一条地震剖面图；
图6为本说明书一种实施方式中x地区滨4油组多属性波形聚类图。
具体实施方式
为了使本技术领域的人员更好地理解本申请中的技术方案，下面将结合本申请实施方式中的附图，对本申请实施方式中的技术方案进行清楚、完整地描述，显然，所描述的实施方式仅仅是本申请一部分实施方式，而不是全部的实施方式。基于本申请中的实施方式，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施方式，都应当属于本申请保护的范围。
针对陆相地层沉积速率横向变化快、沉积相变大、沉积地层厚度变化快，传统的单属性、等时窗波形聚类方法难以满足沉积分类需求。如何有效结合多种地震沉积信息，开展可变时窗的波形聚类是本领域的技术难点。为此，
参考图1所示，本说明书一实施方式提供了频率域多数据体变时窗波形聚类方法，包括以下步骤：
s1、确定目的层段的分析时窗，并获取多种不同的地震属性数据体；
s2、对上述多种不同的地震属性数据体进行降维，获得融合数据体；
s3、将融合数据体从时间域变换为频率域，然后对频率域信号重采样成等长信号，获得频率域数据体；
s4、对频率域数据体进行聚类分析。
在上述实施方式中，确定分析时窗后，先通过提取多种不同属性的数据体，以体现沉积体的不同特性；然后通过降维处理，从多个属性数据体中，提取出与研究对象相关性好的信息并将其融合成一个数据体；再将融合数据体从时间域变换到频率域，并对频率域信号重采样成等长信号，可有效克服时间域地震信号重采样导致的波形畸变；最后通过频率域数据体进行聚类分析，形成一套频率域多数据体变时窗波形聚类方案。该方案可在地层厚度变化大的情况下开展波形聚类分析，并有效结合了多种信息，可以更好地分析沉积类型。
对于s1步骤，在上述实施方式中，目的层段的分析时窗的确定，可以根据本领域的常规方式和指标确定。
在上述实施方式中，本领域中已经有多种可以从单一的地震数据体中提取多种不同的地质属性数据体的方式，对于本实施方式而言，对此并无具体限定，本领域技术人员可以根据实际情况确定合适的方式。
对于s2步骤，在上述实施方式中，本领域技术人员可以采用常规的降维方式对多种不同的地震属性数据体进行处理，从而获得融合数据体。
在本说明书的一些优选实施方式中，采用pca方法(主成分分析法)降维。主成分分析(principalcomponentanalysis,pca)的核心思想就是对数据进行降维。基本原理是用最少的相互独立的变量取代原有的多维变量来概括信息，而每一个独立的变量都代表某一方面的性质。
设合成信号x(x1，x2，…，xn)的主成分f(x1，x2，…，xn)表示为它的线性组合，其中系数矩阵aij。在数学上可以证明，变量x的协方差矩阵的前m个较大的特征值λi，它所对应的特征向量就是主成分f的系数矩阵aij：
主成分分析的步骤如下：
①计算信号x自身协方差矩阵如下式：
②求出协方差矩阵相应特征值λ(λ1xλ2，…，λp)，并进行最大最小的排序
|r-λe|＝0(3)
③主成分的选择，个数m的确定，通过信息方差累计贡献率gm来确定：
当累积贡献率大于85％时，就认为能足够反映原来变量的信息了，对应的m就是抽取的前m个主成分。
④分别求出对应于特征值λi的特征向量αi，主成分fi＝αi′x。
对于s3步骤，在本说明书的一些实施方式中，采用频率域信号插值法将所述融合数据体从时间域变换为频率域。频率域信号重采样是对复函数进行插值。结合信号处理理论，可以推导出频率域信号插值公式：
式(5)中，a(f)为插值函数值；αn为频率域信号；n为离散点号；δf为频率域采样间隔；f为待插值点处频率值。
在本说明书的一些实施方式中，由于精确的sinc插值在现实中很难实现，因为它需要信号无限长，而实际信号都被截断为有限长离散信号。通过推导，可以得到有限长离散信号频率域的离散sinc插值公式：
式(6)中，a(f)为插值函数值；αn为频率域信号；n为离散点号；δf为频率域采样间隔；f为待插值点处频率值。
从上述公式可以看出，离散sincd插值函数要用到信号频率域中所有点，当n较大时会使计算效率降低。因此，实际应用中一般对sincd插值函数进行截断，取少数几项做插值运算。
对于s4步骤，可以采用本领域的常规方式对频率域数据体进行聚类分析。在本说明书的一些优选实施方式中，采用fcm法进行聚类分析。fcm是模糊聚类方法中的经典算法。其方法是：先将n个点分成c类，定义每一个类有一个聚类中心，然后根据点与聚类中心的距离，形成一些具有相同性质的模糊子集，每一个点与聚类中心有一个隶属度。一个类也就是一系列点组成的模糊子集，每个点对于不同的聚类中心有不同隶属度，由此，可以形成在一定隶属度条件下的分类。
在给定参数下，fcm算法试图分割一个有限元素的集合x＝{x1,x2,…,xn}成为一个c均值模糊聚类。假设有一个有限数据集，fcm会返回一个类矢量，如v＝vi,i＝1,2,…,c和一个分割矩阵u＝uij,i＝1,2,…,c,j＝1,2,…,n。其中uij是在[0,1]之间的一个数值，也就是指明xj属于第i个类的隶属度。fcm的目标函数定义为：
式(7)中，
为达到目标函数最小，聚类中心和隶属度的更新如下：
fcm算法流程如下：
(1)初始化：给定聚类类别数c，2≤c≤n，n是样本个数，设定迭代停止阈值ε，初始化聚类原型模式p(0)，设置迭代计数器b＝0。
(2)更新隶属度矩阵u(b)，
隶属度矩阵就是反映不同波形样本与各聚类中心相似度的矩阵。它的计算原则为：样本与某个聚类中心的距离越近，则对该类的隶属程度越大，而两个与某聚类中心距离相同的样本，如果其中之一与其他聚类中心距离更近，则其隶属度相对另一个样本更小。
(3)更新聚类原型模式矩阵p(b+1)：
如果||p(b)-p(b+1)||＜ε，则算法停止并输出划分矩阵u和聚类原型p，否则令b＝b+1，转向(7)。其中||·||为某种合适的矩阵范数。
由以上步骤可以看出，整个计算过程就是反复修改聚类中心和隶属度矩阵的过程。
fcm算法能从任意给定初始点开始沿一个迭代子序列收敛到其目标函数jm(u,p)的局部极小点。
虽然上文描述的过程流程包括以特定顺序出现的多个操作，但是，应当清楚了解，这些过程可以包括更多或更少的操作，这些操作可以顺序执行或并行执行(例如使用并行处理器或多线程环境)。
下面介绍本说明书一示例性实施方式。图2是x地区w1井柱状图，滨4砂组是主要含油层系之一，故以滨4砂组为例说明。
滨4砂组厚度在100-280m之间，在研究区内发育的时间深度范围在2.6-3.5s。岩心观察和测井分析表明，本区滨4砂组为不同岩性地层的互层沉积，北部发育三角洲前缘沉积，西部发育灰云坪，南部发育滨浅湖滩坝沉积(图3)。这为我们判断波形聚类结果提供了依据。
根据井震标定结果，我们在三维地震工区对滨4砂组的顶、底层位做精细解释，并对滨4砂组做常规单数据体时间域波形聚类，图4是波形聚类图。图4的聚类结果暴露出2个问题，一是聚类结果与已有钻井数据不吻合，如w1井与w2井分别为灰云坪相和三角洲前缘相，但由于地震波形较为相似(图5)，因此分类结果将他们分成了同一类。第二是分类结果杂乱而无规律，与已有钻井揭示的沉积规律不吻合。
为此，我们实施以下方案：
s101、确定目的层段的分析时窗，可以得到多种不同的地震属性数据体。
s201、对上述多种不同的地震属性数据体进行pca降维，获得一个融合数据体；
s301、将融合数据体从时间域变换为频率域，在频率域将滨4砂组的地震资料重采样成等长信号，获得频率域数据体；
s401、对频率域数据体进行fcm聚类分析。
采用fcm算法对频率域重采样后的滨4砂组做波形聚类分析，图6是聚类结果。显然，其聚类结果要明显优于常规的时间域单数据体波形聚类结果(图4)。图6将整个工区分为5类，各个分类区域形态完整，分类结果与井资料完全吻合。值得注意的是，三角洲前缘分成了2类，钻井取芯结果表明，二支三角洲分属不同物源，西部三角洲砂体沉积颗粒粗、沉积厚度大，东部三角洲沉积砂体颗粒较细，沉积厚度小于西部三角洲砂体。分类结果很好地体现了2者的差别。
本说明书实施方式提供的频率域多数据体变时窗波形聚类装置包括：
第一模块，第一模块用于确定目的层段的分析时窗，并获取多种不同的地震属性数据体；
第二模块，第二模块用于对上述多种不同的地震属性数据体进行降维，获得融合数据体；
第三模块，第三模块用于将所述融合数据体从时间域变换为频率域，然后对频率域信号重采样成等长信号，获得频率域数据体；
第四模块，第四模块用于对所述频率域数据体进行聚类分析。
在本说明书的一些实施方式中，在第二模块中，采用pca方法降维。
在本说明书的一些实施方式中，在所述第三模块中，采用频率域信号插值法将所述融合数据体从时间域变换为频率域；所述频率域信号插值法中采用以下频率域sinc插值公式：
上述公式中，a(f)为插值函数值；an为频率域信号；n为离散点号；δf为频率域采样间隔；f为待插值点处频率值。
在本说明书的一些实施方式中，频率域sinc插值公式改进后的离散sinc插值公式为：
上述公式中，a(f)为插值函数值；an为频率域信号；n为离散点号；δf为频率域采样间隔；f为待插值点处频率值。
在本说明书的一些实施方式中，在第四模块中，采用fcm法进行聚类分析。
由于该装置解决问题的原理与频率域多数据体变时窗波形聚类方法相似，因此，频率域多数据体变时窗波形聚类装置的实施可以参见频率域多数据体变时窗波形聚类方法的实施，重复之处不再赘述。
为了描述的方便，描述以上装置时以功能分为各种单元分别描述。当然，在实施本说明书时可以把各单元的功能在同一个或多个软件和/或硬件中实现。
本说明书是参照根据本说明书实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
在一个典型的配置中，计算设备包括一个或多个处理器(cpu)、输入/输出接口、网络接口和内存。
内存可能包括计算机可读介质中的非永久性存储器，随机存取存储器(ram)和/或非易失性内存等形式，如只读存储器(rom)或闪存(flashram)。内存是计算机可读介质的示例。
计算机可读介质包括永久性和非永久性、可移动和非可移动媒体可以由任何方法或技术来实现信息存储。信息可以是计算机可读指令、数据结构、程序的模块或其他数据。计算机的存储介质的例子包括，但不限于相变内存(pram)、静态随机存取存储器(sram)、动态随机存取存储器(dram)、其他类型的随机存取存储器(ram)、只读存储器(rom)、电可擦除可编程只读存储器(eeprom)、快闪记忆体或其他内存技术、只读光盘只读存储器(cd-rom)、数字多功能光盘(dvd)或其他光学存储、磁盒式磁带，磁带磁磁盘存储或其他磁性存储设备或任何其他非传输介质，可用于存储可以被计算设备访问的信息。按照本文中的界定，计算机可读介质不包括暂存电脑可读媒体(transitorymedia)，如调制的数据信号和载波。
还需要说明的是，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、商品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、商品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、商品或者设备中还存在另外的相同要素。
本领域技术人员应明白，本说明书的实施例可提供为方法、系统或计算机程序产品。因此，本说明书可采用完全硬件实施例、完全软件实施例或结合软件和硬件方面的实施例的形式。而且，本说明书可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等)上实施的计算机程序产品的形式。
本说明书可以在由计算机执行的计算机可执行指令的一般上下文中描述，例如程序模块。一般地，程序模块包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件、数据结构等等。也可以在分布式计算环境中实践本说明书，在这些分布式计算环境中，由通过通信网络而被连接的远程处理设备来执行任务。在分布式计算环境中，程序模块可以位于包括存储设备在内的本地和远程计算机存储介质中。
本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于系统实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。
以上所述仅为本说明书的实施例而已，并不用于限制本说明书。对于本领域技术人员来说，本说明书可以有各种更改和变化。凡在本说明书的精神和原理之内所作的任何修改、等同替换、改进等，均应包含在本说明书的权利要求范围之内。