地震数据中除有效反射波之外的所有信息都被视作为噪声或者干扰波, 根据噪声信号的不同特征, 可以将噪声分为相干噪声和随机噪声两类^[1]。相干噪声主要包括多次波、面波、海上鸣震、折射波等; 随机噪声主要包括各种环境因素和人为因素引起的噪声^[2]。相干噪声在时空域具有一定的规律性, 在频谱、视速度以及到时等方面与有效信号存在差异; 而随机噪声则缺乏规律性, 在各种域中与有效信号混杂在一起, 存在于整个记录之中, 较难去除^[3]。因此, 要想获取高信噪比的地震资料, 必须进行噪声压制。

目前, 消除地震记录中的随机噪声的方法主要分为时空域方法和变换域方法。时空域方法主要包括多项式拟合滤波^[4-5]、中值滤波^[6]、奇异值分解滤波^[7]、独立分量分析滤波^[8]、双曲滤波^[9]和时频分析法^[10]等; 变换域方法主要包括K-L变换滤波^[11]、F-X预测滤波^[12]、小波变换滤波^[13]和曲波变换滤波^[14]等。上述去噪方法都已得到广泛应用且取得了较好的去噪效果, 但受各自假设条件和方法原理的限制, 对于地震数据中有效信号不连续的边缘(断层)以及突变点(地层或者岩性尖灭点)等位置处的噪声压制效果不甚理想。常规的构造导向滤波方法考虑的是单一窗口内的数据信息, 且窗口大小很难选定合适的范围^[15], 因此, 在地层倾角方位角约束下的去噪效果难以满足需求。而其它常规滤波方法如F-X滤波方法在选取滤波窗口时, 没有很好地考虑地层结构因素, 对地层不连续地质体有效信号的损伤较大。双边滤波方法由于具有简单、保护边缘信息的特点得到了广泛应用, 国内外学者对其进行了大量的研究^[16], 提出了自适应双边滤波方法以及改进型的双边滤波方法。然而由于核函数的建立方式对数据的信噪比和边缘信息较为敏感, 滤波窗口的选取也较为困难, 过大或者过小都容易损伤有效信号, 因而限制其在保持边缘信息情况下压制随机噪声的效果^[17-18], 在不连续地质体(断层、裂缝等)较发育的区域, 上述方法很难取得理想的噪声压制效果。

针对上述情况, 本文采用基于多窗口自适应双边滤波方法对地震信号中的复杂断块区域噪声进行压制。在自适应双边滤波方法基础上, 沿局部地层倾角方位角提取滤波窗口, 进行多窗口边缘保持滤波。首先介绍了计算地质体倾角方位角信息的方法原理, 然后根据地质信息提取沿层多滤波窗口, 再介绍在选定的滤波窗口中执行自适应双边滤波的过程。利用建立的三维模型数据验证了该方法的有效性, 最后在南海L盆地某工区进行了实际应用。

1 方法原理 1.1 局部地层倾角方位角计算

基于多窗口自适应双边滤波方法(简称本文方法), 重点在于沿局部地层倾角方位角提取滤波窗口, 因此, 该方法的第一步在于计算局部地层的法线方向。从梯度结构张量出发^[19], 分别计算三维地震数据3个方向的梯度, 则梯度结构张量T的表达式为:

$ \mathit{\boldsymbol{T}} = \mathit{\boldsymbol{g}}{\mathit{\boldsymbol{g}}^{\rm{T}}} $

(1)

其中,

$ \mathit{\boldsymbol{g}} = \left[ {\begin{array}{*{20}{c}} {\frac{{\partial u}}{{\partial x}}}\\ {\frac{{\partial u}}{{\partial y}}}\\ {\frac{{\partial u}}{{\partial z}}} \end{array}} \right] $

(2)

∂u/∂x, ∂u∂y, ∂u/∂z分别表示地震数据u在3个方向的梯度, 分别记为g₁, g₂, g₃, 则平滑张量T_S为:

