地震随机噪声压缩感知迭代压制方法

刘璐， 刘洋， 刘财， 郑植升

吉林大学地球探测科学与技术学院， 长春 130026

0 引言

随机噪声是地震勘探无法避免的噪声之一，其频率、视速度、方向、振幅均不固定，往往广泛分布于有效信号频带内，严重影响着地震记录的信噪比以及后续的处理和解释工作.随着地震勘探的目标逐渐转向“双复杂”(复杂地表条件和复杂介质)以及深层条件下的资源探查，随机噪声对于复杂信号和弱信号有效识别的影响更加突出，因此，如何有效压制随机噪声，提升地震记录信噪比，恢复更高质量的地震有效信号振幅信息尤为重要.随机噪声主要来源于环境噪声、风动、记录仪器以及检波器与地面耦合差等情况，通常被认为是平稳、高斯随机过程(李庆忠，1993)，而近年来的研究发现地震勘探中的随机噪声并不是传统意义上的平稳随机过程，其平稳性随记录时长增加而降低，且与采集环境也有关系，尤其在山区等地表条件复杂的地区常采集到非平稳随机噪声记录(钟铁等，2017).因此，开发能够保护复杂地震有效信号的非平稳随机噪声压制方法具有重要的现实意义.

目前压制地震随机噪声的方法有多种类型，包括：(1)数据域滤波类的噪声压制方法；(2)f-x域及t-x域预测滤波类方法；(3)基于深度学习的地震数据去噪方法；(4)变换类噪声压制方法；在数据域滤波的地震数据去噪方法中，中值滤波与均值滤波是最常见的方法，Liu等(2006)在传统中值滤波方法的基础上发展了二维多级中值滤波方法，相对于一维中值滤波，该方法能够有效保护信号细节；Bonar 和Sacchi(2012)提出非局部均值滤波方法，利用数据点间的结构相似性以及随机噪声的非相关性进行去噪，具有良好的结构和细节保持能力.预测滤波方法中，Amba和Claerbout(1997)介绍并比较了f-x域及t-x域两种同相轴预测方法，Liu和Li(2018)及国胧予等(2020)在传统f-x域及t-x域预测滤波中引入流式处理框架来解决地震数据随机噪声的压制问题，实现了计算效率的提升.基于深度学习的去噪方法在近年来发展迅速，并在应用中形成了多种网络结构，地震数据去噪中常见的有U-net、Res-net等基于神经网络的去噪方法(Liu et al.，2018；李海山等，2020).基于稀疏变换的地震数据去噪方法主要利用地震数据的稀疏特性进行信噪分离，自20世纪80年代，地震数据的稀疏特性就被应用于解决地震随机噪声的压制问题，小波变换(Mallat and Zhong，1992)克服了傅里叶变换(Canales，1984；Vassiliou and Garossino，1998)在非平稳信号应用中的缺陷，实现了多频率的分解，小波变换至今在地震数据随机噪声的压制中仍是一种应用广泛的变换方法(高静怀等，2006；吴招财和刘天佑，2008)，近年来，针对传统小波变换难以表征高维数据方向性等问题，发展出了多尺度几何分析方法，如ridgelet(Donoho，2000)变换、curvelet变换(Herrmann et al.，2008)、seislet变换、shearlet变换等都在实际应用中取得了良好的效果，其中，seislet变换及shearlet变换以其良好的稀疏性在地震数据的去噪(刘洋等，2009；邓盾等，2019；童思友等，2019)中得到了较多应用，以地震数据稀疏特性为基础的去噪方法是近年来的研究方向之一，其能够提供比传统方法更加有效的弱信号恢复能力.随着数学方法以及信号处理等相关领域的不断发展，地震随机噪声压制方法也在不断地更新，目标是在时空变地震信号保真的前提下提高信噪比，为解决实际应用中的复杂问题提供更多的方法选择.

