樂友喜, 楊杰飛, 陳藝都, 吳佳偉, 楊 濤

(1.中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院,山東青島 266580; 2.黃河勘測規(guī)劃設(shè)計(jì)研究院有限公司,河南鄭州 450003)

隨著勘探技術(shù)的快速發(fā)展,勘探目標(biāo)轉(zhuǎn)向了復(fù)雜的隱蔽油氣藏,對于深部的復(fù)雜地質(zhì)構(gòu)造,地震波在介質(zhì)中經(jīng)過較長距離的傳播、散射以及介質(zhì)非彈性造成的能量衰減,使得深層反射信號(hào)相對較弱,弱信號(hào)常淹沒在背景噪聲中[1-2]。常用的地震弱信號(hào)識(shí)別方法包括基于奇異值分解(SVD)的方法[3]、基于非線性動(dòng)力系統(tǒng)的混沌理論[4]、基于經(jīng)驗(yàn)?zāi)Ｊ椒纸?EMD)的Hilbert-Huang變換[5]、基于獨(dú)立分量分析(ICA)的盲源分離(BSS)技術(shù)[6]以及基于高階統(tǒng)計(jì)量識(shí)別[7]等。Aharon[8]提出了KSVD(K均值奇異值分解)算法,韓文功等[9]提出的SVD分解方法可以較好地進(jìn)行弱信號(hào)的識(shí)別和去噪。但KSVD算法需要多次進(jìn)行SVD分解,效率較慢。Rubinstein等[10]在KSVD的基礎(chǔ)上提出了AKSVD(approximate KSVD)算法,提高了效率,但其效果不甚理想。Irofti 等[11]將正則化引入AKSVD的字典學(xué)習(xí),Dumitrescud等[12]在KSVD基礎(chǔ)上引入了正則化,完善了RAKSVD(正則化近似KSVD)方法,該方法既提高了計(jì)算效率,又可使分解效果基本不受影響。但其中的正則化參數(shù)都是人為給定的固定值,顯然在SVD迭代分解時(shí),誤差項(xiàng)縮小,而系數(shù)誤差項(xiàng)擴(kuò)大,這是不準(zhǔn)確的。筆者對正則化參數(shù)對迭代結(jié)果的影響進(jìn)行研究,通過L曲線進(jìn)行正則化參數(shù)的自動(dòng)優(yōu)選,形成基于正則化參數(shù)自動(dòng)優(yōu)選的RAKSVD去噪方法。

1 RAKSVD方法基本原理

1.1 KSVD算法

信號(hào)可以通過字典進(jìn)行稀疏表示,但初始化固定字典并不能保證已知信號(hào)的稀疏程度,從而不能準(zhǔn)確表達(dá)信號(hào),為了使稀疏分解更準(zhǔn)確,通過KSVD構(gòu)建字典,使已知低維信號(hào)都可以進(jìn)行較為準(zhǔn)確地稀疏表示。KSVD算法是通過訓(xùn)練信號(hào)對初始化字典進(jìn)行逐列更新,用來尋找一個(gè)最佳字典的過程。假定訓(xùn)練的信號(hào)為Y∈Rm×n,其中m為信號(hào)的行數(shù);n為信號(hào)的列數(shù)。每列對應(yīng)一個(gè)訓(xùn)練樣本,對于給定的冗余字典D∈Rm×p,其中p為字典的列數(shù)且m

(1)

先固定字典D,用正交匹配追蹤(OMP)算法對每一個(gè)樣本做稀疏表示,獲取稀疏分解系數(shù)矩陣。然后進(jìn)行字典更新,KSVD對字典原子逐列進(jìn)行。

(2)計(jì)算殘差:

(2)

式中,Ek代表剩余殘差。

(3)

式中,U為左奇異矩陣;V為右奇異矩陣;Λ為奇異值合集;r為非零奇異值的個(gè)數(shù);ui為左奇異矩陣的第i列;σi為第i個(gè)奇異值;νi為右奇異矩陣的第i列。

(4)原子和稀疏系數(shù)應(yīng)更新為

