地震信號線性與非線性時頻分析方法對比

黃昱丞 鄭曉東 欒 奕 楊廷強

(①中國石油勘探開發(fā)研究院油氣地球物理研究所，北京 100083； ②香港中文大學理學院，香港 999077)

1 引言

時頻分析是處理地震信號等非平穩(wěn)信號的常用工具，現(xiàn)已廣泛應(yīng)用于沉積相帶、巖性、儲層和流體的檢測分析。其基本思想是將一維時域地震信號表達為時頻聯(lián)合函數(shù)的形式，在時頻二維空間進行描述，與傳統(tǒng)的Fourier分析相比，在分析信號局部統(tǒng)計特征方面具有明顯的優(yōu)勢。

時頻分析方法在地球物理勘探領(lǐng)域的應(yīng)用可上溯至二十世紀八、九十年代。1982年，Morlet等[1]提出用Gauss包絡(luò)不同波長的零相位余弦子波描述層狀介質(zhì)下地震波傳播與散射機制，通過構(gòu)造復(fù)函數(shù)小波，在時頻域開展地震信號變化特征的研究，其瞬時譜的思想直接觸發(fā)了小波變換(WT)的萌芽。Okaya等[2]利用短時Fourier變換(STFT)對炮集數(shù)據(jù)的地表振動耦合引起的掃頻假象進行時頻濾波。Chakraborty等[3]則將連續(xù)小波變換(CWT)和匹配追蹤分解(MPD)引入地震信號處理，獲得了比STFT更高精度的時頻譜，引起業(yè)界廣泛關(guān)注。Grubb等[4]將離散小波變換用于地震勘探信號的分解與重構(gòu)。李宏兵等[5]基于正演模型將WT用于含氣砂巖的地震衰減檢測。Partyka等[6]將時頻分析正式引入地震資料解釋領(lǐng)域，利用分頻剖面成功地實現(xiàn)了地層厚度估算和可視化，形成了定性與半定量解釋技術(shù)——譜分解(Spectral Decomposition)。Castagna等[7]提出基于地震瞬時譜直接檢測烴類的技術(shù)。Liu等[8-10]則提出了基于Morlet子波和Ricker子波的MPD算法與譜均衡方法，拓寬頻帶，提高了數(shù)據(jù)分辨率，并將其用于河道識別與儲層預(yù)測。Wang[11，12]將MPD應(yīng)用于地震資料精細解釋與烴類檢測，并發(fā)展了多道MPD算法。Sinha等[13]提出了對CWT逆變換做Fourier變換以直接獲得時頻譜的時頻連續(xù)小波變換，省去了CWT中尺度圖與時頻譜之間的繁瑣轉(zhuǎn)化。Pinnegar等[14]、高靜懷等[15]、陳學華等[16，17]分別在S變換(ST)基礎(chǔ)上提出了窗函數(shù)變化可控的參數(shù)化ST方法，統(tǒng)稱為廣義S變換(GST)，并將各自方法應(yīng)用于地震相分析、沉積旋回的識別和油氣檢測。徐陽等[18]利用GST和離散WT進行面波干擾壓制，同時可以減少對有效波的損傷。Li等[19]將平滑的Wigner-Ville分布(SWVD)法用于碳酸鹽巖儲層的地震資料衰減分析，Wu等[20]采用重排的平滑偽Wigner-Ville分布(RSPWVD)方法獲得更為精確的時頻譜用于低頻陰影檢測與儲層展布刻畫。Han等[21]利用第三代Hilbert-Huang變換(HHT)對地震信號進行譜分解，顯示了自適應(yīng)數(shù)據(jù)驅(qū)動類方法在地震數(shù)據(jù)分析中的有效性。曹思遠等[22]利用改進的HHT結(jié)合小波包變換，消除了虛假信號分量，反映了信號的真實頻率特性。薛雅娟等[23]基于集合經(jīng)驗?zāi)B(tài)分解方法(EEMD)結(jié)合WT進行本征模態(tài)函數(shù)的優(yōu)選，研究了地震信號的衰減特征并提出了一套優(yōu)化處理流程。劉晗等[24]應(yīng)用同步擠壓S變換進行地震信號時頻分析，分頻解釋效果優(yōu)于傳統(tǒng)方法。Khonde等[25]較為系統(tǒng)地總結(jié)了譜分解方法在地球物理勘探領(lǐng)域的發(fā)展和應(yīng)用現(xiàn)狀，包括在三維地震數(shù)據(jù)解釋[26]、直接烴類檢測[27]、水合物勘探[28]、煤系地層巖漿巖侵入機理[29]等方面的應(yīng)用。

