基于Adam優(yōu)化的二維自適應(yīng)TFPF地震去噪算法

孟繁磊,范秦寅,穆麗紅

(1.長(zhǎng)春大學(xué) 電子信息工程學(xué)院,長(zhǎng)春 130022；2.大阪大學(xué) 機(jī)械科學(xué)部門,吹田 564-0053， 日本；3.吉林省計(jì)量科學(xué)研究院 幾何量室,長(zhǎng)春 130103)

0 引 言

地震勘探是油氣田礦等資源勘探的主要手段。隨著地震勘探開發(fā)工作的深入,易探易采資源的減少,勘探目標(biāo)向更深、更薄、更不規(guī)則的復(fù)雜地區(qū)發(fā)展的同時(shí),對(duì)勘探的精度及保真度要求越來越高,從事勘探基礎(chǔ)理論研究與應(yīng)用的難度越來越大。一方面對(duì)于淺薄儲(chǔ)集層(如新疆準(zhǔn)噶爾盆地春風(fēng)油田淺薄儲(chǔ)集層)[1]、地層巖性油氣藏(如塔里木盆地碳酸鹽巖油氣),現(xiàn)有的地震資料不能滿足精細(xì)刻畫縫洞體系、準(zhǔn)確預(yù)測(cè)油氣富集區(qū)域的開發(fā)需求[2],因此需要更高精度的新數(shù)據(jù)處理方法。另一方面,無論是為了更好地理解形成和控制陸內(nèi)成礦的深部動(dòng)力學(xué)過程,還是深部找尋資源的現(xiàn)實(shí)需求,勘查走向深部已成為發(fā)展趨勢(shì)[3]。如塔里木盆地、青藏高原腹地、四川盆地、華北地區(qū)及松遼盆地的深層勘探,都需要地震勘探的新理論及交叉學(xué)科的新研究成果解決更多的工程技術(shù)難題。上述全部勘探的基礎(chǔ)是從品質(zhì)好的資料出發(fā)開展后續(xù)解釋工作。壓制地震勘探強(qiáng)隨機(jī)噪聲提高信噪比和分辨率是改善地震記錄品質(zhì)尤為重要的環(huán)節(jié),也是地震信號(hào)研究與處理的難點(diǎn)與熱點(diǎn)之一。

近年來,國(guó)內(nèi)外學(xué)者就消減強(qiáng)隨機(jī)噪聲提高地震勘探資料信噪比問題進(jìn)行了不斷探索,并已取得了實(shí)際意義的成果。如可控方向?yàn)V波器[4]、稀疏去噪[5]、粒子濾波[6]、變分模態(tài)分解[7]、Shearlet[8]以及深度卷積網(wǎng)絡(luò)[9]等方法?；诩哟癢igner-Ville分布(WVD：Wigner-Ville Distribution)的時(shí)頻峰值濾波技術(shù)[10]TFPF(Time-Frequency Peak Filter)能處理極低信噪比情況下(-7 dB)的非平穩(wěn)信號(hào),因此非常適合處理地震數(shù)據(jù)。對(duì)于所有頻率成分,傳統(tǒng) TFPF僅使用固定窗長(zhǎng)進(jìn)行濾波。短窗長(zhǎng)不能有效壓制隨機(jī)噪聲,長(zhǎng)窗長(zhǎng)會(huì)造成有效信號(hào)的衰減?；谶@個(gè)問題,文獻(xiàn)[11-12]提出了徑向道以及自適應(yīng)徑向道TFPF,選擇了更靠近反射同相軸的徑向道方向,而非一維算法沿地震道的方向進(jìn)行TFPF；Tian等[13]提出了變離心率雙曲Trace變換的TFPF；2015年,Zeng等[14]證明TFPF算法可以近似等效于一個(gè)時(shí)不變低通濾波器,并結(jié)合凸優(yōu)化的基本思想,提出了基于濾波器凸集的自適應(yīng)濾波算法,最優(yōu)化求解采用維特比算法。然而,該文獻(xiàn)只是在一維情況下的討論,二維信號(hào)的時(shí)空關(guān)聯(lián)性有可能被忽視甚至被破壞,同時(shí),維特比算法的速度較慢。隨后,曾謙[15]提出了陣列信號(hào)時(shí)空自適應(yīng)濾波算法。在一維自適應(yīng)濾波框架的基礎(chǔ)上選取方向?qū)?shù)作為待優(yōu)化函數(shù)的懲罰項(xiàng),可以兼顧信號(hào)的方向性和時(shí)空相關(guān)性。優(yōu)化函數(shù)的求解使用了傳統(tǒng)的梯度下降法,在每次迭代運(yùn)算中,梯度下降根據(jù)自變量當(dāng)前位置,沿著當(dāng)前位置的梯度更新自變量。然而,自變量的迭代方向僅僅取決于自變量當(dāng)前位置,并且目標(biāo)函數(shù)自變量每個(gè)元素的學(xué)習(xí)率在迭代過程中保持不變,這對(duì)收斂速度與精度都會(huì)造成影響。

