地球物理学报  2021, Vol. 64 Issue (12): 4629-4643   PDF    
地震随机噪声压缩感知迭代压制方法
刘璐, 刘洋, 刘财, 郑植升     
吉林大学地球探测科学与技术学院, 长春 130026
摘要:复杂地表和复杂介质条件下,随机噪声往往严重影响着复杂地震信号的信噪比,同时深层地球物理目标探查中弱地震信号总是被随机噪声所掩盖,如何有效地压制随机噪声干扰、恢复有效地震信号仍然是高精度地震勘探中的关键问题.压缩感知理论突破了奈奎斯特采样定理的限制,利用有效地震信号的可压缩性和稀疏性,提供了从不可压缩随机噪声中进行有效信号分离的数据原理.本文系统分析压缩感知框架下地震随机噪声压制的稀疏优化反问题,提出了基于迭代软阈值算法的"采集-重建-修复"方案对该问题进行求解.在实现高度稀疏表征的基础上进行地震数据的压缩感知随机观测,通过迭代反演对有效地震信号进行重构,有效提高复杂地震数据的信噪比,同时,当求解稀疏优化问题时,如果出现正则化项引起重构信号衰减现象,可以匹配除偏对衰减的有效信号进行修复.通过与工业标准f-x预测滤波方法进行比较,理论模型和实际数据处理的结果表明,压缩感知迭代噪声压制方法对复杂地震数据中的随机噪声有较好的压制效果,可以有效恢复出被较强非平稳随机噪声干扰的时空变同相轴信息.
关键词: 随机噪声压制      压缩感知      基追踪降噪      迭代软阈值      除偏     
Iterative seismic random noise suppression method based on compressive sensing
LIU Lu, LIU Yang, LIU Cai, ZHENG ZhiSheng     
College of Geo-exploration Science and Technology, Jilin University, Changchun 130026, China
Abstract: Random noise can affect the signal-to-noise ratio of complex seismic signals, especially under complex surface and subsurface conditions. Meanwhile weak signals are usually covered by random noise in geophysical exploration at depth. How to effectively suppress such random noise interference and recover the effective seismic signal remains a key problem in high-precision seismic exploration. Compressive sensing breaks though the limitation of the Nyquist sampling theory, providing an approach of effective signal separation from uncompressible random noise based on the compressibility and sparsity of effective signals. In this paper, we analyze the inverse problem of sparse optimization that corresponds to seismic random noise suppression within the compressive sensing framework. An "acquisition-reconstruction-repair" workflow based on an iterative soft threshold is proposed to resolve the inverse problem. With the highly sparse representation of seismic signals, one can obtain random observation to seismic data based on compressive sensing and rebuild only useful signals by using iterative inversion, which improves the signal-to-noise ratio of complex seismic signals. Furthermore, we use a debiasing method to reduce the attenuation of signals due to the presence of the regularization term. This debiasing step improves the proposed method by recovering the attenuated amplitude of effective signals. Compared with the industrial standard FXDECON method, the proposed method shows a better denoising result in model and field data examples and can reasonably recover nonstationary events from noisy data.
Keywords: Random noise suppression    Compressive sensing    Basis pursuit denoising    Iterative soft threshold    Debiasing    
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,1984Vassiliou 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,2008Daubechies et al., 2004)等.在凸优化方法中,迭代收缩阈值方法由于其计算复杂度低,适用于高维大规模数据的性质受到了广泛的应用.尽管压缩感知理论在地震数据随机噪声压制方面有一些应用,但是其去噪数学反问题和匹配的实现方法并不明确.

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

1 基本理论 1.1 压缩感知基本理论回顾

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

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

(1)

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

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

(2)

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

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

(3)

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

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

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

(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)下进行求解:

(5)

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

(6)

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

(7)

式中,λ为正则化参数,p≥1时Lp范数的定义为‖up=,则有L1范数为‖u1=Σi|ui|,δ为非负的实数.公式(7)中第一项中隐含的假设条件是z′在L2范数下最小,当数据d中存在随机噪声时,该条件不能达到最优化求解,此时通过调整第二项的权重可以获得合理的解,其主要原因在于随机噪声不能由稀疏变换基所表征,具有不可压缩性,因此第二项由有效信号的变换域系数所决定.公式(7)所示的基追踪降噪问题在压缩感知框架下有多种解法,称为压缩感知的重构算法.其中,迭代软阈值方法在保持信号连续性和平滑性方面具有良好的特性,是一种应用广泛的凸优化类压缩感知重构算法.

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

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

