三维边缘保持滤波方法在海上地震数据噪声压制中的应用研究——以东海某凹陷为例

孙永壮，李键，秦德文，刘庆文

(中海石油(中国)有限公司 上海分公司研究院，上海 200335)

0 引言

复杂断块区域的构造解释一直是海上地震勘探的难题，受相关地质因素以及地震采集处理技术的影响，该区域地震资料出现信噪比低、断面成像不清等特征，给后续的构造解释以及断块圈闭的识别带来较大困难。因此，提高复杂断块区域的地震资料品质是开展精确构造解释的必由之路。国内外采用多种压制随机噪声用以提高资料的信噪比[1]。目前常用的随机噪声压制方法包括F-X域反褶积[2]、小波变换[3-5]、改进型子波域迭代去噪[6]、聚束滤波[7]、多项式拟合去噪[8]、曲波变换去噪[9]、集成经验模态分解降噪[10]、字典算法去噪[11]、双重稀疏表示法去噪[12]、奇异值分解去噪[13-14]和多手段组合综合去噪[15]等方法。然而上述方法在三维复杂断块区域压制随机噪声时，容易对不连续地质体如断层、岩性尖灭等有效信号造成损伤，致使去噪后数据体水平切片的连续性较差，很难获取理想的去噪效果。

为此，基于各种约束条件的滤波方法被应用于三维地震数据的噪声压制，张恒磊等[16]在反射波各向异性特征和随机噪声各向同性的基础上，提出利用地震资料各向异性特征识别有效反射进行随机噪声的压制，在压制随机噪声的同时保护弱反射波。刘洋等[17]、尹川等[18]在地层倾角等构造条件的约束下进行随机噪声的压制，在保持边界信息的基础上取得较好的噪声压制效果。姚振岸等[19]将曲波变换和各向异性扩散滤波结合起来，提出一种边界和振幅保持的联合去噪方法。然而，上述方法在使用过程中，对于三维数据体的滤波单元窗口形状较难选择，过大或者过小都容易使地震波形变形、同相轴锐化、在断层等边缘处的有效信息容易出现斑块现象。而蔡涵鹏等[20]提出一种基于多窗口相干性的倾角导向主分量滤波方法，结合边缘保持滤波，实现多窗口边界保持和较好的滤波效果。但是，在三维资料中是按照二维时间切片进行滤波，在三维复杂断块区域的滤波效果有待提高，因此，在地质构造约束下选择或建立一种适应性强的滤波算子是一种可选方案。

基于线性相干信号可预测性，F-X滤波方法可以有效预测出相干信号，分离噪声并实现噪声去除的目的。但是常规的F-X滤波方法存在不连续地质体有效反射信息的模糊化缺陷。针对上述问题，本文研究分析F-X滤波和各向异性扩散滤波的特点，综合利用F-X滤波和各向异性扩散滤波的优势，提出一种基于地质构造约束的三维自适应F-X边缘保持滤波方法，新方法在局部地层倾角和方位角的约束下，沿着地层倾角和方位角对三维目标地质体数据提取滤波窗口，在滤波窗口内进行自适应F-X滤波。采用三维模型数据和实际数据进行检验，和常规F-X滤波方法以及各向异性扩散滤波方法相比，新方法在边缘保持和噪声压制方面具有明显的优势和效果。

1 方法原理

各向异性扩散滤波，其优势在于扩散滤波方向沿着构造方向进行，在局部地层结构变化剧烈的地层停止扩散滤波。上述滤波方式能够有效保护不连续有效信号，但是迭代次数对于最后的滤波结果影响较大，迭代次数少难取得较好的去噪效果，迭代次数的增加会给断层等边缘信息带来模糊效应。但是该方法建立的梯度结构张量很好地描述了地质结构特征[21-22]，能够获取地震结构信息，有效识别复杂构造范围，提高地层走向以及倾角估计的空间一致性。常规F-X保边滤波技术虽然能够很好地保持边缘信息，且不受滤波窗口大小影响，但是没有充分考虑目标位置处三维空间的地质构造信息，得到的滤波效果也难以满足后续工作需求。

为综合利用F-X滤波和各向异性张量扩散滤波的优势，本文提出了基于地质构造约束的三维自适应F-X滤波方法，该方法以地质构造导向因素中的地层倾角和方位角为约束，沿着局部地层的倾角和方位角提取沿层滤波窗口，在滤波窗口中进行自适应F-X滤波。