筆者在基于凸集構(gòu)造的二維自適應(yīng)濾波思想基礎(chǔ)上,結(jié)合線性等效后的時(shí)頻峰值濾波,提出了二維自適應(yīng)TFPF算法。使用Adam算法求取自適應(yīng)框架中目標(biāo)函數(shù)的全局最優(yōu)解,Adam算法在傳統(tǒng)梯度下降的基礎(chǔ)上使用了動(dòng)量變量和按元素平方的指數(shù)加權(quán)移動(dòng)平均變量,同時(shí),目標(biāo)函數(shù)自變量中每個(gè)元素都分別擁有自己的學(xué)習(xí)率,這可以較梯度下降法更快速穩(wěn)定地逼近最優(yōu)解。

1 二維自適應(yīng)TFPF

1.1 TFPF原理

TFPF是一種基于瞬時(shí)頻率(IF：Instantaneous Frequency)估計(jì)的非線性信號(hào)平滑算法。利用加窗的WVD能量聚集性好的優(yōu)勢(shì),估計(jì)出經(jīng)過調(diào)頻的含噪信號(hào)的瞬時(shí)頻率進(jìn)行去噪。

假設(shè)一組理想觀測(cè)地震信號(hào)

s=x+g

(1)

其中s∈n×1是含噪觀測(cè)信號(hào),x∈n×1是不含噪的有效信號(hào),g∈n×1是隨機(jī)噪聲。值得注意的是,TFPF是一維去噪算法,需要對(duì)二維含噪地震記錄進(jìn)行逐道濾波。

首先,將s調(diào)制成解析信號(hào)z

(2)

其中μ>0是縮放系數(shù)。顯然,s是z的瞬時(shí)頻率。然后,計(jì)算關(guān)于z的加窗WVD

(3)

其中w(wτ向量形式)為窗函數(shù),且長(zhǎng)度為2L+1。最后,從Wi,f中得到峰值處的瞬時(shí)頻率即為x估計(jì)值

(4)

(5)

其中*表示卷積運(yùn)算,h為TFPF的沖激響應(yīng),且和窗函數(shù)w具有相同的長(zhǎng)度。經(jīng)證明,h是一個(gè)時(shí)不變低通濾波器。由h與w之間的映射關(guān)系易知TFPF的濾波效果被窗函數(shù)唯一決定。

1.2 二維自適應(yīng)濾波框架

下面將構(gòu)造濾波器輸出的凸集,并在其內(nèi)設(shè)計(jì)一個(gè)含有時(shí)空相關(guān)性的二次凸函數(shù),在該函數(shù)關(guān)于某個(gè)指標(biāo)取得最小值時(shí)所對(duì)應(yīng)的濾波器輸出即為最優(yōu)解。這樣,二維自適應(yīng)TFPF等效于帶有箱型約束的凸優(yōu)化問題,并存在全局唯一的最優(yōu)估計(jì)。

現(xiàn)選取由不同窗長(zhǎng)得到的一組沖激響應(yīng)集合,H={h1,h2,…,hK(K≥2)},則H的凸包為

(6)

其中conv(H)是包含H的最小凸集,hc∈conv(H)為不同窗長(zhǎng)參數(shù)下TFPF沖激響應(yīng)的凸組合,合理地選取不同的窗長(zhǎng)參數(shù)使conv(H)可以完全而連續(xù)地覆蓋信號(hào)的頻譜范圍。根據(jù)TFPF輸出與沖激響應(yīng)的線性卷積關(guān)系能方便地得到一組濾波器輸出的凸集{Y1,Y2,…,YK},注意Yi∈n×m是經(jīng)線性TFPF逐道濾波后的二維數(shù)據(jù)。這里使用矩陣向量化(vectorization,vec(·))的表示方式,將n×m的二維地震數(shù)據(jù)的列向量首尾相接轉(zhuǎn)換為nm×1的列向量,所以輸出的凸集可寫為Y={vec(Y1),vec(Y2),…,vec(YK)}∈nm×K,同時(shí)由式(6)可知,利用hc逐道卷積二維含噪地震數(shù)據(jù)S并做向量化處理可得