压缩感知理论由Candès(2006)提出，并在图像处理领域得到了广泛的应用.压缩感知理论能够突破奈奎斯特(Nyquist)采样定理的限制，即使采样频率没有达到信号最高频率的二倍，也可以将原始信号精确地恢复出来，因此，基于压缩感知理论的方法在地震数据的插值、去噪等方面能够取得良好的应用效果.杨冠雨等(2020)在压缩感知中结合shearlet变换，并使用双正则化项进行约束，实现了地震数据的插值重建，南方舟等(2018)将压缩感知与curvelet变换相结合，实现了海底地震仪采集数据的去噪.压缩感知理论下随机噪声的压制方法主要包括两个方面：一是稀疏变换基的选择，二是重构方法的选择.由于地震数据具有特殊性，常规图像处理领域的变换方法往往无法为地震数据提供理想的稀疏表征，因此Fomel和Liu(2010)提出了有效压缩地震数据同相轴的seislet变换，seislet变换属于类小波变换，是在小波提升算法的基础上结合地震数据局部倾角等属性构建的新方法，能够针对地震数据特性实现更优的稀疏表征.刘洋等(2009)将低阶seislet变换扩展到高阶，Liu和Fomel(2010)将seislet变换与波动方程炮检距连续算子结合构造了OC-seislet变换，Liu等(2015,2017)利用动校正方程构造了依赖速度的VD-seislet变换以及广义VD-seislet变换，试图降低复杂地质环境及随机噪声对地震数据属性表征的影响，增强seislet变换对地震数据的稀疏表征能力.有效的稀疏表征是压缩感知去噪效果的前提和保证，由于有效信号具有可压缩性，而随机噪声不具有可压缩性，经稀疏表征后的有效信号和随机噪声可以在重构过程中分离.常见的重构算法有贪婪类方法和凸优化方法.贪婪类方法的基本思想是将全局优化问题转化为分阶段优化，在每一次迭代中只考虑当前阶段的最优解，算法时空复杂度低，但不能保证得到全局最优解，贪婪方法以匹配追踪类(Mallat and Zhang，1993；陈兴飞和孙红梅，2019)为代表.凸优化方法主要包括内点法(Kim et al.，2007)、梯度投影法(Figueiredo et al.，2007)、全变差方法(Chen et al.，2001)、Bregman迭代(Yin et al.，2008)、迭代硬阈值和迭代软阈值(Pope，2008；Daubechies et al.，2004)等.在凸优化方法中，迭代收缩阈值方法由于其计算复杂度低，适用于高维大规模数据的性质受到了广泛的应用.尽管压缩感知理论在地震数据随机噪声压制方面有一些应用，但是其去噪数学反问题和匹配的实现方法并不明确.

本文针对地震数据随机噪声压制问题，建立完善的压缩感知框架下稀疏信号重建数学模型，通过求解基追踪降噪模型下的去噪反演目标函数，提出基于迭代软阈值算法的“采集-重建-修复”方案，对含随机噪声数据中的有效信号进行迭代分离.通过选取和对比不同的稀疏变换基，论证有效信号的稀疏表征能力对于提高地震随机噪声压缩感知迭代压制方法有效性的重要意义，同时，针对稀疏优化问题中正则化项使有效信号衰减的现象，利用除偏算法对损失的振幅进行恢复，进一步提高本文方法的弱信号保护能力.将本文方法与工业标准随机噪声压制FXDECON方法进行比较，模型及实际数据测试结果表明，压缩感知迭代去噪方法可以有效地压制地震勘探记录中的较强振幅随机噪声，有效保护复杂构造的同相轴信息.

1 基本理论

1.1 压缩感知基本理论回顾

假设一组有限长的离散信号x∈RN,在某个变换域内可以线性表示为x=Ψc，式中c是x在Ψ域的变换域系数，Ψ是N×N的稀疏变换基，也用于表示稀疏反变换.当信号在某个变换基Ψ下的系数矩阵中绝大部分元素为零或约等于零时，就称该信号是变换域内的稀疏信号.信号x本身或者变换域系数c具有稀疏性是压缩感知理论应用的基础.

当信号本身是稀疏的，压缩感知的一般采样问题可以描述为(Candès and Wakin，2008)：

y=Φx，

(1)

