基于CEEMDAN的地震信號高分辨率時頻分析方法

尚平萍,李 鵬,楊安琪,陳學國

(1.中國石油化工股份有限公司石油物探技術研究院,江蘇南京211103;2.中國石油化工股份有限公司勝利油田分公司勘探開發研究院,山東東營 257015)

近年來,隨著油氣勘探的難度不斷增大,時頻分辨率低的傳統譜分解技術,越來越難以滿足高分辨率地震數據處理和解釋的需求,高分辨率地震時頻分析技術的研究與應用成為一種趨勢。目前常用的地震時頻分析技術主要包括短時傅里葉變換、小波變換、S變換和魏格納-威利分布(Wigner-Ville distribution,WVD)變換等[1],但這些技術均存在一定不足:短時傅里葉變換的時窗長度固定,確定時窗長度后無法調節時間和頻率分辨率[2-3];連續小波變換克服了短時傅里葉變換時窗固定的限制,具有多分辨率特性,信號特征刻畫效果好,但受Heisenberg不確定性原理影響,時頻分辨率有限[4];S變換是小波變換和短時傅里葉變換的組合,它的窗函數形狀固定,且仍受Heisenberg測不準原理的限制[5-7];魏格納-威利分布變換雖然能夠在時頻域很好地刻畫單一頻率成分的信號,但復雜信號會產生嚴重的交叉干擾,影響WVD變換結果的解釋,盡管相關的平滑改進方法能較好地抑制交叉項,但是以降低信號的局部時頻精度為代價的[8]。可見,上述方法雖然具有簡單快速的優點,但是都存在一定的局限性。

1998年,HUANG等[9]提出了一種適用于非線性非平穩信號的時頻分析新方法—希爾伯特-黃變換(HHT),在信號處理研究中受到廣泛關注,該方法首先對非線性非平穩信號進行EMD,獲得固有模態函數(intrinsic mode function,IMF)分量后,再進行希爾伯特變換(HT),從而得到分辨率較高的Hilbert時頻譜。然而,當信號中存在間斷信號或噪聲干擾時,EMD會產生難以克服的模態混疊現象,破壞了每個IMF分量所蘊含的物理意義,降低了分解的準確性,使得基于EMD的地震數據處理方法難以推廣。

為克服EMD的模態混疊效應,WU等[10]提出了集合經驗模態分解(ensemble empirical mode decomposition,EEMD)方法,使信號在不同尺度上具有連續性,但是該方法是以增加計算成本為代價來提高集合平均次數以降低重構誤差[11]的,所以它不是一個完備的分解方法。針對EEMD方法所存在的問題,TORRES等[12]提出了CEEMDAN方法,該方法可獲得更好的模態頻譜分離結果,而且在較低的計算成本之下可重構出精確的原始信號。

本文以CEEMDAN代替EMD并將其與希爾伯特變換(HT)相結合,形成了一種適應于非線性非平穩信號的高分辨率時頻分析方法。該方法對CEEMDAN得到的每個IMF分量進行HT譜分析,計算其瞬時振幅和瞬時頻率,然后生成瞬時頻譜,用于地震時頻分析。對合成記錄和實際數據,采用不同的時頻分析方法得到的時頻譜進行對比,以證明基于CEEMDAN的時頻分析方法的準確性和穩定性。該方法兼具較高的時間分辨率和頻率分辨率,可獲得更精確的時頻譜,因此作為一種地震資料精細解釋的有力工具,為后續的含油氣檢測提供了有力支撐。

1 方法原理

希爾伯特-黃變換(HHT)主要包括Huang變換和Hilbert變換兩個部分[13],Huang變換又稱作EMD[14-16]。我們先利用EMD獲得數目有限的IMF分量,再利用瞬時頻率方法和HT變換獲得信號的時頻譜。EMD將信號分解為許多單分量信號的和[17],在該過程中,IMF分量頻率與分離的先后順序有關,分離時間越早IMF分量頻率越高,反之分離時間越遲則IMF分量頻率越低[18],原信號的最高頻率成分存在于最早得到的IMF分量中。由于經過了多通帶濾波處理,每個IMF分量都是穩態的。對各IMF分量進行HT變換,可獲得三瞬譜:瞬時相位譜、瞬時振幅譜和瞬時頻率譜。

