基于CEEMDAN的希爾伯特變換在海底天然氣水合物地震探測中的應(yīng)用

夏秋萍，劉懷山,2

(1.中國海洋大學(xué) 海底科學(xué)與探測技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，山東 青島 266100；2.海洋國家實(shí)驗(yàn)室 海洋礦產(chǎn)資源評價(jià)與探測技術(shù)功能實(shí)驗(yàn)室，山東 青島 266071)

1 引 言

天然氣水合物是天然氣與水在高壓低溫條件下形成的類似冰狀的結(jié)晶物質(zhì)，目前已被正式確定為新的礦產(chǎn)資源，又因?yàn)槠渚G色清潔沒有污染，并且在世界上均有分布，因此多個(gè)國家已經(jīng)開始進(jìn)行天然氣水合物的勘探開發(fā)。研究人員應(yīng)用各種高分辨率地震調(diào)查技術(shù)進(jìn)行勘探研究來尋找水合物，目前已經(jīng)發(fā)現(xiàn)了很多標(biāo)志能夠幫助證明水合物的存在，比如BSR(Bottom Simulating Reflector)、極性反轉(zhuǎn)等[1-5]。地震資料的瞬時(shí)屬性有三個(gè)方面，分別是瞬時(shí)振幅、相位、頻率，以上屬性對于處理復(fù)雜地層，能更加全面地反映儲層的特征，在天然氣儲層預(yù)測中得到了廣泛應(yīng)用。

1998年，Huang等[6]提出了一種新的方法，用來進(jìn)行時(shí)頻分析，這個(gè)方法就是希爾伯特-黃變換(HHT)。該方法通過對信號進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解EMD(Empirical Mode Decomposition)，分解出固有模態(tài)函數(shù)IMF(Intrinsic Mode Function)分量之后，對照實(shí)際信息選取分量進(jìn)行Hilbert變換來獲得瞬時(shí)參數(shù)[7-11]。這個(gè)方法改善了時(shí)頻分析的局限性，地震資料的分辨率也得到了提高。但是信號中若存在間斷或者干擾會引起模態(tài)混疊從而導(dǎo)致分解的準(zhǔn)確程度有所降低，因此在地震數(shù)據(jù)中應(yīng)用EMD進(jìn)行處理與分析的方法也受到了一定制約[12]。為了解決EMD存在的問題，Wu等[13]提出了集合經(jīng)驗(yàn)?zāi)B(tài)分解EEMD(Ensemble Empirical Mode Decomposition)方法，然而該方法無法保證重構(gòu)信號的準(zhǔn)確性，同時(shí)會產(chǎn)生很多非必要的迭代過程，大大降低了計(jì)算效率[14，15]。

2011年，Torres等[16]提出了自適應(yīng)噪聲的完備總體集合經(jīng)驗(yàn)?zāi)B(tài)分解(Complete Ensemble Empirical Mode Decomposition of Adaptive Noise，CEEMDAN)，該方法是建立在經(jīng)驗(yàn)?zāi)B(tài)分解的基礎(chǔ)之上的，對信號進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解，從高頻到低頻獲得一系列IMF分量的過程中，將有限次的自適應(yīng)白噪聲加進(jìn)每個(gè)分解階段來消除分解過程產(chǎn)生的模態(tài)混疊現(xiàn)象，能在較少進(jìn)行平均的條件下,對信號進(jìn)行重構(gòu)且?guī)缀鯖]有誤差，能夠得到更好的模態(tài)頻譜分離結(jié)果。

本文將CEEMDAN應(yīng)用于天然氣水合物的OBS數(shù)據(jù)與垂直纜數(shù)據(jù)中，將經(jīng)過CEEMDAN分解得到的IMF分量進(jìn)行Hilbert變換，獲得瞬時(shí)振幅、相位、頻率，將其所形成的瞬時(shí)屬性圖與只進(jìn)行了Hilbert變換獲得的瞬時(shí)屬性圖進(jìn)行對比，分析得到經(jīng)過CEEMDAN處理后的瞬時(shí)屬性圖分辨率更高，可以當(dāng)作地震勘探數(shù)據(jù)精細(xì)解釋的工具，為進(jìn)行儲層分析提供有效支撐。

2 方法原理

2.1 EMD算法

Huang將具有以下兩個(gè)性質(zhì)的信號定義為IMF，第一是該信號有同等數(shù)量的極值點(diǎn)和零交叉點(diǎn)，如有差異，最多相差1個(gè)；第二是由這個(gè)信號繪制的極值上下包絡(luò)線均值為零。EMD算法如下[6]：

給定實(shí)信號x(t)，首先找出x(t)的所有極值點(diǎn)，利用插值函數(shù)擬合所有的局部極大和極小值點(diǎn)，確定上下包絡(luò)線，求上下包絡(luò)線的平均值為h1(t),

