基于Wigner-Ville分布與Chrip-Z變換的高分辨時頻分析方法

李思源 徐天吉

(電子科技大學資源與環境學院，四川成都 611731)

0 引言

隨著勘探開發程度的不斷深化，具有圈閉規模小、儲層薄、烴源豐富卻分布零散等特點的深層微型油氣藏已成為研究熱點[1-2]。然而，受制于油氣藏埋深大、地震信號的分辨能力不足等不利因素，加大了對微型河道、薄砂體等油氣目標的識別難度。時頻分析技術能提供時間域和頻率域的聯合分布信息，故被廣泛應用于地震波等復雜非平穩信號的處理[3]，進而在油氣目標識別中發揮重要作用。

時頻分析技術現已成為研究沉積微相、薄儲層等地質目標的重要手段[4]。馬見青等[5]根據S變換提供的有效地震信息，成功識別了復雜沉積環境； 武迪等[6]結合變分模態分解與包絡導數算子預測碳酸鹽巖溶洞儲層的油氣分布； 朱秋影等[7]利用時頻分析技術預測依拉克構造有利河道砂體的空間展布； 蔣煉等[8]利用分頻技術預測砂體儲層厚度和分布范圍； 張付璦等[9]通過改進窗參數優化S變換從而精細刻畫沉積微相及相應構造； 楊道慶等[10]利用時頻分析及地層切片技術精確描述地下分流河道的空間流向、水平延展、厚度變化等特征。

時頻分析最常用的方法主要包括兩種：一種是以短時傅里葉變換(STFT)[11]、小波變換[12-14]、Gabor變換[15]、廣義S變換[16]為代表的線性時頻分布； 另一種則是以Wigner-Vile分布(WVD)、平滑偽Wigner-Ville分布(SPWVD)為代表的Cohen類廣義雙線性時頻分布[17]。線性時頻分析方法雖然計算速度較快，但存在時頻聚焦性欠佳且分辨率低等不足。WVD具有較高分辨率和較突出的能量集中性，但因其非線性特征，在處理多分量數據時必然會產生大量交叉項從而干擾有效信號，很大程度上限制了其應用范圍[18]。SPWVD屬于Cohen類二次型時頻分布的核函數方法，即基于WVD對信號瞬時自相關函數的時域和頻域方向同時施加窗函數，改變核函數以抑制交叉項的產生。這種核函數法出色地抑制了交叉干擾項的生成，且仍保留了WVD的眾多數學特性，適用于處理復雜地震數據[19]。

地震信號通常具有主頻低、頻帶窄等特點[20]，其時頻分布在有效頻帶內的劃分點數往往較少，難以保證理想的分辨率，使STFT、Gabor變換、廣義S變換、WVD和SPWVD等時頻分析技術在識別埋藏深、厚度小的微型地質目標時的精度不夠高，必然會對油氣藏地下空間分布的認識、儲量估算、鉆井勘探開發過程產生不利影響。

本文在SPWVD基礎上，以線性調頻Z變換(CZT)[21]替代SPWVD中的快速傅里葉變換(FFT)，提出一種提高地震信號時頻分辨率的計算方法(SPWVD-CZT)。利用CZT在三維空間內螺旋采樣的特性，通過調整采樣間隔和采樣點數實現數據插值計算； 在SPWVD基礎上進一步完善地震數據時頻譜的局部細節信息，為地下微型河道、薄儲層等小幅度復雜油氣目標的識別提供方法支撐。

1 SPWVD-CZT方法原理

1.1 SPWVD的定義

定義信號x(t)的Wigner-Ville分布[22]為

(1)

式中：Wx(t，σ)為信號x(t)的WVD，t為時間變量；σ為頻率變量；τ為時間延遲； e-jστ為傅里葉變換的參數因子；x*為x的共軛。定義

(2)

式中Rx(t，τ)為x(t)的瞬時自相關函數。

在時域和頻域同時對WVD進行加窗處理可抑制交叉項的產生，得到SPWVD分布[23]

(3)

式中：Spw(t，σ)為信號x(t)的SPWVD；h(τ)、g(u-τ)分別為時間和頻率方向的窗函數，其中u為頻率延遲。

總之，SPWVD能更清晰地反映信號能量的時頻分布，同時極大程度地抑制了交叉項的產生，被廣泛應用于地震數據的處理。

1.2 CZT的定義