$ {\mathit{\boldsymbol{T}}_{\rm{S}}} = \left[ {\begin{array}{*{20}{c}} {G \cdot {\mathit{\boldsymbol{g}}_1}{\mathit{\boldsymbol{g}}_1}}&{G \cdot {\mathit{\boldsymbol{g}}_1}{\mathit{\boldsymbol{g}}_2}}&{G \cdot {\mathit{\boldsymbol{g}}_1}{\mathit{\boldsymbol{g}}_3}}\\ {G \cdot {\mathit{\boldsymbol{g}}_2}{\mathit{\boldsymbol{g}}_1}}&{G \cdot {\mathit{\boldsymbol{g}}_2}{\mathit{\boldsymbol{g}}_2}}&{G \cdot {\mathit{\boldsymbol{g}}_2}{\mathit{\boldsymbol{g}}_3}}\\ {G \cdot {\mathit{\boldsymbol{g}}_3}{\mathit{\boldsymbol{g}}_1}}&{G \cdot {\mathit{\boldsymbol{g}}_3}{\mathit{\boldsymbol{g}}_2}}&{G \cdot {\mathit{\boldsymbol{g}}_3}{\mathit{\boldsymbol{g}}_3}} \end{array}} \right] $

(3)

其中, 高斯滤波器G的表达式为:

$ G\left( {x,y,z;\sigma } \right) = \exp \left[ {\left( {{x^2} + {y^2} + {z^2}} \right)/\left( {2{\sigma ^2}} \right)} \right] $

(4)

式中:σ表示滤波尺度因子。

对平滑张量进行特征值分解:

$ {\mathit{\boldsymbol{T}}_{\rm{S}}} = \mathit{\boldsymbol{ \boldsymbol{\varLambda} }}{\mathit{\boldsymbol{v}}^{\rm{T}}} $

(5)

式中: Λ表示对角矩阵, 对角线上的值为特征值, 由大到小分别排列为λ₁, λ₂, λ₃; v表示特征向量矩阵, 列向量表示以上特征值对应的特征向量, 分别为v₁, v₂, v₃。局部地层的法线方向就是最大特征值λ₁对应的特征向量v₁。

地层倾角θ和方位角φ表达式如下:

