曹靜杰,楊志權(quán),楊 勇,孫秀麗

(1.河北地質(zhì)大學(xué),河北石家莊050031;2.北京市水科學(xué)技術(shù)研究院,北京100048)

1 方法技術(shù)

含隨機(jī)噪聲的地震數(shù)據(jù)可表示為有效信號(hào)和隨機(jī)噪聲之和:

(1)

式中:d表示有效信號(hào),dobs代表觀測(cè)數(shù)據(jù),n表示加性隨機(jī)噪聲。假設(shè)有效信號(hào)d經(jīng)過(guò)某個(gè)變換后有稀疏的表達(dá),則公式(1)可以表示為:

(2)

式中:Ψ為稀疏變換,x是信號(hào)在變換域的稀疏表示?；趚的稀疏性,當(dāng)噪聲的能量已知時(shí),可以通過(guò)求解模型(3)得到x:

(3)

其中σ表示噪聲能量。由于σ的估計(jì)一般比較困難,因此通常求解與模型(3)等價(jià)的稀疏反演模型:

(4)

其中λ>0是正則參數(shù),用來(lái)平衡擬合誤差‖Ψxk-dobs‖2和x的稀疏性。求解模型(4)得到x后,通過(guò)變換Ψ將x變換到時(shí)間-空間域即可實(shí)現(xiàn)去噪的目的?；谠摲囱菽Ｐ腿ピ氲年P(guān)鍵是選擇合適的λ,當(dāng)隨機(jī)噪聲能量較弱時(shí),應(yīng)選擇較小的λ,當(dāng)隨機(jī)噪聲的能量較強(qiáng)時(shí),應(yīng)選擇較大的λ,只有當(dāng)λ選擇合適時(shí),才能達(dá)到解的精度和去噪效果的最佳平衡,當(dāng)λ選擇不合適時(shí),即使再優(yōu)秀的算法也得不到最優(yōu)解。同時(shí),理論上λ和σ存在一一對(duì)應(yīng)關(guān)系,即對(duì)于某個(gè)σ,存在一個(gè)λ使得模型(3)和模型(4)的解相同。由于不同數(shù)據(jù)的噪聲能量不同,因此需要針對(duì)具體的數(shù)據(jù)來(lái)估計(jì)λ。

在求解模型(4)實(shí)現(xiàn)去噪時(shí),一般先給模型(4)賦予較大的正則參數(shù),然后逐漸減小正則參數(shù),以較大正則參數(shù)確定的模型的解作為較小正則參數(shù)確定的模型的初始解,以此類推直到確定出合適的最小正則參數(shù)。當(dāng)采用迭代軟閾值法求解模型(4)時(shí),正則參數(shù)與閾值相對(duì)應(yīng),基于該方法去噪的關(guān)鍵是獲得合適的最小閾值。當(dāng)噪聲能量未知時(shí),一般通過(guò)人工調(diào)節(jié)獲得較為準(zhǔn)確的最小閾值,這會(huì)耗費(fèi)大量的人力和計(jì)算資源。為了實(shí)現(xiàn)自適應(yīng)去噪,固定某個(gè)正則參數(shù)λk或者閾值τk,解模型(4)得到解xk。定義一個(gè)新的參數(shù):

(5)

隨著正則參數(shù)或者閾值的下降,‖xk‖1逐漸增加,而擬合誤差‖Ψxk-dobs‖2逐漸減小,Pk逐漸增大,直到出現(xiàn)最大極值后逐漸減小。多次數(shù)值試驗(yàn)表明,當(dāng)Pk達(dá)到極大值時(shí),所對(duì)應(yīng)的閾值即為一個(gè)較合適的閾值,利用該閾值能夠獲得較好的去噪結(jié)果[12]?；谏鲜霭l(fā)現(xiàn),本文給出一種自適應(yīng)去噪方法的具體實(shí)現(xiàn)步驟。

步驟1:輸入觀測(cè)數(shù)據(jù)dobs,稀疏變換Ψ,定義迭代次數(shù)N和最小閾值τ,令k=0。

