尺度自适应三维Shearlet变换地震随机噪声压制

程 浩 王德利 王恩德 付建飞 侯振隆

(①东北大学深部金属矿山安全开采教育部重点实验室，辽宁沈阳 110819；②吉林大学地球探测科学与技术学院，吉林长春130026)

0 引言

地震勘探是矿产资源勘查的一种重要手段[1-3]。强随机噪声的干扰直接影响地震资料处理和解释。随机噪声一般表现为无规律振动，由多种不确定因素形成。因此，去除随机噪声是地震数据弱信号提取与增强的重要环节。目前，常用的随机噪声压制方法包括频率域滤波、f-x域去噪、Radon域[4-6]去噪等。20世纪80年代，小波变换[7]作为一种多尺度变换，在图像分析、去噪等方面得到快速发展，能够稀疏表示一维信号，但不能简单地扩展到二维或高维信号。因此Ridgelet变换[8]、Curvelet变换[9-11]和Contourlet变换[12-13]作为多尺度几何分析方法应运而生，在二维和高维信号的稀疏表示和图像处理方面发挥了重要作用，但它们的数学结构复杂，而且变换域的某一尺度的系数矩阵尺度与原始数据不同。

Shearlet变换[14-17]作为一种新的稀疏变换方法，具有多尺度、多方向特性，且能最优地刻画信号的局部特征。与其他多尺度分析方法相比，Shearlet变换具有更简单的数学结构，对尺度和方向表征更细腻，采用更少的系数逼近曲线。对于一个N×N阶的图像，经过Shearlet变换后，各尺度方向的系数仍是N×N阶，且各方向系数仍然保留时间的维度。Shearlet变换可最优地稀疏表示信号，宜于处理包含纹理、轮廓和边缘信息的地震数据。含噪地震数据经Shearlet变换后，通常将有效信号映射到某些特定方向的Shearlet系数上，且对应的Shearlet系数值较大，而随机噪声分布在所有方向的Shearlet系数上，对应的Shearlet系数值较小。因此设置合理的阈值函数可以有效分离有效信号与随机噪声。

刘成明等[18]利用地震数据在二维Shearlet域的稀疏性去除随机噪声，获得了较好效果。程浩等[19]基于Shearlet变换的自适应阈值方法去噪，在提高地震数据信噪比的同时，最大限度地保留有效信号。上述方法都是针对二维单炮记录，并没有考虑到单炮记录间的相关性。

三维Shearlet变换[20]与二维Shearlet变换相比，兼顾了信号的空间相关性，在三维Shearlet域中能更稀疏地表示有效信号。与常用的三维Curvelet变换[21-23]相比，通过三维Shearlet变换得到某一尺度、某一方向的Shearlet系数与原始数据尺度相同，具有更明确的物理意义。段昶等[24]利用三维Shearlet变换融合处理磁共振图像，证明其稀疏性明显优于二维Shearlet变换。Negi等[25]分析了三维Shearlet变换的特点，并用于视频去噪，取得了较好的效果。

本文应用尺度自适应三维Shearlet变换压制多炮地震数据随机噪声，通过将多炮数据变换到三维Shearlet域，充分考虑单炮记录及其间的相关性，在三维Shearlet域更稀疏地表示地震数据。由于有效信号主要分布在低尺度，随机噪声分布在各个尺度，因此在硬阈值的基础上，结合尺度自适应因子压制随机噪声。再通过三维Shearlet反变换，得到去噪地震数据。数值模拟和实际多炮地震数据去噪结果表明，尺度自适应三维Shearlet变换的去噪效果优于二维Shearlet变换、不结合尺度自适应因子的三维Shearlet变换。

1 三维Shearlet变换

1.1 基本理论

Shearlet变换[15]作为一种稀疏变换方法，不仅受尺度和位置变化影响，还受剪切矩阵控制的正交变换影响。当空间维度参数D=3时，利用3个近似锥形区域P(D)和1个中心区域C函数构成Shearlet系统

(1)

