基于自相似性和低秩先驗(yàn)的地震數(shù)據(jù)隨機(jī)噪聲壓制

程文婷,方文倩,付麗華

(中國(guó)地質(zhì)大學(xué)(武漢)數(shù)學(xué)與物理學(xué)院,湖北武漢430074)

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

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

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

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

1 方法理論

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

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

M=X+E

(1)

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

(2)

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

秩函數(shù)是非凸非光滑函數(shù),導(dǎo)致公式(2)的求解是一個(gè)NP難問(wèn)題?？衫煤朔稊?shù)近似秩函數(shù),將公式(2)轉(zhuǎn)化為凸優(yōu)化問(wèn)題:

(3)

1.2 基于截?cái)嗪朔稊?shù)的低秩模型

(4)

因?yàn)楣?4)是非凸的,所以先利用下文定理1對(duì)目標(biāo)函數(shù)進(jìn)行轉(zhuǎn)換,然后再利用APGL算法求解。

1.3 目標(biāo)函數(shù)轉(zhuǎn)換

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

(5)

這樣:

(6)

于是:

(7)

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

(8)

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

1.4 相似塊匹配

相似度定義為歐式距離:

(9)

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

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

(10)

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

1.5 APGL優(yōu)化求解

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

(11)

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

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

(12)

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

(13)

對(duì)于目標(biāo)函數(shù)(公式(10)),令:

g(X)=‖Xi‖*

則:

(14)

對(duì)矩陣X∈Rm×n進(jìn)行SVD分解:

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

(15)

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

(16)

由文獻(xiàn)[20]可知,對(duì)任意τ>0,Y∈Rm×n,有:

(17)

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

(18)

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

(19)

(20)

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

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

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

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

(21)

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

2.1 仿真數(shù)據(jù)實(shí)驗(yàn)

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

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

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

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

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

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

2.2 實(shí)際數(shù)據(jù)應(yīng)用

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

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

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

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

圖9 實(shí)際疊后地震數(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 結(jié)論

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

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

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