0 引 言

断层解释是油气勘探构造解释的关键，解释的精确性和合理性直接影响构造解释的精度.长期以来，研究人员围绕断层系统的精确描述做出了大量的努力，提出并应用了很多描述方法，如常规断层解释方法，断层切片解释方法以及断层的三维数学描述与模拟方法(朱良峰等，2008)等，成功地从三维地震数据中提取了许多突出断层信息的不连续属性体(赵亮，2010).

相干技术是20世纪90年代发展起来的一项十分有效的地震解释技术.基本原理是在偏移后的三维数据体中，对每一地震道分别计算与周围地震数据的相关度(段春节等，2013).当存在断层和地层不连续性变化时，由于地震道波形发生的变化导致了道与道之间的相干性发生变化，引起地震道局部的不连续性(严哲，2010).这种相干技术可以压制连续性，突出不连续性.

第一代相干(Generation-I correlation，C1)算法是由Bahorich和Farmer提出的基于互相关的方法(Bahorichand Farmer，1995)；第二代相干(Generation-Ⅱ correlation，C2)算法是由Marfurt和Kirlin提出的基于多道相似的方法(Marfurt et al.，1998)；第三代相干(Generation-Ⅲ correlation，C3)算法是由Gersztenkorn和Marfurt提出的基于特征结构分析的方法(Gersztenkornand Marfurt，1999).其中，C1算法最简单，计算速度最快，但只用到了相邻的三道数据，受噪声的影响较大；C2和C3算法均为基于协方差矩阵的方法，采用多道数据参与计算，提高了算法抗噪声的能力，但是计算量偏大；C3算法则是类似于基于特征值的主成分分析算法，相比C2算法，进一步提高了相干算法的分辨率和抗噪能力.

高斯过程的高阶累积量(阶次大于2时)等于零，对于非高斯过程，则至少存在大于2的某个阶次的累积量不为零(王书明等，2006).朱培民提出并证明，地震信号的具有对称性和高斯性(朱培民等，2001).因此，可以利用高阶统计量的相干方法可以抑制地震信号以及高斯噪声的影响，突出断层的响应.随着高阶统计量的研究深入，高阶统计量方法的优点逐渐受到人们的关注.高阶统计量分析方法在处理非线性、非高斯性地震信号方面具有独特的优势，已经在地震数据处理和解释中已得到广泛应用.熊晓军等利用高阶统计量比功率谱函数具有对高斯噪声不敏感的优点及地震信号的高阶谱特征，提出了基于高阶统计量的第一代相干(Higher Order Statistics with Generation-I Correlation，HOS-C1)算法，该算法利用地震属性高阶谱特征参数进行油气检测和识别的方法(熊晓军等，2004).王书明等根据高阶统计量有关理论和方法提出了双谱和互双谱的高阶时频分析方法，用于检测碳酸盐岩中由孔、洞、缝等引起的非均匀储集(王书明等，2006).张绍聪利用高阶相干算法提高了相干算法的抗噪能力，用于煤田断层的解释(张绍聪等，2009).杨振武提出了用高阶统计量识别和预测生物礁储集层的新方法，其目的是通过高阶统计量检测出生物礁的非高斯性信息，预测生物礁储层的存在(杨振武等，2010). 综合考虑高阶统计量的优点和已有三代相干技术方法的优缺点，本文从提高相干体算法抑制噪声的能力的角度出发，提出了一种新的基于高阶统计量的相干体方法.该算法以多道相似相干算法为基础，充分考虑了窗口内地震道的能量和相关信息，推导出了类似于第三代相干算法形式的基于高阶统计量特征值相干算法.最终以实际三维地震数据体为例，对提出的相干计算方法进行了验证，实验表明提出的算法有很好的抑制噪声的能力和断层识别的精确度.

1 基于高阶统计量的互相关相干算法 1.1 C1算法

C1算法是以归一化互相关函数为基础，计算相邻测线上地震道之间的相似性.定义纵测线上t时刻相邻的两个地震道u(x_i，y_i)和u(x_i+1，y_i)，在时间延迟为m的互相关函数ρ_x为

(1)

其中，2ω+1为计算相关的时窗长度，u(·)表示地震道记录数据.

同理，横测线方向上相邻的两个地震道u(x_i，y_i)和u(x_i，y_i+1)，在时间延迟为n的互相关函数ρ_y为

(2)

综合考虑式(1)和式(2)，相干值ρ_xy定义为

(3)

其中，maxmρ_x(t，m，x_i，y_i)表示式(1)在不同m延时下的最大值，maxnρ_y(t，n，x_i，y_i)表示式(2)在不同n延时下的最大值.

1.2 HOS-C1算法

