趙品恒,周懷來,王元君,鄔蒙蒙

(1.成都理工大學(xué) 地球物理學(xué)院,成都 610059;2.油氣藏地質(zhì)及開發(fā)工程國家重點(diǎn)實(shí)驗(yàn)室,成都 610059;3.地球勘探與信息技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室(成都理工大學(xué)),成都 610059)

0 引言

隨著時(shí)代發(fā)展,對能源需求增大,對于地震勘探的技術(shù)要求也日趨提高。針對愈發(fā)復(fù)雜的地質(zhì)條件和構(gòu)造模式,對于掌握和精確提取地震資料的三瞬信息就顯得尤為重要[1]。三瞬參數(shù)的提取已廣泛應(yīng)用于石油天然氣儲(chǔ)層預(yù)測及刻畫領(lǐng)域,能有效地反應(yīng)儲(chǔ)層的特征相關(guān)屬性。當(dāng)?shù)卣鸩ù┻^含氣砂巖地層時(shí)高頻能量衰減劇烈,而低頻能量衰減相對較小,因此可通過瞬時(shí)頻率來反映儲(chǔ)層的巖性特征,瞬時(shí)振幅與反射系數(shù)相關(guān)性好,在含油氣區(qū)域地層對應(yīng)的反射強(qiáng)度明顯與非含油氣區(qū)域相異[2]。

近些年,經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)[3]被廣泛研究應(yīng)用在地震信號(hào)時(shí)頻分析處理中,但因模態(tài)混疊和端點(diǎn)效應(yīng)等問題導(dǎo)致該方法不能較好地適用在地震數(shù)據(jù)中。Wu等[4-6]在此算法(EMD)的基礎(chǔ)上,改進(jìn)得到了集合經(jīng)驗(yàn)?zāi)B(tài)分解(EEMD),在一定程度上克服了這些問題。Yeh[7]進(jìn)一步改進(jìn)算法得到互補(bǔ)集合經(jīng)驗(yàn)?zāi)B(tài)分解算法(CEEMD),提高了信號(hào)重構(gòu)精度和計(jì)算效率,但該方法仍有殘余噪聲的問題存在。因此Colominalet等[8]提出改進(jìn)完備經(jīng)驗(yàn)?zāi)B(tài)分解(ICEEMD),該方法所得的模態(tài)分量是由其局部均值來確定,可克服模態(tài)混疊問題[9]。雖然改進(jìn)后的希爾伯特-黃(ICEEMD-Hilbert)在提高地震資料分辨率、消除噪音等方面取得了良好的應(yīng)用效果,但其提取的三瞬參數(shù)仍然會(huì)存在端點(diǎn)效應(yīng)等問題,從而導(dǎo)致解釋結(jié)果誤差較大[10-11]。

TK能量算法[12-14]是一種利用差分運(yùn)算來求取瞬時(shí)屬性的方法,具有操作簡單、有效性高等優(yōu)勢,但缺陷是會(huì)出現(xiàn)負(fù)值。于是O'Toole[15]提出頻率加權(quán)能量算子,即FWEO。相較于TK能量,其不存在負(fù)值,在瞬時(shí)能量追蹤上有更高的精度,具有更強(qiáng)的魯棒性[16-17];FWEO與TK都是只適用于單分量信號(hào)的算子,因此通過ICEEMD方法對地震信號(hào)分解后可以滿足FWEO的應(yīng)用條件。

筆者分別用改進(jìn)的希爾伯特-黃(ICEEMD-Hilbert)和ICEEMD-FWEO算法,在南海某工區(qū)提取瞬時(shí)振幅剖面和瞬時(shí)頻率剖面進(jìn)行比較分析,結(jié)合測井資料對比得出,ICEEMD算法和FWEO能量分離算法組合提取的瞬時(shí)屬性,能夠更加準(zhǔn)確地反映含氣儲(chǔ)層信息。

1 方法原理

1.1 ICEEMD方法

ICEEMD在提取IMF分量時(shí)是加入了特殊噪聲Ek(wi),抑制模態(tài)混疊,增強(qiáng)抗干擾性。其實(shí)現(xiàn)步驟如下:

1)首先對原始信號(hào)加入零均值單位協(xié)方差高斯白噪聲xi=x+ε0ωi,EMD分解得到第1個(gè)殘差:

(1)

2)第1個(gè)模態(tài):

