基于字典學(xué)習(xí)和ADMM的地震數(shù)據(jù)重建

李 慧,韓立國(guó),張 良,賈 帥

(吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院,吉林長(zhǎng)春130026)

重建缺失的地震數(shù)據(jù)是一個(gè)重要的研究課題,傳統(tǒng)的重建方法有波動(dòng)方程法[1]、預(yù)測(cè)濾波法[2]和數(shù)學(xué)變換法[3]等。波動(dòng)方程法通常需要先驗(yàn)參數(shù)信息,這限制了它在實(shí)際工作中的應(yīng)用。預(yù)測(cè)濾波法重建的地震數(shù)據(jù)一般存在誤差較大的問(wèn)題。數(shù)學(xué)變換法包括了采用固定基和字典學(xué)習(xí)兩類(lèi)方法,固定基重建計(jì)算速度較快但是缺乏對(duì)不同的地震數(shù)據(jù)進(jìn)行自適應(yīng)稀疏表示的能力[4-6],字典學(xué)習(xí)能夠?qū)Σ煌牡卣饠?shù)據(jù)進(jìn)行有效的稀疏表示,但是訓(xùn)練字典較為費(fèi)時(shí)[7-9]。

近十年發(fā)展起來(lái)的壓縮感知理論[10-11]對(duì)地震數(shù)據(jù)的重建工作有著重要的指導(dǎo)意義,壓縮感知理論表明:在信號(hào)具有稀疏性的情況下,采集少量的樣本通過(guò)非線性重構(gòu)算法可以重建完整的信號(hào)。目前,壓縮感知理論在地震數(shù)據(jù)處理中已取得了較大進(jìn)展:HERRMANN等[12]提出了基于壓縮感知理論的curvelet變換方法,成功重建了缺失的地震數(shù)據(jù);ZHANG等[13]提出了基于壓縮感知的傅里葉變換重建地震數(shù)據(jù)的方法;LIANG等[14]提出了在壓縮感知理論下采用shearlet變換重建缺失地震數(shù)據(jù)的方法。在壓縮感知中,稀疏基的選擇十分重要,直接影響到地震數(shù)據(jù)的重建精度[15]。字典學(xué)習(xí)能夠根據(jù)不同的地震數(shù)據(jù)訓(xùn)練出相應(yīng)的稀疏表示字典,ZHOU等[7]通過(guò)字典學(xué)習(xí)方法實(shí)現(xiàn)了混采數(shù)據(jù)分離,ZHU等[16]將字典學(xué)習(xí)方法引入地震數(shù)據(jù)去噪中,取得了良好的效果。

BOYD等[17]提出了交替方向乘子算法(ADMM),該方法能夠有效地處理稀疏反演問(wèn)題。STERNFELS等[18]提出了在壓縮感知理論下采用ADMM進(jìn)行奇異值閾值的地震數(shù)據(jù)去噪方法。HARIKUMAR等[19]提出了采用ADMM更新過(guò)完備字典,并在處理變化幅度大的信號(hào)時(shí)取得了良好的效果。ZHOU等[20]將采用ADMM更新過(guò)完備字典的方法成功地應(yīng)用于混采數(shù)據(jù)分離中。

本文在壓縮感知理論下提出了一種基于字典學(xué)習(xí)和ADMM的地震數(shù)據(jù)重建方法,首先對(duì)缺失的地震記錄進(jìn)行字典學(xué)習(xí)得到稀疏表示字典,然后根據(jù)炮記錄中地震道的缺失情況設(shè)計(jì)采樣矩陣,接著建立相應(yīng)的L1范數(shù)約束模型,最后采用ADMM進(jìn)行求解得到重建的地震數(shù)據(jù)。用正演模擬數(shù)據(jù)和實(shí)際地震資料驗(yàn)證了本文方法的有效性。

1 基本原理

1.1 壓縮感知

地震數(shù)據(jù)的重建問(wèn)題可以表示為:

(1)

