時(shí)頻分析方法在形變數(shù)據(jù)中的應(yīng)用研究①

王 寧，吳 云，張 燕

（1.中國(guó)地震局地震研究所地震大地測(cè)量重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430071；2.中國(guó)地震局地殼應(yīng)力研究所武漢科技創(chuàng)新基地，湖北 武漢 430071）

0 引言

傅里葉變換的提出為信號(hào)處理帶來了翻天覆地的變化，其應(yīng)用遍及地震信號(hào)分析［1］、地質(zhì)勘探［2］、機(jī)械振動(dòng)信號(hào)分析［3］等各個(gè)領(lǐng)域。然而，經(jīng)典傅里葉變換不能同時(shí)兼?zhèn)鋾r(shí)間和頻率分辨率，不適宜處理諸如非平穩(wěn)的地震信號(hào)。時(shí)頻分析方法將傳統(tǒng)傅里葉變換的整體譜推廣到局部譜中，對(duì)于非平穩(wěn)地震信號(hào)的分析具有很好的適用性。隨著Gabor變換、小波變換、Hilbert－Huang變換等時(shí)頻分析方法的出現(xiàn)，非平穩(wěn)信號(hào)的時(shí)頻特性得到了很好的展示。周摯等提出了固體潮時(shí)頻分析新方法，即HHT方法，同時(shí)指出了該方法的特點(diǎn)和局限性［4］；周摯等基于HHT提取昆明、下關(guān)重力固體潮的地震前兆信息，設(shè)計(jì)了重力固體潮地震前兆分析的特征參數(shù)，得到震前異常特征參數(shù)為HHT計(jì)算得到的實(shí)際參數(shù)和理論參數(shù)的頻率比值［5］；呂品姬等研究了小波分解及短時(shí)傅里葉變換在地形變觀測(cè)數(shù)據(jù)中的應(yīng)用，得到了不同采樣率的數(shù)據(jù)其時(shí)頻圖的特點(diǎn)及所反應(yīng)的物理事件［6］。然而時(shí)頻分析方法在地形變數(shù)據(jù)分析中尚處于起步階段，還沒有得到普遍應(yīng)用。

本文通過小波變換、Hilbert－Huang變換及S變換在設(shè)計(jì)數(shù)據(jù)以及地形變數(shù)據(jù)分析中的應(yīng)用，分析比較三種方法，以得到這些方法的特點(diǎn)和局限性，為以后方法的利用提供依據(jù)。

1 算法

1.1 小波變換

小波變換最早由法國(guó)的地球物理學(xué)家于上世紀(jì)80年代提出，因其兼具了時(shí)間和頻率的分辨率，被稱為聯(lián)合時(shí)頻分析，簡(jiǎn)稱時(shí)頻分析。信號(hào)h（t）的小波變換定義為

式中，a表示尺度因子；b表示平移因子；ψ（t）為基本小波或母小波。基本小波必須滿足容許性條件，即

小波變換的母小波具有多種形式，本文采用的是分析非平穩(wěn)信號(hào)常用的cmor3－3小波。

1.2 S變換

S變換是以Morlet小波為基本小波的連續(xù)小波變換的延伸［7］，它結(jié)合了短時(shí)傅里葉變換和小波變換的優(yōu)點(diǎn)，解決了短時(shí)傅里葉變換窗口頻率不能調(diào)節(jié)的問題，同時(shí)具備了小波變換多分辨率的優(yōu)勢(shì)，而且基本小波不用滿足容許性條件。信號(hào)h（t）的S變換定義為［8］

S變換不同于小波變換，采用與頻率有關(guān)的可變高斯窗函數(shù)，其時(shí)頻分辨率與頻率有關(guān)，信號(hào)的變換結(jié)果與其傅里葉譜保持直接聯(lián)系。因此通過S變換及其逆變換可以實(shí)現(xiàn)信號(hào)在時(shí)－頻域的無損轉(zhuǎn)換，這些特點(diǎn)使得S變換在分析非平穩(wěn)信號(hào)過程中得到了廣泛的應(yīng)用［9］。

1.3 Hilbert－Huang變換［10］