式中，x∈RN为原始信号，观测矩阵Φ数据维度为M×N，当M≪N时，可以获得关于x的欠采样信号y∈RM，压缩感知理论主要用于实现利用y重建x的过程.

在地震数据中，有效信号本身一般不具有稀疏性，但是可以经过稀疏变换进行稀疏表征，即x=Ψc，此时压缩感知问题为

y=ΦΨc=Θc，

(2)

式中，Θ=ΦΨ称为传感矩阵.

此时，由观测信号y恢复原始信号x的最优化问题可以描述为

min{‖c‖0:y=Φx}，

(3)

信号恢复的过程利用压缩感知中的重构算法.通过压缩感知技术来观测数据，可以节省大量采样空间，且能够准确重建不满足采样定理条件的数据.压缩感知理论不仅可以有效解决稀疏采样的恢复问题，还能够用于压制数据中的随机噪声干扰.

1.2 压缩感知框架下的地震数据随机噪声压制问题

信号中含有随机噪声干扰时，公式(1)变为

y=Φd=Φ(x+z)=Φx+z′，

(4)

式中，d为含随机噪声地震数据，x为不含噪声的纯信号，z为随机噪声，z′为误差项，即压缩感知采样后的噪声，y为压缩感知的观测信号.研究如何从观测信号y中恢复有效信号x成为压缩感知框架下的地震数据随机噪声压制问题(Candès and Wakin，2008).

压缩感知最优化问题可以描述为公式(3)所示的L0范数最小化的最优化问题，但该问题是一种典型的非凸NP-hard问题，在数值计算上难以有效求解，一种可行的方法是用L1范数替代问题中的L0范数.当信号满足在变换域内稀疏，稀疏变换基Ψ和观测矩阵Φ满足非相干性或等距约束性时，最优化问题可以由L0范数最小化转化为L1范数最小化，在基追踪(Basis Pursuit，BP)问题(Chen et al.，2001)下进行求解：

min{‖c‖1∶y=Φx}.

(5)

基追踪准则下的最优解意味着所得解的L1范数最小.在公式(4)数据中含有噪声的条件下，求解压缩感知框架下的随机噪声压制问题需要对基追踪问题(5)的约束条件进行松弛近似，得到基追踪降噪(Basis Pursuit De-Noising，BPDN)问题(Chen et al.，2001)，该问题可以描述为

(6)

以及相应的非约束问题(Figueiredo et al.，2007)：

(7)

迭代软阈值方法由非迭代的软阈值方法发展而来，Donoho等(1995)提出软阈值去噪方法，该方法通过一个简单的阈值步骤将含噪声数据中的有效信号与噪声分离，且估计的去噪结果具有较好的平滑性，但是软阈值去噪的目标函数与基追踪降噪的更一般化目标函数并不相同，因此软阈值方法不能解决基追踪降噪问题.为了适应多种稀疏变换下从含噪声数据中估计有效信号的线性反问题，Daubechies等(2004)提出1≤p≤2时解决Lp问题的迭代软阈值方法，并且给出任意初值时的收敛性证明，为迭代软阈值方法的应用提供了理论基础.

迭代软阈值将软阈值函数与Landweber迭代(Landweber，1951)相结合，在迭代的每一步中应用软阈值函数.迭代软阈值方法可以处理多种Lp范数(1≤p≤2)问题，当取L1范数时，迭代软阈值方法的目标函数与基追踪降噪问题(7)一致，因此，迭代软阈值方法可以用于解决压缩感知框架下的随机噪声压制问题.利用迭代软阈值方法求解压缩感知问题的迭代解如公式(8)所示，迭代初始值任意的情况下公式(8)最终收敛于公式(7)的解：

xn+1=Sλ[xn+ΦT(y-Φxn)]，

(8)

其中：

Sλ[·]=ΨsλΨ-1[·]，

(9)

[·]T表示为矩阵的转置，Ψ-1为稀疏正变换，sλ是对信号中的每一个元素执行软阈值，软阈值函数定义为

(10)

