地球物理学报  2019, Vol. 62 Issue (3): 1181-1192   PDF    
低信噪比地震资料FDOC-seislet变换阈值消噪方法研究
张雅晨1,4, 刘洋1,2,3, 刘财1,2,3,4, 武尚1     
1. 吉林大学地球探测科学与技术学院, 长春 130026;
2. 吉林大学应用地球物理实验教学中心, 长春 130026;
3. 吉林大学地质资源立体探测虚拟仿真实验教学中心, 长春 130026;
4. 国土资源部应用地球物理重点实验室, 长春 130026
摘要:地震数据本质上是时变的,不仅有效同相轴表现出确定性信号的时变特征,而且复杂地表和构造条件以及深部探测环境总是引入时变的非平稳随机噪声.标准的频率-空间域预测滤波只适合压制平面波信号假设下的平稳随机噪声,而处理非平稳地震随机噪声时,需要将数据体分割为小窗口进行分析,但效果不够理想,而传统非预测类随机噪声压制方法往往适应性不高,因此开发能够保护地震信号时变特征的随机噪声压制方法具有重要的工业价值.压缩感知是近年出现的一个新的采样理论,通过开发信号的稀疏特性,已经在地震数据处理中的数据插值以及噪声压制中得到了应用.本文系统地分析了压缩感知理论框架下的地震随机噪声压制问题,建立了阈值消噪的数学反演目标函数;针对时变有效信息具有的可压缩性,利用有限差分算法求解炮检距连续方程,构建有限差分炮检距连续预测算子(FDOC),在seislet变换框架下,提出一种新的快速稀疏变换域——FDOC-seislet变换,实现地震数据的高度稀疏表征;结合非平稳随机噪声不可压缩的特征,提出了一种整形迭代消噪方法,该方法是一种广义的迭代收缩阈值(IST)算法,在无法计算稀疏变换伴随算子的条件下,仍然能够对强噪声环境中的时变有效信息进行有效恢复.通过对模型数据和实际数据的处理,验证了FDOC-seislet稀疏变换域随机噪声迭代压制方法能够在保护复杂构造地震波信息的前提下,有效地衰减原始数据中的强振幅随机噪声干扰.
关键词: 压缩感知理论      炮检距连续方程      FDOC-seislet变换      迭代整形消噪      迭代收缩阈值     
Seismic random noise attenuation using FDOC-seislet transform and threshold for seismic data with low SNR
ZHANG YaChen1,4, LIU Yang1,2,3, LIU Cai1,2,3,4, WU Shang1     
1. College of Geo-exploration Science and Technology, Jilin University, Changchun 130026, China;
2. Central Lab of Applied Geophysics, Changchun 130026, China;
3. Virtual Simulation Experiment Teaching Center for Stereoscopic Exploration of Geological Resources, Changchun 130026, China;
4. Key Laboratory of Applied Geophysics, Ministry of Land and Resources, Changchun 130026, China
Abstract: In seismic exploration, seismic data are essentially time-varying, not only effective events display time-varying characteristics, but also complex surface and subsurface conditions and deep exploration environment always create time-varying and nonstationary random noise.Industrial standard prediction filter in the frequency-space (f-x) domain can only suppress stationary random noise under the assumption of stationary plane waves, for which cutting data into overlapping windows is a common method to handle nonstationary problem. However, such an approach cannot produce good results.On the other hand, random noise attenuation methods without prediction cannot provide good adaptability for different data.It is important to develop nonstationary random noise suppression methods that are suitable for time-varying characteristics of seismic data.Compressive sensing is a new sampling theory, which has found different applications in seismic data processing, for example, missing data interpolation and noise attenuation.In this paper, we analyze the problem of seismic random noise attenuation under compressive sensing framework and establish the objective function for threshold denoising.Aiming at the compressible characteristics of time-varying seismic signal, we solve the offset continuation equation by using finite difference algorithm to construct a finite-difference offset-continuation operato. Then a new fast sparse transform, FDOC-seislet transform, is proposed under seislet transform framework, which provides a highly sparse representation of seismic signals.According to the feature that nonstationary random noise is less compressible, we propose a new iterative shaping denoising method, a kind of generalized iterative shrinkage threshold methods, which is able to recover the time-varying signal under strong random noise environment, even when the adjoint operator of sparse transform cannot be calculated.Compared with conventional denoising methods, the examples of synthetic and field data demonstrate that the random noise attenuation method based on FDOC-seislet transform and iterative shaping can protect complex seismic signal and attenuate random noise even with strong amplitudes.
Keywords: Compressive sensing    Offset continuation equation    FDOC-seislet transform    Iterative shaping denoising    Iterative shrinkage threshold    
0 引言

