0 引言

地震资料的随机噪声通常源于仪器、环境及信号采集时的诸多因素，降低了地震资料的信噪比和分辨率^[1]。压制随机噪声的主要目的是分离有效信号与噪声及保持有效信号的波形^[2]。根据有效信号和噪声在视速度、频率等方面的差异，相继出现了F-X反褶积、Radon变换、小波分析、模态分解等去噪方法^[3-6]。近年来，小波变换以优越的时频聚焦特性迅速从数学、信号处理拓展到各个领域，但其在一维信号的优良特性并不适用于高维，对方向的不敏感使信号表征存在局限^[7-8]。

随着多尺度几何分析方法的不断发展、成熟，逐渐用于地震勘探领域，为高保真去噪提供了丰富的数学基础及新思路。多尺度几何分析也称多分辨率分析^[9]，最早由Rosenfeld等提出，具有多分辨率、局部化、多方向性及各向异性等特征，可更好地表征信号的局部特征，并实现信号的稀疏表示，在去噪和压缩方面具有明显优势。近年来相继出现了Ridgelet^[10]、Curvelet^[11]、Contourlet^[12]、Shearlet^[13]等变换的多尺度几何分析方法，其中Shearlet变换以优越的多分辨率特性在中国发展迅速，已广泛用于图像处理领域，并取得了良好效果^[14-17]。同样，Shearlet变换的多分辨率、多方向性等特性在地震数据去噪中也获得较好的应用效果^[18-20]。

Shearlet变换通过对变换系数设定阈值约束分离信噪，其中阈值函数至关重要。目前常用的阈值函数对多尺度、多方向的Shearlet变换局限性很大，会造成噪声压制不彻底或损失有效信号等问题。为此，本文将信噪比与阈值函数有机关联，在不同尺度下基于信噪比约束自适应求取阈值，实现自适应去噪。在保证有效信号不受损失的情况下，恢复被噪声掩盖的弱信号，有效地改善去噪效果。

1 方法原理 1.1 Shearlet变换

Shearlet变换^[21]是通过特定形式的合成膨胀仿射系统构造的一种多尺度几何分析工具，是小波变换的高维扩展，具有接近最优的非线性误差逼近性能，数学结构严谨，计算复杂度低。Shearlet变换由拉普拉斯金字塔变换和方向滤波器组构成。设

$ \boldsymbol{A}=\left(\begin{array}{ll}{a} & {0} \\ {0} & {\sqrt{a}}\end{array}\right) $

$ \boldsymbol{B}=\left(\begin{array}{cc}{1} & {-s} \\ {0} & {1}\end{array}\right) $

式中：A为尺度矩阵，根据抛物尺度准则划分尺度，其中a＞0为尺度参数；B为剪切矩阵，进行方向剖分，其中s(整数)为剪切参数。令j、l、k分别代表尺度、方向和系数位置序号，则由j、l、k定义的实数域连续可积函数为

$ \psi_{j, l, k}(x)=|\operatorname{det} \boldsymbol{A}|^{j/2} \psi\left(\boldsymbol{B}^{l} \boldsymbol{A}^{j} x-k\right) $

(1)

x为自变量。若使ψ满足Parseval框架，则信号函数f∈L²(R²)的连续Shearlet变换为

$ \mathrm{SH}_{f}(j, k, l)=\left\langle f, \psi_{j, l, k}\right\rangle $

(2)

通过对参量j、l、k的控制实现不同尺度、不同方向的剖分，该过程相当于由平移、伸缩、旋转后的基函数逼近信号函数。在频域中，每一对Shearlet基函数都有一对楔形支撑区间，随尺度因子变小而变窄(图 1)。

1.2 利用Shearlet变换压制随机噪声 1.2.1 Shearlet变换去噪原理