y=vec{hc*S}∈conv(Y)

(7)

這樣,在每個(gè)新的采樣點(diǎn)i=1,2,…,nm上,子集yi是連續(xù)有界的,即

yi∈Ωi=[li,ui]=[minyk,i,maxyk,i]=conv(Yi)

(8)

設(shè)置Ω=Ω1×Ω2×…×Ωnm=conv(Y),便得到經(jīng)向量化后的凸包。進(jìn)而,二維自適應(yīng)TFPF后的所有可能估計(jì)值y均在可行域Ω={y:l≤y≤u}中,其中u=sup{Y}和l=inf{Y}分別為Y的上界和下界。至此,選擇最佳的估計(jì)值y∈conv(Y)等價(jià)于選擇hc∈conv(H),根據(jù)希爾伯特投影理論,在可行域Ω上設(shè)計(jì)一個(gè)關(guān)于濾波輸出的凸函數(shù)J(y),必有全局最小值。可行域Ω的建立過程如圖1所示。

圖1 最優(yōu)解建立凸集過程

1.3 目標(biāo)函數(shù)

在建立了濾波器輸出的某個(gè)凸集Ω后,根據(jù)信號(hào)特點(diǎn),需建立自適應(yīng)濾波的目標(biāo)函數(shù)(凸函數(shù))。假設(shè)不含噪信號(hào)足夠光滑,且至少存在一階絕對(duì)連續(xù)的導(dǎo)數(shù)。因此,將地震勘探隨機(jī)噪聲消減問題表述為凸優(yōu)化問題,也就是尋找一個(gè)定義在Ω上的關(guān)于y的目標(biāo)函數(shù)并使之最小化,擬建立帶懲罰項(xiàng)的最小二乘準(zhǔn)則(PLS：Penalized Least Squares)作為該問題的目標(biāo)函數(shù)[11]

(9)

注意到式(9)采用了矩陣向量化(vec)表示。式(9)第1項(xiàng)是最簡(jiǎn)單的范數(shù)逼近問題----加權(quán)的最小二乘逼近,這項(xiàng)表示濾波結(jié)果y對(duì)含噪信號(hào)s的擬合度,強(qiáng)調(diào)保幅。其中A∈nm×nm是非負(fù)對(duì)角權(quán)矩陣,表示內(nèi)積運(yùn)算。另外,之前假設(shè)不含噪信號(hào)足夠光滑,至少存在一階絕對(duì)連續(xù)的導(dǎo)數(shù),則在最小二乘的基礎(chǔ)上加入了罰函數(shù)(式(9)中第2項(xiàng)與第3項(xiàng)),其中λ稱作平滑因子；D為差分算子,Dr和Dc分別執(zhí)行對(duì)Y的行和列的差分操作,相當(dāng)于對(duì)采樣前的連續(xù)二維信號(hào)Y在兩個(gè)坐標(biāo)方向上的偏導(dǎo)數(shù)。也就是說：Dcy=vec(Dc0Y),Dry=vec(Dr0Y),其中

(10)

參數(shù)|p|=1,可將其理解為cosθ或sinθ,所以第2項(xiàng)與第3項(xiàng)表示了沿著arccosp或arcsinp方向的方向?qū)?shù)的均方值,強(qiáng)調(diào)這一方向的光滑度,也就是強(qiáng)調(diào)去噪。這樣,該目標(biāo)函數(shù)引入了時(shí)空關(guān)聯(lián)性。

采用加權(quán)的PLS為信號(hào)在每個(gè)采樣點(diǎn)選擇符合約束的濾波器。如果權(quán)值是“自適應(yīng)”的,則達(dá)成了在每個(gè)采樣點(diǎn)自適應(yīng)濾波(選擇濾波器)的初衷。因此,乘在最小二乘上的權(quán)值A(chǔ)和乘在懲罰項(xiàng)上的平滑因子λ,共同決定了濾波器輸出的跟蹤性能與平滑性能的折衷。對(duì)A和λ的賦值方法如下