人工地震是地球物理勘探中重要的方法之一,在保证精度和分辨率的情况下,探测深度甚至能够达到莫霍面.目前的石油、天然气等化石能源的确定以地震勘探为主要手段,随着国家重点研发计划“深地资源勘查开采”重点专项的开展,人工地震方法在地球深部资源探查中也成为重要方法.人工地震数据本质上是时变的,其特征表现在地震同相轴的能量、轨迹、时频谱等特征以及随机噪声的统计特征随时间和空间位置不同而变化.一方面,尽管王宏禹和邱天爽(2015)给出了确定性信号的非平稳性定义,但这种概念仍未得到广泛的应用.地震勘探有效同相轴相应的道信号可归类为确定性信号,这种信号是时变而不随机,其也不满足平稳性质,但就数学而言,用平稳性数字特性也可以计算,但不是物理特征,因此本文使用时变性来描述确定性信号的变化特征.另一方面,随机噪声广泛存在于数据中,严重影响着地震资料的分析,在山前带和高陡构造等地区,复杂地下介质容易造成有效信息的能量损失,地表条件等引起的随机干扰会导致低品质数据的采集;另外,随着探测深度的增加,有效信息的能量较弱,随机噪声使得深部地震数据信噪比普遍较低.地震数据中的随机噪声被认为具有平稳随机性(陆基孟,1991李庆忠,1993),随着勘探问题的深入研究,地震随机噪声被发现不是传统意义上的平稳随机过程,实际的地表和构造条件以及深地震观测系统总是会采集到非平稳随机噪声(钟铁等,2017),这种非平稳随机噪声也体现出时变特征.因此,实际地震数据是一个混合信号,即时变确定性信号与时变非平稳随机信号的线性合成信号,在一般数学工具书中都未涉及这类合成信号,物理上用平稳性描述也不合理的,但这个线性混合信号总体表现出时变性.非平稳随机噪声引起的低信噪比问题很难解决,传统的噪声压制方法往往不能有效分离时变确定性地震信号和时变非平稳随机噪声,最终将会引起构造解释精度的下降,甚至导致数据反演的失败.因此,开发有效的非平稳随机噪声压制方法具有重要的实际意义.

许多学者开发了不同类型的地震随机噪声压制方法.(1)利用随机噪声的统计学特征开发的处理方法.包括改进的均值滤波类(Bonar and Sacchi, 2012)和中值滤波类(王伟等, 2012, Liu et al.,2009)方法等.(2)基于各种数学变换的处理方法.包括傅里叶变换(Zhou et al.,2015),余弦变换(陆文凯,2011)以及类小波变换(Wang et al.,2015汪金菊等,2016Chen et al.,2016Zhao et al.,2016)等.(3)基于数学分析的处理方法.例如:经验模式分解(EMD)(Bekara and van der Baan,2009),奇异谱分析(Huang et al.,2016),SVD类方法(刘志鹏等,2012)等.金雷等(2005)将信号处理时频峰值滤波方法引入地震随机噪声压制问题中,其他学者进行了一系列的改进(Liu et al.,2014;林洪波等,2015;Zhang et al.,2016).(4)基于数据轨迹特征的处理方法.包括方向性滤波(黄梅红和李月,2016),线性拟合(Lu and Lu, 2009)以及各项异性分数梯度算子(Wang and Gao, 2014)等.(5)预测滤波类处理方法.地震同相轴具有空间方向连续性,可以利用这种可预测特征从噪声中分辨出有效信息.如果同相轴轨迹是线性分布的,可以使用频率-空间域(f-x域)或时间-空间域(t-x域)的预测技术分析同相轴信息.f-x预测滤波技术由Canales(1984)提出并且由Gulunay(1986)进一步发展,已经成为一种工业界标准方法,被称为“FXDECON”.Naghizadeh和Sacchi(2012)提出多分量f-x矢量自回归算子压制地震随机噪声.Liu等(2012)提出了基于正则化非平稳自回归(RNA)的f-x域预测滤波方法压制随机噪声.预测滤波也可以在t-x域计算,Abma和Claerbout(1995)讨论了利用f-x域和t-x域预测方法从随机噪声中恢复线性同相轴的差异,证明t-x域预测滤波压制随机噪声能力更强.Sacchi和Naghizadeh(2009)提出通过计算时间和空间方向信号变化的预测滤波器压制噪声,尽管上述方法在处理不同地区的随机噪声时能够取得一定效果,但是仍然无法解决低信噪比条件下时变地震信号的准确恢复问题.

