改進的完備經驗模態分解和FWEO能量在南海油氣識別中的應用

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

(1.成都理工大學 地球物理學院,成都 610059;2.油氣藏地質及開發工程國家重點實驗室,成都 610059;3.地球勘探與信息技術教育部重點實驗室(成都理工大學),成都 610059)

0 引言

隨著時代發展,對能源需求增大,對于地震勘探的技術要求也日趨提高。針對愈發復雜的地質條件和構造模式,對于掌握和精確提取地震資料的三瞬信息就顯得尤為重要[1]。三瞬參數的提取已廣泛應用于石油天然氣儲層預測及刻畫領域,能有效地反應儲層的特征相關屬性。當地震波穿過含氣砂巖地層時高頻能量衰減劇烈,而低頻能量衰減相對較小,因此可通過瞬時頻率來反映儲層的巖性特征,瞬時振幅與反射系數相關性好,在含油氣區域地層對應的反射強度明顯與非含油氣區域相異[2]。

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

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

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

1 方法原理

1.1 ICEEMD方法

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

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

(1)

2)第1個模態:

IMF1=x-r1

(2)

3)第2個殘差和模態:

(3)

IMF2=r1-r2

(4)

4)當k= 3、…、K時,IMF分量為:

(5)

IMFk=rk-1-rk

(6)

篩選終止閾值為:

(7)

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

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

上述表達式中,x為原始信號,M(·)代表局部平均值運算符,Ek(·)代表EMD分解出的第K個IMF分量。

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

圖1 模擬信號

圖2 EMD分解的結果

圖3 ICEEMD分解的結果

1.2 頻率加權能量分析

FWEO是計算的信號倒數的包絡,它能通過追蹤信號瞬時能量來消除噪聲分量,并且它從頻率和振幅兩個角度估計信號的瞬時信息,因此能更準確地跟蹤信號的瞬時能量變化,與希爾伯特變換相比更簡單。

信號的瞬時能量通常以包絡形式表征,被定義為:

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

(8)

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

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

(9)

為了得到頻率域信息,在式(9)用導函數引入加權濾波,則算子定義為:

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)分別表示離散的信號和離散的希爾伯特變換。

為了突出FWEO相對于TK的優勢,于是通過對比分析兩種方法跟蹤模擬信號瞬時能量效果,得到如圖4的結果。從圖4可以明顯地看出,TK能量算子計算得到的瞬時能量具有較多沒有實際物理意義的能量負值,而FWEO方法計算得到的瞬時能量沒有能量負值的產生。

圖4 模擬信號及其瞬時能量圖

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

圖5 ICEEMD-FWEO流程圖

2 模型測試

筆者模擬一個如下表達式的復合信號來測試該方法:

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

圖6 模擬信號

圖7 ICEEMD方法分解的IMF分量:IMF 1～IMF 6和一個殘差R

圖8 ICEEMD-Hilbert方法提取的瞬時屬性

圖9 ICEEMD-FWEO方法提取的瞬時屬性

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

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

圖10 ICEEMD方法分解的IMF分量:IMF 1～IMF 6和一個殘差R以及加噪信號

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

3 應用分析

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

表1 過井地震道與其IMF分量的相關性

圖13 連井剖面

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

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

圖15 IMF1分量剖面

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

圖16 瞬時振幅剖面圖

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

圖17 瞬時頻率剖面圖

4 結論及建議

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

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

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