Aa,a=|ua-la-ψ|,ψ=σ(max‖Hk‖F(xiàn)-min‖Hk‖F(xiàn)),λ=(κψ)2,κ>0

(11)

(12)

當(dāng)某時(shí)刻權(quán)值很大時(shí),濾波器會(huì)傾向于選擇較多的擬合度,反之則選擇較多的平滑度,使濾波估計(jì)可以在時(shí)空域某方向上達(dá)到最佳平滑效果。

1.4 最優(yōu)化問題的求解

目標(biāo)函數(shù)J(y)：nm→是一個(gè)實(shí)值函數(shù)。筆者討論的優(yōu)化問題的含義是指能使函數(shù)J達(dá)到最小時(shí),從約束凸集Ω中找出與之相對(duì)應(yīng)的向量y,即arg minJ(y)。由前述可知,在凸集Ω上定義的凸函數(shù)有唯一最優(yōu)解。

筆者采用投影Adam算法求解,Adam是基于梯度下降的優(yōu)化算法[16]

y(k+1)=y(k)-αJ(y(k))

(13)

圖2 以幾何方式演示利用梯度法求二元單值函數(shù)J：2→的極小值點(diǎn)

但給定了α,由于初始點(diǎn)y(0)中的各個(gè)元素位置不同,則在梯度下降迭代y中的nm×1個(gè)元素時(shí),可能導(dǎo)致某些元素越過目標(biāo)函數(shù)最優(yōu)解,這就需要一個(gè)較小的α。然而,這也會(huì)造成自變量朝最優(yōu)解移動(dòng)變慢[18]。

Adam使用了如下的動(dòng)量變量v(k)和按梯度元素平方的指數(shù)加權(quán)移動(dòng)平均變量s(k),并將它們中每個(gè)元素在初始時(shí)刻時(shí)設(shè)置為0[19]。

v(k)=β1v(k-1)+(1-β1)g(k)

(14)

s(k)=β2v(k-1)+(1-β2)g(k)?g(k)

(15)

(16)

(17)

y(k+1)=y(k)-g′(k)

(18)

對(duì)于式(18)的迭代算法,當(dāng)梯度為零時(shí)可作為在理論上的停止規(guī)則。但在實(shí)際應(yīng)用中,很難得到g′(k)=0的一步迭代,因此可以計(jì)算相鄰兩個(gè)迭代點(diǎn)差值的范數(shù)‖y(k+1)-y(k)‖,如果小于預(yù)設(shè)的閾值(如10-5),則停止迭代。Adam優(yōu)化算法的基本流程如圖3所示。

圖3 Adam優(yōu)化算法的基本流程

(19)

點(diǎn)PΩ[y]稱為y到Ω的投影,PΩ為投影算子?？梢钥闯?PΩ是Ω中最接近y的點(diǎn)。

2 仿真結(jié)果與分析

2.1 人工合成地震記錄

首先給出一個(gè)復(fù)雜人工合成地震記錄測(cè)試二維TFPF的去噪能力。如圖4a所示,純凈信號(hào)包含5個(gè)由雷克子波形成的反射同相軸,其形狀有直線軸、雙曲軸、交叉軸與斷軸,其頻率分別為45 Hz,35 Hz,30 Hz,25 Hz和20 Hz。共56道,每道信號(hào)有1 500個(gè)采樣點(diǎn),采樣頻率為 1 000 Hz。為了更好地模擬真實(shí)情況,在純凈信號(hào)中添加了實(shí)際的地震初至前噪聲,使信噪比為-4.766 1 dB,如圖4b所示。從圖4b中可以看到,有效信號(hào)幾乎被完全淹沒。分別采用傳統(tǒng)TFPF、一維和二維自適應(yīng)TFPF處理含噪信號(hào),傳統(tǒng)TFPF的窗長(zhǎng)設(shè)為9點(diǎn),后兩種方法采用的TFPF濾波器集合參數(shù)確保得到一個(gè)合理的對(duì)比：H={h1,h1 499},對(duì)應(yīng)于矩形窗分別為1和1 499點(diǎn)數(shù)據(jù)長(zhǎng)度,該組濾波器的通帶足以覆蓋全部信號(hào)頻譜。其中二維TFPF的方向?qū)?shù)參數(shù)p從0～1取值時(shí),輸出信噪比的變化情況如圖4c所示,可以看出在p=0.1時(shí)取得最大輸出信噪比。濾波結(jié)果分別如圖4d～圖4f所示。對(duì)比3種算法可以看出,TFPF在強(qiáng)實(shí)際地震噪聲背景下的去噪效果并不明顯,仍然有大量的隨機(jī)噪聲殘留；兩個(gè)自適應(yīng)TFPF均較好地恢復(fù)出有效同相軸,但二維算法去噪效果更明顯,背景更干凈。

