0 引言

压制地震记录中的随机噪声是提高地震勘探资料信噪比的重要环节(Warpinski，2009)，固体矿产赋存地带获取的地震资料的质量较差(Leonard and Kennett, 1999；张华等，2011)，含有大量的强随机噪声，信噪比较低(Zhang et al., 1999；马继涛等，2016).传统的滤波方法，如匹配滤波、F-K域消除、双曲Radon域变换、稀疏表示等去噪方法(林红波等，2011)在有效地压制地震资料随机噪声中以往起到较好效果，但是随着地震记录信噪比的降低，对随机噪声压制受到限制.匹配滤波是需要预先得到同相轴作为参考轴并进行地震噪声压制和识别定位的方法，在低信噪比环境下其滤波效果并不显著(Eisner et al., 2008).F-K滤波利用了有效信号与噪声在视速度上的差异进行二维滤波(南方舟等，2018)，从视速度上可以很好地在时空域将有效信号和噪声进行区分，不过有效信号畸变甚至假轴的问题仍需要得到进一步的改进(Song et al., 2014).而稀疏表示的去噪方法(Rodriguez et al., 2011)以及双曲Radon域去噪(Sabbione et al., 2013)也存在着类似的问题，因此非常有必要提出一种新的方法来压制低信噪比地震信号中的随机噪声.相对于以上方法，多尺度集合分析的小波变化在引入到地震信号处理方面以后得到了很好的效果(张华等，2017；邵婕等，2016)，如Curvelet变换(Eslami and Radha, 2003)、Contourlet变换等(Romberg et al., 2003).Shearlet变换是一种结合Curvelet和Contourlet优点的新型多尺度几何分析方法.2005年Labate和Guo(Guo and Labate, 2007)等将复合小波理论和几何多尺度分析相结合，通过特殊形式的具有合成膨胀的仿射系数构造了一种接近最优的多维函数的稀疏表示法——Shearlet变换(Houska and William, 2012), 它比Curvelet和Contourlet方向性更加敏感而且稀疏表示特性更好，运算的速率也更快(Houska and William, 2012；Lim，2010)，非常适合对地震信号进行处理.但是传统的Shearlet变换去噪采用的是硬阈值，在压制随机噪声的同时也滤除了相当一部分的有效信号的成分，使得去噪后的地震资料出现虚假轴，为了解决这个问题，我们提出高阶加权阈值函数.高阶加权阈值函数整体上连续性较好，克服了硬阈值函数存在剧烈变化的缺点以及软阈值在处理较大Shearlet系数总存在恒定偏差(Donoho，1995)的问题，同时保留了传统的软硬阈值函数的优点.实验显示运用本文方法对信噪比为-10 dB的模拟地震资料处理，不仅有效地压制随机噪声，还很好地保留了有效信号的幅度.与Shearlet硬阈值去噪以及小波变换去噪做对比，压制噪声保留有效信号幅度上均具有明显的优势.在本文的实际实验部分，我们选取丘陵地带的实际地震资料作为处理对象，本文的方法也取得较好的去噪效果.

1 Shearlet变换基本原理

与Curvelet不同，Shearlet具有简单的数学结构，它可以通过对一个函数进行伸缩、平移、旋转而生成一组基函数(Kutyniok and Labate, 2007)，这使其可以和多分辨分析关联起来.Shearlet变换不仅具有与Curvelet相同的近似最优的非线性误差逼近阶，而且在频率空间中Shearlet是逐层细分的，这使其具有更好的表示性能.Kutynoik和Labate定义由ψ∈L²(R²)产生的连续Shearlet系统(Zhao et al., 2016)为

(1)

其中，L²(R²)表示二维可积实平面，ψ表示Shearlet母函数，是各向异性膨胀矩阵，S_s= 是剪切矩阵.函数f∈L²(R²)的连续Shearlet变换定义如下：

(2)

连续Shearlet变换将函数f以尺度a、方向s和位置t投影到函数ψ_ast上，这与连续Wavelet变换的经典定义相似(凌振宝等，2016)，其中函数f投影到一个由尺度a和位置t构成的函数ψ_at上，类似连续Wavelet，一个连续的Shearlet的ψ∈L²(R²)必须满足以下条件(Guo and Labate, 2013):

(3)

其中，ξ₁和ξ₂表示二维频域坐标.f∈L²(R²)连续Shearlet反变换定义如下:

(4)

为了方便计算机处理，Kutyniok和Labate提出离散Shearlet系统如下(Kutynoik and Lim, 2011)：