压缩感知的信号重构效果不仅与重构算法的选择有关，还取决于稀疏变换及观测矩阵的选择.选择的稀疏变换应该使有效信号经变换后的变换域系数足够稀疏，实现有效信号和随机噪声的变换域系数振幅差异最大化，保证重构信号的精度.离散小波变换在科学和工程上有着广泛的应用，其主要基于分片平滑假设，在地震数据随机噪声压制方面有着广泛的应用.但传统小波变换在方向性上缺乏灵活性，同时地震数据往往不满足其应用的分片平滑假设，因此，利用图像中方向特性的类小波变换得到了持续的发展，seislet变换(Fomel and Liu，2010)可以根据地震数据的不同特征对地震数据进行更加有效地压缩，因此本文选取离散小波变换和seislet变换作为压缩感知迭代噪声压制方法中的稀疏变换基，并对比处理效果.

观测矩阵的设计是另一个影响信号重建精度的重要环节.压缩感知理论本质上是一个欠定问题，为使观测信号y中包含有足够多原始信号的信息以便进行信号的重构，观测矩阵与稀疏变换基需要满足非相干性，相干性越小，欠采样的观测信号y中包含的有效信息越多，才能够保证信号重建的精确性.在观测矩阵元素的选取方面，Candès和Wakin(2008)给出了随机生成的方案.在地震数据的插值问题中，观测矩阵与观测系统的设计有关，单位矩阵是在全采样的情况下对观测系统的数学表示，欠采样时，插值问题的观测矩阵为对角线随机缺失的单位阵.在去噪问题中，观测矩阵与观测系统无关，主要考虑计算上的优势.高斯随机矩阵作为观测矩阵与正交稀疏变换基满足非相干性，提供了信号恢复的前提，结构简单易于构造且具有普适性，是压缩感知方法中常用的观测矩阵，因此本文选取高斯随机矩阵作为压缩感知迭代噪声压制方法中的观测矩阵，且有Φ～N(0,1/M)，满足均值为0，方差为1/M的高斯分布，其中M为Φ的行维度.

1.3 L1正则化振幅保真的除偏方法

最优化问题(6)和(7)中由于L1正则化项的存在，重构产生的结果往往在幅度上相比于原信号通常会有一定的损失，可以选择除偏(debiasing)算法(Wright et al.，2008)来恢复L1正则项引起的信号幅度损失.

(11)

本文研究表明，当选择的稀疏变换能够提供高效压缩结果时，L1正则化项的影响将变得较小，有效信号的损失也较小，此时除偏步骤可以被省略，因为除偏方法在恢复有效信号的同时会恢复出一定的噪声，降低处理结果的信噪比.除偏步骤主要应用于对信号的保真性要求更高的处理任务中.

1.4 压缩感知迭代噪声压制方法的实现流程

在压缩感知框架中使用高斯随机矩阵与稀疏变换的组合作为压缩感知的传感矩阵，对含噪声信号进行“采集-重建-修复”的压缩感知迭代噪声压制过程可以归纳如下3个步骤：

(1)带有随机噪声的实际数据为原始信号d，根据公式(4)，利用高斯随机矩阵Φ对信号d进行采样得到观测信号y=Φd，实现“采集”数据环节.

具体流程框架如图1所示.本文方法的计算效率不仅取决于Landweber迭代次数，还取决于稀疏变换的计算速度，为了保证噪声压制效果，当选用seislet变换作为压缩感知迭代噪声压制方法中的稀疏变换基时，计算效率比传统的迭代去噪方法更低一些，但是可以通过并行计算改善处理效率.

图1 实现流程框架Fig.1 Processing framework

2 模型测试

首先选取复杂的Sigmoid叠后地震数据模型(Claerbout，2008)测试本文方法，该模型中包括了倾斜地层、断层和不整合面等构造特征，时间方向的采样间隔为4 ms，空间方向的采样间隔为8 m.为了定量分析噪声压制的结果，本文选取全局信噪比作为去噪效果的衡量(刘洋等，2017)：

(12)

