徐天吉 ，程冰潔 ，牛雙晨，秦正曄，王貞貞

（1. 電子科技大學(xué)資源與環(huán)境學(xué)院，成都 611731；2. 電子科技大學(xué)長三角研究院（湖州），浙江湖州 313000；3. 成都理工大學(xué)“油氣藏地質(zhì)及開發(fā)工程”國家重點(diǎn)實(shí)驗(yàn)室，成都 610059；4. 成都理工大學(xué)“地球探測與信息技術(shù)”教育部重點(diǎn)實(shí)驗(yàn)室，成都 610059；5. 四川中成煤田物探工程院有限公司，成都 610072）

0 引言

平直、蛇曲、辮狀等類型多樣的河流相，多年前就已經(jīng)成為油氣領(lǐng)域的研究熱點(diǎn)。在俄羅斯西西伯利亞、美國懷俄明州及中國的鄂爾多斯盆地、四川盆地等區(qū)域，已經(jīng)發(fā)現(xiàn)了許多河流相大型和中型油氣田。但是隨著中國油氣勘探開發(fā)程度的不斷深化，目前河流相資源目標(biāo)已經(jīng)由大型、中型逐漸轉(zhuǎn)向微型油氣藏。微型古河道廣泛分布在陸相與海陸過渡相的淺、中、深等各類地層，油氣儲量豐富，勘探潛力巨大，但具有河道寬度窄、砂體薄、油氣聚集復(fù)雜等不利特征[1-2]。此外，隨著埋藏深度的增加，地震反射表現(xiàn)出能量弱、低頻為主、頻寬窄等不利特征，垂向與橫向分辨率降低，制約了窄古河道、薄儲集層及其他油氣信息的精確提取。微型古河道的識別，已經(jīng)成為增加井位部署、水平井鉆采等勘探開發(fā)風(fēng)險的突出難題。

針對微型古河道的識別難題，目前主要從地震采集、處理和解釋等角度著手，利用高分辨地震數(shù)據(jù)獲取古河道的橫向?qū)挾?、垂向厚度等關(guān)鍵參數(shù)和空間展布特征。就高分辨地震采集而言，盡管在高密度、“兩寬一高”（寬方位、寬頻帶、高密度）等先進(jìn)觀測系統(tǒng)下采集到的地震資料，提升了淺中層油氣藏預(yù)測與描述精度，但是面對深層地質(zhì)目標(biāo)，至今未能突破震源能量衰減快、檢波器難以接收到高頻有效信號等激發(fā)與接收方面的“瓶頸”問題。就高分辨地震資料處理與解釋而言，盡管形成了反Q補(bǔ)償[3]、譜整形[4]、頻寬拓展[5]、薄層反射系數(shù)反演[6]、低頻擴(kuò)展深度學(xué)習(xí)[7]、高分辨阻抗反演[8-9]等方法，可以使地震數(shù)據(jù)“無中生有”地出現(xiàn)更多的反射軸、波阻抗界面增加等現(xiàn)象，但存在信噪比降低、波組有效性評估難、甚至出現(xiàn)假反射界面等特征，必然降低地質(zhì)預(yù)測的可靠性。同時，埋藏越深、規(guī)模越小的微型古河道的識別難度越大；尤其是極窄、流向復(fù)雜的微型古河道，需要從地震波的傳播能量、反射頻率等角度出發(fā)，不斷探索突破地震分辨率極限的有效方法，才能逐步解決微型古河道的識別難題。

識別微型古河道的關(guān)鍵，是基于高分辨時頻分析方法提取地震信號的高分辨時頻響應(yīng)特征，利用地震波的能量、頻率、相位等信息，有效刻畫微型古河道的寬度、流向等特征。目前，反Q補(bǔ)償、頻寬拓展等許多高分辨處理與解釋方法，均以時頻分析為基礎(chǔ)，卻未實(shí)現(xiàn)高分辨時頻響應(yīng)特征提取，致使類似微型古河道的小規(guī)模地質(zhì)目標(biāo)難以實(shí)現(xiàn)高清識別。其實(shí)，自1946年Gabor首次提出了聲波信號的頻譜圖后[10]，迅速形成了線性、雙線性和非線性 3大類核函數(shù)不同的時頻分析方法[11]。其中，Gabor變換、短時傅里葉變換（STFT）、連續(xù)小波變換（CWT）、S變換（ST）和廣義 S變換（GST）等屬于線性時頻分析方法，雖具有計算效率較高的優(yōu)勢，但多數(shù)具有窗口依賴性，且受海森堡（Heisenberg W）不確定性原理控制[12]，時間和頻率的分辨能力不能同時提升。自適應(yīng) STFT[13]、模態(tài)分解[14-15]、同步擠壓小波變換[16]、匹配追蹤（MP）[17]等屬于非線性時頻分析方法，雖具有分辨率較高的優(yōu)勢，但計算量較大，魯棒性不足，難以適應(yīng)低信噪比信號時頻分析。Wigner-Ville分布（WVD）[18]、偽Wigner-Ville分布（PWVD）、平滑偽Wigner-Ville分布（SPWVD）和Choi-Williams分布（CWD）[19]等屬于雙線性（又稱二次型）時頻分析方法，不受海森堡不確定性原理約束[12]，優(yōu)良的數(shù)學(xué)性質(zhì)為高分辨時頻分析奠定了良好的理論基礎(chǔ)。

