自適應秩收斂低秩算法壓制沙漠地震隨機噪聲

李 佳馬海濤李 月

(吉林大學 通信工程學院,長春 130012)

0 引 言

沙漠地震資料中包含大量的復雜噪聲造成地震記錄低信噪比和低分辨率,為后續(xù)的分析研究帶來不便。同時,由于沙漠噪聲的低頻特性使有效信號與噪聲的頻帶產(chǎn)生混疊,造成二者在頻域進行分離時較為困難。

近年來,許多學者提出了一系列的噪聲壓制算法,如f-x反褶積[1]、f-k濾波[2]、小波變換[3]、Shearlet變換[4]、曲波閾值法[5]和非局部均值濾波[6-7]等。Lin等[8]利用時頻峰值濾波壓制地震隨機噪聲。雖然這些算法應用于地震信號分析并取得了一定效果,但在處理復雜多樣的沙漠地震數(shù)據(jù)時表現(xiàn)出明顯的不足。針對更加復雜的地質(zhì)條件和更高要求的勘探目標,傳統(tǒng)地震數(shù)據(jù)去噪技術已難以滿足能源勘探高信噪比、高分辨率、高保真度的需求。

為尋求一種更適合沙漠地震隨機噪聲的去噪方法,筆者提出一種將變分模態(tài)分解(VMD:Variational Mode Decomposition)算法[9]與混合高斯魯棒主成分分析(MoG-RPCA:Mixture of Gauss-Robust Principal Component Analysis)低秩分解[10]相結(jié)合的自適應秩收斂算法,并將其應用于沙漠地區(qū)隨機噪聲消減[11-12]。針對RPCA(Robust Principal Component Analysis)低秩分解算法對噪聲分析過于簡單,提出了加入混合高斯模型的MoG-RPCA低秩分解算法,該方法在視頻圖像觀測中得到廣泛應用并取得良好效果。在利用低秩分解對沙漠地震數(shù)據(jù)進行處理時,秩的數(shù)值大小直接影響分解效果的好壞。為更好地對秩數(shù)進行選擇,在對信號矩陣進行低秩分解后,計算經(jīng)分解產(chǎn)生的低秩和稀疏兩部分矩陣得到分解誤差,當分解誤差滿足預設誤差閾值時完成迭代,從而實現(xiàn)秩數(shù)的自適應選擇[13-14]。通過上述方法,既規(guī)避了由于沙漠地震噪聲復雜特性造成的常規(guī)方法不適用問題,同時又無需針對不同地震記錄多次調(diào)試秩數(shù)大小,進一步保證低秩分解效果。通過模擬實驗以及對實際資料的處理,證明筆者提出的方法可以有效壓制沙漠地震隨機噪聲。

1 隨機噪聲消減算法

筆者選擇MoG-RPCA算法對信號進行低秩分解,該算法針對RPCA算法對噪聲分析過于簡單問題,加入了混合高斯模型,使其更適合處理復雜噪聲數(shù)據(jù),同時為了降低矩陣秩數(shù)選擇對低秩分解效果的影響,采用VMD算法對沙漠地震數(shù)據(jù)進行分解,并將分解得到所有模態(tài)重排成一個新的信號矩陣,利用自適應秩數(shù)選擇的方法,避免對不同地震記錄多次調(diào)試秩數(shù)大小,使去噪算法更具適用性。

1.1 VMD算法原理

VMD算法作為一種新型的自適應模態(tài)分解方法,主要通過求解模態(tài)分量的變分問題確定分量信號的帶寬和中心頻率,并將輸入信號分解出一系列稀疏的分量信號。構造變分量問題時,其基本原理是對分量信號采用希爾伯特變換、混頻等信號處理方法。在各階模態(tài)分量之和等于原信號的前提下,將模態(tài)分量的變分問題轉(zhuǎn)化為尋求估計帶寬最小和的模態(tài)函數(shù)。具體步驟如下[9]。

多分量信號x(t)可表示如下

其中xk(t)為分量信號,k=1,…,N,N為分解模態(tài)的數(shù)量。首先,通過對分量信號xk(t)進行希爾伯特變換得到單邊頻譜,并將各分量信號的頻譜與e-jωkt相乘使其調(diào)整到以預測中心頻率為中心的頻段上。其次,計算混合信號梯度模的平方以估計頻移后分量信號的帶寬。通過上述流程得到如下有約束優(yōu)化問題

為求解式(2)進一步引入一元二次懲罰函數(shù)系數(shù)α和拉格朗日乘子λ,最終得到一個無約束優(yōu)化問題

最后采用交替方向乘子法根據(jù)

其中的傅里葉逆變換實數(shù)部分為分量信號。

交替更新,,λ(n+1)求解式(3)變分問題,直到滿足如下收斂條件

其中ε為收斂值,筆者將其設置為1×10-7。

1.2 MoG-RPCA算法原理

為處理復雜噪聲,在魯棒主成分分析算法[15]的基礎上提出了加入混合高斯(MoG:Mixture of Gauss)的魯棒主成分分析,即MoG-RPCA,該算法更適合對復雜噪聲進行處理,精確度也比RPCA更高[10]。