压缩感知(Donoho,2006)是近年出现的一个新的采样理论,通过开发信号的稀疏特性,在远小于Nyquist采样率的条件下,用随机采样获取信号的离散样本,然后通过非线性重建算法重建信号.该理论一经提出,就引起学术界和工业界的广泛关注,并且已经在地震数据噪声压制中得到了应用.压缩感知理论应用的基础是稀疏变换域和迭代反演算法.近年来,结合Fourier变换(Abma and Kabir, 2006Zwartjes and Sacchi, 2007唐刚和杨慧珠,2010)和curvelet变换(曹静杰等,2012)的稀疏反演方法在地震数据处理中取得了较好的效果,其理论基础都源于压缩感知理论.不过由于这些数学变换都不是特殊针对地震数据开发的方法,并不能为地震数据提供最佳的稀疏表征,因此在应用压缩感知理论解决地震数据处理问题时,存在一定的限制.针对地震数据的特有性质,Fomel和Liu(2010)提出了一种基于地震数据属性特征的seislet类小波变换,刘洋等(2009)开发了seislet变换的高阶形式,能够对沿同相轴方向的复杂信号(例如带有AVO效应)进行有效地压缩.然而,复杂地震波场(尤其含散射波场)的动力学表征很难被早期的seislet变换所解决.Liu和Fomel(2010)进一步扩展seislet变换,通过结合炮检距连续算子(Deregowski and Rocca, 1981Bolondi et al.,1982Salvador and Savelli, 1982Fomel,2003a)开发了傅里叶炮检距连续seislet变换(FTOC-seislet)方法,尽管该方法对于复杂地震波场有较好的压缩效果,但是,其利用全波场信息,计算复杂性和计算量都比较大,不适合用于解决叠前大规模地震数据处理任务.另一方面,在迭代反演算法中,常见有基追踪算法、匹配追踪算法、正交匹配追踪法等,其中迭代阈值算法(Daubechies et al.,2004)是最常见的一种方法,但是该方法需要计算稀疏变换算子的伴随矩阵,在一定程度上限制了该方法的应用.

本文是在国家自然科学基金重点项目“大兴安岭西盆地群域构造与地球物理场”研究中,解决低品质地震数据的信号提取问题而开展的深入研究.针对时变有效信息具有可压缩性而非平稳随机噪声不可压缩的特征,利用有限差分算法求解炮检距连续方程,在seislet变换框架下,构建新的能够准确表征强随机噪声干扰下非全波场地震信息的FDOC-seislet稀疏变换,扩展传统迭代收缩阈值方法,提出一种广义迭代整形消噪算法,对强噪声环境中的时变复杂有效信息进行准确恢复.通过对模型数据和实际数据的处理,验证了本次研究方法能够在压制原始数据中强振幅随机噪声的同时最大程度地保护复杂构造下的地震波场信息.

1 理论基础 1.1 压缩感知框架下的随机噪声压制问题

随机噪声压制问题广泛存在于各个领域的数据处理任务,在早期的研究中,随机噪声压制问题可以通过稀疏域阈值去噪模型进行求解(Mallat and Hwang, 1992Donoho,1995).例如,在地震数据中,随机噪声一般为随机分布,如果存在一个稀疏变换域A,能够将地震信号压缩在有限的几个系数中,而不可压缩的噪声能量分布于整个稀疏变换域内,可以找到一个合适的阈值,当变换域系数小于该阈值时,认为这时的系数主要是由噪声引起的;当变换域系数大于该阈值时,则认为其主要是由地震信号引起的.当对变换域系数进行阈值处理后,反变换就可以达到去除噪声而保留有用信号的目的.经典的软阈值去噪方法主要解决如下最优化问题:

(1)

其中,y=Ab为含噪声数据的变换域系数,b为含噪数据,A为稀疏变换,x为待求的去噪数据的变换域系数,方程(1)假设在变换域中随机噪声为高斯分布,当随机噪声满足数据域内高斯分布时,需要建立不同的目标函数.

本质上,地震数据随机噪声压制可以由线性估计问题引出.假设b为观测数据,m为待求解的模型,利用正向模型算子L可以获得模型和观测数据之间的方程

(2)

当正向模型算子L为单位矩阵I时,并且如果模型空间是可压缩的(即具有极少的非零元素),那么压缩感知理论(Donoho,2006)可以将求解公式(2)转化为如下问题:

(3)

其中,‖ ‖0表示L0范数,即数据的非零个数.数学非凸优化问题公式(3)(Natarajan,1995)往往很难得到全局最优解,因此在数学上往往把公式(3)转化为凸优化问题:

(4)

其中,‖ ‖1表示L1范数,Elad and Bruckstein(2002)证明了公式(3)和(4)的等效条件.

将约束问题(公式(4))进一步转化为非约束问题:

(5)

其中,μ为惩罚参数.此时公式(5)描述了随机噪声消减的数学问题.然而,消噪后的地震数据m很少存在可直接压缩(只有少数元素为非零)的情况.为了将反问题(2)在压缩感知理论框架下进行求解,需要对m进行变量代换.如果m在某种变换域内具有可压缩的特性,即

(6)

其中,A-1为稀疏反变换,x为变换域系数,同时x的非零数据个数远远小于m的非零数据个数.利用x的稀疏特性以及公式(6)对公式(3)—(5)进行修改,可得基追踪降噪问题的目标函数(Chen et al.,2001)

(7)

公式(7)与公式(1)具有本质的不同,即随机噪声分别在数据域和稀疏变换域内具有高斯分布.在求解公式(7)时,关键在于稀疏变换A的设计和迭代算法的选取.

1.2 FDOC-Seislet变换