CZT是對FFT改進得到的，其原理是在Z平面內用一條螺旋線進行等間隔采樣，采樣點在螺旋線上呈現等角分布[24]。該方法計算功率譜和時頻分布十分有效，且具有計算快速、不受數據序列長度限制、靈活性強等優點。

長度為N的離散信號序列x(n)的Z變換為

(4)

式中：n為離散時間變量；z為Z變換采樣因子。

令zr=AW-r，A=A0ejθ0，W=W0ejφ0，則有

zr=A0ejθ0W0e-jφ0r

(5)

式中：zr為Z變換的采樣因子，代替式(4)中的z,其中r為空間步長； ejθ0為步長因子； e-jφ0r為幅角因子；A0、W0均為任意正實數。給定A0、W0、θ0、φ0，當r=0，1，…，∞時，可得到Z平面上的一系列點z0、z1、…、z∞，再對其做Z變換

(6)

式中：A-n為采樣步長；Wnr為采樣幅角。

對信號做頻譜分析時，應在單位圓上實現CZT[25]。A0和W0都應取為1，設離散信號X(n)的長度為n=0，1，…，N-1，變換的長度r=0，1，…，M-1，式(6)就變為

(7)

由Bluestein公式

(8)

式(7)又可寫為

(9)

(10)

1.3 離散SPWVD-CZT的實現

將CZT與SPWVD結合，利用CZT的特性，通過在Z平面上對信號采樣重構增加頻帶內的劃分點數，從而改善信號的局部細節，有效提升信號的時頻分辨率。實現過程如下：

將連續信號離散化x(t)?x(n)，設時間采樣點數為N的解析信號為z(n)，有

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

(11)

式中：解析信號z(n)的實部是x(n)，虛部H[x(n)]為x(n)的Hilbert變換。z(n)的WVD定義為

(12)

式中：Wz(n，σ)為z(n)的WVD；l為步長；z*為z的共軛； e-j2lσ為傅里葉變換因子。

在頻域方向加上長度為N/4的矩形窗g(l)，當|l|>N/8時，g(l)=0。可得頻域方向加窗后的離散瞬時自相關函數為

R(n,l)=z(n+l)z*(n-l)g(l)g*(-l)

(13)

式中g*為g的共軛。

再對R(n，l)的時域方向加上長度為N/2的矩形窗h(l)，當|l|>N/4時，h(l)=0。可得時域方向加窗后的離散瞬時自相關函數為

z*(n-l)g(l)g*(-l)

(14)

此時，K(n，l)也被稱為SPWVD的核函數。

由于CZT和FFT的計算域為非負數(l=0，1，…，N-1)，還需對核函數K(n，l)做循環位移，使得l位于0～(N-1)內，得到重新排序后的核函數

(15)

式中：K(n，l)為排序前正頻率方向的核函數；K(n，-l)為排序前負頻率方向的核函數。

設待分析頻段的點數為M，根據FFT的快速計算原理，取最接近N+M-1的二次冪Q作為計算時的序列長度。由于WVD的周期為π，還應將CZT中的相角θ0和幅角φ0擴展為原來的兩倍，即A=A0ej2θ0，W=W0ej2φ0。構造序列

(16)

單位響應序列δ(l)的長度為Q，對SPWVD的序列補0，使得兩個卷積的序列長度相同。再構造序列

(17)

式中：f(l)為重排的序列；f*(l)=KN(n,l)A-l×Wl2/2，0≤l≤N-1；KN(n，l)為式(15)中SPWVD循環位移后得到的核函數。

然后分別對δ(l)、f(l)做傅里葉變換

(18)

再將Δ(r)與F(r)相乘，得到

Y(r)=F(r)Δ(r)

(19)

進一步作逆變換

y(l)=IFFT[Y(r)]

(20)

將y(l)乘上加權系數

(21)

式中SPWVD-CZT(zl)為信號的SPWVD-CZT分布結果。

2 方法測試

2.1 模擬信號時頻分析

利用雷克子波可合成與實際地震信號近似的模擬信號，以驗證SPWVD-CZT的時頻分析能力。采用主頻分別為40、20、10Hz的雷克子波，由下式

(22)

合成如圖1所示的模擬地震信號，可見其振幅峰值分別出現在0.1、0.3、0.5s位置。

圖1 模擬地震信號波形

將模擬地震信號做Hilbert變換，對加窗后的核函數做FFT獲得SPWVD結果。與實際地震信號一樣，模擬地震信號也具有主頻低、頻段窄的特點。選擇頻段0～100Hz，觀察信號完整的時頻分布。

