0 引言

由于表层激发条件的限制以及其他各种不可抗干扰因素的影响，仅在野外采集阶段完全压制随机噪声不太现实，需要在室内用有效的数学去噪方法提高叠前地震资料的信噪比，以满足不同处理环节的需求^[1-4]。目前，有许多行之有效的随机噪声压制方法，其中稀疏变换应用较为广泛，包括傅里叶变换^[5-8]、小波变换^[9-10]、Radon变换^[11-12]、Contourlet变换^[13]、曲波变换^[14]以及Seislet变换^[15]等。该类方法根据随机噪声和有效波在稀疏变换域中能量分配的差异，对稀疏系数采用适当的阈值去除随机噪声的影响。为了取得较好的去噪效果，稀疏基需要尽可能地捕获地震波的有效信息，并且少数较大的稀疏系数能够代表信号的主要特征，而大部分较小、被滤除的系数不影响原始数据的主要特征^[16]。曲波变换能够有效地表示地震数据的尺度性和方向性，更加稀疏地表示地震波场局部细节特征，所以许多学者选用曲波变换方法去噪，且效果显著^[17-20]。

现有曲波变换去噪方法的前提是地震数据为均匀网格采样，而对于非均匀网格采样的含噪数据则效果不佳。由于野外地形条件及施工环境的制约，在很多情况下，地震数据常为非均匀采样^[21]。如果在去噪过程中将非均匀采样视为均匀采样进行处理，就得不到连续的地震波场，自然也不能充分压制噪声干扰，从而会影响到后续其他处理方法的应用。Hennenfent等^[22]提出基于非均匀采样的曲波变换重建方法，但未讨论该方法在噪声压制中的应用。

本文在前人研究的基础上，针对非均匀采样地震数据，在传统曲波变换的基础上，引入非均匀傅里叶变换，建立均匀曲波系数与空间非均匀采样地震道之间的规则化反演算子；然后选择合适的噪声估计值，使用线性Bregman方法反演，在迭代过程中采用软阈值对曲波系数去噪，由反演得到无噪声的均匀曲波系数；再进行常规曲波反变换，得到去噪后的地震数据，由此形成基于非均匀曲波变换和线性Bregman方法的随机噪声压制方法。模型数据和实际数据的实验结果表明了本文方法的有效性。

1 非均匀曲波变换

Candès等^[23]提出了第二代曲波变换，从而使曲波变换更容易被理解、运算效率更高、实现更简单，能够为地震信号提供最优的稀疏表示方式。实际上，实现快速离散曲波变换主要包括两个步骤：①对地震数据应用二维傅里叶变换，得到频率—波数域系数；②在频率—波数域形成角度楔形，将每一个楔形围绕到原点重新装配，并对每一个装配好的楔形应用二维傅里叶反变换，得到离散曲波系数。参考文献[22]定义曲波正变换算子为

$ \mathit{\boldsymbol{A}} = \mathit{\boldsymbol{TF}} $

(1)

式中：F实现了上述离散曲波变换第①步：T则实现了第②步。

定义曲波反变换算子为

$ {\mathit{\boldsymbol{A}}^{\rm{H}}} = {\mathit{\boldsymbol{F}}^{\rm{H}}}{\mathit{\boldsymbol{T}}^{\rm{H}}} $

(2)

式中：上标“H”表示共轭转置；F^H表示二维傅里叶反变换，将频率—波数域转换到时间—空间域中；T^H表示曲波平铺算子，即将曲波系数变换到频率—波数域的过程。

式(1)中，由于傅里叶变换要求数据均匀采样，所以常规的二维曲波变换不能处理非均匀采样数据，否则会严重扭曲原始地震数据的真实位置，使地震同相轴不连续，得不到正确的频率—波数域系数。用F_x^H代表沿着空间轴的一维傅里叶反变换算子，F_t^H表示沿着时间轴的一维傅里叶反变换算子，则二维傅里叶反变换可以表示为两个维度变换的克罗内克乘积，即

$ {\mathit{\boldsymbol{F}}^{\rm{H}}} = \mathit{\boldsymbol{F}}_x^{\rm{H}} \otimes \mathit{\boldsymbol{F}}_t^{\rm{H}} $

(3)

式中“⊗”表示克罗内克乘积。在这个分解变换中，时间方向是均匀理想采样，不需要重建。因此，可以用非均匀快速傅里叶反变换N_x代替一维傅里叶反变换算子F_x^H，则非均匀曲波反变换算子可以定义为

$ \mathit{\boldsymbol{C}} = ({\mathit{\boldsymbol{N}}_x} \otimes \mathit{\boldsymbol{F}}_t^{\rm{H}}){\mathit{\boldsymbol{T}}^{\rm{H}}} $

