一种基于矩阵降维的特征结构相干体快速实现方法与流程

本发明属于勘探地球物理领域，特别涉及一种基于矩阵降维的特征相干体快速实现方法。

背景技术：

石油天然气是一种重要的战略资源，在国民经济及国家安全中占有举足轻重的地位。地震属性是对地震勘探数据的一种度量，能够在视觉上增强具有解释意义的特征，对地层特征以及储层特性有直接反映。因此，地震属性可在油气藏勘探开发中发挥重要的作用。几何属性是地震资料解释中最为重要的地震属性之一，主要用于增强和显示地震层位的几何形态，主要包括倾角\方位角、连续性等。其中刻画连续性的重要方位为相干体算法，其结果不但可直接用于勘探工区的断层解释，还可为还可用于监控地震资料处理的质量。然而，应用效果非常明显的特征结构相干体运算量包含两个部分：构建协方差矩阵及对协方差矩阵进行特征分解，这极大地限制了它在工业界的广泛应用。因此，提高特征结构相干体算法的计算效率具有极为重要的应用价值。

目前已有多种技术用于提高特征结构相干体的计算效率，主要技术如下所示：

现有技术1：基于空间/时间递推的特征结构相干体快速实现方法

该类方法瞄准于如何减少构建协方差矩阵的运算量。这类方法主要利用相邻分析窗之间的时间及空间相关性，避免在当前分析窗内构建协方差矩阵，采用对前一个分析窗内构建的协方差矩阵进行简单改造以获得当前分析窗内的协方差矩阵的策略减少运算量。

现有技术1的特点：

优点：可大幅减少构建协方差矩阵的运算量。

缺点：1、只能减少特征结构相干体中构建协方差的运算量；2、无法减少特征结构相干体中矩阵特征分解的运算量。

现有技术2：基于乘幂法的特征结构相干体快速实现方法

该类方法瞄准于如何减少协方差矩阵分解的运算量。这类方法避免求取协方差矩阵的全部特征值，转而利用乘幂法来求取协方差矩阵的最大特征值来达到减少计算量的目的。

现有技术2的特点：

优点：实现简单、可以较大幅度地减少矩阵特征分解的计算量。

缺点：随着矩阵分析窗尺寸的增加，矩阵特征分解的计算量大幅上涨。

技术实现要素：

本发明的目的在于提供一种基于矩阵降维的特征结构相干体快速实现方法，通过矩阵降维，能够更为高效地计算三维叠后地震勘探数据的相干体，在有限的时间内为断层解释提供基础数据体。

为了实现上述目的，本发明采用如下技术方案：

一种基于矩阵降维的特征结构相干体快速实现方法，包括以下步骤：

步骤01：采集得到三维地震数据体，并设定二维空间分析窗大小和一维时间分析窗大小；

步骤02：依照预先设定的二维空间分析窗大小从步骤01采集的三维地震数据体中切取得到一个小三维数据体；

步骤03：在小三维数据体中，计算每一时刻x方向相干度量和y方向相干度量；

步骤04：在小三维数据体内将所有时刻的x方向相干度量和y方向相干度量相乘得到最终的相干度量；

步骤05：判断x及y是否遍历完毕，若完毕输出三维相干体，否则改变x及y并返回步骤02。

进一步的，步骤01中，采集得到三维地震数据体，并设定二维空间分析窗大小和一维时间分析窗大小，具体包括：

用f(x,y,t)表示大小为nx×ny×nt的三维地震数据，x和y表示三维地震数据在两个空间方向的指标，t表示三维地震数据在时间方向的指标，其中x和y的变化范围分别为[1,nx]和[1,ny]，t的变化范围为[1,nt]；空间分析窗尺寸为mx×my×mt，其中mx、my及mt为三维分析窗沿x、y及t方向的大小，且其同时均为奇数；mx＝2wx+1，my＝2wy+1，mt＝2wt+1，其中wx、wy及wt为x、y及t方向的半窗大小。

进一步的，步骤02中，依照预先设定的二维空间分析窗大小从步骤01采集的三维地震数据体中切取得到一个小三维数据体，具体包括：

以x及y为中心可从三维地震数据中划取得到一个小三维地震数据体(其大小为mx×my×nt，其在原始三维数据体中的指标范围为[x-wx,x+wx]、[y-wy,y+wy]及[1,nt])。同时假设小三维数据体内沿x方向、y方向的局部指标为mx及my，两者的变化范围为[1,mx]及[1,my]。