Hilbert－Huang變換是一種基于經(jīng)驗(yàn)的數(shù)據(jù)分析方法。其核心思想是把非線性非平穩(wěn)信號(hào)看成是由若干的本征函數(shù)構(gòu)成，通過經(jīng)驗(yàn)?zāi)B(tài)分解（EMD）將其分解出來，然后利用Hilbert變換構(gòu)造解析信號(hào)，得到數(shù)據(jù)的瞬時(shí)屬性，進(jìn)而得到其時(shí)頻譜。

1.3.1 經(jīng)驗(yàn)?zāi)B(tài)分解（EMD）

經(jīng)驗(yàn)?zāi)B(tài)分解建立的前提是假設(shè)任何數(shù)據(jù)都是由一系列的本征模態(tài)函數(shù)（IMF）構(gòu)成，IMF需滿足以下兩個(gè)條件：（1）整個(gè)數(shù)據(jù)中，極值點(diǎn)和過零點(diǎn)的數(shù)量應(yīng)保持一致或最多相差一個(gè)；（2）在任意點(diǎn)處，由極大值和極小值定義的包絡(luò)線的均值應(yīng)該為0。

EMD方法包括以下過程：首先，找出原始數(shù)據(jù)序列x（t）的所有極大值和極小值，并用三次樣條曲線擬合其上下包絡(luò)線；用原始數(shù)據(jù)x（t）減去上下包絡(luò)線的均值m1（t）得到第一個(gè)分量h1（t），一般情況下h1（t）不滿足IMF條件，需要重復(fù)以上處理過程；經(jīng)過k次處理后產(chǎn)生第一個(gè)IMF分量c1（t），即

其次，將原始數(shù)據(jù)序列x（t）減去c1（t），得到余量r2（t）；重復(fù)以上步驟得到rn（t）。最后，可知原始數(shù)據(jù)可以由以上IMF分量及余量疊加而成。即

經(jīng)EMD分解變換得到的IMF序列是直接從原始時(shí)序數(shù)據(jù)中分離出來的，事先無需確定分解階次，不受人為因素影響，不存在機(jī)械分解。因此IMF序列能更好反映原始數(shù)據(jù)固有的物理特性，其分解是客觀的、內(nèi)在的和自適應(yīng)的［11］。

1.3.2 Hilbert變換及瞬時(shí)屬性提取

將每一個(gè)IMF進(jìn)行Hilbert變換，進(jìn)而計(jì)算出信號(hào)的瞬時(shí)頻率，通過合成得到的就是HHT時(shí)頻譜。

2 實(shí)驗(yàn)數(shù)據(jù)分析

實(shí)驗(yàn)將設(shè)計(jì)的仿真信號(hào)分別進(jìn)行S變換、Morlet小波變換以及HHT變換，得到這三種方法的時(shí)頻圖，研究它們?cè)诜治鰯?shù)據(jù)上的適用范圍。

2.1 單頻率的正弦波

圖1 正弦曲線的時(shí)頻圖Fig.1 Time－frequency spectrum of sinusoid

地形變數(shù)據(jù)具有典型的波動(dòng)特性，圖1給出了頻率為1Hz、出現(xiàn)在40～50s期間的正弦波曲線。圖1（d）中為顯示HHT時(shí)頻圖的細(xì)節(jié)信息，將其40～55s處圖形放大。從圖可以看出，3種時(shí)頻分析方法在40～50s時(shí)刻都能很好地顯示正弦信號(hào)的頻率信息，且都具有很好的時(shí)間和頻率分辨率。其中，S變換其頻率分辨率較低，圖上顯示為圖譜較寬；小波變換低頻處存在虛假頻率；HHT方法低頻擾動(dòng)和頻率分辨率都較好，但其端部效應(yīng)明顯。

2.2 不連續(xù)的正弦波

為驗(yàn)證三種時(shí)頻方法對(duì)突變頻率的檢驗(yàn)，圖2采用兩個(gè)不同頻率（f1＝1Hz，f2＝2Hz）信號(hào)合成的信號(hào)進(jìn)行對(duì)比分析，比較在頻率變化處三種方法的優(yōu)劣。

從圖2可以看出，三種方法都能很好地在頻率突變處對(duì)不同頻率進(jìn)行區(qū)分，其中小波變換和S變換得到的效果圖由于存在時(shí)間和頻率分辨率的問題，顯然在兩個(gè)不同的頻率處的分辨率不同，高頻部分的頻率帶變寬，尤其S變換高頻處的頻率分辨率更低，HHT則不存在分辨率不統(tǒng)一的問題。但三種方法在頻率的突變處都存在交叉影響，端部效應(yīng)也比較明顯。