图1 频率空间锥形区域P1(左)、P2(中)、P3(右)

(2)

并且使

(3)

满足

(4)

当D=1，2，3，l=(l1,l2)∈Z2,与三维Shearlet系统相关的锥形区域P(D)可以定义为

(5)

其中

(6)

(7)

1.2 Shearlet变换的稀疏性

基于信号分解的随机去噪方法的去噪效果依赖于方法对有效信号的逼近程度。为了展现Shearlet变换的稀疏特性，将其与Curvelet变换和Contourlet变换进行稀疏度对比。首先，利用不同的变换方法分解不含噪地震数据。然后，将变换域中的系数从大到小排列，保留不同百分比的系数重构，计算重构结果与原始数据的误差。图2为重构误差曲线。由图可见，保留的最大系数百分比从0.5%～30.0%，当保留4.0%的最大系数时，Shearlet变换的重构误差趋近于零，而Curvelet变换和Contourlet变换保留10.0%的最大系数时，重构误差才近乎为零，说明Shearlet变换具有更好的稀疏特性。

图2 重构误差曲线

1.3 Shearlet变换的局部刻画

Shearlet变换通过多尺度、多角度的特性，可以分解地震数据。由于地震同相轴是连续的，在局部同相轴具有一定的方向性，当同相轴的方向与某一Shearlet基的方向大致相同时，在Shearlet域对应较大Shearlet系数，同相轴的局部方向不同，在Shearlet域可被分解到不同的方向。由于随机噪声不具有方向性，在Shearlet域对应较小的Shearlet系数，因此通过设置合理的阈值可去除随机噪声。

图3为第50、125、200炮模拟单炮记录。利用二维Shearlet变换处理单炮记录，仅利用共炮点道集自身的特点压制随机噪声。但对于多炮记录，还可以利用共检波点方向的相关性。图4为多炮地震数据的空间切片。由图可见，各单炮记录间存在一定联系，类似于三维数据的单炮记录，即多炮地震数据可以看成一个伪三维数据。在变换域随机噪声压制过程中，不仅考虑单炮记录自身的特点，还兼顾单炮记录间的相关性，可获得更好的效果。

图5为三维Shearlet变换分解多炮地震数据的3组系数，可见由于Shearlet变换的方向性，地震同相轴根据自身的方向匹配到Shearlet域的不同方向，说明Shearlet变换可将地震数据分解到不同方向刻画。

图3 第50(左)、125(中)、200炮(右)模拟单炮记录

图4 多炮地震数据的空间切片

图5 三维Shearlet变换分解多炮地震数据的3组系数

2 三维Shearlet变换随机噪声压制

在随机噪声压制过程中，假设待处理的含噪地震数据为

f(t)=s(t)+n(t)

(8)

式中：s(t)为有效信号；n(t)为随机噪声。

三维Shearlet变换随机噪声压制过程与二维Shearlet变换类似，但后者仅考虑共炮点道集(单炮记录)中有效信号的分布特点，而前者不仅考虑共炮点道集有效信号的分布特点，还兼顾共检波点或共炮检距方向有效信号的分布特点，可更好地去除随机噪声，保护有效信号。三维Shearlet去噪主要包括以下步骤。

(1)通过三维Shearlet变换，对含噪地震数据进行多尺度、多角度剖分，获取对应的系数

(9)

式中：C(j,l,k)为参数尺度j、方向l、位置k对应的三维Shearlet系数；SH(·)为三维Shearlet正变换。

(2)在步骤(1)中，当信号的某一方向与对应某一角度的Shearlet基函数方向大致相同时，将获得较大的Shearlet系数，相反，当所选Shearlet基函数方向与信号方向偏差较大时，将获得较小的Shearlet系数。因为地震同相轴都具有一定方向性，所以在Shearlet域对应尺度、角度上将获得较大的Shearlet系数，而随机噪声不具有方向性，将获得相对较小的Shearlet系数。通过阈值函数对Shearlet系数进行阈值化处理，保留大于阈值的Shearlet系数，去除小于阈值的Shearlet系数，达到去除随机噪声的目的。常用硬阈值函数为