(8)

其中:

(9)

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

(10)

将该迭代过程

最终得到的估计解记为,对应的变换域系数记为,且有.Landweber迭代的收敛速度有限(Bredies and Lorenz, 2008),但其主要优势在于其迭代收敛的稳定性和结构简单,因此,本文选取该方法作为基本计算框架进行问题的求解,已有学者针对IST的计算速度进行研究,如快速阈值迭代法(Beck and Teboulle, 2009)、冷却阈值迭代法(Herrmann and Hennenfent, 2008)等,对于大规模数据快速计算时,也可以参考本文的设计方案进行改进.压缩感知迭代噪声压制方法的有效性基于:地震勘探有效信号具有可预测的特征,可以通过稀疏变换进行有效压缩,而随机噪声具有不可压缩性,因此,随机噪声经过稀疏变换后仍然是随机噪声,以较小的值随机分布于整个变换域内,而有效信号所包含的信息则集中于少数幅值较大的系数,通过设置阈值函数可以分离变换域中的信号和噪声.

压缩感知的信号重构效果不仅与重构算法的选择有关,还取决于稀疏变换及观测矩阵的选择.选择的稀疏变换应该使有效信号经变换后的变换域系数足够稀疏,实现有效信号和随机噪声的变换域系数振幅差异最大化,保证重构信号的精度.离散小波变换在科学和工程上有着广泛的应用,其主要基于分片平滑假设,在地震数据随机噪声压制方面有着广泛的应用.但传统小波变换在方向性上缺乏灵活性,同时地震数据往往不满足其应用的分片平滑假设,因此,利用图像中方向特性的类小波变换得到了持续的发展,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正则项引起的信号幅度损失.

除偏算法是针对变换域系数进行的修复,主要计算最小二乘目标函数f=(1/2)‖Ac-c′‖22,其中c′=-1d为含噪声数据变换域系数的随机观测值,A为高斯随机矩阵.与普通最小二乘问题不同的是,除偏算法将重构结果中的零元素固定为零,在剩余非零项的基础上对目标函数进行最小化求解,可以表示为

(11)

由于在迭代的过程中零元素始终固定不变,除偏过程相当于仅求解非零元素,式(11)下标集合I表示重构结果中非零元素的位置,AI表示观测矩阵A与重构结果中非零元素对应相乘的子列,为非零元素位置上的待求未知数组成的向量.共轭梯度法求解上述最小二乘问题(11)等价于求解方程组.使用共轭梯度法进行求解,最终得到除偏后的变换域系数,再对其进行反变换得到最终除偏结果.

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

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

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

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

(2) 确定软阈值函数的阈值λ,给定迭代的初始值x0和终止条件,根据公式(8)—(10)执行迭代软阈值算法,当满足终止条件时退出迭代,得到最终迭代结果及其变换域系数,其中为未进行除偏的去噪信号,实现“重建”数据环节.

(3) 将迭代最终结果的变换域系数按照除偏算法进行部分系数的更新,反变换得到除偏后的结果,实现有效信号的“修复”环节.

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

图 1 实现流程框架 Fig. 1 Processing framework
2 模型测试

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

(12)

式中,x为不含噪声的信号,x为含噪声信号或经噪声压制后的信号.

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

(13)

图 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.

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

在模型数据的测试中,为了对比压缩感知迭代噪声压制方法的效果,使用工业标准的f-x预测滤波(FXDECON)方法、f-x域流式预测滤波方法(国胧予等,2020)以及f-x域正则化非平稳自回归滤波(RNA)方法(Liu et al., 2012)处理结果进行比较.图 4a图 4b分别为f-x预测滤波方法的去噪结果和压制的噪声剖面,从图中可以看到该方法仅能去除一部分随机噪声,并且衰减了部分有效信号,这是由于本模型中有效信号并不符合平面波,平稳的f-x预测滤波方法对有效信号尤其是弯曲同相轴损失较大,同时在沿地震剖面同相轴的方向产生了虚假同相轴,该问题常见于频率域预测滤波方法,经处理后数据的信噪比提升至7.4 dB.接下来,首先使用离散小波变换作为稀疏变换基,在压缩感知理论的框架下使用迭代阈值方法进行随机噪声的压制,本文中高斯随机矩阵对原始信号的采样率取为90%,迭代软阈值中采用百分比阈值,通过在迭代中选取阈值的固定排序百分比,实现阈值的逐渐递减,取迭代误差‖xn-xn-122/‖xn-122<1×10-6时认为迭代收敛.图 4c图 4d为不使用振幅除偏修复的去噪结果及压制的噪声剖面,选取百分比阈值为16%,迭代次数为6,处理结果信噪比为8.0 dB,图中可以看到去噪效果不理想,其处理效果与平稳f-x预测滤波相近,主要是小波变换对地震数据中的有效信号无法进行有效的稀疏表征,在小波变换域中有效信号与随机噪声的振幅值区分不明显,进行阈值处理时会损失较多有效波信息而且还会保留一部分具有较大幅值的非平稳随机噪声.当使用除偏方法后,处理效果得到了一定的提升(图 4e),剖面中有效波振幅较除偏前明显增强,同相轴更清晰,在去除的噪声剖面中(图 4f),有效波信息的能量较图 4d有所减小,也说明除偏具有有效的振幅恢复能力,但是除偏也会恢复一部分随机噪声,使得除偏结果的信噪比降低至7.6 dB.

