0 引言

利用地震数据识别断层^[1]、裂缝和不整合面等不连续性结构对隐蔽油气藏、地下水勘探等具有重要的意义，这些地质构造的地震响应为反射波同相轴不连续。地震数据中的不连续性识别方法从简单到复杂，经过了几十年的发展，已经逐步趋向成熟。目前，相干体技术、方差体技术等很多方法均能实现地震数据不连续性识别。Bahorich等^[2]于1995年提出了在三维断层解释技术中占据重要地位的基于互相关的第一代相干算法。Marfurt等^[3]在1998年提出了基于多道相似性的第二代相干算法。Gersztenkorn等^[4]在1999年提出了基于特征结构的第三代相干算法。方差体技术^[5]是在误差分析理论的基础上发展起来的，该技术使用局部方差来反映反射振幅信号的横向不连续性，但在倾斜地层或者起伏地层区域不能有效识别地层不连续性。Randen等^[6]在2000年将混沌体技术应用于地震横向不连续性识别，该技术依据地震振幅数据的规律性和混乱性来识别不连续变化。相干值、方差和混沌值都是基于滑动窗口计算的，对处理数据信息中的局部非平稳特征有一定难度^[7]，参数的选取往往会不同程度地影响最终的应用效果，并且对大倾角地层的地震数据应用效果不理想。

为了提高地震数据不连续性识别的辨识度和稳定性，本文将在指纹图像处理中衡量局部方向场信息各向异性强度的度量“一致性”，与基于正交多项式曲面拟合的Facet模型梯度算子进行结合，用于检测地震数据的不连续性。通过与C₃相干算法和方差算法的识别结果进行对比，验证该方法的可行性和有效性。

1 方法原理 1.1 一致性属性

叠后地震数据实际是地下地质构造形态的反射波成像，地震记录中的同相轴类似指纹图像中的脊谷线，所以可以借鉴指纹识别技术对地震数据进行处理。在指纹图像中，灰度沿指纹脊谷线走向变化最慢，而沿指纹脊谷线的垂向方向变化最快，表现为方向各向异性特征^[8]，该特征可用于求取脊谷线走向信息。将衡量各向异性强度的度量称为“一致性”^[8]。对地震数据求取一致性属性的过程如下。

在直角坐标系下，二维数据体中数据点(x, y)的梯度向量表示为

式中：G_x(x, y)、G_y(x, y)分别为点(x, y)沿x、y方向上的梯度值; I(x, y)表示点(x, y)的振幅。将[G_x(x, y), G_y(x, y)]^T简记为[G_x, G_y]^T。具体计算梯度时，本文选用了基于Facet模型的曲面拟合算子，在1.2节作详细介绍。

在极坐标系下，公式(1) 可表示成为

式中：G_ρ为极径；G_φ为极角。则

定义数据点(x, y)的平方梯度向量表达式为

提取一个以数据点(x, y)为中心、包含n×n个数据点的窗口B，将该数据块的平方梯度向量表示为

其中，

公式(5) 计算的目的是将目标点局部邻域范围内与同相轴走向垂直的平方梯度向量求和，以作为目标点的梯度向量输出，从而降低噪声影响，进而获取更为准确的各向异性特征信息，以提高求取地层走向信息的精度。

地震同相轴与指纹脊谷线相似，具有两侧边界线。理想情况下，若按公式(1) 计算，则会在两侧边界线上出现大小相等、指向相反的梯度向量。若在局部邻域范围内对各梯度向量求和，指向相反的向量就会相互抵消，由此不能得到有用的各向异性特征信息，所以定义如公式(4) 所示的平方梯度向量。平方梯度向量的角度是梯度向量的两倍，这样指向相反的向量变成了相同指向的向量，按公式(5) 求和，各向异性特征更加突出。

为了衡量数据块B各向异性的强度，Kass等^[8]提出一致性概念：

定义Coh_B为数据块B中心数据点(x, y)的一致性，0≤Coh_B≤1。

如果数据块B内的地层连续，则相应各点平方梯度向量的指向就会相同，各向异性指示性强，此时，所有向量模的和等于所有向量和的模，Coh_B等于1；如果数据块B内地层存在横向不连续或强随机噪声，那么各点平方梯度向量的指向不尽相同，各向异性特征不明显，此时，所有向量模的和远大于所有向量和的模，Coh_B接近于0或等于0；其他情况，Coh_B介于0和1之间。由此可见，在地层不连续范围内，可利用相邻地震道之间一致性值的突变表征地层的不连续性。

对于进行边缘保护滤波后的地震数据，仍会存在一定的残余噪声，而在残余噪声区及弱反射地震信号区的一致性值也低于1，这会干扰解释人员对横向不连续性的识别。经过大量实验发现，由残余噪声或弱反射地震信号产生的低一致性值与地层不连续位置处的近0值差值较大。依据此数据特性，本文对一致性数据Coh_B作进一步阈值化处理，以压制伪不连续性信息的干扰：