(4)

该算子描述了非均匀采样下离散曲波系数与地震数据之间的关系。

由于快速离散曲波变换具有紧支撑性，正、反曲波变换前、后能量无损失，即

$ \mathit{\boldsymbol{A}}\;{\mathit{\boldsymbol{A}}^{\rm{H}}} = \mathit{\boldsymbol{I}} $

(5)

式中I为单位矩阵。对于新的非均匀曲波反变换算子C，由于包含非均匀快速傅里叶变换，那么算子C就不满足式(5)。为此定义非均匀曲波正变换算子C⁺为

$ \begin{array}{l} {\mathit{\boldsymbol{C}}^ + }:\mathit{\boldsymbol{d}} \to \mathit{\boldsymbol{\tilde x}} = \arg \min \left\| {{x}} \right\|_1\\ \;\;\;\;\;\;\;满足\;\;\;\;\;\;\;{\left\| {\mathit{\boldsymbol{d}} - \mathit{\boldsymbol{Cx}}} \right\|_2} \le \sigma \end{array} $

(6)

式中：d表示非均匀含噪地震数据；σ为大于零的噪声均方根估计值，主要根据噪声能量计算该值的大小，该值的选择尤为关键，直接决定去噪效果。如果σ=0，则表示仅对非均匀含噪数据进行规则化，不进行噪声压制。从式(6)可知，C⁺将非均匀含噪数据d非线性映射到曲波系数向量${\mathit{\boldsymbol{\tilde x}}} $，有无限多向量x满足d=Cx，${\mathit{\boldsymbol{\tilde x}}} $是满足最小L₁范数中的一个，并且该式两边相等保证了从空间非均匀采样地震数据获得的曲波系数通过反变换算子C，可以再次精确地产生相同位置的非均匀采样地震数据。Hennenfent等^[22]采用谱梯度投影算法求解式(6)，但该算法参数较多，运行速度较慢，本文选用线性Bregman方法求解该L₁范数最优化问题^[24]。

2 线性Bregman方法

式(6)最优化问题可转为求解下述基追踪规则化问题^[25]

$ \begin{array}{l} \mathop {\min }\limits_x \left( {\lambda {{\left\| \mathit{\boldsymbol{x}} \right\|}_1} + \frac{1}{2}\left\| \mathit{\boldsymbol{x}} \right\|_2^2} \right)\\ \;\;\;\;\;\;满足\frac{1}{2}\left\| {\mathit{\boldsymbol{d}} - \mathit{\boldsymbol{Cx}}} \right\|_2^2 \le {\sigma ^2} \end{array} $

(7)

式中λ是一个阈值权衡因子，平衡L₁范数和L₂范数的比重，可以在第一次迭代过程计算。线性Bregman方法求解迭代式为

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{z}}_{k + 1}} = {\mathit{\boldsymbol{z}}_k} - {l_k}\mathit{\boldsymbol{C}}_k^{\rm{H}}({\mathit{\boldsymbol{C}}_k}{\mathit{\boldsymbol{x}}_k} - {\mathit{\boldsymbol{d}}_k})\\ {\mathit{\boldsymbol{x}}_{k + 1}} = {S_\lambda }({\mathit{\boldsymbol{z}}_{k + 1}}) \end{array} \right. $

(8)

式中：z_k表示第k次迭代得到的曲波系数向量；l_k为动态步长，定义为

$ {l_k} = \frac{{\left\| {{\mathit{\boldsymbol{C}}_k}{\mathit{\boldsymbol{x}}_k} - \mathit{\boldsymbol{d}}} \right\|_2^2}}{{\left\| {\mathit{\boldsymbol{C}}_k^{\rm{H}}({\mathit{\boldsymbol{C}}_k}{\mathit{\boldsymbol{x}}_k} - \mathit{\boldsymbol{d}})} \right\|_2^2}} $

(9)

软阈值函数为

$ {S_\lambda }\left( \mathit{\boldsymbol{x}} \right) = {\rm{sign}}\left( \mathit{\boldsymbol{x}} \right)\max \left( {\left| \mathit{\boldsymbol{x}} \right| - \lambda , 0} \right) $

(10)

