三維地震數據頻域無監督隨機噪聲壓制方法

薛亞茹，蘇軍利*，馮璐瑜，張程，梁琪

（1.中國石油大學（北京）信息科學與工程學院，北京 102249；2.中國石油大學（北京）油氣資源與工程全國重點實驗室，北京 102249）

0 引言

地震數據往往受自然環境、采集儀器等因素的影響而被不同程度噪聲干擾，從而降低了信噪比和保真度。因此，對采集的地震數據去噪是不可或缺的重要環節。

依據地震數據特征和不同先驗假設提出了各種去噪理論和實現方法，主要可分為以下幾類。①基于預測的濾波方法，如t-x域預測濾波[1]、f-x域反褶積濾波[2]、非平穩性預測濾波[3]等。該類方法基于地震道之間的可預測性，適用于線性同相軸。②基于稀疏變換的濾波方法，如傅里葉變換[4]、小波變換[5]、曲波變換[6]等。此類方法假設地震數據在變換域能被稀疏表征，通過抑制較小的參數值達到去噪目的，在一定程度上可以獲得良好的去噪效果，但變換域噪聲壓制的閾值等超參數難以確定。 ③矩陣降秩方法。該類方法假設地震數據可以表征為一個低秩矩陣，隨機噪聲的存在會增加矩陣的秩，因此可以將去噪問題轉化為低秩矩陣逼近問題。如Cadzow 濾波[7]、奇異值分解[8]、多通道奇異譜分析(Multichannel Singular Spectrum Analysis，MSSA)[9]等方法。字典學習本質也屬于降秩方法，它能夠從地震數據中自適應地學習基函數，從而對數據進行稀疏表示，以便將信號與噪聲進行分離，達到好的去噪效果[10]，然而其計算成本較高，而且算法在字典更新步驟中會破壞稀疏系數的結構。模態分解也是降秩方法之一，通過提取有限的模態分量達到信噪分離的目的，目前使用最為廣泛的有經驗模態分解[11]和變分模態分解[12]。

近年來，隨著計算機性能的迅速提高及大數據時代的來臨，深度學習研究備受青睞。不少學者將深度學習引入到地球物理領域實現對地震數據去噪處理[13-15]。最初的研究多是采用有監督深度學習模型，相較于依賴先驗信息的傳統算法，有監督模型在地震數據去噪和插值方面表現出了優異的性能[16-18]，但是海量高質量標簽數據的獲取，無疑增加了人工和計算成本。因此，無監督學習受到了廣泛關注?；跓o監督學習的自編碼網絡模型通過對輸入數據編碼壓縮和解碼重構達到去噪的目的[19-20]。Zhang 等[21]通過對網絡中間層加入稀疏約束提高網絡去噪性能；Yang等[22]提出了一種全連接型深度去噪無監督網絡用于三維地震數據噪聲壓制；Saad等[23]將無監督網絡用于增強微地震信號和壓制隨機噪聲，通過對原始數據進行切片操作，初步減小輸入數據的維度，之后利用自編碼網絡實現信噪分離。為了強化網絡對地震數據高階特征的學習能力，在編碼器和解碼器之間加入跳躍連接，使得解碼器在獲得全局特征的同時可以融入編碼器提取的局部特征[24]?？紤]到疊前地震資料經動校正(NMO)后具有較高的自相似性，Liu 等[25]提出了一種基于無監督深度學習壓制面波的方法，將經過動校正后的地震數據作為神經網絡的輸入，利用UNet 型卷積網絡提取數據特征，從而達到信噪分離的目的。

上述基于神經網絡的地震數據去噪方法大多在時空域完成。將地震數據剖分成小的時空切片[24]，利用切片之間的相似性[25]，通過神經網絡提取有效信號特征，實現信噪分離。此類方法將數據切分，未曾考慮切片之間的連續性，也未考慮三維數據體之間的相關性。地震數據在頻域是可預測的，MSSA 法利用該特性實現噪聲壓制和數據重構[9]，并能夠較好地保持橫向連續性，但該方法需要構建大型Hankel 矩陣，計算量大；同時基于奇異值分解理論的壓縮去噪方法提取的地震數據特征是相互正交的，對信號進行稀疏表征的程度有限，而且有效信號秩的優選也比較困難。