由高阶统计量的理论可知，高斯噪声的三阶或三阶以上的累积量为零.考虑到高阶统计量的这个性质和互相关相干算法的优点，熊晓军等提出了HOS-C1算法.HOS-C1算法较C1算法提高了抗噪能力，具体计算过程如下：

(4)

其中，τ₁为(x_i+1，y_i)点地震道相对于(x_i，y_i)点地震道的时间延迟，τ₂为(x_i，y_i+1)点地震道相对于(x_i，y_i)点地震道的时间延迟.HOS-C1算法类似于C1算法中的纵向与横向上的相干系数相乘，该算法构造出了四个地震道u(·)的延时乘积，故称为高阶统计量互相关算法.于是定义

(5)

为最终目标地震道的相干值.

通过比较C1算法与HOS-C1算法，可以看出两种算法都需要计算不同延时情况下相邻三个地震道的相似性.不同之处是，C1算法计算的是2个地震道u(·)的乘积，HOS-C1算法计算的是4个地震道u(·)的乘积.设M和N分别表示ρ_x和ρ_y的最大延时，S表示一道数据的相关计算量，那么C1算法的时间复杂度为O((M+N)·S)，HOS-C1算法的计算复杂度为O(M·N·S).

2 基于高阶统计量的特征值相干算法 2.1 C3算法

C1算法是沿着横、纵测线方向的3道数据进行相干处理，显然这种方法没有充分利用中心地震道附近的地震道信息，C2算法则对上面的想法进行了考虑.该分析方法则是以中心地震道为中心建立一个椭圆形或者矩形分析窗，假设窗口内包含了J道地震数据，进行相干计算.以t时刻的(x_i，y_i)点为中心，对分析窗内的数据进行相似系数的判定，相似系数定义为

(6)

其中，u(t-px_i-qy_i，x_i，y_i)为中心地震道，u(t-px_j-qy_j，x_j，y_j)为中心地震道周围的J条地震道，p和q分别为x和y方向的视倾角.上式即为分析窗内平均道的能量与各道平均能量的比值.

为了保证计算的稳定性，考虑2ω+1为计算相关的时窗长度，即

(7)

令t=nΔt，t为时间，n为正整数，Δt为采样间隔，那么有

(8)

那么，C2算法的相干值为

(9)

其中，u_jm=u(mΔt-px_i-qy_i，x_i，y_i).将相似系数写成协方差矩阵形式

(10)

假设λ_j(j=1，2，…，J)是协方差矩阵Cov_C3(t，p，q，x_i，y_i)的J个特征值.参考主成分分析方法，认为最大特征值对应信号的信息，其余特征值对应的是噪声，则C3算法的相干值定义为

(11)

2.2 基于高阶统计量的特征值相干算法

通过以上高阶统计量在C1算法中的应用，本文提出了基于高阶统计量的特征值相干(Higher Order Statistics with Generation-Ⅲ Correlation，HOS-C3)算法.C2算法是基于目标地震道相邻地震道的幅度展开的，只考虑了周围地震道的幅度大小而忽略了它们与中心目标地震道之间的关系.考虑到以上关系，本文提出同时考虑目标地震道和周围地震道的相关性和能量关系的HOS-C3算法.

定义HOS-C3算法下的相关函数为

(12)

其中，u(t-px_i-qy_i，x_i，y_i)为中心地震道，u(t-px_j-qy_j，x_j，y_j)为中心地震道周围的地震道，p和q分别为x和y方向的视倾角.式(12)同时考虑了分析窗内的能量与各道的相关性.为了保证计算的稳定性，考虑2ω+1为计算相关的时窗长度，即

(13)

对于时间中心为t时，设t=nΔt，视倾角为(p，q)时(13)式对应的协方差矩阵为

(14)

其中，u′_jm=u(mΔt-px_j-qy_j，x_i，y_i)u(mΔt-px_j-qy_j，x_j，y_j)为在时间t=mΔt-px_i-qy_i时中心节点和相邻节点插值的乘积.假设λ_j(j=1，2，…，J)是协方差矩阵Cov_HOS(t，p，q，x_i，y_i)的J个特征值.于是，参考C3算法的相干体算法定义方式，定义

(15)

为HOS-C3算法的相干值.

通过比较C3算法与HOS-C3算法，可以看出两种算法在计算协方差矩阵过程中都需要计算延时情况下相邻地震道的插值的乘积，C3算法计算的是2个地震道u(·)插值的乘积，HOS-C3算法计算的是4个地震道u(·)插值的乘积. 设J表示窗口内地震道数，R表示协方差矩阵计算特征值的计算量，W表示相关计算的延时数，那么C3算法的计算复杂度为O(4J·W·R)，HOS-C3算法的计算复杂度为O(4J·W·R).