根据Shearlet变换的多尺度、多方向性特点，地震数据经过Shearlet变换会得到一系列不同尺度、不同方向的Shearlet系数。若Shearlet基函数的方向越逼近有效信号，则Shearlet系数越大；若基函数的方向与信号方向偏差越大，则Shearlet系数越小。由于随机噪声不具有方向性，所以经Shearlet变换后所得系数较小。利用阈值函数去掉较小的Shearlet系数，保留较大的部分，就能压制随机噪声，再进行Shearlet反变换得到去噪后的记录。获得Shearlet系数的流程(图 2)为：①利用拉普拉斯金字塔变换多尺度剖分信号，得到低、高频子带；②进行多方向分解，并在高频子带的伪极坐标进行傅里叶变换；③用Meyer函数滤波步骤②得到的结果；④对低频子带循环进行步骤①~③；⑤在笛卡尔坐标系对滤波数据进行快速傅里叶反变换，得到Shearlet变换系数。

1.2.2 自适应阈值函数

在阈值类去噪方法中，阈值选取和阈值函数设计直接影响去噪效果。传统的阈值选取方法是对所有变换域系数使用统一阈值，但对Shearlet变换而言，各尺度、各方向的有效信号和噪声均存在差异，因此全局硬阈值存在一定局限性^[22-23]。局部阈值则根据一定范围内的系数分布情况确定。本文在局部阈值的基础上改进贝叶斯阈值，形成一种适用于Shearlet变换的自适应阈值函数。

贝叶斯阈值是Chang等^[24]利用小波系数的广义高斯分布提出的一种具有自适应性的阈值函数，其公式为

$ T=\frac{\sigma_\text{n}^{2}}{\sigma_\text{s}} $

(3)

式中：T为阈值；σ_n²为噪声方差；σ_s为信号标准差。在式(3)引入一个与各尺度记录信噪比相关的权值系数λ，形成新的阈值函数

$ T_{\lambda}=\lambda \frac{\sigma_\text{n}^{2}}{\sigma_\text{s}} $

(4)

式中σ_n由中值估计法得到

$ \sigma_{\mathrm{n}}(j, l)=\frac{f_{\mathrm{Median}}(|S(j, l, k)|)}{0.6745} $

(5)

f_Median代表求取中值，σ_n和σ_s满足σ_d²=σ_s²+σ_n²，且

$ \sigma_{\mathrm{s}}(j, l)=\sqrt{\max \left[\sigma_{\mathrm{d}}^{2}(j, l)-\sigma_{\mathrm{n}}^{2}(j, l), 0\right]} $

(6)

$ \sigma_{\mathrm{d}}^{2}(j, l)=\frac{1}{n} \sum\limits_{k=1}^{n} S^{2}(j, l, k) $

(7)

σ_d为Shearlet系数的均方根；S(j, l, k)为Shearlet系数；n为某尺度、某方向Shearlet系数的个数。式(4)由不同尺度、不同方向的σ_n、σ_s及λ_j共同确定，λ_j与信噪比相关

$ \lambda_{j}=\frac{1}{\operatorname{SNR}(j)} $

(8)

SNR(j)为不同尺度Shearlet系数单独构成的地震信噪比，由相邻地震道互相关求得

$ \operatorname{SNR}(j)=\frac{N \sum\limits_{i=1}^{N-1} Q_{i,i+1}}{(N-1) \sum\limits_{i=1}^{N} Q_{i, i}(0)-N \sum\limits_{i=1}^{N-1} Q_{i, i+1}} $

(9)

式中：Q_{i, i+1}为第i道和第i+1道互相关函数最大值；Q_{i, i}(0)为第i道自相关函数最大值；N为道数。

确定阈值后即可调节Shearlet系数

