一种基于曲波变换的自适应地震随机噪声消除方法
2018-02-27曹静杰杨志权孙秀丽
曹静杰,杨志权,杨 勇,孙秀丽
(1.河北地质大学,河北石家庄050031;2.北京市水科学技术研究院,北京100048)
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
步骤5:输出最终解dfinal=Ψxk。
由于线性阈值下降方法的阈值下降太快,为防止阈值变化剧烈错过最合适的阈值,本文采用指数下降方法来降低阈值,使得迭代后期阈值下降缓慢。设定在迭代N次内下降到最小阈值,实际则依靠步骤4终止迭代。本文算法针对每个τk对模型(4) 进行一次迭代软阈值运算,求出的结果作为阈值为τk+1时模型(4)的初始解,这样会显著地减小计算量,提高计算效率。为了使得变换域的系数更加稀疏,本文采用曲波变换[12]为稀疏变换。曲波变换不需要对数据进行分块,比小波变换和傅里叶变换有更加稀疏的表达,并且对曲波变换进行转置即可实现求逆运算,因此是较理想的稀疏变换。关于曲波变换更加详细的内容见参考文献[12]。
2 数值实验
2.1 模拟数据
利用图1所示模拟数据对本文方法进行了测试。
首先采用软阈值方法去噪。令阈值从变换域的最大系数开始,按照指数方式下降到一个很小的阈值,由于原始数据已知,因此可以求出每个阈值所对应的去噪结果的信噪比。令总的阈值下降次数为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曲线曲率最大的点确定正则参数,需要计算出所有的正则参数确定的反演模型才能得到合适的正则参数,因此计算量较大。
图1 模拟数据实验a 原始单炮数据; b 含随机噪声的单炮数据
2.2 实际地震数据
采用实际地震数据验证了本文方法的有效性。图5a 是我国西北某地区含随机噪声的叠后地震数据,由于噪声的存在,同相轴模糊,断层划分困难。
图2 模拟数据实验a 信噪比随迭代次数的变化; b 模型(4)的L曲线
图3 模拟数据实验a 本文方法去噪结果; b 滤掉的噪声
图5b是利用本文方法对其去噪的结果,设定的最大迭代次数为60次,最小阈值为0.004。经过34次迭代后算法自动终止,得到的阈值为0.097。由图5b可见,利用本文方法去噪后同相轴变得清晰连贯,能够清楚地确定断层位置。图5c是图5b与原始含噪声数据(图5a)的差值剖面,该剖面没有明显的有效信号,说明本文去噪方法去掉的基本上为随机噪声。图6a,图6b和图6c分别是该地区实际地震数据去噪过程中拟合误差‖Ψxk-dobs‖2、解的1范数‖xk‖1和参数Pk随迭代次数的变化情况,可见本文提出的自适应去噪方法能够自动获得阈值,进而有效去除随机噪声。
图4 模拟数据实验a 拟合误差‖Ψxk-dobs‖2随迭代次数的变化; b 解的1范数‖xk‖1随迭代次数的变化; c 参数Pk随迭代次数的变化
图5 实际数据实验(西北某地区)a 实际含噪声数据; b 本文方法的去噪结果; c 本文方法去噪结果与原始数据的差
图7a为我国东部某地区含噪声地震数据,同相轴连续性差,随机噪声强。图7b是本文方法对该数据去噪的结果,设定的最大迭代次数为60,最小阈值
为0.006。经过29次迭代后算法自动终止,确定的阈值为0.1756。由图7b可见,利用本文方法去噪后,一些不连续和不稳定的同相轴可以得到清楚的判识,说明本文方法取得了较好的去噪效果。图7c是图7a 与图7b之间的差值剖面,几乎看不到有效信号,说明本文方法能在去噪的同时,不损伤有效信号。图8a,图8b和图8c分别是该地区实际地震数据去噪过程中拟合误差‖Ψxk-dobs‖2、解的1范数‖xk‖1和参数Pk随迭代次数的变化情况,在Pk开始下降时迭代停止,再次说明本文提出的自适应去噪方法能够自动获得阈值,有效去除随机噪声。
图6 实际数据实验(西北某地区)a 拟合误差‖Ψxk-dobs‖2随迭代次数的变化; b 解的1范数‖xk‖1随迭代次数的变化; c 参数Pk随迭代次数的变化
图7 实际数据实验(东部某地区)a 原始含噪声地震数据; b 本文方法的去噪结果; c 本文方法去噪结果与原始数据的差
图8 实际数据实验(东部某地区)a 拟合误差‖Ψxk-dobs‖2随迭代次数的变化; b 解的1范数‖xk‖1随迭代次数的变化; c 参数Pk(c)随迭代次数的变化
3 结束语
本文提出一种自适应地震数据随机噪声消除方法,采用曲波变换作为稀疏变换,利用解的稀疏性和拟合误差之间的内在关系自动获得合适的正则参数(阈值),只需要选择初始阈值和迭代次数就能够实现迭代去噪过程的自动终止,实现噪声的有效去除。由于不需要对所有固定正则参数的模型进行求解就能够获得较好的去噪结果,因此计算成本大大减少。
需要注意的是,为了防止初始迭代时解的不稳定,前几次迭代结果不能作为判断终止迭代的标准。另外,若采用三维曲波变换为稀疏变换,则本文算法可用于三维地震数据去噪。
致谢:衷心感谢长安大学包乾宗老师和中国石油大学(北京)袁三一老师对本文研究的指导和帮助。
[1] 王华忠,冯波,王雄文,等.压缩感知及其在地震勘探中的应用[J].石油物探,2016,55(4):467-474
WANG H Z,FENG B,WANG X W,et al.Compressed sensing and its application in seismic exploration[J].Geophysical Prospecting for Petroleum,2016,55(4):467-474
[2] 康治,于承业,贾卧,等.f-x域去噪方法研究[J].石油地球物理勘探,2003,38(2):136-138
KANG Z,YU C Y,JIA W,et al.A study on noise-suppression method inf-xdomain[J].Oil Geophysical Prospecting,2003,38(2):136-138
[3] 谢凤兰,赵改善.利用SVD法重建地震剖面[J].石油物探,1990,29(4):33-40
XIE F L,ZHAO G S.Seismic section reconstruction by SVD technique[J].Geophysical Prospecting for Petroleum,1990,29(4):33-40
[4] JONES I F,LEVY S.Signal-to-noise ratio enhancement in multichannel seismic data via the Karhunen-Loeve transform[J].Geophysical Prospecting,1987,35(1):12-32
[5] 崔树果,朱凌燕,王建花.f-x域Cadzow技术分块压制随机噪声及其应用[J].石油物探,2012,51(1):43-50
CUI S G,ZHU L Y,WANG J H.f-xdomain Cadzow blocking technique for random noise suppression and its application[J].Geophysical Prospecting for Petroleum,2012,51(1):43-50
[6] 王维强,杨国权.基于EMD与ICA的地震信号去噪技术研究[J].石油物探,2012,51(1):19-29
WANG W Q,YANG G Q.The study of seismic denoising based on EMD and ICA[J].Geophysical Prospecting for Petroleum,2012,51(1):19-29
[7] 孙成禹,邵婕,蓝阳,等.基于独立分量分析基的地震随机噪声压制[J].石油物探,2016,55(2):196-204
SUN C Y,SHAO J,LAN Y,et al.Seismic random noise suppression based on independent component a-
nalysis basis functions[J].Geophysical Prospecting for Petroleum,2016,55(2):196-204
[8] YUAN S Y,WANG S X,LI G F.Random noise reduction using Bayesian inversion[J].Journal of Geophysics and Engineering,2012,9(1):60-68
[9] CAO J J,ZHAO J T,HU Z Y.3D seismic denoising based on a low-redundancy Curvelet transform[J].Journal of Geophysics and Engineering,2015,12(4):566-576
[10] BECKOUCHE S,Ma J W.Simultaneous dictionary learning and denoising for seismic data[J].Geophysics,2014,79(3):A27-A31
[11] 陈香朋,曹思远.第二代小波变换及其在地震信号去噪中的应用[J].石油物探,2004,43(6):547-550
CHEN X P,CAO S Y.The second generation wavelet transform and its application in seismic signal denoising[J].Geophysical Prospecting for Petroleum,2004,43(6):547-550
[12] 曹静杰,王彦飞,杨长春.地震数据压缩重构的正则化与零范数稀疏最优化方法[J].地球物理学报,2012,55(2):596-607
CAO J J,WANG Y F,YANG C C.Seismic data restoration based on compresive sensing using the regularization and zero-norm sparse optimization[J].Chinese Journal of Geophysics,2012,55(2):596-607
[13] 巩向博,韩立国,王恩利,等.压制噪声的高分辨率Radon变换法[J].吉林大学学报(地球科学版),2009,39(1):152-157
GONG X B,HAN L G,WANG E L,et al.Denoising via high resolution Radon transform[J].Journal of Jilin University:Earth Science Edition,2009,39(1):152-157
[14] ZHANG R F,ULRYCH T J.Physical wavelet frame denoising[J].Geophysics,2003,68(1):225-231
[15] PAROLAI S.Denoising of seismograms using the S transform[J].Bulletin of the Seismology Society of America,2009,99(1):226-234
[16] DAUBECHIES I,DEFRISE M,MOL C.An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[J].Communications on Pure and Applied Mathematics,2004,57(11):1413-1457
[17] 曹静杰,王本锋.基于一种改进凸集投影方法的地震数据同时插值和去噪[J].地球物理学报,2015,58(8):2935-2947
CAO J J,WANG B F.An improved projection onto convex set method for simultaneous interpolation and denoising[J].Chinese Journal of Geophysics,2015,58(8):2935-2947
[18] YAGOLA A G,LEONOV A S,TITARENKO V N.Data errors and an error estimation for ill-posed problems[J].Inverse Problems in Science and Engineering,2002,10(2):117-129
[19] MONTEFUSCO L B,PAPI S.A parameter selection method for wavelet shrinkage denoising[J].BIT Numerical Mathematics,2003,43(3):611-626
