尚平萍,李 鵬,楊安琪,陳學(xué)國

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

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

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

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

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

1 方法原理

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

1.1 HHT存在的模態(tài)混疊問題及改進方法

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

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

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

(1)

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

1.2 自適應(yīng)噪聲的總體集合經(jīng)驗?zāi)B(tài)分解方法原理

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

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

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

(2)

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

(3)

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

(4)

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

(5)

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

(6)

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

(7)

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

(8)

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

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

(9)

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

(10)

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

(11)

2 合成信號時頻譜分析

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

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

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

從圖2中可以看出,采用CEEMDAN方法重構(gòu)得到的信號誤差小,幾乎可以忽略?？梢?CEEMDAN方法可以精確地重構(gòu)原始信號。

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

圖2 采用CEEMDAN方法得到的重構(gòu)信號

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

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

3 實際應(yīng)用

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

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

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

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

圖6 某工區(qū)1井、2井儲層對比

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

圖8 某工區(qū)連井線吸收屬性剖面

4 結(jié)束語

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

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