本文基于信號最大熵準(zhǔn)則（Maximum Entropy Principle），在克服 WVD交叉項(xiàng)干擾的前提下，研究增強(qiáng) WVD頻譜聚焦特性和提升地震信號時頻分辨率的方法，以突出窄細(xì)古河道響應(yīng)特征。并采用仿真信號、正演數(shù)值模擬等理論數(shù)據(jù)開展試驗(yàn)與論證分析，以形成微型古河道識別新方法，為井位部署、水平鉆井等提供支撐。

1 方法原理

1.1 最大熵準(zhǔn)則

熵（Entropy），是指離散隨機(jī)事件或特定信息的出現(xiàn)概率，由Shannon C E在1948年首次定義，主要用于表征隨機(jī)變量或特征信息的不確定性[20]。1957年，Jaynes E T 提出了最大熵準(zhǔn)則，即在先驗(yàn)信息準(zhǔn)確的前提下，最能代表當(dāng)前信息狀態(tài)的概率分布，是熵最大的概率分布，換言之，當(dāng)符合先驗(yàn)信息的概率分布具有不確定性時，熵最大的信息分布是合理的[21]。

1967年，Burg將最大熵準(zhǔn)則引入信號頻譜分析中，形成了最大熵時頻分析Burg算法[22]。Burg認(rèn)為，熵是信號X的分布函數(shù)，僅與X的取值概率有關(guān)，而與實(shí)際取值無關(guān)。當(dāng)信號X的取值為x時，分布密度函數(shù)為p(x)，則其熵H(x)為：

此時，Burg定義信號X的功率譜熵為：

當(dāng)利用（2）式求取功率譜時，需要在H[P(f)]最大的條件下，遞推求解信號X的自相關(guān)函數(shù)，使功率譜誤差最小，分辨率獲得提升。其實(shí)，Burg算法等價于線性自回歸AR（Auto Regressive）預(yù)測模型，最大熵功率譜P(f)的計算公式[22]為：

1.2 離散Wigner-Ville分布與功率譜

1992年Boashash定義[23]連續(xù)時序信號x(t)的解析信號為z(t)時，Wigner-Ville分布（WVD）為：

然而，實(shí)測的信號并非連續(xù)信號，實(shí)際應(yīng)用時，需要對W(t,f)進(jìn)行離散化。設(shè)x(t)、z(t)的離散表達(dá)式分別為x(n)和z(n)（0 ≤n≤N-1），對x(n)作希爾伯特（Hilbert）變換，在計算出z(n)后即可獲得W(t,f)的離散表述：

據(jù)WVD的性質(zhì)可知，離散WVD是對瞬時自相關(guān)函數(shù)進(jìn)行離散傅立葉變換。離散信號的瞬時自相關(guān)函數(shù)為：

對于每一個樣點(diǎn)n，kn(l)將組成一個卷積核序列kz(n)，如（7）式。當(dāng)n=0和n=N- 1 時，卷積核序列kz(n)=kz(0)和kz(n) =kz(N- 1)最短；當(dāng)時，最長。以kn(0)為中心，對kz(n)左右補(bǔ)充0值，可將kz(n)擴(kuò)充為長度為N的卷積核序列。

對第n個采樣點(diǎn)的卷積核序列kz(n)，進(jìn)行離散傅立葉變換，即可獲得WVD瞬時功率譜，如下：

其中，Wz(n,f)的子項(xiàng)可表述為：

通過計算每一個采樣點(diǎn)n對應(yīng)的瞬時功率譜，最終可以得到離散信號x(n)的Wigner-Ville功率譜，即：

當(dāng)然，受 WVD雙線性的影響，在利用瞬時自相關(guān)函數(shù)kn(l)前后數(shù)據(jù)相乘時，將產(chǎn)生交叉項(xiàng)干擾，不利于高精度功率譜特征分析。