图2a是不含噪声的模型数据，图2b是加入较强振幅的非平稳随机噪声数据，含噪声数据信噪比为4.4 dB，与噪声方差为7.5×10-7的平稳噪声信噪比相当.模型数据中的非平稳噪声被定义为时间和空间方向上性质不稳定的高斯白噪.将时间方向上平稳的高斯随机噪声序列记为n1(t)(图3a)，对平稳随机噪声序列中的部分取值进行放大或改变取值的分布特征可以得到非平稳高斯随机噪声(钟铁，2016)，则通过平稳的噪声序列n1(t)可以构造并表示出非平稳的噪声序列(图3b)，非平稳噪声添加在时间方向上，记为n2(t)：

(13)

其中k为噪声序列中取值放大的开始位置，取值放大的长度为50，放大倍数取为2倍，不同地震道上噪声放大的起始时间k随机.将全频带的非平稳噪声进行带通滤波，仅保留合理频率范围内的噪声，如图3b所示，构造出的随机噪声具有非平稳的时变特征，更接近于真实的随机噪声.

图2 模型数据(a) Sigmoid模型； (b) 含噪声数据.Fig.2 Synthetic model data(a) Sigmoid model； (b) Noisy data.

图3 不同平稳性的随机噪声(a) 平稳随机噪声； (b) 非平稳随机噪声.Fig.3 Random noise with different stationarities(a) Stationary random noise； (b) Non-stationary noise.

图4 不同去噪方法结果对比(a) f-x预测滤波去噪结果； (b) f-x预测滤波压制的噪声； (c) 本文迭代方法小波变换下去噪结果(未除偏)； (d) 本文迭代方法小波变换下压制的噪声(未除偏)； (e) 本文迭代方法小波变换下去噪结果(除偏)； (f) 本文迭代方法小波变换下压制的噪声(除偏).Fig.4 Comparison of denoised results by different methods(a) Result using f-x prediction filtering； (b) The noise removed by f-x prediction filtering； (c) Result using the proposed method with wavelet transform(before debiasing)； (d) The noise removed by the proposed method with wavelet transform(before debiasing)； (e) Result using the proposed method with wavelet transform(after debiasing)； (f) The noise removed by the proposed method with wavelet transform (after debiasing).

为了进一步提升本文方法的有效性，使用seislet变换作为稀疏变换基，图5a为seislet稀疏变换在本文迭代方法框架下的去噪结果，为了和小波变换进行对比，选取百分比阈值为16%，迭代次数为10，从去噪结果中可以看出本文方法的能够有效压制大部分随机噪声，尤其有效弱信号损失的恢复能力增强，处理之后的同相轴连续性提高，信噪比提升至11.7 dB(图5a).在差剖面(图5b)中可以看出，本文方法压制的主要为随机噪声和少部分的断层信息.由于seislet变换对有效信号的压缩能力更强，因此L1正则化项的影响很小，故“修复”步骤被省略，此时除偏修复的幅度主要来自噪声，会降低处理效果，因此“修复”步骤主要适用于压缩能力一般的稀疏变换基.与f-x域流式预测滤波方法(图5c)相比，本文方法计算结果信噪比更高，f-x域流式预测滤波方法处理结果信噪比为8.5 dB，差剖面(图5d)中可以看出本文方法去掉的噪声更多且损失有效信号更少，但流式方法的计算量要低于本文方法；f-x域RNA方法(图5e)处理结果信噪比为10.7 dB，对比可见本文方法的信噪比更高，本文方法差剖面(图5f)中去掉的有效信息更少，而且没有频率域预测滤波方法中常见的虚假同相轴问题.

图5 本文迭代方法处理结果(a) seislet变换下去噪结果； (b) seislet变换下压制的噪声; (c) f-x域流式预测滤波方法去噪结果； (d) f-x域流式预测滤波方法压制的噪声； (e) f-x域RNA方法去噪结果； (f) f-x域RNA方法压制的噪声.Fig.5 Denoising results using the proposed method(a) The denoised result using the proposed method with seislet transform ； (b) The noise removed by the proposed method with seislet transform； (c) The denoised result using f-x streaming prediction filter； (d) The noise removed by f-x streaming prediction filter； (e) The denoised result using f-x RNA； (f) The noise removed by f-x RNA.