總體而言，時頻分析方法根據(jù)核函數(shù)表達形式的不同可分為線性、雙線性和非線性三大類。線性方法包括Gabor展開[30]、STFT[31]、CWT[32]、ST[33]和GST[34-36]等。線性方法計算速度快，受限于不確定準則，時頻分辨率有限，李振春等[37]曾對線性方法做過系統(tǒng)總結(jié)。雙線性方法，如Wigner-Ville分布(WVD)[38]、偽Wigner-Ville分布(PWVD)、平滑偽Wigner-Ville分布(SPWVD)[39]等，歸于Cohen類時頻分布范疇。其中，WVD具有所有時頻分析方法中最高的時頻聚集性能，但受交叉項干擾，分辨率極差，后續(xù)發(fā)展的加窗改進方法則以犧牲部分分辨率為代價壓制交叉項。除以上兩類方法之外的都歸為非線性方法，包括基于數(shù)據(jù)驅(qū)動的自適應(yīng)分解方法，如自適應(yīng)STFT[40，41]，根據(jù)信號特性靈活變化時窗長度，計算效率很高，分辨率較STFT有很大提高。Huang等[42]提出的經(jīng)驗?zāi)B(tài)分解(EMD)及其改進版本的集合經(jīng)驗?zāi)B(tài)分解(EEMD)、完全集合經(jīng)驗?zāi)B(tài)分解(CEEMD)[43，44]的核心均在于經(jīng)驗?zāi)B(tài)分解，盡管時頻分辨率非常高，并采用加噪提純技術(shù)消除模態(tài)混疊現(xiàn)象，但仍然存在邊界效應(yīng)等問題，且抗噪性較差，計算效率也不高。近年提出的融入EMD和CWT的同步擠壓小波變換[45]，并以此發(fā)展出的一大類基于不同線性方法(STFT、ST)的“后處理”方法——同步擠壓變換(SST)[46，47]，兼具較高的分辨率和計算效率并可實現(xiàn)信號重構(gòu)，然而對低信噪比信號分析魯棒性較差?；谛盘栂∈璞硎镜姆蔷€性方法，如包括基于Gabor原子的匹配追蹤分解[48，49]、基于Chirplet變換(CT)[50，51]的自適應(yīng)分解(Chirplet Decomposition，CD)[52]、指數(shù)調(diào)頻小波變換[53]等，這類方法若匹配適當可以實現(xiàn)信號的稀疏分解，分辨率很高且免于交叉項干擾，但結(jié)果常常受控于時頻原子庫的架構(gòu)，易出現(xiàn)原子彌散問題，參數(shù)較多，不易控制，且計算效率極低。還有衍生自CT的多項式調(diào)頻小波變換、樣條調(diào)頻小波變換、泛諧波調(diào)頻小波變換[54，55]、廣義線性調(diào)頻小波變換(GLCT)[56]等，這些新方法或通過擠壓時頻譜帶寬，或通過調(diào)頻核函數(shù)與時頻脊線迭代擬合的方式，往往能夠獲得較高精度的信號時頻表征，但一定程度上也受窗函數(shù)影響。近十年來還發(fā)展了如約束最小二乘譜分析[57]、基于反褶積的STFT[58]、自適應(yīng)稀疏時頻分解[59]、最大熵Wigner-Ville分布聯(lián)合算法[60]等方法。

時頻分析是一種將一維時域信號映射到二維時頻空間描述和表達的方法，在信號局部或瞬態(tài)特征刻畫方面具有獨特的優(yōu)勢。目前使用的時頻分析方法眾多，使用條件和適用范圍存在差異。Cas-tagna等[61]認為在地震解釋領(lǐng)域時頻分析方法沒有對錯之分，而是針對具體應(yīng)用合不合適的問題。因此，比較它們之間的異同，分析參數(shù)選擇的影響，對合理應(yīng)用時頻分析方法有實際價值。因此，系統(tǒng)地比較不同時頻分析方法有助于改善地震時頻分析在儲層預(yù)測中的應(yīng)用效果。

2 基本方法分析與參數(shù)選擇

