姚家駿，楊立明，馮建剛

（1.青海省地震局，青海 西寧 810000；2.中國(guó)地震局蘭州地震研究所，甘肅 蘭州 730000）

0 引言

對(duì)于頻率等特征量不隨時(shí)間變化的平穩(wěn)信號(hào)，傳統(tǒng)的傅立葉變換就能很好地建立時(shí)域與頻域的映射關(guān)系。而對(duì)于地震波等非平穩(wěn)信號(hào)，由于其頻率等特征量隨時(shí)間變化而變化，傳統(tǒng)的傅立葉變換因?yàn)椴荒芊从承盘?hào)的時(shí)頻局域性信息不再適用。用時(shí)頻分析法建立起時(shí)間和頻率的二維函數(shù)將一維的時(shí)間域信號(hào)（如數(shù)字化地震波記錄）擴(kuò)展到二維的時(shí)頻平面上，直觀而全面地反映信號(hào)頻率隨時(shí)間的變化關(guān)系，使我們同時(shí)掌握信號(hào)的時(shí)域及頻域局域性信息。

時(shí)頻分析方法可分為線性與非線性兩種［1］。其中線性時(shí)頻表示是由傅立葉變換轉(zhuǎn)化而來，典型的線性變換形式有短時(shí)傅立葉變換（STFT）、小波變換（WT）及S變換。常用的非線性時(shí)頻表示有維格納－威爾分布（WVD）、喬伊－威廉斯分布（CWD）、趙－阿特拉斯－馬克斯時(shí)頻分布（ZAM）等。非平穩(wěn)信號(hào)分析與處理的研究工作主要是從20世紀(jì)40年代開始的。1946年Gabor［2］首次發(fā)現(xiàn)在信號(hào)處理的過程中時(shí)頻分析的重要性，隨后他便提出用加窗傅立葉變換（STFT）法去找出聲波頻率的確切位置；Wigner于1932年首先提出了Wigner分布的概念，并把它用于量子力學(xué)領(lǐng)域；1948年Ville［3］首先把它應(yīng)用于信號(hào)分析，因此Wigner分布又稱Wigner－Ville分布，簡(jiǎn)稱為 WVD；1966年Cohen［4］給出了各種雙線性時(shí)頻分布的統(tǒng)一表示形式；1996年Stockwell提出一種介于STFT和小波變換之間的S變換［5］，吸取了兩者的優(yōu)點(diǎn)。近幾年來人們利用時(shí)頻分析技術(shù)識(shí)別弱震相、區(qū)分人工地震與天然地震、研究大震前小震的頻譜變化等已經(jīng)做了一些工作［6－11］。本文試用短時(shí)傅立葉變換（STFT）、S變換、CWD分布以及人們較少關(guān)注的ZAM分布分析數(shù)字地震波P波信號(hào)、S波信號(hào)的時(shí)頻分布，分析它們?cè)诜直娴卣鸩〞r(shí)的應(yīng)用效果，總結(jié)其優(yōu)缺點(diǎn)。

1 時(shí)頻分析的算法原理

1.1 短時(shí)傅立葉變換（STFT）

短時(shí)傅立葉變換（STFT）思路依托于處理平穩(wěn)信號(hào)的經(jīng)典傅立葉變換，其基本思想［1］是把信號(hào)劃分成許多小的時(shí)間間隔，用傅立葉變換分析每一個(gè)時(shí)間間隔，以便確定每個(gè)時(shí)間間隔所存在的頻率，沿時(shí)間發(fā)展的這些頻率的總體就表現(xiàn)頻譜在時(shí)間上的變化過程。短時(shí)傅立葉變換概念直接，算法簡(jiǎn)單，是其他時(shí)頻方法的基礎(chǔ)。

信號(hào)s（t）的STFT 定義［12］為

其中h（t）是窗函數(shù)，反變換的公式為