圖2 3種方法對(duì)突變頻率的檢驗(yàn)Fig.2 Three methods for the inspection of mutation frequency

2.3 多頻率成分?jǐn)?shù)據(jù)

地形變數(shù)據(jù)同一時(shí)間往往具有多個(gè)頻率成分。為驗(yàn)證這一特點(diǎn)，圖3采用頻率隨時(shí)間變化較大的信號(hào)，也是經(jīng)典的S變換采用的實(shí)例信號(hào)進(jìn)行比較。

其中l(wèi)en為一常量，表示采樣的數(shù)據(jù)點(diǎn)數(shù)。

從對(duì)具有多頻段的模擬信號(hào)的分析來看，S變換和小波變換可以顯示出信號(hào)在某段時(shí)間內(nèi)的頻率及其能量的變化趨勢(shì)，而HHT由于每一個(gè)時(shí)刻由不同的點(diǎn)表示不同的頻率，因而更適合分析信號(hào)的瞬時(shí)特性。S變換在低頻部分的分辨率明顯優(yōu)于小波變換的分辨率，小波變換容易產(chǎn)生低頻的虛假信號(hào)。

總體來看，在做連續(xù)時(shí)間序列的時(shí)頻分析中，HHT方法不適合做相關(guān)分析。小波變換與S變換均可作為時(shí)頻分析的有效工具，應(yīng)根據(jù)所分析對(duì)象包含的頻率信息選擇適當(dāng)方法，低頻信息較多地采用S變換，高頻信息較多且頻率成分較多地采用小波分析方法；在分析高低頻率混雜的數(shù)據(jù)時(shí)，可先做 小波分析，進(jìn)一步針對(duì)某一低頻事件做S分析。

圖3 多頻率成分的時(shí)頻分析比較Fig.3 Comparison of multiple frequency data＇s time－frequency analysis

3 地形變觀測(cè)數(shù)據(jù)時(shí)頻分析

在對(duì)上述三種時(shí)頻分析方法各自的優(yōu)越性及適用范圍進(jìn)行對(duì)比分析的基礎(chǔ)上，本文采用部分地形變數(shù)據(jù)作進(jìn)一步的分析。

3.1 寬頻帶數(shù)據(jù)應(yīng)用

研究資料為汶川地震前前1小時(shí)距離較近的郫縣臺(tái)和安康臺(tái)的寬頻帶數(shù)據(jù)（數(shù)據(jù)采樣頻率為2 Hz）。汶川地震震中經(jīng)度為103.4°，緯度為31°。選擇的兩個(gè)臺(tái)站的位置和汶川地震的震中距如表1。

表1 臺(tái)站的震中距Table 1 Position of the stations

對(duì)兩個(gè)臺(tái)站的寬頻帶數(shù)據(jù)分別使用三種時(shí)頻分析方法進(jìn)行處理，結(jié)果見圖4、5。看出在分析寬頻帶數(shù)據(jù)時(shí)，HHT不能揭示信號(hào)隨時(shí)間的頻率變化趨勢(shì)，S變換和小波變換顯示了地震信號(hào)存在0.1～0.3Hz的頻率信息。還可以看出，S變換較小波變換在高頻部分具有很好的時(shí)間分辨率，低頻部分具有很好的頻率分辨率，且產(chǎn)生的虛假頻率較少，但其高頻部分的頻率分辨率略低于小波變換。

秒采樣數(shù)據(jù)量較大，在做時(shí)頻分析時(shí)得到的是震前一小時(shí)的圖譜，可以看出通過秒采樣得到的圖譜震前異常不明顯，與臺(tái)站距離地震的遠(yuǎn)近關(guān)系不大。震前異常尚有待進(jìn)一步的研究。

3.2 洞體應(yīng)變數(shù)據(jù)的應(yīng)用