IMF1=x-r1

(2)

3)第2個(gè)殘差和模態(tài):

(3)

IMF2=r1-r2

(4)

4)當(dāng)k= 3、…、K時(shí),IMF分量為:

(5)

IMFk=rk-1-rk

(6)

篩選終止閾值為:

(7)

式(7)中,經(jīng)過多次試驗(yàn),s的值一般在0.05～0.1之間效果最好。

對于ICEEMD方法的邊界條件的做法是:①當(dāng)信號(hào)經(jīng)過樣條插值時(shí),對于樣值點(diǎn)的第一個(gè)和最后一個(gè)點(diǎn),通常采取式(1)將它們同時(shí)作為最大值和最小值;②式(2)根據(jù)最近的極值點(diǎn)將它們作為最大值或者最小值,旨在保證最大值和最小值之間的交替性。

上述表達(dá)式中,x為原始信號(hào),M(·)代表局部平均值運(yùn)算符,Ek(·)代表EMD分解出的第K個(gè)IMF分量。

為了測試ICEEMD方法在分解信號(hào)方面的機(jī)能,模擬一個(gè)復(fù)合信號(hào)來測試。合成的復(fù)合信號(hào)由x1=(2+0.5cos(20πt))*cos(400πt+10sin(10πt))和x2=cos(100πt)組成,如圖1(c)所示。分別對其作EMD和ICEEMD處理,結(jié)果如圖2和圖3所示。從圖2可知,分解結(jié)果中7個(gè)IMF分量和1個(gè)殘差R,且可以明顯看出從IMF2分量就開始出現(xiàn)模態(tài)混疊現(xiàn)象;然后從IMF3～I(xiàn)MF7這5個(gè)分量都受到IMF2的影響。從圖3可知,ICEEMD方法可以完全把合成信號(hào)分離出來,從殘差R上可以看到,分離信號(hào)沒有任何的損失,因此證明ICEEMD方法比常規(guī)的經(jīng)驗(yàn)?zāi)B(tài)分解更好。

圖1 模擬信號(hào)

圖2 EMD分解的結(jié)果

圖3 ICEEMD分解的結(jié)果

1.2 頻率加權(quán)能量分析

FWEO是計(jì)算的信號(hào)倒數(shù)的包絡(luò),它能通過追蹤信號(hào)瞬時(shí)能量來消除噪聲分量,并且它從頻率和振幅兩個(gè)角度估計(jì)信號(hào)的瞬時(shí)信息,因此能更準(zhǔn)確地跟蹤信號(hào)的瞬時(shí)能量變化,與希爾伯特變換相比更簡單。

信號(hào)的瞬時(shí)能量通常以包絡(luò)形式表征,被定義為:

s[x(t)]=|x(t)+jH[x(t)]|2

(8)

其中:x(t)分別表示信號(hào);H[·]表示希爾伯特變換。對于信號(hào)x(t)應(yīng)該是一個(gè)簡單諧波,即具有以下形式:

x(t)=Acos(ω0t+φ)

(9)

為了得到頻率域信息,在式(9)用導(dǎo)函數(shù)引入加權(quán)濾波,則算子定義為:

E[x(t)]=|x'(t)+jH[x'(t)]|2=

x'2(t)+H[x'(t)]2

(10)

其離散形式為:

h2(n+1)+h2(n-1)]-

h(n+1)h(n-1)]

(11)

其中:x(n)和h(n)分別表示離散的信號(hào)和離散的希爾伯特變換。

為了突出FWEO相對于TK的優(yōu)勢,于是通過對比分析兩種方法跟蹤模擬信號(hào)瞬時(shí)能量效果,得到如圖4的結(jié)果。從圖4可以明顯地看出,TK能量算子計(jì)算得到的瞬時(shí)能量具有較多沒有實(shí)際物理意義的能量負(fù)值,而FWEO方法計(jì)算得到的瞬時(shí)能量沒有能量負(fù)值的產(chǎn)生。

圖4 模擬信號(hào)及其瞬時(shí)能量圖

筆者結(jié)合ICEEMD和FWEO能量算子,首先將地震剖面逐道進(jìn)行ICEEMD分解,提取反映烴類信息最多的一個(gè)或幾個(gè)IMF分量剖面進(jìn)行分析;然后利用FWEO算子分別計(jì)算其瞬時(shí)頻率ω(t)和瞬時(shí)振幅a(t);最后為了更好地表征分量特性及烴類信息,對ICEEMD-FWEO得到的瞬時(shí)譜進(jìn)行了高斯平滑。具體的流程如圖5所示。

