變分模態(tài)分解與包絡(luò)導數(shù)算子結(jié)合的時頻分析方法及溶洞儲層預測

武 迪 宋維琪 劉 軍 陳禮豪 陳俊安 楊子鵬

(①中國石油大學(華東)地球科學與技術(shù)學院，山東青島 266580；②中國石化西北油田分公司勘探開發(fā)研究院，新疆烏魯木齊 830011)

0 引言

深層縫洞型碳酸鹽巖儲層中縫洞空間分布散亂，不同縫洞的規(guī)模、形態(tài)和內(nèi)部結(jié)構(gòu)等特征差異明顯，并且由于儲層埋深較大，造成地震信號信噪比低。因此由縫洞引起的“串珠狀”地震響應(yīng)與背景分離度差，儲層識別困難[1-2]。碳酸鹽巖地層中的裂縫、溶洞往往是良好的油氣聚集空間[3-4]，因此縫洞體儲層預測具有十分重要的意義。

油氣的存在導致地震信號頻率和能量異常，通過譜分解技術(shù)，容易發(fā)現(xiàn)隱藏的特定頻段的異常信息，目前已經(jīng)被廣泛用于儲層預測、烴類檢測等方面[5-8]。近年來，在經(jīng)典譜分析方法的基礎(chǔ)上，涌現(xiàn)出經(jīng)驗小波變換[9]、稀疏時頻分析[10]、同步壓縮變換[11]等多種有效算法。作為一種時頻分析手段，聯(lián)合經(jīng)驗模態(tài)分解(Empirical Mode Decomposition，EMD)與希爾伯特變換的希爾伯特—黃變換(HHT)方法由Huang等[12]提出，能夠自適應(yīng)地將信號分解為一系列的固有模態(tài)函數(shù)(Intrinsic Mode Function，IMF)，繼而分析各IMF的瞬時頻率、瞬時振幅等[13]。EMD具有更高的時頻分辨能力，但也存在模態(tài)混疊、端點效應(yīng)等問題[14]。為了解決模態(tài)混疊問題，人們相繼提出引入噪聲輔助的集合經(jīng)驗模態(tài)分解(Ensemble Empirical Mode Decomposition，EEMD)、完備集合經(jīng)驗模態(tài)分解(Complete Ensemble Empirical Mode Decomposition，CEEMD)等方法[15-16]，在一定程度上消除了模態(tài)混疊，已被用于地震資料處理[17]及烴類檢測[18-20]，但由于其本質(zhì)上缺乏完備的理論支撐，在實際應(yīng)用中存在一定局限性。Dragomiretskiy等[21]提出了變分模態(tài)分解(Variational Mode Decomposition，VMD)方法，具有完善的理論基礎(chǔ)[22]，能夠?qū)⒍喾至啃盘栕赃m應(yīng)、非遞歸地分解為一系列具有帶限性質(zhì)的IMF，有效控制了模態(tài)混疊問題。由于VMD的本質(zhì)是維納濾波的推廣，因而具有良好的抗噪性，在地震資料去噪[23-24]及時頻分析領(lǐng)域[25]取得了良好的效果。

Teager等[26]提出了Teager能量算子，Kaiser[27]給出了其離散形式——Teager-Kaiser(TK)算子，因具有良好的局部分析能力及考慮了頻率特性，de Matos等[28]最早將其引入地震勘探鄰域，用于分析碳酸鹽巖儲層。陳學華等[29]、唐湘蓉等[30]在時頻域進行儲層預測及流體識別。Xue等[31]在EMD的基礎(chǔ)上，利用TK算子代替希爾伯特變換獲得瞬時屬性進行聯(lián)合時頻分析，解決了希爾伯特變換的端點效應(yīng)問題，然而TK算子對噪聲敏感的特性限制了其應(yīng)用。O′Toole等[32]提出了具有良好非負特性的包絡(luò)導數(shù)算子(Envelope Derivative Operator，EDO)，與TK算子相比，同為具有良好的追蹤信號瞬時變化能力的非線性能量算子，但EDO抗噪性和穩(wěn)定性更好[33]，已經(jīng)被用于軸承故障檢測等方面[34]。在本文中，為了滿足該類非線性能量算子的應(yīng)用條件，利用VMD技術(shù)將地震信號分解為包含低頻、高頻有效信息的IMF，再引入EDO獲取瞬時振幅、瞬時頻率，克服了希爾伯特變換的端點效應(yīng)且具有更高的時頻分辨率，最終形成一種高精度時頻分析方法，并應(yīng)用于中國西部N區(qū)縫洞型碳酸鹽巖儲層預測，取得了良好的效果。