本文選取了包括線性方法中的STFT、CWT、ST、GST，雙線性方法中的WVD、PWVD和SPWVD，以及非線性方法中的HHT、基于STFT的SST、CD和GLCT，總計三類11種方法，梳理各種方法的適用條件，研究參數(shù)選取的原則。上述各種時頻分析方法的基本公式與控制參數(shù)如表1所示。

2.1 線性方法

線性方法是指信號變換滿足線性疊加原理的時頻分析方法[31]，包括STFT、WT、ST、GST等。

STFT[62]的本質(zhì)是對信號加窗分段進行滑動譜分析，并假定信號在時窗內(nèi)平穩(wěn)。其時窗函數(shù)常用Gauss窗、Hamming窗等。采用Gauss窗可以獲得信號最小的時寬和帶寬乘積，此時STFT亦稱作Gabor變換[30]。STFT的概念與定義所反映的物理意義簡單明確，計算效率很高，但受固定窗函數(shù)影響導(dǎo)致分辨率有限，時窗長度是唯一制約時頻表征的參數(shù)。實際應(yīng)用中，根據(jù)輸入信號特點和不同的分析需求權(quán)衡時間分辨率和頻率分辨率。

表1 時頻分析方法歸類及定義匯總

WT在STFT基礎(chǔ)上對時頻網(wǎng)格進行改造，具有多分辨率特性[63]。由于地震信號分析中常用Gabor小波或Morlet小波，這類母小波不存在尺度函數(shù)，故一般以存在信息冗余的連續(xù)小波變換(CWT)的形式出現(xiàn)。在CWT表達式中，大尺度對應(yīng)低頻端，時間分辨率低，頻率分辨率高；反之，小尺度對應(yīng)高頻端，時間分辨率高，頻率分辨率低。為分析的直觀性，一般需要將尺度轉(zhuǎn)換為頻率，它反映的是一個頻段而非頻點。良好的時頻局域特性是CWT的優(yōu)勢所在，但若母小波選擇不當，包括類型和峰值頻率等因素，應(yīng)用效果會大受影響。

ST繼承和發(fā)展了STFT和CWT局部化的優(yōu)點[33]，信號的ST是一種非嚴格意義上的CWT，若將ST看做采用特定窗函數(shù)的STFT，則相比Gauss窗函數(shù)，ST的窗函數(shù)相當于將尺度參數(shù)σ(標準差)表達為頻率絕對值的倒數(shù)σ=1/|f|。與CWT類似，既滿足了多分辨特性，同時又不需要從尺度到頻率的轉(zhuǎn)換，物理意義更加直觀明了，計算效率很高，同時存在逆變換，可以實現(xiàn)信號濾波與重構(gòu)。

雖然ST是小波變換的改進，但對信號局部時頻特性的分析能力仍有待提升，因為ST采用固定的母小波，其時間窗形態(tài)隨頻率線性變化，使該方法在實際應(yīng)用中仍受到一定限制。因此，眾多學者提出了若干形態(tài)各異的時窗作為ST的基函數(shù)，統(tǒng)稱為GST。本文采用陳學華[64]等提出的GST，即引入兩個參數(shù)λ和p對ST的核函數(shù)進行改造，在ST的基礎(chǔ)上，改寫尺度參數(shù)σ=1/(λ|f|p)，其中λ>0、p>0。參數(shù)λ和p調(diào)節(jié)窗函數(shù)形狀和變化趨勢。λ和p越大，GST高頻端時間分辨率越高，頻率分辨率越低，且變化趨勢對指數(shù)參數(shù)p更敏感，但調(diào)節(jié)過程中應(yīng)避免高頻失真的情況。因此，GST通過改變窗函數(shù)隨頻率變化的形態(tài)進一步改善分辨率，在實際應(yīng)用中更為靈活。

總之，線性方法均受Heisenberg不確定性原理限制，無法同時在時域和頻域獲得最高的能量聚集性能。

2.2 雙線性方法

雙線性方法屬于非線性方法范疇，滿足二次疊加原理[31]。這類方法從能量角度刻畫信號的時頻分布，相比STFT的譜圖或CWT的尺度圖等線性方法取模值平方的方式，雙線性方法更適于描述地震信號能量隨時間頻率的變化。為了避免負頻率帶來的交叉項干擾，雙線性方法一般是把基于實信號Hilbert變換獲得的解析信號作為輸入項以獲取單邊時頻譜。

