黃昱丞 鄭曉東 欒 奕 楊廷強(qiáng)

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

1 引言

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

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

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

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

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

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

2.1 線性方法

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

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

表1 時(shí)頻分析方法歸類及定義匯總

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

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

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

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

2.2 雙線性方法

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

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

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

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

2.3 非線性方法

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

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

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

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

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

3 時(shí)頻分辨率分析

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

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

(1)

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

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

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

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

(2)

信號(hào)的瞬時(shí)頻率為

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

(3)

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

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

方法STFTCWTSTGSTWVDPWVDSPWVDHHTSSTCDGLCT熵4.98855.10514.89125.07694.72914.52154.62670.04621.71213.22114.7242

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

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

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

(4)

信號(hào)的兩個(gè)瞬時(shí)頻率分量為

(5)

兩分量在0.6s前相互緊貼，隨后瞬時(shí)頻率變化趨勢(shì)逐漸分開(kāi)，信號(hào)后半段頻率變化劇烈。

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

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

方法STFTCWTSTGSTWVDPWVDSPWVDHHTSSTCDGLCT熵3.87403.49374.01574.05104.79724.11414.02880.99062.26422.75364.3537

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

3.3 模型3：合成地震信號(hào)

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

圖3 合成地震信號(hào)、單頻子波以及不同方法獲得的時(shí)頻譜

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

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

表4 合成地震信號(hào)不同方法時(shí)頻譜的Rényi熵

4 計(jì)算效率分析

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

圖4 不同時(shí)頻分析方法三個(gè)模型信號(hào)處理效率對(duì)比

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

5 抗噪性能分析

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

(6)

式中fc為信號(hào)的Nyquist頻率。通過(guò)與理論瞬時(shí)頻率曲線IF(t)相比較可計(jì)算IFe(t)的均方誤差

(7)

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

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

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

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

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

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

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

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

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

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

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

7 結(jié)論

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

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

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

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

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