圖4 復(fù)雜人工合成地震記錄例子

另外二維方法使用Adam優(yōu)化,圖5a給出了學(xué)習(xí)率為0.005時(shí)梯度下降和Adam算法的迭代次數(shù)與損失函數(shù)值的關(guān)系,從中可以看出Adam比梯度下降算法有更快的收斂速度。在圖5b中,當(dāng)學(xué)習(xí)率設(shè)置為0.1時(shí),梯度下降算法由于學(xué)習(xí)率過大導(dǎo)致梯度彌散,而Adam只需20次迭代即可穩(wěn)定。

a 學(xué)習(xí)率為0.005 b 學(xué)習(xí)率為0.1

接著,為了更加深入地對(duì)比,對(duì)上述純凈、含噪、傳統(tǒng)TFPF、一維與二維自適應(yīng)TFPF 5個(gè)記錄分別比較了F-K譜圖,如圖6a～圖6e所示,可以明顯地看出,二維自適應(yīng)TFPF 的F-K譜比一維的方法更接近純凈信號(hào),這說明二維自適應(yīng)方法在復(fù)雜地形與強(qiáng)噪聲背景下依然可以有效恢復(fù)信號(hào)。

圖6 各記錄F-K譜圖

最后利用

(20)

分別計(jì)算不同記錄的輸出信噪比定量分析實(shí)驗(yàn)結(jié)果,式(20)中x0是無噪信號(hào),x是濾波后的信號(hào),M和N是矩陣的維度。實(shí)驗(yàn)分析結(jié)果列入表1中,由表1可看出,二維自適應(yīng)TFPF的輸出信噪比最高。

表1 各記錄的SNR

2.2 實(shí)際地震記錄

圖7 實(shí)際地震記錄對(duì)比

筆者通過處理實(shí)際地震資料驗(yàn)證了二維自適應(yīng)TFPF方法壓制地震勘探隨機(jī)噪聲的適用性。選取兩幅實(shí)際地震記錄。圖7a給出了一個(gè)168道的共炮點(diǎn)地震記錄,每道包含2 025個(gè)采樣點(diǎn),記錄的采樣頻率為250 Hz。從圖7a中可以發(fā)現(xiàn),記錄中同相軸的形狀十分不規(guī)則,很多有效軸的清晰度和連續(xù)性都被隨機(jī)噪聲所干擾。采用一維和二維自適應(yīng)TFPF分別處理這兩幅地震資料,其濾波器集合參數(shù)均為H={h1,h2 025},二維方法使用Adam優(yōu)化,初始學(xué)習(xí)率為0.1。結(jié)果分別如圖 7b、圖7c所示。其中一維方法處理后的記錄,隨機(jī)噪聲并沒有完全去除；而二維方法的噪聲壓制得更加徹底,幾乎沒有殘留,同相軸清晰可見。有代表性的區(qū)域已用白色矩形框標(biāo)出。

3 結(jié) 語

筆者針對(duì)一維自適應(yīng)TFPF沒用利用空間相關(guān)信息的問題,研究了基于濾波器凸集的二維自適應(yīng)TFPF算法。在一維自適應(yīng)目標(biāo)函數(shù)框架下加入了方向?qū)?shù)參數(shù),同時(shí)使用投影的Adam算法對(duì)目標(biāo)函數(shù)進(jìn)行優(yōu)化,從而實(shí)現(xiàn)了快速穩(wěn)定的二維濾波。通過人工合成實(shí)驗(yàn)和實(shí)際資料處理結(jié)果表明,筆者提出的算法充分利用了地震同相軸的橫向相關(guān)性,能更清晰地恢復(fù)出有效同相軸,并且對(duì)噪聲有較好的壓制作用。