式(8)只能进行规则化重建，而不能够同时压制噪声。在式(8)中引入新的投影函数$ \mathop \prod \limits_\sigma ({\mathit{\boldsymbol{C}}_k}{\mathit{\boldsymbol{x}}_k} - {\mathit{\boldsymbol{d}}_k})$处理地震数据中的随机噪声^[26]，形成新的迭代公式

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{z}}_{k + 1}} = {\mathit{\boldsymbol{z}}_k} - {l_k}\mathit{\boldsymbol{C}}_k^{\rm{H}}\mathop \prod \limits_\sigma ({\mathit{\boldsymbol{C}}_k}{\mathit{\boldsymbol{x}}_k} - {\mathit{\boldsymbol{d}}_k})\\ {\mathit{\boldsymbol{x}}_{k + 1}} = {S_\lambda }({\mathit{\boldsymbol{z}}_{k + 1}}) \end{array} \right. $

(11)

投影函数$ \mathop \prod \limits_\sigma ({\mathit{\boldsymbol{C}}_k}{\mathit{\boldsymbol{x}}_k} - {\mathit{\boldsymbol{d}}_k})$定义为

$ \begin{array}{l} \mathop \prod \limits_\sigma ({\mathit{\boldsymbol{C}}_k}{\mathit{\boldsymbol{x}}_k} - {\mathit{\boldsymbol{d}}_k})\\ = \max \left( {0, 1 - \frac{\sigma }{{{{\left\| {{\mathit{\boldsymbol{C}}_k}{\mathit{\boldsymbol{x}}_k} - {\mathit{\boldsymbol{d}}_k}} \right\|}_2}}}} \right)({\mathit{\boldsymbol{C}}_k}{\mathit{\boldsymbol{x}}_k} - {\mathit{\boldsymbol{d}}_k}) \end{array} $

(12)

该式可以保证在每次迭代过程中对含噪地震数据进行噪声压制，从中也可以看出线性Bregman方法实现非常简单，没有过多的调节参数。

3 噪声压制原理

一个含噪声二维地震信号模型可以表示成为

$ d(i, j) = s(i, j) + e(i, j) $

(13)

式中：s为真实信号；e为噪声。去噪过程就是从非均匀含噪地震信号d中提取真实信号s，去除噪声干扰信号e(图 1)。

非均匀曲波变换依赖于式(6)中的曲波系数。对于随机干扰噪声的去除，在采用线性Bregman方法求解式(6)的过程中，由于较大曲波系数代表有效信号本身，通过对曲波系数进行软阈值处理^[27-28]，经过多次迭代就可以得到无噪声的均匀曲波系数估计值${\mathit{\boldsymbol{\tilde x}}} $。因此，去噪后的地震波场${\mathit{\boldsymbol{\tilde s}}} $可以通过下式得到

$ \mathit{\boldsymbol{\tilde s}} = {\mathit{\boldsymbol{A}}^{\rm{H}}}\mathit{\boldsymbol{\tilde x}} $

(14)

式中A^H是常规的曲波反变换算子。通过对反演后的无噪声均匀曲波系数应用常规曲波反变换，可得到均匀采样后的去噪地震数据。

4 数值模拟

信噪比定义为^[29]

$ R = - 10\lg \frac{{{{\left\| {\mathit{\boldsymbol{s}} - {\mathit{\boldsymbol{s}}_0}} \right\|}_2}}}{{{{\left\| {{\mathit{\boldsymbol{s}}_0}} \right\|}_2}}} $

(15)

式中：s₀表示原始均匀不含噪声数据；s表示去噪后的均匀地震数据。R越高，表示去噪效果越好。

图 2a为合成的256道地震数据，有4个反射波同相轴，各反射同相轴能量有差异，采样间隔为1ms，道间距为5m，每道1024个采样点。加入标准差为0.05的高斯白噪声，如图 2b所示，其信噪比为-5.11dB。对含噪数据进行空间均匀傅里叶变换，然后再进行空间非均匀傅里叶反变换，得到新的空间非均匀采样下的256道地震数据，如图 2c所示，信噪比为-5.25dB，道间距范围为0~10m。如果直接用均匀曲波变换法对图 2c数据去噪，则扭曲的同相轴得不到校正，导致去噪后的地震记录误差较大。因此首先采用非均匀曲波变换方法对图 2c进行规则化处理，此时噪声估计值σ=0，曲波变换所选择的尺度数为5，在第二个最粗尺度上的角度数为32。规则化结果如图 2d所示，可以看见规则化后的地震波场非常连续，几乎没有视觉上的差异，并且从误差剖面(图 2e)也可以看出规则化前后几乎没有误差，说明非均匀曲波变换方法规则化效果好。

