Shearlet域尺度角度自适应深反射地震数据随机噪声压制方法

2022-06-16王通刘建勋王兴宇李广才田密

物探与化探 2022年3期

王通,刘建勋,王兴宇,李广才，田密

(1.国家现代地质勘查工程技术研究中心，河北廊坊 065000；2.中国地质科学院地球物理地球化学勘查研究所，河北廊坊 065000)

0 引言

在深地探测领域中，深反射地震技术被认为是探测岩石圈、解决深部地质问题最为准确有效的方法之一[1]，深反射地震探测方法采用大药量激发，长时间观测，地震波能够带来莫霍面及以下的深部信息，深反射地震数据能够对深部地质结构实现有效的成像。

深反射地震接收来自地层深部的反射信号，由于大地滤波作用使得反射信号能量很弱，受背景随机噪声干扰严重，一些微弱的构造响应会被随机噪声淹没，对深反射地震数据处理和解释带来很大影响。深反射地震数据的处理结果直接关系到地震剖面的整体成像质量，高分辨率地震数据处理技术对深反射地震探测至关重要[2-3]。朱小三等[4]指出怎样除去深部信号中噪声的同时又能最大限度地保留有效信号的能量是深反射地震资料去噪的关键。

通过对深反射地震数据的处理发现，常规处理(异常振幅噪声衰减、局部异常振幅压制、相干噪声和非相干噪声压制、球面扩散补偿、地表一致性补偿、地震子波反褶积等)后深反射地震数据中大部分噪声能够得到很好的压制，但是来自莫霍面及以下构造的弱反射信号不能有效地从背景噪声中呈现出来。传统的随机噪声去除方法(频率域滤波、f-x域去噪、Radon域去噪等)不能满足深部弱信号高精度成像的需求。

20世纪80年代，小波变换(Wavelet)[5-6]作为一种多尺度变换工具，在去噪和图形处理领域得到了快速的发展，利用小波变换方法分尺度对地震数据中随机噪声进行压制取得了很好的应用效果。程浩等[7]采用基于小波变换的自适应阈值，利用分层自适应因子，对矿山微震信号中所包含的随机噪声进行压制。罗红梅等[8]采用基于小波变换的经验模态EMD分解的地震记录重建方法，较好地提高了弱信号的能量信息。陈亮等[9]提出一种基于神经网络改进小波的地震数据随机噪声去除方法，采用级联BP神经网络模型提取出多类别随机噪声信号，实现地震数据的随机信号压制。基于小波变换的随机噪声压制方法实现了噪声阈值的尺度自适应，但是还不满足地震数据多方向性的需求[10-12]。

在Wavelet多尺度变换的基础上，多尺度多方向的Curvelet变换应运而生，有效的应用到地震数据随机噪声压制中。刘磊等[13]、金丹等[14]、何柯等[15]利用Curvelet域的多种阈值方法对二维地震数据实现了弱信号提取和随机噪声压制，张华等[16]采用曲波变换对地震数据抽取时间切片，实现了三维地震数据随机噪声压制。孙苗苗等[17]采用基于曲波域稀疏约束的OVT域地震数据去噪方法，有效提取了OVT域地震数据的中、深反射层的弱有效信号，在压制强随机噪声的同时减少了弱有效信号的损失。但是，Curvelet变换的数学结构复杂，且在变换的某一尺度角度域中系数矩阵尺度与原始数据不一致。

相较于Curvelet变换，Shearlet变换[18]数学结构更为简单，同样具有多尺度、多方向特性，利用更少的系数逼近曲线。对于一个N×N阶的图像，在Shearlet域中各尺度方向的数据体仍是N×N阶，且各方向系数仍然保留时间的维度。程浩等[19]采用尺度自适应三维Shearlet变换对地震数据中随机噪声进行压制。童思友等[20]在Shearlet域局部阈值的基础上改进贝叶斯阈值，将信噪比与阈值函数有机关联，并将信噪比作为阈值设定的因素，实现了地震资料中随机噪声的压制。薛琳等[21]结合Shearlet变换多尺度、多方向特性，提出一种随尺度和方向同时自适应变化的阈值实现地震数据随机噪声的压制。以上对Shearlet域内随机噪声估计的参数过于单一，对于深层弱反射信号下随机噪声压制效果欠佳。

