基于時(shí)頻域的變分模態(tài)分解油氣檢測(cè)技術(shù)及應(yīng)用

張 舜， 周懷來(lái),2， 王元君,2， 周玉毅， 韋 豪

(1.成都理工大學(xué) 地球物理學(xué)院，成都 610059； 2.“油氣地質(zhì)及開(kāi)發(fā)工程”國(guó)家重點(diǎn)實(shí)驗(yàn)室，成都 610059)

0 引言

地震波在地下介質(zhì)中傳播時(shí)，受外界摩擦力和地層巖性變化等原因，會(huì)伴隨著相應(yīng)的能量與頻率的損失和衰減。近年來(lái)，大量研究者依照這種衰減屬性，結(jié)合時(shí)頻分析技術(shù)提取“衰減梯度”進(jìn)行含氣性儲(chǔ)層檢測(cè)的方法，取得了有效的應(yīng)用效果[1-5]。

Mitchell[6]提出了EAA方法，通過(guò)計(jì)算高頻段擬合斜率來(lái)計(jì)算衰減吸收系數(shù)，但該計(jì)算方法僅適用于高信噪比或平穩(wěn)的頻譜，而地震信號(hào)高頻段中往往摻雜著噪音的干擾，限制了該方法的應(yīng)用；張固瀾等[7]利用改進(jìn)的廣義S變換對(duì)低頻段進(jìn)行衰減梯度的提取，并進(jìn)行了驗(yàn)證，但該方法需要避免強(qiáng)干擾信號(hào)的影響，需要事先掌握干擾信號(hào)的頻率特征；熊曉軍[8]使用三參數(shù)小波變換獲得更高精度的時(shí)頻譜結(jié)合曲線非線性擬合提取衰減梯度，提高了儲(chǔ)層預(yù)測(cè)的可靠性；薛雅娟等[9]將EEMD與小波變換集合，通過(guò)提取特征剖面，進(jìn)一步提取更高精度的衰減梯度，為含氣性檢測(cè)提供了一種新思路；于婭[10]提出了匹配追蹤提取高頻衰減梯度的方法，其對(duì)淺層油氣檢測(cè)更為敏感準(zhǔn)確。

Huang[11]提出了經(jīng)驗(yàn)?zāi)Ｊ椒纸?EMD)方法，該方法可以用于對(duì)非線性和非平穩(wěn)信號(hào)的分解，但存在端點(diǎn)效應(yīng)和模態(tài)混疊效應(yīng)；Wu等[12]提出集合經(jīng)驗(yàn)?zāi)B(tài)分解(EEMD)作為EMD的一種改進(jìn)算法，但分解后的殘留噪聲導(dǎo)致重構(gòu)誤差較大，這也限制了其應(yīng)用；Yeh[13]進(jìn)一步改進(jìn)提出互補(bǔ)集合經(jīng)驗(yàn)?zāi)B(tài)分解(CEEMD)，該方法不僅改善了模態(tài)混疊而且保證了重構(gòu)的精確性,但受遞歸分解方式的影響，不能完全抑制模態(tài)混疊現(xiàn)象；Dragomiretskiy[14]提出了變分模態(tài)分解(VMD)，VMD彌補(bǔ)了以往模態(tài)分解沒(méi)有數(shù)學(xué)理論作支撐的弊端，它有著堅(jiān)實(shí)的數(shù)學(xué)基礎(chǔ)做支撐，并通過(guò)非遞歸的方式解決了模態(tài)混疊問(wèn)題，在機(jī)械故障、橋梁和地球物理等領(lǐng)域得到了廣范應(yīng)用[15-16]；Zhou[17]應(yīng)用VMD結(jié)合奇異值法對(duì)地震數(shù)據(jù)進(jìn)行去噪，比起其他方法，該方法去噪效果更高，有效信號(hào)保留相對(duì)完整；龍丹[18]利用VMD改進(jìn)的HHT方法，顯示了其在突出地質(zhì)信息方面的時(shí)頻聚焦的優(yōu)越性。

筆者受前人研究的啟發(fā)，提出了利用VMD與高時(shí)頻分辨率的廣義S變換結(jié)合來(lái)加強(qiáng)含氣儲(chǔ)層的弱反射信號(hào)，從而獲得更精確的衰減梯度進(jìn)行含氣檢測(cè)的方法。通過(guò)對(duì)實(shí)際工區(qū)資料的處理，結(jié)合實(shí)際已知含氣性資料及測(cè)井?dāng)?shù)據(jù)對(duì)本文方法進(jìn)行了驗(yàn)證。結(jié)果表明，本文方法對(duì)含氣目的層有很好地指示效果和應(yīng)用價(jià)值，有望在以后的油氣檢測(cè)應(yīng)用中進(jìn)一步發(fā)揮作用和推廣。