由式（1）可知，短時(shí)傅立葉變換就是在傅立葉變換的基礎(chǔ)上通過固定長(zhǎng)度的滑動(dòng)窗來計(jì)算信號(hào)頻譜的。但該滑動(dòng)窗口與頻率無關(guān)，它的時(shí)間分辨率和頻率分辨率受Heisenberg測(cè)不準(zhǔn)原理約束，利用短窗口，時(shí)間分辨率比較好，但是頻率分辨能力差；當(dāng)利用長(zhǎng)窗口時(shí)，有較高的頻率分辨率，但時(shí)間分辨能力就弱。

1.2 S變換

信號(hào)s（t）的S變換（ST）的定義［5］為

式中τ、f分別表示時(shí)間和頻率，均為實(shí)數(shù)。反變換為

它是以Morlet小波為基本小波的連續(xù)小波變換的延伸?；拘〔ㄊ怯珊?jiǎn)諧波與高斯函數(shù)的乘積構(gòu)成的，定義為

基本小波中的簡(jiǎn)諧波在時(shí)間域僅作伸縮變換，而高斯函數(shù)則進(jìn)行伸縮和平移。由式（2）可知，S變換也可以看作是在傅立葉變換的基礎(chǔ)上通過滑動(dòng)窗來計(jì)算信號(hào)頻譜的，但滑動(dòng)窗口與頻率有關(guān)。S變換采用寬度可變的高斯窗函數(shù)，其時(shí)窗寬度隨頻率呈反比變化，在低頻段的時(shí)窗較寬，可獲得較高的頻率分辨率；而高頻段的時(shí)窗較窄，可獲得很高的時(shí)間分辨率。因此，較之STFT，S變換具有更好的時(shí)頻分辨率，特別是在高頻部分S變換的時(shí)間分辨率優(yōu)勢(shì)更加明顯。而在連續(xù)小波變換中簡(jiǎn)諧波與高斯函數(shù)進(jìn)行同樣的伸縮和平移。因此與小波變換、短時(shí)Fourier變換等時(shí)—頻域方法相比，S變換有其獨(dú)特的優(yōu)點(diǎn)，如信號(hào)的S變換的時(shí)頻譜分辨率與頻率（即尺度）有關(guān)，且與其Fourier譜保持直接的聯(lián)系，基本小波不必滿足容許性條件［5］等，這些特點(diǎn)在實(shí)際應(yīng)用中是非常有用的。

1.3 喬伊－威廉斯分布（CWD）

非線性時(shí)頻變換中最基本的的Wigner－Ville分布（WVD）［5］為

其中＊表示解析信號(hào)s（t）的復(fù)共軛，s（t）在信號(hào)分時(shí)出現(xiàn)兩次，因此稱其為雙線性變換。由于式中不包含任何的窗函數(shù)，避免了線性時(shí)頻表示中時(shí)間、頻率分辨率相互牽制的矛盾。雖然WVD的時(shí)頻聚集度最高，能滿足的數(shù)學(xué)性質(zhì)最多，但多信號(hào)的WVD會(huì)出現(xiàn)交叉項(xiàng)，在不少場(chǎng)合會(huì)限制其效果。后來人們又提出多種改進(jìn)形式，如指數(shù)分布、廣義指數(shù)分布、廣義雙線性時(shí)頻分布（Cohen類時(shí)頻分布）等，因此 WVD被譽(yù)為 “所有時(shí)頻分布之母［13］”。

20世紀(jì)60年代中期，Cohen將眾多的雙線性時(shí)頻分布統(tǒng)一為一種形式，習(xí)慣稱之為Cohen類時(shí)頻分布，定義［14］為

其中φ（θ，τ）叫做核函數(shù)。

選取不同的核函數(shù)可得到不同的時(shí)頻表示.但核函數(shù)在減少交叉項(xiàng)的同時(shí)也會(huì)帶來負(fù)面作用：信號(hào)自身項(xiàng)展寬，信號(hào)分辨率的降低。抑制交叉項(xiàng)與提高信號(hào)分辨率，如同短時(shí)傅立葉變換中時(shí)域和頻域分辨率一樣，不可能同時(shí)做到最好，只能折中最優(yōu)。目前大多抑制交叉項(xiàng)的方法是以降低分辨率為代價(jià)的。分析地震信號(hào)時(shí)頻分布，需要尋找一個(gè)適應(yīng)于地震信號(hào)性質(zhì)的濾波核函數(shù)，最大地抑制交叉項(xiàng)同時(shí)盡可能小地降低分辨率。對(duì)于喬伊－威廉斯時(shí)頻分布，取核函數(shù)