式中:y∈RM表示原始未缺失的M維地震數(shù)據(jù);x∈RN表示N(N

采用訓(xùn)練得到的過(guò)完備字典D(D∈RM×N且M

(2)

公式(2)是一個(gè)欠定問(wèn)題,壓縮感知理論指出,當(dāng)信號(hào)的稀疏度為r,并且φD滿足有限等距性質(zhì)(restricted isometry property,RIP)的常數(shù)為νr∈(0,1)時(shí):

(3)

公式(2)有唯一解,但是通常情況下很難通過(guò)公式(3)進(jìn)行嚴(yán)格驗(yàn)證,而B(niǎo)ARANIUK[21]放寬了這一條件,當(dāng)φ與D不相關(guān)時(shí),φD滿足公式(3)中條件的概率很大,φD的相關(guān)性定義為:

(4)

式中:i,j分別表示φD中的行和列。GLEICHMAN等[22]驗(yàn)證了過(guò)完備字典和0,1組成的對(duì)角矩陣的不相關(guān)性,此時(shí),公式(2)可以描述為如下數(shù)學(xué)問(wèn)題:

(5)

(6)

式中:λ為正則化參數(shù)。在滿足公式(4)的條件下,公式(6)可以改寫(xiě)為:

(7)

(8)

1.2 字典學(xué)習(xí)

本文采用了K-奇異值分解(K-singular value decomposition,K-SVD)來(lái)獲取公式(2)中的過(guò)完備字典D,K-SVD是一個(gè)字典學(xué)習(xí)的優(yōu)化方法[23],它的優(yōu)化過(guò)程可以表示為:

(9)

式中:F∈RM×N表示輸入的訓(xùn)練樣本;A為樣本在字典中的稀疏系數(shù);T表示非零元素個(gè)數(shù)上限;al為A的第l行。首先需要設(shè)定一個(gè)初始字典,本文的初始字典選取了過(guò)完備的DCT字典,如圖1所示。

然后,固定初始字典D,輸入訓(xùn)練樣本F,求出稀疏系數(shù)al:

(10)

圖1 DCT過(guò)完備字典

式中:fl為F的列。在求得稀疏系數(shù)al的同時(shí),固定al,對(duì)字典D進(jìn)行逐列更新:

(11)

(12)

1.3 交替方向乘子算法

與K-SVD采用正交匹配追蹤算法求解公式(6)不同,為了提高計(jì)算精度,通過(guò)對(duì)L0約束問(wèn)題向L1約束問(wèn)題轉(zhuǎn)化,本文采用ADMM對(duì)公式(7)進(jìn)行求解。

公式(7)是一個(gè)凸優(yōu)化問(wèn)題,ADMM引入了輔助變量p,h來(lái)求解(7)式:

(13)

式中:f(h)(h∈R)要求為凸函數(shù)并且光滑;g(p)(p∈R)要求為凸函數(shù)。將公式(7)寫(xiě)成ADMM形式:

(14)

其增廣目標(biāo)函數(shù)為:

(15)

式中:τ為增廣拉格朗日參數(shù)并且大于0,它控制著迭代步長(zhǎng)。將公式(15)寫(xiě)成增廣拉格朗日函數(shù):

(16)

式中:Λ為對(duì)偶變量。ADMM的第i次迭代形式為:

(17)

其中,I表示單位矩陣,soft為軟閾值函數(shù),定義:

(18)

2 技術(shù)流程

本文技術(shù)流程實(shí)現(xiàn)步驟如下。

1) 輸入不完整的地震數(shù)據(jù)作為訓(xùn)練樣本,初始字典選為DCT過(guò)完備字典,采用的字典學(xué)習(xí)方法為K-SVD,輸出相應(yīng)的稀疏字典D。

2) 利用D對(duì)地震數(shù)據(jù)進(jìn)行稀疏表示,采用ADMM求解基于壓縮感知理論建立的數(shù)學(xué)模型(公式(7))。

3 模擬數(shù)據(jù)實(shí)驗(yàn)

本文采用信噪比(RSN)作為地震數(shù)據(jù)重建性能的判斷標(biāo)準(zhǔn):

(19)

圖2 模擬數(shù)據(jù)(a)和隨機(jī)缺失50%后的數(shù)據(jù)(b)

對(duì)圖2b中缺失50%后的模擬數(shù)據(jù)進(jìn)行K-SVD訓(xùn)練,可以得到相應(yīng)的稀疏字典,如圖3所示。

對(duì)圖2b中缺失50%的地震數(shù)據(jù)進(jìn)行重建。分別采用curvelet稀疏基和OMP、curvelet稀疏基和ADMM、字典學(xué)習(xí)和OMP、字典學(xué)習(xí)和ADMM重建模擬數(shù)據(jù)的結(jié)果及對(duì)應(yīng)的重建誤差如圖4至圖7所示。這4種重建方法的RSN分別為17.7138,21.8057,21.9478,24.9869dB,耗時(shí)分別為2.56,4.32,132.48,271.09s。從計(jì)算效率看,當(dāng)需要處理的數(shù)據(jù)量較大時(shí),會(huì)面臨耗時(shí)過(guò)長(zhǎng)的問(wèn)題,若數(shù)據(jù)量較小且需要高精度時(shí),選取本文方法比較有利。其中,ADMM的正則化參數(shù)λ取1,控制步長(zhǎng)的參數(shù)τ取1,ADMM中的初始a,p,Λ設(shè)為0。