xk=σ1ν1,dk=u1.

(4)

1.2 AKSVD算法

由于KSVD算法分解效率較慢,而且每次更新一列都要SVD分解,導(dǎo)致算法較慢,實(shí)際應(yīng)用不大。對此,Rubinstein等[10]提出了一種近似KSVD(AKSVD)算法。AKSVD算法的原子和相關(guān)系數(shù)分別按如下公式更新:

(5)

(6)

1.3 RAKSVD算法

AKSVD算法相對于KSVD算法,明顯提高了計(jì)算效率,但由于近似使SVD分解結(jié)果產(chǎn)生了一定的誤差,因而影響了分解效果。針對這一問題,Irofti等[11]提出了正則化近似KSVD(RAKSVD)算法,通過引入正則化參數(shù)對學(xué)習(xí)模型進(jìn)行優(yōu)化,建立嶺估計(jì)模型,解決稀疏分解中的病態(tài)問題,從而使SVD分解結(jié)果的誤差得到很大程度上的控制,提高了稀疏分解效果。

地震信號(hào)的稀疏分解是一個(gè)典型的病態(tài)問題,而正則化方法一直是解決病態(tài)問題的一種有效方法。正則化字典模型問題可轉(zhuǎn)化為

(7)

式中,μ>0為正則化參數(shù);第一項(xiàng)用來衡量樣本的預(yù)測值和真實(shí)值之間的誤差,第二項(xiàng)通過對參數(shù)加一個(gè)系數(shù)來調(diào)整其權(quán)值。

(8)

根據(jù)式(8),從KSVD算法的正則化[11]可以得出,正則化約束后的更新原子方法和原算法一致,而正則化參數(shù)是應(yīng)用到稀疏系數(shù)的更新,使系數(shù)變?yōu)?/(1+μ)倍,通過給稀疏系數(shù)一個(gè)權(quán)值,每次通過這個(gè)權(quán)值影響后續(xù)的原子更新,最終完成字典的更新。所以RAKSVD的原子和稀疏系數(shù)更新為

(9)

(10)

2 RAKSVD算法的正則化參數(shù)自動(dòng)優(yōu)選

正則化參數(shù)μ對正則化效果有較大的影響,而在前人的研究文獻(xiàn)中,μ一般取為常數(shù),這樣雖對模型進(jìn)行了正則化處理,卻得不到最優(yōu)結(jié)果。因此正則化參數(shù)的自動(dòng)優(yōu)選是極其重要的。

2.1 方法原理

為了求得L-corner的參數(shù)值,通常引入曲率,曲率最大時(shí)則對應(yīng)此處的值,曲率表達(dá)式一般為

(11)

其中

式中,με為ε次迭代所對應(yīng)的正則化參數(shù)變量;ε為迭代次數(shù);δ為相應(yīng)的曲率函數(shù);η′、ρ′為一階導(dǎo)數(shù);η″、ρ″為二階導(dǎo)數(shù)。

對于RAKSVD的正則化參數(shù)選取,理想的效果應(yīng)為每次循環(huán)都會(huì)利用L曲線法求取μ值,但L曲線計(jì)算需要先對字典D進(jìn)行奇異值分解,影響分解效率,所以需要進(jìn)一步改進(jìn)算法,將整個(gè)過程插入若干個(gè)控制節(jié)點(diǎn),對于RAKSVD每一次更新字典時(shí)的正則化參數(shù)估算值經(jīng)數(shù)值擬合得到:

(12)

式中,u為分段優(yōu)選正則化參數(shù)的控制節(jié)點(diǎn)序號(hào);μu、Su為調(diào)節(jié)系數(shù)(μu、Su>0);μ1為通過L曲線求取的第一次迭代的μ值;S1為給定的初始值。

調(diào)節(jié)系數(shù)Su和μu通過梯度下降法求得,對應(yīng)的損失函數(shù)為

(13)

(14)

(15)

2.2 稀疏分解效果對比

