自適應變分模態分解同態反褶積方法

衛澤，潘樹林*，程祎，茍其勇，王暢

(1.西南石油大學地球科學與技術學院，四川成都 610500；2.中國石油吐哈油田公司勘探開發研究院，新疆哈密 839009；3.中國石油西南油氣田公司頁巖氣研究院，四川成都 610056)

0 引言

反褶積作為一種提高分辨率的處理方法，一直是研究的熱點。Ramachandran等[1]給出了褶積的一般原則，為反褶積方法的發展奠定了基礎。Ulrych[2]首次將同態變換應用于地震子波分離。烏爾里克等[3]進一步將同態理論應用于反褶積，取得了較好的效果。同態反褶積不需要假設子波為最小相位，也不需要假設反射系數具有白噪特征，因此受到了廣泛關注。凌云等[4]提出了零相位同態反褶積算法；郭向宇等[5]利用同態反褶積法對混合相位子波進行了相位估算和矯正；王君等[6]將統計法同態反褶積與幾種已有的方法結合起來，使混合相位反褶積向實用化更進了一步。

同態變換由于特有的優勢得到了廣泛應用。Liao等[7]將混合同態反褶積方法應用于實際斷裂模型分析；王丹等[8]利用同態反褶積提取子波，以定位管道泄露；De Macedo等[9]結合同態反褶積和確定性估計的優勢，結合測井數據更準確地預估了子波。同時，同態反褶積在醫學信號處理領域也發揮了重要作用，Hernndez[10]將同態反褶積算法應用于胎兒心率信號的分析；同態反褶積在聲源定位中也能起到關鍵的作用[11]。

雖然同態反褶積具有較好的適用性，但是其抗噪能力較差，影響了應用效果。Jin 等[12]首次從理論和實驗兩個方面分析了噪聲對同態反褶積處理結果的影響，并對方法進行了改進。在實際數據中由于噪聲的存在，很難準確確定子波和反射系數在同態譜中的界限，限制了方法的使用。因此，目前同態反褶積研究的關鍵問題是如何在同態域分離子波與反射系數。

模態分解方法是一種有效分析不同特征信號的方法。Huang等[13]提出了一種適用于非線性平穩信號分析的經驗模態分解(EMD)方法，在信號處理研究中起到了重要的作用。Wu等[14]針對模態混疊現象提出了改進的EMD算法——集合經驗模態分解方法(EEMD)；Yeh等[15]對EEMD算法做了進一步改進，提出了一種互補集合經驗模態分解方法，蔡俊雄等[16]將該方法用于地震數據去噪。Jegadeeshwaran等[17]提出的變分模態分解 (Variational Mode Decomposition，VMD)方法，能夠自適應地實現信號的頻域劃分及各分量的有效分離，表現出更強的抗噪能力。方江雄等[18]提出了基于頻率域內全局自適應的VMD的地震隨機噪聲壓制方法，實現了信號頻帶的高分辨率、自適應剖分；武迪等[19]提出了一種聯合VMD與包絡導數能量算子的高精度時頻分析方法，能夠將地震信號自適應、非遞歸地分解為一系列具有帶限性質的固有模態函數(IMF)，最終獲取的時頻分布能夠追蹤地震信號的能量變化，從而預測縫洞型儲層。胡瑞卿等[20]將VMD方法應用于地震資料去噪，有效壓制了隨機背景噪聲，同時對陡傾角線性干擾有明顯的壓制效果。針對VMD方法參數確定較難的問題，李華等[21]提出了基于信息熵的參數確定方法，馬洪斌[22]采用包絡熵與蛙跳算法相結合的優化方式，驗證了以包絡熵為判斷標準確定最優預設參數的可行性。目前VMD方法主要應用于地震資料去噪處理，未見其他方面的應用。

本文利用VMD方法在信號分離方面的優勢，結合同態變換特點，提出了一種基于VMD的同態反褶積方法，在同態域采用VMD算法自動分離子波與反射系數，以達到提高地震資料分辨率的目的。

1 方法原理

1.1 同態反褶積

根據褶積模型，地震記錄可以表示為

x(t)=b(t)*r(t)

(1)

式中：b(t)為地震子波；r(t)表示反射系數。將式(1)變換到頻率域，為