所有的雙線性方法中，WVD[38]因具有良好的邊緣特性和極高的時頻聚集性能而成為最基本、應(yīng)用也最廣泛的時頻分布。解析信號z(t)的WVD定義為信號瞬時相關(guān)函數(shù)關(guān)于時延參數(shù)τ的Fourier變換。WVD對單分量平穩(wěn)信號或者線性調(diào)頻信號具有最高的時頻分辨率，但其嚴重缺陷在于分析多分量信號時面臨交叉項的干擾，也正由此派生了眾多圍繞抑制交叉項的雙線性方法。

WVD可看做是窗函數(shù)為脈沖函數(shù)δ(t)的雙線性分布，為了抑制交叉項，對時延參數(shù)τ加窗，實現(xiàn)對頻率軸的平滑，由此衍生出PWVD[39]。本質(zhì)上PWVD是以犧牲一定的時頻聚集性能為代價，對WVD的結(jié)果進行了低通濾波。時窗越長，交叉項抑制越明顯，但(頻率)分辨率也越低。若τ時窗的窗寬為1，則退化為WVD。

事實上，若考慮同時對時間變量t和時延參數(shù)τ加窗平滑，亦即同時在時間和頻率域進行平滑，可進一步抑制交叉項，從而出現(xiàn)了SPWVD[39]，但從一重積分變?yōu)槎胤e分的代價是計算量明顯增加。具體而言，時間變量窗口越寬，每次循環(huán)參與計算的樣點數(shù)越多，耗時越長。若時間變量窗口寬度為1，則退化為PWVD。

2.3 非線性方法

非線性方法包括：基于數(shù)據(jù)驅(qū)動的HHT、時頻原子匹配追蹤自適應(yīng)分解以及同步擠壓變換等參數(shù)化方法，它們都屬于廣義非線性方法范疇。

HHT是Huang等[42]提出的基于數(shù)據(jù)驅(qū)動的時頻分析方法，其算法實現(xiàn)分為兩步：①經(jīng)驗?zāi)B(tài)分解(EMD)，將數(shù)據(jù)信號分解為有限個本征模態(tài)函數(shù)(IMF)和背景趨勢值或殘差；②Hilbert譜分析，對這些IMF進行Hilbert變換求取瞬時頻率和瞬時振幅，對瞬時頻率重排，將信號映射為二維時頻譜，亦稱Hilbert譜。HHT可以獲得高精度線譜，但邊界效應(yīng)、模態(tài)混疊和薄弱的數(shù)理基礎(chǔ)成為HHT推廣應(yīng)用的困難。為避免模態(tài)混疊問題，本文采用集合經(jīng)驗?zāi)B(tài)分解(EEMD)[43]算法篩選IMF，該算法參數(shù)主要包括參與計算的白噪聲集合數(shù)目NE和白噪聲標準差Nstd。NE越大，理論上IMF越穩(wěn)定，但耗時越長。對強時變信號，取較大的Nstd為宜，但一般不超過30%。

SST源自于Daubechies等[45]提出的一種“后處理”時頻分析方法——同步擠壓小波變換(Synchro-squeezed Wavelet Transform，SWT)。SWT的基本思想借鑒了HHT的EMD方法，將信號看做多個本征函數(shù)的疊加，在CWT基礎(chǔ)上，進一步“擠壓”瞬時頻帶，獲得更精確的時頻譜?；谛〔ㄏ禂?shù)對時間求導(dǎo)，初步估計瞬時頻率。這樣，將時間—尺度域轉(zhuǎn)化為時間—頻率域。沿用這一思路，將擠壓操作與STFT、ST等方法結(jié)合，相繼出現(xiàn)了同步擠壓短時Fourier變換[65]、同步擠壓S變換[46]等方法。將這類方法統(tǒng)稱為“同步擠壓變換”(Synchro-squeezing Transform，SST)。SST具有極高的時頻聚集性，運算效率較高，但基于微分算子的擠壓操作會引起數(shù)值不穩(wěn)定，抗噪性低。本文采用基于STFT的SST進行信號分析，其優(yōu)勢在于時頻網(wǎng)格均勻性和中低頻段的信號時頻特征刻畫的精度和穩(wěn)定性[66]。