1.3 最大熵準(zhǔn)則約束下的Wigner-Ville功率譜計算方法

當(dāng)瞬時自相關(guān)函數(shù)kn(l)的一階導(dǎo)數(shù)為0時，求得的Wigner-Ville分布即為最大熵譜。由于一維最大熵譜與自回歸模型AR譜是等效的，故只需要求得AR模型參數(shù)、計算出 AR模型譜，就能獲得最大熵譜[24]。這樣，可以避免利用（10）式計算功率譜，繞開瞬時自相關(guān)函數(shù)前后數(shù)據(jù)相乘時產(chǎn)生的交叉項(xiàng)干擾。

基于 Burg算法求取最大熵功率譜P(f)時，（3）式與（10）式是等效的，只需基于Levinson-Durbin遞推規(guī)則[25]，在求取預(yù)測誤差E、自回歸系數(shù)a等 AR模型的參數(shù)后，代入（3）式就能獲得最大熵準(zhǔn)則約束下的WVD功率譜。Burg算法通過前向預(yù)測誤差和后向預(yù)測誤差，遞推 AR模型的自回歸系數(shù)[26-27]。將設(shè)為初始條件，遞推公式如下：

通過離散信號x(n)及其解析信號z(n)，可以計算預(yù)測誤差的初始功率，即

利用m階前向預(yù)測誤差和后向預(yù)測誤差，可以計算預(yù)測誤差的平均功率，即

據(jù)Levinson-Durbin遞推約束規(guī)則，設(shè)am(j)為M階AR模型在階次序號為m時的第j個自回歸系數(shù)，Em為m階次時的最小預(yù)測誤差，則其遞推關(guān)系為

這樣，將am( 1)，am( 2)，…，am(m- 1)和Em代入（3）式，就能計算出 WVD最大熵功率譜。WVD最大熵功率譜在時間域和頻率域的分布，即為Wigner-Ville最大熵時頻譜（簡稱MEWVD）。

概括起來，Burg算法的WVD最大熵功率譜的計算步驟[28]主要包括：①設(shè)m=1，計算P0(f)、EF,0(n)與EB,0(n)；②計算am(m)；③計算階次為m時的第j個自回歸系數(shù)④計算Pm(f)；⑤計算階次為m時的EF,m(n)和EB,m(n)；⑥設(shè)m=m+1，重復(fù)②—⑤，直到Pm(f)不再明顯減小時，利用Levinson-Durbin遞推關(guān)系，可以獲得AR模型的參數(shù)，最終計算出WVD最大熵功率譜。

在時間域，MEWVD并未將瞬時自相關(guān)函數(shù)前后數(shù)據(jù)相乘，而僅使用當(dāng)前n點(diǎn)附近的數(shù)據(jù)，避免了交叉項(xiàng)干擾；在頻率域，受最大熵準(zhǔn)則約束，使每個n點(diǎn)對應(yīng)的功率譜更精確，能量更聚焦。同時，在最大熵的前提下，采用Levinson-Durbin約束遞推未知瞬時自相關(guān)函數(shù)，使計算的功率譜誤差減小，提升了信號時頻分辨能力。

需要指出的是，最大熵準(zhǔn)則在數(shù)學(xué)、經(jīng)濟(jì)、信息等領(lǐng)域被廣泛應(yīng)用，主要優(yōu)勢在于其不僅保留了信息分布的不確性時，同時將預(yù)測風(fēng)險降到了最低。在油氣領(lǐng)域，常常面對橫向較窄、縱向較薄的古河道、砂壩、點(diǎn)礁等微型地質(zhì)體，地震波衰減、調(diào)諧、噪聲等因素增加了微型地質(zhì)體的預(yù)測難度，計算高精度、高聚焦性的地震功率譜，是識別微型地質(zhì)體的重要基礎(chǔ)。MEWVD在熵最大時計算最能代表微型地質(zhì)體的地震頻譜響應(yīng)特征，不僅能實(shí)現(xiàn)時頻特征的分類取優(yōu)，而且可以在噪聲環(huán)境下提高信號分析的魯棒性[29]，有利于提升微型地質(zhì)體識別的可靠性。

2 方法實(shí)驗(yàn)

2.1 仿真地震信號的MEWVD實(shí)驗(yàn)及定量分析

地震波在地下傳播時，受地層埋深、孔隙度、填充流體等因素影響，將產(chǎn)生能量衰減，引起動力學(xué)、運(yùn)動學(xué)和幾何學(xué)等特征發(fā)生變化。采用振幅、吸收系數(shù)、頻率、相位和走時等參數(shù)，可以仿真描述地震波傳播過程，如下：