從圖4至圖7可以看出,對(duì)于正演模擬數(shù)據(jù),基于壓縮感知理論的固定基和字典學(xué)習(xí)方法都可以重建地震數(shù)據(jù),而字典學(xué)習(xí)重建的性能比固定基重建的性能更好,重建誤差也小。在壓縮感知的重建算法中,ADMM重建的性能比OMP重建的性能更優(yōu)越,無(wú)論是在固定基中還是在字典學(xué)習(xí)中,ADMM重建的誤差均小于OMP重建的誤差。圖8是采用4種方法對(duì)不同缺失率的地震數(shù)據(jù)插值后的信噪比曲線,可以看出,在缺失10%～50%地震數(shù)據(jù)的情況下,采用本文方法重建的效果比采用其它幾種方法的效果好,且有較強(qiáng)的魯棒性。

圖3 模擬數(shù)據(jù)的K-SVD字典

圖4 采用curvelet稀疏基和OMP重建模擬數(shù)據(jù)的結(jié)果(a)及重建誤差(b)

圖5 采用curvelet稀疏基和ADMM重建模擬數(shù)據(jù)的結(jié)果(a)及重建誤差(b)

圖6 采用字典學(xué)習(xí)和OMP重建模擬數(shù)據(jù)的結(jié)果(a)及重建誤差(b)

圖7 采用字典學(xué)習(xí)和ADMM重建模擬數(shù)據(jù)的結(jié)果(a)及重建誤差(b)

圖8 采用4種方法對(duì)不同缺失率的地震數(shù)據(jù)插值后的信噪比曲線

4 實(shí)際數(shù)據(jù)實(shí)驗(yàn)

圖9a為缺失50%的實(shí)際數(shù)據(jù)。選取100道數(shù)據(jù),采樣點(diǎn)數(shù)為400個(gè),時(shí)間采樣率為0.001s,道間距為10m,將缺失的實(shí)際數(shù)據(jù)作為訓(xùn)練樣本,得到的稀疏字典如圖9b所示。

分別采用curvelet稀疏基和OMP、curvelet稀疏基和ADMM、字典學(xué)習(xí)和OMP、字典學(xué)習(xí)和ADMM 4種方法對(duì)圖9a中缺失50%的實(shí)際數(shù)據(jù)進(jìn)行重建,其結(jié)果及對(duì)應(yīng)的插值部分如圖10至圖13所示。這4種方法的耗時(shí)分別為4.15,7.52,182.07,411.30s。其中ADMM中的正則化參數(shù)λ取1,控制步長(zhǎng)的參數(shù)τ取1,ADMM中的初始a,p,Λ設(shè)為0。

圖9 實(shí)際數(shù)據(jù)(a)和相應(yīng)的字典(b)

圖10 采用curvelet稀疏基和OMP重建實(shí)際數(shù)據(jù)的結(jié)果(a)及對(duì)應(yīng)的插值部分(b)

圖11 采用curvelet稀疏基和ADMM重建實(shí)際數(shù)據(jù)的結(jié)果(a)及對(duì)應(yīng)的插值部分(b)

圖12 采用字典學(xué)習(xí)和OMP重建實(shí)際數(shù)據(jù)的結(jié)果(a)及對(duì)應(yīng)的插值部分(b)

圖13 采用字典學(xué)習(xí)和ADMM重建實(shí)際數(shù)據(jù)的結(jié)果(a)及對(duì)應(yīng)的插值部分(b)

從圖10至圖13可以看出,壓縮感知理論下的固定基和字典學(xué)習(xí)方法都可以很好地重建實(shí)際地震數(shù)據(jù),其中,采用字典學(xué)習(xí)比采用固定基的效果更好,可以對(duì)一些細(xì)節(jié)部分進(jìn)行較好的重建,而壓縮感知理論下重建算法中的ADMM的重建性能比OMP的重建性能更好,無(wú)論是能量較強(qiáng)的同相軸還是能量較弱的細(xì)節(jié)部分都可以得到有效重建。

5 結(jié)論

本文提出了一種壓縮感知下基于K-SVD字典學(xué)習(xí)和ADMM的地震數(shù)據(jù)重建方法,并建立了技術(shù)實(shí)現(xiàn)流程。正演模擬數(shù)據(jù)和實(shí)際地震資料測(cè)試結(jié)果表明,本文方法與固定基地震數(shù)據(jù)插值方法以及字典學(xué)習(xí)中常用的OMP插值方法相比,重建精度有很大提高。并且,在地震數(shù)據(jù)缺失10%～50%時(shí),采用本文方法可以較好地對(duì)地震數(shù)據(jù)進(jìn)行重建,重建精度較高。但本文方法有耗時(shí)過(guò)長(zhǎng)的缺陷,當(dāng)面臨較大數(shù)據(jù)量的插值任務(wù)時(shí),這一缺陷會(huì)更大,如何減少耗時(shí)、提高算法的運(yùn)行效率是下一步的研究方向。