(10)

式中：M为常数变量，根据不同工区环境的噪声级别选取不同值，同一工区采集的多炮地震数据，由于所处环境基本一致，所以所选值相同；λ为某一尺度和角度的Shearlet系数均方差；σ为噪声标准差。

由于有效地震信号主要分布在相对低频的尺度，随机噪声则分布在各个尺度，因此为了尽可能地保留有效信号，在硬阈值函数中加入尺度自适应因子ε， 将式(10)变为尺度自适应阈值，即

(11)

其中

(12)

式中J为所选的尺度值。

(13)

为了体现所提方法的有效性，分别计算地震数据去噪前、后的信噪比(SNR)、峰值信噪比(Peak signal to noise ratio,PSNR)和均方差(Mean square error，MSE)进行对比[17]。

3 理论数据测试

利用模拟的多炮地震数据作为伪三维数据，验证尺度自适应三维Shearlet变换去噪方法。图6为数值模型及第125炮原始单炮记录、加噪单炮记录。由图可见：第125炮原始单炮记录浅部有效信号能量较强，底部有效信号能量较弱，对应两个能量较弱的反射同相轴(图6b箭头所示位置)；在模拟记录中添加一定程度的随机噪声，仍可以有效识别浅部能量较强的信号，但底部能量较弱的信号被随机噪声完全淹没(图6c箭头所示位置)而无法识别。

图6 数值模型(a)及第125炮原始单炮记录(b)、加噪单炮记录(c)利用有限差分法进行模拟，共251炮模拟地震数据，每炮记录包含251道和1000个样点，采样间隔为0.002s，采样时间为2s。经计算可知，添加随机噪声地震记录的SNR、PSNR和MSE分别为-1.2772dB、74.1516dB和2.5×10-3

分别利用二维(图7a)、三维(图7b)Shearlet变换和尺度自适应三维Shearlet变换(图7c)处理加噪多炮地震记录。可见：①二维Shearlet变换硬阈值去噪结果较充分地去除了随机噪声，图像清晰，但在去噪的同时损害了有效信号(图7a箭头位置)，还产生一定的边界效应；对应的差剖面具有清晰可见的有效信号(图8a)，SNR、PSNR和MSE分别为15.1687dB、89.9425dB和6.589×10-5(表1)，说明二维Shearlet硬阈值去噪方法具有一定的随机噪声压制能力，但一定程度上损害了能量较弱的有效信号。②三维Shearlet变换硬阈值去噪结果在去除随机噪声的同时，保留了底部能量较弱的有效信号，且由变换产生的边界效应也较弱，不影响数据的成像效果(图7b)；对应的差剖面基本不残留有效信号(图8b)，SNR、PSNR和MSE分别为19.4816dB、94.9103dB和2.0992×10-5(表1)。因此与二维Shearlet变换相比，三维Shearlet变换不但考虑了单炮记录的特点，还考虑了单炮记录间的相关性，在三维Shearlet域地震数据更稀疏，对随机噪声的去除更彻底。③与二维(图7a)、三维(图7b)Shearlet变换去噪结果相比，尺度自适应三维Shearlet变换阈值去噪结果对随机噪声的去除效果更好，充分保留了底层能量较弱的有效信号，由变换产生的边界效应更小，反射同相轴连续性更好、更清晰(图7c箭头所示位置)；对应的差剖面不存在有效信号(图8c)，SNR、PSNR和MSE分别为20.1058dB、95.5345dB和1.8181×10-5(表1)，在信噪比和保真度方面得到进一步提高。因此尺度自适应三维Shearlet变换去噪方法在去除随机噪声的同时，最大限度地保留了有效信号。

图7 二维(a)、三维(b)Shearlet变换和尺度自适应三维Shearlet变换(c)对图6c数据去噪结果所选尺度参数为4，对应的方向剖分参数为[1 1 2 2]

图8 图6c数据去噪结果与原始单炮记录的差剖面