Seislet变换是一种利用地震数据模式特征压缩数据的稀疏变换方法;以小波提升算法为基础,与地震数据模式识别相结合,使用不同的模式对地震数据进行预测,能够得到不同类型的seislet变换.例如,沿着地震同相轴倾角方向表征地震数据,可以得到基本的seislet变换的定义(Fomel and Liu, 2010刘洋等,2009),然而,单一数据点位置上的交叉同相轴(如散射波存在情况)很难被早期的seislet变换准确表征.地震数据有别于一般的图像信息,其在固体介质中传播的物理过程基于地震波动力学理论,并且采用特殊的获取手段(地面接收的观测系统),波动方程能够描述复杂地震波场的传递关系.这种波场的可预测特性与seislet变换框架相结合,可以用来定义新的FDOC-seislet稀疏变换方法A.

Seislet变换是由离散小波变换和地震数据模式构成的,其中离散小波变换基于提升算法(Sweldens,1995).提升算法的核心是更新算法和预测算法,通过预测算子得到高频信息,而通过更新算子得到低频信息,正变换的步骤如下:

(1) 将数据分解为奇数序列o和偶数序列e.

(2) 计算奇数序列与偶数序列预测值之间残差r以及数据的近似值c

(8)

其中,P是预测算子,U是更新算子.

(3) 近似值c成为新的数据,重复以上步骤得到下一级数的变换系数.

设计不同的预测算子和更新算子可以产生不同的小波变换,当选取两个相邻数据样点的线性插值为预测算子,以及保留信号的连续平均值作为更新算子时,可以得到提升方案下Cohen-Daubechies-Feauvea(CDF)5/3双正交小波变换算子:

(9)

(10)

其中,k是奇数元素的位置.反变换可以通过反向计算以上步骤来恢复原始数据.

对小波提升方案中离散小波变换的预测算子和更新算子进行修改,通过结合炮检距连续算子进行预测运算.同样以CDF 5/3双正交小波为例,可以得到FDOC-seislet变换的更新算子和预测算子:

(11)

(12)

其中:表示频率-空间域数据在沿炮检距坐标轴方向的偶数位置数据向量;表示频率-空间域数据在沿炮检距坐标轴方向的偶数位置数据向量预测值与奇数位置数据向量之间的残差;是有限差分炮检距连续(FDOC)预测算子,+和-分别表示沿共炮检距方向从左、右两侧对频率域数据进行预测,k是频率-空间域数据在沿炮检距坐标轴h方向的奇数数据位置.接下来给出FDOC预测算子的具体形式.

在炮检距连续微分理论中,炮检距连续偏微分方程有如下形式(Fomel,2003b):

(13)

其中,U(x, h, tn)是地震数据,xhtn分别为中心点坐标、炮检距坐标以及NMO(正常时差校正)后的时间轴,n为该时间轴上的采样点,h=(r-s)/2是炮检距坐标,x=(r+s)/2是共中心点坐标,rs分别为检波点和炮点位置.利用全通滤波器的谱能量不变性,在Z变换下求解方程(13),有如下解的形式:

(14)

公式(14)代表一个隐式有限差分格式,其中:

(15)

(16)

(17)

(18)

Zx为共中心点坐标xZ变换.最终,h1h2之间的有限差分炮检距连续预测算子为

(19)

与傅里叶炮检距连续预测算子(Liu and Fomel, 2010)同时依赖于炮检距和共中心点坐标的计算方式不同,有限差分炮检距连续预测算子的计算在共中心点坐标方向上互相独立,因此能够更好地实现中心点方向数据分割处理,为后期的三维大规模地震数据处理应用提供可行性.

1.3 整形迭代方法

为了求解非约束问题(公式(7)),可以从公式

(20)

出发,如果算子A-H可以看作A的近似算子,即级联算子A-HA-1约等于单位算子I,那么存在简单的迭代方法

(21)

公式(21)收敛于如下方程的解

(22)

公式(21)被Goldin(1986)称为“R方法”,其收敛的充分条件为公式(22)右侧表达式是稀疏的.当A-HA-1的伴随算子时,公式(21)也称为“Landweber迭代”(Landweber,1951).由于公式(21)往往不满足所期待的模型解,Fomel(2008)提出了一种非线性整形方法来约束模型空间的属性,通过定义一个整形算子T,公式(21)被修改为

(23)

其中,整形算子T可以选取为非线性阈值算子(Donoho and Johnstone, 1994).当公式中A-H满足A-1伴随算子的条件时,整形迭代公式(23)等价为Daubechies等(2004)提出的经典收缩迭代方法,即迭代收缩阈值方法(IST),公式(23)最终收敛于公式(7)的解.该方法假设数据域内的随机噪声为高斯白噪,即满足噪声在数据域内L2范数下最小,最终的消噪结果为m=A-1xNN为迭代次数.在本文中,FDOC-seislet变换A为线性变换,其反变换A-1的伴随算子A-H可以构建,因此,在数据测试中,迭代整形消噪算法退化为经典的迭代收缩阈值算法.

2 模型测试 2.1 正弦反射界面模型