式中，T是基于一致性数据Coh_B计算得到的阈值。公式(7) 的计算，在一定程度上消除了因噪声或弱反射地震信号造成的伪地层不连续性干扰，使不连续性信息更加直观、清晰。其中，阈值T的选取要依据不同数据而定。T太大，对背景干扰的压制作用不明显；T太小，会容易将有效表征地层构造不连续特征的信息去掉。本文选取Tian等^[9]提出的自适应阈值选取方法计算T值。该方法计算的阈值可以在消除伪地层不连续性信息的同时，最大程度地保留真实的地层不连续性信息。

为了能更精确地检测到地层不连续的具体位置，也为了能进一步加强不连续性特征信息的连通性，本文在公式(7) 的基础上，再利用数学形态学方法中的腐蚀、膨胀、细化3种算法作进一步处理，便可以得到能精确指示地层构造不连续性的线条。

以上过程便是本文提出的利用一致性检测地层不连续性的过程，并将最终的检测结果定义为Coh。

1.2 Facet模型曲面拟合算子

Haralick等^[10]在1984年首次提出了基于多项式曲面拟合的边缘检测算子，称为Facet模型梯度算子。Facet模型将二维灰度图像转化为空间的曲面，为图像分析提供了更加丰富而准确的信息，改善了传统梯度算子的准确性和精度^[11]。

Facet模型梯度算子的基本思想：对目标点局部邻域的数值进行最小二乘拟合，然后对拟合函数进行一阶求导，便可得到梯度算子，这样做可以在提高边缘定位精度的同时实现对噪声的压制^[12]，使梯度算子更易进行扩展。姜杭毅等^[13]对该方法进行了详细的数学推导，Facet模型梯度算子的计算过程如下。

设S是以目标点(x₀, y₀)为坐标原点、大小为R×R(R为奇数)的二维局部对称邻域，F(x, y)为该局部邻域内某点的数值，利用一个二元三次多项式f(x, y)对F(x, y)进行拟合：

式中，k_i为拟合系数(i=1, 2, …, 10)。在邻域S内，二元三次多项式函数的一组二维离散正交多项式基函数系为

式中，

将拟合函数f(x, y)用上述离散正交多项式组表示成如下形式：

式中，K_i为正交多项式组的拟合系数。对比公式(8) 和公式(9) 可见，k_i可用K_i线性表示。利用正交多项式对S邻域内的数据F(x, y)进行最小二乘曲面拟合，K_i值应满足公式(10) 取值最小：

根据正交多项式的性质，有

通过推导可求得K_i为

可见，每个K_i都可以单独由F(x, y)的线性组合表示，其对应的权值模板W_i为

公式(8) 在目标点(x₀, y₀)处的一阶偏导数为

因此，若求目标点(x₀, y₀)的梯度，只需要求得k₂和k₃。因为k_i可用K_i线性表示，那么k_i的权值模板ω_i也可以由K_i的权值模板W_i线性表示。由公式(14) 可知，k₂和k₃的权值模板ω₂和ω₃可进一步表示为ω_x和ω_y，则

至此，曲面拟合方法求取梯度最终转化成利用两个垂直方向的权值模板ω_x和ω_y对原始图像的局部邻域进行模板卷积的过程。

卷积模板的大小会直接影响地层不连续位置的定位精度及对噪声的压制效果。卷积模板越大，对噪声的压制效果越强，但所检测到地层不连续位置宽度也越大，同时会降低地层不连续特征信息的连通性；卷积模板越小，对地层不连续位置的精度定位越准确，但同时对噪声的压制效果会变差。所以，针对高质量地震数据，可以选用较小的卷积模板，而对信噪比低的地震数据则可选用较大的卷积模板。

针对本文地震数据，在求取一致性的过程中，Facet模型梯度算子的卷积模板选用5×5大小。

2 理论模型试算

为了验证本文方法对地震数据不连续性识别的可靠性，选取理论合成地震记录“sigmoid”模型^[14]进行测试(图 1a)。该模型包含顶部的单向倾斜地层、中部的向斜构造、底部的背斜构造、地层间的角度不整合接触和一个正断层。该模型的时间深度为200个采样点，采样间隔为4 ms；测线长度为200道，道间距为0.01 km。

用基于Facet模型梯度算子的一致性算法对图 1a进行地震数据不连续性信息识别，其中参数n=7、T=0.89，结果如图 1b所示。可以看出：在地层连续处，一致性数值接近1；而在断层处及地层不整合接触位置，一致性值接近于0，从而突出地层的横向不连续性。该结果证明了本文方法用于地震数据不连续性识别的有效性。为了验证本文方法的优越性，选取传统的C₃相干算法(图 1c)和传统的方差算法(图 1d)对图 1a进行地震数据不连续性识别。由图 1c、d可见，因为褶皱构造的存在，地层倾角在短距离范围内变化较大，而这两种算法对地层倾角很敏感，虽然能识别出地层不整合接触的位置，但是并不能清晰地识别出断层信息。为了进一步比较，分别对这两种算法进行倾角调向^[15]，识别结果如图 1e、f所示。由图 1e、f可见，考虑倾角的C₃相干算法和方差算法虽能在一定程度上压制陡倾地层造成的背景干扰，但识别出的地层不连续信息仍旧比较模糊，辨识度较差。本文提出的一致性算法能很精确地指示出地层不连续的具体位置，并且表征不连续性线条的连通性要高于前面两种算法。

