基于WVD-MEM的高分辨河道識別方法

牛雙晨 程冰潔

(①成都理工大學油氣藏地質及開發工程國家重點實驗室，四川成都 610059；②成都理工大學地球勘探與信息技術教育部重點實驗室，四川成都 610059；③四川中成煤田物探工程院有限公司，四川成都 610072)

0 引言

眾所周知，河道砂體往往是油氣藏的重要載體，隨著勘探開發不斷深化，埋深更深、寬度更窄、厚度更薄的河道逐漸成為油氣勘探領域的研究熱點。針對河道空間展布特征識別的難題，業界廣泛運用振幅、頻率、相位等地震屬性和速度、密度、阻抗等反演參數開展研究。然而，受河道埋藏深度大、縱橫向分布復雜和地震資料主頻偏低、頻帶窄、空間分辨率不足等因素的制約，傳統方法對河道的識別精度往往不高，識別效果欠佳，不利于氣藏空間分布預測、儲量估算、水平井鉆井等勘探開發和工程作業。

近幾年來，時頻分析技術在地震勘探領域發揮著越來越重要的作用[1]，被廣泛應用于主振幅和主頻率等參數計算[2]、地震與測井信號的多分辨率分析[3]、儲層預測[4]和吸收補償[5]等。目前，地震信號時頻分析技術也成為河道識別的重要手段，應用越來越廣泛，針對地下河道沉積相的空間延展、厚度和寬度變化等特征，基于地震波在河道介質中的傳播特性，利用地震信號開展時頻分析，在時間域和頻率域同時研究振幅、功率譜、能量、相位等地震屬性。如，陳杰等[6]分析了短時傅里葉變換(STFT)的適用條件，并應用于河道砂體識別；楊道慶等[7]通過離散傅里葉變換，利用分頻技術識別河道砂體；馮金義[8]應用時頻分析技術提取了不同頻率的瞬時振幅沿層切片，以檢測河道砂體的空間分布。

1946年，Gabor[9]最早提出了利用聲波信號頻譜圖等時頻分析技術。根據核函數的不同，把時頻分析法分為線性、雙線性和非線性三大類[10]。其中，線性方法包括Gabor頻譜圖、STFT[11]、連續小波變換(CWT)[12]、S變換(ST)[13]和廣義S變換(GST)[14-16]等，該類方法計算速度較快，但多具有窗口依賴性，因而無法同時提高時間分辨率與頻率分辨率。雙線性方法又稱為二次型時頻分布，包括WVD[17]、平滑偽Wigner-Ville分布(SPWVD)[18]和Choi-Williams分布(CWD)[19]等，以WVD為例，本質屬于信號自相關系數的傅里葉變換，不受Heisenberg測不準原理[20]的制約，具有理論上優良的數學性質[21]。非線性方法包括基于數據驅動的自適應分解方法，如自適應STFT[22-23]，相比STFT窗口長度靈活可變；Huang等[24]提出的經驗模態分解(EMD)、集合經驗模態分解(EEMD)以及完全集合經驗模態分解(CEEMD)[25-26]，時頻分辨率高，但仍存在受信噪比影響大和計算效率低等問題；融合CWT與EMD的同步擠壓小波變換[27]，以及融合STFT、ST的同步擠壓變換(SST)[28-29]，分辨率和計算效率都很高，缺點是對低信噪比信號分析魯棒性不足，且嚴格意義上同步擠壓技術只是對時頻變換結果的譜圖重排處理。此外，還有基于信號稀疏表示的非線性方法，如基于Gabor原子的匹配追蹤分解[30]等。