步驟2:定義初始解為x0=0,初始?xì)埐顬閞0=dobs,初始閾值為τ0=max‖ΨHdobs‖∞,ΨH為Ψ的共軛轉(zhuǎn)置或者逆變換,令P0=lg(‖x0‖1)·lg(‖r0‖2)。

步驟3:通過(guò)軟閾值運(yùn)算更新變換域的系數(shù),即xk+1=Hτk[ΨH(rk+Ψxk)],其中Hτk(·)為軟閾值運(yùn)算。定義rk+1=dobs-Ψxk+1,令

(6)

步驟4:如果Pk-1Pk+1,則轉(zhuǎn)步驟5;否則通過(guò)指數(shù)法將閾值降低為τk+1,令k=k+1,轉(zhuǎn)步驟3。

步驟5:輸出最終解dfinal=Ψxk。

由于線性閾值下降方法的閾值下降太快,為防止閾值變化劇烈錯(cuò)過(guò)最合適的閾值,本文采用指數(shù)下降方法來(lái)降低閾值,使得迭代后期閾值下降緩慢。設(shè)定在迭代N次內(nèi)下降到最小閾值,實(shí)際則依靠步驟4終止迭代。本文算法針對(duì)每個(gè)τk對(duì)模型(4) 進(jìn)行一次迭代軟閾值運(yùn)算,求出的結(jié)果作為閾值為τk+1時(shí)模型(4)的初始解,這樣會(huì)顯著地減小計(jì)算量,提高計(jì)算效率。為了使得變換域的系數(shù)更加稀疏,本文采用曲波變換[12]為稀疏變換。曲波變換不需要對(duì)數(shù)據(jù)進(jìn)行分塊,比小波變換和傅里葉變換有更加稀疏的表達(dá),并且對(duì)曲波變換進(jìn)行轉(zhuǎn)置即可實(shí)現(xiàn)求逆運(yùn)算,因此是較理想的稀疏變換。關(guān)于曲波變換更加詳細(xì)的內(nèi)容見(jiàn)參考文獻(xiàn)[12]。

2 數(shù)值實(shí)驗(yàn)

2.1 模擬數(shù)據(jù)

利用圖1所示模擬數(shù)據(jù)對(duì)本文方法進(jìn)行了測(cè)試。

首先采用軟閾值方法去噪。令閾值從變換域的最大系數(shù)開(kāi)始,按照指數(shù)方式下降到一個(gè)很小的閾值,由于原始數(shù)據(jù)已知,因此可以求出每個(gè)閾值所對(duì)應(yīng)的去噪結(jié)果的信噪比。令總的閾值下降次數(shù)為70,信噪比隨閾值的變化如圖2a所示,可以看出在第35次迭代(對(duì)應(yīng)的閾值為0.153)時(shí),得到的信噪比最高。模型(4)的L曲線如圖2b所示,此曲線與基于2范數(shù)正則化的L曲線差異很大,因此按照L曲線法不能獲得適合模型(4)的正則參數(shù)。圖3a是本文自適應(yīng)去噪方法獲得的結(jié)果,信噪比為12.666dB,經(jīng)過(guò)36次迭代后算法停止(對(duì)應(yīng)的閾值為0.1543),確定的閾值與圖2a中信噪比最大時(shí)的閾值吻合。由此可見(jiàn),采用本文算法能夠得到一個(gè)合適的閾值,并且自動(dòng)終止迭代,獲得較好的去噪結(jié)果。圖3b是圖1b和圖3a之間的差值剖面,表示濾掉的噪聲。圖3b中除含有隨機(jī)噪聲外,還存在一定程度的有效信號(hào),原因是基于稀疏變換和閾值運(yùn)算的去噪方法將變換域中小于閾值的系數(shù)都作為噪聲,當(dāng)變換域有效信號(hào)的系數(shù)小于閾值時(shí),此方法不僅去掉了噪聲,而且去掉了一部分有效信號(hào)。這是基于稀疏變換的閾值去噪方法自身存在的問(wèn)題,需要研究出更加合適的稀疏表達(dá)方式來(lái)解決。圖4a,圖4b和圖4c分別是模擬數(shù)據(jù)去噪過(guò)程中擬合誤差‖Ψxk-dobs‖2、解的1范數(shù)‖xk‖1和參數(shù)Pk隨迭代次數(shù)的變化情況。隨著閾值的下降,‖Ψxk-dobs‖2逐漸減小,而‖xk‖1逐漸增大,Pk在迭代開(kāi)始階段逐漸增大,在第36次迭代時(shí)開(kāi)始降低,因此取此時(shí)的閾值最為合適。MONTE-FUSCO等[19]提出了H曲線方法,該方法通過(guò)計(jì)算H曲線曲率最大的點(diǎn)確定正則參數(shù),需要計(jì)算出所有的正則參數(shù)確定的反演模型才能得到合適的正則參數(shù),因此計(jì)算量較大。