CD源于Mann等[50]提出的Chirplet變換(CT)，在Gabor原子基礎(chǔ)上加入線性調(diào)頻相位項，實現(xiàn)在時頻空間對非平穩(wěn)信號的一階逼近?；贛PD[48]的思想，利用預(yù)定義的歸一化四參數(shù)Gaussian Chirplet原子庫[67]，對實信號進行Hilbert變換得到解析信號，通過不斷迭代匹配搜索最佳參數(shù)集合，滿足時頻原子在信號殘差方向上投影最大化。經(jīng)過若干次迭代滿足判定條件，匹配終止并構(gòu)建無交叉項WVD譜。CD在時頻原子庫選擇恰當?shù)那闆r下，可以實現(xiàn)信號的稀疏匹配，時頻聚集性能極高，并免于交叉項干擾。本文使用的算法[68]需要預(yù)定義匹配分解的Chirplet原子個數(shù)NA，原子分解個數(shù)M控制CD分解擬合的精細程度，算法耗時正比于NA，故分量越多、變化越復(fù)雜的信號，處理耗時越長。

GLCT由Yu等[56]提出，其本質(zhì)是在線性調(diào)頻小波變換(Linear Chirplet Transform，LCT)的基礎(chǔ)上進行調(diào)頻斜率旋轉(zhuǎn)優(yōu)化。通過在每個時頻樣點上計算不同的調(diào)頻斜率值對應(yīng)的LCT，挑選其投影能量最大值作為該點的時頻響應(yīng)。實際操作中，需要對時頻空間離散化，Chirplet原子旋轉(zhuǎn)角度劃分的稠密度N是關(guān)鍵參數(shù)。理論上N值越大，擬合局部時頻特征越準確。GLCT在處理強時變、多分量信號方面具有一定優(yōu)勢，與傳統(tǒng)方法相比，GLCT可以自適應(yīng)擬合多分量非平穩(wěn)信號時頻特征，具有較高的時頻分辨率，計算效率優(yōu)于稀疏匹配類算法，且表征多分量信號時頻譜時不受交叉項干擾。不足之處在于其本身受窗函數(shù)影響(默認采用Gauss窗)，且N的選取存在主觀性，隨著N值增大，計算量也顯著增加，一般不宜超過5。

3 時頻分辨率分析

分別采用上述方法對單分量諧波調(diào)頻信號、多分量指數(shù)調(diào)頻信號以及Ricker子波合成地震信號進行分析，前兩類信號給出基于瞬時相位變化率定義的瞬時頻率理論值作為參考，信號長度均為1024樣點，采樣間隔1ms。關(guān)于時頻分辨率的優(yōu)劣一般都是來自時頻譜視覺效果上的直接比對，定量分析指標很少，但基于信息熵的時頻聚集性檢測則可以為此提供一個參考的定量指標。

Baraniuk等[69]提出利用Rényi熵估計信號所含信息量與時頻譜復(fù)雜程度的思路。信號時頻譜TFRs(t,f)的Rényi熵定義為

(1)

式中：α為Rényi熵的階次，一般取3； 對數(shù)底b一般取2即可。Rényi熵是Shannon熵的廣義形式，Shannon熵是Rényi熵在α→1時的極限，Rényi熵更具普適性。Rényi熵值越小，說明信號在時頻域的能量復(fù)雜度越低，時頻聚集性越高；反之，則復(fù)雜度越高，能量分布越均勻，聚集性越差。需要注意的是，聚集性和分辨率是兩個概念，采用三階Rényi熵僅作為刻畫時頻分析方法能量聚集性高低的一種定量手段，聚集性高并不說明分辨率一定高，也不一定反映信號真實的局部頻率特征。

3.1 模型1：諧波調(diào)頻信號

諧波調(diào)頻信號的瞬時頻率具有周期振蕩特征(圖1)，時頻譜上呈現(xiàn)諧波形態(tài)，是一類典型的非平穩(wěn)信號，即

s(t)=cos{2π[sin(30t)+60t]t}

(2)

信號的瞬時頻率為

IF(t)=30cos(30t)+60

(3)

利用該信號考察各類方法對振蕩調(diào)頻瞬時頻率曲線的時頻聚集性能。由于瞬時頻率振蕩變化較為劇烈，故采用較短時窗處理。由圖1可見，線性方法的分辨率普遍較低，WVD因交叉項的干擾出現(xiàn)頻率假象。STFT、PWVD、SPWVD和GLCT等方法可以較準確地反映信號瞬時頻率變化特征，但聚集性略差。Chirplet原子分解呈現(xiàn)時頻特征的一階逼近，故在擬合信號瞬時頻率曲線峰谷點處誤差較大。無噪聲情況下，對于這種固定周期震蕩調(diào)頻信號的時頻表征，HHT和SST具有最高的時頻聚集性能。

