基于Adam優化的二維自適應TFPF地震去噪算法

孟繁磊,范秦寅,穆麗紅

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

0 引 言

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

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

筆者在基于凸集構造的二維自適應濾波思想基礎上,結合線性等效后的時頻峰值濾波,提出了二維自適應TFPF算法。使用Adam算法求取自適應框架中目標函數的全局最優解,Adam算法在傳統梯度下降的基礎上使用了動量變量和按元素平方的指數加權移動平均變量,同時,目標函數自變量中每個元素都分別擁有自己的學習率,這可以較梯度下降法更快速穩定地逼近最優解。

1 二維自適應TFPF

1.1 TFPF原理

TFPF是一種基于瞬時頻率(IF：Instantaneous Frequency)估計的非線性信號平滑算法。利用加窗的WVD能量聚集性好的優勢,估計出經過調頻的含噪信號的瞬時頻率進行去噪。

假設一組理想觀測地震信號

s=x+g

(1)

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

首先,將s調制成解析信號z

(2)

其中μ>0是縮放系數。顯然,s是z的瞬時頻率。然后,計算關于z的加窗WVD

(3)

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

(4)

(5)

其中*表示卷積運算,h為TFPF的沖激響應,且和窗函數w具有相同的長度。經證明,h是一個時不變低通濾波器。由h與w之間的映射關系易知TFPF的濾波效果被窗函數唯一決定。

1.2 二維自適應濾波框架

下面將構造濾波器輸出的凸集,并在其內設計一個含有時空相關性的二次凸函數,在該函數關于某個指標取得最小值時所對應的濾波器輸出即為最優解。這樣,二維自適應TFPF等效于帶有箱型約束的凸優化問題,并存在全局唯一的最優估計。

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

(6)

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

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

(7)

這樣,在每個新的采樣點i=1,2,…,nm上,子集yi是連續有界的,即

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

(8)

設置Ω=Ω1×Ω2×…×Ωnm=conv(Y),便得到經向量化后的凸包。進而,二維自適應TFPF后的所有可能估計值y均在可行域Ω={y:l≤y≤u}中,其中u=sup{Y}和l=inf{Y}分別為Y的上界和下界。至此,選擇最佳的估計值y∈conv(Y)等價于選擇hc∈conv(H),根據希爾伯特投影理論,在可行域Ω上設計一個關于濾波輸出的凸函數J(y),必有全局最小值?？尚杏颚傅慕⑦^程如圖1所示。

圖1 最優解建立凸集過程

1.3 目標函數

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

(9)

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

(10)

參數|p|=1,可將其理解為cosθ或sinθ,所以第2項與第3項表示了沿著arccosp或arcsinp方向的方向導數的均方值,強調這一方向的光滑度,也就是強調去噪。這樣,該目標函數引入了時空關聯性。

采用加權的PLS為信號在每個采樣點選擇符合約束的濾波器。如果權值是“自適應”的,則達成了在每個采樣點自適應濾波(選擇濾波器)的初衷。因此,乘在最小二乘上的權值A和乘在懲罰項上的平滑因子λ,共同決定了濾波器輸出的跟蹤性能與平滑性能的折衷。對A和λ的賦值方法如下

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

(11)

(12)

當某時刻權值很大時,濾波器會傾向于選擇較多的擬合度,反之則選擇較多的平滑度,使濾波估計可以在時空域某方向上達到最佳平滑效果。

1.4 最優化問題的求解

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

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

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

(13)

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

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

Adam使用了如下的動量變量v(k)和按梯度元素平方的指數加權移動平均變量s(k),并將它們中每個元素在初始時刻時設置為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)

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

圖3 Adam優化算法的基本流程

(19)

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

2 仿真結果與分析

2.1 人工合成地震記錄

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

圖4 復雜人工合成地震記錄例子

另外二維方法使用Adam優化,圖5a給出了學習率為0.005時梯度下降和Adam算法的迭代次數與損失函數值的關系,從中可以看出Adam比梯度下降算法有更快的收斂速度。在圖5b中,當學習率設置為0.1時,梯度下降算法由于學習率過大導致梯度彌散,而Adam只需20次迭代即可穩定。

a 學習率為0.005 b 學習率為0.1

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

圖6 各記錄F-K譜圖

最后利用

(20)

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

表1 各記錄的SNR

2.2 實際地震記錄

圖7 實際地震記錄對比

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

3 結 語

筆者針對一維自適應TFPF沒用利用空間相關信息的問題,研究了基于濾波器凸集的二維自適應TFPF算法。在一維自適應目標函數框架下加入了方向導數參數,同時使用投影的Adam算法對目標函數進行優化,從而實現了快速穩定的二維濾波。通過人工合成實驗和實際資料處理結果表明,筆者提出的算法充分利用了地震同相軸的橫向相關性,能更清晰地恢復出有效同相軸,并且對噪聲有較好的壓制作用。