本文提出一種頻域降秩的深度學習去噪方法。首先分析三維地震數據體的低秩特性，利用奇異值分解理論指導自編碼網絡的建立；針對頻域噪聲分布特征，引入K-L(Kullback-Leibler)散度優化損失函數。合成記錄和實際地震資料處理結果表明，本文方法在隨機噪聲壓制效果和處理效率等方面均優于MSSA和K-SVD(K-奇異值分解)方法。

1 方法原理

1.1 三維地震數據的低秩表征

地震勘探通過多炮激發實現地下構造反射成像。為提高成像精度，震源在探區多次激發，構成了炮點坐標(xs，ys)、檢波點坐標(xr，yr)及時間t的五維地震數據，可表示為d(xs，ys，xr，yr，t)。為方便解釋，本文以一條測線的炮集為例進行分析，但同樣適用其他三維數據體。假設地層是水平層狀介質，第k炮數據可表示為dk(xs，xr，t)，則第k+1炮記錄可表示為

式中：Δxs為炮間距；p為射線參數。變換到頻率域，有

由上式可知，相鄰炮記錄之間的頻譜僅有線性相位差s=exp(-i2πfpΔxs)。對某個頻率而言，各炮數據構成低秩矩陣，可以表示為

式中：n為檢波點個數；m為炮數。由此可以看到該矩陣的秩是1，也就是說各炮記錄是相似的。但是考慮到實際地層并非水平，也并非各向同性，上述矩陣秩不一定是1，但仍滿足低秩假設，因此對于三維地震數據, 其頻率切片滿足低秩特征。

M的奇異值分解可表示為

式中：r為矩陣M的秩；U是矩陣MMT的特征向量矩陣，由各炮f-xr譜特征向量組成，且相互正交；V是矩陣MTM的特征向量矩陣，由各檢波點的f-xs譜特征向量組成；S表示各特征值構成的對角矩陣；σk為第k個特征值；uk、vk為第k個特征值對應的特征向量。

根據上述分析，頻率切片會呈現一定的低秩特征(或稀疏性)，對任一炮的某一頻率分量可以表示為特征向量的稀疏表示，即

式中：mj表示第j炮頻率分量；cj=S(VT)j為稀疏表征參數。該表示方式也可以理解為稀疏字典概念，其中U也稱為字典。通過適當的降秩，可以達到壓制噪聲的目的。

奇異值分解方法獲得的各特征向量是相互正交的，這限制了數據稀疏表征程度。為了解決該問題，提出了獨立分量分析、基于數據驅動的冗余字典等方法。本文將采用數據驅動的神經網絡框架，在頻域實現地震數據壓縮去噪。

1.2 無監督自編碼網絡模型

自編碼器是Lecun[26]提出的一種具有對稱結構的無監督型神經網絡。對于輸入數據為X，編碼器函數可表示為

該過程期望獲得數據X的稀疏編碼。式中：h為隱層神經元的輸出；w為網絡輸入的權重；b為網絡輸入的偏置。與式(5)相比，實質就是尋找稀疏編碼系數c；編碼器各神經元的權重系數w近似為式(5)中的特征向量矩陣U的轉置，但沒有U各特征向量的正交性限制。自編碼網絡可以通過神經網絡的訓練過程優化各特征向量，打破正交約束，大大提高數據特征的稀疏性。

然后，網絡通過將來自隱層的低維空間特征進行解壓縮達到數據重構和還原的目的，這一過程也被稱為解碼。解碼器函數可以表示為

該過程即為式(6)的逆過程，利用編碼器的稀疏系數重構數據。式中：Y為神經網絡的輸出；w′為網絡輸出的權重；b′為網絡輸出的偏置。在整個過程中由于壓縮作用，會將不相干的噪聲壓制，而只保留相干性較強的有效信號。

