复杂地表和地下地质条件下, 地震资料的信噪比问题是地震勘探面临的核心问题之一。例如, 我国南方山地工区由于地表和地下地质条件双复杂, 地震波在地表和地下传播时复杂多变, 噪声干扰使单炮记录中有效地震波场难以呈现清晰的反射同相轴, 降低了地震资料的信噪比, 增加了地震资料的处理难度。地震噪声主要包括相干噪声和随机噪声。相干噪声主要是指具有一定运动学规律的干扰信号, 如多次波、面波、直达波、鬼波等。随机噪声是指没有规律可循的环境噪声和次生散射噪声等。从时间域看, 在整个地震记录中都有可能存在随机噪声, 而且其传播速度不同, 方向不定; 从频率域看, 既存在低频随机噪声也存在高频随机噪声。一般来说, 如果地震数据中的信号成分或者噪声成分可以用某种变换的一个或少数的基函数来表示, 则在分离信号和噪声时这种变换将发挥很大作用^[1]。这是因为这些成分通常在变换域中稀疏性强, 因而能较好地分离信号和噪声, 如奇异值分解方法^[2]。

地震有效信号可以根据其横向连续性进行识别, 而这种连续性被用来区分背景随机噪声和地震有效信号^[3]。在地震资料处理中用于随机噪声压制的方法技术有很多, 其中以f-x域预测误差滤波方法最为常用^[4-6], 该方法又称为三维随机噪声衰减技术(3DRNA), 最初由CANALES^[7]提出, 后来由GULUNAY^[8]根据TREITEL^[9]的复数序列预测发展起来。其基本假设条件是如果时间域局部区域内地震信号的同相轴是线性的, 则频率空间域同相轴信号是可以预测出来的^[10-11], 而随机噪声不是相干信号, 所以不能有效预测。根据此假设条件, 就可以设计不同频率的预测滤波算子, 然后用不同频率的预测滤波算子预测出地震数据中的相干信号, 进而去除不可预测的随机噪声。f-x域预测滤波方法已经普遍应用于地震资料的噪声压制, 是目前有效压制随机噪声、提高地震资料信噪比的方法之一, 但是该方法在实际应用中仍存在模糊有效不连续信息(如断层)的问题, 从而对后续的地震属性分析, 特别是地质不连续性分析带来影响^[12-14]。通常只能采取人工解释断层位置后, 再进行非断层区噪声压制以保证断层信息不受损失, 工作效率低, 效果受人工解释的精度影响较大。针对常规预测误差滤波器存在的模糊有效不连续信息的不足和人工干预影响大的问题, 本文方法利用不同方向预测滤波器对地震信号中的不连续性影响的差异, 提出一种自适应智能化双向滤波策略, 在有效压制随机噪声的同时实现断层保护。

1 方法原理

定义有效信号、含噪声地震信号、噪声分别为s(i, j)、x(i, j)和n(i, j), 其中i是地震道数, j是时间。含噪声地震信号可以表示为:

$ \mathit{\boldsymbol{x}}(i,j) = \mathit{\boldsymbol{s}}(i,j) + \mathit{\boldsymbol{n}}(i,j) $

(1)

沿时间方向进行傅里叶变换(FFT), 方程(1)可以改写成:

$ \mathit{\boldsymbol{X}}(i,\omega ) = \mathit{\boldsymbol{S}}(i,\omega ) + \mathit{\boldsymbol{N}}(i,\omega ) $

(2)

式中:X(i, ω)、S(i, ω)和N(i, ω)分别为f-x域中的含噪声地震信号、有效信号和噪声; ω为频率。对于一个单独的频率成分k, 沿地震道i提取一个一维信号如下:

$ Y(i) = X(i,k) $

(3)

由于线性地震同相轴在f-x域是可预测的, 因此, 为了保证地震同相轴是局部线性的, 需要先将地震数据分成若干小块, 然后将每个小块地震数据变换到f-x域再进行预测滤波。设每个小块包含M个地震道, 每道有I个采样点。对于频率k, 可以得到一个L阶的预测滤波器f:

$ \mathit{\boldsymbol{Af}} = \mathit{\boldsymbol{d}} $

(4)

式中:A为输入数据矩阵; d为预测的有效信号。表达式分别为:

$ \mathit{\boldsymbol{f}} = [{f_1},{f_2},{f_3}, \cdots ,{f_L}] $

(5)

$ \mathit{\boldsymbol{d}} = [Y(L + 1),Y(L + 2), \cdots ,Y(M)] $

(6)

$ \mathit{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} {Y(1)}&{Y(2)}& \cdots &{Y(L)}\\ {Y(2)}&{Y(3)}& \cdots &{Y(L + 1)}\\ \vdots & \vdots &{}& \vdots \\ {Y(M - L)}&{Y(M - L + 1)}& \cdots &{Y(M - 1)} \end{array}} \right] $

(7)

方程(4)的最小二乘解为:

$ \mathit{\boldsymbol{f}} = {({\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}})^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{d}} $

(8)

得到预测滤波器f后, 可以得出预测的有效信号$ \boldsymbol{\hat d}$:

$ \mathit{\boldsymbol{\hat d}} = \mathit{\boldsymbol{Af}} $

(9)

传统的f-x域预测滤波, 需要进行前向和反向滤波, 其滤波算子分别记为f₁和f₂, 滤波结果分别记为$\boldsymbol{\hat d}_1$和$\boldsymbol{\hat d}_2$, 变换到t-x域后, 传统的f-x域滤波结果$\boldsymbol{\hat x}$表示为:

$ \mathit{\boldsymbol{\hat x}}(i,j) = 0.5{{\mathit{\boldsymbol{\hat d}}}_1}(i,j) + 0.5{{\mathit{\boldsymbol{\hat d}}}_2}(i,j) $

(10)

可以看出, 传统f-x域预测滤波方法采用前向和反向滤波器的预测结果相同, 所以对于连续无间断的采样点可以取得较好的滤波结果。但对那些存在间断的采样点, 一般来说, 被预测的采样点和用来进行预测的采样点分布在间断的同一侧时, 可以得到较好的滤波结果; 反之, 被预测的采样点和用来进行预测的采样点分布在间断的两侧时, 会存在较大的误差, 因此产生模糊边界效应。在本文方法中, 考虑到不同方向的预测滤波器受地震信号中间断区的影响不同以及预测精度的差异, 提出一种自适应智能化双向滤波的策略:

$ \mathit{\boldsymbol{\hat x}}(i,j) = a(i,j){{\mathit{\boldsymbol{\hat d}}}_1}(i,j) + [1 - a(i,j)]{{\mathit{\boldsymbol{\hat d}}}_2}(i,j) $

(11)

式中:a(i, j)为加权值, 其取值范围为[0, 1]。

假设地震断层是几乎垂直于地面的, 一般来说, 用前向滤波器来预测断层左侧的采样点时, 其预测误差会远小于利用反向滤波器预测该采样点, 而对断层右侧的采样点进行预测时则相反。当采样点周围不存在断层时, 则前向和反向滤波器的预测误差相差不大。基于上述特征, 本文提出可以自动识别地震数据中断层位置信息, 进而获得加权函数的方法。也可以针对不同的问题(如强反射边界的间断区等)设定不同的加权函数。本文选择的加权函数a(i, j)为:

$ \begin{array}{*{20}{c}} {a(i,j) = {\rm{tanh}} [\frac{{{e_2}(i,j)}}{{{e_1}(i,j) + {e_2}(i,j)}} \times 4 - 2] \times }\\ {0.5 + 0.5} \end{array} $

(12)

式中: $e_{1}(i, j)=\left| \boldsymbol{\hat{d}}_{1}(i, j)- \boldsymbol{x}(i, j)\right|, e_{2}(i, j)=\left|\hat{\boldsymbol{d}}_{2}(i, j)-\boldsymbol{x}(i, j)\right|$, 分别对应前向和反向滤波器的预测误差。a(i, j)越接近于0.5, 表明双向的预测误差越小, 说明(i, j)位置不属于断层区域, 如果a(i, j)远离0.5, 表明双向的预测误差越大, 说明(i, j)位置属于断层区域。