首先,建立如图 1所示的正弦反射界面模型,通过该模型来说明FDOC-seislet的地震数据压缩特性以及信噪分离过程.利用克希霍夫方法正演地震数据,如图 2a所示.图 2a所示立方体数据的三个表面分别代表立方体内部的三个切片(即上表面代表 0.516 s的时间切片,前表面代表半炮检距0.256 km位置的共炮检距剖面,右表面代表中心点1.6 km位置的共中心点道集),从图中可以看到,正弦反射界面产生了复杂的反射波场,许多位置出现了反射波同相轴交叉的情况,反射界面向斜构造甚至会引起远炮检距位置三条同相轴同时汇聚的特征,这种复杂的波场使传统的噪声压制方法很难获得理想的处理效果.为了测试FDOC-seislet的随机噪声压制能力,在图 2a中加入反射波最大振幅1.6倍左右的高斯随机噪声,图 2b为加入噪声后的数据体,如此极端的强噪声情况在常规地震勘探中较少出现,此处仅为了测试本文方法的有效性.此时反射波已经淹没在随机噪声中,原有的波场特征无法被分辨.利用FDOC-seislet变换方法对数据体(图 2b)进行压缩处理,图 3a仅显示频率域内沿炮检距方向进行压缩的变换域系数实部,可以看到,有效反射波信息被压缩在第1个级数范围内,而随机噪声分布在所有变换域系数位置,由于压缩后的有效波信息变换域系数振幅要明显大于随机噪声的变换域系数,因此可以利用软阈值在变换域内进行信噪分离.软阈值处理后的变换域系数如图 3b所示,图 3b为处理后仅保留振幅值较大的有效波FDOC-seislet变换系数.对图 3b进行反变换,可以得到噪声压制后的有效波场信息,如图 4a所示.对比图 4a图 2b,不难看到被淹没的反射波信息被有效地恢复,尤其是交叉同相轴也被较好地保护.通过图 4a图 2a中有效波振幅的关系,可以看到噪声压制后的有效波振幅要小于原始信号,主要是因为软阈值无法区分变换域内弱信号和噪声,部分有效波能量也被当作噪声被衰减,但是在分离的噪声中(图 4b),未发现明显的反射波信息,说明被衰减的有效信号能量较小.

图 1 正弦反射界面模型 Fig. 1 Model with sinusoid reflector
图 2 正弦模型数据 (a)理想数据体;(b)加入随机噪声后的数据体. Fig. 2 Model data (a) Ideal data; (b) Data added random noise.
图 3 FDOC-seislet变换域系数 (a)含噪声数据的变换域系数;(b)软阈值处理后的变换域系数. Fig. 3 Coefficients in FDOC-seislet transform domain (a) Coefficients of noisy data in transform domain; (b) Coefficients of processed data in transform domain.
图 4 软阈值处理结果 (a)噪声压制后的信号;(b)分离的噪声成分. Fig. 4 Processed data with soft thresholding (a) The denoised result; (b) The removed noise.
2.2 French模型

为了对比FDOC-seislet变换迭代去噪方法与传统方法的处理效果,选取标准French模型(French,1974)(图 5a)的二维切片(图 5b),通过克希霍夫正演建立一个时间-共中心点-炮检距域二维叠前数据体(图 6a),三个剖面分别选取为:0.6 s的时间切片、1.0 km的共炮检距剖面以及0.2 km的共中心点剖面.图 6a的地震波同相轴信息十分复杂,在构造拐点处有较强的散射波存在,并且复杂构造使得散射波与反射波同相轴有互相叠加的情况.对图 6a数据加上较强随机噪声,结果如图 6b所示.选取的高斯随机噪声最大振幅值略大于有效信号的最大振幅,这种幅度的随机噪声与实际数据情况更加接近.但是噪声仍然对地震资料的中深层弱反射波和散射波产生较大的影响,该噪声模型的处理难点在于同时处理多维、振幅时变以及交叉同相轴.在此复杂波场条件下,很难用简单的消噪方法实现对该类型噪声的压制.为了比较,本文使用工业标准三维FXDECON方法压制随机噪声,结果如图 7a所示,三维FXDECON方法能够消除大部分的随机干扰,但是,该方法对弯曲同相轴产生了比较大的损害,尤其在共炮检距道集中最下层的同相轴连续性被破坏,这是因为FXDECON方法是一种非时变滤波方法,很难处理弯曲的同相轴情况,分时窗处理能够改善这种现象,但是此时需要大量的人工参数选取,并且窗口参数无法进行自适应调节,处理效果仍然有限.图 7b为应用本文方法对噪声数据的处理结果,FDOC-seislet变换阈值消噪方法能够压制更多的随机噪声,同时复杂的有效波场信息也被很好地保留.

图 5 标准French模型 (a)三维模型;(b)二维模型切片. Fig. 5 Benchmark French model (a) 3D model; (b) 2D slice of model.
图 6 French模型数据 (a)理想数据体;(b)加入随机噪声后的数据体. Fig. 6 Synthetic data of French model (a) Ideal data; (b) Noisy data.
图 7 不同方法噪声压制对比 (a)三维FXDECON方法;(b)本文方法. Fig. 7 The denoised results using different methods (a) 3D FXDECON; (b) The proposed method.
3 实际数据测试