如图6所示为Sigmoid模型不同方法处理结果的f-k谱.图6a为不含噪声模型的f-k谱，图6b为含噪声模型的f-k谱，对比各方法f-k谱可见，基于seislet变换的本文方法(图6f)去噪效果最佳，基于wavelet变换的本文方法(图6d)去除随机噪声不够彻底，除偏(图6e)还会引入少量噪声，而f-x预测滤波方法(图6c)处理后仍然存在部分高波数的随机噪声干扰.f-x域流式预测滤波方法(图6g)处理后在有效信号的频带范围内残留噪声较多，f-x域RNA方法(图6h)去除噪声较为彻底，仅在高波数部分残留少量噪声.各方法间去噪性能对比如表1所示.基于seislet变换的本文方法去噪结果信噪比(SNR)最高，均方误差(MSE)最小，虽然迭代次数稍多于基于wavelet变换的本文方法，但其去噪效果好于后者，f-x预测滤波方法信噪比最低，均方误差最大.一些其他的稀疏变换基，如shearlet(童思友等，2019)、contourlet(Li and Gao，2013)也可以用于本文去噪方法框架中，理论上也会实现较好的效果.

图6 模型数据f-k谱(a) 原始信号f-k谱； (b) 含噪声数据f-k谱； (c) f-x预测滤波结果f-k谱； (d) 本文迭代方法小波变换下未除偏结果f-k谱； (e) 本文迭代方法小波变换下除偏结果f-k谱； (f) 本文迭代方法seislet变换下去噪结果f-k谱； (g) f-x域流式预测滤波方法去噪结果f-k谱； (h) f-x域RNA方法去噪结果f-k谱.Fig.6 Comparison of different f-k spectras(a) Original signal； (b) Noisy data； (c) Denoised result using f-x prediction filtering； (d) Denoised result using the proposed method with wavelet before debiasing； (e) Denoised result using the proposed method with wavelet after debiasing； (f) Denoised result using the proposed method with seislet； (g) Denoised result using f-x streaming prediction filter； (h) Denoised result using f-x RNA method.

表1 各方法去噪性能对比Table 1 Comparison of denoising performance in different methods

3 实际数据处理

3.1 实际叠后数据处理

接下来，利用某地区实际叠后时间偏移数据对本文方法进行测试(图7a)，该数据在时间方向上为700个采样点，采样间隔为1 ms，空间方向310个地震道.该数据浅层存在近平行的同相轴，在深层存在比较复杂的具有不同倾角的倾斜同相轴，数据中的较强随机噪声严重影响着复杂同相轴的连续性，进而影响高精度的地层解释.图7b为f-x预测滤波方法的去噪结果，噪声水平得到一定程度的衰减，同时浅层同相轴连续性增强.将滤波后的结果与原始数据相减得到滤除的噪声剖面，如图7c所示.在噪声剖面中可以看到标准f-x预测滤波方法对有效波振幅的保护效果差，即使浅层比较平稳的近水平同相轴也被部分衰减.图7d为使用seislet变换作为稀疏变换的压缩感知迭代噪声压制方法的去噪结果，迭代次数为10，百分比阈值选取为18%.从图7b和图7d对比可以看到，本文方法的处理结果优于f-x预测滤波方法，处理后数据同相轴清晰，深层复杂构造得到更好的保护，而f-x预测滤波的处理结果中同相轴被一定程度地模糊化，有效信号损失较大.对比两种方法的差剖面(图7c和图7e)，本文方法去除的主要是随机噪声，很难直观看到有效波的损失，而f-x预测滤波损失了大量有效波同相轴振幅.最后，通过处理结果的f-k谱来进一步对比两种方法的处理效果，如图8所示.图8c表明本文能够更好地压制原始数据低频和高波数范围内(图8a)的随机噪声干扰，而f-x预测滤波方法处理后仍然存在部分高波数的随机噪声干扰(图8b).