在某些地质构造特别复杂或地震资料信噪比很低的探区, 分块地震数据不能满足f-x域预测滤波方法的线性同相轴假设条件, 即使采用上述方法压制噪声仍会导致有效信号的损失。针对此问题, 在上述方法的基础上做进一步优化, 提出了噪声回归的处理策略, 将损失的有效信号补偿回来。首先得到所消除的噪声为:

$ \mathit{\boldsymbol{e}}(i,j) = \mathit{\boldsymbol{x}}(i,j) - \mathit{\boldsymbol{\hat x}}(i,j) $

(13)

估计得到|e(i, j)|的均值, 记为m_e= $\frac{\sum\limits_{i=1, j=1}^{M , I}| \boldsymbol{e}(i, j)|}{M I}$, 构造如下加权系数:

$ w(i,j) = \left\{ {\begin{array}{*{20}{c}} 1&{|\mathit{\boldsymbol{e}}(i,j)| < {\sigma _1}{m_e}}\\ {1 - \frac{{|\mathit{\boldsymbol{e}}(i,j)| - {\sigma _1}{m_e}}}{{{\sigma _2}{m_e} - {\sigma _1}{m_e}}} + \frac{1}{{{\sigma _2}}}}&{{\sigma _1}{m_e} \le |\mathit{\boldsymbol{e}}(i,j)| \le {\sigma _2}{m_e}}\\ {1/{\sigma _2}}&{|\mathit{\boldsymbol{e}}(i,j)| > {\sigma _2}{m_e}} \end{array}} \right. $

(14)

其中, σ₁, σ₂是用户给定的噪声水平参数, 且σ₁ < σ₂。最终输出的滤波结果为:

$ \mathit{\boldsymbol{\hat s}}(i,j) = \mathit{\boldsymbol{x}}(i,j) - \mathit{\boldsymbol{e}}(i,j)w(i,j) $

(15)

从公式(14)中可见, 当|e(i, j)| < σ₁m_e时, 表示噪声能量偏低, 此时加权值为1, 即视为无间断区; 当|e(i, j)|>σ₂m_e时, 表示噪声能量偏高, 加权值为反比例函数, 即视为间断区需要保护; 当σ₁m_e≤|e(i, j)|≤σ₂m_e时, 表示无间断和间断的过渡区间。用户可以根据实际情况通过调整σ₁, σ₂参数来调节噪声压制的水平, 从而满足实际生产的需要。

根据上述原理设计如下去噪流程。

1) 根据公式(9)在原始数据上逐点进行双向预测算子计算, 并形成双向的预测结果$\boldsymbol{\hat d}_1$和$\boldsymbol{\hat d}_2$;

2) 通过公式(12)获得加权函数a(i, j), 可根据加权函数的大小自动识别间断区在剖面中的位置, 并通过公式(11)获得双向滤波后的去噪结果$\boldsymbol{\hat x}$;

3) 通过公式(13)获得噪声剖面e;

4) 根据步骤2)识别的间断区位置, 再给定噪声水平参数σ₁, σ₂, 构造加权系数w(i, j), 根据公式(15)将噪声剖面中的有效信号补充回信号剖面, 即可实现断层保护的随机噪声压制。

在整个去噪流程中, 只需人为给定两个噪声水平参数σ₁、σ₂, 按照去噪流程设计程序, 即可完成随机噪声的压制和断层附近有效信号的保护。由于有效信号具有线性特征, 因此其在噪声剖面上表现为连续的同相轴和强振幅, 只要识别强振幅的范围就可以确定间断区间的范围。在实际应用中, 可以自动统计整个噪声剖面e的能量范围, 然后设置σ₁和σ₂两个参数为整个噪声剖面相对能量的百分比即可, 例如σ₁=80%, σ₂=90%, 其含义为90%最大能量以上的视为间断区必须完全保护, 80%~90%最大能量的为过渡区进行部分保护, 80%最大能量以下的可以不保护。

2 模拟数据测试