圖1 模擬數(shù)據(jù)實(shí)驗(yàn)a 原始單炮數(shù)據(jù); b 含隨機(jī)噪聲的單炮數(shù)據(jù)

2.2 實(shí)際地震數(shù)據(jù)

采用實(shí)際地震數(shù)據(jù)驗(yàn)證了本文方法的有效性。圖5a 是我國(guó)西北某地區(qū)含隨機(jī)噪聲的疊后地震數(shù)據(jù),由于噪聲的存在,同相軸模糊,斷層劃分困難。

圖2 模擬數(shù)據(jù)實(shí)驗(yàn)a 信噪比隨迭代次數(shù)的變化; b 模型(4)的L曲線

圖3 模擬數(shù)據(jù)實(shí)驗(yàn)a 本文方法去噪結(jié)果; b 濾掉的噪聲

圖5b是利用本文方法對(duì)其去噪的結(jié)果,設(shè)定的最大迭代次數(shù)為60次,最小閾值為0.004。經(jīng)過(guò)34次迭代后算法自動(dòng)終止,得到的閾值為0.097。由圖5b可見(jiàn),利用本文方法去噪后同相軸變得清晰連貫,能夠清楚地確定斷層位置。圖5c是圖5b與原始含噪聲數(shù)據(jù)(圖5a)的差值剖面,該剖面沒(méi)有明顯的有效信號(hào),說(shuō)明本文去噪方法去掉的基本上為隨機(jī)噪聲。圖6a,圖6b和圖6c分別是該地區(qū)實(shí)際地震數(shù)據(jù)去噪過(guò)程中擬合誤差‖Ψxk-dobs‖2、解的1范數(shù)‖xk‖1和參數(shù)Pk隨迭代次數(shù)的變化情況,可見(jiàn)本文提出的自適應(yīng)去噪方法能夠自動(dòng)獲得閾值,進(jìn)而有效去除隨機(jī)噪聲。

圖4 模擬數(shù)據(jù)實(shí)驗(yàn)a 擬合誤差‖Ψxk-dobs‖2隨迭代次數(shù)的變化; b 解的1范數(shù)‖xk‖1隨迭代次數(shù)的變化; c 參數(shù)Pk隨迭代次數(shù)的變化

圖5 實(shí)際數(shù)據(jù)實(shí)驗(yàn)(西北某地區(qū))a 實(shí)際含噪聲數(shù)據(jù); b 本文方法的去噪結(jié)果; c 本文方法去噪結(jié)果與原始數(shù)據(jù)的差

圖7a為我國(guó)東部某地區(qū)含噪聲地震數(shù)據(jù),同相軸連續(xù)性差,隨機(jī)噪聲強(qiáng)。圖7b是本文方法對(duì)該數(shù)據(jù)去噪的結(jié)果,設(shè)定的最大迭代次數(shù)為60,最小閾值