1 基本原理

1.1 變分模態(tài)分解原理

傳統(tǒng)的EMD通過(guò)遞歸的方式對(duì)非線性信號(hào)進(jìn)行分解，主要依靠對(duì)極值點(diǎn)的需找方法，缺少數(shù)學(xué)理論的支撐，存在嚴(yán)重的模態(tài)混疊現(xiàn)象。CEEMD是一種噪聲輔助的改進(jìn)經(jīng)驗(yàn)?zāi)Ｊ椒纸馑惴āＫ鼘⒁欢ǚ鹊恼?fù)成對(duì)白噪聲添加到時(shí)間序列，然后由EMD分解并保存結(jié)果，改善了EMD的模態(tài)混疊現(xiàn)象，但沒(méi)有改變EMD的遞歸分解模式的本質(zhì)。

變分模態(tài)分解(Variational Mode Decomposition，VMD)不同于EMD和CEEMD,可以將其看作是多個(gè)維納濾波器的組合，是一種尋找最優(yōu)解的變分問(wèn)題。VMD具有堅(jiān)實(shí)的數(shù)學(xué)基礎(chǔ)，它改變了以往經(jīng)驗(yàn)?zāi)B(tài)分解類(lèi)方法遞歸分解的特征，是一種完全非遞歸的變分問(wèn)題分解模型，克制了模態(tài)混疊現(xiàn)象。VMD的分解過(guò)程是假設(shè)每一個(gè)分量函數(shù)都圍繞一個(gè)固定的中心頻率ωk，通過(guò)尋找最優(yōu)解，圍繞中心頻率擾動(dòng)，不斷更迭每個(gè)模態(tài)函數(shù)和中心頻率，從而使得最終分解的每個(gè)模態(tài)函數(shù)為具有窄帶特性的調(diào)幅-調(diào)頻(AM-FM)信號(hào)。

調(diào)幅-調(diào)頻信號(hào)表示為：

uk(t)=Ak(t)cos(φk(t))

(1)

式中：Ak(t)、φk(t)分別表示各模態(tài)函數(shù)的瞬時(shí)振幅和瞬時(shí)相位。

變分模態(tài)分解(VMD)主要對(duì)信號(hào)進(jìn)行以下變換計(jì)算處理[14]：

VMD是一個(gè)約束變分問(wèn)題，由式(2)、式(3)表示。

(2)

(3)

式中：δ(t)為狄里克萊函數(shù)；uk(t)為第K個(gè)模態(tài)函數(shù)。

每個(gè)模式的帶寬由解析信號(hào)的L2范數(shù)決定。引入二次罰項(xiàng)α和拉格朗日乘子λ，將式(2)的約束問(wèn)題轉(zhuǎn)化為非約束問(wèn)題：

(4)

式中：α為懲罰因子；f(t)為輸入信號(hào)；λ(t)為L(zhǎng)agrange函數(shù)。

(5)

(6)

(7)

5) 重復(fù)步驟2)～步驟4)，對(duì)于給定判別精度e>0，若

(8)

則停止迭代。最終，得到具有AM-FM信號(hào)特征的K個(gè)IMF分量信號(hào)。

1.2 廣義S變換原理

設(shè)h(t)為平方可積的函數(shù)，則可以定義h(t)的S變換為：

(9)

式中：τ為時(shí)間；f為頻率；S(τ,f)為h(t)的S變換。

其中,S變換的基本小波函數(shù)定義為：

gf(t)exp(-i2πft)

(10)

式中：gf(t)為高斯窗函數(shù)，可定義為：

(11)

由于S變換的基本小波窗函數(shù)形態(tài)固定，使其不能夠依據(jù)實(shí)際情況來(lái)調(diào)整窗口的大小，在高頻部分分辨率低，導(dǎo)致在實(shí)際應(yīng)用中不能得到較好的效果。

廣義S變換是在S變換上的改進(jìn)算法，對(duì)于S變換現(xiàn)在已經(jīng)有多種改進(jìn)方法，本文采用陳學(xué)華提出的二參數(shù)廣義S變換法[19]，通過(guò)加入兩個(gè)參數(shù)rgs和ρ來(lái)調(diào)節(jié)高斯窗函數(shù)(式(11))的大小，見(jiàn)式(12)。

(12)

從而得到信號(hào)h(t)的二參數(shù)廣義S變換為