本文提出的一致性算法、C₃相干算法和方差算法的计算结果均在0和1之间，具有相同的指示地震数据不连续性的数据性质；因此，分别从图 1b、e、f中抽取各地震道在同一时间深度(0.352 s)所对应的检测结果值(图 2)，经分析可见：在同一地层不连续处，考虑倾角的C₃相干算法和方差算法的计算结果值较大，与相邻连续区域的计算结果值区别不大，不连续性不利于被识别；而本文方法计算的一致性数据接近于0值，与连续区域的1值差别很大，故在相同色标显示下，能更清晰地指示出地层的构造不连续位置。

为了进一步验证本文方法的优越性，分别利用本文所提的一致性算法、基于倾角调向的C₃相干算法和基于倾角调向的方差算法对三维地震数据体^[16](图 3a)沿其横测线方向进行地层横向不连续性识别。该三维地震数据体包含顶部的水平地层、中部的陡倾地层、底部的倾斜地层以及不同类型地层的不整合接触关系和一个地堑式构造。图 3a中的3个侧面分别代表三维数据体内部标识线位置处的切片(前侧面为横测线0.45 km位置的纵测线剖面，右侧面为纵测线0.95 km位置的横测线剖面，顶侧面为3.3 s位置处的时间切片)。从断层检测结果(图 3b、c、d)对比来看：本文所提的一致性算法不仅能很清晰、精细地指示出断层的位置，还能较完整且准确地指示出地层角度不整合接触的位置，最主要的是该算法不受地层倾角的影响，背景很干净，没有多余的干扰存在；然而，基于倾角调向的C₃相干算法和基于倾角调向的方差算法均不能完全消除掉陡倾地层所产生的背景干扰，并且这两种方法均不能清晰地指示出地层角度不整合接触的具体位置。

综上可知，C₃相干算法和方差算法对陡倾地层非常敏感，计算结果明显受陡倾地层的干扰，将这两种算法经倾角调向后的计算结果也不能完全消除掉地层倾角所产生的低值干扰，而且该两种方法并不能检测出由角度不整合接触所造成的地层横向不连续性。与这两种方法对比，本文所提方法面对陡倾地层无需倾角调向环节，计算结果丝毫不会受地层倾角的影响，而且该方法还能精确、完整地指示出断层、角度不整合接触等地层不连续的位置。图 2和图 3的对比结果验证了本文提出的一致性算法检测地层不连续性的可行性和优越性。

3 实际数据分析

选取墨西哥湾某地区的二维时间偏移剖面对本文方法进行测试^[14]，该数据已利用Liu等^[17]的滤波方法进行了去噪(图 4a)。该地震数据的时间深度范围为0.40~1.48 s，有271个采样点，采样间隔为4 ms，测线长度为125道。用基于Facet模型梯度算子的一致性算法对图 4a进行地震数据不连续性信息的识别(图 4b)，参数n=9、T=0.81。为了与其对比，同时给出基于倾角调向的C₃相干算法(图 4c)和基于倾角调向的方差算法(图 4d)的识别结果。经过对比分析可见：后两种方法无法完全消除弱反射地震信号及残余噪声造成的背景干扰，识别到的不连续性位置比较模糊；而本文方法对不连续性的辨识度更高，能更完整地展示不连续性的垂向延伸长度，对伪不连续性信息的压制效果最佳。

同样，为了进一步测试本文所提方法在三维实际数据中应用的可行性，分别利用本文所提一致性算法、基于倾角调向的C₃相干算法和基于倾角调向的方差算法对墨西哥湾某地区三维叠后数据体^[18](图 5a)沿横测线方向进行处理，3种算法的识别结果分别如图 5b、c、d所示。通过对比可见，本文所提算法抗噪性最高，而且其所识别到的断层线在连续性和纵向延伸完整度方面及对断层位置的指示精度方面要明显优于其他两种算法。

综上，无论是对二维实际地震数据还是三维实际地震数据，本文方法均能有效识别其中的不连续性信息，且识别效果要明显优于其他两种算法，具有一定的实际应用价值。

4 结论

1) 为了实现对地震数据不连续性的识别，本文将一致性作为不连续性识别的新方法应用于地震数据处理，利用基于Facet模型梯度算子计算一致性，并对一致性数据作阈值化及数学形态学处理，实现了地震数据不连续性的自动识别。

2) 利用基于Facet模型梯度算子的一致性、基于倾角调向的C₃相干算法和基于倾角调向的方差算法分别对理论模型和实际地震数据进行处理，结果证明基于Facet模型梯度算子的一致性在地震数据不连续性识别方面不受地层倾角影响，无需倾角调向环节，具有更高的辨识度，能更好地压制伪不连续性信息。因此，该方法为地震数据不连续性识别提供了一种有效的可选工具。