設(shè)計(jì)弱信號(hào)的褶積模型(圖1(a)),該模型共有3層,上層為強(qiáng)軸,中間層為傾斜層,由上到下逐漸變?nèi)?并與上下兩層斜交,下層為弱信號(hào),共25道,每道300個(gè)采樣點(diǎn),其中雷克子波主頻為25 Hz,模型中存在嚴(yán)重高斯白噪聲(圖1(b)),峰值信噪比為1.2 dB,用來進(jìn)行字典訓(xùn)練。

圖1 褶積噪音模型

本次實(shí)驗(yàn)硬件環(huán)境(i5-8250U, 1.6 GHz,8 GB)算法均在MATLAB2016a中實(shí)現(xiàn),用RMSE代表均方根誤差,其中字典大小為64×64,將地震體分成8×8的塊,對加噪模型進(jìn)行處理。迭代次數(shù)為30,為了求取準(zhǔn)確,循環(huán)10次求平均值,圖2是經(jīng)優(yōu)選得到的每次迭代所對應(yīng)的正則化參數(shù)取值結(jié)果,其中μ為無量綱參數(shù),圖中5個(gè)圓點(diǎn)代表了根據(jù)L曲線優(yōu)化得到的控制節(jié)點(diǎn)處的μ值,從第二個(gè)控制點(diǎn)開始對調(diào)節(jié)系數(shù)μu、Su進(jìn)行修正,直到最大迭代次數(shù)結(jié)束。不同正則化參數(shù)所對應(yīng)的均方根誤差(RMSE)隨迭代次數(shù)的變化關(guān)系如圖3所示。

圖2 正則化參數(shù)與迭代次數(shù)的關(guān)系

由圖3可以看出,當(dāng)訓(xùn)練字典在25次迭代后RMSE基本趨于穩(wěn)定,分別對RAKSVD的正則化參數(shù)μ在0～0.3之間取固定值進(jìn)行了測試:當(dāng)μ為0時(shí),訓(xùn)練誤差最差,RMSE為1.440,隨著μ增大,RMSE逐漸下降;當(dāng)μ取0.2時(shí),RMSE降為1.403,第一次迭代通過L曲線計(jì)算得到的μ為0.22。 但若μ取固定值0.22時(shí),其分解效果并沒有比0.2好,RMSE增大到1.405,且當(dāng)μ為0.3時(shí)誤差變得更大,說明整個(gè)迭代過程中較合適的μ在0.1～0.22,μ過大或過小分解效果都不好,當(dāng)使用正則化參數(shù)計(jì)算公式(12)每次迭代更新μ進(jìn)行字典訓(xùn)練時(shí),就形成了正則化參數(shù)自動(dòng)優(yōu)選的RAKSVD方法(μ自動(dòng)優(yōu)選),可以避免μ取固定值不準(zhǔn)確所帶來的誤差,從而實(shí)現(xiàn)了每次迭代自動(dòng)選取μ值的目的,此時(shí)RMSE為1.344,達(dá)到了最優(yōu)稀疏分解效果。

圖3 不同正則化參數(shù)取值的誤差對比

3 正則化參數(shù)自動(dòng)優(yōu)選的RAKSVD去噪方法

一般來講,地震資料數(shù)據(jù)量較大,首先將其分割成若干個(gè)較小且大小相等的地震子塊,在對實(shí)際地震資料去噪處理時(shí),從整體稀疏分解問題上考慮,在式(7)中加入一個(gè)約束項(xiàng),目標(biāo)函數(shù)轉(zhuǎn)變?yōu)槿缦聝?yōu)化問題:

(16)

式中,λ為拉格朗日參數(shù);N為重構(gòu)信號(hào);Y為輸入信號(hào),即含噪的地震剖面;D為需要更新的字典;i、j為分塊矩陣的編號(hào);Lij為選取Y的位置矩陣;Xij為相應(yīng)的分塊系數(shù)矩陣。