X(ω)=B(ω)·R(ω)

(2)

再對兩邊取對數，可得

lnX(ω)=lnB(ω)+lnR(ω)

(3)

(4)

對式(4)進行反傅里葉變換，可得

(5)

(6)

(7)

式中L為截止時間，一般取較小值。然后進行同態反變換就可獲得時間域的子波和反射系數，即

(8)

(9)

(10)

以上處理過程只是一種理想狀況。通常情況下，子波和反射系數的同態譜難以通過一個簡單的濾波器進行分離，特別是當數據中存在噪聲時，同態反褶積方法在實際應用中效果不理想。

VMD能夠自適應地實現信號的頻域劃分以及各分量的有效分離，并且抗噪性、對信號的敏感度都強于其他模態分解方法[23]。如果可以利用VMD方法更合理地分離同態域子波與反射系數，則可以有效提高同態反褶積的應用效果。

1.2 基于重構能量差的自適應VMD分解

VMD方法在應用中，分解層數的確定是一個復雜的過程，通常需要進行多次實驗。針對此問題，本文提出一種基于分解前、后能量差的自適應變分模態分解算法。如果分解層數過多，處理后信號與原始信號的能量差值會發生突變，因此通過判斷能量差值是否發生突變以自適應確定分解層數。

信號的能量可表示為

(11)

式中：μ為信號序列；N為序列長度。當分解層數為k時，定義各模態分量能量之和與原始信號能量的相對能量差及其變化率分別為

(12)

(13)

式中：Ej為第j個模態分量的能量；Es為原始信號的能量。ρ′最大值對應了能量差值最不穩定的分層參數，可以將其作為分層參數選擇的標志。

數據在重構過程中存在能量損耗現象，分解的層數越多損耗越嚴重，通過對ρ求導數，找到求導后的最大值，就可以確定最佳分解層數。

為驗證上述方法的可行性，采用三個單頻諧波合成的信號(圖1)進行實驗。單頻信號頻率分別為10、50、100 Hz，合成信號的能量為508.7350。

根據式(12)、式(13)計算能量差和能量差的變化率，可以自動確定最佳分解層數。對信號進行1層分解，相當于信號沒有變化，因此不對1層分解進行分析。將2～5層分解結果進行統計，結果如表1所示。

表1 合成信號不同層數分解結果統計

由表1可知，在分解層數由3變化到4時，信號能量差具有較大值，差值變化率為最大值。表明分解4層時，重構信號與原始信號能量差值發生突變，因此圖1a數據分解三層最合適。

圖2和圖3分別為圖1a信號的3層和4層模態(IMF)分解結果。在按照4層進行分解時，信號出現了能量畸變，分解出的信號不是期望信號(圖3)。因此當分解層數過多時，會造成分解過度、模態混疊，出現信號雜亂現象(如圖3中IMF3)，不合適的分解層數將無法保證分解結果的可靠性。如按照本文方法確定的層數進行分解，能得到最優結果。

圖1 合成信號(左)及其頻譜(右)

圖2 3層VMD結果

圖3 4層VMD結果

1.3 基于自適應VMD的同態反褶積方法

基于VMD方法分離同態域的子波與反射系數的同態反褶積方法實現步驟如圖4所示。

圖4 基于VMD的同態反褶積流程

為了更好地說明本文方法的基本步驟，采用主頻為15 Hz的零相位Ricker子波合成采樣間隔為1 ms、樣點數為500的單道地震記錄進行反褶積處理(圖5)。

首先對圖5數據進行傅里葉變換，得到數據的頻譜，再對頻譜分別取對數，將子波與反射系數在頻率域的相乘關系轉換為相加關系(圖6)。由圖6可以看出，圖6a的子波信息為低頻信息，而圖6b的反射系數為高頻信息。

圖5 合成的單道地震記錄

對圖6c的對數譜做傅里葉反變換轉換到同態域(圖7)，圖7c對應的地震數據的同態譜結果可以看作子波同態譜(圖7a)與反射系數同態譜(圖7b)相加的結果。

圖6 圖5數據的對數譜

