變分模態分解與包絡導數算子結合的時頻分析方法及溶洞儲層預測

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

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

0 引言

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

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

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

1 方法原理

1.1 VMD

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

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

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

(1)

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

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

(2)

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

(3)

(4)

式中n為迭代次數。對應的中心角頻率為

(5)

通過

(6)

(7)

則迭代結束，得到K個IMF。

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

1.2 包絡導數算子

1.2.1 TK算子

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

(8)

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

(9)

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

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

(10)

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

1.2.2 EDO

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

(11)

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

(12)

則包絡導數算子為[32]

(13)

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

其中

(14)

將TK算子、EDO應用于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能量計算時出現負值。圖1為信號y(t)及TK算子、EDO能量計算結果。由圖可見：對于能量的變化特征，TK算子、EDO均取得了良好的追蹤結果，但在部分極值點處TK算子出現負值，這將在進一步分離能量時出現奇異值；EDO則具有良好的非負特性。

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

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，原始數據中的有效信息將主要由一個或幾個IMF所反映，選取包含儲層信息最多的IMF進行分析。

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

(19)

(20)

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

(21)

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

2 模型驗證

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

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

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

表1 模型參數

3 實際資料應用分析

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

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

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

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

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

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

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

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

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

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

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

4 結論與認識

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

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

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