$ S_{\mathrm{new}}=\left\{\begin{array}{ll}{S} & {|S| \geqslant T_{\lambda}} \\ {0} & {|S|<T_{\lambda}}\end{array}\right. $

(10)

式中S、S_new分别为阈值处理前、后的Shearlet系数。自适应阈值函数方法通过将信噪比作为阈值设定的因素，即不同信噪比的权值系数不同，可以自适应求取不同尺度阈值，最大限度地改善去噪效果，避免有效信号损失。

2 模型试算与实际资料应用 2.1 模型试算

为测试Shearlet变换自适应阈值方法(下称本文方法)的去噪效果，首先进行模拟数据试算。图 3展示了模拟单炮记录去噪效果。由图可见：采用全局硬阈值方式的Curvelet变换方法虽然去除了大部分高斯随机噪声，但仍有部分残留，效果不佳(图 3c)；Shearlet变换全局阈值方法去噪效果较明显，基本压制了高斯随机噪声(图 3d)；本文方法去噪结果的信噪比进一步提高(图 3e)。

表 1展示了不同方法去噪前、后信噪比，可以看出本文方法的去噪效果更好。图 4为不同方法去噪前、后的残差剖面。由图可见：采用Curvelet变换全局阈值(图 4a)、Shearlet变换全局阈值(图 4b)去噪方法在深度较大的位置(红圈区域)造成有效信号损失，信号的保真度较低；采用本文方法去噪的残差剖面上不存在有效信号(图 4c)，去噪结果(图 3e)与模拟单炮记录(图 3a)基本一致。

2.2 实际资料应用

为进一步检验本文方法的有效性，选取A区陆地三维地震资料测试随机噪声压制效果。该区的偏移地震数据虽经过初步处理，但仍含有较强随机噪声，致使目标层局部有效弱信号被噪声淹没。图 5为Inline 100629剖面Shearlet变换不同尺度的地震记录。由图可见，不同尺度的有效信号与噪声分布存在差异(即信噪比不同)，因此采用传统的阈值方法去噪存在局限性。截取Inline 100629剖面中深部区域(图 6a的黑框区域)进行去噪，图 6为实际地震记录去噪效果。由图可见：经各种方法去噪后随机噪声得到有效压制、信噪比得到提高，其中Shearlet变换全局阈值(图 6d)较Curvelet变换全局阈值(图 6c)的去噪效果好，局部弱信号得到恢复(红圈区域)，但部分有效信号仍受噪声干扰，去噪效果不佳；经本文方法去噪后，弱信号更清晰、连续，明显提高了信噪比(图 6e)。图 7为Shearlet变换信噪比改善程度—阈值曲线。由图可见：传统全局硬阈值求取方法采用黑色竖线与各尺度曲线交点处的值作为阈值，没有考虑各尺度的信噪特征；图中4条曲线的标注点对应的纵坐标数值为本文方法确定的阈值。

本文方法通过求取各尺度分量的信噪比，并将其作为约束条件，获得符合多尺度、多分辨率特征的阈值，去噪记录保真度更高。

为了更好地体现本文方法的自适应性，图 8给出了Inline 100629剖面不同方法去噪前、后的残差剖面。由图可见：采用Curvelet变换全局阈值(图 8a)、Shearlet变换全局阈值(图 8b)去噪方法虽然有效地压制了噪声，但由于受各尺度阈值设定的单一性限制，去噪程度基本一致；本文方法处理的残差剖面中出现较明显的弱能量区(红圈区域)，这是由于目标层有效信号区域的初始信噪比较高(红圈区域中与有效信号混在一起的随机噪声已得到初步压制)所致，意味着采用本文方法的去噪效果不明显，充分体现了去噪的自适应特点。

表 2为不同区域去噪前、后信噪比。由表可见，原始剖面中信噪比不同的区域使用本文方法去噪后信噪比基本相同，这也说明自适应阈值函数通过计算局部信噪比约束去噪效果，可避免有效信号的损失，使剖面质量明显提高。

3 结束语

Shearlet变换作为一种多尺度几何分析新方法，较小波变换等传统信号分析方法具有更多优良特性，且运算效率高。另外，针对地震数据去噪过程中传统全局阈值方法的局限性，本文通过改进自适应阈值函数压制随机噪声，在模型试算与实际资料应用中取得了满意的去噪效果，验证了方法的可行性和有效性。

尚需指出，目前Shearlet变换自适应阈值方法仅利用了Shearlet变换的多方向、多分辨率特性，其优良的稀疏特性同样值得进一步探索并用于地震数据处理领域，如尝试将Shearlet变换稀疏性与压缩感知理论结合，以期得到更好的去噪或拓频方法等。

在模型试算与实际地震资料处理中使用了中国石油集团东方地球物理公司提供的GeoEast软件，在此表示感谢。