上式中，仿真地震信號xF(t)分別由衰減分量x1(t)和平穩(wěn)分量x2(t)疊加而成。其中，x1(t)包含了吸收系數(shù)、振幅、頻率、相位和走時等參數(shù)，x2(t)不包含吸收系數(shù)，傳播能量將不會產(chǎn)生衰減。

為了對比分析地震波傳播過程中的時頻特征，利用xF(t)的參數(shù)，構(gòu)造出了關(guān)聯(lián)特征較強(qiáng)的高、低頻能量強(qiáng)弱不一樣的仿真地震信號。如表1，信號xF1和xF2的起始能量A1和傳播時間t相同，xF1中A2、f2、j2為0，表示xF1僅由衰減分量構(gòu)成，且其A1、f1和Q均較大，類似強(qiáng)能量、高頻和強(qiáng)衰減的地震信號。xF2包含了xF1，僅通過改變f2，在xF2的基礎(chǔ)上增加了低頻、弱振幅的平穩(wěn)分量，意味著xF2在t時間內(nèi)的傳播能量、振幅等將隨頻率變化而改變。如圖 1a和圖 1b中第 1列，分別顯示了xF1和xF2的波形特征。可見，xF1隨著t的變化，振幅減弱，能量衰減很快；xF2隨著t的變化，振幅減弱，能量也在衰減，但衰減分量的控制作用越來越弱，在信號傳播約0.25 s之后，信號主要受平穩(wěn)分量控制，且由于f2差異，使xF1振蕩周期短、波形“瘦”，xF2振蕩周期長、波形“胖”。

表1 仿真地震信號的參數(shù)描述

同時，對比不同的時頻分析方法，發(fā)現(xiàn) MEWVD具有非常優(yōu)良的時頻聚焦特性和分辨率。如圖1a和圖1b，第 2列顯示了仿真地震信號的快速傅里葉變換（FFT）頻譜，揭示信號xF1的能量峰值聚集在80 Hz，xF2的能量峰值聚集在10 Hz和80 Hz。第3至10列，分別顯示了 STFT、GST、CWT、WVD、PWVD、SPWVD、CWD和 MEWVD等方法的時頻效果。對比分析，盡管xF1和xF2的最強(qiáng)能量均聚集在 80 Hz高頻位置，但是，不同的時頻分析方法顯示的能量聚焦程度與頻率分辨能力存在差異，且MEWVD聚焦效果與分辨率最優(yōu)、WVD次之、STFT最弱。GST和CWT的時頻分辨能力相似，對80 Hz高頻強(qiáng)振幅與10 Hz低頻弱振幅均能有效識別，但能量聚焦程度、時間分辨能力等較MEWVD差。此外，WVD時頻特征還存在交叉項(xiàng)干擾，PWVD雖在一定程度上壓制了交叉項(xiàng)干擾，但能量聚焦性減弱，降低了xF2的低頻、弱振幅信息的分辨能力；SPWVD和CWD基本壓制了交叉項(xiàng)干擾，但也幾乎無法有效識別xF2中低頻、弱振幅響應(yīng)，較 PWVD的分辨率損失更嚴(yán)重。

圖1 仿真地震信號及不同時頻分析方法呈現(xiàn)出的頻譜響應(yīng)差異

2.2 窄薄數(shù)值模型MEWVD實(shí)驗(yàn)及效果分析

為了進(jìn)一步證實(shí) Wigner-Ville最大熵時頻譜（MEWVD）的頻譜聚焦性和時頻高分辨能力，采用表2所示參數(shù)，設(shè)計出二維數(shù)值模型，包含了埋深、寬度和厚度各異的窄薄地質(zhì)體，其“透鏡狀”的地質(zhì)特征與微型古河道的橫截面類似?；诙S地震波動方程正演數(shù)值模擬方法，以35 Hz的Rick子波作為點(diǎn)震源激發(fā)，以10 m采樣間距、1 ms采樣率接收，獲得了正演模擬記錄。模型④—⑨的砂體厚度逐漸變薄、寬度逐漸變窄；模型④—⑥較寬、厚度大于 1/4波長（約30 m），頂?shù)捉缑娴姆瓷渫噍S明顯，而⑦—⑨較窄、厚度小于1/4波長，地震記錄表現(xiàn)為諧波反射，且⑦和⑧表現(xiàn)出強(qiáng)調(diào)諧反射，極薄、極窄的⑨則表現(xiàn)為弱調(diào)諧反射（見圖2）。

表2 數(shù)值模型參數(shù)描述

圖2 窄薄各異模型的單道模擬地震信號及MEWVD等頻譜特征