與非線性方法等其他時頻分析方法相比，WVD具有極高的時頻聚焦性，然而WVD存在交叉項干擾的缺陷。目前，最常用的消除交叉項的方法是核函數平滑法，如偽Wigner-Ville分布(PWVD)、SPWVD以及CWD等，然而這些方法均以犧牲時頻分辨率為代價，且對交叉項干擾抑制不徹底。有學者也提出了許多消除交叉項的新方法，劉希康等[31]提出聯合STFT和WVD消除交叉項，結合STFT不存在交叉項干擾和WVD時頻聚焦性高的優點，既抑制了交叉項干擾，又保留了高分辨率特性，但仍存在依賴窗函數和地震信噪比等問題。孔慧芳等[32]提出基于分數階傅里葉變換(FRFT)抑制WVD的交叉項，在有效抑制交叉項干擾的同時具有一定的抗噪性能，但僅限于線性調頻信號。隨后，紀永禎等[33]提出了結合貝葉斯學習算法和WVD消除交叉項，該方法在分辨率和保真性能方面具有一定優勢，但依賴雷克子波庫。與上述消除交叉項的方法相比，Zoukaneri等[34]提出的WVD-MEM具有極高的時頻分辨率，計算簡便、抗噪性能好，不依賴雷克子波等各種子波庫，特別是對交叉項干擾有良好的抑制效果，且適用于包括地震信號在內的非平穩隨機信號。

本文在前人研究基礎上，利用非線性調頻(FM)信號、合成地震信號和實際地震數據研究WVD-MEM消除交叉項干擾的方法。在對比分析目前廣泛應用的STFT、CWT、GST、SPWVD、CWD等時頻分析方法的基礎上，針對調頻信號、仿真信號、單道實測地震信號和楔狀正演模型，開展WVD-MEM薄層分辨率定量分析和抗噪性能測試，并成功應用于四川盆地侏羅系沙溪廟組氣藏的河道識別，基于WVD-MEM方法實現了對實測地震信號功率譜、瞬時振幅等高分辨地震屬性的提取，獲取了河道空間展布、寬度、厚度等重要信息，為氣藏勘探開發提供了方法支撐。

1 方法原理

1.1 Wigner-Ville分布(WVD)

Wigner-Ville給出的WVD定義為

(1)

式中：z(t)為信號x(t)的解析信號(解析信號包括原始信號和希爾伯特變換)；上標“*”表示共軛轉置；t是時間；f是頻率；τ是時延；j為虛數單位。

在時間域進行離散處理，信號x(n)、z(n)的關系被定義為

z(n)=x(n)+jH[x(n)]

(2)

式中：H[x(n)]為信號x(n)的希爾伯特變換；n為采樣點整數，且n∈[0，N-1]；N為最大采樣點數。

利用解析信號z(n)可以生成信號自相關系數協方差矩陣C

C=zzT

(3)

式中：C為N×N階矩陣；上標“T”表示復共軛轉置。

協方差矩陣C具有沿主對角線兩端對稱及復數實部相等、虛部相反的性質，求解C的一半即可完整求出矩陣。矩陣C沿每條交叉對角線的序列，構成了離散WVD的核函數，可以表示為

K(n)={kn(-l),…,kn(0),…,kn(l)}

(4)

式中：kn(0)表示矩陣C的主對角線元素；l∈[0，N-1]。定義kn(l)為

(5)

當l=0時，中心項kn(0)也是矩陣C的主對角線元素，有

kn(0)=z(n)z*(n)n=1，…，N

(6)

W(n)={Wn[-(N-1)],…,

Wn(0),…,Wn(N-1)}

(7)

1.2 基于最大熵(MEM)原理的功率譜計算方法

最大熵(Maximum Entropy Method，MEM)原理指出，信號熵值最大的概率分布是最合理的分布。這樣，基于最大熵原理的功率譜求解與自回歸(autoregressive,AR)模型功率譜的求解是等價的[35]，基本方程可以表示為

(8)

式中：P(f)為功率譜；p代表最大熵AR模型的階數；m表示階數p的序數，取值范圍為[0，p-1]；cm是序數為m階時的自回歸系數，當m=0時c=1；Ep-1對應p-1階時的預測誤差能量；Δt表示采樣間隔。

在其他變量已知的前提下，方程求解的關鍵在于自回歸系數c和預測誤差能量E的求取。利用式(6)可得到E0，利用Burg算法[35]可直接從z(n)中求出自回歸系數c

(9)

c(m,i)=c(m-1,i)+c(m,m)c*(m-1,m-i)

i=1,…,m-1

(10)