m1(t)=x(t)-h1(t)

(1)

把m1(t)作為新信號，重復(fù)式(1)過程，得到上下包絡(luò)的均值h2(t)

m2(t)=m1(t)-h2(t)

(2)

重復(fù)k次，直到mk(t)滿足IMF條件，得到給定信號的第一階IMF分量

mk(t)=mk-1(t)-hk(t)

(3)

第一階IMF分量是頻率最高的一階分量，記為c1(t)，從給定信號x(t)中減去c1(t),獲得殘余量r1(t),

r1(t)=x(t)-c1(t)

(4)

把r1(t)作為新的給定信號，重復(fù)上述過程以獲得全部的IMF分量，當(dāng)rn(t)比預(yù)定值小或者rn(t)變?yōu)閱握{(diào)函數(shù)的時(shí)候，分解停止。最終得到給定信號:

(5)

2.2 CEEMDAN算法

從本質(zhì)上來說，自適應(yīng)噪聲的總體集合經(jīng)驗(yàn)?zāi)B(tài)分解[17]依然是以EMD為基礎(chǔ)進(jìn)行改進(jìn)得到的方法[18]。類似于EEMD，對信號進(jìn)行噪聲輔助分析的方法，與EMD和EEMD不同的是，CEEMDAN克服了模態(tài)混疊，在提高計(jì)算效率的同時(shí)精確重構(gòu)了原始信號，得到了不錯(cuò)的模態(tài)分離譜[19]。

首先設(shè)Ej(·)為所要分析的信號經(jīng)過EMD獲得的第j階模態(tài)分量；ωi(n)為高斯白噪聲，i=1，2，…，I；εk為各階段的信噪比，k=0，…，K。設(shè)x(t)為目標(biāo)信號， CEEMDAN的算法如下[20]。

IMF1(t)表示為分解出的第一個(gè)固有模態(tài)分量，分解方法等同于EMD，以不同的噪聲為基礎(chǔ)，將EMD分解過程重復(fù)I次，對集合求平均值，得到x(t)的IMF1(t)：

(6)

k=1時(shí)，求一階殘差r1(t)得：

r1(t)=x(t)-IMF1(t)

(7)

將通過EMD方法獲得的噪聲分量E1[ωi(t)]代入到式(7)中后，再次做EMD，將獲得的集合平均值定為第二個(gè)固有模態(tài)分量IMF2(t)：

(8)

k=2，…，K時(shí)(K表示模態(tài)總數(shù))，求k階殘差rk(t)得：

rk(t)=rk-1(t)-IMFk(t)

(9)

將rk(t)+εkEk[ωi(t)]中的IMF分量提出并求取集合平均值，獲得x(t)的第k+1個(gè)固有模態(tài)分量IMF(k+1)(t)：

(10)

將以上步驟反復(fù)至殘差無法再進(jìn)行分解的時(shí)候，獲得結(jié)果就是最終殘差R(t):

(11)

可得到目標(biāo)信號x(t):

(12)

從式(12)可知，CEEMDAN方法能準(zhǔn)確地重構(gòu)原始信號。

設(shè)e為希望得到的信號分解相對誤差，由相關(guān)準(zhǔn)則將白噪聲[17]確定為：

(13)

式(13)中：ε=σn/σs,表示的是白噪聲和原信號的幅值標(biāo)準(zhǔn)差的比；ξ=σh/σs,表示的是高頻信息和原信號幅值標(biāo)準(zhǔn)差的比。所以，式(13)可寫為

(14)

通過Huang 等[14]總結(jié)出的規(guī)律計(jì)算集合平均次數(shù)N，設(shè)置好的e值范圍則為：

(15)

3 實(shí)際數(shù)據(jù)應(yīng)用

本文將基于CEEMDAN的希爾伯特變換用于南海北部神狐海域的OBS(海底地震儀)和垂直纜數(shù)據(jù)中。OBS是一種置于海底接收地震信號的儀器，其接收到的是共接收點(diǎn)記錄，因?yàn)樵诠ぷ髌陂g，OBS一直位于海底接收不同炮間距的地震信號，因此一般情況下，OBS不僅能接收反射記錄，還對深部地層界面的折射波有所記錄。不同于傳統(tǒng)的檢波器在橫向布設(shè)以接收信號，垂直纜是一種垂向布設(shè)檢波器并在底部懸掛海底地震儀的接收裝置。理論上垂直纜在海水中是垂直布設(shè)的，水聽器記錄來自震源激發(fā)后經(jīng)反射界面?zhèn)鞑セ貋淼牡卣鸩ā?/p>