CWD可以有效的抑制由不同時(shí)間和不同頻率中心引起的交叉干擾，但是它不能減少同一頻率不同時(shí)間中心或是同一時(shí)間不同頻率中心之間的干擾。而且，CWD的計(jì)算機(jī)處理速度較慢，交叉項(xiàng)的出現(xiàn)常常導(dǎo)致時(shí)頻平面上出現(xiàn)偽影現(xiàn)象。

1.4 趙－阿特拉斯－馬克斯時(shí)頻分布（ZAM）

若取核函數(shù)

則Cohen類時(shí)頻分布變?yōu)橼w－阿特拉斯－馬克斯時(shí)頻分布（ZAM）。在趙－阿特拉斯－馬克斯原著中g(shù)（τ）＝1，α＝1／2。研究表明［15］類似形如式（8）的錐形核函數(shù)很好的適合于地震數(shù)據(jù)分析。這個(gè)核函數(shù)不像頻譜圖那樣消除脈沖，但可完全地抑制它們，從而得到一個(gè)有效和準(zhǔn)確的頻率局部化表示。

2 算例分析

選用四個(gè)基本高斯元構(gòu)成一個(gè)合成信號(hào)，時(shí)間長(zhǎng)度為128s，其特征見表1。

表1 四個(gè)基本高斯元的特征情況

為了較好地在不同TFD方法下考察信號(hào)間的相干擾程度，設(shè)計(jì)出一對(duì)信號(hào)在同一時(shí)間中心，另一對(duì)在同一頻率中心，并且相鄰（圖1）。其中圖1（a）為合成信號(hào)的時(shí)域波形，圖1（b）為合成信號(hào)理想情況的時(shí)頻分布圖。

圖1 合成信號(hào)及其理論上的時(shí)頻分布Fig.1 Synthesis signal and its academic TFD.

合成信號(hào)的短時(shí)傅立葉變換（STFT）、S變換、喬伊－威廉斯時(shí)頻分布（CWD）、趙－阿特拉斯－馬克斯時(shí)頻分布（ZAM）如圖2所示。由于STFT的時(shí)頻分辨率與窗函數(shù)長(zhǎng)度有很大的關(guān)系，我們經(jīng)過實(shí)踐，選取合理的長(zhǎng)度，以便同時(shí)在時(shí)間分辨率、頻率分辨率上取個(gè)折中。由圖2（a）、圖2（b）的色彩顯示可以看出，對(duì)于STFT來說，信號(hào)頻譜圖在頻率方向上對(duì)應(yīng)的寬度是不變的，即在低頻和高頻的頻率分辨率是不變的，不利于低頻、高頻信號(hào)的分辨。對(duì)于S變換來說，信號(hào)頻譜圖在頻率方向上對(duì)應(yīng)的寬度是不同的，可看到在信號(hào)頻率變化處有較明顯的邊界，能把高頻信號(hào)分離出，并且不存在交叉干擾。S變換在低頻部分具有較高的頻率分辨率，高頻部分具有較高的時(shí)間分辨率。恰恰能滿足法國(guó)地球物理學(xué)家Morlet于20世紀(jì)80年代初在分析地震信號(hào)時(shí)提出的“地震信號(hào)的低頻端應(yīng)該具有很高的頻率分辨率，而在高頻端頻率分辨率可以較低”的要求。但信號(hào)的時(shí)頻聚集性較差。圖2（c）所示CWD時(shí)頻分布盡管具有較好的時(shí)頻聚集性，有明顯的交叉項(xiàng)干擾。信號(hào)的ZAM時(shí)頻分布（圖2（d））同樣具有較好的時(shí)頻聚集性，而且和已知仿真信號(hào)時(shí)間起始點(diǎn)和頻率點(diǎn)相比，交叉項(xiàng)干擾要小于CWD時(shí)頻分布結(jié)果，但時(shí)頻聚集性稍微不如CWD。