為0.006。經(jīng)過(guò)29次迭代后算法自動(dòng)終止,確定的閾值為0.1756。由圖7b可見(jiàn),利用本文方法去噪后,一些不連續(xù)和不穩(wěn)定的同相軸可以得到清楚的判識(shí),說(shuō)明本文方法取得了較好的去噪效果。圖7c是圖7a 與圖7b之間的差值剖面,幾乎看不到有效信號(hào),說(shuō)明本文方法能在去噪的同時(shí),不損傷有效信號(hào)。圖8a,圖8b和圖8c分別是該地區(qū)實(shí)際地震數(shù)據(jù)去噪過(guò)程中擬合誤差‖Ψxk-dobs‖2、解的1范數(shù)‖xk‖1和參數(shù)Pk隨迭代次數(shù)的變化情況,在Pk開(kāi)始下降時(shí)迭代停止,再次說(shuō)明本文提出的自適應(yīng)去噪方法能夠自動(dòng)獲得閾值,有效去除隨機(jī)噪聲。

圖6 實(shí)際數(shù)據(jù)實(shí)驗(yàn)(西北某地區(qū))a 擬合誤差‖Ψxk-dobs‖2隨迭代次數(shù)的變化; b 解的1范數(shù)‖xk‖1隨迭代次數(shù)的變化; c 參數(shù)Pk隨迭代次數(shù)的變化

圖7 實(shí)際數(shù)據(jù)實(shí)驗(yàn)(東部某地區(qū))a 原始含噪聲地震數(shù)據(jù); b 本文方法的去噪結(jié)果; c 本文方法去噪結(jié)果與原始數(shù)據(jù)的差

圖8 實(shí)際數(shù)據(jù)實(shí)驗(yàn)(東部某地區(qū))a 擬合誤差‖Ψxk-dobs‖2隨迭代次數(shù)的變化; b 解的1范數(shù)‖xk‖1隨迭代次數(shù)的變化; c 參數(shù)Pk(c)隨迭代次數(shù)的變化

3 結(jié)束語(yǔ)

本文提出一種自適應(yīng)地震數(shù)據(jù)隨機(jī)噪聲消除方法,采用曲波變換作為稀疏變換,利用解的稀疏性和擬合誤差之間的內(nèi)在關(guān)系自動(dòng)獲得合適的正則參數(shù)(閾值),只需要選擇初始閾值和迭代次數(shù)就能夠?qū)崿F(xiàn)迭代去噪過(guò)程的自動(dòng)終止,實(shí)現(xiàn)噪聲的有效去除。由于不需要對(duì)所有固定正則參數(shù)的模型進(jìn)行求解就能夠獲得較好的去噪結(jié)果,因此計(jì)算成本大大減少。

需要注意的是,為了防止初始迭代時(shí)解的不穩(wěn)定,前幾次迭代結(jié)果不能作為判斷終止迭代的標(biāo)準(zhǔn)。另外,若采用三維曲波變換為稀疏變換,則本文算法可用于三維地震數(shù)據(jù)去噪。

致謝:衷心感謝長(zhǎng)安大學(xué)包乾宗老師和中國(guó)石油大學(xué)(北京)袁三一老師對(duì)本文研究的指導(dǎo)和幫助。

[1] 王華忠,馮波,王雄文,等.壓縮感知及其在地震勘探中的應(yīng)用[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] 康治,于承業(yè),賈臥,等.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] 崔樹(shù)果,朱凌燕,王建花.f-x域Cadzow技術(shù)分塊壓制隨機(jī)噪聲及其應(yīng)用[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] 王維強(qiáng),楊國(guó)權(quán).基于EMD與ICA的地震信號(hào)去噪技術(shù)研究[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] 孫成禹,邵婕,藍(lán)陽(yáng),等.基于獨(dú)立分量分析基的地震隨機(jī)噪聲壓制[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] 陳香朋,曹思遠(yuǎn).第二代小波變換及其在地震信號(hào)去噪中的應(yīng)用[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] 曹靜杰,王彥飛,楊長(zhǎng)春.地震數(shù)據(jù)壓縮重構(gòu)的正則化與零范數(shù)稀疏最優(yōu)化方法[J].地球物理學(xué)報(bào),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] 鞏向博,韓立國(guó),王恩利,等.壓制噪聲的高分辨率Radon變換法[J].吉林大學(xué)學(xué)報(bào)(地球科學(xué)版),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ìn)凸集投影方法的地震數(shù)據(jù)同時(shí)插值和去噪[J].地球物理學(xué)報(bào),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