0 引言

实际地震勘探中采集的数据包含多种干扰波，有效信号受其影响而被扭曲、畸变甚至湮没。多年来，为了消除干扰噪声以突出有效信号，人们展开一系列探究并提出许多行之有效的方法。

随机噪声是地震数据中最常见的背景噪声，呈高斯性、平稳性^[1-3]分布。随机噪声压制方法的研究从未停滞，从20世纪小波滤波、F-X预测滤波、K-L(Karhunen-Loeve)变换等方法的兴起，到21世纪被广泛应用于地震数据处理领域，并被不断改进，同时又推出许多新方法，且大多取得较显著成效。

这些方法常用的有：通过估计信号特征设计的滤波方法，如F-X预测滤波^[4-6]。该方法能很好地提取有效信号特征，但由于在实际应用中需假设一些条件并存在不可预测的因素，因此该方法所能达到的滤波效果很有限。对于同相轴为直线型或近似直线型，且相邻两地震道的间距不变的地震记录，该方法具有较好的预测效果，但对于含有曲线型同相轴的记录，其去噪效果较差。

有利用相关性增强有效信号的滤波方法，如K-L变换法^[7-8]、SVD(Singular value decomposition)滤波^[9]等。K-L变换是基于统计特性的一种变换，其突出优点是去相关性好，在均方误差意义下具有最佳性能。但使用该方法需首先知道信源的协方差矩阵并求出其特征值，而这些在维数较高时不易获知和求解。然后，通过保留相关性较强的信号特征值，舍弃相关性较小的噪声特征值，再进行重构以实现信噪分离。实际应用中，该方法对水平或倾角较小的同相轴增强效果明显，而对于倾角较大的同相轴跟踪能力不足，因此增强效果很有限。SVD滤波通过对含噪信号分解后所得特征值矩阵中奇异值的合理划分，将表征有效信号的奇异值保留用以重构，得到去噪后的信号。由于奇异值矩阵中的元素值是从大到小排列的，且奇异值减小得很快，一般选取前10%甚至1%的元素即可捕捉有效信号的特征。但该方法分解出的矩阵解释性不够强，奇异值取值个数对计算精度有较大影响。因此，单纯采用SVD滤波往往难以得到理想去噪效果。

有基于时频变换的滤波方法，如基于WVD (Wigner-Ville分布)和PWVD(Pseudo Wigner-Ville分布)的时频峰值滤波^[10-11]。PWVD是WVD的加窗形式，更适用于非线性、非平稳信号的处理。在实际应用中，时频峰值滤波法多采用PWVD，以尽可能满足信号局部线性化的无偏估计条件。因此，也使该方法存在窗长选择与噪声压制之间的矛盾^[12-16]。

还有基于多尺度分析的滤波方法，如小波滤波^[17-20]、曲波滤波^[21-23]、Shearlet滤波^[24-25]等。该类方法通过对含噪数据进行多尺度分解，根据有效信号和噪声分布在不同尺度系数上的特点，采用阈值划分等手段实现信噪分离。但这些方法存在无方向性，或不能达到最优的稀疏逼近，或对阈值函数的选择较苛刻等不足，其滤波效果仍有较大提升空间。

波原子变换是由Demanet等^[26]提出的一种新型多尺度几何分析工具，可看成是二维小波包变换的变体。该变换具有良好的方向特性，每个波包的振动周期和支撑尺寸满足抛物尺度关系，即波长约等于支撑尺寸的平方^[26-29]，在此意义下，可简单地将波原子理解为方向小波与Gabor原子的插值。

波原子变换对于纹理模型具有最优的稀疏表示。一条地震记录可看作为一幅纹理图像，利用波原子变换可以更好地捕捉其特征^{[27, 29]}。在此基础上，需采用适宜的划分方法使波原子分解系数中的信号成分与噪声成分有效分离，舍弃噪声系数，然后通过波原子重构得到去噪后的地震记录。

本文采用对波原子分解所得系数进行自适应Wiener滤波^[30-36]的方案，设计出性能更优且易于实现的滤波方法压制地震勘探随机噪声以突出有效信号。通过实验验证，该方法较之于阈值分离方法可达到更优的效果。

1 方法原理 1.1 波原子变换原理

相比于小波变换、曲波变换等，波原子变换对具有振荡性的函数或具有丰富方向纹理特征的信号能呈现最优稀疏表示^{[26, 28]}。定义波原子为φ_μ(x)，μ=(j, m, n)，其中m=(m₁, m₂)，n=(n₁, n₂)。五个参量j、m₁、m₂、n₁、n₂是整数，且有j≥0，m≥0, n∈Z。定义相空间中一点(x_μ, ω_μ)为