分別抽取代表模型④—⑨的單道地震信號，計算MEWVD后，可分析模型的時頻差異。如圖2，顯示了模型④—⑨的單道地震信號及MEWVD頻譜特征，可見，窄薄各異的數(shù)值模型的反射信號，MEWVD均表現(xiàn)出了較強(qiáng)的頻譜聚焦性和高分辨特征。其中，模型④—⑥較寬、較厚，單道地震信號表現(xiàn)出明確的頂?shù)捉绶瓷?，MEWVD在頻率軸方向表現(xiàn)出較窄的頻譜分布。模型⑦—⑨較窄、較薄，單道地震信號僅顯示出了頂界面的反射，底界面受波場調(diào)諧作用而難以識別。MEWVD在頻率軸方向表現(xiàn)出較寬的頻譜分布，且厚度越薄、頻譜分布范圍越寬。

在單頻剖面上，MEWVD的頻譜聚焦性非常顯著，不同頻率的功率譜展現(xiàn)出了不同的時頻分辨率。如圖3，顯示了正演地震信號的瞬時振幅及25，35，45 Hz的MEWVD頻譜特征。其中，圖3a顯示瞬時振幅受時頻分辨能力的局限，僅較厚的模型④表現(xiàn)出了頂?shù)捉缑媲逦摹半p軸”強(qiáng)瞬時振幅特征，而厚度逐漸減薄的模型⑤—⑨卻表現(xiàn)為“單軸”強(qiáng)瞬時振幅分布，連厚度大于1/4波長的模型⑤和⑥的頂、底界面也難以有效識別。圖3b—圖3c分別顯示了3種頻率的MEWVD頻譜分布，且均比瞬時振幅具有更高的分辨率，但不同頻率的 MEWVD頻譜響應(yīng)各異。其中，45 Hz的MEWVD頻譜分辨率最高，可以準(zhǔn)確識別窄薄模型⑤和⑥的頂?shù)捉缑婕皺M向?qū)挾取?5 Hz的MEWVD頻譜分辨率次之，也能夠有效識別⑤和⑥的頂?shù)捉缑婕皺M向?qū)挾龋瑢O窄、極薄模型⑧和⑨的頻譜聚焦響應(yīng)優(yōu)良，且受不同寬度的影響，二者的頻譜分布厚度相近、寬度差異顯著。進(jìn)一步分析極窄、極薄模型⑧和⑨，二者厚度相同、寬度差 5倍，瞬時振幅橫向差異不明顯，MEWVD頻譜異常顯著。尤其在圖3d中，45 Hz高頻聚焦性極好，顯示出了厚度一樣、寬度差 5倍的MEWVD頻譜異常特征。同時，雖然25 Hz的MEWVD頻譜分辨率總體不如35 Hz和45 Hz的頻譜，但是對極窄、極薄模型⑧和⑨的頻譜聚焦響應(yīng)卻有明顯差異。其中，35 Hz和45 Hz的頻譜表現(xiàn)為“團(tuán)狀”現(xiàn)象，而25 Hz的頻譜卻表現(xiàn)出“層狀”頻譜聚焦現(xiàn)象，表明Wigner-Ville最大熵譜與地震信號的強(qiáng)弱無關(guān)，當(dāng)調(diào)諧作用使模型頂?shù)捉缑娴墓β首V分布不確定時，熵最大的功率譜即為最合理的分布。

圖3 正演地震信號的瞬時振幅及25，35，45 Hz的MEWVD頻譜特征

總之，仿真地震信號和窄薄正演數(shù)值模型實(shí)驗(yàn)表明，MEWVD具有良好的頻譜聚焦性和時頻高分辨能力。MEWVD能夠精確表征信號的時頻分布，且不同頻率的MEWVD能夠準(zhǔn)確獲取不同尺度窄薄模型的頻譜響應(yīng)；同時，由于Wigner-Ville最大熵譜僅與頻譜分布概率相關(guān)，當(dāng)?shù)刭|(zhì)體較窄或薄度小于1/4波長而發(fā)生調(diào)諧作用時，利用地震反射同相軸不能直接識別，熵最大的功率譜即為地質(zhì)體的最合理的分布空間。因此，綜合不同頻率的MEWVD頻譜特征，可以識別類似窄薄模型的微型古河道的空間分布。

3 應(yīng)用案例