選取了2008年3月到6月黃梅洞體應(yīng)變NS分量的數(shù)據(jù)及黃梅臺(tái)水管傾斜NS分量的數(shù)據(jù)，采樣周期為小時(shí)。三種時(shí)頻分析方法的結(jié)果見圖6、7。可以看出，S變換和小波變換時(shí)頻圖上的固體潮日波，半日波較明顯，且與傅里葉譜的頻率值對(duì)應(yīng)；HHT譜不能用在揭示這種長(zhǎng)周期性連續(xù)性事件上；S變換在低頻部分的分辨率顯然高于小波變換，產(chǎn)生的虛假頻率較少。從小波和S變換圖中可以看到5月12日前產(chǎn)生了一段低于日波頻率的信號(hào)，且信號(hào)強(qiáng)度在臨震時(shí)逐漸加強(qiáng)，這是否與2012年5月12日汶川地震有關(guān)，有待進(jìn)一步的證實(shí)。

4 結(jié)果與討論

（1）從仿真信號(hào)的時(shí)頻圖分析，小波變換、S變換以及Hilbert－Huang變換都能很好地展示信號(hào)頻率隨時(shí)間的變化。

圖4 安康臺(tái)寬頻帶數(shù)據(jù)的時(shí)頻分析比較Fig.4 Comparison of time－frequency analysis based on the broadband data from Ankang station

圖5 郫縣臺(tái)寬頻帶數(shù)據(jù)的時(shí)頻分析比較Fig.5 Comparison of time－frequency analysis based on the broadband data from Pixian station

（2）由地形變數(shù)據(jù)的小波變換結(jié)果看：小波變換的頻譜在低頻部分會(huì)產(chǎn)生虛假頻率信息，但在高頻部分的頻率分辨率較高，表現(xiàn)為中心頻率集中；S變換時(shí)間和頻率的分辨率都比較好，尤其是在低頻部分具有很好的分辨率，但其在高頻部分的頻率分辨率較低，需要采用必要的方法進(jìn)行改進(jìn)，對(duì)低頻成分比較多的信號(hào)S變換優(yōu)勢(shì)明顯；Hilbert－Huang變換從數(shù)據(jù)本身進(jìn)行分解得到，保存了信號(hào)本身的物理特性，同時(shí)也具有很好的時(shí)間和頻率分辨率，其頻率的分辨率統(tǒng)一，不存在不同頻率處分辨率不同的影響，但存在一定的邊緣效應(yīng)，且其和S變換一樣受各種條件限制，目前處理的數(shù)據(jù)量十分有限，同時(shí)在揭示趨勢(shì)變化上存在缺陷，還有待進(jìn)一步的改進(jìn)。目前提出的EEMD方法是在Hilbert－Huang變換基礎(chǔ)上發(fā)展起來的，其使用優(yōu)勢(shì)還有待進(jìn)一步的驗(yàn)證［12］。因此，在對(duì)形變數(shù)據(jù)分析時(shí)，針對(duì)高頻部分可以采用小波分析的方法，低頻部分用S變換進(jìn)行細(xì)節(jié)分析。分析信號(hào)組成成分特征時(shí)，可以采用HHT方法將信號(hào)進(jìn)行分解，依次分析各分量的物理特性。

圖6 黃梅臺(tái)應(yīng)變儀數(shù)據(jù)時(shí)頻圖比較Fig.6 Comparison of time－frequency analysis based on the strainmeter＇s data from Huangmei station

圖7 黃梅臺(tái)傾斜儀數(shù)據(jù)時(shí)頻圖比較Fig.7 Comparison of time－frequency analysis based on the tiltmeter＇s data from Huangmei station

（3）在針對(duì)整時(shí)數(shù)據(jù)的震前后時(shí)頻分析中，可以看到存在一定的頻率異常信號(hào)，但其是否屬于震前異常，需要進(jìn)一步的證明。

（References）

［1］楊培杰，印興耀，張廣智.希爾伯特—黃變換地震信號(hào)時(shí)頻分析與屬性提取［J］.地球物理學(xué)進(jìn)展，2007，22（5）：1585－1590.YANG Pei－jie，YIN Xing－yao，ZHANG Guang－zhi.Seismic Signal Time－frequency Analysis and Attributes Extraction Based on HHT［J］.Progress in Geophysics，2007，22（5）：1585－1590.（in Chinese）

［2］楊志強(qiáng)，單娜琳，劉占興，等.S變換在巖溶區(qū)地震映像資料處理中的應(yīng)用［J］.工程地球物理學(xué)報(bào)，2012，9（2）：227－230.YANG Zhi－qiang，DAN Na－lin，LIU Zhan－xing，et al.Application of S－Transform to Seismic Lmaging Data Processing in Karst Areas［J］.Chinese Journal of Engineering Geophysics，2012，9（2）：227－230.（in Chinese）