1 方法原理

1.1 VMD

VMD的目標是將輸入信號分解為一系列子信號，分解后得到的模態(tài)分量具有稀疏特性并且可以重建原始信號。與EMD相比，VMD具有完備的理論基礎(chǔ)，從本質(zhì)上看，VMD是維納濾波的推廣。

VMD非遞歸地將多分量信號自適應(yīng)地分解為具有特定稀疏屬性的帶限模態(tài)分量，首先定義具有中心頻率的有限帶寬的IMF

u(t)=A(t)cosφ(t)

(1)

式中：信號的相位φ(t)需滿足其對時間t的一階偏導數(shù)φ′(t)≥0；A(t)為包絡(luò)。

定義了式(1)的IMF之后，VMD過程包括變分問題的構(gòu)造和求解，若將分析信號分解成K個IMF，則對應(yīng)的約束變分模型為

(2)

為了求解式(2)的最優(yōu)解，引入二次懲罰因子α和Lagrange乘法算子λ(t)，則約束變分問題轉(zhuǎn)換為非約束變分問題。增廣Lagrange公式為

(3)

(4)

式中n為迭代次數(shù)。對應(yīng)的中心角頻率為

(5)

通過

(6)

(7)

則迭代結(jié)束，得到K個IMF。

由VMD分解信號時，需要預先設(shè)定IMF的個數(shù)K，在測試模擬信號時發(fā)現(xiàn)K取l+1(l為合成信號分量的個數(shù))時效果最好。由于實際信號復雜多變，吳文軒等[35]、劉尚坤[36]分別利用峭度及互信息準則等方法確定K值，文中利用經(jīng)驗方法，確定當K=3時可取得最佳的分解效果。

1.2 包絡(luò)導數(shù)算子

1.2.1 TK算子

對于連續(xù)信號x(t)，TK算子被定義為二階微分方程的形式[26]

(8)

x(m)=a(m)cosφ(m)

(9)

式中a(m)為調(diào)幅信號，m為采樣點號。x(m)的離散形式的TK算子為[27]

ψd[x(m)]=x2(m)-x(m-1)x(m+1)

(10)

TK算子是僅利用差分運算計算瞬時能量的非線性算子，考慮了信號的頻率特性，最明顯的優(yōu)勢是良好的局域特性及計算的簡潔、高效性。

1.2.2 EDO

對于信號瞬時能量的求取，典型的方法是以振幅值的平方進行量化，或者以包絡(luò)的形式

(11)

表征。式中H[·]為希爾伯特變換。為了引入頻率域的信息，可以在式(11)中以導函數(shù)的形式引入加權(quán)濾波器。根據(jù)傅里葉變換(FT)的性質(zhì)，有

(12)

則包絡(luò)導數(shù)算子為[32]

(13)

式(13)結(jié)合了信號的頻率變化信息與隨時間變化的包絡(luò)，稱為包絡(luò)導數(shù)算子(EDO)。EDO的定義形式與TK算子較相似，只有第二項存在區(qū)別。O′Toole等[32]詳細對比了TK算子、EDO用于x(t)=Acos(ω0t+φ)(A為振幅，ω0為初始角頻率，φ為相位)或x(t)=Aertcos(ω0t+φ)(r為調(diào)節(jié)系數(shù)，用于描述振幅的變化)形式的簡單信號時的特性，均取得較好效果。對于線性組合信號y(t)=x1(t)+x2(t)

其中

(14)

將TK算子、EDO應(yīng)用于y(t)，分別得到

ψ[y(t)]=ψ[x1(t)]+ψ[x2(t)]+

cos[(ω1-ω2)t+φ1-φ2]}

(15)

Γ[y(t)]=Γ[x1(t)]+Γ[x2(t)]+

a[sin(ω1t+φ1)sin(ω2t+φ2)+

cos(ω1t+φ1)cos(ω2t+φ2)]

=Γ[x1(t)]+Γ[x2(t)]+

acos[(ω1-ω2)t+φ1-φ2]

