程文婷,方文倩,付麗華

(中國地質大學(武漢)數(shù)學與物理學院,湖北武漢430074)

地震數(shù)據(jù)隨機噪聲壓制方法較多,主要有基于稀疏變換的方法,如Fourier變換[1]、Wavelet變換[2-3]、Curvelet變換[4]、Radon變換[5]、Shearlet變換[6]、Seislet變換[7]等;基于濾波理論的方法,如中值濾波[8-9]、f-x反褶積[10]、聚束濾波[11]等;基于波動方程的方法[12]和基于字典學習的方法[13-14]等。近年來,矩陣降秩作為一種重要的數(shù)據(jù)處理方法,成為地震信號處理領域的研究熱點[15-17]。

基于降秩理論的地震數(shù)據(jù)去噪原理是:不含噪聲的地震數(shù)據(jù)具有低秩結構,而隨機噪聲的存在會增加數(shù)據(jù)矩陣的秩。因此,地震數(shù)據(jù)去噪問題可以轉化為矩陣降秩問題。奇異譜分析(singular specturm analysis,SSA)[18]是經(jīng)典的降秩方法,通過對地震數(shù)據(jù)頻率切片構造的Hankel矩陣進行降秩達到去噪目的。由于奇異值分解(singular value decomposition,SVD)計算時間長,OROPEZA等[19]以隨機SVD代替SVD,提高了計算效率。但是,此類方法需要事先確定秩的大小,而秩的選取將直接影響數(shù)據(jù)的去噪效果。核范數(shù)最小化(nuclear norm minimization,NNM)方法是另一類低秩矩陣近似方法。核范數(shù)作為秩函數(shù)的凸松弛函數(shù),在其基礎上構建的目標函數(shù)是凸問題,可使用凸優(yōu)化方法求解,如:奇異值閾值法(singular value threshold,SVT)[20]、加速近端梯度法[21]、交替乘子方向法(alternating direction method of multipliers,ADMM)[22]、不動點迭代法(fixed point continuation,FPC)[23]等。然而NNM方法是將所有奇異值的和最小化,求出的解往往只是原始秩最小化問題的次優(yōu)解。為了解決該問題,HU 等[24]提出截斷核范數(shù)方法,對[min(m,n)-r]個較小奇異值的和進行最小化處理,其中m,n為矩陣的尺度,r為矩陣的秩。相較于核范數(shù),截斷核范數(shù)能更好地逼近秩。

近年來,非局部自相似先驗在圖像去噪領域取得了巨大的應用成果,如BUADES等[25]提出的非局部均值(non-local means,NLM)算法以及DABOV等[26]提出的三維塊匹配(block matching and 3D collaborative,BM3D)算法等。事實上,地震數(shù)據(jù)相鄰道之間同相軸的時差變化規(guī)律相近,地震記錄在縱向和橫向上存在較強的相似性,自相似塊構造的矩陣具有更好的低秩性能。BONAR等[27]首次利用NLM算法壓制地震數(shù)據(jù)隨機噪聲。SHANG等[28]提出一種自適應濾波參數(shù)估計的NLM算法,提高了自相似塊的搜索效率。張巖等[29]提出非局部自相似與字典學習相結合的算法,通過對基于多道相似組的地震數(shù)據(jù)構建字典和稀疏編碼,提高稀疏性壓制隨機噪聲。ZHANG等[30]將二維相似塊組排列成三維順序結構,增強自相似塊的相干特性,進而利用濾波對地震數(shù)據(jù)隨機噪聲進行壓制。

本文提出基于自相似性和低秩先驗的地震數(shù)據(jù)去噪方法,通過塊匹配算法搜索地震數(shù)據(jù)的自相似塊,然后以“自相似塊組”為單元,利用低秩約束進行隨機噪聲壓制。TNNR比經(jīng)典的核范數(shù)更加接近原始的秩函數(shù),但由于該問題是非凸問題,因此,我們將目標函數(shù)進行轉換,然后利用APGL算法進行求解。