原始數(shù)據(jù)(圖1)經(jīng)過CEEMDAN之后會分解出多個(gè)IMF分量。如圖2和圖3所示截取了OBS數(shù)據(jù)和垂直纜數(shù)據(jù)的前3個(gè)IMF分量，第一階IMF形成的剖面,在地震數(shù)據(jù)中是頻率最高的一階，會包含一些細(xì)節(jié)，也存在高頻噪音，所以文中分析不選擇IMF1。隨著IMF階數(shù)增高,信號幅度會迅速下降,主要成分在前幾階分量上能夠得到較好的反應(yīng)，在IMF3之后的分量形成的剖面都比較模糊，因此本文中選擇以數(shù)據(jù)的第二和第三階分量為基礎(chǔ)來進(jìn)行分析。

圖1 原始數(shù)據(jù)地震記錄Fig.1 Seismic record of original data

圖2 OBS數(shù)據(jù)經(jīng)CEEMDAN處理的前三階分量Fig.2 The first three order components of OBS data processed by CEEMDAN

3.1 瞬時(shí)振幅

瞬時(shí)振幅反映的是振幅包絡(luò)，一種對于反射強(qiáng)度的度量，反映能量的變化，地層反射的強(qiáng)弱表現(xiàn)得更加明顯，在瞬時(shí)振幅剖面上，強(qiáng)反射更強(qiáng)，弱反射更弱。一般情況下，來自水合物穩(wěn)定帶底面的反射也大致與海底平行，稱為BSR，在地震剖面中表現(xiàn)為與海底輪廓相似的強(qiáng)反射層[21],地層中含有天然氣水合物的波阻抗要高于下伏地層中含有游離氣地層的波阻抗，在剖面上BSR的特征是振幅較強(qiáng)，含有水合物地層因?yàn)橛兄容^均勻的密度，所以特征是振幅相對較弱，瞬時(shí)振幅信息比常規(guī)地震剖面更加突出BSR的強(qiáng)反射特征。

圖4和圖5分別為南海神狐海域天然氣水合物的OBS數(shù)據(jù)的和垂直纜數(shù)據(jù)的瞬時(shí)振幅，其中圖4(a)和圖5(a)中的為原始數(shù)據(jù)經(jīng)過希爾伯特變換得到的瞬時(shí)振幅圖，圖4(b)、圖4(c)分別為OBS數(shù)據(jù)經(jīng)過CEEMDAN獲得的IMF2和IMF3通過希爾伯特變換得到的瞬時(shí)振幅圖圖5 (b)、圖5 (c)。從圖4和圖5 中均可看到BSR，用該方法所成的瞬時(shí)振幅圖，不論是OBS還是垂直纜數(shù)據(jù)，分辨率都較圖4(a)、圖5(a)中有所提高，但是可以看到圖4(b)和圖4(c)中BSR海底反射不太連續(xù)，推測是因?yàn)镺BS直接沉放在海底的，在一定程度上會受到接收的直達(dá)波和海底表層反射得到的反射波的干涉效應(yīng)；在圖5中，垂直纜接收信號的直達(dá)波與海底一次反射波能清晰分開，水合物儲層BSR效應(yīng)也很明顯，圖5(b)和圖5(c)顯示的BSR清晰，與圖5(a)中直接經(jīng)過希爾伯特變換得到的比較，分辨率更高。

圖4 OBS數(shù)據(jù)瞬時(shí)振幅Fig.4 Instantaneous amplitude of OBS data

圖5 垂直纜數(shù)據(jù)瞬時(shí)振幅Fig.5 Instantaneous amplitude of vertical cable data

圖6 OBS數(shù)據(jù)瞬時(shí)相位Fig.6 Instantaneous phase of OBS data

圖7 垂直纜數(shù)據(jù)瞬時(shí)相位Fig.7 Instantaneous phase of vertical cable data

3.2 瞬時(shí)相位

瞬時(shí)相位是同相軸連續(xù)性的度量，一種不考慮振幅強(qiáng)度的相位剖面，地震波穿過有巖性差異的地層時(shí)，瞬時(shí)剖面上能反映出地震波的相位變化。圖6和圖7分別為南海神狐海域天然氣水合物的OBS和垂直纜數(shù)據(jù)的瞬時(shí)相位圖，其中圖6(a)、圖7(a)為原始數(shù)據(jù)經(jīng)過希爾伯特變換得到的剖面，圖6(b)、圖6(c)和圖7 (b)、圖7(c)分別為原始數(shù)據(jù)經(jīng)過CEEMDAN得到的IMF2和IMF3分量經(jīng)過希爾伯特變換得到的剖面。從圖6和圖7中可以看到，原始數(shù)據(jù)經(jīng)過希爾伯特變換得到的圖6(a)和圖7(a)中，相位屬性信息復(fù)雜，難以分辨有效信息，而經(jīng)過CEEMDAN處理后的進(jìn)行希爾伯特變換得到的圖6(b)和圖6(c)中，也沒有明顯的追蹤到BSR，在瞬時(shí)相位屬性圖中，如果BSR不連續(xù)，是能在圖中找到，對后續(xù)解釋有很大幫助；而如果BSR界面平行于地層或相交的角度比較小時(shí)，在圖中效果不明顯。因此，判斷在本文的屬性圖中沒有有效識別到BSR的原因，應(yīng)是本文使用的數(shù)據(jù)得到的BSR與地層大致平行或交角較小的緣故。