(5)

其中ψ∈L²(R²).式(5)可以由式(1)通过在离散集上对连续参数a、s和t采样而得.

2 高阶加权阈值函数原理

传统的Shearlet变换去噪的实质是将含噪地震信号经过Shearlet变换得到相应的Shearlet系数，然后通过设置该系数的阈值去消减噪声.

2.1 对阈值加权处理方法

当Shearlet基函数方向与地震信号方向大致相同的时候，经过变换得到较大的Shearlet系数，反之，则得到较小的Shearlet系数.在实际的地震资料中，随机噪声不具有方向性，因此将含噪信号变换到Shearlet域，有效信号的Shearlet系数较大，而噪声的Shearlet系数比较小.这样我们设定一个阈值，采用硬阈值函数来除去较小的Shearlet系数，再对新的系数进行Shearlet反变换得到去噪后的地震信号.

传统的阈值函数主要包括硬阈值函数与软阈值函数(曹静杰等，2015)，分别如式(6)和式(7)所示：

(6)

(7)

(8)

c(i, j)表示Shearlet系数，c′(i, j)表示经过阈值函数处理的Shearlet系数，sgn()表示符号函数，σ表示噪声的标准差，N表示地震资料的大小.

图 1中可以看出硬阈值函数存在剧烈的变化，处理后的Shearlet系数相关性较差，因此反变换后得到的地震信号会产生震荡.而软阈值函数得到的Shearlet系数c′(i, j)总是小于原来的Shearlet系数c(i, j)，存在固有的误差，因此会影响重构后的信号与真实信号的逼近程度.两种传统的阈值在处理低信噪比的地震资料已经不能取得令人满意的效果.为了克服硬阈值函数存在剧烈的变化和软阈值存在固有误差的问题，本文提出高阶加权阈值函数.首先我们引入加权因子W.

(9)

其中，c(i, j)表示待处理的Shearlet系数，σ表示噪声的标准差，N表示地震资料的长度，n为加权系数，可以取1到+∞的任何数.我们利用加权因子构建高阶加权阈值函数，如式(10)所示：

(10)

式(10)中c(i, j)表示Shearlet系数，c′(i, j)表示经过高阶加权阈值函数处理后的Shearlet系数，w表示上文提到的加权因子.由图 1可以看出高阶加权阈值函数连续性比传统的硬阈值函数要好.

2.2 误差理论分析

对于高阶加权阈值函数的误差，以c(i, j)≥λ为例，由公式(7)可以看出软阈值函数的恒定误差为λ，无论Shearlet系数大小是多少，误差总为λ.现在我们将高阶加权阈值函数化简得到公式(11)：

(11)

由公式(11)可以看出高阶加权阈值函数的误差是，由于c(i, j)≥λ，所以1，故高阶加权阈值函数的误差小于软阈值函数的误差.而且当c(i, j)越大时，高阶加权阈值函数的误差越小，这就巧妙的解决了软阈值函数在处理较大Shearlet系数存在较大恒定误差的问题.往往有效信号的Shearlet系数较大，这样运用高阶加权阈值函数就可以大大减小有效信号Shearlet系数的误差，使得反变换后得到的去噪信号更加接近真实的信号.下面我们来说明n如何取值适合处理地震信号的Shearlet系数.可以看出当n=1的时候，高阶加权阈值函数就是传统的软阈值函数；当n=+∞的时候，高阶加权阈值函数接近传统的硬阈值函数.可以看出经过高阶加权阈值函数处理的Shearlet系数相比软阈值函数处理后的固有误差较小，而且曲线的连续性更好，克服了硬阈值函数的缺点.为了更加直观我们给出n=2, 4, 6时的高阶加权阈值函数的示意图，如图 1所示.可以看出n越大，高阶加权阈值函数趋近于传统的硬阈值函数，连续性较差；n越小，高阶加权阈值函数趋近于传统的软阈值函数，固有误差较大.经过我们多次试验，当n=4时的高阶加权阈值函数最适合对地震信号的Shearlet系数进行处理.

2.3 本文算法大致流程

(a) 分解：对需要处理的低信噪比地震信号进行Shearlet变换.

(b) 新的阈值函数处理：在Shearlet域利用高阶加权阈值函数对系数进行阈值收缩处理.

(c) 反变换：对处理后的Shearlet系数进行Shearlet反变换，得到去噪后的地震资料.

3 仿真实验及结果 3.1 模拟实验