1 方法理論

1.1 基于降秩理論的地震數(shù)據(jù)去噪模型

地震數(shù)據(jù)去噪模型可表示為:

M=X+E

(1)

式中:M和X分別表示含噪地震數(shù)據(jù)和待恢復的無噪數(shù)據(jù);E為隨機噪聲?；诮抵壤碚摰牡卣饠?shù)據(jù)去噪模型采用如下目標函數(shù):

(2)

式中:rank(X)為X的秩,表示X非零奇異值的個數(shù);‖·‖F(xiàn)為Frobenious范數(shù);λ為正則化參數(shù)。

秩函數(shù)是非凸非光滑函數(shù),導致公式(2)的求解是一個NP難問題。可利用核范數(shù)近似秩函數(shù),將公式(2)轉化為凸優(yōu)化問題:

(3)

1.2 基于截斷核范數(shù)的低秩模型

(4)

因為公式(4)是非凸的,所以先利用下文定理1對目標函數(shù)進行轉換,然后再利用APGL算法求解。

1.3 目標函數(shù)轉換

定理1[24]:已知矩陣X∈Rm×n,存在矩陣A∈Rr×m,B∈Rr×n,滿足AAT=Ir×r,BBT=Ir×r,對任意非負整數(shù)r(r≤min(m,n)),有:

(5)

這樣:

(6)

于是:

(7)

對公式(7)進行迭代求解。令X1=M,第l次迭代時,對Xl進行SVD分解得到Al和Bl;然后固定Al和Bl,更新Xl+1。目標函數(shù)((4)式)可轉化為如下優(yōu)化問題:

(8)

因此,可采用凸優(yōu)化方法求解。APGL算法是經(jīng)典的凸優(yōu)化方法,其優(yōu)勢為收斂速度快,對噪聲具有較強的魯棒性[21]。本文采用APGL算法進行優(yōu)化求解。

1.4 相似塊匹配

相似度定義為歐式距離:

(9)

式中:mi和mj分別表示目標塊和搜索的相似塊,di,j越小,則mi和mj越相似。由于地震數(shù)據(jù)量較大,在全局搜索相似塊計算較為耗時,因此,設置大小為Q×Q的搜索窗,在搜索窗內按照公式(9)所示的相似性度量尋找相似塊。相似塊大小q,相似塊個數(shù)H及搜索窗Q是相似塊匹配中較為重要的參數(shù)。實驗中可以通過多次試誤調節(jié),在計算復雜度和去噪效果之間折衷選擇參數(shù)。

由相似塊組成的矩陣Mi是本文需要恢復的矩陣,因此,公式(8)轉化為基于相似塊矩陣形式的去噪模型:

(10)

圖1 自相似塊匹配與去噪示意

1.5 APGL優(yōu)化求解

APGL算法是一種有效且穩(wěn)定的凸優(yōu)化方法,可以解決如下無約束的非光滑凸問題:

(11)

式中:g為閉凸函數(shù);f為可微凸函數(shù)。

APGL算法將公式(11)寫成如下F(X)在某一定點Y處的近似形式:

(12)

對任意t>0,迭代更新X,Y,t從而優(yōu)化公式(11)。在第k次迭代,通過最小化G(X,Yk)更新Xk+1:

(13)

對于目標函數(shù)(公式(10)),令:

g(X)=‖Xi‖*

則:

(14)

對矩陣X∈Rm×n進行SVD分解:

X=UΣVTΣ=diag({σi}1≤i≤min(m,n))

(15)

矩陣奇異值收縮算子Dτ定義:

(16)

由文獻[20]可知,對任意τ>0,Y∈Rm×n,有:

(17)

因此,公式(14)可轉化為:

(18)

接著,更新ti,k+1和Yi,k+1:

(19)

(20)

本文所提的SP-TNNR算法具體實現(xiàn)過程見表1。

表1 基于自相似性和低秩先驗的地震數(shù)據(jù)隨機噪聲壓制算法實現(xiàn)過程