式(16)中第一部分是一個(gè)約束項(xiàng),促使含噪信號(hào)與稀疏去噪后的信號(hào)能夠盡可能地趨近;第二部分是正則化,目的在于解決模型稀疏分解的病態(tài)問題;第三部分確保每個(gè)地震子塊經(jīng)稀疏分解后都有很好的稀疏性。

實(shí)現(xiàn)步驟如下:

(1)假設(shè)初始值。N=Y,D0設(shè)定為冗余DCT字典;

(2)迭代循環(huán)。

a) 稀疏編碼階段通過正交匹配追蹤(OMP)算法求取每個(gè)子塊在字典下的稀疏分解系數(shù),然后問題轉(zhuǎn)化為

b)利用L曲線求出μ1,對該次迭代進(jìn)行正則化字典與系數(shù)更新。

c)按照公式(12)更新正則化系數(shù),更換字典D中原子,求出Xij。

d)達(dá)到最大迭代次數(shù)停止循環(huán)。

(3)去噪后地震數(shù)據(jù)的重構(gòu)。

前面已經(jīng)通過正則項(xiàng)更新了字典,求解了正則化參數(shù),故可以固定正則項(xiàng),使問題簡化為

在已經(jīng)更新完畢的學(xué)習(xí)字典下,應(yīng)用稀疏算法對每個(gè)子塊求解相應(yīng)的表示系數(shù)Xij,利用公式進(jìn)行地震數(shù)據(jù)重構(gòu),得到去噪后的結(jié)果。公式為

式中,I為單位矩陣。

4 弱信號(hào)模型測試

圖1(b)中存在較強(qiáng)的高斯白噪聲,強(qiáng)軸雖能明顯看到,但弱軸淹沒在噪聲中不易識(shí)別,圖4為不同的正則化參數(shù)對去噪效果圖。同時(shí),將本方法與KSVD算法和AKSVD算法進(jìn)行對比。表1為不同方法去噪效果和運(yùn)行時(shí)間,RAKSVD方法指的是參數(shù)自動(dòng)優(yōu)選的RAKSVD方法。

從圖4中可看到:對于低信噪比弱地震信號(hào),當(dāng)μ為0時(shí)(此時(shí)為AKSVD算法),去噪結(jié)果雖能保留弱信號(hào),但去噪效果較差,殘余很多毛刺,此時(shí)峰值信噪比為12.03 dB;如果RAKSVD算法中正則化參數(shù)取固定值,當(dāng)μ取0.2時(shí),去噪效果明顯提高,此時(shí)峰值信噪比為12.76 dB;如果根據(jù)正則化參數(shù)自動(dòng)優(yōu)選取值時(shí),隨著每次迭代,其值是相應(yīng)調(diào)整的,此時(shí)峰值信噪比為13.49 dB,能夠達(dá)到最佳去噪效果,說明正則化參數(shù)自動(dòng)優(yōu)選的RAKSVD方法能提高弱信號(hào)的去噪效果和識(shí)別能力。

圖4 RAKSVD算法不同正則化參數(shù)去噪結(jié)果

從表1可以看出:AKSVD相較KSVD,其計(jì)算時(shí)間縮減了67%,去噪效果有所下降;正則化參數(shù)自動(dòng)優(yōu)選的RAKSVD方法,其計(jì)算時(shí)間與AKSVD差不多,但去噪效果明顯提升。在實(shí)際應(yīng)用中,選用本文中提出的正則化參數(shù)自動(dòng)優(yōu)選的RAKSVD去噪方法,去噪后既能保證地震弱信號(hào)不發(fā)生畸變,達(dá)到較好的去噪效果,同時(shí)計(jì)算效率還可得到明顯提升。

表1 不同方法的去噪效果和運(yùn)行時(shí)間對比

5 實(shí)際資料處理

為了進(jìn)一步測試本文的去噪方法,將其應(yīng)用于實(shí)際資料的去噪處理。圖5為實(shí)際地震剖面,從圖5中可以看出隨機(jī)噪聲對弱同相軸的影響較大,強(qiáng)軸之間的弱信號(hào)同相軸不清晰、不連續(xù)或錯(cuò)斷。

圖5 實(shí)際地震剖面