1.1 HHT存在的模態混疊問題及改進方法

盡管與傳統時頻分析方法相比HHT在求取地震記錄的瞬時參數時,可獲得更準確的結果,但當信號不連續或者存在噪聲時,EMD會導致模態混疊現象,這基本是不可避免的,此時每個IMF分量所蘊含的物理意義都受到影響,導致求取的瞬時參數不理想,解釋結果存在誤差,制約了其實際應用范圍。

為了克服或者降低模態混疊效應,一些學者提出了不同的解決方案。FLANDRIN等[19]設定了截止頻率以分離間歇信號,但是門檻值的選擇具有很強的主觀性,故難以獲得較好的效果。EEMD是對EMD的一種改進,其核心思想在于分解原始信號之前加入頻譜分布均勻的高斯白噪聲[11]。EEMD的基本思想是先在原始信號中添加高斯白噪聲,再對其整體進行EMD,由于高斯白噪聲的頻譜分布均勻,這使得不同尺度的信號成分自動映射到相應的參考尺度上,對分解出的分量進行平均以此抵消前期加入的白噪聲。EEMD算法中有兩個重要的參數需要預先設定:添加高斯白噪聲的幅值和集合平均次數,HUANG等[9]依靠經驗確定這兩個參數,既不方便也不客觀。作為一種噪聲輔助分析方法,EEMD在一定程度上克服了EMD的模態混疊效應,改善了分解效果,但同時引入了其它問題,如殘余噪聲、分量解非唯一、白噪聲幅值和集合平均次數不能自適應地獲取等。針對上述問題,YEH等[20]提出了一種改進的EEMD方法,即互補集合經驗模態分解(complete ensemble empirical mode decomposition,CEEMD)方法,可消除IMF分量的殘余噪聲。

在CEEMD方法中,原始信號中加入成對的白噪聲,可產生兩個IMF分量集合。加入噪聲后的信號如下:

(1)

式中:S為原信號;N為白噪聲;M1、M2分別為加入正、負成對噪聲后的信號。

1.2 自適應噪聲的總體集合經驗模態分解方法原理

自適應噪聲的總體集合經驗模態分解(CEEMDAN)[12]本質上也屬于EMD,是對其改進后得到的一種變換形式[21]。與EEMD方法類似,該方法也屬于噪聲輔助分析方法。但不同之處在于:CEEMDAN方法克服了模態混疊效應,可精確地重構出原始信號,并獲得了較好的模態分離譜,同時提高了計算效率。

首先定義算子Ej(·)為給定信號通過EMD得到的第j階模態分量;ωi(n)表示高斯白噪聲,其中i=1,2,…,I;εk系數表示在每個階段的信噪比,其中k=0,…,K。設x(t)為目標信號,CEEMDAN方法包括以下6個步驟[22]。

1) 第一階固有模態分量IMF1(t)的求取方法與EMD方法相同,基于不同的噪聲,重復EMD分解過程I次,計算集合平均值,目標信號x(t)的IMF1(t):

(2)

2) 當k=1時,計算一階殘差r1(t):

(3)

3) 對r1(t)添加經EMD分解后的噪聲分量E1[ωi(t)],再進行一次EMD分解,并定義集合平均值為第二階固有模態分量IMF2(t):

(4)

4) 當k=2,…,K時(其中,K是模態總體數量),計算k階殘差rk(t):

(5)

5) 提取rk(t)+εkEk[ωi(t)]的IMF分量,計算它們的集合平均值,得到目標信號的第k+1階固有模態分量IMF(k+1)(t):

(6)

6) 重復上述步驟,當殘差不能再被分解時,得到的結果即為最終殘差R(t):

(7)

目標信號x(t):

(8)

上式表明,CEEMDAN方法可以精確地重構原始信號。