为了验证本方法的有效性,在实际数据测试中使用波兰某区域的陆上实际数据来进行测试.该数据具有较低的信噪比,图 8a为去除大部分面波后的共中心点道集,强随机噪声使有效反射波的特征不可见.尽管如此,常规速度扫描仍然能够用于得到比较合理的动校正速度(图 8b),速度场将用于本文方法的数据预处理.试验证明,本文方法对动校正速度并不敏感,因此,即使动校正速度有一定的误差,对处理的效果也产生较小的影响.利用本文提出的方法对该数据体进行随机噪声压制,由于该数据体的覆盖次数较低,因此在FDOC-seislet变换域内,部分能量较弱的有效信号的变换域系数与随机噪声的变换域系数接近,采用软阈值处理则会衰减有效信息,根据有效波信息压缩后的能量集中在较小级数尺度上的特点,在实际数据中使用级数方向的线性窗函数代替公式(23)中的软阈值,仅保留前4个级数尺度上的信号.通过迭代处理后,得到随机噪声压制后的数据(图 9b),从共炮检距道集中可以看到,数据的信噪比得到提高,能够看到比较明显的反射同相轴的特征,同时在时间切面上也能观察到明显的同相轴信息.

图 8 波兰某地区实际陆上数据 (a) CMP道集;(b)速度谱. Fig. 8 Field land data in Poland (a) CMP gathers; (b) Velocity spectra.

为了验证本文方法对实际数据的消噪效果,选取工业标准FXDECON方法与其进行对比(图 9a),可以观察到,FXDECON方法强随机噪声的压制能力有限,仅能够压制一些近炮检距的干扰,消噪后有效信号能量仍然较弱,主要是频率-空间域预测滤波类方法要求数据的信噪比不能过低,否则会引起滤波器无法对有效信息进行估计的情况.为了进一步对比,对两种消噪方法的处理结果进行叠加,图 10a图 10b分别为两种方法处理后的叠加剖面,图中可以发现:传统FXDECON方法变换得到的叠加剖面信噪比低,中深层的有效波信息不能被有效识别;而在FDOC-seislet变换迭代消噪结果对应的叠加剖面中,整体剖面的质量有很大的提高,尤其是0.5~2.0 s之间的地层连续性被提高,验证了本文方法可以达到压制随机噪声的效果,尤其恢复低信噪比条件下有效信号的能量.

图 9 不同方法的处理结果比较 (a) FXDECON;(b)本文方法. Fig. 9 The denoised results using different methods (a) FXDECON; (b) The proposed method.
图 10 不同处理方法的叠加剖面比较 (a) FXDECON;(b)本文方法. Fig. 10 The stacking results using different methods (a) FXDECON; (b) The proposed method.
4 结论

本文介绍了一种基于FDOC-seislet变换的随机噪声迭代阈值衰减方法.该方法以波动方程动力学关系为基础,能够准确地表征低信噪比条件下复杂地震波场,有效地压缩复杂地震波同相轴,在压缩感知理论框架下,通过整形迭代算法改善地震资料的信噪比.理论模型测试结果表明,本文方法能够有效压制含散射波复杂波场的强随机噪声,证明了方法的有效性.实际地震资料处理结果表明本文方法能够有效压制随机噪声,恢复低信噪比条件下有效波场特征,为后续的处理工作提供基础资料.