Em=Em-1[1-c(m,m)c*(m,m)]

鋁股和芯棒在穩定狀態下的應力水平仿真結果如表9所示。由圖11和表9可看出：導線鋁股自外向內應力逐漸減小，導線內層鋁股與芯棒接觸位置應力最大。仿真結果證明導線外層所施加的徑向載荷能夠很好地傳遞給芯棒。

(11)

式中kz(m)為輸入信號自相關函數。聯合式(9)～式(11)，進入Levinson遞歸求解完整的自回歸系數c和預測誤差能量E，再由式(8)求得AR模型的最大熵功率譜。

1.3 基于WVD-MEM的高分辨時頻屬性計算方法

離散WVD核函數的累計求和，是輸入信號自相關函數的偶數項序列，表示為

(12)

Burg在Levinson遞歸方程中推導出自相關函數的計算公式

c(m,m)Em-1

(13)

聯合式(12)和式(13)，得到WVD-MEN的核函數

c(m,m)Em-1

(14)

與WVD核函數求取方法類似，根據協方差矩陣的共軛性質，只需求取右半部分的值即可求出完整的WVD-MEM的核函數。之后對其進行快速傅里葉變換，可求出WVD的最大熵瞬時功率譜，即WVD-MEM時頻分布。由于時頻分布的峰值對應瞬時振幅的平方，且峰值在時頻面的投影為瞬時頻率，因而利用峰值檢測法[36]可以獲得瞬時頻率和瞬時振幅等屬性參數。

WVD-MEM高分辨時頻屬性的主要計算步驟如下：

(1)由式(2)對原始輸入數據x(n)做希爾伯特變換，獲得解析信號z(n)。

(2)為了抑制數據變換結果的端點效應(即時頻表示平面上兩端數據隨窗口增大而消失的現象)，需對z(n)進行加窗口擴展，在z(n)首尾各添一段數據，窗口長度均為L，設定值均為0。

(15)

式中：L為窗口長度，取奇數；zL(n)為解析信號z(n)擴展后的新數據，以z(n)為中心。

(3)設定最大熵AR模型階數。

(4)由式(6)計算最大熵模型的預測誤差初值(即WVD-MEM核函數的中心項)。

(5)利用式(9)～式(11)求自回歸系數c。

(6)利用式(14)及自回歸系數c，求WVD-MEM的核函數knWM(m)。

(7)對knWM(m)做快速傅里葉變換，獲得WVD-MEM高分辨時頻屬性。

需要說明的是，由于第(2)步擴展數據時首尾加上了窗口長度L的數據，因此需對求得的WVD-MEM時頻屬性結果矩陣邊緣長度為L的數據做舍去處理。

采用上述步驟不僅可以求取地震信號的功率譜，還可以計算瞬時頻率、瞬時振幅等地震屬性。由于WVD-MEM時頻屬性計算方法結合了WVD的時頻高聚焦特性和MEM的最大概率分布兩項優勢，因此獲得的功率譜、瞬時頻率、瞬時振幅等地震屬性具有顯著的高分辨特征，可為薄儲層預測、窄河道識別及其他微型地質體研究提供方法支撐。

2 方法試驗與分辨率定量分析

2.1 調頻信號定量試驗

利用非線性調頻(FM)信號研究WVD-MEM的時頻分辨能力。圖1a顯示了正弦調頻信號x1、雙曲線調頻信號x2和二者線性疊加信號x的時間域波形特征。由圖可見，x1的振幅隨著時延增大,曲線變密，x2的振幅隨著時延減小，x的振幅隨著時延變化并不明顯，但曲線前半部分較為稀疏，后半段較密。圖1b顯示x1和x2在頻率域的變化特征，其中，x1的頻率呈高低起伏的“波浪狀”，x2的頻率隨著時延逐漸降低，且降低的幅度越來越小。

圖1 時間域和頻率域的非線性調頻信號