$ \begin{aligned} &\boldsymbol{x_{\mu}}=2^{-j} \boldsymbol{n} \quad \boldsymbol{\omega_{\mu}}={\rm{ \mathit{ π} }} 2^{j} \boldsymbol{m} \\ &C_{1} 2^{j} \leqslant \max \limits_{i=1,2}\left|m_{i}\right| \leqslant C_{2} 2^{j} \end{aligned} $

(1)

式中：C₁、C₂是两个正常数，根据实际应用情况进行合理选择，可都取为1；位置向量x_μ和波向量ω_μ分别表征φ_μ(x)的空域和频域中心。

波包{φ_μ}的框架元称为波原子，波原子围绕相空间点(x_μ, ω_μ)遵循局部化条件。即对于任意M＞0，有

$ \left|\boldsymbol{\varphi_{\mu}}(x)\right| \leqslant C_{M} \cdot 2^{j}\left(1+2^{j}\left|\boldsymbol{x}-\boldsymbol{x_{\mu}}\right|\right)^{-M} $

(2)

$ C_{M} \cdot 2^{-j}\left(1+2^{-j}\left|\boldsymbol{\omega}+\boldsymbol{\omega_{\mu}}\right|\right)^{-M} $

(3)

对于任意f(x)∈L²R²，在2^-j尺度上的空间域一维波原子系数为

$ {C_{j,\mathit{\boldsymbol{m}},\mathit{\boldsymbol{n}}}} = \int {\mathit{\boldsymbol{\psi }}_\mathit{\boldsymbol{m}}^j} \left( {\mathit{\boldsymbol{x}} - {2^{ - j}}\mathit{\boldsymbol{n}}} \right)f(\mathit{\boldsymbol{x}}){\rm{d}}\mathit{\boldsymbol{x}} $

(4)

式中ψ_m^j(x-2^-jn)是实值波包ψ_m，n^j(x)在2^-j上的抽样卷积。根据Plancherel关系，得到频域一维波原子系数为

$ {C_{j,\mathit{\boldsymbol{m}},\mathit{\boldsymbol{n}}}} = \frac{1}{{2{\rm{ \mathit{ π} }} }}\int {\exp } \left( {{\rm{i}}{2^{ - j\mathit{\boldsymbol{n\omega }}}}} \right)\overline {\mathit{\boldsymbol{\hat \psi }}_\mathit{\boldsymbol{m}}^j(\mathit{\boldsymbol{\omega }})} \hat f(\mathit{\boldsymbol{\omega }}){\rm{d}}\mathit{\boldsymbol{\omega }} $

(5)

式中频域 $\overline {\mathit{\boldsymbol{\hat \psi }}_\mathit{\boldsymbol{m}}^j(\mathit{\boldsymbol{\omega }})} $为 ${\mathit{\boldsymbol{\hat \psi }}_\mathit{\boldsymbol{m}}^j(\mathit{\boldsymbol{\omega }})}$的共轭。

Hilbert变换记为H，在二维情形下，定义正交基及对偶正交基分别为

$ \mathit{\boldsymbol{\varphi }}_\mathit{\boldsymbol{\mu }}^ + \left( {{x_1},{x_2}} \right) = \mathit{\boldsymbol{\psi }}_{{m_1}}^j\left( {{x_1} - {2^{ - j}}{n_1}} \right)\mathit{\boldsymbol{\psi }}_{{m_2}}^j\left( {{x_2} - {2^{ - j}}{n_2}} \right) $

(6)

$ \mathit{\boldsymbol{\varphi }}_\mathit{\boldsymbol{\mu }}^ - \left( {{x_1},{x_2}} \right) \equiv {\rm{H}}\mathit{\boldsymbol{\psi }}_{{m_1}}^j\left( {{x_1} - {2^{ - j}}{n_1}} \right){\rm{H}}\mathit{\boldsymbol{\psi }}_{{m_2}}^j\left( {{x_2} - {2^{ - j}}{n_2}} \right) $

(7)