目前对于计算数据倾角方位角算法较多，主要包括离散扫描算法[23]、复数道分析方法[24]和梯度结构张量算法(GST)[25]等。本文采用GST算法计算地层倾角和方位角，该方法采用梯度张量矩阵来描述地震数据的构造特征，能够有效度量地层产状以及断层、砂体、河道等特殊构造的内部纹理特征。

1.1 沿地层倾角方位角的滤波多窗口提取

从P-M扩散方程出发[26]，建立三维局部梯度结构张量J：

(1)

式中：Uσ表示需要滤波的图像。结构矩阵J为半正定矩阵，其特征值为不小于零的实数，通过矩阵分解，结构矩阵变为：

(2)

求解式(2)可得其特征值λ以及对应的特征向量v，当3个特征值λ表现不同的数值关系时，对应的特征向量v表示地震资料不同的地质意义，具体如表1所示。求解式(2)之后，得到局部地层法线方向v1，具体如图1所示。将其在x、y、z方向的分量求解出来：

则局部地层倾角dip为：

(3)

表1 特征值、特征向量对应的地质结构Table 1 The geological structure table corresponding to eigenvalues and eigenvectors

方位角azimuth为：

azimuth=atan2(v1y,v1x),

(4)

atan2是一个函数，用于计算方位角，根据倾角方位角，计算x、y方向的延迟时间，进而插值出x、y方向邻道数值，用以提取目标点位置处的局部沿层滤波窗口，进行自适应F-X保边滤波。x、y方向的延迟时间计算公式如下：

(5)

选定滤波窗口之后，在滤波窗口内划分多边形，多边形的划分如图2所示。

图1 根据结构张量矩阵估计得到局部地层的法线方向v1Fig.1 Estimate the normal direction v1 of the local strata according to the structure tensor matrix

a—六边形滤波窗口;b—五边形滤波窗口;c—正方形滤波窗口a—hexagonal filter window;b—pentagonal filter window;c—square filter window图2 多窗口滤波窗口示意Fig.2 Schematic diagram of multi-window filter window

如图2所示，在包含局部地层中心点(绿色实心圆点)的邻域内(图1所示)采用9个小多边形选中数据，然后计算各多边形内的方差，选择方差最小的多边形作为最终滤波窗口，在该窗口内执行自适应F-X保边滤波，将滤波结果作为中心点输出结果。

1.2 三维自适应F-X保边滤波算法

常规F-X滤波通常在二维剖面中沿着时间轴提取一维滤波窗口，然后在一维滤波窗口内执行双向滤波，将正反向滤波结果的均值作为输出结果[27]。其公式如下：

F=0.5F1+0.5F2,

(6)

式中：F1、F2分别表示正向滤波和反向滤波结果，F表示当前窗口下最终的滤波输出结果。由于不同方向的滤波对于地震有效信号的影响不同，常规的F-X方法容易损伤断层、砂体边界等边缘不连续反射信息。为进一步提升地震数据的信噪比，从三维地震数据中提取二维沿层滤波窗口，具体如图2所示。然后在二维窗口内进行不同方向的F-X滤波，从而有效提升F-X的滤波效果。因此，本文对双向滤波结果的合并方式进行改进，利用地层信息自适应设置正反向的滤波系数，则式(6)变为：

(7)

式中：α为自适应权系数，其取值范围为0～1；F1i、F2i为每道正向滤波和反向滤波结果，F1、F2表示图2所示的沿地层倾角方位角提取的多窗口中各道正向滤波之和和反向滤波之和，其中自适应权重的计算方式如下：

(8)

根据图2中提取的多窗口，对窗口内的地震数据进行常规的F-X正、反向滤波，根据滤波结果确定权重系数α。上述权重确定方式是数据驱动，不依赖任何模型分布，具有自适应性，因此，更加符合地震数据地质特征。

2 模型测试