(13)

通過(guò)調(diào)整rgs和ρ，就可以獲得高的時(shí)頻分辨率。

1.3 衰減梯度求取

“低頻伴影”技術(shù)已經(jīng)成為含氣性檢測(cè)一種常用方法，它主要依據(jù)當(dāng)遇到油氣時(shí)“低頻能量增強(qiáng)，高頻能量減弱”的特性。其需要對(duì)不同頻段的剖面進(jìn)行對(duì)比，且需要注意頻率的選取范圍。衰減梯度已經(jīng)成為油氣檢測(cè)的一個(gè)重要屬性,衰減梯度主要依據(jù)在含氣目的層高頻能量快速衰減的原理，其僅需一個(gè)剖面即可反映我們所需信息，且不需人為尋找頻段。衰減梯度值大的地方往往表示此處存在強(qiáng)振幅異常，高頻衰減較快，排除其他因素影響一般表現(xiàn)為含油氣反應(yīng)。衰減梯度的提取離不開(kāi)時(shí)頻分析方法，筆者采用具有高時(shí)頻精度的廣義S變換來(lái)進(jìn)行衰減梯度的提取。

衰減梯度通常采用EAA方法求取，根據(jù)粘彈性理論，依據(jù)均勻不完全彈性介質(zhì)產(chǎn)生的吸收效應(yīng),使地震波振幅隨地震波傳播距離的增加呈指數(shù)衰減的特性，擬合出頻譜中高頻衰減異常部分，這里我們使用廣義S變換來(lái)獲得相應(yīng)的頻譜信息，即對(duì)于地震信號(hào)h(t)的廣義S變換時(shí)頻分布Xω(τ，f)，它與衰減梯度a的擬合關(guān)系可以表示為:

Xω(τ,f)=ce-af

(14)

式中的c為常數(shù),然后對(duì)式(14)中兩側(cè)分別取對(duì)數(shù)：

ln[Xω(τ,f)]=lnc-af

(15)

然后，通過(guò)計(jì)算總能量85%與65%對(duì)應(yīng)的位置頻率與所對(duì)應(yīng)的地震波能量，通過(guò)擬合斜率，對(duì)式(15)進(jìn)行變換，可得式(16)。

(16)

EAA方法示意圖見(jiàn)圖1。

圖1 EAA方法原理Fig.1 EAA schematic

應(yīng)用VMD結(jié)合廣義S變換提取衰減梯度的方法步驟：

1) 對(duì)地震數(shù)據(jù)剖面進(jìn)行VMD分解。對(duì)分解后的模態(tài)函數(shù)進(jìn)行相關(guān)性計(jì)算，從相關(guān)系數(shù)中選取具有強(qiáng)相關(guān)性的分量作為特征分量，將選取的特征分量合成重構(gòu)地震剖面。

2) 衰減梯度提取。傳統(tǒng)的EAA技術(shù)只適用于高信噪比的平滑頻譜，利用VMD分解后，獲得高信噪比的信號(hào)，再對(duì)重構(gòu)后的地震剖面逐道進(jìn)行二參數(shù)廣義S變換，利用最小二乘法來(lái)擬合求取頻譜斜率，并將其作為該時(shí)間點(diǎn)的衰減梯度值，逐道逐點(diǎn)計(jì)算衰減梯度值，最終得到整個(gè)衰減梯度剖面。

3)通過(guò)衰減梯度剖面，結(jié)合實(shí)際工區(qū)已知的相關(guān)資料，對(duì)含氣儲(chǔ)層進(jìn)行分析。

具體流程如圖2所示。

圖2 VMD與廣義S變換提取衰減梯度流程Fig.2 Schematic diagram of VMD and generalized S transform extraction attenuation gradient

需要值得注意的是VMD分解的參數(shù)對(duì)結(jié)果影響較大,K值過(guò)大會(huì)導(dǎo)致分解過(guò)度，產(chǎn)生虛假分量；K值過(guò)小，則不利于各模態(tài)分量的中心頻率的捕獲。K值過(guò)大或過(guò)小都會(huì)直接影響到結(jié)果的準(zhǔn)確性，在多次試驗(yàn)中發(fā)現(xiàn)，K值一般取值在3～5之間時(shí)較為合適。

2 模型測(cè)試