首先把本文的方法应用到合成的地震记录中，运用雷克子波构建模拟的地震记录(Fletcher et al., 2002)，模拟含噪地震记录如图 2a所示，其采样频率为500 Hz, 模拟信号中含有三条主频为25 Hz的同相轴.我们在模拟的原始信号中加入高斯白噪声，使该地震记录的信噪比为-10 dB，如图 2b所示.从图 2b可以看出地震信号中的同相轴被强噪声淹没，很难识别出有效信号.实验中，我们使用传统的基于小波变换的去噪算法和Shearlet硬阈值算法作为对比实验，结果如图 2c、图 2d和图 2e所示.从图 2c和图 2d的去噪结果中可以看到基于小波变换的去噪算法和Shearlet硬阈值算法, 可以滤除大部分的随机噪声，但是噪声去除的不够彻底，同相轴的幅值也严重衰减，因此效果并不理想.相比之下，图 2e基于本文方法的去噪结果地震有效信号同相轴恢复的更加彻底，噪声也被更加有效地移除.尤其在图中的蓝色标识区，可以很明显的看出本文方法优于另外的两种方法.

为了更好的观察文中各种方法的去噪结果，我们随机取出第18道进行时域波形和频谱对比，对比的结果如图 3a和图 3b所示.通过对比图 3中的处理结果不难发现，在波形上，本文方法在子波的波峰和波谷都可以进行有效的恢复；在频谱上，本文方法在低频追踪和高频压制上都更加接近原始纯净信号.

最后我们定量的验证本文算法的有效性，我们在图 2a所示的原始纯净地震信号中加入不同强度的高斯白噪声，通过对比信噪比来验证本文算法的有效性.我们用三种方法处理了信噪比为-5 dB，-10 dB和-15 dB的含噪信号，去噪后的信噪比(SNR)和均方误差(MSE)如表 1、表 2所示.从表中的数据可以看到在不同的信噪比水平下，三种方法去噪后的信噪比都有所提升.其中基于小波的去噪算法去噪后信噪比最低，而本文的方法去噪后具有较高的信噪比.尤其是在处理信噪比较低的地震信号时，本文算法的优势更加明显.同时相对于其他两种算法，本文方法去噪后的地震数据均方误差也更小.从SNR和MSE这两个指标综合比较后可以看出，本文算法对地震信号的去噪效果更好.

3.2 野外地震记录处理

为了验证本文方法的实际应用效果，我们将该方法应用于云南某山地丘陵地带实际地震记录中，野外数据采集60次覆盖，道间距40 m，采样间隔1 ms.我们对共80道实际地震记录中截取1024点进行处理.由于丘陵地带的地形复杂、树木较多、风变化大等因素，丘陵地带的随机噪声具有横向变化不稳定特性，如图 4a右侧红框所示.同时丘陵地带的随机噪声较强，这些都给丘陵地带的随机噪声的压制带了巨大的困难.我们分别采用基于小波变换的阈值去噪算法、Shearlet硬阈值算法和本文算法对含噪信号进行滤波处理，去噪结果如图 4b—d所示.通过对比图中红色标识区内的地震信号可以发现，本文算法处理后的实际地震信号同相轴恢复的更加清晰连续，同时也更加易于识别.

为了更好的对比去噪结果，我们将图 4中平行四边形框部分区域进行了局部放大，如图 5所示，从图中可以看出，本文算法噪声压制效果理想，而且原来杂乱无章的信号变得清晰连贯，信号幅度保持较好.也验证了本文算法在处理实际地震资料时，比其他两种方法更加有效，可以更好的去除噪声，并恢复同相轴.

4 结论与讨论

本文提出了一种基于高阶加权阈值函数的Shearlet变换去噪算法，提出的高阶加权阈值函数整体连续性较好，克服了传统阈值的缺点.经过仿真模拟与野外地震资料的处理，结果表明，用本文提出算法在压制随机噪声的同时也很好的保留了有效信号的幅度，而且信号同相轴更清晰连续，原本淹没的同相轴也显现出来，增加了可以从记录中获得的信息；和小波变换去噪算法、Shearlet硬阈值去噪算法的对比，在压制低信噪比地震信号中的随机噪声、保留有效信号的幅度均具有明显的优势.

本文就山地丘陵地带随机噪声进行的实际消噪处理，取得明显的效果；一种好的滤波方法技术需要有广泛的适用性，这需要把本文方法应用于沙漠地带、黄土塬地带等不同地面环境的去噪分析研究与处理；进一步得到高阶加权阈值的广泛适用性.