四川盆地是一個多旋回富氣沉積盆地，截至2020年，在震旦系、寒武系、志留系、石炭系、三疊系、侏羅系等27套層系中，已經(jīng)發(fā)現(xiàn)了構(gòu)造、構(gòu)造-巖性、巖性-構(gòu)造和巖性4類氣藏，形成了四川盆地東北部、西部、南部和北部 4大工業(yè)氣區(qū)[30-31]。其中，川西氣區(qū)的侏羅系沙溪廟組氣藏，主要分布在西鄰龍門山推覆構(gòu)造帶、東接川中隆起區(qū)的川西坳陷。在川西東斜坡地區(qū)，發(fā)現(xiàn)了埋深小于3 100 m的大型巖性-構(gòu)造氣田——中江氣田。

中江氣田面積為2 350 km2，區(qū)內(nèi)侏羅系沙溪廟組以淺水三角洲沉積為主，在三角洲平原亞相中發(fā)育分流河道、分流間灣、河口壩、決口扇等沉積微相。其中，分流河道主要發(fā)育細(xì)—中粒巖屑長石砂巖、長石巖屑砂巖、巖屑砂巖和巖屑石英砂巖等砂巖儲集層，平均孔隙度達(dá) 8.94%，平均滲透率達(dá)0.5×10-3μm2。目前，針對資源豐富的沙溪廟組氣藏，已投產(chǎn)工業(yè)氣井約150口，日產(chǎn)氣超過 220′104m3。然而，隨著開發(fā)程度的深化，微型分流河道氣藏更加隱蔽，儲集層以條帶狀疊置薄層河道砂體為主，橫向?qū)挾日v向厚度薄，精細(xì)刻畫難度極大。需要探索窄細(xì)微型古河道的精確識別方法，為氣藏描述、水平鉆井、水力壓裂等提供支撐。

中江氣田侏羅系中統(tǒng)沙溪廟組劃分為上（J2s1）、中（J2s2）、下（J2s3）共 3段；其中，J2s3段由 J2s31、J2s32和 J2s33組成。依據(jù)古河道和砂巖沉積特征，J2s33層被細(xì)分為J2s33-1和J2s33-2小層。目前，J2s33-2小層是主力產(chǎn)氣層，層內(nèi)網(wǎng)狀和條帶狀古河道主要沿北東流向南西，延伸長度18～58 km，寬度0.2～3.5 km，砂厚6.5～45.2 m。儲集層呈席狀、箱狀和鐘型分布，表現(xiàn)為低頻、“亮點(diǎn)”、“空白”等地震反射特征，沿古河道方向，地震反射總體較連續(xù)、穩(wěn)定，垂直河道方向，呈“透鏡狀”或弱反射異常。

近幾年，隨著中江氣田 J2s33-2小層氣藏的深入開發(fā)，寬度小于500 m、厚度小于35 m的窄薄古河道砂巖成為重點(diǎn)目標(biāo)，微型古河道的精確識別也成為氣藏開發(fā)的關(guān)鍵環(huán)節(jié)。然而，受埋藏深度、寬度、厚度、巖性、物性等多種因素的影響，不同的微型古河道呈現(xiàn)出明顯的地震反射差異。如圖4a圖顯示J2s33-2小層的地震主頻約30 Hz，頻帶8～70 Hz；圖4b顯示E井J2s33-2小層古河道出現(xiàn)砂體GR低值異常，地震反射為局部雙波峰“透鏡狀”異常，預(yù)測河道寬約354 m、厚約35 m；圖4c顯示D井J2s33-2小層古河道出現(xiàn)GR低值異常，預(yù)測河道砂巖寬約310 m、厚約30 m，較E井河道更窄、更薄，“透鏡狀”反射更弱；圖4d顯示C井J2s33-2小層古河道出現(xiàn)GR低值異常，預(yù)測河道砂巖寬約248 m、厚約30 m，與D、E井相比較，“透鏡狀”異常消失，呈弱反射響應(yīng)。顯然，通過對比圖4中的C、D、E井微型古河道地震反射特征，隨著寬度變窄、厚度減薄，河道隱藏性增強(qiáng)，地震反射越弱，識別難度越大。

圖4 目標(biāo)地層地震振幅譜及典型井微型分流河道地震反射剖面