將MoG-RPCA模型視為一個合成模型,則可表示為

其中L為低秩分量,E為稀疏分量。

針對稀疏分量部分,假設E中的每個eij都遵循MoG分配,即

其中πk為混合比例且為高斯分組數(shù),為平均值為μ和精度為τ的高斯分布。通過引入指標變量zijk可將式(9)等價為

由于自動關聯(lián)決策(ARD:Automatic Relevance Determination)具有快速和可擴展性良好的優(yōu)點,采用其為低秩分量部分進行建模。秩為l≤min(m,n)的低秩矩陣L∈R m×n作為矩陣U∈R m×R和V∈R n×R經(jīng)奇異值分解(SVD:Singular Value Decomposition)得到的結(jié)果,公式如下

其中R>l,u·r(v·r)為矩陣U(V)的第r列。為實現(xiàn)U和V中的列稀疏,讓U和V中某些列向量趨近于零,進一步保證矩陣L的低秩特性。

為求解式(11)引入如下的先驗公式到U和V中

其中I m為m×m的單位矩陣。每個精度變量γr的同源先驗可表示為

記矩陣U和V中的每列對u·r,v·r具有以共同精度變量γr為特征的相似稀疏度曲線,經(jīng)驗證明此模型會產(chǎn)生較大的精度值γr,使矩陣L具有較好的低秩估計[16]。結(jié)合式(8)和式(10)~式(13)建立基于混合高斯模型的魯棒主成分分析,即MoG-RPCA,則目標轉(zhuǎn)換為求解所有相關變量的后驗分布

其中Z={zij};μ=(μ1,μ2,…,μK);τ=(τ1,τ2,…,τK);γ=(γ1,γ2,…,γR)。

為從目標矩陣Y中恢復低秩矩陣L,采用變分貝葉斯(VB:Variational Bayes)理論[17]求解U,V和γ。VB通過求解變分優(yōu)化得到真實后驗為觀測數(shù)據(jù))的近似分布q(x),則式(14)所示后驗分布可用如下因子化公式近似表示

其中u i·(v j·)為矩陣U(V)的第i(j)行。對矩陣U中每行u i·,經(jīng)式(15)計算得

其中μu i·為均值為協(xié)方差,

同理可求得矩陣V中每行v j·的近似分布q(v j·)。對控制矩陣L秩大小的變量γ,根據(jù)式(13)其近似分布可表示為

1.3 自適應迭代收斂算法原理

在利用低秩分解算法進行沙漠地震數(shù)據(jù)處理時,秩數(shù)大小直接影響分解效果好壞。為更好地對秩數(shù)進行選擇,利用低秩分解產(chǎn)生的低秩和稀疏兩部分矩陣計算分解誤差,當分解誤差滿足預設誤差閾值時完成迭代,從而實現(xiàn)秩數(shù)的自適應選擇,具體步驟如下:

1)根據(jù)式(8)信號矩陣Y經(jīng)低秩分解得到低秩部分L和稀疏部分E,記低秩分解誤差為Eerror,利用

計算分解誤差。

2)第i次低秩分解完成后,判斷分解誤差Eerrori是否滿足誤差閾值,如果不滿足則需要進行第i+1次迭代,同時秩參數(shù)減1;

3)當Eerror達到閾值條件時迭代停止得到最終的低秩分解矩陣L。

1.4 低頻噪聲消減算法

筆者將VMD算法與MoG-RPCA低秩分解相結(jié)合,并采用自適應迭代方式對沙漠地區(qū)低頻噪聲進行壓制,具體流程如下:

1)設原始地震記錄為x m×n(t),利用VMD對其進行模態(tài)分解重排,將每道信號分解出N個模態(tài),并將所有模態(tài)合成為一個新的信號矩陣X m×nN(t);

2)對X m×nN(t)進行MoG-RPCA分解,得到低秩部分L和稀疏部分S;

3)根據(jù)式(18)計算分解誤差Eerror;

4)在第i次MoG-RPCA分解完成后,判斷分解誤差Eerrori是否滿足誤差閾值(文中設置Eerror<0.02),如果不滿足條件則需要進行第i+1次迭代,同時秩參數(shù)減1,當Eerror達到閾值條件時迭代停止得到最終的低秩矩陣L;

5)將4)中所得低秩矩陣L每N列進行疊加并與原始記錄作差得到最終的去噪結(jié)果。

2.1 模擬地震記錄

為檢測該方法去噪性能,利用筆者算法處理一幅50道模擬地震記錄。如圖1a所示,其采樣頻率為500 Hz,采樣點個數(shù)為700個,采樣時間為1 400 ms,信號主頻為30 Hz,幅值為1.0。由圖1a、圖1b和圖1c可見,噪聲基本上將有效信號同相軸淹沒。對比圖1d和圖1e的傳統(tǒng)去噪方法處理結(jié)果,筆者方法的整體去噪效果更好(見圖1f)。小波變換以及MoG-RPCA雖然也具備一定的去噪效果,但其噪聲消減沒有筆者算法徹底,去噪后有效信號同相軸不夠連續(xù)清晰。因此,筆者算法針對沙漠低頻噪聲具有更好的壓制效果。