3 实际应用

为了验证本文提出相干算法的效果，用实际三维地震数据体进行了断层识别，并与其他相干算法的结果进行了比较.实验采用的实际三维地震数据体是某工区的三维地震数据，选取的数据为横测线方向上120~400共281道，纵测线方向上350~600共251道，数据的采样时间是1400~1700 ms，采样间隔是4 ms，共76个采样点.实验过程中选取的是1484 ms时的水平时间切片数据，如图 1所示.

从图中可以看出，在这个切片中可以看出存在多个横向断层，对于不明显的部分及小断层的存在及位置不能确认.利用本文中提到的C1算法、HOS-C1算法和C3算法及本文提出的HOS-C3算法，分别对以上切片进行了断层识别.为了统一比较不同相干算法的性能，所有画图结果的相干范围是，具体结果如图 2-图 5所示.

经过相干算法对原始数据切片的处理，全部相干算法能够将原始切片中明显存在的断层检测出来.经过传统的C1算法、C3算法和HOS-C1算法、HOS-C3算法比较，可以明显看出，经过高阶统计量处理的相干算法结果都有了明显的提升.

经过C1算法和HOS-C1算法结果图的比较，可以看出，C1算法可以检测出断层的位置，包含噪声也被尽可能多的剔除，但是检测精度不高，不能很好的区分出断层和背景信息，不利于后期对断层的精细检测；HOS-C1算法能够尽可能大的区分开了断层和噪声，并且断层能够明显的显示出来，不足之处是背景信息的噪声较大，对后期的处理也有较大的影响.C3算法较C1算法有了很大的改进，能够比较明显的检测出大的断层的存在，不足之处是检测结果存在较大的噪声影响，噪声已经明显影响到了小断层的检测；HOS-C1算法则是在增强断层检测结果的同时也增加了结果中的噪声.经过C3算法和HOS-C3算法结果图的比较，可以看出，C3算法已经在C1算法基础上增强了断层结果的检测，但由于C3算法利用主特征值的方法导致相邻地震道的连续性变得很差，结果变得模糊不清；HOS-C3算法较C3算法更能检测出断层的所在，在大断层以外的小断层也能有较好的检出效果，显示出的断层更清晰，并且更容易识别.HOS-C3算法则是在C3基础上改进而来，既保留了C3断层结果的优良效果，又没有过多的受到噪声和主特征值导致的不连续的影响.HOS-C3算法和HOS-C1算法的结果比较也可以看出，HOS-C3算法比HOS-C1算法的结果也有了很大的提升，断层位置的检测精度和清晰度有了提高，非断层位置的结果也比较清楚.

总体而言，经过高阶统计量处理后的相干算法性能都有了较大的提升，无论是在检测效果还是抗噪声方面.HOS-C3算法的起点是C3算法，将高阶统计量方法良好的抗噪性和C3相干算法的高精度相结合，提取出了更加明显的相干数据切片，取得了很好的效果.

经过对不同相干算法计算过程中复杂度的分析，在程序实现过程中对不同相干算法的运行时间进行比较.不同相干算法的运行时间如表 1所示，表中显示的时间是程序运行时间，包括数据的读取和相干过程的计算.计算过程中，C1算法中的延时分别取为和，C2算法中的窗口长度为7.程序运行时间也显示，HOS-C1算法较C1算法的计算时间明显增加，HOS-C3算法在计算时间上较C3算法的计算时间大约增加了15%.以上程序运行的时间与分析得出的相干算法的时间复杂度是一致的.虽然HOS-C3算法比C3算法运行时间较长，但从取得的效果说明HOS-C3算法是可以接受的.

4 结 论

本文提出了HOS-C3算法，基本原理是将高阶统计量的过程与C3算法结合.通过实验结果的比较，本文提出的HOS-C3算法能够较C1算法、C3算法和HOS-C1算法显示出更清晰的断层，说明提出的方法比这些算法有更好的抗噪性.断层结果同时说明HOS-C3算法在断层的精度上更是比以往的算法有很大的提升，这将为后续的断层识别提供保证.在计算复杂度方面考虑，高阶统计量的方法当和C1算法结合时，具有比C1算法更高的计算复杂度；当高阶统计量方法和C3算法结合时，仅有与C3算法相似的计算复杂度.虽然新的算法在计算量有了一定的增加，但是计算精度和结果的准确性有了提高.后期工作的重点是从相干体中提取断层的展示.

致 谢 感谢教育部新世纪人才支持计划(NCET-13-1011)和国家自然科学基金(61102103，61302138)对本文的全力支持.感谢张世晖老师提供的三维地震数据.