圖2是分別采用STFT、CWT、GST、WVD、SPWVD、CWD和WVD-MEM等7種時頻分析方法獲得的疊加信號x的時頻譜。圖2a顯示疊加信號x的實際時頻譜是x1和x2兩個分量的疊加，在時頻域被分離成兩部分，同時具有x1和x2的頻譜特征。圖2b～圖2d為采用線性時頻分析方法STFT、CWT和GST處理獲得的時頻譜，雖不存在交叉項干擾，對分量x2的頻譜聚焦性也較好(黃色箭頭所示)，但對分量x1的頻譜聚焦性相對較差(紅色箭頭所示)。WVD方法雖然時頻聚焦性好，但存在嚴重的交叉項干擾(圖2e中綠色箭頭所示)。SPWVD(圖2f)和CWD(圖2g)處理雖然消除了WVD產生的交叉項干擾，但降低了頻率分辨率，尤其是CWD，存在時窗平滑效應產生的帶狀干擾(圖2g白色箭頭所示)。相比之下，WVD-MEM方法(圖2h)不僅完整消除了WVD存在的交叉項干擾，而且保持了與實際頻率基本一致的時頻分布，清楚地顯示出x1和x2頻譜分別由低到高和由高向低的變化特征。

圖2 采用不同方法獲得的疊加信號x時頻分布特征

總之，試驗表明這7種方法均能分辨疊加信號x中的x1、x2分量，且線性時頻分析方法中，GST的聚焦性優于STFT和CWT；二次型時頻分析方法中SPWVD、CWD和WVD-MEM的聚焦性整體優于線性時頻分析法。其中WVD-MEM方法具有更顯著的高分辨率優勢。

2.2 單道仿真地震信號——非平穩信號定量試驗

利用單道仿真地震信號定量檢測WVD-MEM方法對非平穩信號的高分辨時頻分析能力，并為實際地震信號功率譜、瞬時振幅等高分辨時頻屬性的提取和應用奠定基礎。圖3a為使用5個雷克子波分量合成的單道仿真地震信號，子波分量的主頻分別為：①20Hz；②30Hz和60Hz；③40Hz、30Hz和60Hz；④60Hz、40Hz和60Hz；⑤60Hz、60Hz和60Hz。最大振幅分別為：①30dB；②20dB和-10dB；③15dB、-10dB和20dB；④30dB、20dB和-30dB；⑤20dB、40dB和20dB。信號長度(采樣總數)分別為①100；②100；③100；④100；⑤101。具體參數選取和計算方程如式(16)所示。

(16)

式中：r(f,t)為雷克子波表達式；f1～f5為5個雷克子波分量頻率。

圖3b～圖3h分別為采用STFT、CWT、GST、WVD、SPWVD、CWD和WVD-MEM等7種方法獲得的仿真信號時頻分布。由圖可見，7種方法均在0.1、0.3、0.5、0.7及0.9s附近檢測到相應的時頻能量響應，但能量分布具有明顯的差異。三種線性時頻分析方法中，GST聚焦性和分辨率最好，CWT次之，STFT最差；相比之下，二次型時頻分布的聚焦性均優于GST，但WVD法在0.2、0.4、0.6和0.8s附近出現了嚴重的干擾項(圖3e中綠色箭頭所示)；SPWVD和CWD雖然消除了圖3e中的交叉項干擾，但對于0.5、0.7和0.9s附近的交叉項干擾卻無法完全抑制(黃色箭頭所示)。而WVD-MEM則完全消除了WVD中存在的所有交叉項干擾，表現出最優的時頻聚焦性和頻率定位精度，清楚地顯示了20、30、40、60Hz等頻率位置的功率譜強弱變化。例如，在0.9s附近的仿真信號由3個60Hz主頻、中間振幅較高(40dB)、兩側振幅較低且相同(20dB)的雷克子波合成，WVD-MEM時頻分析方法不僅準確展示出子波能量變化特征，而且精確定位在60Hz頻率位置，其他方法均未能達到如此高的分辨率。

圖3 單道仿真地震信號不同方法時頻分析效果對比

2.3 單道實測地震信號——非平穩信號試驗