(16)

式中a=2A1A2ω1ω2。由式(15)、式(16)可見：Γ[y(t)]包含一個等于信號y(t)兩個分量頻率之差的項；ψ[y(t)]則包含了額外的附加項，這將在TK能量計算時出現(xiàn)負值。圖1為信號y(t)及TK算子、EDO能量計算結(jié)果。由圖可見：對于能量的變化特征，TK算子、EDO均取得了良好的追蹤結(jié)果，但在部分極值點處TK算子出現(xiàn)負值，這將在進一步分離能量時出現(xiàn)奇異值；EDO則具有良好的非負特性。

圖1 信號y(t)(上)及TK算子、EDO能量(下)計算結(jié)果

TK算子采用向前差分法獲得離散形式，為了得到EDO的離散形式，定義以下中心差分形式

(17)

因此，離散形式的EDO為

h2(m+1)+h2(m-1)]+

h(m+1)+h(m-1)]

(18)

式中h(·)=Η[x(·)]。

1.2.3 基于VMD-EDO的時頻分析方法

基于VMD-EDO的時頻分析方法如圖2所示。

圖2 基于VMD-EDO的譜分解算法

首先，將地震信號進行VMD，原始數(shù)據(jù)中的有效信息將主要由一個或幾個IMF所反映，選取包含儲層信息最多的IMF進行分析。

然后，采用能量分離(ESA)算法獲得各分量信號EDO能量的瞬時頻率和瞬時振幅[37]。由于ESA算法只適用于單分量信號，不能直接用于地震信號，而VMD方法將地震信號分解為一系列窄帶信號，近似滿足非線性能量算子的計算條件，可獲得EDO能量。通過三點對稱差分運算，得

(19)

(20)

這里，瞬時頻率ω(m)和瞬時振幅A(m)都是時間的函數(shù)，因此可以定義一個三維空間[t，A(t)，ω(t)][38]

(21)

最終可以獲得IMF的時頻分布，進一步可以分離得到不同頻段的高精度瞬時譜。

2 模型驗證

為了驗證VMD-EDO時頻分析方法的效果，建立了地質(zhì)模型(圖3a、表1)。采用主頻30Hz的雷克

圖3 基于模型的VMD-EDO譜分解

子波作為震源進行模擬，并添加7%的隨機噪聲。基于模型的VMD-EDO譜分解結(jié)果表明：在合成地震記錄上縫洞體呈“串珠狀”響應(yīng)(圖3b)；在VMD-EDO低頻剖面上縫洞體呈強能量體特征(圖3c)；相對于地層信息，在VMD-EDO高頻剖面上縫洞體能量衰減明顯(圖3d)。因此VMD-EDO譜分解結(jié)果總體體現(xiàn)了“低頻能量加強、高頻能量衰減”的特征，證明了方法的有效性。

表1 模型參數(shù)

3 實際資料應(yīng)用分析

為進一步驗證VMD-EDO時頻分析方法的效果，選取中國西部N區(qū)的疊后地震資料進行測試。N區(qū)位于順托果勒低隆的北部，處于阿瓦提、滿加爾坳陷與沙雅隆起的結(jié)合部，緊鄰南部海相烴源巖灶，同時發(fā)育本地寒武系烴源巖，構(gòu)造位置有利，油源充足，是油氣長期運移、聚集的有利區(qū)。目前研究表明，N區(qū)奧陶系油氣成藏條件優(yōu)越，儲層地震響應(yīng)主要以同相軸錯斷、異常彎曲和“串珠狀”為主。由于縫洞儲層自身的規(guī)模、形態(tài)和內(nèi)部結(jié)構(gòu)差異明顯以及埋深較大，對于時頻分析精度提出了更高的要求。

圖4為A井井旁地震道CDP105及VMD結(jié)果，圖5為IMF2振幅譜及由TK算子、EDO得到的瞬時振幅。由圖可見：①IMF1為超低頻信息(圖4b上)，而中低頻及高頻有效信息主要分別包含在IMF2(圖4b中)與IMF3(圖4b下)中。②整體來看，對于4600～4700ms層段的油氣異常特征，分別采用TK算子(圖5b上)與EDO(圖5b下)得到的瞬時振幅相近，并且均避免了希爾伯特變換的端點問題。③由于TK能量計算時出現(xiàn)負值，在部分極值點附近瞬時振幅出現(xiàn)部分奇異值，需要采用取絕對值才能避免此問題(圖5b上)；EDO計算結(jié)果更平滑，且EDO良好的非負特性有效避免了奇異值(圖5b下)。