2 數(shù)據(jù)實驗

利用仿真數(shù)據(jù)和實際地震數(shù)據(jù)驗證SP-TNNR算法的性能,并與傳統(tǒng)的Curvelet變換、f-x反褶積、TNNR算法以及基于自相似塊匹配的核范數(shù)最小化方法(SP-NNM)進行對比分析,用信噪比(RSN)評價算法性能:

(21)

式中:S與S*分別表示不含噪數(shù)據(jù)和去噪后的數(shù)據(jù)。

2.1 仿真數(shù)據(jù)實驗

仿真數(shù)據(jù)實驗選取的地震數(shù)據(jù)共256道,每道256個時間采樣點,采樣間隔為2ms。在綜合考慮算法時間復雜度、計算復雜度及去噪效果等因素的基礎上,SP-TNNR算法實現(xiàn)中,設置搜索窗Q=30,相似塊大小q=8,相似塊個數(shù)H=150,超參數(shù)λ=0.95。f-x反褶積頻帶為1～100Hz,濾波長度為12。而在TNNR方法中秩取15。在仿真和實際數(shù)據(jù)實驗中Curvelet變換尺度參數(shù)均為5,閾值大小取決于去噪數(shù)據(jù)的信噪比。SP-NNM方法中搜索窗、相似塊大小以及相似塊個數(shù)設置均與SP-TNNR方法一致。圖2a為原始仿真疊前地震數(shù)據(jù),圖2b為加入了隨機噪聲(均值為0,標準差為50的高斯白噪)的仿真疊前地震數(shù)據(jù)(RSN=5.1)。圖3和圖4分別為5種方法去噪結果及對應的殘差剖面。Curvelet變換、f-x反褶積、TNNR、SP-NNM以及SP-TNNR方法去噪后的信噪比分別為23.0,20.2,12.5,26.3,27.5,耗時分別為2.8,3.4,860.1,931.7,932.3s。從去噪結果和殘差剖面中可以看出,Curvelet變換去除了絕大部分噪聲,但同相軸有模糊現(xiàn)象且同相軸附近還殘留部分噪聲。f-x反褶積方法去噪結果中仍然可以明顯看到噪聲殘留。TNNR去噪結果中同相軸附近殘留噪聲湮沒了有效信息,去噪效果較差。SP-NNM和SP-TNNR方法去噪后噪聲基本消除,但SP-NNM去噪后同相軸周圍邊緣過于光滑,在去噪的同時去掉了部分有效信號,而SP-TNNR去噪結果更干凈,去噪后同相軸十分清晰,去噪效果最佳。

圖2 原始仿真疊前地震數(shù)據(jù)(a)及加噪數(shù)據(jù)(b)

圖3 5種方法去噪結果a Curvelet變換; b f-x反褶積; c TNNR;d SP-NNM; e SP-TNNR

圖4 應用不同方法對仿真疊前地震數(shù)據(jù)去噪后的殘差剖面a Curvelet變換; b f-x反褶積; c TNNR; d SP-NNM; e SP-TNNR

圖5給出了不同噪聲水平下5種方法去噪結果的信噪比?？梢钥闯?隨著噪聲水平的增加,5種方法的信噪比逐漸降低,而本文提出的SP-TNNR方法去噪效果最優(yōu)。

圖5 不同噪聲水平下5種方法去噪結果的信噪比

2.2 實際數(shù)據(jù)應用