进一步的，步骤03中具体包括：

步骤031：依据其前一时刻的x方向协方差矩阵对当前时刻的x方向协方差矩阵进行更新，以获得每一时刻的x方向协方差矩阵；

步骤032：在小三维数据体中对每一时刻，对其对应的x方向协方差矩阵进行特征分解以获得最大特征值，进而获得每一时刻x方向相干度量；

步骤033：在小三维数据体中对每一时刻，依据其前一时刻的y方向协方差矩阵对当前时刻的y方向协方差矩阵进行更新，以获得每一时刻的y方向协方差矩阵；

步骤034：在小三维数据体中在每一时刻，对其对应的y方向协方差矩阵进行特征分解以获得最大特征值，进而获得每一时刻y方向相干度量。

进一步的，步骤031中在小三维数据体中对每一时刻,依据其前一时刻的x方向协方差矩阵对当前时刻的x方向协方差矩阵进行更新,以获得每一时刻的x方向协方差矩阵，具体包括：

x方向协方差矩阵的大小为mx×mx，通过如下方法初始化x方向协方差矩阵第i行第j个元素值cx(x,y,0)ij：

若在t≥1且t≤wt+1时，用如下方法更新x方向协方差矩阵第i行第j个元素值cx(x,y,t)ij：

若在t＞wt+1且t≤nt-wt时，用如下方法更新x方向协方差矩阵第i行第j个元素值cx(x,y,t)ij：

若在t＞nt-wt时，用如下方法更新x方向协方差矩阵第i行第j个元素值cx(x,y,t)ij：

步骤032中，在小三维数据体中对每一时刻，对其对应的x方向协方差矩阵进行特征分解以获得最大特征值，进而获得每一时刻x方向相干度量，具体包括：

对t时刻的x方向协方差矩阵cx(x,y,t)进行特征分解：

其中λx1、λx2…λxmx分别为cx(x,y,t)的特征值，且满足λx1≥λx2≥…≥λxmx,u1、u2…umx为对应的特征向量，则t时刻x方向相干度量ax(x,y,t)由最大特征值λx1与所有mx个特征值之和的比构建：

步骤033中，在小三维数据体中对每一时刻,依据其前一时刻的y方向协方差矩阵对当前时刻的y方向协方差矩阵进行更新,以获得每一时刻的y方向协方差矩阵，具体包括：

y方向协方差矩阵的大小为my×my，通过如下方法初始化y方向协方差矩阵第i行第j个元素值cy(x,y,0)ij：

若在t≥1且t≤wt+1时，用如下方法更新y方向协方差矩阵第i行第j个元素值cy(x,y,t)ij：

若在t＞wt+1且t≤nt-wt时，用如下方法更新y方向协方差矩阵第i行第j个元素值cy(x,y,t)ij：

若在t＞nt-wt时，用如下方法更新y方向协方差矩阵第i行第j个元素值cy(x,y,t)ij：

步骤034，在小三维数据体中每一时刻，对其对应的y方向协方差矩阵进行特征分解以获得最大特征值，进而获得每一时刻y方向相干度量，具体包括：

对cy(x,y,t)进行特征分解：

其中λy1、λy2…λymy分别为cy(x,y,t)的特征值，且满足λy1≥λy2≥…≥λymy,u1、u2…umy为对应的特征向量，则t时刻y方向相干度量ay(x,y,t)由最大特征值λy1与所有my个特征值之和的比来构建：

进一步的，步骤04中，在小三维数据体内将所有时刻的x方向相干度量和y方向相干度量相乘得到最终的相干度量，具体包括：

在小三维数据体内的将任一时刻的x方向相干度量ax(x,y,t)和对应时刻的y方向相干度量ay(x,y,t)相乘得到最终的相干度量a(x,y,t)：

a(x,y,t)＝ax(x,y,t)ay(x,y,t)。

进一步的，步骤05中，判断x及y是否遍历完毕，若完毕输出三维相干体，否则改变x及y并返回步骤02，具体包括：

判断x的取值是否在[1,nx]范围内遍历完毕，判断y的取值是否在[1,ny]范围内遍历完毕；若两者都遍历完毕说明三维数据体已经处理完毕，则输出三维相干度量数据体，否则改变x及y的取值并返回步骤02。

