自適應(yīng)變分模態(tài)分解同態(tài)反褶積方法

衛(wèi)澤，潘樹林*，程祎，茍其勇，王暢

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

0 引言

反褶積作為一種提高分辨率的處理方法，一直是研究的熱點(diǎn)。Ramachandran等[1]給出了褶積的一般原則，為反褶積方法的發(fā)展奠定了基礎(chǔ)。Ulrych[2]首次將同態(tài)變換應(yīng)用于地震子波分離。烏爾里克等[3]進(jìn)一步將同態(tài)理論應(yīng)用于反褶積，取得了較好的效果。同態(tài)反褶積不需要假設(shè)子波為最小相位，也不需要假設(shè)反射系數(shù)具有白噪特征，因此受到了廣泛關(guān)注。凌云等[4]提出了零相位同態(tài)反褶積算法；郭向宇等[5]利用同態(tài)反褶積法對(duì)混合相位子波進(jìn)行了相位估算和矯正；王君等[6]將統(tǒng)計(jì)法同態(tài)反褶積與幾種已有的方法結(jié)合起來(lái)，使混合相位反褶積向?qū)嵱没M(jìn)了一步。

同態(tài)變換由于特有的優(yōu)勢(shì)得到了廣泛應(yīng)用。Liao等[7]將混合同態(tài)反褶積方法應(yīng)用于實(shí)際斷裂模型分析；王丹等[8]利用同態(tài)反褶積提取子波，以定位管道泄露；De Macedo等[9]結(jié)合同態(tài)反褶積和確定性估計(jì)的優(yōu)勢(shì)，結(jié)合測(cè)井?dāng)?shù)據(jù)更準(zhǔn)確地預(yù)估了子波。同時(shí)，同態(tài)反褶積在醫(yī)學(xué)信號(hào)處理領(lǐng)域也發(fā)揮了重要作用，Hernndez[10]將同態(tài)反褶積算法應(yīng)用于胎兒心率信號(hào)的分析；同態(tài)反褶積在聲源定位中也能起到關(guān)鍵的作用[11]。

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

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

本文利用VMD方法在信號(hào)分離方面的優(yōu)勢(shì)，結(jié)合同態(tài)變換特點(diǎn)，提出了一種基于VMD的同態(tài)反褶積方法，在同態(tài)域采用VMD算法自動(dòng)分離子波與反射系數(shù)，以達(dá)到提高地震資料分辨率的目的。

1 方法原理

1.1 同態(tài)反褶積

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

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

(1)

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

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

(2)

再對(duì)兩邊取對(duì)數(shù)，可得

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

(3)

(4)

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

(5)

(6)

(7)

式中L為截止時(shí)間，一般取較小值。然后進(jìn)行同態(tài)反變換就可獲得時(shí)間域的子波和反射系數(shù)，即

(8)

(9)

(10)

以上處理過程只是一種理想狀況。通常情況下，子波和反射系數(shù)的同態(tài)譜難以通過一個(gè)簡(jiǎn)單的濾波器進(jìn)行分離，特別是當(dāng)數(shù)據(jù)中存在噪聲時(shí)，同態(tài)反褶積方法在實(shí)際應(yīng)用中效果不理想。

VMD能夠自適應(yīng)地實(shí)現(xiàn)信號(hào)的頻域劃分以及各分量的有效分離，并且抗噪性、對(duì)信號(hào)的敏感度都強(qiáng)于其他模態(tài)分解方法[23]。如果可以利用VMD方法更合理地分離同態(tài)域子波與反射系數(shù)，則可以有效提高同態(tài)反褶積的應(yīng)用效果。

1.2 基于重構(gòu)能量差的自適應(yīng)VMD分解

VMD方法在應(yīng)用中，分解層數(shù)的確定是一個(gè)復(fù)雜的過程，通常需要進(jìn)行多次實(shí)驗(yàn)。針對(duì)此問題，本文提出一種基于分解前、后能量差的自適應(yīng)變分模態(tài)分解算法。如果分解層數(shù)過多，處理后信號(hào)與原始信號(hào)的能量差值會(huì)發(fā)生突變，因此通過判斷能量差值是否發(fā)生突變以自適應(yīng)確定分解層數(shù)。

信號(hào)的能量可表示為

(11)

式中：μ為信號(hào)序列；N為序列長(zhǎng)度。當(dāng)分解層數(shù)為k時(shí)，定義各模態(tài)分量能量之和與原始信號(hào)能量的相對(duì)能量差及其變化率分別為

(12)

(13)

式中：Ej為第j個(gè)模態(tài)分量的能量；Es為原始信號(hào)的能量。ρ′最大值對(duì)應(yīng)了能量差值最不穩(wěn)定的分層參數(shù)，可以將其作為分層參數(shù)選擇的標(biāo)志。

數(shù)據(jù)在重構(gòu)過程中存在能量損耗現(xiàn)象，分解的層數(shù)越多損耗越嚴(yán)重，通過對(duì)ρ求導(dǎo)數(shù)，找到求導(dǎo)后的最大值，就可以確定最佳分解層數(shù)。