圖1 諧波調(diào)頻信號、真實瞬時頻率以及不同方法獲得的時頻譜

方法STFTCWTSTGSTWVDPWVDSPWVDHHTSSTCDGLCT熵4.98855.10514.89125.07694.72914.52154.62670.04621.71213.22114.7242

表2表明HHT和SST兩種方法的時頻聚集性是最高的，同時，其余雙線性、非線性方法的時頻聚集性能也都高于線性方法，這與時頻譜時頻分辨率一致。但WVD時頻譜明顯受交叉項干擾，Rényi熵卻低于所有線性方法，也說明Rényi熵測度的定量判據(jù)肯定了WVD的高聚集性，但忽視了其交叉項干擾導(dǎo)致的低分辨率事實。

3.2 模型2：指數(shù)調(diào)頻信號

參考Sinha等[13]的模型，信號由兩個具有不同指數(shù)調(diào)頻因子的掃頻信號疊加而成，即

(4)

信號的兩個瞬時頻率分量為

(5)

兩分量在0.6s前相互緊貼，隨后瞬時頻率變化趨勢逐漸分開，信號后半段頻率變化劇烈。

利用該信號有效頻帶范圍遠離Nyquist頻率，考察各類方法的多分辨性能，并檢驗低頻端時域分辨率，對各類方法性能提出了較高要求(圖2)。由圖2可見線性方法中只有GST可以較好地區(qū)分開兩個分量并兼顧高低頻端的分辨率。雙線性方法受交叉項干擾嚴重，或者也因窗函數(shù)導(dǎo)致高頻段分辨率降低。SST受窗函數(shù)影響，長時窗導(dǎo)致高頻端分辨率較低。Chirplet原子分解因為原子匹配度不高，一階擬合信號時頻分量存在較大誤差。HHT和GLCT兩種方法可以高分辨地顯示信號的兩個掃頻分量，但HHT在高頻段存在IMF線譜振蕩失穩(wěn)的現(xiàn)象，而GLCT存在Chirplet時頻原子旋轉(zhuǎn)交疊的問題?？傮w而言，GST和HHT對該信號時頻刻畫效果最佳。

圖2 指數(shù)調(diào)頻信號、真實瞬時頻率以及不同方法獲得的時頻譜

方法STFTCWTSTGSTWVDPWVDSPWVDHHTSSTCDGLCT熵3.87403.49374.01574.05104.79724.11414.02880.99062.26422.75364.3537

從Rényi熵角度看，HHT表現(xiàn)出了最高的時頻聚集性能，GLCT算法本身在時頻面每個樣點能量都是非零的，故而Rényi熵較高，也說明了GLCT處理此類指數(shù)調(diào)頻信號存在不足。

3.3 模型3：合成地震信號

參考Castagna等[7]的模型，給出不同峰值頻率和時延的零相位Ricker子波合成的地震信號。圖3給出了各單頻子波及合成信號，可知有十二個零相位Ricker子波參與信號合成，包括兩個40Hz子波、七個30Hz子波、兩個20Hz子波和一個10Hz子波。該信號考察各類方法的多分辨性能，檢驗高頻端頻域分辨率亦即薄層(相鄰子波)分辨率。

圖3 合成地震信號、單頻子波以及不同方法獲得的時頻譜

由圖3可見，線性方法中STFT不能兼顧高低頻，CWT、ST、GST等方法在低頻端和高頻端總體聚集性不高。相對而言，CWT、GST高頻端可以區(qū)分相鄰Ricker子波。SST呈現(xiàn)稀疏表示的特征，但對薄層的分辨率并不高。雙線性方法中SPWVD表現(xiàn)最好，高低頻分辨率大致保持，并具有較高的薄層分辨率，WVD、PWVD受交叉項干擾。HHT時頻聚集性能極高，但從圖中不能清晰區(qū)分子波各分量，分辨率不高。參數(shù)選取恰當時，CD可以提取出大部分子波分量，對低頻的分辨能力較為突出，但薄層仍不能分辨。相比之下，GLCT則具有類似SPWVD的優(yōu)異表現(xiàn)。