图 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个地震道.该数据浅层存在近平行的同相轴,在深层存在比较复杂的具有不同倾角的倾斜同相轴,数据中的较强随机噪声严重影响着复杂同相轴的连续性,进而影响高精度的地层解释.图 7bf-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道集,可以看到该数据信噪比较低,强随机噪声的存在严重影响着该数据中的反射同相轴的连续性,有效信号的能量较弱.本文方法可以直接扩展为三维方式,三维处理方法能够进一步提高效果,但是将会增加计算负担,所以将三维数据按照二维切片进行独立处理.图 9bf-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及超完备冗余字典可能会提供更加优秀的数据表征.

致谢  感谢同济大学王本锋研究员的建议和讨论,感谢两位匿名审稿专家提出的宝贵意见.
References
Amba R, Claerbout J. 1997. Lateral prediction for noise attenuation by t-x and f-x techniques. Geophysics, 60(6): 1887-1896.
Beck A, Teboulle M. 2009. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM Journal on Imaging Sciences, 2(1): 183-202. DOI:10.1137/080716542
Bonar D, Sacchi M. 2012. Denoising seismic data using the nonlocal means algorithm. Geophysics, 77(1): A5-A8. DOI:10.1190/geo2011-0235.1
Bredies K, Lorenz D A. 2008. Linear convergence of iterative soft-thresholding. Journal of Fourier Analysis and Applications, 14(5): 813-837.
Canales L L. 1984. Random noise reduction. //54th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 525-527.
Candès E J. 2006. Compressive sampling. //Proceedings of the International Congress of Mathematicians. Madrid, Spain: European Mathematicial Society Publishing House, 3: 1433-1452.
Candès E J, Wakin M B. 2008. An introduction to compressive sampling. IEEE Signal Processing Magazine, 25(2): 21-30. DOI:10.1109/MSP.2007.914731
Chen S S, Donoho D L, Saunders M A. 2001. Atomic decomposition by basis pursuit. SIAM Review, 43(1): 129-159. DOI:10.1137/S003614450037906X
Chen X F, Sun H M. 2019. Seismic data denoising method based on compressed sensing theory. Progress in Geophysics (in Chinese), 34(3): 1025-1031. DOI:10.6038/pg2019CC0208
Claerbout J F. 2008. Basic Earth imaging. Stanford Exploration Project, http://sepwww.stanford.edu/sep/prof/.
Daubechies I, Defrise M, De Mol C. 2004. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. Communications on Pure and Applied Mathematics, 57(11): 1413-1457. DOI:10.1002/cpa.20042
Deng D, Pei J X, Deng Y, et al. 2019. Key techniques research for improving mid-deep layer structure imaging in deep water seismic data of South China Sea. Progress in Geophysics (in Chinese), 34(2): 702-708. DOI:10.6038/pg2019BB0551
Donoho D L. 1995. De-noising by Soft-Thresholding. IEEE Transactions on Information Theory, 41(3): 613-627. DOI:10.1109/18.382009
Donoho D L. 2000. Orthonormal ridgelets and linear singularities. SIAM Journal on Mathematical Analysis, 31(5): 1062-1099. DOI:10.1137/S0036141098344403
Figueiredo M A T, Nowak R D, Wright S J. 2007. Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems. IEEE Journal of Selected Topics in Signal Processing, 1(4): 586-597. DOI:10.1109/JSTSP.2007.910281
Fomel S, Liu Y. 2010. Seislet transform and seislet frame. Geophysics, 75(3): V25-V38. DOI:10.1190/1.3380591
Gao J H, Mao J, Man W S, et al. 2006. On the denoising method of prestack seismic data in wavelet domain. Chinese Journal of Geophysics (in Chinese), 49(4): 1155-1163.
Gong X B, Han L G, Wang E L, et al. 2009. Denoising via high resolution Radon transform. Journal of Jilin University (Earth Science Edition) (in Chinese), 39(1): 152-157.
Guo L Y, Liu C, Liu Y, et al. 2020. Seismic random noise attenuation based on streaming prediction filter in the f-x domain. Chinese Journal of Geophysics (in Chinese), 63(1): 329-338. DOI:10.6038/cjg2020N0030
Herrmann F J, Hennenfent G. 2008. Non-parametric seismic data recovery with curvelet frames. Geophysical Journal of the Royal Astronomical Society, 173(1): 233-248. DOI:10.1111/j.1365-246X.2007.03698.x
Herrmann F J, Wang D L, Hennenfent G, et al. 2008. Curvelet-based seismic data processing: A multiscale and nonlinear approach. Geophysics, 73(1): A1-A5.
Kim S J, Koh K, Lustig M, et al. 2007. An interior-point method for large-scale L1 regularized least-squares. IEEE Journal of Selected Topics in Signal Processing, 1(4): 606-617. DOI:10.1109/JSTSP.2007.910971
Landweber L. 1951. An iteration formula for Fredholm integral equations of the first kind. American Journal of Mathematics, 73(3): 615-624. DOI:10.2307/2372313
Li H S, Chen D W, Wu J, et al. 2020. Pre-stack random noise suppression with deep residual network. Oil Geophysical Prospecting (in Chinese), 55(3): 493-503.
Li Q, Gao J H. 2013. Contourlet based seismic reflection data non-local noise suppression. Journal of Applied Geophysics, 95: 16-22. DOI:10.1016/j.jappgeo.2013.05.002
Li Q Z. 1993. The way to Obtain A Better Resolution in Seismic Prospecting (in Chinese). Beijing: Petroleum Industry Press: 112-114.
Liu C, Liu Y, Yang B J, et al. 2006. A 2D multistage median filter to reduce random seismic noise. Geophysics, 71(5): V105-V110. DOI:10.1190/1.2236003
Liu G C, Chen X H, Du J, et al. 2012. Random noise attenuation using f-x regularized nonstationary autoregression. Geophysics, 77(2): V61-V69. DOI:10.1190/geo2011-0117.1
Liu J, Lu W, Zhang P. 2018. Random noise attenuation using convolutional neural networks. //88th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts.
Liu Y, Fomel S. 2010. OC-seislet: seislet transform construction with differential offset continuation. Geophysics, 75(6): WB235-WB245. DOI:10.1190/1.3479554
Liu Y, Fomel S, Liu C. 2015. Signal and noise separation in prestack seismic data using velocity-dependent seislet transform. Geophysics, 80(6): WD117-WD128. DOI:10.1190/geo2014-0234.1
Liu Y, Fomel S, Liu C, et al. 2009. High-order seislet transform and its application of random noise attenuation. Chinese Journal of Geophysics (in Chinese), 52(8): 2142-2151. DOI:10.3969/j.issn.0001-5733.2009.08.024
Liu Y, Li B X. 2018. Streaming orthogonal prediction filter in the t-x domain for random noise attenuation. Geophysics, 83(4): F41-F48. DOI:10.1190/geo2017-0322.1
Liu Y, Li B X, Wang D, et al. 2017. Local SNR estimation method based on regularization for seismic data. Chinese Journal of Geophysics (in Chinese), 60(5): 1979-1987. DOI:10.6038/cjg20170529
Liu Y, Zhang P, Liu C. 2017. Seismic data interpolation using generalised velocity-dependent seislet transform. Geophysical Prospecting, 65(S1): 82-93.
Mallat S, Zhong S. 1992. Characterization of signals from multiscale edges. IEEE Transactions on Pattern Analysis and Machine Intelligence, 14(7): 710-732. DOI:10.1109/34.142909
Mallat S G, Zhang Z F. 1993. Matching pursuits with time-frequency dictionaries. IEEE Transactions on Signal Processing, 41(12): 3397-3415. DOI:10.1109/78.258082
Nan F Z, Xu Y, Liu W, et al. 2018. Denoising methods of OBS data based on sparse representation. Chinese Journal of Geophysics (in Chinese), 61(4): 1519-1528. DOI:10.6038/cjg2018L0130
Pope G. 2008. Compressive Sensing: A Summary of Reconstruction Algorithms. Zürich: Eidgen ssische Technische Hochschule.
Tong S Y, Gao H, Liu R, et al. 2019. Seismic random noise adaptive suppression based on the Shearlet transform. Oil Geophysical Prospecting (in Chinese), 54(4): 744-750.
Vassiliou A A, Garossino P. 1998-12-15. Time-frequency processing and analysis of seismic data using very short-time Fourier transforms: US, 5850622.
Wright S J, Nowak R D, Figueiredo M A T. 2008. Sparse reconstruction by separable approximation. //2008 IEEE International Conference on Acoustics, Speech and Signal Processing. Las Vegas, NV, USA: IEEE, 3373-3376.
Wu Z C, Liu T Y. 2008. Wavelet transform methods in seismic data noise attenuation. Progress in Geophysics (in Chinese), 23(2): 493-499.
Yang G Y, Luan X W, Meng F S, et al. 2020. Seismic data reconstruction based on Shearlet transform and total generalized variation regularization. Chinese Journal of Geophysics (in Chinese), 63(9): 3465-3477. DOI:10.6038/cjg2020M0424
Yin W T, Osher S, Goldfarb D, et al. 2008. Bregman iterative algorithms for L1-minimization with applications to compressed sensing. SIAM Journal on Imaging Sciences, 1(1): 143-168. DOI:10.1137/070703983
Zhong T. 2016. A study on the characteristics of land-seismic-prospecting random noise based on modern statistical theory[Ph. D. thesis]. Changchun: Jilin University.
Zhong T, Li Y, Yang B J, et al. 2017. Statistical features of the random noise in land seismic prospecting. Chinese Journal of Geophysics (in Chinese), 60(2): 655-664. DOI:10.6038/cjg20170219
陈兴飞, 孙红梅. 2019. 基于压缩感知理论的地震数据降噪方法. 地球物理学进展, 34(3): 1025-1031. DOI:10.6038/pg2019CC0208
邓盾, 裴健翔, 邓勇, 等. 2019. 南海深水盆地改善中深层地震成像关键技术研究. 地球物理学进展, 34(2): 702-708. DOI:10.6038/pg2019BB0551
高静怀, 毛剑, 满蔚仕, 等. 2006. 叠前地震资料噪声衰减的小波域方法研究. 地球物理学报, 49(4): 1155-1163. DOI:10.3321/j.issn:0001-5733.2006.04.030
巩向博, 韩立国, 王恩利, 等. 2009. 压制噪声的高分辨率Radon变换法. 吉林大学学报(地球科学版), 39(1): 152-157.
国胧予, 刘财, 刘洋, 等. 2020. 基于f-x域流式预测滤波器的地震随机噪声衰减方法. 地球物理学报, 63(1): 329-338. DOI:10.6038/cjg2020N0030
李海山, 陈德武, 吴杰, 等. 2020. 叠前随机噪声深度残差网络压制方法. 石油地球物理勘探, 55(3): 493-503.
李庆忠. 1993. 走向精确勘探的道路. 北京: 石油工业出版社: 112-114.
刘洋, Fomel S, 刘财, 等. 2009. 高阶seislet变换及其在随机噪声消除中的应用. 地球物理学报, 52(8): 2142-2151. DOI:10.3969/j.issn.0001-5733.2009.08.024
刘洋, 李炳秀, 王典, 等. 2017. 基于正则化条件的地震数据局部信噪比估计方法. 地球物理学报, 60(5): 1979-1987. DOI:10.6038/cjg20170529
南方舟, 徐亚, 刘伟, 等. 2018. 基于稀疏表达的OBS去噪方法. 地球物理学报, 61(4): 1519-1528. DOI:10.6038/cjg2018L0130
童思友, 高航, 刘锐, 等. 2019. 基于Shearlet变换的自适应地震资料随机噪声压制. 石油地球物理勘探, 54(4): 744-750.
吴招才, 刘天佑. 2008. 地震数据去噪中的小波方法. 地球物理学进展, 23(2): 493-499.
杨冠雨, 栾锡武, 孟凡顺, 等. 2020. 基于Shearlet变换和广义全变分正则化的地震数据重建. 地球物理学报, 63(9): 3465-3477. DOI:10.6038/cjg2020M0424
钟铁. 2016. 基于现代统计理论的陆地地震勘探随机噪声特征研究[博士论文]. 长春: 吉林大学.
钟铁, 李月, 杨宝俊, 等. 2017. 陆地地震勘探随机噪声统计特性. 地球物理学报, 60(2): 655-664. DOI:10.6038/cjg20170219