［3］楊世錫，胡勁松，吳昭同，等.旋轉(zhuǎn)機(jī)械振動(dòng)信號(hào)基于EMD的希爾伯特變換和小波變換時(shí)頻分析比較［J］.中國(guó)電機(jī)工程學(xué)報(bào)，2003，23（6）：102－107.YANG Shi－xi，HU Jin－song，WU Zhao－tong，et al.The Comparison of Vibration Signals’Time－frequency Analysis Between EMD－Based HT and WT Method in Rotating Machinery［J］.Proceedings of the CSEE，2003，23（6）：102－107.（in Chinese）

［4］周摯，山秀明，梁虹，等.固體潮時(shí)頻分析新方法［J］.地震研究，2005，28（4）：334－339.ZHOU Zhi，SHAN Xiu－ming，LIANG Hong，et al.A New Approach of Earth Tides Time－Frequency Analysis［J］.Journal of Seismological Research，2005，28（4）：334－339.（in Chinese）

［5］周摯，山秀明，張立，等.基于HHT提取昆明，下關(guān)重力固體潮的地震前兆信息［J］.地球物理學(xué)報(bào)，2008，51（3）：836－844.ZHOU Zhi，SHAN Xiu－ming，ZHANG Li，et al.The Gravity Tide of Kunming ＆ Xiaguan Based on the HHT［J］.Chinese Journal of Geophysics，2008，51（3）：836－844.（in Chinese）

［6］呂品姬，趙斌，陳志遙，等.小波分解－STFT方法在地形變觀測(cè)數(shù)據(jù)中的應(yīng)用［J］.大地測(cè)量與地球動(dòng)力學(xué)，2011，31（5）：136－140.LV Pin－ji，ZHAO Bin，CHEN Zhi－yao，et al.Application of Wavelet－Decomposition and STFT Method in Continuous Deformation Observation Analysis［J］.Journal of Geodesy and Geodynamics.2011，31（5）：136－140.（in Chinese）

［7］姚家駿，楊立明，馮建剛.常用時(shí)頻分析方法在數(shù)字地震波特征量分析中的應(yīng)用［J］.西北地震學(xué)報(bào)，2011，33（2）：105－110.YAO Jia－jun，YANG Li－ming，F(xiàn)ENG Jian－gang，et al.Application of Common Time－Frequency Analysis Methods in Analyzing Characteristic Quantity of Digital Seismic Wave［J］.Northwestern Seismological Journal，2011，33（2）：105－110.（in Chinese）

［8］R G Stockwell，L Mansinha，R P Lowe.Localization of the Complex Spectrum：the S Transform［J］.Signal Processing，IEEE Transactions on，1996，44（4）：998－1001.

［9］周竹生，陳友良.含可變因子的廣義S變換及其時(shí)頻濾波［J］.煤田地質(zhì)與勘探，2011，39（6）：63－66.ZHOU Zhu－sheng，CHEN Yong－liang.Generalized S－transform with Variable－factor and its Time－frequency filtering［J］.Coal Geology ＆ Exploration，2011，39（6）：63－66.（in Chinese）

［10］Huang N E.Lntroduction to the Hilbert－Huang Transform and its Related Mathematical Problems［J］.Lnterdisciplinary Mathematics，2005，（5）：1－26.

［11］武安緒，吳培稚，蘭從欣，等.Hilbert－Huang變換與地震信號(hào)的時(shí)頻分析［J］.中國(guó)地震，2005，21（2）：207－215.WU An－xu，WU Pei－zhi，LAN Cong－xin，et al.Hilbert－Huang Transform and Time－Frequency Analysis of Seismic Signal［J］.Earthquake Research in China，2005，21（2）：207－215.（in Chinese）

［12］李大虎，梁明劍，黎小剛，等.2011年四川爐霍MS5.3地震加速度記錄的時(shí)頻分析與能量計(jì)算［J］.西北地震學(xué)報(bào)，2013，34（4）：335－341.LI Da－h(huán)u，LIANG Ming－jian，LI Xiao－gang，et al.Time－frequency Analysis and Energy Calculation for Acceleration Records of Luhuo MS5.3Earthquake in Sichuan Province in 2011［J］.Northwestern Seismological Journal，2013，34（4）：335－341.（in Chinese）