圖5 ICEEMD-FWEO流程圖

2 模型測試

筆者模擬一個(gè)如下表達(dá)式的復(fù)合信號(hào)來測試該方法:

將y1與y2相加得到如圖6所示,然后對其使用ICEEMD方法分解(圖7),分解為IMF1～I(xiàn)MF6以及一個(gè)殘差R,其中迭代次數(shù)為200。于是分別使用本文方法和ICEEMD-Hilbert方法,對IMF1分量求取各個(gè)的瞬時(shí)振幅和瞬時(shí)頻率,得到的結(jié)果分別如圖8、圖9所示。

圖6 模擬信號(hào)

圖7 ICEEMD方法分解的IMF分量:IMF 1～I(xiàn)MF 6和一個(gè)殘差R

圖8 ICEEMD-Hilbert方法提取的瞬時(shí)屬性

圖9 ICEEMD-FWEO方法提取的瞬時(shí)屬性

圖8(a)、圖9(a)中虛線為IMF1分量波形圖。在圖8(a)中,紅色虛線橢圓中,也就是0 ms～80 ms處,出現(xiàn)明顯的端點(diǎn)效應(yīng)。同樣地,在圖9(b)中,0 ms～100 ms的首段以及990 ms～1 000 ms末端也發(fā)生了端點(diǎn)效應(yīng)。這是由于Hilbert變換窗口效應(yīng)導(dǎo)致的瞬時(shí)振幅出現(xiàn)突變現(xiàn)象。反觀圖9,用本文方法提取的瞬時(shí)振幅和瞬時(shí)頻率都沒有在端點(diǎn)處發(fā)生突變。

為了測試方法的抗噪性,在圖6(c)上加入了SNR=20的高斯白噪聲,對其進(jìn)行ICEEMD分解后得到圖10所示。分別使用本文方法和ICEEMD-Hilbert方法,對IMF4分量求取各個(gè)的瞬時(shí)振幅和瞬時(shí)頻率,得到的結(jié)果見圖11和圖12。

圖10 ICEEMD方法分解的IMF分量:IMF 1～I(xiàn)MF 6和一個(gè)殘差R以及加噪信號(hào)

圖11(a)、圖12(a)中虛線為IMF4分量波形圖。從圖12中可以清晰地看到,筆者使用ICEEMD-FWEO方法提取的瞬時(shí)屬性很好地壓制了噪聲,并且端點(diǎn)處也沒有發(fā)生突變。反觀圖11使用ICEEMD-Hilbert方法提取的瞬時(shí)振幅和瞬時(shí)頻率被噪聲干擾,具有較差的魯棒性。因此,測試信號(hào)處理結(jié)果表明,使用FWEO算子代替Hilbert,在追蹤信號(hào)的瞬時(shí)能量時(shí)具有更高的實(shí)用性。

3 應(yīng)用分析

為了測試本文方法的實(shí)際有效性,選取中國南海某工區(qū)作為研究對象(圖13)。工區(qū)某連井剖面共683道,采樣間隔為2 ms。圖13中包含2口井,曲線為密度曲線,紅色虛線橢圓代表著儲(chǔ)層位置。抽取兩條井旁地震道,分別進(jìn)行ICEEMD分解,結(jié)果如圖14所示。分別求取井旁地震道與各自的IMF分量的相關(guān)系數(shù),如表1所示。

表1 過井地震道與其IMF分量的相關(guān)性

圖13 連井剖面

圖14 兩條過井地震道及其ICEEMD分解結(jié)果

對比圖14以及表1可以發(fā)現(xiàn),IMF1與原始井旁地震道相關(guān)性最大,于是抽取每道的IMF1分量來代替原始地震剖面進(jìn)行屬性分析,其中IMF1分量剖面如圖15所示。在圖15中我們可以清楚地發(fā)現(xiàn),儲(chǔ)層位置的細(xì)節(jié)信息要比原始地震剖面更豐富,而且紅色虛線橢圓處的儲(chǔ)層地震響應(yīng)特征要更強(qiáng)。于是,分別使用ICEEMD-Hilbert方法以及本文方法來提取其各自的瞬時(shí)屬性。