通過上述分析可見，自編碼網絡工作原理與奇異值分解相似，這為自編碼網絡提供了工作機理，克服了神經網絡的黑箱特性。圖1解釋了自編碼網絡與奇異值分解在結構和原理上的對應關系。

圖1 自編碼網絡與奇異值分解原理映射圖

基于上述自編碼網絡去噪原理，設計了無監督網絡模型(圖2)。首先，應用快速傅里葉變換(FFT)將三維數據變換到頻率域，由于頻率域數據是復數，將某一切片數據實部和虛部串接作為神經網絡的輸入。然后，采用自編碼網絡實現噪聲壓制，將處理過的所有頻率切片重構回三維數據頻譜，進行快速傅里葉逆變換(IFFT)得到去噪后數據。

圖2 無監督去噪模型框架

在該處理流程中，自編碼網絡結構是影響去噪的關鍵因素。一是網絡不同層之間的連接方式。本文在隱層添加跳躍連接塊以增強網絡對地震數據特征的學習，如此，解碼層的輸入包含兩部分，即來自隱層的輸出和來自跳躍連接塊的輸出。跳躍連接塊的使用可以加速無監督模型收斂，增強網絡提取高階特征的能力。解碼器重建抽象特征，然后將這些抽象特征映射回原始地震數據空間，以實現隨機噪聲衰減。通過迭代學習抽象特征表示，盡可能減少編碼器和解碼器的差異。

二是網絡中間層的維度決定了數據壓縮程度。維度太小，將導致信號壓縮過度，損失有效信號；維度太大，壓縮不到位，噪聲去除效果不佳。采用深度學習的目的是希望獲得比SVD 方法更好的壓縮編碼效果，因此利用SVD 分解后特征值能量占比為90%時對應的階數，作為中間層維度的參考。此策略保證自編碼網絡的稀疏度低于SVD 方法，能學習到更稀疏的數據特征。

1.3 基于K-L 散度的優化模型

均方誤差函數是回歸問題中較為常見的損失函數。根據Bayes 反演理論，最小化均方誤差函數是以零均值高斯白噪聲為假設條件。但是由于本文方法是在頻率域完成，白噪聲的存在導致頻率切片的均值偏移，因此對去噪結果進行均值約束，有利于數據歸位于其真值。K-L 散度可以度量兩個分布之間的距離，定義為

式中：Pde=mean(Y)，表示對去噪后數據Y取均值；Pref為參考均值，理論上應該等于干凈信號的均值。式（8）度量了去噪數據均值與參考均值的距離，當Pde=Pref時，該值達到最小。

將上述K-L 散度約束與均方誤差函數相結合，定義新的損失函數為

正則化項的加入可以驅動去噪結果均值接近參考值。

2 實驗結果分析

使用一組合成數據和兩組海洋實際數據驗證本文方法的有效性，并根據信噪比(Signal-Noise Ratio，SNR)、局部相似性對方法進行評價。

區域(包括本礦區)開展過1∶20萬、1∶5萬水系沉積物地球化學測量以及1∶1萬土壤地球化學測量工作。1∶5萬水系沉積物測量成果顯示，本區元素異常以Mo、W為主，二者強度高、規模大,且與1∶20萬水系沉積物地球化學鉬異?？傮w形態相似，套合較好。1∶1萬土壤地球化學測量Mo、W、Ag、Cu異常強度高、規模較大，是礦區的主要成礦元素和伴生元素。Mo異常集中分布在上棋盤北側和南側，以及大尖山一帶。上棋盤北側和南側的Mo異常，總體展布方向呈北北西—南南東向，大尖山一帶的Mo異常，總體展布方向呈北西—南東向，與區域上Mo異常成礦帶基本一致。經普查地質工作，化探鉬異常多為礦致異常。

SNR作為全局衡量標準，定義為

2.1 合成數據