由于信號的稀疏性，HHT和SST的Rényi熵值均小于零，故這兩種方法的時頻聚集性最高，然而，兩者對合成地震信號時頻譜的刻畫均出現(xiàn)較大程度的失真。

表4 合成地震信號不同方法時頻譜的Rényi熵

4 計算效率分析

在XW9400型HP工作站(64位、AMD Opte-ron 2.8GHz雙核處理器、RAM 8GB)上處理以上三個模型信號，計算時間如圖4所示。所有的線性方法以及雙線性方法中的WVD和PWVD處理1024個樣點信號的耗時都在10-2～10-1s數(shù)量級，其余大部分非線性方法耗時都在1s以上，高出一到兩個數(shù)量級。大部分方法計算效率受窗函數(shù)長度影響顯著，如SPWVD在處理模型2時采用較寬的雙重時窗上處理時間較多。此外，SST耗時主要體現(xiàn)在頻帶擠壓操作中，但在非線性方法中效率相對較高；HHT在加噪提純與IMF的篩選中耗時過長；CD處理時間正比于匹配的原子個數(shù)，若一階擬合三個模型信號，分別預(yù)定義了10個、8個和12個原子，這種基于稀疏匹配的貪婪算法成為最耗時方法；GLCT耗時正比于原子旋轉(zhuǎn)角度劃分的稠密度N?？傮w而言，線性方法的計算效率遠高于大部分非線性方法，非線性方法中諸如CD等方法是基于MPD算法思想，雖然對某些信號的時頻估計效果較好，但在實際應(yīng)用中計算成本較高。

圖4 不同時頻分析方法三個模型信號處理效率對比

一般而言，地震信號記錄長度很有限，所以在單道試算時，時間頻度表達式中被忽略的不同系數(shù)最高次冪和各低次冪項造成的時間消耗不容忽視?？傮w而言，各類方法的計算效率優(yōu)劣對比明顯(一到兩個數(shù)量級)，實際計算耗時會因為海量的地震道數(shù)(如涉及疊前或三維地震數(shù)據(jù)的計算)而陡增，影響算法的實際應(yīng)用，這也是應(yīng)當考慮的現(xiàn)實問題。

5 抗噪性能分析

參考Yu等[56]的抗噪性能分析，對不同方法獲得的時頻譜|TFR(t,f)|進行瞬時峰值頻率檢測。對每一個樣點t，瞬時峰值頻率為

(6)

式中fc為信號的Nyquist頻率。通過與理論瞬時頻率曲線IF(t)相比較可計算IFe(t)的均方誤差

(7)

以模型1諧波調(diào)頻信號為例，對信號加入不同級別(無噪聲，10dB，0dB)高斯白噪聲，通過對比各類方法估計的瞬時峰值頻率的均方誤差考察其穩(wěn)健性。WVD因交叉項干擾嚴重，而CD本身的不穩(wěn)定性較高，兩種方法檢測結(jié)果均明顯偏離真實值，故未參與測評。

不同的方法隨著噪聲級別增大(圖5a)，瞬時頻率估計值的標準偏差都不同程度地增大，其中CWT、SPWVD和GLCT的穩(wěn)健性最高，是低信噪比情況下抗噪效果較好的方法；由不同方法的抗噪性對比(圖5b)可見STFT、ST、GST、SST、PWVD和HHT隨噪聲級別增大估計偏差迅速增加，尤以GST(9.33%)、ST(7.61%)、HHT(6.82%)和SST(6.56%)受噪聲干擾嚴重。對高信噪比(無噪聲、加10dB白噪)數(shù)據(jù)，線性方法中的STFT(0.08%)和非線性方法中的GLCT(0.05%～0.09%)具有最高的穩(wěn)健性；對低信噪比(0dB)數(shù)據(jù)，CWT(0.44%)、SPWVD(0.57%)和GLCT(1.84%)是首選方法。

綜上所述，各類方法性能表現(xiàn)按“最低<低<一般<較高<高”排列分級，總結(jié)如表5所示。

圖5 不同噪聲級別下不同方法估計偏差對比(a)和不同方法在不同噪聲級別下的穩(wěn)健性變化(b)

