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

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

0 引言

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

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

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

1 方法原理

1.1 VMD

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

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

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

(1)

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

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

(2)

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

(3)

(4)

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

(5)

通過

(6)

(7)

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

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

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

1.2.1 TK算子

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

(8)

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

(9)

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

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

(10)

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

1.2.2 EDO

對(duì)于信號(hào)瞬時(shí)能量的求取，典型的方法是以振幅值的平方進(jìn)行量化，或者以包絡(luò)的形式

(11)

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

(12)

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

(13)

式(13)結(jié)合了信號(hào)的頻率變化信息與隨時(shí)間變化的包絡(luò)，稱為包絡(luò)導(dǎo)數(shù)算子(EDO)。EDO的定義形式與TK算子較相似，只有第二項(xiàng)存在區(qū)別。O′Toole等[32]詳細(xì)對(duì)比了TK算子、EDO用于x(t)=Acos(ω0t+φ)(A為振幅，ω0為初始角頻率，φ為相位)或x(t)=Aertcos(ω0t+φ)(r為調(diào)節(jié)系數(shù)，用于描述振幅的變化)形式的簡(jiǎn)單信號(hào)時(shí)的特性，均取得較好效果。對(duì)于線性組合信號(hào)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)]包含一個(gè)等于信號(hào)y(t)兩個(gè)分量頻率之差的項(xiàng)；ψ[y(t)]則包含了額外的附加項(xiàng)，這將在TK能量計(jì)算時(shí)出現(xiàn)負(fù)值。圖1為信號(hào)y(t)及TK算子、EDO能量計(jì)算結(jié)果。由圖可見：對(duì)于能量的變化特征，TK算子、EDO均取得了良好的追蹤結(jié)果，但在部分極值點(diǎn)處TK算子出現(xiàn)負(fù)值，這將在進(jìn)一步分離能量時(shí)出現(xiàn)奇異值；EDO則具有良好的非負(fù)特性。

圖1 信號(hào)y(t)(上)及TK算子、EDO能量(下)計(jì)算結(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的時(shí)頻分析方法

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

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

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

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

(19)

(20)

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

(21)

最終可以獲得IMF的時(shí)頻分布，進(jìn)一步可以分離得到不同頻段的高精度瞬時(shí)譜。

2 模型驗(yàn)證

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

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

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

表1 模型參數(shù)

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論與認(rèn)識(shí)

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

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

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