為驗(yàn)證上述方法的可行性，采用三個(gè)單頻諧波合成的信號(hào)(圖1)進(jìn)行實(shí)驗(yàn)。單頻信號(hào)頻率分別為10、50、100 Hz，合成信號(hào)的能量為508.7350。

根據(jù)式(12)、式(13)計(jì)算能量差和能量差的變化率，可以自動(dòng)確定最佳分解層數(shù)。對(duì)信號(hào)進(jìn)行1層分解，相當(dāng)于信號(hào)沒有變化，因此不對(duì)1層分解進(jìn)行分析。將2～5層分解結(jié)果進(jìn)行統(tǒng)計(jì)，結(jié)果如表1所示。

表1 合成信號(hào)不同層數(shù)分解結(jié)果統(tǒng)計(jì)

由表1可知，在分解層數(shù)由3變化到4時(shí)，信號(hào)能量差具有較大值，差值變化率為最大值。表明分解4層時(shí)，重構(gòu)信號(hào)與原始信號(hào)能量差值發(fā)生突變，因此圖1a數(shù)據(jù)分解三層最合適。

圖2和圖3分別為圖1a信號(hào)的3層和4層模態(tài)(IMF)分解結(jié)果。在按照4層進(jìn)行分解時(shí)，信號(hào)出現(xiàn)了能量畸變，分解出的信號(hào)不是期望信號(hào)(圖3)。因此當(dāng)分解層數(shù)過多時(shí)，會(huì)造成分解過度、模態(tài)混疊，出現(xiàn)信號(hào)雜亂現(xiàn)象(如圖3中IMF3)，不合適的分解層數(shù)將無(wú)法保證分解結(jié)果的可靠性。如按照本文方法確定的層數(shù)進(jìn)行分解，能得到最優(yōu)結(jié)果。

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

圖2 3層VMD結(jié)果

圖3 4層VMD結(jié)果

1.3 基于自適應(yīng)VMD的同態(tài)反褶積方法

基于VMD方法分離同態(tài)域的子波與反射系數(shù)的同態(tài)反褶積方法實(shí)現(xiàn)步驟如圖4所示。

圖4 基于VMD的同態(tài)反褶積流程

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

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

圖5 合成的單道地震記錄

對(duì)圖6c的對(duì)數(shù)譜做傅里葉反變換轉(zhuǎn)換到同態(tài)域(圖7)，圖7c對(duì)應(yīng)的地震數(shù)據(jù)的同態(tài)譜結(jié)果可以看作子波同態(tài)譜(圖7a)與反射系數(shù)同態(tài)譜(圖7b)相加的結(jié)果。

圖6 圖5數(shù)據(jù)的對(duì)數(shù)譜

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

圖7 圖5數(shù)據(jù)的同態(tài)譜

圖8 圖7c數(shù)據(jù)VMD法分離結(jié)果

圖5～圖9使用了一個(gè)簡(jiǎn)單的合成記錄對(duì)基于VMD的同態(tài)反褶積方法實(shí)現(xiàn)過程進(jìn)行了說(shuō)明。在反褶積過程中，子波隨時(shí)間和空間變化會(huì)對(duì)反褶積效果有很大的影響。

圖9 本文方法與常規(guī)同態(tài)反褶積分離出子波(a)和反射系數(shù)(b)對(duì)比

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

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

圖10 時(shí)變子波合成記錄的本文方法處理結(jié)果

圖11 含噪時(shí)變子波合成記錄本文方法處理結(jié)果

對(duì)比本文方法與預(yù)測(cè)反褶積在處理時(shí)變子波記錄的結(jié)果(圖12)可見：預(yù)測(cè)反褶積對(duì)實(shí)際數(shù)據(jù)分辨率有所提高，但無(wú)法準(zhǔn)確恢復(fù)反射系數(shù)的位置，處理效果較差；本文的基于VMD的同態(tài)反褶積方法可以較準(zhǔn)確地突出反射系數(shù)位置，有效壓制噪聲，提高分辨率。

圖12 常規(guī)預(yù)測(cè)反褶積與本文方法處理結(jié)果的對(duì)比

2 實(shí)際資料應(yīng)用

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

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

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

圖13 本文方法實(shí)際數(shù)據(jù)處理結(jié)果

圖14 圖13井位處的放大顯示對(duì)比

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

圖16 原始地震剖面(藍(lán)色)與本文方法(紅色)處理結(jié)果的頻譜對(duì)比

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

3 結(jié)論

本文提出了一種使用VMD算法對(duì)同態(tài)域地震信號(hào)進(jìn)行高分辨率處理的方法。該方法能夠較好地分離地震子波與反射系數(shù)，有效提高地震資料的分辨率。得出以下幾點(diǎn)結(jié)論：

(1)針對(duì)VMD算法分解層數(shù)不易確定的問題，采用了重構(gòu)后數(shù)據(jù)與原數(shù)據(jù)能量差作為標(biāo)準(zhǔn)，實(shí)現(xiàn)了分解層數(shù)的自動(dòng)確定；

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

(3)應(yīng)用本文方法對(duì)實(shí)際資料處理，結(jié)果與測(cè)井?dāng)?shù)據(jù)合成記錄一致，提高了地震分辨率，有利于后續(xù)薄層、斷層的解釋。