本文应用结合信噪比(SNR)、L2范数、随机噪声残差的Shearlet域尺度角度自适应阈值方法对深反射地震数据中随机噪声进行压制，将地震数据变换到Shearlet域中，由于地震数据具有方向性，而随机噪声是无规律分布的，所以同一尺度下不同角度上有效信号与随机噪声比存在微弱差异，在此基础上利用信噪比(SNR)、随机噪声残差和Shearlet系数L2范数对每个角度域数据进行自适应阈值求取，实现更为精准的随尺度角度自适应变化的随机噪声压制方法。

1 Shearlet变换

Shearlet变换是基于合成小波理论，结合仿射系统和多尺度几何分析，带有合成膨胀系统的几何变换方法，是一种多尺度几何分析工具，具有最优的非线性误差逼近性能，数学结构严谨，计算复杂度低。连续的仿射系统在L2(Rn)具有以下形式[18]：

{TtDMψ,M∈G,t∈Rn}

(1)

式中ψ∈L2(Rn)，Tt是平移算子，并且定义：Ttf(x)=f(x-t)，DM表示伸缩算子，定义：DMf(x)=|detM|-1/2f(M-1x)，G⊂GLn(R)是一组矩阵，当n=2，G是含有2个参数的扩张组：

(2)

矩阵Mas可以表示为Mas=SsAa，式中S为剪切矩阵，A是尺度矩阵：

(3)

式中：s为剪切参数，a为尺度参数。令i,j,k分别代表尺度、方向和系数位置序号，则由i,j,k定义的实数域连续可积函数为：

φi,j,k(x)=|detA|i/2φ(SjAkx-k) ，

(4)

式中:x为自变量。若φ满足Parserval框架，则对于任何一个平方可积的函数f的Shearlet变换表示为：

(5)

式中：C为Shearlet系数矩阵；SH{·}为Shearlet变换；〈·,·〉表示内积；f为含噪信号。

同样可以对Shearlet系数矩阵C进行反变换重构原始函数：

f=SH-1{C} ，

(6)

式中，SH-1表示Shearlet逆变换。

2 Shearlet域自适应阈值

假设地震记录中有效信号为P0(t)，随机噪声为n(t)，则观测到的包含随机噪声的地震数据为：

P(t)=P0(t)+n(t) 。

(7)

传统阈值法压制随机噪声是通过计算地震数据，筛选出一个随机噪声强度相应的阈值参数Th，这里Th是一个常数。利用阈值参数Th对地震数据分选，大于Th的数据被视为有效信号进行保留，小于Th的数据被视为随机噪声进行压制，最终获得去除随机噪声后的地震数据。

Shearlet是一种稀疏变换方法，地震数据经过Shearlet变换后得到分尺度、分方向表示稀疏的Shearlet系数。Shearlet变换具有多尺度多角度特性，在特定的尺度角度中，当地震数据同相轴走向与Shearlet变换基函数的方向大致相同时，地震信号就对应着较大的Shearlet系数；当地震数据同相轴走向与Shearlet变换基函数的方向相差较大时，相应的Shearlet系数较小。在实际地震数据中，有效的地震信号通常具有方向性，对应着较大的Shearlet系数，而随机噪声并无规律而言，通常对应着较小的Shearlet系数。对含噪声的地震数据P(t)进行Shearlet稀疏变换，获得对应的Shearlet系数，如下：

Cp(j,l,k)=SH{P(t)}=〈P(t),φj,l,k〉。

(8)

利用稀疏变换后的Shearlet系数估计随机噪声的方差，再通过阈值估计算子获得对应的阈值参数：

(9)

(10)

式中：N表示地震数据的总采样点数；median(·)表示对数据矩阵中所有元素取中值。

对稀疏变换后的Shearlet系数进行阈值处理

(11)

式中：CTh(j,l,k)表示阈值处理后的Shearlet系数矩阵，保留较大的Shearlet系数，去除包含随机噪声较小的Shearlet系数，就实现了Shearlet变换阈值法去除随机噪声的目的[21]。

式(10)利用随机噪声的残差及尺度参量对阈值大小进行估计，所以随机噪声的阈值仅随尺度的不同而发生改变，实现了阈值随尺度变化的自适应。

3 Shearlet域尺度角度自适应阈值