應用本文方法分離的同態域反射系數和子波如圖8所示。與圖7a、圖7b對比可以看出，本文方法分離出的子波形態與實際子波形態基本一致，反射系數位置也比較準確，證明了VMD 方法可在同態域有效分離子波與反射系數信息。與常規同態反褶積相比，經過本文方法處理后，原始記錄的分辨率大幅提高(圖9)。

圖7 圖5數據的同態譜

圖8 圖7c數據VMD法分離結果

圖5～圖9使用了一個簡單的合成記錄對基于VMD的同態反褶積方法實現過程進行了說明。在反褶積過程中，子波隨時間和空間變化會對反褶積效果有很大的影響。

圖9 本文方法與常規同態反褶積分離出子波(a)和反射系數(b)對比

使用一維合成記錄驗證本文方法對時變子波的適用性。采用最小相位時變子波(圖10a)與反射系數褶積合成地震記錄(101～200 ms，主頻為40 Hz；201～300 ms，主頻為35 Hz；301～400 ms，主頻為30 Hz；401～500 ms，主頻為25 Hz)，結果如圖10b所示。圖10c為本文方法處理結果與實際反射系數對比，可以看出，分離出反射系數與實際反射系數位置基本吻合，而常規反褶積方法難以做到。

為測試方法的抗噪聲能力，在圖10b合成記錄中加入高斯噪聲，合成了信噪比為5的記錄(圖11a)。本文方法對含噪記錄的處理結果如圖11b所示?？梢钥闯?，噪聲雖然對處理結果有所影響，但影響很小，不但分離出了統計子波，分離的反射系數與實際位置基本吻合，而且提高了地震記錄的分辨率。

圖10 時變子波合成記錄的本文方法處理結果

圖11 含噪時變子波合成記錄本文方法處理結果

對比本文方法與預測反褶積在處理時變子波記錄的結果(圖12)可見：預測反褶積對實際數據分辨率有所提高，但無法準確恢復反射系數的位置，處理效果較差；本文的基于VMD的同態反褶積方法可以較準確地突出反射系數位置，有效壓制噪聲，提高分辨率。

圖12 常規預測反褶積與本文方法處理結果的對比

2 實際資料應用

為了驗證本文的基于VMD的同態反褶積方法的實用性，對實際疊后數據進行處理，結果如圖13所示。由圖可以看出，經本文方法處理后的剖面，同相軸明顯變細、數量增多、同相軸連續性增強，同時目的層位處理前后振幅無明顯變化，能夠顯著提升地震記錄剖面的分辨率。對比圖13的井位處的放大顯示(圖14)可見：本文方法處理結果可以清晰追蹤原始剖面中無法分辨的薄層，與測井數據合成記錄保持一致，證明了本文方法高分辨處理結果的正確性，為后續解釋提供了資料基礎。

對比本文方法實際數據處理前、后時間切片(圖15)可以看出，橫向分辨率大大提高，斷點更加清晰，斷層更加連續，尤其紅色橢圓內，原本分辨不出的斷層也能清晰刻畫，這證明了本文方法是一種適用性強的高分辨率處理方法。

圖16為圖13剖面目的層的頻譜曲線，可以看出，經本文方法處理，在一定程度上拓寬了頻帶，提升了主頻。提取本文方法處理前(圖13a)、后(圖13b)目的層振幅曲線(圖17)可見，本文方法在拓寬地震數據頻帶的基礎上，處理前、后振幅曲線基本一致。

圖13 本文方法實際數據處理結果

圖14 圖13井位處的放大顯示對比

圖15 本文方法處理前(左)、后(右)的時間切片對比

圖16 原始地震剖面(藍色)與本文方法(紅色)處理結果的頻譜對比

圖17 本文方法處理前、后沿目的層的振幅曲線對比

3 結論

本文提出了一種使用VMD算法對同態域地震信號進行高分辨率處理的方法。該方法能夠較好地分離地震子波與反射系數，有效提高地震資料的分辨率。得出以下幾點結論：

(1)針對VMD算法分解層數不易確定的問題，采用了重構后數據與原數據能量差作為標準，實現了分解層數的自動確定；

(2)對時變子波合成記錄可以較好地分離子波與反射系數，表明本文方法對時變子波具有較好的適應能力；

(3)應用本文方法對實際資料處理，結果與測井數據合成記錄一致，提高了地震分辨率，有利于后續薄層、斷層的解釋。