為驗(yàn)證VMD的優(yōu)越性，我們用主頻分別為30 Hz、45 Hz和50 Hz的單頻信號(hào)y1、y2、y3合成信號(hào)y分別進(jìn)行CEEMD分解和VMD分解進(jìn)行比較。其中分量信號(hào)表達(dá)式為y1=sin(2π*30t),t∈[0,1000](圖3(a));y2=cos(2π*45t),t∈[600,1000](圖3(b));y3=sin(2π*50t+15π),t∈[0,400](圖3(c))；y(圖3(e))是將y1、y2、y3疊加后再加入白噪聲(圖3(d))后的合成信號(hào)。

圖3 合成信號(hào)及分量信號(hào)Fig.3 Synthesized signal and component signal(a)分量信號(hào)y1;(b)分量信號(hào)y2;(c)分量信號(hào)y3;(d)噪聲信號(hào);(e)合成信號(hào)y

為了展示VMD分解方法的優(yōu)勢(shì)，對(duì)圖3中合成信號(hào)y進(jìn)行分解。圖4為CEEMD分解后得到的模態(tài)分量及其對(duì)應(yīng)的振幅譜。由圖4發(fā)現(xiàn)，CEEMD分解中IMF1分量主要為高頻噪聲分量；IMF3為30 Hz分量，出現(xiàn)輕微模態(tài)混疊現(xiàn)象； IMF2與IMF4出現(xiàn)了嚴(yán)重的模態(tài)混疊現(xiàn)象，這主要是因?yàn)镮MF1分解不完全，及殘留噪聲等因素影響，通過(guò)遞歸傳遞向下傳遞分解所導(dǎo)致的。

圖4 CEEMD分解模態(tài)分量及其振幅譜Fig.4 CEEMD decomposed modal components and the amplitude spectrum(a)分解模態(tài)分量;(b)分解振幅譜

圖5為VMD分解的分量及其對(duì)應(yīng)的振幅譜。由圖5中，可以清楚地看到，合成信號(hào)經(jīng)過(guò)VMD分解后，完美的將30 Hz、45 Hz、50 Hz三個(gè)不同頻率成分的分量信號(hào)分解了出來(lái)，在IMF4中主要為添加的噪聲信號(hào)，有效克服了模態(tài)混疊問(wèn)題，體現(xiàn)了VMD的優(yōu)越性。這是由于VMD的非遞歸分解完全改變了以往的遞歸模態(tài)分解，且有著堅(jiān)實(shí)的數(shù)學(xué)基礎(chǔ)做約束，從而其分解分量都為窄帶信號(hào)且約束在其中心頻率附近。與CEEMD分解得到的8個(gè)IMF分量相比，VMD只需分解4個(gè)分量就可以分解出每個(gè)子信號(hào)，分解后的模態(tài)分量更加精確,且VMD可以有效對(duì)噪聲信號(hào)進(jìn)行分離，這也是很多研究者選擇使用VMD去噪的原因之一。

圖5 VMD分解模態(tài)分量及其振幅譜Fig.5 VMD decomposition mode components and the amplitude spectrum(a)分解模態(tài)分量;(b)振幅譜

3 實(shí)際工區(qū)應(yīng)用

本文所用數(shù)據(jù)為南海某實(shí)際工區(qū)的疊后地震數(shù)據(jù)中截取的一個(gè)連井剖面如圖6(a)，截取時(shí)間為2.4 s～2.85 s，受采集方式和外界其他因素等影響，儲(chǔ)層較周邊反射層響應(yīng)較弱。其中包含了井1、井2和井3三口井，測(cè)井曲線為密度曲線，結(jié)合測(cè)井曲線與已知相關(guān)資料，得知井1為高產(chǎn)井，井2為含水氣井(含氣量少)，井3為弱含氣井，筆者主要對(duì)高產(chǎn)井井1含氣性進(jìn)行研究。由實(shí)際工區(qū)資料及測(cè)井解釋結(jié)果得知井1中包含油氣層分別位于如圖所指位置(圖6(b))。

圖6 實(shí)際工區(qū)Fig.6 Orignal seismic section(a)實(shí)際工區(qū)連井剖面；(b)井1放大剖面

首先對(duì)井1的井旁地震道作VMD分解見(jiàn)圖7，其中參數(shù)為K=3。并將各分量與原始地震道進(jìn)行相關(guān)性計(jì)算(表1)，從表1中可以得到IMF2相關(guān)系數(shù)最大，表現(xiàn)為強(qiáng)相關(guān)。對(duì)原始地震剖面進(jìn)行VMD，得到IMF1～I(xiàn)MF3三個(gè)分量剖面(圖8(a))。從圖8中可以看到，相對(duì)于原始剖面，IMF2剖面很好的包含了主要目的層信息。原始剖面目的反射層能量受巖性等其他因素影響顯示較弱，IMF2剖面對(duì)目的層的反射明顯得到了增強(qiáng)，且經(jīng)過(guò)VMD分解后的信號(hào)為窄帶信號(hào)，具有更高信噪比與頻譜更加平穩(wěn)的優(yōu)點(diǎn)，更有利于后期對(duì)衰減梯度的提取。因此選取IMF2剖面作為后期研究的特征剖面。