上述仿真信號是由與實際地震信號近似的雷克子波組成，已知子波頻率、相位、振幅等信息，而實測地震信號中這些信息是未知的，且更加復雜，二者存在明顯差異。因此，有必要針對實測地震信號開展試驗，進一步分析該方法對實際地震信號的適用性，為后續利用功率譜、瞬時頻率、瞬時振幅等地震屬性進行高分辨河道識別提供參考。圖4為單道實測地震信號及采用7種方法的時頻分析結果對比。同樣，在STFT、CWT和GST等3種線性時頻分析方法中，GST的時頻分辨率和聚焦性最優；WVD方法受到交叉項的干擾，有效的時頻特征被完全淹沒，無法進行信號分析；SPWVD和CWD方法的時頻聚焦性優于線性時頻分析方法，但交叉項干擾抑制不徹底(黃色箭頭所指)。WVD-MEM則徹底消除了WVD存在的嚴重交叉項干擾，在保持高度時頻聚焦性的同時，具有明顯更高的時頻分辨特性，對能量隨頻率和時間變化的細節刻畫得更精確。可見，針對更復雜的實際地震信號，WVD-MEM方法仍然能獲得高聚焦、高時頻分辨率的處理效果。

圖4 單道實測地震信號時頻分析效果對比

2.4 多道正演數值模擬數據的WVD-MEM時頻分辨能力定量分析

河道識別本質上是地震極限分辨率條件下的薄層識別問題。通過提高瞬時功率譜、瞬時頻率、瞬時振幅等地震屬性的分辨率，可以提高地震信號的薄層分辨率，從而實現河道識別。設計一個楔狀模型正演數值模擬地震數據(圖5a)，震源為主頻35Hz的雷克子波，檢波器的采樣間距為5m，接收長度為1.5s，共201道。在地層變薄位置(紅色與藍色圈內)，可以觀察到薄層調諧效應產生的明顯復合地震反射現象。

圖5b～圖5d分別顯示了利用GST、SPWVD和WVD-MEM三種方法提取的35Hz單頻剖面。由圖可見，在紅色圈標識的薄層調諧位置，GST、SPWVD及WVD-MEM方法的最大識別位置分別為第66、54、41道；在藍色圈標識的薄層調諧位置，GST、SPWVD及WVD-MEM方法的最大識別位置分別為第132、144、161道。顯然，WVD-MEM表現出更高的能譜聚焦性，頻譜被高度細化，增強了時頻分辨能力，能更好地消除薄層調諧效應的影響，更有利于薄層識別。

圖5 正演數值模擬記錄及采用不同方法獲得的35Hz單頻剖面

2.5 WVD-MEM的抗噪性能測試分析

實際地震信號通常含有一定程度的噪聲，因而，有必要在不同信噪比(SNR)條件下測試WVD-MEM方法的抗噪性能，以進一步明確其實用價值。圖6顯示單道實測地震信號(圖6a)及SNR分別為30、20和10噪聲條件下，應用WVD-MEM獲得的瞬時振幅譜特征。可見，隨著SNR由高向低轉變，相應的WVD-MEM處理效果與原始地震信號處理結果相比較并未發生明顯改變。由此可以證實，WVD-MEM具有良好的抗噪性能，能夠適用于實測地震信號的時頻處理。

圖6 不同SNR實測地震信號的WVD-MEM處理結果對比

綜上所述，在非線性調頻信號、合成地震信號、單道實測地震信號及多波正演數值模擬地震記錄的定量測試分析中，WVD-MEM時頻分析方法都表現出較強的能譜聚焦和譜細化特性，時頻分辨能力顯著，且抗噪性能良好。因此，當地震波在地下介質中傳播，遇到目標地質體時可能出現功率譜、瞬時振幅、瞬時頻率等時頻屬性異常，借助WVD-MEM方法可以更精確地突出時頻屬性異常特征，進而有效識別薄儲層、窄河道及其他微型地質體。

3 應用實例

圖7 過JS205井(左)和JS301井(右)河道的地震反射特征

地震數據頻譜分析對提取瞬時振幅及分頻屬性具有指導作用。首先對研究區地震數據進行頻譜分析(圖8)，可以看出，工區地震數據主頻約為35Hz，頻帶范圍為5～80Hz，優勢頻帶范圍為7～65Hz。結合前文方法與試驗研究，采用WVD-MEM方法處理工區地震數據并提取瞬時頻率、瞬時振幅和不同頻率功率譜等地震屬性，從中發現時頻屬性異常，以多視角分析出河道的異常響應特征，實現河道空間展布特征識別。