對比模擬信號的SPWVD(圖2a)和SPWVD- CZT(圖2b)處理結果，可見模擬地震信號在0.1、0.3、0.5s處的能量分布分別集中于40、20、10Hz處； SPWVD-CZT保留了SPWVD完整的時頻分布和數學特性。由于模擬地震信號較平滑且時頻域范圍較大，難以直觀地看到信號分辨率的提升，因此選取0.2～0.4s時段、10～40Hz頻段做局部細節放大顯示(圖2c和圖2d)，對比可見SPWVD-CZT算法利用采樣插值特性增加細化點數方式平滑了選定頻段的時頻分布(圖2d)，進一步增強了模擬信號時頻聚焦能力，同時完善了局部邊緣細節，提升了信號的時頻分辨率。

圖2 模擬地震信號的時頻分析效果對比(a)SPWVD； (b)SPWVD-CZT； (c)圖a局部放大顯示； (d)圖b局部放大顯示

2.2 實際單道地震記錄時頻分析

相較于模擬地震信號，實際地震信號的頻率、相位、振幅等信息更復雜。通過對實際單道地震記錄進行SPWVD-CZT分析，進一步驗證該時頻分析方法的有效性。

圖3為實際單道地震信號，采樣率為2ms，時間采樣點數(N)為1001。對其分別于0、100Hz頻段做SPWVD-CZT處理(圖4b)，細化點數(M)為600。

圖3 實際單道地震信號波形

與實際單道地震信號的SPWVD處理結果(圖4a)對比，可見SPWVD-CZT算法(圖4b)對復雜實際單道地震信號仍具有完整的時頻分析能力，不僅準確刻畫出能量在時域和頻域內的分布，而且保留了SPWVD優良的數學特性。選取0.5～1.5s時段、50～70Hz頻段，同樣對上述處理結果做局部細節放大顯示(圖4c、圖4d)，對比可知SPWVD-CZT算法(圖4d)通過在各采樣點之間進行插值，增加有效頻段內的劃分點數，在保持高度時頻聚焦性的同時平滑了地震信號的能量分布，更精確地刻畫出功率譜隨時間和頻率的變化。

圖4 實際單道地震信號的時頻分析對比(a)SPWVD； (b)SPWVD-CZT； (c)圖a局部放大顯示； (d)圖b局部放大顯示

3 應用案例

川西中侏羅統沙溪廟組屬于淺水三角洲沉積體系，油氣資源豐富，是目前四川盆地天然氣產能建設的重要層系之一。依據沉積特征和地層厚度，沙溪廟組被劃分為上沙溪廟組和下沙溪廟組，其中上沙溪廟組又進一步細分為沙溪一段和沙溪二段，下沙溪廟組為沙溪三段。上沙溪廟組厚度為450～700m，是由棕(褐)色泥巖、粉砂質泥巖與褐灰、淺綠灰色細砂巖組成的厚互層； 下沙溪廟組厚度為134～250m，是由暗(褐)紫色含粉砂質泥巖為主、夾淺(綠)灰色細粒灰質巖屑長石砂巖組成的薄互層。

圖5 原始剖面(a)與20(b)、35(c)、50Hz(d)的SPWVD-CZT單頻剖面

圖三小層的兩種算法的50Hz功率譜層的層的層的SPWVD；層的層的層的SPWVD-CZT

4 結論

通過SPWVD-CZT方法研究，利用模擬地震信號、實際單道地震信號驗證該算法對信號時頻分辨率的提升效果，并將其應用于川西中江地區分流河道的識別。獲得以下認識和結論。

(1)SPWVD-CZT具有優良的時頻聚焦性，能完整地反映信號能量在時頻域內的變化，它利用采樣插值提升特性，增加有效頻段劃分點數，在SPWVD基礎上進一步提升了信號時頻分辨率。

(2)利用SPWVD-CZT提取的模擬地震信號和實際單道記錄的瞬時功率譜的能量分布更平滑、更均勻，時頻譜圖的邊緣細節也更清晰。

(3)針對埋藏深、隱蔽性強的分流河道儲層，SPWVD-CZT能在一定程度上克服地震信號主頻低、頻帶窄導致的分辨率不足問題，突出地下河道的橫向寬度、縱向厚度、沿層走向等空間分布特征，提升油氣資源勘探精度。該方法也為微型地質體識別提供了一種新思路，且拓展了傳統時頻分析理論。