3.3 瞬時(shí)頻率

瞬時(shí)頻率表示瞬時(shí)相位的時(shí)間變化率。地層是否含有氣可以通過反射波某一時(shí)刻的即時(shí)頻率即瞬時(shí)頻率來反映，在AVO屬性中，游離氣的賦存狀態(tài)主要通過振幅隨偏移距的變化情況來描述，而瞬時(shí)頻率則通過地震波傳播到不同地層的不同頻率響應(yīng)來描述，含游離氣的地層在瞬時(shí)頻率剖面中表現(xiàn)為低頻特征(圖8、圖9)。

圖8 OBS數(shù)據(jù)瞬時(shí)頻率Fig.8 Transient frequency of OBS data

圖9 垂直纜數(shù)據(jù)瞬時(shí)頻率 Fig.9 Instantaneous frequency of vertical cable data

如圖8、圖9所示，分別為南海神狐海域OBS數(shù)據(jù)和垂直纜數(shù)據(jù)的瞬時(shí)頻率圖，由于希爾伯特變換是全通濾波，極其容易受到高頻隨機(jī)噪聲的干擾，如果接收到的地震信號含有高頻噪聲，噪聲會把經(jīng)過希爾伯特變換獲得的瞬時(shí)特點(diǎn)掩蓋。在圖8(a)和圖9(a)中，不論是OBS數(shù)據(jù)還是垂直纜數(shù)據(jù)以希爾伯特變換為基礎(chǔ)得到的瞬時(shí)頻率圖，無法清楚辨別出高低頻信息，難以分辨各層的位置，且無法在圖8(a)、圖9(a)中識別到BSR，而圖8(b)、圖8(c)和圖9(b)、圖9(c)是經(jīng)過CEEMDAN的IMF2和IMF3再做希爾伯特變換得到的瞬時(shí)頻率圖。相較于圖8(a)和圖9(a)，雖然依然含有噪音，但是都能識別到BSR。如圖8和9所指示的反射能量較強(qiáng)的地方，也就是水合物和游離氣共存產(chǎn)生 BSR 相應(yīng)的位置。 BSR以上表示水合物儲集區(qū)(該區(qū)水合物儲層厚度在 30 m 左右)。同時(shí)垂直纜與OBS相比，由于垂直纜接收數(shù)據(jù)具有較高的信噪比和相對較寬的頻帶，對比圖8與圖9可知，圖9中垂直纜數(shù)據(jù)的剖面分辨率更高。

4 結(jié) 論

通過對天然氣水合物實(shí)際數(shù)據(jù)應(yīng)用的分析可知，瞬時(shí)振幅、頻率、相位可以很好地識別出天然氣水合物的賦存特征，基于CEEMDAN的希爾伯特變換可以用在海洋地震勘探數(shù)據(jù)OBS和垂直纜數(shù)據(jù)的分析中，能獲得更加清楚的層間信息，為減少地質(zhì)解釋的多解性提供了可參考的依據(jù)。

1)瞬時(shí)振幅中，垂直纜數(shù)據(jù)經(jīng)過處理后分辨率提高，而OBS數(shù)據(jù)分辨率雖有提高，但BSR有些間斷，是因?yàn)镺BS中包含了縱波分量和橫波分量，橫波地震記錄中水合物的BSR效應(yīng)不明顯，而垂直纜豎直投放在海水中，環(huán)境噪音相對較小，因此垂直纜的分辨率要高于OBS；

2)瞬時(shí)相位中，無論是只經(jīng)過希爾伯特變換處理得到的瞬時(shí)屬性圖還是經(jīng)過CEEMDAN后再進(jìn)行希爾伯特變換的瞬時(shí)屬性圖，都沒有有效追蹤到BSR，判斷BSR與地層大致平行或交角較?。?/p>

3)瞬時(shí)頻率中，經(jīng)過處理后的剖面比只經(jīng)過希爾伯特變換獲得的剖面分辨率更高，能識別到BSR的存在；

4)在獲得的地震瞬時(shí)屬性中，對比垂直纜數(shù)據(jù)和OBS數(shù)據(jù)，垂直纜數(shù)據(jù)能獲得分辨率更高的信息。