時頻域變分模態分解地震資料去噪方法

胡瑞卿 何俊杰 李華飛 張曉莉 裴家定 劉億偉

(①東方地球物理公司研究院庫爾勒分院，新疆庫爾勒 841000；②西安石油大學地球科學與工程學院，陜西西安 710065)

0 引言

目前地震資料隨機噪聲常規壓制方法主要包括空間域壓制和變換域壓制。前者有中值濾波、各向異性擴散濾波等；后者有小波閾值法、曲波變換閾值法等。其本質都是利用在特定域中有效信號與噪聲具有可區分的特征進行信噪分離。

時頻分析類方法主要根據時頻域中有效信號與噪聲的能量分布，進行噪聲的分析和定位及剔除或壓制[1-3]，主要有基于短時傅里葉變換、連續小波變換和S變換等方法[4]。Morlet[5]提出的改進小波變換可在時頻域中尋找目標信號的局部特征。但基于小波變換的一系列延伸方法對具有復雜同相軸的地震資料，效果并不理想。針對該問題，Candès等[6]提出了具有多方向、多尺度的曲波變換。包乾宗等[7]將曲波變換與奇異值分解方法相結合，對地震資料中的不同成分進行分離；Neelamani等[8]將曲波變換用于三維地震數據的隨機噪聲的壓制；孟閣閣等[9]基于曲波變換，設計了一種自適應閾值去噪方法，效果比較明顯；張恒磊等[10-11]采用非線性閾值法基于曲波變換對疊前資料進行了去噪處理，隨后提出了適用于曲波域和空間域的傾角掃描去噪技術，并在復雜區塊資料中取得了較好的效果；張之涵等[12]提出一種基于三維曲波變換的噪聲壓制方法，通過改進非線性閾值，更好地保護了有效信號；姚振岸等[13]將曲波變換與各向異性擴散濾波相結合，改善了常規曲波變換降噪方法的過度平滑與環繞效應，實現了數據邊界和有效信號振幅的保護。

近幾年，諸如奇異值分解、稀疏表征等新方法也逐步應用于地震資料的噪聲壓制。它們通過構建適用于目標資料的字典，替代小波變換、曲波變換中的確定性基函數，其本質是將去噪的過程轉化為優化問題。基于該思路，Sanchis等[14]通過獨立分量分析算法構建基函數，但該方法需要有效成分的統計估計；李海山等[15]采用形態分量分析方法進行字典設計，通過求解表征目標函數達到去噪的效果，但測試結果顯示仍存在一定的偽影現象；Tang等[16]通過引入學習型過完備字典實現最優稀疏表征，并采用字典訓練與總變差量的最小化共同抑制偽吉布斯效應，避免了偽影的干擾，但該方法需要已知噪聲的統計信息，在實際應用中需要前期統計與實驗；徐彥凱等[17]提出了基于自適應奇異值分解的局部正交壓制方法，通過將殘余信號從隨機噪聲中的二次分離，提高了資料的信噪比。

以上時頻分析方法受不確定準則限制，只能在時間分辨率與頻率分辨率之間采取折中處理[18]。針對這一缺陷，有學者提出了用于處理時序非穩態信號的經驗模態分解(EMD)方法[19-20]。基于該方法的時頻分析可以在時間、頻率尺度上均達到較高的分辨率。EMD方法在地震資料噪聲壓制和屬性分析中均獲得了較好效果[21-23]。而對于復雜的時序非平穩信號，僅在時間域進行分解往往效果不佳，或出現假頻等現象，而采用EMD方法，將單道地震記錄或剖面在合適的域內分解為一系列本征模態分量(IMF)，可以表征不同頻段內的局部能量異常。傳統EMD方法的適用條件可參見文獻[21-24]。Bekara等[23]率先將EMD方法用于頻率—空間域中的各頻率切片數據，將分解得到的第一模態作為噪聲剔除，重構之后得到噪聲壓制后的信號，取得一定效果。