使用褶積模型合成三維數據體。該模型共128×128 道，每道256 個采樣點，時間采樣間隔為2 ms。對合成數據添加方差為0.3 的高斯噪聲，信噪比低至-5.96 dB(圖3a)，分別采用MSSA、K-SVD 和本文方法去噪。K-SVD 方法去噪效果最差，信噪比僅提升至2.72dB，同相軸附近殘留噪聲較多，殘差道集也可以看出同相軸附近痕跡明顯，說明該方法去除噪聲并不徹底(圖3b)；MSSA 方法較K-SVD 方法去噪效果有所改善，信噪比提高至4.52 dB，但殘差道集有同相軸殘留(圖3c)，說明該方法對有效信號造成了一定程度的損傷；本文方法去噪較為徹底、效果最好，信噪比提升至6.55 dB，同相軸刻畫較為清晰，殘差數據也基本看不到同相軸痕跡(圖3d)，說明本文方法在壓制隨機噪聲的同時不會對有效信號造成損傷。

圖3 合成數據三種方法去噪效果三維對比

為方便觀察各方法的數據處理細節，抽取該三維數據Inline 方向第11 個剖面進行對比(圖4)。與原始合成數據(圖4a 左)相比，加噪數據的同相軸被淹沒在強噪聲中，信噪比低至-5.98 dB(圖4a 右)。K-SVD方法去噪結果的信噪比僅提高到2.47 dB(圖4b 左)，殘差剖面中同相軸痕跡明顯(圖4b右)。MSSA 方法去噪結果的信噪比達到4.37 dB，但同相軸處仍有明顯的噪聲殘留(圖4c 左)，殘差剖面中也留有同相軸痕跡(圖4c右)。本文方法去噪后，同相軸更連續，信噪比提高到6.37 dB(圖4d 左)，殘差剖面中同相軸痕跡不明顯(圖4d 右)。由此可以看出，本文方法壓制噪聲效果優于其他兩種方法。

圖4 合成數據三種方法去噪效果剖面對比

進一步抽取主頻為30 Hz 的頻率剖面，測試重構階數或中間層維度對SVD 算法和無監督網絡去噪結果的影響，即比較編碼稀疏度對去噪結果的影響(圖5)。當中間層維度取19 時，無監督網絡去噪結果的SNR 最高、效果最好(圖5中綠色虛直線②所示)，恰是該頻率切片奇異值能量占比為90%所對應的維度。從圖中不難發現，無監督網絡去噪結果的SNR 一直高于SVD 方法，即當兩者有相同的稀疏度時，無監督網絡性能更優(圖5 中綠色虛直線①所示)。即使在更小的稀疏參數下，無監督網絡也可以獲得比SVD 方法更佳的去噪效果。

圖5 中間層維度對去噪性能的影響

提取出無監督網絡第二編碼層的特征向量，并計算其協方差矩陣，即各特征向量之間的相關譜(圖6)，除了主對角線上，其余位置仍然有不少非零值，可見其特征向量并未完全正交，即采用無監督網絡打破了傳統SVD 方法特征向量正交這一約束，有利于提高數據表征的稀疏性。

圖6 第二編碼層特征向量相關譜

去噪方法的時間成本同樣是一個值得關注的問題。對于合成三維數據體，MSSA、K-SVD 和本文方法處理時間分別為792、1358、276 s。由于K-SVD 方法每次迭代都需要奇異值分解運算，耗時最大；MSSA 也需要大規模矩陣的奇異值分解，耗時也較大；而本文方法采用的數據結構簡單，自編碼網絡通過反向傳播優化網絡，計算簡單，可以快速實現。

2.2 實際數據一

為了定量評價本文方法的實際數據處理能力，對實際海洋地震數據進行去噪實驗。該數據共有400炮，每一炮包含400道數據，每道數據有800個時間采樣點，采樣間隔為2 ms。圖7a是原始地震道集及其加入高斯白噪聲的數據，信噪比降至-7.30 dB，原始數據中弱信號同相軸被噪聲淹沒。圖7b～圖7d 分別是K-SVD、MSSA 和本文方法的處理結果及其殘差道集，可以看出本文方法的殘差道集中有效信號能量最弱。對圖7 紅框區域進行放大顯示，如圖8所示?？梢悦黠@看出，MSSA 和K-SVD 方法去噪后的數據同相軸附近仍殘留大量噪聲，有效信號能量較弱，而本文方法去噪后同相軸最清晰。K-SVD、MSSA 和本文方法處理結果的SNR 分別為2.76、4.73 和5.91 dB，顯然本文方法的去噪效果最好。

