一种基于曲波变换的自适应地震随机噪声消除方法
曹静杰;杨志权;杨勇;孙秀丽
【摘 要】The conventional sparse inversion-based random noise elimination utilizes a thresholding operation to conduct denoising, on the basis that seismic signals are sparsely expressed in a transform domain.This would produce effective denoising when thresh-old values match noise energy.However,owing to variety in the noise energy of different data,the reasonable threshold is usually obtained by manual adjustment,which is time- and labor-consuming.This paper proposes an adaptive random noise elimination method that does not rely on noise energy estimation.The method uses a Curvelet transform as the sparse transform and chooses a proper threshold value through the inner relationship between solution sparsity and fitting error,thus terminating the iteration au-tomatically.Testing on both synthetic data and field data demonstrate that the proposed method can eliminate random noise while preserving effective signal.%基于稀疏反演的随机噪声消除方法需要估计一个与噪声能量相匹配的阈值才能获得可靠的去噪结果.由于不同数据的噪声能量不同,因此通常采用人
工调节的方法获得合理的阈值估计,这会耗费大量的计算资源和人力成本.为此提出一种自适应的随机噪声消除方法,以曲波变换为稀疏变换,通过迭代过程中解的稀疏性与拟合误差之问的内在关系确定合适的阈值,并且自动终止迭代,因而不依赖于对噪声能量的估计就能实现对噪声的消除.利用理论模型数据及两个地区实际地震数据验证了方法的有效性.
【期刊名称】《石油物探》
【年(卷),期】正则化参数的自适应估计2018(057)001
【总页数】8页(P72-78,121)
【关键词】阈值估计;去噪;自适应;稀疏反演;曲波变换;稀疏性;拟合误差;正则参数
【作 者】曹静杰;杨志权;杨勇;孙秀丽
【作者单位】河北地质大学,河北石家庄 050031;河北地质大学,河北石家庄 050031;北京市水科学技术研究院,北京 100048;河北地质大学,河北石家庄 050031
【正文语种】中 文
【中图分类】P631
野外地震数据不可避免地含有各种噪声,尤其是随机噪声,会影响偏移、反褶积、插值、多次波消除等处理的结果[1],消除随机噪声、获得高信噪比的地震数据,对于后期地震数据处理、解释和反演非常重要。随机噪声消除方法有f-x域去噪方法[2]、局部奇异值分解法[3]、Karhunen-Loeve(K-L)变换法[4]、Cadzow滤波法[5]、经验模态分解法[6]、基于独立分量分析的去噪方法[7]、基于稀疏反演的去噪方法[8-9]等等。近年来字典学习方法[10]也被用于地震数据去噪。这些方法大都假设地震记录的相邻道含有随机噪声,并且仅适合线性同相轴记录,基于这些方法去噪需要对数据进行分块处理。基于稀疏反演的方法利用解的某种先验信息建立反演模型,通过求解该模型实现去噪的目的。先验信息有全变差极小[5]、地震数据在某个变换域的稀疏性约束[8-9]等。常用的变换有小波变换[11]、曲波变换[12]、Radon变换[13]、物理小波框架[14]、S变换[15]等。当采用曲波变换[9]等多尺度变换时,基于稀疏反演的方法可以解决非线性同相轴记录的去噪问题。稀疏反演模型确定后,还需估计一个合适的正则参数(阈值参数、噪声能量)才能得到可靠的解[16-17]。而估计一个合适的正则参数往往需要进行多次人工调试,这会耗费很多计算成本和时间。目前反演领域中获得正则参数(阈值)的方法主要有广义交叉检验准则和L曲线法。广义交叉检验准则定义函数其中δ为变换域的阈值,ωδ为
经过阈值处理后的变换域的系数,ω为阈值运算之前的变换域的系数,M为变换域系数的个数,M0为经过阈值处理后被置为0的变换域的系数个数。调整阈值δ的大小使得GCV函数取极小时,经过阈值去噪后的信号均方根误差逼近最小,从而获得合适的正则参数。L曲线方法以目标函数和相应的约束项函数的对数为坐标得到一条参数曲线,该曲线一般是L形状,在曲率最大点处对应的正则参数一般被当作最优正则参数。但是上述两种方法都是针对2范数为正则项的反演模型提出的,并且在某些情况下会失效[18]。本文基于稀疏变换的去噪原理,针对地震数据去噪问题的1范数正则化模型,提出一种适合该模型的自适应随机噪声消除方法,利用理论模型数据和实际地震数据对方法进行了测试。
1 方法技术
含随机噪声的地震数据可表示为有效信号和随机噪声之和:
(1)
式中:d表示有效信号,dobs代表观测数据,n表示加性随机噪声。假设有效信号d经过某个变换后有稀疏的表达,则公式(1)可以表示为:
(2)
式中:Ψ为稀疏变换,x是信号在变换域的稀疏表示。基于x的稀疏性,当噪声的能量已知时,可以通过求解模型(3)得到x:
(3)
其中σ表示噪声能量。由于σ的估计一般比较困难,因此通常求解与模型(3)等价的稀疏反演模型:
(4)
其中λ>0是正则参数,用来平衡拟合误差‖Ψxk-dobs‖2和x的稀疏性。求解模型(4)得到x后,通过变换Ψ将x变换到时间-空间域即可实现去噪的目的。基于该反演模型去噪的关键是选择合适的λ,当随机噪声能量较弱时,应选择较小的λ,当随机噪声的能量较强时,应选择较大的λ,只有当λ选择合适时,才能达到解的精度和去噪效果的最佳平衡,当λ选择不合适时,即使再优秀的算法也得不到最优解。同时,理论上λ和σ存在一一对应关系,即对于某个σ,存在一个λ使得模型(3)和模型(4)的解相同。由于不同数据的噪声能量不同,因此需要针对具体的数据来估计λ。
在求解模型(4)实现去噪时,一般先给模型(4)赋予较大的正则参数,然后逐渐减小正则参数,以较大正则参数确定的模型的解作为较小正则参数确定的模型的初始解,以此类推直到确定出合适的最小正则参数。当采用迭代软阈值法求解模型(4)时,正则参数与阈值相对应,基于该方法去噪的关键是获得合适的最小阈值。当噪声能量未知时,一般通过人工调节获得较为准确的最小阈值,这会耗费大量的人力和计算资源。为了实现自适应去噪,固定某个正则参数λk或者阈值τk,解模型(4)得到解xk。定义一个新的参数:
(5)
随着正则参数或者阈值的下降,‖xk‖1逐渐增加,而拟合误差‖Ψxk-dobs‖2逐渐减小,Pk逐渐增大,直到出现最大极值后逐渐减小。多次数值试验表明,当Pk达到极大值时,所对应的阈值即为一个较合适的阈值,利用该阈值能够获得较好的去噪结果[12]。基于上述发现,本文给出一种自适应去噪方法的具体实现步骤。
步骤1:输入观测数据dobs,稀疏变换Ψ,定义迭代次数N和最小阈值τ,令k=0。
步骤2:定义初始解为x0=0,初始残差为r0=dobs,初始阈值为τ0=max‖ΨHdobs‖∞,ΨH为Ψ的共轭转置或者逆变换,令P0=lg(‖x0‖1)·lg(‖r0‖2)。
步骤3:通过软阈值运算更新变换域的系数,即xk+1=Hτk[ΨH(rk+Ψxk)],其中Hτk(·)为软阈值运算。定义rk+1=dobs-Ψxk+1,令
(6)
步骤4:如果Pk-1<Pk>Pk+1,则转步骤5;否则通过指数法将阈值降低为τk+1,令k=k+1,转步骤3。
步骤5:输出最终解dfinal=Ψxk。
由于线性阈值下降方法的阈值下降太快,为防止阈值变化剧烈错过最合适的阈值,本文采用指数下降方法来降低阈值,使得迭代后期阈值下降缓慢。设定在迭代N次内下降到最小阈值,实际则依靠步骤4终止迭代。本文算法针对每个τk对模型(4) 进行一次迭代软阈值运算,求出的结果作为阈值为τk+1时模型(4)的初始解,这样会显著地减小计算量,提高计算效率。为了使得变换域的系数更加稀疏,本文采用曲波变换[12]为稀疏变换。曲波变换不需要对数据进行分块,比小波变换和傅里叶变换有更加稀疏的表达,并且对曲波变换进行转置即可实现求逆运算,因此是较理想的稀疏变换。关于曲波变换更加详细的内容见参考文献[12]。
2 数值实验
2.1 模拟数据
利用图1所示模拟数据对本文方法进行了测试。
图1a 是一个模拟地震单炮数据,共128道,时间采样点数为256,时间采样率为4ms;图1b为图1a加入随机噪声的结果,信噪比为6.7325dB。信噪比的定义为RSN=10lg(‖x/‖x-xnosie),其中xnoise表示含有噪声的数据,x表示不含噪声的数据。
首先采用软阈值方法去噪。令阈值从变换域的最大系数开始,按照指数方式下降到一个很小的阈值,由于原始数据已知,因此可以求出每个阈值所对应的去噪结果的信噪比。令总的阈值下降次数为70,信噪比随阈值的变化如图2a所示,可以看出在第35次迭代(对应的阈值为0.153)时,得到的信噪比最高。模型(4)的L曲线如图2b所示,此曲线与基于2范数正则化的L曲线差异很大,因此按照L曲线法不能获得适合模型(4)的正则参数。图3a是本文自适应去噪方法获得的结果,信噪比为12.666dB,经过36次迭代后算法停止(对应的阈值为0.1543),确定的阈值与图2a中信噪比最大时的阈值吻合。由此可见,采用本文算法能够得到一个合适的阈值,并且自动终止
迭代,获得较好的去噪结果。图3b是图1b和图3a之间的差值剖面,表示滤掉的噪声。图3b中除含有随机噪声外,还存在一定程度的有效信号,原因是基于稀疏变换和阈值运算的去噪方法将变换域中小于阈值的系数都作为噪声,当变换域有效信号的系数小于阈值时,此方法不仅去掉了噪声,而且去掉了一部分有效信号。这是基于稀疏变换的阈值去噪方法自身存在的问题,需要研究出更加合适的稀疏表达方式来解决。图4a,图4b和图4c分别是模拟数据去噪过程中拟合误差‖Ψxk-dobs‖2、解的1范数‖xk‖1和参数Pk随迭代次数的变化情况。随着阈值的下降,‖Ψxk-dobs‖2逐渐减小,而‖xk‖1逐渐增大,Pk在迭代开始阶段逐渐增大,在第36次迭代时开始降低,因此取此时的阈值最为合适。MONTE-FUSCO等[19]提出了H曲线方法,该方法通过计算H曲线曲率最大的点确定正则参数,需要计算出所有的正则参数确定的反演模型才能得到合适的正则参数,因此计算量较大。

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。