若记 $\mathit{\boldsymbol{\varphi }}_\mathit{\boldsymbol{\mu }}^{(1)} = \frac{{\mathit{\boldsymbol{\varphi }}_\mathit{\boldsymbol{\mu }}^ + + \mathit{\boldsymbol{\varphi }}_\mathit{\boldsymbol{\mu }}^ - }}{2}, \mathit{\boldsymbol{\varphi }}_\mathit{\boldsymbol{\mu }}^{(2)} = \frac{{\mathit{\boldsymbol{\varphi }}_\mathit{\boldsymbol{\mu }}^ + - \mathit{\boldsymbol{\varphi }}_\mathit{\boldsymbol{\mu }}^ - }}{2}$，则{φ_μ}=φ_μ⁽¹⁾, φ_μ⁽²⁾就形成了波原子紧支框架，且有

$ \sum\limits_\mathit{\boldsymbol{\mu }} {{{\left| {\left\langle {\mathit{\boldsymbol{\varphi }}_\mathit{\boldsymbol{\mu }}^{(1)},f} \right\rangle } \right|}^2}} + \sum\limits_\mathit{\boldsymbol{\mu }} {{{\left| {\left\langle {\mathit{\boldsymbol{\varphi }}_\mathit{\boldsymbol{\mu }}^{(2)},f} \right\rangle } \right|}^2}} = {\left\| f \right\|^2} $

(8)

那么，二维波原子变换系数为

$ {C_\mathit{\boldsymbol{\mu }}} = \left\langle {f,\mathit{\boldsymbol{\varphi }}_\mathit{\boldsymbol{\mu }}^{(1)}} \right\rangle + \left\langle {f,\mathit{\boldsymbol{\varphi }}_\mathit{\boldsymbol{\mu }}^{(2)}} \right\rangle $

(9)

该系数参量是一个关于j、m、n三参数的三维矩阵。j为该矩阵包含的系数矩阵组数，这里j分别取1和2，得到两组系数矩阵，每组系数矩阵都是二维矩阵，其维度分别为m₁×n₁、m₂×n₂。在后续滤波中，需对这两组系数分别进行处理；在最后的重构过程中，还需将处理过的两组系数进行组合。

对多道地震信号做二维滤波能很好地保持有效信号之间的相关性，这也是其优于一维滤波之处。本文选择二维波原子变换正是充分考虑地震记录中各道有效信号之间的相关性，使有效信号得到更大程度的凸显。

1.2 自适应Wiener滤波

维纳滤波器^[30-31]是Wiener在二十世纪四十年代提出的一种基于最小均方误差准则的自适应最佳线性滤波器。该方法被应用于众多领域，如地震勘探^[32-33]、生物医学^[34]、遥感图像^[35]以及语音信号处理^[36]等。

方法中需处理二维波原子分解系数，因此需采用二维自适应Wiener滤波。假设一个不含噪信号为ξ(l, k)，加性随机噪声为η(l, k)，那么含噪信号为

$ s(l, k)=\xi(l, k)+\eta(l, k) $

(10)

式中l、k分别为二维信号的行、列元素序号。

滤波的目的是从含噪信号中去除随机噪声，得到有效信号的估计值 $\hat \xi {\rm{(}}l, k{\rm{)}}$。而Wiener滤波所遵循的最小均方误差准则为

$ \operatorname{MSE}(\hat{\xi})=\sum\limits_{l=1}^{L} \sum\limits_{k=1}^{K}[\hat{\xi}(l, k)-\xi(l, k)]^{2} $

(11)

式中L、K分别为二维信号的行、列元素个数。二维滤波时，需选取滤波局域窗(掩模)。那么，在局域窗中采样点的局部均值和方差可分别表示为

$ \lambda=\frac{1}{P Q} \sum\limits_{p=l-r_{1} / 2}^{l+r_{1} / 2} \sum\limits_{q=k-r_{2} / 2}^{k+r_{2} / 2} s(p, q) $

(12)

$ \sigma^{2}=\frac{1}{P Q} \sum\limits_{p=l-r_{1} / 2}^{l+r_{1} / 2} \sum\limits_{q=k-r_{2} / 2}^{k+r_{2} / 2} s^{2}(p, q)-\lambda^{2} $

(13)

式中：r₁×r₂为局域窗维度；P、Q分别为局部区域内行、列元素个数，p、q为相应元素序号。对于该区域内的采样点，其二维自适应Wiener滤波可表示为

$ y(p, q)=\lambda+\frac{\sigma^{2}-\nu^{2}}{\sigma^{2}}[s(p, q)-\lambda] $

(14)

式中ν²为噪声方差。