EMD方法不需要設置針對性的參數，具有較強的泛化能力。但如果在頻率—空間域的頻率切片體上進行EMD處理，邊界效應等現象常常導致分解過程精度低或不穩定[25]。在EMD方法的基礎上發展的變分模態分解(VMD)方法，通過將地震數據進行解析化，在時頻域內將地震信號分解為各個模態分量。對于時頻域中帶限頻率切片的處理，VMD方法[26]明顯優于傳統EMD方法，是一種自適應、完全非遞歸的模態變分方法。基于奇異值分解的VMD方法通過確定模態分解個數，自適應地匹配每種模態的最優中心頻率和有限帶寬，最終得到變分問題的最優解，實現IMF的高效分離。該方法克服了傳統EMD的端點效應與模態分量混疊問題，求解過程具有更高的穩定性。

本文通過希爾伯特-黃變換(HHT)將地震數據解析化，在時頻域進行變分模態分解。通過分頻分解的思路，分析每個模態分量各頻段上的噪聲和有效信號的能量分布，設置合適的分量數，分頻段壓制噪聲，再將有效信號模態分量重構回時頻域。實際資料測試結果表明，相較于常規RNA隨機噪聲壓制方法，該方法有效地壓制了隨機背景噪聲，同時對陡傾角的線性干擾有明顯的壓制效果。

1 算法原理

對于地震剖面中的線性同相軸，可將其描述為空間坐標x和時間t的函數

(1)

式中：w為地震子波；c為波的視速度。在頻率域可由Fourier變換表示為

(2)

式中W(f)為子波的頻譜。

(3)

(4)

(5)

對上式進行傅里葉變換后，可得頻率域解

(6)

(7)

算法基本流程如下：

(1)應用合適時窗選取部分目標數據，進行VMD參數初始化；

(2)應用HHT將目標數據從時間域變換到時頻域；

(3)對各個頻率切片進行VMD分解，確定各模態分量的中心頻率(圖1)；

(4)按能量強弱對VMD分解后得各模態分量(記為VMDIMF)排序，保留前若干個有效VMD-IMF，重構得到濾波后的頻率切片；

(5)將濾波后的所有頻率切片反變換回時空域；

(6)移動時窗，重復上述步驟，直到完成所有數據的處理。

對于三維數據體，僅增加一項正交空間坐標基即可，整體算法流程不變。

對地震數據進行HHT，可分析含噪非平穩信號在不同頻率段內能量沿時間軸的分布。圖2為原始單道地震記錄及其前半段(1～500ms)時頻譜，其中IMF1～IMF6依次對應信號的從高頻至低頻信息(對于單道數據，頻率切片簡化為IMF)。由圖2可見，高頻端和低頻端主要由噪聲構成，有效信息集中在IMF3～IMF5。以IMF5的時頻譜為例，對其進行VMD處理，分解得到各階VMDIMF，如圖3所示。圖中將各VMDIMF的能量進行歸一化，方便對比分析。若分量中存在異常值，可對圖2進行綜合分析，決定是將其作為異常噪聲進行剔除，還是作為局部有效信號予以保留。

用2、25、288Hz諧波合成信號對比VMD法與傳統EMD法分解效果(圖4和圖5)。在分解效率方面，由于VMD遞歸確定最優中心頻率，僅3個VMD-IMF即分解出三個頻率成分，而傳統EMD方法需要6個IMF才分解出288Hz成分；在精度方面，VMD法分解得到的VMDIMF2、VMDIMF3、VMDIMF4均準確定位了合成信號的各頻率成分(圖4)，而EMD法分解結果中IMF2、IMF3與IMF5均為無效分量(圖5)。可見，對諧波合成信號的分解，VMD法具有明顯的效率和精度優勢。

圖1 VMD流程

圖2 單道地震記錄及其各階IMF的時頻譜

圖3 圖2中IMF5的VMDIMF結果

圖4 諧波合成信號的VMD法分解結果(左)及其頻譜(右)