EEMD方法中所添加的高斯白噪聲的幅值和集合平均次數不能自適應地獲取,針對這一問題,CEEMDAN方法首先設置期望的信號分解相對誤差e,考慮到EMD方法得到的第一個IMF分量不一定能夠準確描述信號中的高頻信息,故可以通過抑制低頻信號來獲取高頻信號,該過程是通過添加白噪聲實現的。根據相關準則[12]確定白噪聲的幅值范圍:

(9)

式中:ε=σn/σs為加入白噪聲幅值標準差σn與原信號幅值標準差σs的比值;ξ=σh/σs為信號中高頻成分幅值標準差σh與原始信號幅值標準差σs的比值。因此,公式(9)可表示為:

(10)

由公式(10)可知加入的白噪聲的幅值范圍。設定e值之后,再根據HUANG等[9]提出的統計規律可得到集合平均次數N:

(11)

2 合成信號時頻譜分析

分別采用EMD、EEMD及CEEMDAN 3種方法對合成信號進行測試,以驗證CEEMDAN方法的優勢。圖1a是各組分信號及其合成信號,自上而下分別是10Hz(0～1s)和20Hz(1～2s)的余弦信號、100Hz的Morlet子波(1.8s處)、50Hz的間歇性余弦信號以及合成信號。圖1b是采用EMD方法分解合成信號得到的結果,可以發現EMD方法將合成信號分解成6個IMF分量,分量存在明顯的模態混疊效應。

圖1c是添加10%的高斯白噪聲且集合平均次數為500時利用EEMD方法分解合成信號得到的結果(圖中只列出前6個IMF分量),可以看出,EEMD方法已在很大程度上降低了模態混疊效應,IMF1和IMF2兩個分量中分別恢復出了100Hz的高頻Morlet子波和50Hz的余弦信號;然而IMF2中摻雜了20Hz的余弦信號,并且IMF3也混合了20Hz和10Hz的余弦信號,這說明仍存在模態混疊效應,但是相比EMD方法得到的結果,其模態混疊效應已經有了明顯降低。

CEEMDAN方法采用了與EEMD方法相同的高斯白噪聲,但集合平均次數為100,分解結果如圖1d所示。圖中只列出了前6個IMF分量,可以看出,IMF1完全恢復了100Hz的高頻子波,基于IMF2、IMF3和IMF4可完全提取50Hz、20Hz和10Hz的簡諧波分量。采用CEEMDAN方法得到的結果幾乎完全克服了模態混疊效應的影響。相較于EEMD方法,CEEMDAN方法得到的分解結果更準確。

從圖2中可以看出,采用CEEMDAN方法重構得到的信號誤差小,幾乎可以忽略。可見,CEEMDAN方法可以精確地重構原始信號。

圖1 各組分信號及合成信號和應用EMD、EEMD、CEEMDAN方法分解結果a 各組分信號及合成信號; b EMD方法分解結果; c EEMD方法分解結果(前6個分量); d CEEMDAN方法分解結果(前6個分量)

圖2 采用CEEMDAN方法得到的重構信號

為了驗證CEEMDAN方法的時頻高分辨率特性,基于圖1a中的合成信號進行了試算,采用不同時頻分析方法獲得的時頻譜如圖3所示。圖3a、圖3b和圖3c分別為采用STFT(short-time Fourier trans-form)、CWT(continuous wavelet transform)和廣義S變換方法得到的時頻譜,可以看出,采用STFT方法得到的時頻譜分辨率單一,能量聚焦性低;采用CWT和廣義S變換方法得到的時頻譜在低頻端頻率分辨率高,在高頻端時間分辨率高,但因受Heisenberg測不準原理限制,時間和頻率分辨率不能同時達到最佳。圖3d、圖3e和圖3f分別為采用EMD、EEMD和CEEMDAN方法得到的時頻譜,這3種方法都克服了傳統方法的不足,得到了較為精確的結果,但是采用EMD方法得到的時頻譜出現了嚴重的模態混疊現象;采用EEMD方法得到的時頻譜模態混疊現象雖有改善,但仍存在;采用CEEMDAN方法得到的時頻譜幾乎不存在模態混疊現象,各頻率分量都清晰可辨。可見CEEMDAN方法得到的時頻譜最為精細和準確。