3 數(shù)字化地震波記錄時(shí)頻分析

甘肅省平?jīng)鰯?shù)字地震臺(tái)記錄的2009年11月28日四川成都彭州市MS5.0地震，震相清晰無干擾。儀器的采樣率為100。截取其UD道P波、S波，如圖3、圖4，用上述四種方法做TFD。由圖3（a）可以看出，該P(yáng)波信號(hào)的11.8s到13.5s、18.5～19.5s之間是振幅較大的區(qū)域。由圖3（b）的傅立葉譜可以看出，該P(yáng)波信號(hào)的主頻率為1.9Hz，主要分布在1.5～3Hz之間。由圖4（a）可以看出，該S波信號(hào)的5.5～7.1s、17.5s到20s之間是振幅較大的區(qū)域。由圖4（b）的傅立葉譜可以看出，該S波信號(hào)的主頻率為0.7～0.9Hz，主要分布在1.7 Hz之內(nèi)。圖5可知：四種方法都能較準(zhǔn)確地檢測(cè)出該信號(hào)的主頻率的時(shí)頻位置，但S變換能檢8Hz的高頻弱信號(hào)，其余三種方法只能檢測(cè)到4～5Hz左右的高頻信號(hào)；非線性的CWD法與ZAM法的時(shí)頻聚集性強(qiáng)，主信號(hào)的時(shí)頻分布情況一目了然；CWD時(shí)頻分布盡管具有較好的時(shí)頻聚集性，但交叉項(xiàng)干擾比較明顯，容易出現(xiàn)偽影現(xiàn)象，產(chǎn)生虛假頻率；ZAM法雖時(shí)頻聚集性稍遜于CWD時(shí)頻分布，但在消除交叉項(xiàng)的干擾的優(yōu)勢(shì)上比較明顯。

圖2 合成信號(hào)的時(shí)頻分布Fig.2 TFD of the synthesis signal（（a）STFT；（b）ST；（c）CWD；（d）ZAM）.

圖3 地震波P波信號(hào)及其傅立葉譜Fig.3 P wave（a）and its Fourier spectrum （b）.

圖4 S波及其傅立葉譜Fig.4 S wave（a）and its Fourier spectrum （b）.

圖5 P波的時(shí)頻分布Fig.5 TFD of P wave（（a）STFT；（b）ST；（c）CWD；（d）ZAM）.

圖6為該S波信號(hào)經(jīng)短時(shí)傅立葉變換（STFT）、S變換、喬伊－威廉斯時(shí)頻分布（CWD）、趙－阿特拉斯－馬克斯時(shí)頻分布（ZAM）所做得時(shí)頻分布。由圖6可知：四種方法都能較準(zhǔn)確地檢測(cè)出該信號(hào)的主頻率的時(shí)頻位置，但S變換能檢2～3 Hz的高頻弱信號(hào)，其余三種方法只能檢測(cè)到2s左右出現(xiàn)的1.7Hz左右的高頻主信號(hào)；非線性的CWD法與ZAM法的時(shí)頻聚集性的優(yōu)勢(shì)依然很明顯；但仍然知道CWD時(shí)頻分布雖然具有較好的時(shí)頻聚集性，可交叉項(xiàng)干擾比較明顯的；相對(duì)于CWD法，ZAM法具有消除交叉項(xiàng)的干擾的優(yōu)勢(shì)。

4 結(jié)論

通過對(duì)模擬信號(hào)以及平?jīng)鰯?shù)字地震臺(tái)記錄的2009年11月28日在四川彭州市MS5.0地震P波、S波的分析可知：

（1）四種方法都能比較精確地檢測(cè)出地震信號(hào)的主頻率的時(shí)頻位置，二維時(shí)頻分布圖能充分顯示出地震波信號(hào)的非穩(wěn)定性特征；

（2）線性的STFT法運(yùn)算速度最快，但時(shí)頻分辨率表現(xiàn)最差；S變換在檢測(cè)高頻弱信號(hào)方面有獨(dú)特的優(yōu)勢(shì)，但時(shí)頻聚集性較差；