致谢  感谢吉林大学杨宝俊教授的讨论和建议。感谢匿名审稿专家提出的宝贵意见.
References
Abma R, Claerbout J. 1995. Lateral prediction for noise attenuation by t-x and f-x techniques. Geophysics, 60(6): 1887-1896. DOI:10.1190/1.1443920
Abma R, Kabir N. 2006. 3D interpolation of irregular data with a POCS algorithm. Geophysics, 71(6): E91-E97. DOI:10.1190/1.2356088
Bekara M, van der Baan M. 2009. Random and coherent noise attenuation by empirical mode decomposition. Geophysics, 74(5): V89-V98. DOI:10.1190/1.3157244
Bolondi G, Loinger E, Rocca F. 1982. Offset continuation of seismic sections. Geophysical Prospecting, 30(6): 813-828. DOI:10.1111/gpr.1982.30.issue-6
Bonar D, Sacchi M. 2012. Denoising seismic data using the nonlocal means algorithm. Geophysics, 77(1): A5-A8.
Canales L L. 1984. Random noise reduction.//54th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 525-527.
Cao J J, Wang Y F, Yang C C. 2012. Seismic data restoration based on compressive sensing using the regularization and zero-norm sparse optimization. Chinese Journal of Geophysics (in Chinese), 55(2): 596-607. DOI:10.6038/j.issn.0001-5733.2012.02.022
Chen S S, Donoho D L, Saunders M A. 2001. Atomic decomposition by basis pursuit. SIAM Review, 43(1): 29-159.
Chen Y K, Ma J W, Fomel S. 2016. Double sparsity dictionary for seismic noise attenuation. Geophysics, 81(2): V193-V206.
Daubechies I, Defries 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/(ISSN)1097-0312
Deregowski S M, Rocca F. 1981. Geometrical optics and wave theory of constant offset sections in layered median. Geophysical Prospecting, 29(3): 374-406. DOI:10.1111/gpr.1981.29.issue-3
Donoho D L, Johnstone I M. 1994. Ideal spatial adaptation by wavelet shrinkage. Biometrika, 81(3): 425-455. DOI:10.1093/biomet/81.3.425
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. 2006. Compressed sensing. IEEE Transactions on Information Theory, 52(4): 1289-1306. DOI:10.1109/TIT.2006.871582
Elad M, Bruckstein A M. 2002. A generalized uncertainty principle and sparse representation in Pairs of bases. IEEE Transactions on Information Theory, 48(9): 2558-2567. DOI:10.1109/TIT.2002.801410
Fomel S. 2003a. Theory of differential offset continuation. Geophysics, 68(2): 718-732.
Fomel S. 2003b. Seismic reflection data interpolation with differential offset and shot continuation. Geophysics, 68(2): 733-744.
Fomel S. 2008. Nonlinear shaping regularization in geophysical inverse problems.//78th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 2046-2051.
Fomel S, Liu Y. 2010. Seislet transform and seislet frame. Geophysics, 75(3): V25-V38.
French W S. 1974. Two-dimensional and three-dimensional migration of model-experiment reflection profiles. Geophysics, 39(3): 265-277.
Goldin S V. 1986. Seismic traveltime inversion.SEG.
Gulunay N. 1986. FXDECON and complex Wiener prediction filter.//56th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts.
Huang M H, Li Y. 2016. Suppression of seismic random noise based on steerable filter. Chinese Journal of Geophysics (in Chinese), 59(5): 1815-1823. DOI:10.6038/cjg20160524
Huang W L, Wang R Q, Chen Y K, et al. 2016. Damped multichannel singular spectrum analysis for 3D random noise attenuation. Geophysics, 80(4): V261-V270.
Jin L, Li Y, Yang B J. 2005. Reduction of random noise for seismic data by time-frequency peak filtering. Progress in Geophysics (in Chinese), 20(3): 724-728.
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 Q Z. 1993. A Systematical Analysis of High Resolution Seismic Exploration Seismic Exploration (in Chinese). Beijing: Petroleum Industry Press: 112-114.
Lin H B, Ma H T, Xu L P. 2015. A radial time-frequency peak filtering based on ROAD for suppressing spatially nonstationary random noise in seismic data. Chinese Journal of Geophysics (in Chinese), 58(7): 2546-2555. DOI:10.6038/cjg20150729
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 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, Liu C, Wang D. 2009. A 1D time-varying median filter for seismic random, spike-like noise elimination. Geophysics, 74(1): V17-V24.
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 P, Dang B, Li Y, et al. 2014. Local spatiotemporal time-frequency peak filtering method for seismic random noise reduction. Journal of Applied Geophysics, 111: 76-85. DOI:10.1016/j.jappgeo.2014.09.018
Liu Z P, Zhao W, Chen X H, et al. 2012. Local SVD for random noise suppression of seismic data in frequency domain. Oil Geophysical Prospecting (in Chinese), 47(2): 202-206.
Lu J M. 1991. Principles and Data Interpretation of Seismic Prospecting (in Chinese). Beijing: Petroleum Industry Press: 140-142.
Lu W K. 2011. Seismic random noise suppression based on the discrete cosine transform. Oil Geophysical Prospecting (in Chinese), 46(2): 202-206.
Lu Y H, Lu W K. 2009. Edge-preserving polynomial fitting method to suppress random seismic noise. Geophysics, 74(4): V69-V73. DOI:10.1190/1.3129907
Mallat S, Hwang W L. 1992. Singularity detection and processing with wavelets. IEEE Transactions on Information Theory, 38(2): 617-643. DOI:10.1109/18.119727
Naghizadeh M, Sacchi M. 2012. Multicomponent f-x seismic random noise attenuation via vector autoregressive operators. Geophysics, 77(2): V91-V99. DOI:10.1190/geo2011-0198.1
Natarajan B K. 1995. Sparse approximate solutions to linear systems. SIAM Journal on Computing, 24(2): 227-234. DOI:10.1137/S0097539792240406
Sacchi M, Naghizadeh M. 2009. Adaptive linear prediction filtering for random noise attenuation.//79th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 3347-3351.
Salvador L, Savelli S. 1982. Offset continuation for seismic stacking. Geophysical Prospecting, 30(6): 829-849. DOI:10.1111/gpr.1982.30.issue-6
Sweldens W. 1995. Lifting scheme:A new philosophy in biorthogonal wavelet constructions.//Proceedings of the SPIE 2569, Wavelet Applications in Signal and Image Processing Ⅲ. San Diego, CA, United States:SPIE, 2569:68-79. http://spie.org/Publications/Proceedings/Paper/10.1117/12.217619
Tang G, Yang H Z. 2010. Seismic data compression and reconstruction based on Poisson Disk sampling. Chinese Journal of Geophysics (in Chinese), 53(9): 2181-2188. DOI:10.3969/j.issn.0001-5733.2010.09.018
Wang D H, Gao J H. 2014. A new method for random noise attenuation in seismic data based on anisotropic fractional-gradient operators. Journal of Applied Geophysics, 110: 135-143. DOI:10.1016/j.jappgeo.2014.09.011
Wang H Y, Qiu T S. 2015. Unified classification methods for determinate nonstationary signals and random nonstationary signals. Journal on Communications (in Chinese), 36(2): 2015028.
Wang J J, Yuan L, Liu W R, et al. 2016. Dual-tree complex wavelet domain bivariate method for seismic signal random noise attenuation. Chinese Journal of Geophysics (in Chinese), 59(8): 3046-3055. DOI:10.6038/cjg20160827
Wang W, Gao J H, Chen W C, et al. 2012. Random seismic noise suppression via structure-adaptive median filter. Chinese Journal of Geophysics (in Chinese), 55(5): 1732-1741. DOI:10.6038/j.issn.0001-5733.2012.05.030
Wang X K, Gao J H, Chen W C, et al. 2015. The seismic random noise attenuation method based on enhanced bandelet transform. Journal of Applied Geophysics, 116: 146-155. DOI:10.1016/j.jappgeo.2015.03.002
Zhang C, Li Y, Lin H B, et al. 2016. Seismic directional random noise suppression by radial-trace time-frequency peak filtering using the Hurst exponent statistic. Geophysical Prospecting, 64(3): 571-580. DOI:10.1111/gpr.2016.64.issue-3
Zhao X, Li Y, Zhuang G H, et al. 2016. 2-D TFPF based on Contourlet transform for seismic random noise attenuation. Journal of Applied Geophysics, 129: 158-166. DOI:10.1016/j.jappgeo.2016.03.030
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
Zhou J X, Lu W K, He J W, et al. 2015. A data-dependent Fourier filter based on image segmentation for random seismic noise attenuation. Journal of Applied Geophysics, 114: 224-231. DOI:10.1016/j.jappgeo.2015.01.020
Zwartjes P M, Sacchi M D. 2007. Fourier reconstruction of nonuniformly sampled, aliased seismic data. Geophysics, 72(1): V21-V32. DOI:10.1190/1.2399442
曹静杰, 王彦飞, 杨长春. 2012. 地震数据压缩重构的正则化与零范数稀疏最优化方法. 地球物理学报, 55(2): 596-607. DOI:10.6038/j.issn.0001-5733.2012.02.022
黄梅红, 李月. 2016. 基于方向可控滤波的地震勘探随机噪声压制. 地球物理学报, 59(5): 1815-1823. DOI:10.6038/cjg20160524
金雷, 李月, 杨宝俊. 2005. 用时频峰值滤波方法消减地震勘探资料中随机噪声的初步研究. 地球物理学进展, 20(3): 724-728. DOI:10.3969/j.issn.1004-2903.2005.03.024
李庆忠. 1993. 走向精确勘探的道路. 北京:石油工业出版社, 112: 112-114.
林红波, 马海涛, 许丽萍. 2015. 压制空间非平稳地震勘探随机噪声的ROAD径向时频峰值滤波算法. 地球物理学报, 58(7): 2546-2555. DOI:10.6038/cjg20150729
刘洋, Fomel S, 刘财, 等. 2009. 高阶seislet变换及其在随机噪声消除中的应用. 地球物理学报, 52(8): 2142-2151. DOI:10.3969/j.issn.0001-5733.2009.08.024
刘志鹏, 赵伟, 陈小宏, 等. 2012. 局部频率域SVD压制随机噪声方法. 石油地球物理勘探, 47(2): 202-206.
陆基孟. 1991. 地震勘探原理及资料解释. 北京:石油工业出版社: 140-142.
陆文凯. 2011. 基于离散余弦变换的地震随机噪声压制技术. 石油地球物理勘探, 46(2): 202-206.
唐刚, 杨慧珠. 2010. 基于泊松碟采样的地震数据压缩重建. 地球物理学报, 53(9): 2181-2188. DOI:10.3969/j.issn.0001-5733.2010.09.018
王宏禹, 邱天爽. 2015. 非平稳确定性信号与非平稳随机信号统一分类法的探讨. 通信学报, 36(2): 2015028.
汪金菊, 袁力, 刘婉如, 等. 2016. 地震信号随机噪声压制的双树复小波域双变量方法. 地球物理学报, 59(8): 3046-3055. DOI:10.6038/cjg20160827
王伟, 高静怀, 陈文超, 等. 2012. 基于结构自适应中值滤波器的随机噪声衰减方法. 地球物理学报, 55(5): 1732-1741. DOI:10.6038/j.issn.0001-5733.2012.05.030
钟铁, 李月, 杨宝俊, 等. 2017. 陆地地震勘探随机噪声统计特性. 地球物理学报, 60(2): 655-664. DOI:10.6038/cjg20170219