$ \theta = \left\{ \begin{array}{l} 0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{v_{1x}} = 0\;且\;{v_{1y}} = 0\\ \frac{{\rm{ \mathsf{ π} }}}{2} - \left| {\arctan \frac{{{v_{1z}}}}{{\sqrt {v_{1x}^2 + v_{1y}^2} }}} \right|\;\;\;{v_{1x}} \ne 0\;或\;{v_{1y}} \ne 0 \end{array} \right. $

(6)

$ \varphi = \left\{ \begin{array}{l} \arctan \frac{{{v_{1x}}}}{{{v_{1y}}}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{v_{1x}} \ne 0\;且\;{v_{1z}}{v_{1y}} < 0\\ {\rm{ \mathsf{ π} }}{\mathop{\rm sgn}} {v_{1x}} + \arctan \frac{{{v_{1x}}}}{{{v_{1y}}}}\;\;\;\;\;{v_{1x}} \ne 0\;且\;{v_{1z}}{v_{1y}} > 0\\ \frac{{\rm{ \mathsf{ π} }}}{2}{\mathop{\rm sgn}} {v_{1x}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{v_{1x}} \ne 0\;且\;{v_{1y}} = 0\\ \frac{{\rm{ \mathsf{ π} }}}{2} + \frac{{\rm{ \mathsf{ π} }}}{2}{\mathop{\rm sgn}} {v_{1y}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{v_{1x}} = 0 \end{array} \right. $

(7)

式中:v_1x, v_1y, v_1z分别表示v₁在x, y, z方向的分量。地层倾角和方位角计算完毕之后, 按公式(8)和公式(9)分别估算x, y方向的延迟时间T_x, T_y。

$ {T_x} = \tan \theta \sin \varphi $

(8)

$ {T_y} = \tan \theta \cos \varphi $

(9)

根据两个方向的延迟时间来获取每一个目标点位置处的局部沿层多滤波窗口, 然后执行多窗口自适应双边滤波。

1.2 沿层提取多滤波窗口

为了进一步压制地震数据中的随机噪声, 提高信噪比, LUO等^[20]提出一种多窗口保边滤波方法。该方法在每一个目标点位置处设置一个滑动窗口, 以寻找所有包含该目标点的窗口。然后计算每个窗口内所有样点信息的标准方差, 计算方差最小窗口内的所有样点的平均值, 并将该平均值作为目标点的输出结果^[21], 能够取得较好的去噪效果。为利用多窗口保边滤波的优势, 我们将多窗口滤波的多窗口提取思路引入本文滤波技术流程。具体为:根据上一步计算的局部地层的倾角、方位角, 沿层提取目标点周围多窗口, 该多窗口中包含不同空间位置处的数据, 多窗口的提取示意如图 1所示。图 1中, 红色区域表示提取的多窗口范围, 红色区域上下虚线框表示具有倾角θ和方位角φ的地层, V表示地层的法线方向向量。然后计算目标点周围多窗口的标准方差, 找到方差最小的窗口, 在该窗口内进行双边滤波。

多窗口的设计参考WANG等^[22]的思想, 采用窗口多维划分方法, 以二维数据的总窗口5×5为例, 划分为9个小多边形, 分为4个五边形(蓝色虚线框), 4个六边形(绿色实线框)和1个矩形(红色实线框), 具体如图 2所示。

1.3 双边滤波方法

双边滤波方法最初用于图像信息处理, 它结合给定的窗口中数据信息对噪声进行处理, 其公式如下^[23]:

$ N\left( {x,y} \right) = \frac{1}{{{w_{\rm{p}}}}}\sum\limits_{{{\rm{p}}_{i,j}} \in \mathit{\Omega }} {{w_{\rm{s}}}\left( {i,j} \right){w_{\rm{r}}}\left( {i,j} \right)I\left( {i,j} \right)} $

(10)

$ {w_{\rm{s}}}\left( {\xi ,c} \right) = \exp \left[ { - \frac{{{d^2}\left( {\xi - c} \right)}}{{2\sigma _{\rm{s}}^2}}} \right] $

(11)

$ {w_{\rm{r}}}\left( {\xi ,c} \right) = \exp \left[ { - \frac{{{\delta ^2}\left( {\xi - c} \right)}}{{2\sigma _{\rm{r}}^2}}} \right] $

(12)

$ {w_{\rm{p}}} = \sum\limits_{i,j \in \mathit{\Omega }} {{w_{\rm{s}}}\left( {i,j} \right){w_{\rm{r}}}\left( {i,j} \right)} $

(13)

式中:N(x, y)表示去噪后图像; w_p为归一化参数; I(i, j)表示原始图像; Ω表示滤波的窗口范围; w_s(i, j)表示高斯距离权值; w_r(i, j)表示高斯像素相似度权值; d和δ分别表示两个像素ξ和c的空间距离差和像素差; σ_s和σ_r分别表示相似度方差和空间方差。公式(10)中, 主要通过处理目标点位置像素与窗口内其余像素的距离和像素亮度信息作为权重进行加权平均之后, 得到目标位置处的像素。

自适应双边滤波方法根据图像的局部特征自适应地设置空间参数, 通过计算图像各个像素点的目标尺度得到目标点像素周围的平滑区域范围, 用来控制空间距离权重中的相似度权重^[24]。在边缘和不连续信息变化区域, 目标尺度小, 则相似度权重越小; 在平滑连续区域, 目标尺度大, 则相似度权重越大。相似度方差表达式^[25-26]为:

$ {\sigma _{\rm{s}}} = \arg \max \left\{ {\frac{{\sum {\exp \left\{ { - \frac{{{{\left[ {I\left( {x,y} \right) - I\left( {i,j} \right)} \right]}^2}}}{{2\sigma _{\rm{u}}^2}}} \right\}} }}{{{B_{ij}}\left( \mathit{\Omega } \right)}}} \right\} $

(14)

式中:I(x, y)为目标位置像素值; σ_u表示图像梯度分布的统计参数; B_ij(Ω)表示滤波窗口的区域范围。根据前文提取的多窗口, 在滤波窗口内执行双边滤波。

2 模型测试

为对比常规构造导向滤波方法、F-X滤波方法和本文方法, 设计一个三维模型, 该模型浅部为连续沉积的倾斜地层, 中深部为一个包含断层的连续沉积波状地层(图 3a), 以此为基础, 添加一定的随机噪声, 含噪模型如图 3b所示。以含噪模型为处理数据, 分别采用上述3种方法进行滤波处理。选取处理后数据的纵横测线剖面和时间切片(分别见图 3中的虚线、实线位置)进行对比, 3种方法的处理结果如图 4所示。

处理结果表明, 3种方法都能在一定程度上压制噪声, 但是不同的去噪方法其压制噪声效果不同。采用F-X滤波方法能有效抑制噪声, 但是对断层边缘信息有所损伤, 使得断层边界变得模糊, 断面成像不干脆(图 4a, 图 4b); 且去除的噪声中包含明显的断层有效信息(图 4c和图 4d)。采用常规构造导向滤波方法的结果中断层和不整合面的边缘保持以及噪声压制都优于F-X滤波方法, 但是断面边缘仍然存在一定的模糊(图 4e, 图 4f); 去除的噪声中仍然包含一定的断面和不整合面的有效反射(图 4g, 图 4h)。采用本文方法去除噪声的效果最好, 断面清晰, 边缘保持较好(图 4i, 图 4j); 去除的噪声中断面和不整合面的有效反射信息不明显(图 4k, 图 4l)。

3 实际应用

采用南海L盆地某工区的实际资料进行测试, 该地震资料除包含与自由表面相关的多次波和绕射多次波外, 在断层、裂缝等不连续地质体发育的中深层, 随机噪声较为严重^[27-28], 具体表现为:①该区块属于浅海相沉积, 水深变化剧烈, 水道砂、席状砂以及各种裂隙、裂缝和断裂较为发育, 断面反射不清晰, 地层连续性较差; ②浅层含气范围广, 中深层属于高温高压区, 波阻抗差异小, 各向异性较强, 吸收衰减严重, 模糊区较为发育, 地震资料信噪比低。由于该区随机噪声严重, 对地质体边缘的解释工作容易形成误判, 而常规去噪方法很难解决复杂构造情况下的噪声压制问题, 且易伤害地层和断面的有效反射信号, 因而影响后续地震属性体的品质, 不利于圈闭和储层预测。

分别采用F-X滤波方法、常规构造导向滤波方法和本文方法去除噪声。采用F-X滤波方法去除噪声的结果如图 5和图 6所示, 从图 5c中可以看出, 去除的噪声中出现有效反射信息(图中椭圆处), 且在去除的噪声水平切片中也出现了断面以及地层的边缘信息(图 6c), 表明F-X滤波方法损害了有效反射信息。采用常规构造导向滤波方法去除噪声的结果如图 7和图 8所示, 可以看出, 在去除的噪声剖面中也出现断层有效反射信息(图 7c), 且在去除的噪声水平切片中也出现了断面以及地层的边缘信息(图 8c), 表明常规构造导向滤波方法同样损害了有效反射信息。采用本文方法去除噪声的结果如图 9和图 10所示。对比图 5c, 图 7c和图 9c可以看出, 采用本文方法去除噪声后原始数据中断面信息基本没有受到损伤(图 9c椭圆位置); 对比图 6c, 图 8c和图 10c可以看出, 采用本文方法去除的是不规则无规律的信号, 表明滤掉的是随机噪声, 且没有损伤地层和断面的有效反射信号(图 10c)。以上结果充分验证了本文方法的有效性及合理性。

4 结论

本文提出的多窗口自适应双边滤波方法, 改进了常规双边滤波方法的单窗口技术流程, 采用多窗口滤波流程, 多窗口的选择根据地质构造信息进行, 减小其核函数对地震数据信噪比和地质体边缘信息的依赖性。与常规构造导向滤波方法和F-X滤波方法相比, 本文方法能够在有效压制随机噪声的基础上, 更好地保持边缘信息。经模型数据和实际数据验证, 该方法能够更好地实现边缘保持, 同时有效压制地震信号中的随机噪声, 去除噪声后的资料信噪比得到明显提高, 地震剖面反射同相轴连续, 保真性好, 具有较强的实用性。

由于该方法是在多窗口中进行双边滤波, 寻找最小方差的窗口作为滤波窗口进行滤波并输出结果, 与常规滤波方法相比, 不足之处在于计算时间较长。因此, 下一步工作是对该方法进行优化, 在实现保边和去噪效果的基础上提高其计算效率。