圖15 IMF1分量剖面

瞬時(shí)振幅可以反映信號(hào)能量強(qiáng)度。實(shí)際地震記錄的瞬時(shí)振幅通常用于表征地層連續(xù)性和儲(chǔ)層流體和巖性的橫向變化。強(qiáng)振幅通常與儲(chǔ)層相對應(yīng),所以據(jù)此屬性可觀察到儲(chǔ)層的位置和形狀。從圖16中可以明顯看出,使用FWEO算子比Hilbert有更高的時(shí)間分辨率,尤其是在儲(chǔ)層的位置,即紅色虛線橢圓處,儲(chǔ)層處的地震響應(yīng)更加強(qiáng)烈。同時(shí)在圖16(a)中,黑色箭頭所指的地方有很明顯的端點(diǎn)效應(yīng),反觀本文方法提取的瞬時(shí)振幅剖面上則沒有。因此筆者利用FWEO算子代替Hilbert變換在實(shí)際地震數(shù)據(jù)求取瞬時(shí)振幅有更高的時(shí)間分辨率,可以更好地應(yīng)用于儲(chǔ)層的識(shí)別與刻畫。

圖16 瞬時(shí)振幅剖面圖

圖17為兩種方法提取的瞬時(shí)頻率剖面,瞬時(shí)頻率作為識(shí)別油氣的原理是,瞬時(shí)頻率反映的是地震目的層段瞬時(shí)主頻隨著巖性的變化而改變。當(dāng)瞬時(shí)頻率經(jīng)過含油氣的儲(chǔ)層時(shí),由于高頻能量被吸收,導(dǎo)致地震波經(jīng)過儲(chǔ)層位置時(shí),其瞬時(shí)頻率會(huì)變小,圖17中儲(chǔ)層位置處,紅色代表低瞬時(shí)頻率值。我們可以從圖17中清楚地發(fā)現(xiàn),本文方法提取的瞬時(shí)頻率剖面在2.553 s～2.558 s、2.56 s～2.579 s以及2.615 s～2.62 s左右分別為gas1、gas2和gas3,這與圖中左側(cè)的鈣質(zhì)曲線相吻合,而且相較于ICEEMD-Hilbert方法有更高的分辨率。反觀圖17(a),使用ICEEMD-Hilbert方法只能識(shí)別一套氣層,分辨率不高,無法準(zhǔn)確識(shí)別儲(chǔ)層的精細(xì)構(gòu)造,而且在圖17(a)中黑色箭頭中可以發(fā)現(xiàn)有明顯的端點(diǎn)效應(yīng)。在圖16和圖17中存在掛面現(xiàn)象,橫向變化不自然,其原因可能是地震數(shù)據(jù)在經(jīng)過ICEEMD方法分解時(shí),由于終止閾值為0.2,才會(huì)出現(xiàn)部分IMF分量有端點(diǎn)效應(yīng)造成的。因此,筆者使用ICEEMD方法結(jié)合FWEO算子,在實(shí)際應(yīng)用中提取瞬時(shí)屬性具有一定的實(shí)用性,可以更好地用于儲(chǔ)層預(yù)測。

圖17 瞬時(shí)頻率剖面圖

4 結(jié)論及建議

1)筆者使用ICEEMD方法將復(fù)雜的地震信號(hào),分解為一系列從高頻到低頻的IMF分量,并且選取相關(guān)性大于0.8以上的IMF分量代替原始地震信號(hào)進(jìn)行提取其瞬時(shí)屬性,因此正確選擇ICEEMD分解后的IMF分量,是基于ICEEMD的分析方法進(jìn)行烴類檢測的前提。

2)相較于Hilbert變換,采用FWEO算子在時(shí)間上有更高的分辨率,不僅解決了Hilbert變換由于窗口效應(yīng)導(dǎo)致的端點(diǎn)效應(yīng),而且在瞬時(shí)能量追蹤上克服了負(fù)值的出現(xiàn)。相比于希爾伯特變換,計(jì)算速度更快,更適合數(shù)量龐大的地震數(shù)據(jù)分析。

3)筆者使用ICEEMD結(jié)合FWEO算子,在實(shí)際地震數(shù)據(jù)中提取瞬時(shí)屬性相較于Hilbert變換,可以更好地刻畫儲(chǔ)層位置。