如果噪声方差未知，则可用所有局部估计方差的均值代替。根据式(14)求出各局部区域的信号估计值，再将这些估计值合并为一个整体矩阵即为最终估计值。Wiener滤波的此种实现方法有别于其他传统实现方法。传统实现方法中，需求解Wiener滤波器的响应函数，然后用输入信号卷积(时域实现)或乘积(频域实现)此响应函数而得到输出信号。以最小均方误差准则(式(11))为滤波误差约束条件。

但是，该误差公式用到信号的期望值，而该期望值在实际中很难得到。式(14)所示的Wiener滤波实现过程不依赖于期望信号，主要是利用局部均值与方差得到有效信号的估计值，是一种更实用的自适应滤波方法。

2 在地震数据处理中的应用 2.1 处理过程概述

前文述及，对波原子分解系数进行滤波处理，需去除不必要的系数，即表征噪声的系数，保留表征有效信号的系数，然后对这些有效信号系数进行波原子重构即可得到滤波后的信号。本文采用自适应Wiener滤波对分解系数进行处理，是利用其在最小均方误差意义下可达到最佳线性滤波的特点。

对地震记录进行处理时，一般采样点数与道数不是2的整数次幂，因此首先要将地震数据通过补零得到行、列元素个数均为2的整数次幂的矩阵。在分解过程中，得到如式(9)中的φ_μ⁽¹⁾和φ_μ⁽²⁾，依照式(12)~式(14)分别对这两组系数进行二维自适应Wiener滤波。滤波时，需根据记录中的同相轴形态与分布特征，合理选择局域窗，即对于同相轴不对称分布的记录一般选用矩形窗，而对于同相轴呈近似对称分布的记录一般选用正方形局域窗。对于双曲同相轴记录与线性同相轴记录选取的窗函数的维度也有所区别，同等噪声条件下，前者滤波窗的r₂一般小于后者。最后将滤波后的波原子系数进行重构，即可得到能量增强后的地震记录。

2.2 模拟地震记录处理

首先针对含有线性同相轴的模拟地震记录进行实验。该记录为2ms采样，道间距为10m，共有50道。其中包含四个线性同相轴，每个同相轴由主频为30Hz的Ricker子波构成。对该不含噪模拟地震记录(图 1a)加入高斯白噪声，使其信噪比约为-5dB(图 1b)。

SVD阈值(方法1)、二维多尺度小波分解阈值(方法2)和二维波原子分解阈值(方法3)三种滤波方法中，阈值的选取须综合考虑随机噪声压制与有效信号保幅，并历经多次实验测试。三种方法滤波后重构的地震记录依次如图 1c~图 1e所示。

对含噪记录进行二维波原子分解系数自适应Wiener滤波(本文方法)时，考虑记录中有效同相轴分布情况，选取的局域窗为长宽不等的矩形窗，即r₁=5、r₂=187。当噪声强度改变时，可适当调整r₁和r₂值。由此，本文方法滤波后重构的地震记录如图 1f所示。

从滤波结果可见：方法1对背景噪声的去除及对有效信号的保留效果均较差(图 1c中矩形框)；方法2对有效同相轴的损伤较严重且在同相轴边缘产生阴影和模糊现象(图 1d中椭圆和矩形框)；方法3滤波效果较前二者有所提高，但模拟记录两边地震道有效信号损失较严重(图 1e中矩形框)；本文方法处理效果最好，在随机噪声压制和有效信号保持方面均表现较优(图 1f)。据图中右侧能量值色标也能看出图 1f中有效信号能量显著强于图 1c~图 1e。经定量计算，得到四种方法滤波后记录的信噪比分别为0.1686、7.6282、7.8338和8.7075dB。

对上述纯净记录加入不同强度随机噪声，得到不同信噪比含噪记录，分别采用四种方法做滤波实验，滤波前后地震记录信噪比数据如表 1所示。

从表 1可见：方法1所能达到的信噪比很低，对原含噪记录信噪比的提高较少；方法2、方法3的滤波后记录能达到较高信噪比，其中方法3信噪比提升幅度稍高于方法2，但二者都存在同相轴局部损失较多的问题，且方法2还会导致有效同相轴边缘不清晰；本文方法相对于其他三种方法优势更明显，对有效同相轴的增强及边缘保留效果最好，对含噪记录信噪比的提升能力最强。

本次实验采用一台高性能便携式工作站。在相同计算机条件下，上述四种滤波方法处理图 1b所示含噪记录的时间依次为：0.0071、0.7044、3.4747、3.9913s。可见方法1虽用时最短，但效果不佳；方法2用时次短，但效果也欠理想；方法3和本文方法的用时虽较前二者长，但本文方法处理效果更好。