采用Wigner-Ville最大熵時頻譜計算方法，提取了中江氣田沙溪廟組的MEWVD頻譜特征，發(fā)現(xiàn)不同寬度和厚度的微型古河道呈現(xiàn)出了不同的MEWVD異常，揭示了窄細(xì)古河道的空間展布。如圖5a顯示的瞬時振幅剖面上，E井J2s33-2小層古河道出現(xiàn)瞬時振幅強(qiáng)異常，異常范圍較大，聚焦性較差，分辨率較低；圖5b—圖5d顯示的MEWVD頻譜剖面上，E井J2s33-2小層古河道出現(xiàn)強(qiáng)異常，異常范圍隨著頻率的增加而減小，聚焦能力增強(qiáng)，分辨率提高。比較瞬時振幅與25 Hz、35 Hz的MEWVD頻譜特征，在45 Hz的MEWVD剖面上，分辨能力更高，縱向頻譜被分開，橫向頻譜收窄，頻譜異常由“團(tuán)狀”變化為“雙層狀”，與砂體GR低值異常更加吻合，有效地揭示出了E井J2s33-2小層35 m薄古河道砂巖的頂?shù)捉缑妗？梢姡肕EWVD頻譜異常，能夠有效刻畫窄細(xì)古河道的縱向分布。

圖5 E井瞬時振幅及25，35，45 Hz的MEWVD頻譜特征

同時，利用MEWVD頻譜異常，還能刻畫微型古河道的橫向展布。如圖6a瞬時振幅高值異?？坍嫵隽讼鄬^寬的Ⅰ—Ⅳ、Ⅵ、Ⅷ號主河流的沿層展布；但在黑色框示區(qū)域，河道很窄、河道砂巖很薄，地震反射能量較強(qiáng)，分辨率低，利用瞬時振幅無法識別窄細(xì)古河道的空間位置、寬度和流向。在圖 6b中，利用MEWVD方法獲得的25 Hz頻譜出現(xiàn)高值異常，不僅刻畫出了較寬的Ⅰ—Ⅳ、Ⅵ、Ⅷ號主河流的沿層展布，而且中低值異常刻畫出了更窄的Ⅴ和Ⅶ號河道的空間展布；尤其在左側(cè)黑色框內(nèi)，Ⅴ和Ⅶ號河流的寬度和流向非常清晰。與Ⅰ—Ⅷ微型河道比，Ⅸ號河道更窄，且流向彎曲變化較大，隱蔽性更強(qiáng)，利用瞬時振幅和25 Hz頻譜難以識別，而在35 Hz和45 Hz的MEWVD頻譜特征中，Ⅸ號曲流河道的沿層特征卻非常清晰。由圖6可見，MEWVD較瞬時振幅分辨更高，且不同頻率的MEWVD頻譜特征，可以識別不同寬度的窄細(xì)河道。當(dāng)然，由于任何物質(zhì)都有相應(yīng)的敏感頻段，河道也不例外，Ⅴ和Ⅶ號河道對低頻更敏感，故在35 Hz和45 Hz的MEWVD頻譜特征中反而不如25 Hz刻畫更清晰（見圖6黑色方框），但整體規(guī)律仍然是高頻頻譜刻畫效果更佳。

圖6 瞬時振幅及25，35，45 Hz的MEWVD分別刻畫的J2s33-2小層微型古河道的展布

采用RGB融合（紅、綠、藍(lán)混合生成其他色彩）顯示方式，將不同頻率的MEWVD頻譜進(jìn)行融合，能更清晰的展示不同寬度和厚度河道的空間分布。如圖7a顯示了采用 RGB融合方式，將 25，35，45 Hz的MEWVD頻譜進(jìn)行融合，有效刻畫出J2s33-2小層Ⅰ—Ⅸ號古河道的展布；尤其是具有極窄特征的Ⅴ、Ⅶ和Ⅸ號河道，以及Ⅷ號河道以東、Ⅵ號河道以西的更多微型古河道（紅色箭頭所指），也精確地揭示出了其寬度和流向等空間信息。當(dāng)然，結(jié)合區(qū)內(nèi)沉積構(gòu)造背景，利用巖心分析與測井解釋資料，可以證實(shí)采用MEWVD

頻譜異常識別微型古河道的可靠性。如圖7b所示，在縱向上，D井主要發(fā)育分流河道、分流間彎、河口壩等三角洲前緣沉積微相。其中，J2s33-2小層以分流河道沉積微相為主，聲波時差和伽馬測井曲線表現(xiàn)出低值響應(yīng)特征，深側(cè)向和淺側(cè)向電阻率曲線表現(xiàn)出高值異常。在分流河道中部2 869.57～2 871.49 m井段的鉆井取心顯示，古河道水平層理發(fā)育，在較強(qiáng)的水動力作用下，沉積形成了厚度約30 m的致密砂巖儲集層。受巖性、物性、厚度等因素影響，儲集層與圍巖的阻抗差異不大，在圖4c所示的“透鏡狀”河道反射較弱，但在圖6和圖7卻具有較突出的MEWVD頻譜異常?？梢?，綜合地震反射、鉆井取心、測井響應(yīng)等信息，利用MEWVD頻譜異?？坍嬑⑿凸藕拥谰哂休^高可靠性。