相对于现有技术，本发明具有以下有益效果：本发明为针对三维叠后地震数据计算特征结构相干体以解释断层的快速实现方法，为一种快速数据处理方法；本发明首先在三维叠后数据体上依据预先设定的分析窗尺寸针对每一点构建沿x方向协方差矩阵及沿y方向协方差矩阵，然后分别对x方向的协方差矩阵及y方向协方差矩阵进行特征分解得到x方向相关度量及y方向的相干度量，最终将x方向及y方向的相干度量进行综合以得到最终的特征结构相干体，以更为快捷的速度为后续地震资料解释人员提供可靠的相干体属性。本发明可通过矩阵降维更加快速地得到特征结构相干体属性，现对于现有技术，能够节约处理时间50％以上，分析窗尺寸越大的时候，节约的时间更多，最多可以节约处理时间90％以上。

附图表说明

图1为本发明流程图；

图2为常规方法计算特征结构相干体结果的剖面展示；

图3为本发明方法计算特征结构相干体结果的对应剖面展示。

具体实施方式

下面结合附图及表格对本发明做进一步详细的说明。

本发明为针对三维叠后地震数据计算特征结构相干体以解释断层的快速实现方法，为一种快速数据处理方法。本发明首先在三维叠后数据体上依据预先设定的分析窗尺寸针对每一点构建沿x方向协方差矩阵及沿y方向协方差矩阵，然后分别对x方向的协方差矩阵及y方向协方差矩阵进行特征分解得到x方向相关度量及y方向的相干度量，最终将x方向及y方向的相干度量进行综合以得到最终的特征结构相干体，以更为快捷的速度为后续地震资料解释人员提供可靠的相干体属性。

请参阅图1所示，本发明提供一种基于矩阵降维的特征结构相干体快速实现方法，主要包括：

步骤01：采集得到三维地震数据体，并设定二维空间分析窗大小和一维时间分析窗大小；

步骤02：依照预先设定的二维空间分析窗大小从步骤01采集的三维地震数据体中切取得到一个小三维数据体；

步骤03：在小三维数据体中，计算每一时刻x方向相干度量和y方向相干度量；

步骤04：在小三维数据体内将所有时刻的x方向相干度量和y方向相干度量相乘得到最终的相干度量；

步骤05：判断x及y是否遍历完毕，若完毕输出三维相干体，否则改变x及y并返回步骤02。

步骤03具体包括：

步骤031：依据其前一时刻的x方向协方差矩阵对当前时刻的x方向协方差矩阵进行更新，以获得每一时刻的x方向协方差矩阵；

步骤032：在小三维数据体中对每一时刻，对其对应的x方向协方差矩阵进行特征分解以获得最大特征值，进而获得每一时刻x方向相干度量；

步骤033：在小三维数据体中对每一时刻，依据其前一时刻的y方向协方差矩阵对当前时刻的y方向协方差矩阵进行更新，以获得每一时刻的y方向协方差矩阵；

步骤034：在小三维数据体中在每一时刻，对其对应的y方向协方差矩阵进行特征分解以获得最大特征值，进而获得每一时刻y方向相干度量。

步骤01中采集得到三维地震数据体，并设定二维空间分析窗大小和一维时间分析窗大小，具体包括：

用f(x,y,t)表示大小为nx×ny×nt的三维地震数据，x和y表示三维地震数据在两个空间方向的指标，t表示三维地震数据在时间方向的指标，其中x和y的变化范围分别为[1,nx]和[1,ny]，t的变化范围为[1,nt]。设定的分析窗大小为mx×my×mt(其中mx与my为二维空间分析窗的大小，mt为时间分析窗的大小，且三者均为奇数；由于mx、my及mt均为奇数，可将其表示为mx＝2wx+1，my＝2wy+1，mt＝2wt+1，其中wx、wy及wt为x、y及t方向的半窗大小)。

步骤02中依照预先设定的二维空间分析窗尺寸mx×my并以x及y为中心从三维地震数据中切取得到一个小三维数据体，具体包括：

假设目前关注三维数据体的x指标及y指标，以x及y为中心可从三维地震数据中划取得到一个小三维地震数据体(其大小为mx×my×nt，其在原始三维数据体中的指标范围为[x-wx,x+wx]、[y-wy,y+wy]及[1,nt])。同时假设小三维数据体内沿x方向、y方向的局部指标为mx及my，两者的变化范围为[1,mx]及[1,my]。