方法窗函數(shù)交叉項聚集性低頻分辨率薄層分辨率計算效率抗噪性STFT固定無低低低高較高CWT多分辨無低較高較高高高ST多分辨無低較高一般高低GST多分辨無低高高高低WVD無有高低低高-PWVD固定有較高高一般高一般SPWVD固定有較高高高一般高HHT無無高較高一般低低SST固定無高高較高較高低CD無無較高高低最低-GLCT固定無較高較高高一般高

6 應(yīng)用實例分析

本文實際數(shù)據(jù)所在研究區(qū)位于塔里木盆地塔中地區(qū)，目的層段是上奧陶統(tǒng)良里塔格組碳酸鹽巖，埋深約5000m，溶蝕孔隙極為發(fā)育，非均質(zhì)性強[19]。工區(qū)三維地震數(shù)據(jù)采樣間隔4ms，截取3000～4500ms時窗，層位mfs與sb3為目的層段頂、底界面，資料信噪比較低，主頻不足20Hz，常規(guī)剖面(圖6)難以識別儲層形態(tài)與含油氣特征。

選取STFT、CWT、SPWVD和GLCT作為代表，進行效果對比。圖7給出了四種方法5Hz分頻剖面，STFT剖面中反射軸不連續(xù)，也不能區(qū)分目的層頂、底；CWT低頻端時窗較寬，因而時間分辨率很低，難以分辨目的層響應(yīng)；SPWVD是基于能量的雙線性分布，其淺層強反射掩蓋了目的層信息，頂、底區(qū)分不開；GLCT顯示了沉積層較多細節(jié)，目的層臺地內(nèi)頂?shù)酌娣瓷浔粎^(qū)分開并可以連續(xù)追蹤，淺層地層反射連續(xù)性也更好。圖8為四種方法獲得的35Hz層內(nèi)均方根振幅切片，在GLCT和SPWVD切片上臺地邊界(黃色箭頭所示)以及潮汐水道(虛線圓框所示)更為清晰，臺地內(nèi)部水道(白色箭頭所示)和瀉湖沉積區(qū)前緣與內(nèi)部斷裂體系(黑色箭頭所示)細節(jié)也更豐富，而STFT、CWT的切片效果相當，不能明顯觀察到這些地質(zhì)特征。總體上，盡管信噪比較低，兩類較穩(wěn)健的非線性時頻分析方法仍可以挖掘數(shù)據(jù)體中豐富的低頻信息和平面細節(jié)，儲層幾何特征與空間展布也得到精細刻畫。

圖6 實際數(shù)據(jù)剖面

圖7 實際數(shù)據(jù)5Hz不同方法分頻剖面

圖8 層間均方根振幅35Hz不同方法分頻切片

7 結(jié)論

本文系統(tǒng)分析了用于地震勘探領(lǐng)域的11種時頻分析方法，基于三個模型信號和碳酸鹽巖儲層實際地震數(shù)據(jù)比較了各類方法的時頻分辨率、計算效率和抗噪性能。結(jié)論如下：

(1)線性類方法受限于不確定準則，時頻聚集性能普遍較低，但計算效率很高，且不存在交叉項干擾；另外，STFT、CWT的穩(wěn)健性較高，GST、GLCT的薄層分辨率較高； 非線性類方法時頻聚集性能普遍較高，除SPWVD和GLCT外，抗噪性普遍較差；雙線性方法中WVD、PWVD計算速度快但受交叉項干擾較嚴重，其余非線性類方法計算效率普遍不高。

(2)除了信號本身特征外，窗函數(shù)普遍影響各種方法的時頻分辨率，而特定方法、特定的參數(shù)集合也會影響時頻分辨率、計算效率和抗噪性能，如HHT(EEMD)參與計算的白噪信號的數(shù)目、Chirplet原子分解法中預(yù)定義的Chirplet原子數(shù)目、GLCT中時頻面角度劃分稠密度等參數(shù)對時頻分辨率和計算耗時影響較大，前兩種方法的實際應(yīng)用中計算成本較高。

(3)碳酸鹽巖儲層實際地震資料分析顯示，穩(wěn)健性較高的非線性方法獲得的分頻剖面上可以更方便地追蹤儲層頂、底反射界面，分頻切片上可以更清晰地反映潮汐水道和礁灘相帶展布等特征。

(4)總體而言，低信噪比數(shù)據(jù)分析宜采用線性類方法，高信噪比數(shù)據(jù)宜采用時頻分辨率更高的非線性類方法，有可能獲得高精度的地層阻抗與局部頻率變化信息。