圖8 研究區頻譜分析

圖9 研究區WVD-MEM提取的不同頻段頻譜剖面

圖10 利用WVD-MEM方法提取的時頻屬性識別層河道分布

采用分頻融合技術對WVD-MEM計算的低、中、高三個頻段的功率譜特征開展多頻融合增強顯示，進一步增強河道識別效果。圖11a為RGB融合后識別效果，圖11b為CMY融合后的河道識別效果。由圖可見，經過多頻融合后，不同尺度、不同流向的河道得到進一步清晰展示，對一些較寬河道刻畫效果更加明顯(紅色箭頭所示)，對于原始振幅切片上無法識別的兩條不明顯窄河道(綠色箭頭所示)刻畫效果明顯改善。因此，利用WVD-MEM方法提取低、中、高三個頻段的功率譜屬性并融合顯示，能夠突出共性、弱化差異，有效改善單頻段頻譜在河道識別方面的不足，增強對工區河道的高分辨識別效果。

圖11 研究區WVD-MEM不同頻率功率譜融合層位切片特征

這里需要進一步指出的是，基于時頻分析方法進行河道識別，是將地震數據從時間域轉換到時頻域，即將地震數據從一個時間采樣點對應一個原始振幅值轉換到一個時間采樣點對應不同頻率下的多個值(無論振幅譜值還是功率譜值)，本質上不同于從時間域的角度分析地震信號，屬于在時頻域對原始地震信號的直接應用，不會破壞原始振幅的連續性，因而識別結果具有可靠性。而影響河道高分辨識別效果的因素取決于時頻分析方法的能譜聚焦性和分辨率。由于WVD-MEM方法結合了Wigner-Ville分布和最大熵原理的非平穩信號分析優勢，不僅能夠有效消除WVD信號分析中存在的交叉項，而且可以保持與瞬時頻率基本一致的高度時頻聚焦性，使計算的瞬時振幅、功率譜等地震屬性具有高分辨率特性。

該方法在川西地區侏羅系沙溪廟組氣藏河道識別的應用實例表明，WVD-MEM提取的功率譜、瞬時振幅等地震屬性，能從垂向和橫向較好地刻畫河道的沉積厚度、延展寬度和流動方向等地質細節，可以在其他類似的油氣探區開展推廣應用。

4 結論與認識

通過WVD-MEM方法研究和對非線性調頻信號、仿真地震信號、單道實測地震(非平穩)信號、楔狀模型正演模擬信號等定量試驗，以及在川西侏羅系沙溪廟組氣藏河道識別的實際應用，得到以下認識：

(1)相比于STFT、CWT、GST、WVD、SPWVD、CWD等時頻分析方法，WVD-MEM時頻分析方法結合了WVD的高度時頻聚焦性以及MEM的最大概率分布優勢，不僅完整消除了WVD中存在的嚴重交叉項，且具有顯著的頻譜聚焦性和高分辨率優勢，可以應用于薄儲層、窄河道及其他地下微型地質體的識別研究；

(2)針對埋藏較深、寬度較窄、厚度較薄、流向變化的河道高分辨識別難題，WVD-MEM方法能夠在一定程度上解決地震資料主頻不夠高、頻帶不夠寬等不利因素導致的空間分辨率不足的問題，利用WVD-MEM獲得的瞬時振幅、分頻功率譜等高分辨地震屬性以及多頻融合等手段，能夠有效突出河道的厚度、寬度、流向等空間展布細節，進一步提高河道的識別精度；

(3)在對WVD-MEM實測地震信號分析過程中，就整個工區而言，階數p一般選擇定值。如何在每道實測地震信號計算前自適應調整p的取值、優化WVD-MEM計算效率和程序設計以及方法的應用推廣，將是下一步研究的方向。

感謝中國石化西南油氣分公司為此研究提供實驗數據。