圖5 諧波合成信號的EMD法分解結果(左)及其頻譜(右)

2 模型測試

使用VMD法對地震資料進行處理時，首先需要指定總的模態個數K，該參數直接影響去噪效果。在實際應用前，首先對小塊資料(比如200×200個樣點)進行測試。對模型數據(圖6a)加入頻帶為8～55Hz的高斯隨機噪聲，噪聲平均能量大于有效信號能量的50%，其信噪比為3.09dB(圖6b)。該加噪數據不同K值的VMD法去噪結果如圖7所示。當K=3時背景噪聲壓制效果最好，但有效成分出現較為明顯的割裂痕跡(圖7a紅色箭頭所示)，是由時窗截取的邊緣效應被過低K值放大所致，可根據實際情況，調整時窗長度。在K為5或7時，去噪結果的割裂痕跡明顯得到改善；隨著K值增大，更多的背景噪聲被保留，當K值過大時(該例K大于7時)，出現縱向高傾角處理痕跡(圖7d～圖7f黃框所示)。

圖8為圖7去噪后數據的時頻譜，圖9為去噪結果信噪比(SNR)隨K值的變化曲線，可見：K值較小時，處理結果的信噪比更高；隨著K值的增大，更多噪聲成分被保留，信噪比逐漸下降。

圖6 模型(a)及其加噪(b)數據

圖7 含噪模型數據不同K時的VMD法去噪結果(a)K=3； (b)K=5； (c)K=7； (d)K=9； (e)K=11； (f)K=13

圖8 含噪模型數據不同K時的VMD法去噪后的時頻譜(a)K=3； (b)K=5； (c)K=7； (d)K=9； (e)K=11； (f)K=13

圖9 模態數K與處理結果信噪比相關趨勢

3 實際資料應用

實際地震剖面如圖10a所示，采樣間隔為2ms，記錄長度為8s，道間距為12.5m，具有較強的隨機噪聲干擾，還有部分線性相關噪聲，中深層同相軸連續性較差。本文方法去噪結果如圖10b所示，與原始剖面相比，淺層水平同相軸和深層低伏度同相軸的質量得到明顯提升，連續性得到明顯加強，構造形態更加清晰。與常規RNA去噪方法的處理結果(圖10c)相比，本文方法去噪結果未出現輕微類似斜干擾的人為痕跡(紅色橢圓所示)。通過局部放大對比(圖11、圖12)，處理前淺層水平同相軸背景有殘留的面波成分(圖11a)經處理后得到明顯的壓制(圖11b)；中深層的低幅度構造，在處理前受到線性噪聲和高頻背景噪聲干擾(圖12a)，處理后線性噪聲痕跡幾乎消失，背景隨機噪聲能量得到較強壓制，與有效同相軸的能量不在同一數量級(圖12b)。由處理前、后的頻譜(圖13)可見，由于壓制了低頻面波和高頻背景噪聲能量，低頻能量和高頻能量有所減弱。

圖10 實際地震數據去噪結果(a)原始剖面； (b)本文方法處理剖面； (c)RNA去噪方法處理剖面

圖11 實際地震資料本文方法處理前(a)、后(b)淺層局部放大對比

圖12 實際地震資料本文方法處理前(a)、后(b)深層局部放大對比

圖13 實際地震資料處理前(藍線)、后(紅線)頻譜對比

4 結束語

本文通過HHT得到地震資料在各頻率切片上的能量分布，然后通過VMD對各頻率切片進行分解；分析有效信號與噪聲在時頻切片上的能量分布，在此基礎上優選出有效信號模態分量重構時頻切片；最后反變換回時空域，達到噪聲壓制的目的。分頻處理可以在不同頻率段內對二維信號，甚至單道信號進行有效信息的篩選，從而對局部異常噪聲進行更精細的壓制。模型數據測試和實際資料應用結果表明，與常規RNA方法相比，本文方法可有效壓制線性干擾和隨機背景噪聲，處理后資料信噪比更高。