時頻域變分模態(tài)分解地震資料去噪方法

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

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

0 引言

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

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

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

以上時頻分析方法受不確定準(zhǔn)則限制，只能在時間分辨率與頻率分辨率之間采取折中處理[18]。針對這一缺陷，有學(xué)者提出了用于處理時序非穩(wěn)態(tài)信號的經(jīng)驗?zāi)B(tài)分解(EMD)方法[19-20]?；谠摲椒ǖ臅r頻分析可以在時間、頻率尺度上均達(dá)到較高的分辨率。EMD方法在地震資料噪聲壓制和屬性分析中均獲得了較好效果[21-23]。而對于復(fù)雜的時序非平穩(wěn)信號，僅在時間域進(jìn)行分解往往效果不佳，或出現(xiàn)假頻等現(xiàn)象，而采用EMD方法，將單道地震記錄或剖面在合適的域內(nèi)分解為一系列本征模態(tài)分量(IMF)，可以表征不同頻段內(nèi)的局部能量異常。傳統(tǒng)EMD方法的適用條件可參見文獻(xiàn)[21-24]。Bekara等[23]率先將EMD方法用于頻率—空間域中的各頻率切片數(shù)據(jù)，將分解得到的第一模態(tài)作為噪聲剔除，重構(gòu)之后得到噪聲壓制后的信號，取得一定效果。

EMD方法不需要設(shè)置針對性的參數(shù)，具有較強的泛化能力。但如果在頻率—空間域的頻率切片體上進(jìn)行EMD處理，邊界效應(yīng)等現(xiàn)象常常導(dǎo)致分解過程精度低或不穩(wěn)定[25]。在EMD方法的基礎(chǔ)上發(fā)展的變分模態(tài)分解(VMD)方法，通過將地震數(shù)據(jù)進(jìn)行解析化，在時頻域內(nèi)將地震信號分解為各個模態(tài)分量。對于時頻域中帶限頻率切片的處理，VMD方法[26]明顯優(yōu)于傳統(tǒng)EMD方法，是一種自適應(yīng)、完全非遞歸的模態(tài)變分方法?；谄娈愔捣纸獾腣MD方法通過確定模態(tài)分解個數(shù)，自適應(yīng)地匹配每種模態(tài)的最優(yōu)中心頻率和有限帶寬，最終得到變分問題的最優(yōu)解，實現(xiàn)IMF的高效分離。該方法克服了傳統(tǒng)EMD的端點效應(yīng)與模態(tài)分量混疊問題，求解過程具有更高的穩(wěn)定性。

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

1 算法原理

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

(1)

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

(2)

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

(3)

(4)

(5)

對上式進(jìn)行傅里葉變換后，可得頻率域解

(6)

(7)

算法基本流程如下：

(1)應(yīng)用合適時窗選取部分目標(biāo)數(shù)據(jù)，進(jìn)行VMD參數(shù)初始化；

(2)應(yīng)用HHT將目標(biāo)數(shù)據(jù)從時間域變換到時頻域；

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

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

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

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

對于三維數(shù)據(jù)體，僅增加一項正交空間坐標(biāo)基即可，整體算法流程不變。

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

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

圖1 VMD流程

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

圖3 圖2中IMF5的VMDIMF結(jié)果

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

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

2 模型測試

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

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

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

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

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

圖9 模態(tài)數(shù)K與處理結(jié)果信噪比相關(guān)趨勢

3 實際資料應(yīng)用

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

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

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

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

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

4 結(jié)束語

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