為了消除隨機(jī)噪聲的影響,分別采用KSVD、AKSVD以及不同正則化參數(shù)的RAKSVD方法對圖5進(jìn)行了去噪處理,圖6為不同方法去噪后的地震剖面,表2為不同方法的去噪效果和運(yùn)行時(shí)間,其中RAKSVD方法指的是參數(shù)自動(dòng)優(yōu)選的RAKSVD方法。

從圖6可以看出,正則化參數(shù)的選取對去噪效果存在著重要影響,其中μ為0時(shí)去噪效果最差,峰值信噪比為10.658 dB;一旦加入正則項(xiàng),去噪效果就會(huì)有所提升。如果正則化參數(shù)取固定值,當(dāng)μ取0.1時(shí),峰值信噪比為12.064 dB,能夠取得較好的去噪效果;而當(dāng)μ分別取0.3和0.01時(shí),去噪效果不如μ取0.1,說明能夠達(dá)到最佳去噪效果的正則化參數(shù)μ介于0.01～0.3,采用正則化參數(shù)自動(dòng)優(yōu)選的RAKSVD方法,可避免由于μ取固定值不準(zhǔn)確對去噪效果的影響,每次迭代自動(dòng)選取μ值,迭代結(jié)束后峰值信噪比提升到12.823 dB,達(dá)到了最佳去噪效果。

從表2可以看出,AKSVD方法相較KSVD方法,其計(jì)算效率得到了明顯提升,但去噪效果下降明顯,此時(shí)峰值信噪比為10.658 dB;RAKSVD方法在提升計(jì)算效率的同時(shí),還提升了去噪效果。對比圖6綠色方框中的強(qiáng)信號(hào)和青色方框中的弱信號(hào)的去噪結(jié)果可以看出,KSVD方法會(huì)產(chǎn)生過度去噪問題,對強(qiáng)軸的影響不是很大,但弱信號(hào)發(fā)生了明顯的畸變,因此該方法對弱地震信號(hào)產(chǎn)生較大影響;AKSVD方法去噪后的峰值信噪比相比于KSVD方法有所降低,因此總體去噪效果不如KSVD方法,但能夠保留弱信號(hào)特征,這是因?yàn)锳KSVD使用了所有奇異值分量;正則化參數(shù)自動(dòng)優(yōu)選的RAKSVD去噪方法在達(dá)到最佳去噪效果的同時(shí),能夠使弱信號(hào)產(chǎn)生畸變的程度降到最低,從而有利于對弱信號(hào)的提取和識(shí)別。

圖6 實(shí)際地震剖面去噪結(jié)果對比

表2 不同方法的去噪效果和運(yùn)行時(shí)間對比

6 結(jié)束語

在對低信噪比地震弱信號(hào)進(jìn)行去噪處理時(shí),KSVD方法有其獨(dú)特的優(yōu)勢,但含噪信號(hào)存在病態(tài)問題,對于弱信號(hào)容易產(chǎn)生去噪過度問題,且計(jì)算效率較差;AKSVD通過修改更新的原子和相關(guān)參數(shù),極大地提高了計(jì)算效率,但其去噪效果有所下降;將正則化方法引入到AKSVD之后,其去噪效果有明顯提升,但正則化參數(shù)的選取對去噪效果存在著重要影響,通過L曲線擬合自動(dòng)優(yōu)化選取正則化參數(shù),能顯著提高弱地震信號(hào)的去噪效果。在實(shí)際應(yīng)用中可采用正則化參數(shù)自動(dòng)優(yōu)選的RAKSVD方法,不僅能使強(qiáng)反射信號(hào)達(dá)到較好的去噪效果,而且加強(qiáng)了對弱信號(hào)的保護(hù),去噪后地震弱信號(hào)基本不發(fā)生畸變,從而較好地解決了傳統(tǒng)去噪方法中弱信號(hào)保真度不高、容易發(fā)生畸變的問題,在地震弱信號(hào)的提取和識(shí)別中具有廣闊的應(yīng)用前景。