圖1 模擬記錄處理結(jié)果Fig.1 The processing results of synthetic desert seismic record

通過如圖2所示的純噪聲記錄與去噪結(jié)果差值記錄對比,可更直觀地展示出幾種去噪方法的保幅情況。

圖2 去噪差值結(jié)果Fig.2 The results of denoising residual

從圖2可以看出,小波變換和MoG-RPCA在差值記錄中均有明顯的有效信號殘留,而筆者算法則沒有出現(xiàn)此問題。小波變換在彎軸部分保幅效果較差,對噪聲壓制也不徹底,而MoG-RPCA在直軸部分有較多有效信號殘留。因此,筆者算法對有效信號保留最為完整。

為進一步展示幾種去噪算法對有效信號的保幅效果,將模擬含噪記錄和各算法去噪結(jié)果中第20道信號進行對比,結(jié)果如圖3所示。從圖3可以看出,小波變換與MoG-RPCA對信號的保幅效果較差,而經(jīng)筆者算法處理后有效信號恢復較為完整,信號幅度保持效果能達到85%以上。

圖3 第20道信號振幅比較Fig.3 Amplitude comparison of signals of the 20th trace

通過以上對比實驗可知,筆者算法在噪聲壓制以及有效信號保留均有較好的效果。進一步對去噪結(jié)果進行量化分析,表1給出在不同噪聲水平下,幾種去噪方法處理前后的信噪比比較,使用的信噪比計算公式為

其中s(t)和x(t)分別為原始純凈信號和去噪后信號。由表1可知,經(jīng)筆者算法去噪處理后的結(jié)果信噪比提升最明顯。

表1 不同噪聲水平下的信噪比對比Tab.1 Comparison of SNR with different desert noise levels (dB)

從頻域角度驗證筆者算法性能,圖4為整幅模擬含噪記錄以及各個算法去噪結(jié)果的頻率-波數(shù)譜圖(FK譜圖),圖5為去噪差值FK譜圖。

圖4 模擬記錄FK譜圖Fig.4 FK spectrum of synthetic desert seismic record

從圖4d、圖4e可以看出,小波變換對低頻噪聲壓制并不徹底,且對有效信號造成了明顯的損失,MoG-RPCA雖然對有效信號的恢復效果有所提升,但也有較多噪聲殘留;而從圖4f可以看出,筆者算法對低頻噪聲壓制較為徹底,基本沒有造成有效信號能量損失。

圖5進一步驗證了筆者算法在噪聲壓制以及有效信號恢復的優(yōu)越性。從圖5b、圖5c可看出,小波變換對有效信號損害嚴重,并且對有效信號與噪聲混疊部分壓制效果較差;MoG-RPCA對噪聲具備一定壓制效果但是并不徹底,同時也對有效信號造成一定能量損失。對比可知,筆者算法在噪聲壓制以及有效信號保留均有較好的效果。

圖5 去噪差值FK譜圖Fig.5 FK spectrum of denoising residual

2.2 實際地震記錄

實際地震記錄處理結(jié)果如圖6所示。

圖6 實際記錄處理結(jié)果Fig.6 The processing results of real desert seismic record

為了驗證筆者算法在實際地震勘探中的有效性,對圖6a所示的一幅210道實際地震記錄進行處理,其采樣頻率為500 Hz,采樣時間為2 400 ms。

由圖6中方框標出部分可知,圖6b所示小波變換具有一定去噪效果,但對面波壓制不夠徹底,且有較多噪聲殘留。圖6c所示MoG-RPCA算法同樣存在面波壓制不徹底的問題,并且同相軸不連貫,有效信號保幅較差。經(jīng)筆者算法處理后結(jié)果如圖6d所示,該記錄中的初至噪聲以及面波均得到有效壓制,同相軸清晰且連續(xù),證明筆者算法在噪聲壓制以及有效信號恢復均有較好效果。

3 結(jié) 語

針對沙漠地震噪聲復雜特性,筆者提出一種將VMD算法與MoG-RPCA算法結(jié)合的自適應秩收斂低秩算法,實現(xiàn)沙漠地震低頻噪聲消減。該方法自適應性強,具有較好特征識別效果。相比于RPCA算法,MoG-RPCA可以對更復雜的噪聲進行處理,且與VMD算法結(jié)合改進了因頻譜混疊導致MoG-RPCA分解所得低秩矩陣包含大量信號的問題,同時采用自適應迭代算法對秩數(shù)進行自動選擇,使低秩分解更精確,提高有效信號保幅效果。模擬記錄處理結(jié)果表明,筆者算法對沙漠地震記錄中隨機噪聲具有較好消減效果,信噪比提升明顯且有效信號幅度衰減較小;實際地震記錄處理結(jié)果表明,筆者算法在降噪同時保留了有效頻率成分,同相軸清晰連續(xù)。綜上可得,筆者提出的基于VMD與MoG-RPCA相結(jié)合的自適應迭代收斂低秩算法在沙漠地震去噪方面實現(xiàn)簡單,具有較好的實用性。