實際疊后地震數(shù)據(jù)來自網(wǎng)站https:∥github.com/sevenysw/MathGeo2018,如圖6a所示,包含256道,每道256個時間采樣點。加入隨機噪聲(均值為0,標準差為50的高斯白噪)的結果見圖6b(RSN=9.0),可以看出加噪數(shù)據(jù)的同相軸邊緣信息模糊。經(jīng)過多次試誤調節(jié),實驗中搜索窗Q=30,相似塊大小q=9,相似塊個數(shù)H=150,正則化參數(shù)λ=0.99。f-x反褶積頻帶為1～100Hz,濾波長度為14。TNNR方法中秩取18。圖7為實際疊后地震數(shù)據(jù)去噪結果。圖7a為Curvelet變換去噪結果,該方法去掉了大部分噪聲,但是同相軸過于光滑,有效信號損失較大;圖7b為f-x反褶積方法去噪結果,部分有效信號和噪聲同時被去掉;圖7c為TNNR方法去噪結果,效果極差;圖7d為SP-NNM方法去噪結果,去掉大部分噪聲的同時也去掉了部分有效信號;圖7e 為本文提出的SP-TNNR方法去噪結果,可以看出噪聲被去掉,同相軸清晰可見且連續(xù)性較好,具有較高的保真度。5種方法去噪后的信噪比依次為19.8,18.9,13.1,20.7,21.9,耗時分別為2.3,3.2,883.2,906.3,912.8s。圖8為圖7對應的殘差剖面。可以看出,本文提出的SP-TNNR方法去噪后有效信息得到了最大程度的保留,信噪比高于另外4種算法。圖9 至圖11為對應的頻率-波數(shù)圖。可以看出,SP-TNNR方法去噪后頻率比較集中,其它4種方法去噪后依然存在頻散現(xiàn)象。

圖6 實際疊后地震數(shù)據(jù)(a)及加噪數(shù)據(jù)(b)

圖7 實際疊后地震數(shù)據(jù)去噪結果a Curvelet變換; b f-x反褶積; c TNNR; d SP-NNM; e SP-TNNR

圖8 實際疊后地震數(shù)據(jù)去噪后的殘差剖面a Curvelet變換; b f-x反褶積; c TNNR; d SP-NNM; e SP-TNNR

圖9 實際疊后地震數(shù)據(jù)(a)和加噪數(shù)據(jù)(b)頻率-波數(shù)圖

圖10 利用 Curvelet變換(a)和 f-x反褶積(b)方法去噪后的地震數(shù)據(jù)頻率-波數(shù)圖

圖11 利用 TNNR(a)、 SP-NNM(b)和SP-TNNR(c)方法去噪后的地震數(shù)據(jù)頻率-波數(shù)圖

3 結論

本文在TNNR基礎上引入非局部自相似先驗,使用相似塊匹配方法構造低秩矩陣,提出了基于自相似性和低秩先驗的地震數(shù)據(jù)隨機噪聲壓制方法,即SP-TNNR算法。與經(jīng)典的核范數(shù)方法相比,本文方法優(yōu)點為:TNNR算法保留最大的幾個奇異值僅懲罰較小奇異值,可以更加準確地逼近秩;與截斷核范數(shù)方法相比,本文方法優(yōu)點為:相似塊匹配方法構造低秩矩陣充分利用地震數(shù)據(jù)的自相似性,通過相似塊構建的矩陣具有更好的低秩性能,加入自相似先驗的截斷核范數(shù)方法比沒有增添自相似先驗的截斷核范數(shù)方法具有更好的隨機噪聲壓制能力;與經(jīng)典的Curvelet變換和f-x反褶積方法相比,本文方法優(yōu)點為:SP-TNNR方法去噪后的數(shù)據(jù)信噪比更高,去噪后的地震數(shù)據(jù)紋理細節(jié)更加完整,有效信號充分保留。

仿真和實際地震數(shù)據(jù)實驗結果均表明本文提出的SP-TNNR方法相較于其它方法去噪效果更好,具有更高的信噪比。但本文方法也存在著不足,相似塊搜索耗時較長,導致本文方法計算效率較低?？焖俚乃阉魉惴ㄒ约跋嗨茐K匹配中參數(shù)的科學選擇方法將是未來的研究方向;另外,TNNR算法能更加逼近秩,但求解過程中需要兩層SVD分解,計算復雜度較高,如何提高計算效率也是下一步的工作。

致謝:本文在完成過程中采用了哈爾濱工業(yè)大學數(shù)學學院馬堅偉教授團隊的MathGeo2018工具包,在此表示感謝!