（3）在對(duì)信號(hào)的時(shí)頻聚集性方面，非線性的CWD法與ZAM法遠(yuǎn)比線的STFT、S變換方法要強(qiáng)；

（4）CWD時(shí)頻分布盡管具有較好的時(shí)頻聚集性，但交叉項(xiàng)干擾比較明顯，容易出現(xiàn)比較明顯的偽影現(xiàn)象，產(chǎn)生虛假頻率，ZAM法雖時(shí)頻聚集性稍遜于CWD時(shí)頻分布，但在消除交叉項(xiàng)的干擾方面要遠(yuǎn)比CWD法出色。

數(shù)字化地震觀測(cè)技術(shù)的普及與提高，使現(xiàn)代信號(hào)處理技術(shù)能得以很好的比較與應(yīng)用。時(shí)頻分析可以在時(shí)頻二維平面很好地解釋地震波的非平穩(wěn)性，可以更豐富地顯示地震波特征，因此時(shí)頻分析在地震震相識(shí)別、模式識(shí)別以及利用小震頻譜特征預(yù)報(bào)大地震等方面有廣泛的應(yīng)用前景。

圖6 S波的時(shí)頻分布Fig.6 TFD of S wave（（a）STFT；（b）ST；（c）CWD；（d）ZAM）.

［1］L科恩（美）.時(shí)頻分析－理論與應(yīng)用［M］.西安：西安交通大學(xué)出版社，1998.

［2］Gabor D.Theory of communication［J］.Journal of the IEEE，1946，93：488－497.

［3］Ville J.Theorie et application de la notion de signal analytique［J］.Cables et transmissions，1948，2A：61－74.

［4］Cohen L.Generalized phase－space distribution functions［J］.Jour.Math.Phys.，1966，7：781－786.

［5］R Stockwell，L Mansinha，R Lowe.Localization of the Complex Spectrum：The S Transform［J］.IEEE Trans.on Signal Process，1996，44（4）：998－1001.

［6］許康生.地震信號(hào)的時(shí)頻分析［J］.西北地震學(xué)報(bào)，2000，22（4）：479－482.

［7］劉希強(qiáng)，周惠蘭，李紅.基于小波包變換的地震數(shù)據(jù)時(shí)頻分析方法［J］.西北地震學(xué)報(bào)，2000，22（2）：143－146.

［8］劉希強(qiáng)，沈萍，山長(zhǎng)侖，等.數(shù)字化地震波形資料的時(shí)頻分析方法及應(yīng)用［J］.西北地震學(xué)報(bào)，2004，26（2）：118－125.

［9］熊曉軍，賀振華，黃德濟(jì)，等.廣義S變換在地震高分辨處理中的應(yīng)用［J］.勘探地球物理進(jìn)展，2006，29（6）：415－418.

［10］陳學(xué)華，賀振華.改進(jìn)的S變換及在地震信號(hào)處理中的應(yīng)用［J］.數(shù)據(jù)采集與處理，2005，20（4）：449－453.

［11］嚴(yán)鋒，靳平，范廣超.地震臺(tái)陣上信號(hào)方位角和慢度的時(shí)、頻域估計(jì)方法比較［J］.西北地震學(xué)報(bào)，2006，28（4）：227－230.

［12］趙淑紅.時(shí)頻分析方法及其在地震數(shù)據(jù)處理中的應(yīng)用［D］.西安：長(zhǎng)安大學(xué)，2006.

［13］張建，沈民奮，宋驥.三種時(shí)頻分析方法在心音信號(hào)分析中的應(yīng)用［J］.汕頭大學(xué)學(xué)報(bào)（自然科學(xué)版），2003，18（2）：62－69.

［14］王紅萍.非平穩(wěn)信號(hào)時(shí)頻分析方法性能研究［D］.南京：南京航空航天大學(xué)，2008.

［15］張曄.信號(hào)時(shí)頻分析及應(yīng)用［M］.哈爾濱：哈爾濱工業(yè)大學(xué)出版社，2006.