圖7 實際海洋數據三種方法去噪效果對比

圖8 圖7 紅框區域的放大顯示

為了驗證本文方法保持數據橫向連續性的優勢，抽取第10 個共炮檢距道集0.4～1.0 s 部分(圖9)。由圖9 可以看出，本文方法在強噪聲背景下仍較好地恢復了同相軸連續性，而且最清晰，尤其是在圖中箭頭處，效果明顯優于其他兩種方法。

圖9 實際海洋數據三種方法去噪結果的共炮檢距道集對比

對于實際數據一，MSSA、K-SVD 和本文方法處理時間分別為6524、1315、1160s。由于實際數據維度較大，因此MSSA 的計算效率低于K-SVD?？梢?，無論是合成數據還是實際數據，本文方法的計算效率最高。

本文方法中參數Pref的取值會影響各個子網絡對頻率切片的去噪性能。選取不同的Pref值對30 Hz含噪頻率切片數據進行去噪處理并計算其SNR(圖10)。含噪切片數據的SNR 為-6.59dB，當Pref取值0.03 時，去噪后SNR 最高，這也恰是無噪數據的均值。隨著Pref逐漸偏離真值，去噪后SNR 也會隨之下降。由此可以看出，K-L 散度約束對去噪結果有一定影響。

圖10 參數Pref對去噪效果的影響

2.3 實際數據二

為了評價本文方法的實用性，對實際含噪數據二進行處理。實際數據一中的隨機噪聲是人為添加，而實際數據二中是真實噪聲。圖11a 是原始含噪道集，噪聲導致同相軸的清晰度降低；圖11b 為K-SVD 方法去噪結果及殘差道集，可以看出，噪聲去除并不徹底，尤其是遠炮檢距仍然可以看到大量噪聲存在；圖11c 為MSSA 方法去噪結果及殘差剖面，噪聲得到了有效壓制，但是同相軸的連續性并未得到恢復，并且殘差道集中有明顯的有效信號痕跡；圖11d是本文方法去噪結果及殘差剖面，噪聲壓制比較徹底、同相軸更清晰，殘差道集中有效信號不明顯。需承認的是，本文方法在去噪的同時對弱信號會造成輕微的損傷，如圖中4.8 s 以后弱信號難以恢復。這是由于在頻率域，淺層強信號和深層弱信號耦合，經過自編碼網絡的壓縮會壓制弱信號。

圖11 實際地震數據二三種方法去噪效果對比

去噪結果與去除的噪聲之間的局部相似性能評價信號泄漏程度[27]，高相似異常就表明去噪結果在相應的位置有信號泄漏。圖12 展示了K-SVD、MSSA和本文方法去噪結果與去除的噪聲之間的相似性。由圖可見：K-SVD 方法處理后在2.4～3.2 s 大炮檢距處，殘差信號與去噪結果存在較高的相似性，說明去噪結果中仍然保留了部分噪聲，壓制不徹底；經過MSSA 方法去噪后的數據在1.8～3.4 s除了大炮檢距處噪聲殘留，對有效信號損傷也較大，這在圖11c的殘差道集上也很明顯；而本文方法去除的噪聲與保留的信號相干性較小，說明在去噪時沒有損傷有效信號，具有較高的保真度。

圖12 三種方法去噪結果與殘差的局部相似性對比

3 結論

本文利用三維數據頻域相似性，建立了一種奇異值分解理論指導的自編碼無監督去噪網絡模型。通過將奇異值分解映射為自編碼網絡，克服了奇異值分解各特征向量正交的約束，提高了數據表示的稀疏性。根據頻率域噪聲分布特點在損失函數中引入K-L散度約束，驅動數據恢復其均值分布。本文方法實現了一種運算速度快、又保留傳統方法物理可解釋性的無監督學習算法模型。合成和實際數據處理結果均表明，本文方法在去噪的同時能夠保留信號細節，具有較高的保真度。