表1 信噪比、峰值信噪比、均方差对比

分别抽取第125炮第40道(图9)和第90道(图10)的1.2～2.0s的记录进行振幅分析。由图9可见：原始有效信号淹没在含噪数据中(图9a)；三维Shearlet变换去噪方法大致在1.7s前去噪效果与保幅性较好，振幅基本与原始信号吻合，在1.7s以后振幅基本归零，与有效信号差异较大(图9b)；尺度自适应三维Shearlet变换去噪方法的整体去噪效果与保幅性都很好，振幅与原始信号吻合(图9c)。由图10可见：原始有效信号淹没在含噪数据中(图10a)；三维Shearlet变换去噪方法在1.5s前振幅拟合相对较好，但是幅值存在一定差异，大致在1.75s又出现一定振幅响应，差异较大，其他幅值基本为零(图10b)；尺度自适应三维Shearlet变换去噪方法在1.5s前振幅拟合较好，幅值差异相对较小，较好地拟合了1.65～1.80s的有效振幅响应(图10c)。由图9、图10可知：由于存在炮检距变化，第40道的有效信号幅值相对较大，去噪效果较好，振幅拟合程度高，保幅性好(图9)；第90道有效信号的幅值相对较小，去噪效果和振幅拟合程度受到影响(图10)。

图9 第125炮第40道数据振幅分析

图10 第125炮第90道数据振幅分析

通过上述分析可知，尺度自适应三维Shearlet变换去噪方法不但考虑了单炮记录自身的特点，还顾及了单炮记录间的相关性，充分地利用了地震数据在三维Shearlet域的稀疏性。同时，结合有效信号与随机噪声在三维Shearlet域的分布特点，利用尺度自适应因子最大限度地保留有效信号，提高信噪比和保真度。

4 实际数据测试

选取A区实际地震数据测试尺度自适应三维Shearlet变换的去噪效果。图11为第123炮记录及去噪结果。由图可见：①单炮记录上部同相轴较清晰，能量最强，连续性好，底部有效信号几乎被随机噪声淹没，同相轴相对能量较弱，连续性差，不易识别，影响层位和断层解释(图11a)。②二维Shearlet变换较好地压制了随机噪声，去噪结果上部同相轴连续性变好，整体清晰，但底部有效信号被损害，同相轴连续性差，严重影响底部层位解释(图11b)；在对应的差剖面上可较清晰地看到有效信号的趋势(图12a)。③三维Shearlet变换充分地压制了随机噪声，去噪结果上部反射同相轴连续性更好，去噪结果较好地保留了底部信号，底部反射同相轴连续性较好，并保留了部分能量较弱的有效信号(图11c)；在对应的差剖面上几乎不残留有效信号(图12b)。④与前两种去噪方法相比，尺度自适应三维Shearlet变换去噪结果上部反射同相轴连续性最好、最清晰，能量最强，并最大限度地保留了底部有效信号，有效改善了反射同相轴的连续性(图11d)；在对应的差剖面上基本上不残留有效信号(图12c)。

图11 第123炮记录及去噪结果

图12 图11a数据去噪结果与图11a的差剖面

5 结束语

在传统稀疏变换去噪方法的基础上，笔者应用尺度自适应三维Shearlet变换去除地震数据噪声。尺度自适应三维Shearlet变换去噪方法不但考虑了单炮记录自身的特点，还顾及了单炮记录间的相关性，充分地利用了地震数据在三维Shearlet域的稀疏性。同时，结合有效信号与随机噪声在三维Shearlet域的分布特点，利用尺度自适应因子，最大限度地保留有效信号，提高信噪比和保真度。

尚需指出，由于对多炮地震数据进行三维Shearlet变换，对应的任意尺度、方向上的系数矩阵规模均与原始数据一致，所以本文方法对服务器的内存有一定要求。此外，通过对比理论数据单道振幅可知，由于炮检距的变化，有效信号的幅值在横向产生变化，因此所提方法可能对幅值相对较小的有效信号产生损害。