图7 实际数据处理结果(a) 实际数据； (b) f-x预测滤波方法的去噪结果； (c) f-x预测滤波方法压制的噪声； (d) 本文方法的去噪结果； (e) 本文方法压制的噪声.Fig.7 Processing result of real data(a) Real data; (b) Denoised result using f-x prediction filtering; (c) Noise removed by f-x prediction filtering； (d) Denoised result using the proposed method with seislet transform； (e) Noise removed by the proposed method with seislet transform.

图8 处理结果f-k谱对比(a) 叠后实际数据f-k谱； (b) f-x预测滤波方法的去噪结果f-k谱； (c) 本文方法的去噪结果f-k谱.Fig.8 Comparison of f-k spectra for processing results(a) Poststack real data； (b) Denoised result using f-x prediction filtering； (c) Denoised result using the proposed method.

3.2 实际叠前数据处理

利用某地区陆上实际叠前数据对本文方法进行进一步的测试，图9a为去面波后的CMP道集，可以看到该数据信噪比较低，强随机噪声的存在严重影响着该数据中的反射同相轴的连续性，有效信号的能量较弱.本文方法可以直接扩展为三维方式，三维处理方法能够进一步提高效果，但是将会增加计算负担，所以将三维数据按照二维切片进行独立处理.图9b为f-x预测滤波方法的去噪结果，可见由于强噪声影响，滤波后的结果中也几乎不可见同相轴的存在，数据信噪比仍然较低.图9c为使用VD-seislet变换作为稀疏变换的压缩感知迭代噪声压制方法的去噪结果，其中seislet变换的倾角模式由均方根速度计算(Liu et al.，2015)，迭代次数为11，百分比阈值选取为10%.从图9c中可以看到，经本文方法处理后，结果中同相轴连续性明显增强，同相轴形态恢复良好，数据信噪比得到提升.对比两种方法处理后的叠加剖面(图10a和图10b)，本文方法的叠加剖面中同相轴连续性明显增强，尤其中浅层结构恢复较为明显，f-x预测滤波方法的叠加剖面中同相轴能量较弱，连续性较差，随机噪声的存在仍严重影响着该剖面的信噪比.

图9 叠前实际数据处理结果(a) 叠前实际数据； (b) f-x预测滤波方法的去噪结果； (c) 本文方法的去噪结果.Fig.9 Processing results of pre-stack real data(a) Pre-stack real data; (b) Denoised result using f-x prediction filtering; (c) Denoised result using the proposed method with seislet transform.

图10 不同方法处理结果叠加剖面(a) f-x预测滤波方法去噪结果叠加剖面； (b) 本文方法的去噪结果叠加剖面.Fig.10 Stacked profiles using different methods(a) f-x prediction filtering； (b) The proposed method with seislet transform.

4 结论与展望

本文系统地论述了压缩感知框架下的随机噪声压制问题，建立了压缩感知、软阈值去噪和基追踪降噪的理论关联，提出了基于迭代软阈值算法的“采集-重建-修复”技术方案，为提升本文方法的时空变信号保护及非平稳随机噪声压制能力，选用对地震信号具有更高压缩能力的seislet变换作为稀疏变换，同时利用匹配的除偏方法，对最优化问题中由于较差压缩能力的稀疏变换在L1正则化条件中所引起的振幅衰减实现有效的恢复.通过与工业标准的f-x预测滤波方法进行对比，本文方法能够更好地处理同相轴倾角的时空变化，避免频率域预测滤波类方法的虚假同相轴问题，能够压制强随机噪声且更好地平衡时空变有效信号保护和随机噪声压制的关系，同时能够适应随机噪声的非平稳分布特征，提供一种有效的随机噪声压制方法.理论模型和实际数据测试结果表明本文提出方法能够在保护时空变有效信息的前提下压制较强幅度的非平稳随机噪声，提高同相轴连续性，为后续高精度的解释提供数据基础.为提升本文方法的性能，仍需要进一步研究，如寻求对数据更优的稀疏表示方法，例如Shearlet、Contourlet及超完备冗余字典可能会提供更加优秀的数据表征.

致谢感谢同济大学王本锋研究员的建议和讨论，感谢两位匿名审稿专家提出的宝贵意见.