圖3 對合成信號采用不同時頻分析方法得到的時頻譜a SFFT; b CWT; c 廣義S變換; d EMD; e EEMD; f CEEMDAN

3 實際應用

將CEEMDAN方法應用于某工區的實際地震資料處理,以檢驗該方法的實際應用效果,尤其是該方法的時頻高分辨率特性。圖4對比了分別采用STFT、CWT、廣義S變換及CEEMDAN方法得到的時頻譜結果,可以看出相較于其它3種時頻分析方法,采用CEEMDAN方法得到的時頻譜時間分辨率和頻率分辨率明顯提升,時頻譜主頻的變化特征清晰,分辨率高,這有助于描述地質細節,也證明了該方法的有效性。從圖5可以看出,該工區應用CEEMDAN方法后得到的地震剖面的分辨率明顯提高,地質細節刻畫清晰。

圖4 某工區實際單道地震記錄及采用不同時頻分析方法得到的時頻譜a 實際單道地震記錄; b STFT; c CWT; d 廣義S變換; e CEEMDAN

圖5 某工區應用CEEMDAN方法前(a)、后(b)的地震剖面

研究表明,地震波的吸收衰減性質主要取決于巖石骨架、孔隙度和孔隙流體性質,當地層含油氣時,地層對高頻成分吸收增強,導致高頻成分快速衰減,利用這一特性可以預測地層的含油氣性。如圖6的鉆測井數據所示,紅色表示含油層,橙色代表連井的砂體分布,1井2砂組油層厚度8.6m,1井1砂組含油水層厚度15.2m;2井1砂組油層厚度14.5m，有效厚度6.7m。圖7中黑色表示頻譜變化斜率屬性高值,紅色表示頻譜變化斜率屬性較高值,黃色表示頻譜變化斜率屬性較低值,綠色表示頻譜變化斜率屬性低值,對比圖7a和圖7b可以看出,應用CEEMDAN方法后砂組預測結果與測井數據更匹配,高頻信息更豐富。將給定頻率低頻分量的頻譜變化斜率減去應用CEEMDAN方法后的頻譜變化斜率,得到的吸收屬性剖面如圖8所示,利用該屬性可較好地預測儲層的含油氣性。圖8中黑色表示吸收屬性高值,紅色表示該屬性的較高值,黃色表示該屬性的較低值,綠色表示該屬性的低值,可以看出1井1砂組含油氣性較差,2井1砂組含油氣性較好,這與測井數據吻合,應用CEEMDAN方法后得到的剖面高頻信息更豐富,具有較高的分辨率,有利于提升薄油層的預測效果。

圖6 某工區1井、2井儲層對比

圖7 某工區原始地震剖面頻譜變化斜率屬性剖面(a)及應用CEEMDAN方法后頻譜變化斜率屬性剖面(b)

圖8 某工區連井線吸收屬性剖面

4 結束語

EMD作為HHT的核心技術,存在嚴重的模態混疊效應,信號分析的精度不高。基于EMD的改進算法CEEMDAN,不僅克服了模態混疊效應,而且能夠精確地重構原始信號,同時能夠自適應地確定EEMD算法中需預先確定的兩個重要參數:添加高斯白噪聲的幅值和集合平均次數。

本文以CEEMDAN代替EMD并與HT相結合,提出了一種適用于非線性非平穩信號的高分辨率時頻分析方法。信號頻譜分析和實際應用結果表明,與STFT、CWT、廣義S變換方法等傳統的時頻分析方法相比,本文方法能獲得較高的時間和頻率分辨率,得到更精準的信號頻譜。因此基于CEEMDAN的時頻分析方法可作為地震資料精細解釋的有力工具,應用于精細沉積微相的分析以及薄層儲集體的刻畫。