步骤031中在小三维数据体中对每一时刻,依据其前一时刻的x方向协方差矩阵对当前时刻的x方向协方差矩阵进行更新,以获得每一时刻的x方向协方差矩阵，具体包括：

x方向协方差矩阵的大小为mx×mx，可通过如下方法初始化x方向协方差矩阵第i行第j个元素值cx(x,y,0)ij：

若在t≥1且t≤wt+1时，用如下方法更新x方向协方差矩阵第i行第j个元素值cx(x,y,t)ij：

若在t＞wt+1且t≤nt-wt时，用如下方法更新x方向协方差矩阵第i行第j个元素值cx(x,y,t)ij：

若在t＞nt-wt时，用如下方法更新x方向协方差矩阵第i行第j个元素值cx(x,y,t)ij：

步骤032中，在小三维数据体中对每一时刻，对其对应的x方向协方差矩阵进行特征分解以获得最大特征值，进而获得每一时刻x方向相干度量，具体包括：

对t时刻的x方向协方差矩阵cx(x,y,t)进行特征分解：

其中λx1、λx2…λxmx分别为cx(x,y,t)的特征值，且满足λx1≥λx2≥…≥λxmx,u1、u2…umx为对应的特征向量，则t时刻x方向相干度量ax(x,y,t)可由最大特征值λx1与所有mx个特征值之和的比来构建：

步骤033中，在小三维数据体中对每一时刻,依据其前一时刻的y方向协方差矩阵对当前时刻的y方向协方差矩阵进行更新,以获得每一时刻的y方向协方差矩阵，具体包括：

y方向协方差矩阵的大小为my×my，可通过如下方法初始化y方向协方差矩阵第i行第j个元素值cy(x,y,0)ij：

若在t≥1且t≤wt+1时，用如下方法更新y方向协方差矩阵第i行第j个元素值cy(x,y,t)ij：

若在t＞wt+1且t≤nt-wt时，用如下方法更新y方向协方差矩阵第i行第j个元素值cy(x,y,t)ij：

若在t＞nt-wt时，用如下方法更新y方向协方差矩阵第i行第j个元素值cy(x,y,t)ij：

步骤034，在小三维数据体中每一时刻，对其对应的y方向协方差矩阵进行特征分解以获得最大特征值，进而获得每一时刻y方向相干度量，具体包括：

对cy(x,y,t)进行特征分解：

其中λy1、λy2…λymy分别为cy(x,y,t)的特征值，且满足λy1≥λy2≥…≥λymy,u1、u2…umy为对应的特征向量，则t时刻y方向相干度量ay(x,y,t)可由最大特征值λy1与所有my个特征值之和的比来构建：

步骤04中，在小三维数据体内将所有时刻的x方向相干度量和y方向相干度量相乘得到最终的相干度量，具体包括：

在小三维数据体内的将任一时刻的x方向相干度量ax(x,y,t)和对应时刻的y方向相干度量ay(x,y,t)相乘得到最终的相干度量a(x,y,t)：

a(x,y,t)＝ax(x,y,t)ay(x,y,t)。

步骤05中，判断x及y是否遍历完毕，若完毕输出三维相干体，否则改变x及y并返回步骤02，具体包括：

判断x的取值是否在[1,nx]范围内遍历完毕，判断y的取值是否在[1,ny]范围内遍历完毕。若两者都遍历完毕说明三维数据体已经处理完毕，则输出三维相干度量数据体，否则改变x及y的取值并返回步骤02。

以某油田的三维叠后地震数据(数据大小为600*200*500)为例，分别采用常规实现方法及本发明方法计算特征结构相干体。常规方法计算得到的特征结构相干体剖面如图2所示，利用本发明方法计算得到的特征结构相干体剖面如图3所示，本发明的结果与常规实现方法基本一致。表1为不同分析窗大小时采用常规实现方法和本专利方法运算时间对比，由图2及图3可以看出，本发明方法可获得与常规实现方法基本相同的效果，但运算时间大幅减少，说明采用本方法可以更加快速的计算三维地震数据体的特征结构相干体。

表1为常规方法和本发明方法计算相同大小数据体的特征结构相干体所用时间对比；

最后需要说明的是，以上算例对本发明的目的，技术方案以及有益效果提供了进一步的验证，这仅属于本发明的具体实施算例，并不用于限定本发明的保护范围，在本发明的精神和原则之内，所做的任何修改，改进或等同替换等，均应在本发明的保护范围内。