为测试算法的有效性，建立三维断层模型，断层的走向成三角函数关系变化，然后在原始无噪模型中添加随机噪声(5 dB)，原始无噪模型和含噪模型分别如图3a、b以及图4所示。采用常规多窗口滤波方法、各向异性扩散滤波方法和本文方法分别对含噪模型进行去噪，去噪结果如图5、图6和图7所示。对比3种去噪结果，可以看到常规方法去噪的结果在断面下部的地层同相轴边缘位置出现台阶状模糊效应。对应的剖面和水平切片上，不仅地层边缘有所损伤，断面边缘也有一定的斑状模糊效应，具体如图5中线圈及箭头位置处所示。各向异性扩散滤波方法能够较好压制随机噪声，且能够较好地保护边缘，但需多次试验确定合适的参数，噪声压制和边缘保持之间很难同时取得较好效果，边缘保持较好的同时往往噪声压制不够理想，具体如图6所示。而本文的去噪方法不仅能有效压制随机噪声，同时对于地层以及断面的有效反射信息保护较好，具体如图7所示。同时，抽取XLINE剖面中的第55道数据的局部进行对比分析，具体如图8所示。可以看到：各向异性(绿线)可以很好地消除噪声，但由于迭代次数的选择，导致边缘信息有所损伤；常规多窗口滤波(浅蓝细线)噪声压制存在一定的不足，且边缘信息也有损伤；而本文新方法不仅压噪效果好，同时对于边缘信息保护更加有利。

a—原始无噪三维模型;b—含5 dB噪声的三维模型a—raw model without noise;b—the model with noise of 5 dB图3 原始三维断层模型及含噪三维断层模型Fig.3 Three-dimensional fault model with and without noise

a—原始无噪模型主测线方向;b—原始无噪模型横测线方向;c—原始无噪模型水平切片;d—噪声模型主测线方向;e—噪声模型横测线方向;f—噪声模型水平切片a—raw model inline;b—raw model xline;c—raw model time slice;d—model with noise inline;e—model with noise xline;f—model with noise time slice图4 原始模型及噪声模型切片Fig.4 Slices of raw model and model with noise

a—去噪数据主测线方向;b—去噪数据横测线方向;c—去噪数据水平切片a—de-noising model inline;b—de-noising model xline;c—de-noising model time slice图5 常规F-X方法去噪结果Fig.5 De-noising result by conventional F-X method

a—去噪数据主测线方向;b—去噪数据横测线方向;c—去噪数据水平切片a—de-noising model inline;b—de-noising model xline;c—de-noising model time slice图6 各向异性扩散滤波去噪结果Fig.6 De-noising result by anisotropic diffusion filter method

a—去噪数据主测线方向;b—去噪数据横测线方向;c—去噪数据水平切片a—de-noising model inline;b—de-noising model xline;c—de-noising model time slice图7 本文方法去噪结果Fig.7 De-noising result of the method proposed in this paper

1—无噪声数据；2—本文方法去噪数据；3—各向异性去噪数据；4—噪声信号；5—常规F-X去噪数据1—raw data;2—de-noising result by our new approach;3—de-noising result by the anisotropic approach;4—noise data;5—de-noising result by traditional F-X approach图8 原始无噪数据和含噪数据以及3种去噪方法结果单道数据对比显示Fig.8 Original data,noisy data and the results of three de-noising methods in single trace comparison

3 实际数据验证

采用东海某工区的实际三维地震资料进行验证。工区位于凹陷西部斜坡带，断层呈顺向断阶展布，对整个凹陷的构造、沉积、成藏起主要的控制作用。在张扭性应力作用下，地层倾角变大，方位角多变，地层倾角方位角如图9b、c所示，受复杂断块的影响，现有地震资料信噪比低，随机噪声严重，严重制约后续的解释和反演工作，如图9a所示。为此，对该工区的数据，利用倾角、方位角信息作为约束引导，开展了基于地质构造导向的多窗口三维滤波方法处理，同时和常规多窗口滤波方法，各向异性滤波方法进行对比。滤波结果如图10、图11所示。抽取时间切片的某条主测线的数据进行对比分析，具体如图12所示。