圖7 中江氣田J2s33-2小層RGB融合MEWVD刻畫河道分布及D井沉積相測井解釋

此外，MEWVD頻譜異常還可以刻畫斷層分布特征，但斷層與古河道具有顯著的MEWVD差異。如圖6d和圖7a顯示，中江氣田沉積構(gòu)造環(huán)境整體比較穩(wěn)定，在試驗(yàn)工區(qū)南部僅發(fā)育少量近東西向的斷層，裂縫不發(fā)育。斷層的MEWVD頻譜異常表現(xiàn)為上盤與下盤邊界頻譜較強(qiáng)、中間的頻譜較弱。斷層與古河道雖皆呈細(xì)長條帶狀，但古河道邊界頻譜較弱、中間頻譜較強(qiáng)。可見，斷層與古河道的MEWVD頻譜異常雖有相似之處，卻也存在顯著差異。本質(zhì)上，二者的共性與差異皆由地質(zhì)環(huán)境決定，這些地震響應(yīng)被MEWVD異常精確的刻畫了出來。

總之，MEWVD頻譜聚焦性較強(qiáng)，時頻分辨率較高，且不同頻率的MEWVD對窄細(xì)不同的河道具有不同程度的敏感性，有效的揭示了中江氣田沙溪廟組微型古河道的寬度、砂巖厚度和流向等空間信息。

4 結(jié)論

熵最大時的頻譜特征，最能代表地震信號的頻譜分布狀態(tài)，MEWVD不僅避免了 WVD產(chǎn)生的交叉項(xiàng)干擾，而且最大程度的增強(qiáng)了頻譜聚焦性，有效提升了地震信號的時頻分辨率。

仿真地震信號和窄薄模型正演模擬信號試驗(yàn)揭示，較STFT、GST、CWT、WVD、PWVD、SPWVD、CWD等方法，MEWVD能更加精確的提取信號頻譜特征；尤其是針對窄薄各異的數(shù)值模型，由于最大熵的差異而表現(xiàn)出了特殊的頻譜響應(yīng)，采用不同頻率的MEWVD能夠準(zhǔn)確的識別尺度各異的窄薄地質(zhì)體。

基于最大熵準(zhǔn)則與 Wigner-Ville分布的微型古河道識別方法，充分發(fā)揮了MEWVD頻譜聚焦性較強(qiáng)、時頻分辨率較高等優(yōu)勢，有效的識別出了中江氣田不同尺度的微型古河道，且不同頻率的MEWVD對窄細(xì)各異的河道具有不同程度的敏感性，準(zhǔn)確地刻畫了微型古河道的寬度、砂巖厚度和流向等空間信息。

符號注釋：

a——自回歸系數(shù)，又稱為預(yù)測誤差因子，無因次；am——序號為m的自回歸系數(shù)，無因次；1A——衰減分量的振幅，dB；A2——平穩(wěn)分量的振幅，dB；E——預(yù)測誤差，W/Hz；EB，m——序號為m的后向預(yù)測誤差，W/Hz；EF，m——序號為m的前向預(yù)測誤差，W/Hz；f——信號的頻率，Hz；f1——衰減分量的頻率，Hz；f2——平穩(wěn)分量的頻率，Hz；GR——自然伽馬，API；H——熵函數(shù)，無因次；i——虛數(shù)單位；j——自回歸系數(shù)am的序號，無因次；kn——自相關(guān)函數(shù)，無因次；kz——kn的卷積核，無因次；l——離散信號擴(kuò)展采樣點(diǎn)序號，無因次；m——am、EF，m、EB，m、wn等函數(shù)的序號，無因次；M——AR模型的階數(shù)，無因次；n——離散信號采樣點(diǎn)序號，無因次；N——離散信號長度，無因次；p——分布密度函數(shù)，無因次；P——功率譜，W/Hz；P0——初始功率譜，W/Hz；Pm——平均功率譜，W/Hz；Q——吸收系數(shù)，無因次；Rlld——深側(cè)向電阻率，W·m；Rlls——淺側(cè)向電阻率，W·m；t——時間，s；W——連續(xù)WVD函數(shù)，W/Hz；Wz——離散WVD函數(shù)，W/Hz；wn——Wz的子項(xiàng)，W/Hz；X——信號，無因次；x——某時段信號，dB；x1——衰減分量，dB；x2——平穩(wěn)分量，dB；Dt——時間采樣率，s；xF——仿真信號，dB；j1——衰減分量的相位，rad；j2——平穩(wěn)分量的相位，rad；t——時延，s；z——離散信號，dB。