利用如图 1所示的模拟数据验证本文方法压制随机噪声的效果。图 1a为无噪声原始地震剖面, 剖面中间有一个断层, 表现为间断的地震同相轴。图 1b为在图 1a所示的地震剖面中加入10 dB的随机噪声的地震剖面, 可见, 剖面中间的断层被模糊化, 水平地震信号也受到了随机噪声的污染。分别利用常规f-x域随机噪声压制方法和本文方法对图 1b所示地震数据进行随机噪声压制(图 2, 图 3), 并对比经噪声压制后的信号剖面和被压制的噪声剖面, 从而分析两种噪声压制方法的效果。

利用常规f-x域随机噪声压制方法对含噪声数据(图 1b)进行处理后的结果如图 2a所示。从图 2a中可见, 中间的断层位置被模糊化了, 这会导致后续的地震资料解释不能准确地确定断层的位置。从噪声剖面(图 2b)上可清晰地看到反映断层位置不连续性的有效信号被当作噪声压制了。由此可见, 常规f-x域随机噪声压制方法的局限性对断层的准确成像影响较大。

利用本文方法对含噪声数据(图 1b)进行随机噪声压制后的处理结果如图 3a所示。从图 3a中可以看到, 断层位置清晰, 表明在利用本文方法压制随机噪声时没有损伤有效信号, 使得断层得到了准确成像。从图 3b所示的噪声剖面中也可以看出, 相对于常规的f-x域随机噪声压制方法, 本文方法能更好地保护断层信息, 有利于后续地震属性分析和断层检测。

3 实际数据处理

选择某实际工区的地震资料的成像结果进行噪声压制的测试。原始数据中包含有效反射信号和噪声, 其中包含强反射界面同样对本文方法进行验证。

3.1 实例1

首先利用图 4所示的地震剖面进行测试, 该数据为偏移成像后的叠加数据。分别采用常规f-x域随机噪声压制方法和本文方法对该数据进行处理, 结果如图 5和图 6所示。由图 4可见, 该数据中存在非常强的基底反射界面和断层, 强反射界面在常规的f-x域去噪中能量会有损失, 在图 5b所示的噪声剖面中残留大量的基底反射能量。由于基底的振幅大, 类似于振幅间断点, 采用本文方法同样能够识别强反射边界的位置, 并且在噪声压制过程中对基底有效信号进行保护, 如图 6b所示的噪声剖面中无基底有效信号的损伤。对比图 5a和图 6a可以看出, 采用本文方法的噪声压制结果较好, 表明采用本文方法能够在保护有效信号的同时进行随机噪声压制。

3.2 实例2

为了进一步验证本文方法的有效性, 选取断层位置存在明显的断面波反射的剖面(图 7)进行方法的测试分析。从图 7所示的剖面中可以清晰看到断面波的位置和形态, 剖面中存在随机噪声干扰, 并且不易识别, 很容易在噪声压制的同时损失有效信号的成像振幅。分别采用常规f-x域随机噪声压制方法和本文方法对该剖面进行处理, 结果如图 8和图 9所示。对比图 8a和图 9a可以看出, 常规f-x域随机噪声压制方法与本文方法整体去噪效果相当, 本文方法的结果在断层位置处效果更好。从图 8b和图 9b可以看出, 常规f-x域随机噪声压制方法对于断层信息的损失是显而易见的, 而本文方法对于断层的信息保护很好, 无明显的有效信号损失。可见, 与常规f-x域随机噪声压制方法相比, 本文方法对间断区域的识别及保护断层有效信号具有明显的优势。

4 结论

本文从随机噪声和有效信号的基本特征差异出发, 引入了预测误差加权函数, 通过间断点双向预测误差大的特点智能化地确定了间断区间, 然后引入加权系数对间断区间进行有效信号保护, 进而避免了常规f-x域预测滤波的线性假设造成的间断信号的损失问题, 因此能够实现在保护断层有效信号的前提下, 进行随机噪声压制。本文方法对于间断区间的识别完全是自动完成, 无需人工干预, 可以实现智能化的随机噪声压制。同时, 还提供了人工干预的接口参数可以调整模块的噪声压制效果, 应用方便, 有利于工业应用时的效率提升。