圖6為A井井旁道 CDP109時頻譜。由圖可見，有效信息主要分布在15～40Hz頻段，在4650ms附近出現(xiàn)強能量異常，VMD-EDO時頻譜(圖6a)的時頻分辨率高于連續(xù)小波變換(CWT)時頻譜(圖6b)。

撒利明等[39]采用“低頻共振、高頻衰減”特征檢測油氣，薛雅娟等[18]分別討論了低頻(14～18Hz)、高頻(26～30Hz)的時頻特征，并預測了鄂爾多斯盆地海相碳酸鹽巖儲層(頻帶范圍為12～30Hz)。圖7為過A井地震剖面及VMD結(jié)果。由圖可見：在過A井地震剖面的CDP105、4650ms附近鉆遇油氣，呈小“串珠狀”地震響應(yīng)(圖7a)；IMF2(圖7c)、IMF3(圖7d)分別反映了中低頻、高頻有效信息。圖8、圖9分別為VMD-EDO、CWT-EDO的低頻、高頻剖面。由圖可見：①在含油氣區(qū)域(紅圈位置)，兩種方法都取得了良好的檢測結(jié)果，具體表現(xiàn)為：在低頻剖面上呈不連續(xù)的強亮點及周邊存在弱能量團(圖8a、圖9a)；隨著頻率升高，在高頻剖面上含油氣區(qū)域(紅圈位置)的能量衰減明顯(圖8b、圖9b)。②VMD-EDO方法檢測結(jié)果的分辨率更高，能更好地區(qū)分不同深度的能量異常。圖10為對應(yīng)圖7c、圖7d的VMD-EDO低頻、高頻目的層沿層切片。由圖可見，A井、B井均為良好的產(chǎn)油井，兩口井所在區(qū)域均明顯表現(xiàn)為低頻段強能量異常(圖10a)、高頻能量衰減(圖10b)特征，進一步驗證了方法的有效性。因此，在縫洞型碳酸鹽巖儲層中，基于VMD-EDO的高精度時頻分析算法的儲層預測效果較好。

圖4 A井井旁地震道CDP105(a)及VMD結(jié)果(b)A井為良好的產(chǎn)油井，4600～4700ms(紅線圍成的區(qū)域)為油層

圖5 IMF2振幅譜(a)及由TK算子(上)、EDO(下)得到的瞬時振幅(b)

圖6 A井井旁道 CDP109時頻譜(a)VMD-EDO； (b)CWT

圖7 過A井地震剖面及VMD結(jié)果(a)過A井地震剖面； (b)IMF1； (c)IMF2； (d)IMF3

圖8 對應(yīng)圖7c、圖7d的VMD-EDO低頻(19～22Hz)(a)、高頻(33～35Hz)(b)剖面

圖9 對應(yīng)圖7c、圖7d的CWT-EDO低頻(19～22Hz)(a)、高頻(33～35Hz)(b)剖面

圖10 對應(yīng)圖7c、圖7d的VMD-EDO低頻(19～22Hz)(a)、高頻(33～35Hz)(b)目的層沿層切片

4 結(jié)論與認識

本文在HHT的理論基礎(chǔ)上，提出了一種結(jié)合VMD與EDO能量算子的高精度時頻分析方法，并將其應(yīng)用于中國西部N區(qū)的縫洞儲層預測，獲得了較好效果。

在VMD獲得IMF的基礎(chǔ)上，運用TK/EDO方法求取瞬時振幅、瞬時頻率可以避免希爾伯特變換引起的端點效應(yīng)，其中EDO的結(jié)果更平滑，且具有良好的抗噪性及穩(wěn)定性，性能優(yōu)于TK算子。

實際資料應(yīng)用效果表明，與經(jīng)典時頻分析方法相比，基于VMD-EDO的時頻分析方法具有更高的時頻分辨率，能有效識別隱藏在寬頻地震數(shù)據(jù)中的能量異常。結(jié)合含油氣儲層“低頻能量加強、高頻能量衰減”的特點，VMD-EDO的時頻分析方法具備良好的油氣檢測能力。