圖7 井1井旁地震道及其VMD分解分量Fig.7 Seismic track beside well 1(a)原始信號(hào)；(b)IMF1；(c)IMF2;(d)IMF3

圖8 原始剖面及VMD分解剖面Fig.8 Original section and VMD decomposition section(a)原始剖面;(b)IMF1剖面;(c)IMF2剖面;(d)IMF3剖面

表1 井1井旁地震道與其各IMF分量的相關(guān)系數(shù)

為突出本文方法的可行性與有效性，首先對(duì)圖8(a)進(jìn)行分頻處理，并利用“低頻伴影”的方法來(lái)進(jìn)行含氣性的檢測(cè)，所得結(jié)果如圖9所示。從圖9中可以看到，通過(guò)分頻得到了兩個(gè)不同頻段的分頻剖面，在高頻段(圖9(b))中含油氣層比低頻段(圖9(a))的能量有所減弱，表現(xiàn)出“低頻伴影”效應(yīng)。但從圖9中可以觀察到分頻剖面分辨率低，無(wú)法準(zhǔn)確指示出目的層，且需要對(duì)兩圖進(jìn)行比較才可以得出結(jié)論。

圖9 分頻剖面Fig.9 Frequency-division section(a)10 Hz～20 Hz分頻剖面; (b)20 Hz～30 Hz分頻剖面

圖10為應(yīng)用不同方法提取的衰減梯度剖面。由圖10觀察到，傳統(tǒng)基于原始剖面直接使用廣義S變換提取的衰減梯度(圖10(a))，目的層能量較弱，反應(yīng)不明顯。筆者經(jīng)VMD處理后再使用廣義S變換獲得的梯度剖面，衰減梯度強(qiáng)能量明顯，分辨率高。對(duì)于強(qiáng)含氣井井1，經(jīng)這里VMD結(jié)合廣義S變換獲得的衰減梯度剖面(圖10(b))，在目標(biāo)區(qū)域位置都表現(xiàn)了強(qiáng)衰減梯度反映，表明此處存在振幅異常和頻率異常特征，可能為含油氣儲(chǔ)層。對(duì)比上述所選用的分頻方法，本文方法僅需一個(gè)剖面即可，更為直觀，分辨率更高。為進(jìn)一步確定本文方法有效性，將圖10(b)中井1旁的衰減梯度剖面進(jìn)行放大，并加入測(cè)井曲線(所選測(cè)井為密度曲線)，見(jiàn)圖11。依據(jù)密度曲線含油氣性處通常顯示為低值的特性與已知實(shí)際工區(qū)資料結(jié)合，從圖11可以看到，使用本文方法所得衰減梯度剖面很好地吻合了含油氣部位。綜上所述，證實(shí)了本文方法對(duì)研究工區(qū)儲(chǔ)層含油氣分布分析的可行性與有效性。

圖10 不同方法得到的衰減梯度剖面Fig.10 Attenuation gradient profiles obtained by different methods(a)僅使用廣義S變換的方法;(b)本文方法

圖11 井1衰減梯度剖面Fig.11 Well 1 attenuation gradient profil

4 結(jié)論

由于海上數(shù)據(jù)采集中受采集方式等多方面干擾，造成工區(qū)數(shù)據(jù)中含油氣反射層信息較弱，筆者將VMD結(jié)合二參數(shù)廣義S變換，通過(guò)VMD分解維納濾波的特性和自身分解方式的優(yōu)勢(shì)降低了噪聲的干擾，并加強(qiáng)了地震信號(hào)中弱含氣信號(hào)的響應(yīng)特征。然后，通過(guò)二參數(shù)廣義S變換的靈活性和對(duì)時(shí)頻分辨能力的適應(yīng)性，使得在高頻段得以提取更精準(zhǔn)的衰減梯度，從而獲得更高的分辨力。使用本文方法,可以直觀精確刻畫(huà)寬帶地震數(shù)據(jù)中的特定頻帶含氣儲(chǔ)層的強(qiáng)振幅能量異常與高頻的快速衰減。通過(guò)對(duì)實(shí)際工區(qū)的應(yīng)用處理，進(jìn)一步證明了本文方法在含油氣性檢測(cè)中的有效性和優(yōu)越性。