采用均匀曲波变换方法对非均匀含噪数据(图 2c)直接去噪，对分解后的曲波系数进行软阈值处理，将保留下来的曲波系数进行曲波反变换，从而得到最终的去噪结果，如图 3a所示，其信噪比为4.45dB。图 3b为去除的噪声干扰，可以看出，其中含有部分有效波信号，主要原因是常规曲波变换不能直接处理非均匀采样数据，而是把非均匀采样数据当成均匀采样数据进行处理。根据噪声能量计算出噪声估计值σ=2.5，采用本文方法在非均匀采样数据规则化过程中同时去噪，结果如图 3c所示，其信噪比为5.11dB。去除的噪声如图 3d所示。可以看出，本文方法不仅可以将非均匀采样地震数据内插为均匀采样数据，还可以同时有效地去除随机噪声干扰，去噪后的地震同相轴更加连续、清晰，大幅度地提高了信噪比，且去除的噪声中几乎不含有效波能量，表明对于非均匀采样地震数据本文方法具有较好的去噪效果。

为了详细对比去噪效果，将原始理论地震数据、非均匀地震数据、含噪非均匀地震数据、常规曲波变换去噪结果和本文方法去噪结果分别进行局部放大，如图 4所示。可以看到非均匀地震数据显示在均匀网格上使有效波同相轴扭曲错动明显，再加上噪声的影响，这种扭曲现象更为严重，使有效波能量和噪声能量相互缠绕在一起，降低了整个地震数据的信噪比。从去噪结果的局部放大显示来看，由于常规曲波变换的去噪前提条件是均匀采样，所以在压制非均匀地震数据噪声时不能有效调整非均匀采样点位置，使去噪后的有效波同相轴仍然扭曲错动，并且噪声压制也不彻底。而本文方法去噪后的有效波与原始地震数据有效波非常接近，去噪后的同相轴光滑连续，这进一步说明本文方法在去除随机噪声的同时，也可以将非均匀地震数据调整为均匀采样数据。

为了对比不同噪声水平下本文方法的去噪效果，对原始理论数据加入标准差分别为0.075和0.100的高斯白噪声，然后采用非均匀傅里叶反变换获得非均匀地震数据，如图 5a和图 5b所示，信噪比分别为-6.99dB和-8.23dB。可以看出有效波信号被噪声淹没，难以识别，并且有效波同相轴局部扭曲，模糊不清，尤其是图 5b，几乎看不出有效信号。采用本文方法同时进行数据规则化和噪声压制，噪声估计值分别为σ=3.8和σ=5.1，结果如图 5c和图 5d所示。去噪后的信噪比分别为3.34dB和1.99dB，可以看出本文方法能够有效地去除非均匀采样数据中的随机噪声，并且将其内插为均匀采样，极大地提高了原始含噪数据的信噪比，但随着噪声能量大幅度增强，有效波损伤也会增加。

5 应用实例

图 6a为海上单炮地震数据，数据采集时设计道距为12.5m，然而由于海上拖缆的羽状漂移导致采集道距不均匀，道距范围为6~20m。图 6b为其局部放大，可以看出随机噪声较为发育，并且非均匀采样扭曲了部分有效波同相轴。采用本文方法进行去噪，分解的尺度为5，在第二个最粗尺度上的角度数为16，噪声估计值σ=8.5。图 6c为本文方法去噪结果，图 6d为局部放大显示，可以看出绝大部分噪声能量得到了有效压制，去噪后的地震数据同相轴较为连续，几乎没有损失有效波。图 7a和图 7b分别为原始单炮数据和本文方法去噪结果所对应的二维频谱，也可以看出本文方法去噪相对彻底，并且将非均匀采样地震数据内插为均匀采样。从去噪前、后振幅谱差(图 7c)可以看出，去除了大部分随机噪声，有效波的损伤较小，大幅提高了信噪比。

6 结论与认识

本文在多尺度多方向二维曲波变换基础上，提出了基于非均匀二维曲波变换和线性Bregman方法的地震数据随机噪声压制方法。由于非均匀曲波变换在正、反变换过程中具有能量无损性质，而这种性质特别适合在曲波域处理非均匀采样下的含噪地震数据，并且将其内插为均匀采样，在此过程中能有效地压制随机噪声。理论和实际资料的处理结果表明，该方法可以有效地压制地震数据中的随机噪声，且尽可能地保护了微弱的有效波信号，使反射波同相轴更加连续、清晰。

本文方法采用非均匀采样快速傅里叶变换，需要进行褶积和反褶积运算，计算速度远远低于均匀采样快速傅里叶变换，而且曲波变换实现过程复杂、冗余度高，因此本文所提出的非均匀曲波变换去噪方法计算效率相对较低。本文只进行了二维非均匀曲波变换，如推广到三维数据，则计算时间会进一步增加，需要发展其他快速算法提高本文方法的计算效率。