从图中可以发现，通过3种方法的滤波处理，原始地震剖面的随机噪声得到很好的压制。对于常规F-X滤波结果而言，随机噪声得到有效压制，地层同相轴连续性增强，但是断层边缘变得模糊，具体如图10b。各向异性滤波结果相对常规F-X滤波结果，其噪声压制更好，但是对于断层边缘存在一定程度的平滑。本文方法的滤波结果中，不仅随机噪声得到压制，在增强地层反射同相轴连续性的同时，断层的边缘信息也得到有效保护，如图10a所示。分别提取3种滤波方法的水平切片，可以看出，各向异性滤波结果的水平切片平滑性和连续性都较好，常规F-X去噪方法在压制噪声的同时，使水平切片的平滑性和连续性较差，如图11b。进一步分析水平切片的单道记录如图12所示，3种方法的滤波结果都较平滑，但是各向异性滤波结果(绿线)对边缘平滑更大，常规F-X滤波结果(蓝线)的局部仍然存在信号跳跃，相对而言，新方法滤波结果(红线)信号平滑好，边缘保护更好。

a—原始地震数据;b—计算的构造倾角;c—计算的构造方位角a—raw seismic data;b—structural dip;c—structural azimuth图9 原始地震数据及对应的构造倾角、方位角剖面Fig.9 Raw seismic data and corresponding structural dip and azimuth profiles

a—新方法去噪结果剖面;b—常规F-X方法去噪结果剖面;c—各向异性扩散滤波结果剖面a—de-noising result with new filtering approach;b—de-noising result with traditional F-X filtering approach;c—de-noising result with anisotropic diffusion filtering approach图10 不同的去噪结果剖面Fig.10 Profiles of different filtering approaches

a—新方法去噪结果水平切片;b—常规F-X方法去噪结果水平切片;c—各向异性扩散滤波结果水平切片a—de-noising time slice with new filtering approach;b—de-noising time slice with traditional F-X filtering approach;c—de-noising time slice with anisotropic diffusion filtering approach图11 不同去噪结果水平切片Fig.11 Time slices of different method

1—原始数据；2—本文方法滤波结果；3—各向异性滤波结果；4—常规F-X滤波结果1—raw data;2—de-noising result by our new approach;3—de-noising result by the anisotropic approach;4—de-noising result by traditional F-X approach图12 数据压噪前后对比示意Fig.12 De-noising result before and after using different filter approaches

压制的噪声剖面和压制的噪声时间切片中，常规F-X滤波对地震反射的有效同相轴有所损伤，具体如图13b和图14b。各向异性滤波结果相对而言，能够压制部分噪声，对边缘保护较好，但是仍然对有效反射信号有所损伤，具体如图13c和14c所示。而本文的方法在压制随机噪声的同时，水平切片同相轴连续性和平滑性得到较好的提高，如图11a。压制的噪声剖面和切片中，基本不包含地震有效反射，对有效信号的保护更好，具体如图13a和14a所示。

a—新方法压制的噪声剖面;b—常规F-X方法压制的噪声剖面;c—各向异性扩散滤波压制的噪声剖面a—de-noising error with new filtering approach;b—de-noising error with traditional F-X filtering approach;c—de-noising error with anisotropic diffusion filtering approach图13 不同去噪方法压制的噪声剖面Fig.13 Profiles of de-noising error with different filtering approaches

a—新方法压制的噪声水平切片;b—常规F-X方法压制的噪声水平切片;c—各向异性扩散滤波压制的噪声水平切片a—time slice of de-noising error with new filtering approach;b—time slice of de-noising error with traditional F-X filtering approach;c—time slice of de-noising error with anisotropic diffusion filtering approach图14 不同去噪方法压制的噪声水平切片Fig.14 Time slices of de-noising error with different filtering approaches

4 结论

本文提出的基于地质构造约束的三维自适应F-X滤波方法，采用地质构造因素中的地层倾角和方位角为约束，实现边界保持滤波。实际资料表明，该方法能够有效提高地震数据的信噪比，有助于地层的边界保护及增强同相轴连续性，而有效的不连续信号(断层等)也得到增强，更加清晰反映地质体的反射结构和分布特征。与常规F-X滤波以及各向异性扩散滤波方法相比，本文方法具有以下特点：

1)相对于常规F-X滤波算法而言，新方法选取目标点周围的三维空间窗口，充分利用目标点的空间地质信息，但是计算量较大，计算效率主要和扩散张量的构建以及多窗口参数的设置相关。

2)与常规的各向异性扩散滤波算法相比，新方法提取目标点的滤波多窗口时构建了扩散张量，排除了迭代次数对边缘信息的破坏。而基于地质构造约束的三维F-X随机噪声压制方法可以有效压制随机噪声，在提高复杂断块区域地震资料信噪比的同时，有效保护断层等不连续地质体的有效反射信号。