再对含有双曲同相轴的地震记录(图 2)做滤波实验。该数据采样间隔为1ms，道间距为20m，采用主频为35Hz的Ricker子波，共有76道。层速度分别为2300、2500、2900m/s。对不含噪记录(图 2a)加入高斯白噪声使其信噪比为-5dB(图 2b)，三种方处理后的地震记录分别如图 2c~图 2e所示。

由于图 2a同相轴的形态和分布与图 1a不同，因此需调整Wiener滤波窗函数。调整时须兼顾随机噪声压制与有效信号保持。此时，选取r₁=5、r₂=117，即可得到本文方法处理后的地震记录(图 2f)。

观察四种滤波方法所得滤波记录(图 2)可知：方法1对背景噪声的消除能力很有限且对有效同相轴损失较多(图 2c中矩形和椭圆框)；方法2和方法3能较有效地压制随机噪声，但对有效信号的损失不容忽视，前者会造成同相轴边缘模糊(图 2d中矩形和椭圆框)，后者会使双曲同相轴拱起部分损失较严重(图 2e中椭圆框)；而本文方法能在随机噪声压制和有效信号保持方面做到很好权衡(图 2f)。四种方法处理后记录的信噪比分别为0.3301、7.1616、7.4222、9.1073dB。

同样，对上述纯净记录加入不同强度随机噪声，得到不同信噪比含噪记录，分别采用四种方法进行滤波实验，滤波前、后的信噪比数据如表 2所示。分析对比该表数据，也能得到与表 1相似的结论。

2.3 实际地震数据处理

选取实际共炮点记录做滤波对比测试。该共炮点记录共有168道，采样间隔为1ms，仅截取记录的上部(0~2048ms)。首先分别采用方法2、方法3和本文方法进行滤波处理(图 3)。此次处理选取的Wiener滤波窗为正方形，即r₁=r₂=95。

从滤波结果(图 3)可见：方法2(图 3b)和方法3(图 3c)滤波效果明显不理想，对有效同相轴造成较大程度的损失及畸变；而本文方法使有效同相轴更清晰、连续，能量得到显著增强(图 3d中梯形和矩形框)。

对该实际记录还应用方法1和一维时频峰值滤波法进行处理。方法1处理时，若选取过少特征值，就会导致有效同相轴损失太多；若选取较多特征值，就几乎没有去噪效果。对于此次实际记录，选取分解后奇异值矩阵中的前60个特征值进行重构较适宜(图 4b)。但该方法对背景噪声的压制效果非常有限，且对有效同相轴造成了一定程度的畸变。

进行时频峰值滤波时，滤波窗长不宜太大或太小，若太大虽能较好地压制随机噪声，但有效信号损失较大；若太小则不能很好地压制随机噪声。对于该实际记录，选取滤波窗长为13点(图 4c)。与本文方法所得结果(图 4d)对比，可见一维时频峰值滤波处理结果虽然背景噪声有所压制，但同相轴的清晰度和连续性提升效果欠理想，而本文方法处理后的记录中同相轴的清晰度更高、连续性更好。

3 结论

本文提出二维波原子分解系数自适应Wiener滤波方法，对地震资料进行随机噪声消减处理以增强有效同相轴。其实现过程是，先对地震记录进行二维波原子分解，再对所得系数进行二维自适应Wiener滤波，最后对滤波后的系数进行波原子重构得到同相轴增强的记录。通过对人工合成记录与实际记录的滤波处理，得到以下认识和结论：

(1) 通过对模拟地震记录的滤波实验，并与SVD阈值滤波、二维多尺度小波分解阈值滤波及二维波原子分解阈值滤波结果作比较，验证本文所提滤波方法在地震随机噪声压制及有效信号保持方面具有优越性。

(2) 采用二维多尺度小波分解阈值滤波、二维波原子分解阈值滤波、SVD阈值滤波、一维时频峰值滤波及二维波原子分解系数自适应Wiener滤波五种方法对实际共炮点记录进行滤波处理，并将处理后的结果作对比，发现本文本文方法能够较彻底地压制记录中的随机噪声，同时使有效同相轴更清晰、连续性更好。

因此，应用本文方法对地震数据进行处理，可为后续的地震成像及综合解释提供可靠资料，对提高地震资料解释的准确性具有重要意义。