对一单炮地震数据进行Shearlet变换，分别求取不同尺度角度下地震数据的信噪比和L2范数，如图1和2所示，分别展示了第1、2、3、4尺度下各个角度的信噪比和L2范数柱状图。通过对柱状图分析可知，在不同的尺度角度域中，数据的信噪比、L2范数相差很大。

单一的利用尺度参量进行随机噪声阈值估计不能最大程度地利用Shearlet的稀疏特性，同时也不能达到最优的去噪效果。地震数据同相轴具有方向性，而随机噪声比较均一，在Shearlet变换的某些方向上，地震数据的分布将会很少，这时Shearlet系数绝大部分表示随机噪声，用此部分Shearlet系数就能估计出某一尺度下不同角度随机噪声的能值范围，即将有效信号分布最小的尺度和方向上分布的噪声作为该尺度的噪声，并将其剔除，得到分布在各尺度和方向上的有效地震数据。随尺度和方向自适应变化的阈值函数如下：

图2 Shearlet域不同尺度角度下地震数据的L2范数Fig.2 L2 norm of seismic data at different scale angles in Shearlet domain

(12)

式中：ej,l为尺度j、方向l上的稀疏变换系数矩阵的L2范数，j=1,2,…,J，J为分解总尺度；σ为噪声标准差； min(·)为最小值函数；λj与信噪比相关，

(13)

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

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

自适应阈值函数将信噪比、随机噪声残差、Shearlet系数L2范数作为阈值设定的因素，即不同尺度、角度的阈值估计权值系数不同，可以自适应求取相应的阈值，以达到最优的随机噪声压制效果，避免有效信号的损伤。

4 理论模型试算

为验证基于Shearlet稀疏变换自适应阈值法对随机噪声压制的有效性，截取部分Sigsbee速度模型(图3a)采用有线差分进行正演模拟，采用道距20 m、炮距20 m，采样间隔1 ms，8 s接收，加入随机噪声的单炮记录如图3b所示。随机噪声对深层弱反射信号干扰很大，部分微弱的反射信号在随机噪声的背景下已经很难用肉眼识别。采用传统Shearlet域阈值法和Shearlet域尺度角度自适应阈值法对此单炮记录进行随机噪声压制，结果如图4。对比图3和图4可以看出，基于Shearlet变换的传统阈值方法和本文提出的尺度角度自适应阈值方法都能够起到压制单炮记录中随机噪声的作用，主要的反射同相轴成像效果得到了提升。但是，图4a的局部弱反射同相轴并没有得到提升，有些本该存在的同相轴受到了损伤，相应的图4c中随机噪声波场含有有效反射波同相轴信息。图4b的随机噪声去除效果明显优于图4a，并且通过观察去除的随机噪声波场(图4d)未发现有效反射同相轴信息，压制随机噪声的同时，有效保护了中深层的弱反射信息。图4b中构造边缘的弱反射同相轴成像效果得到了提升，单炮记录上对构造细节的呈现优于图4a，表明本文提出的Shearlet域尺度角度自适应阈值方法在随机噪声压制和弱信号保护上是有效的。图4e、f分别为传统阈值方法去噪的第三、第四尺度结果，图4g、h分别为本文方法去噪的第三、第四尺度结果，在单独的尺度域中，本文提出的Shearlet域尺度角度自适应阈值方法在随机噪声压制和弱信号保护上也是优于传统Shearlet域阈值方法的。

图3 Sigsbee模型(a)及包含随机噪声的单炮记录(b)Fig.3 Sigsbee model(a) and single shot record containing random noise(b)

a—传统Shearlet域阈值方法;b—Shearlet域尺度角度自适应阈值方法;c—传统Shearlet域阈值方法去除噪声;d—Shearlet域尺度角度自适应阈值方法去除的噪声;e—传统Shearlet域阈值方法第三尺度结果;f—传统Shearlet域阈值方法第四尺度结果;g—Shearlet域尺度角度自适应阈值方法第三尺度结果;h—Shearlet域尺度角度自适应阈值方法第四尺度结果a—traditional Shearlet domain threshold method;b—Shearlet domain scale angle adaptive threshold method;c—random noise by traditional Shearlet domain threshold method;d—random noise by Shearlet domain scale angle adaptive threshold method;e—the 3rd scale result of traditional Shearlet domain threshold method;f—the 4th scale result of traditional Shearlet domain threshold method;g—the 3rd scale result of Shearlet domain scale angle adaptive threshold method;h—the 4th scale result of Shearlet domain scale angle adaptive threshold method图4 单炮记录随机噪声去除效果Fig.4 Shearlet domain threshold method for the removal of random noise

