正則化參數(shù)自動優(yōu)選的RAKSVD方法在地震弱信號去噪中的應(yīng)用

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

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

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

1 RAKSVD方法基本原理

1.1 KSVD算法

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

(1)

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

(2)計算殘差:

(2)

式中,Ek代表剩余殘差。

(3)

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

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

xk=σ1ν1,dk=u1.

(4)

1.2 AKSVD算法

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

(5)

(6)

1.3 RAKSVD算法

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

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

(7)

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

(8)

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

(9)

(10)

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

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

2.1 方法原理

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

(11)

其中

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

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

(12)

式中,u為分段優(yōu)選正則化參數(shù)的控制節(jié)點序號;μ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è)計弱信號的褶積模型(圖1(a)),該模型共有3層,上層為強軸,中間層為傾斜層,由上到下逐漸變?nèi)?并與上下兩層斜交,下層為弱信號,共25道,每道300個采樣點,其中雷克子波主頻為25 Hz,模型中存在嚴重高斯白噪聲(圖1(b)),峰值信噪比為1.2 dB,用來進行字典訓練。

圖1 褶積噪音模型

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

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

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

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

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

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

(16)

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

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

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

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

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

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

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

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

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

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

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

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

式中,I為單位矩陣。

4 弱信號模型測試

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

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

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

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

表1 不同方法的去噪效果和運行時間對比

5 實際資料處理

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

圖5 實際地震剖面

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

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

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

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

表2 不同方法的去噪效果和運行時間對比

6 結(jié)束語

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