为验证本文提出的方法对叠后数据的适用性，抽取含随机噪声模拟数据的零偏移距数据，如图5，理论上零偏移距数据相当于一次叠加的叠后剖面。分别采用传统Shearlet域阈值法和本文提出的Shearlet域尺度角度自适应阈值法对此零偏移距数据进行随机噪声压制，结果如图6所示。对比图6a和图6b，两种方法对随机噪声都起到了很好的去除效果，但是，在交互薄层、构造边缘及深部构造上，图6b保留的有效信息更全，更能适应高精度高保真度的地震数据处理要求。通过在零偏移距地震数据上测试可知，Shearlet稀疏变换自适应阈值法的随机噪声压制方法同样适用于叠后地震数据，能够取得理想的处理效果。

图5 包含随机噪声的零偏移距数据Fig.5 Zero-offset record containing random noise

a—传统Shearlet域阈值方法;b—Shearlet域尺度角度自适应阈值方法;c—传统Shearlet域阈值方法去除的噪声;d—Shearlet域尺度角度自适应阈值方法去除的噪声a—traditional Shearlet domain threshold method;b—Shearlet domain scale angle adaptive threshold method;c—Random noise by traditional Shearlet domain threshold method;d—Random noise by Shearlet domain scale angle adaptive threshold method图6 零偏移距数据随机噪声去除效果Fig.6 Shearlet domain threshold method for the removal of random noise

5 实际数据处理

为探明松辽盆地基底构造结构及盆山构造演化过程，本单位近几年在松辽盆地部署了多条深反射地震测线。图7a为松辽盆地西部的一条深反射叠加剖面，观测系统为道距25 m、炮距100 m，960道接收，采样间隔1 ms，25 s记录。经过常规地震数据处理流程处理后(包含静校正、异常振幅压制、面波压制、先行干扰波压制、补偿、反褶积等)，剖面能够呈现出主要的地质构造走势。但是受随机噪声的影响，剖面整体信噪比并不高，中浅层构造细节刻画不清晰，中深层莫霍面上下反射同相轴不连续，给后期的解释工作带来一定干扰。

将本文提出的Shearlet域尺度角度自适应阈值随机噪声压制方法应用到松辽盆地深反射地震数据中。图7b为经过Shearlet域自适应阈值随机噪声压制后的叠加剖面，对比图7a和图7b可以发现，经过随机噪声压制后，剖面的成像质量明显提升，中深层的构造细节和莫霍面的走势都得到清晰呈现。图7c为去除的随机噪声，可以发现本文方法在实际数据应用中只是去除了背景场的随机噪声部分，未对中深层的弱有效反射信号带来损伤。图8a、c对应图7a中的黄色和蓝色方框圈定的位置，图8b、d对应图7b中的黄色和蓝色方框圈定的位置，通过局部放大进行对比可以看出，经过Shearlet域尺度角度自适应阈值随机噪声压制方法处理后，微弱的地震同相轴和构造边缘位置都得到了清晰的成像，剖面的精细成像质量得到了很大的提升。此外，通过图9去噪前后频谱对比分析发现，去噪后10～65 Hz主频带能量没有损耗，证明了该方法对地震剖面实现了保幅处理。

a—松辽盆地深反射地震数据叠加剖面;b—Shearlet域尺度角度自适应阈值去噪后剖面;c—去除的随机噪声a—superimposed profile of deep reflection seismic data in Songliao Basin;b—Shearlet domain scale angle adaptive threshold denoised profiles;c—random noise图7 松辽盆地地震数据随机噪声去除效果Fig.7 Random noise removal effect of field seismic data in Songliao Basin

a—图7a中黄色方框放大;b—图7b中黄色方框放大;c—图7a中蓝色方框放大;d—图7b中蓝色方框放大a—enlarge the yellow box area in figure 7a;b—enlarge the yellow box area in figure 7b;c—enlarge the blue box area in figure 7a;d—enlarge the blue box area in figure 7b图8 图7局部放大对比Fig.8 Zoom in of Fig.7