2011年四川爐霍MS5．3地震加速度記錄的時頻分析與能量計算①

李大虎，梁明劍，黎小剛，顧勤平

（1．四川省地震局，四川 成都 610041；2．中國地震局地球物理研究所，北京 100081；3．江蘇省地震局，江蘇 南京 210014）

0 引言

2011年4月10日17時02分，四川省甘孜藏族自治州爐霍縣境內(nèi)發(fā)生MS5．3地震（東經(jīng)100．9°，北緯31．3°），震中位于鮮水河斷裂爐霍段上。四川省數(shù)字強震動觀測臺網(wǎng)在震區(qū)附近的8個臺站獲取了主震的三分向加速度記錄，其中在震中距40km范圍內(nèi)獲得了3條140cm／s2以上的記錄，最大水平峰值加速度記錄為爐霍地辦臺記錄到的EW向分量403．5cm／s2。臺站分布位置如圖1所示。利用P波初動方向資料和CAP方法求解出的震源機制解反映此次地震為純走滑型的斷層錯動性質(zhì)。

為了提取到地震發(fā)生后更多關于強震記錄的時頻特性和動力特征，為地震波的傳播規(guī)律及抗震設計等方面的研究奠定基礎［1］，長期以來地震工程界對強震加速度記錄最常用的處理方法就是傅里葉變換和小波分析，抑或S變換、CWD分布、ZAM 分布等時頻分析方法［2］，而這些方法大都是建立在信號為穩(wěn)態(tài)的基礎上，并且由于Fourier變換無法進行時頻局部化分析和小波變換存在基函數(shù)選取不定等問題［3］。因此，近年來很多專家與學者采用了具備精確分解優(yōu)勢和非線性動態(tài)數(shù)據(jù)刻劃能力的希爾伯特黃變換（HHT）來對強震記錄進行時頻分析［4－5］。但是由于經(jīng)驗模態(tài)分解（EMD）過程中存在的模態(tài)混疊問題，這不但會使得單獨的固有模態(tài)函數(shù)（IMF）分量缺乏明確的物理意義，還會造成EMD算法運行過程中的不穩(wěn)定。為了抑制模態(tài)混疊效應同時提取到關于地震動更多的時頻特性及地震波的動力學特征，本文采用一種基于聚類經(jīng)驗模態(tài)分解（EEMD）的HHT計算方法，對爐霍地震中獲取的不同震中距的臺站記錄（爐霍專業(yè)臺、爐霍地辦臺以及甘孜地辦臺）進行時頻分析和能量計算，以得到不同臺站獲取的加速度記錄信號的邊際譜、相位譜、Hilbert時頻幅值譜及能量譜，量化提取其平均周期、中心頻率和Hilbert能量在不同時頻段上的分布等特性，進而與快速傅里葉變換（FFT）、小波變換（WT）進行對比研究，并對所選取的強震臺站中各頻帶的能量進行定量計算，最后從瞬時能量譜的角度分析其與逆—走滑型汶川8．0級地震的強震記錄的不同，以期從中獲得一些新的認識。

圖1 區(qū)域地震構(gòu)造與臺站分布圖Fig．1 Distribution of regional seismic tectonics and seismic stations．

1 EEMD方法原理與譜分析

本文之所以選擇基于聚類經(jīng)驗模態(tài)分解（EEMD）的HHT進行加速度記錄的時頻分析在于其不但能夠有效抑制以往EMD分解過程中出現(xiàn)的模態(tài)混疊問題，而且使得分解出的IMF分量的物理意義更加清晰。該方法主要由EEMD分解和Hilbert變換兩部分組成。

1．1 EEMD分解原理及算法流程

本文所采用聚類經(jīng)驗模態(tài)分解原理為：當附加的白噪聲均勻分布在整個時頻空間時，該時頻空間就由濾波器組分割成的不同尺度成分組成，那么，對于每次EMD分解，添加的白噪聲在整個時頻空間是均勻分布的，信號的不同頻率尺度被自動投影到由白噪聲所建立的均勻時頻空間的相應頻率尺度上。由于每次EMD分解添加不同的白噪聲，噪聲之間不相關，因此對所有EMD分解的相應IMF求整體平均后，人為添加的噪聲被抵消掉，全體的均值最后將會被認為是真正的結(jié)果。唯一持久穩(wěn)固的部分是信號本身，所加入的多次測試是為了消除附加的噪聲［6－7］。對于將EEMD分解應用于在提取強震動記錄時頻特性方面，筆者已對汶川地震不同斷層距的強震記錄做過相應的對比研究［8］。其算法流程如下：

（1）初始化EMD執(zhí)行總次數(shù)M，白噪聲的幅值系數(shù)k，m＝1；

（2）執(zhí)行第m次EMD分解試驗；

①在輸入信號x（t）中加入一組隨機高斯白噪聲序列nm（t），得到加噪的待處理信號xm（t）：xm（t）＝x（t）＋k nm（t）；

②用EMD分解xm（t），得I個IMFcjm（j＝1，…，I），cjm表示第m次計算分解出的第j個IMF；

③若m＜M，返回步驟（2），m＝m＋1；

（3）對M次試驗的每個IMF計算均值 （j＝1，2，…I；m＝1，2，…M）；

（4）輸出 作為EEMD分解得到的第j個IMF，j＝1，2，…I。

1．2 Hilbert變換與譜分析

對地震動加速度信號的各IMFci（t）做Hilbert變換［9－10］，得到

上面兩式明確表達了瞬時振幅和瞬時相位，很好地反映了信號的瞬時特征。

在此基礎上再定義瞬時頻率為

于是可以得到

進而得到

這就稱作Hilbert譜：

再定義Hilbert邊際譜：

另外作為Hilbert邊際譜的附加結(jié)果，可以定義Hilbert瞬時能量譜：

瞬時能量譜提供了信號能量隨時間的變換情況。

如果振幅的平方對時間積分，可以得到Hilbert能量譜：

2 加速度記錄的時頻分析

2．1 EEMD分解結(jié)果

選取不同震中距的臺站獲取的加速度記錄進行處理和計算，包括爐霍專業(yè)臺（震中距14．2km）NS向分量和爐霍地辦臺（震中距24．0km）NS向分量。對爐霍專業(yè)臺NS向分量進行EEMD分解結(jié)果的分量及其功率譜分別如圖2所示，從圖2中可以看出，原信號已被分解成8個IMF分量（C1～C8）和1個殘余項R（res）。并且隨著分解的進行，所得到的IMF分量頻率逐漸降低，直到分解出頻率已經(jīng)很低的最后一個趨勢項R。其中各個IMF分量包含了不同的尺度特征。

關于每個IMF分量具體的時頻信息，經(jīng)計算后列于表1。

圖2 基于EEMD分解得到的爐霍專業(yè)臺各IMF分量Fig．2 The IMFcomponents from decomposition of EEMD recorded at Luhuo professional station．

表1 IMF分量統(tǒng)計特征值

除了以上基于EEMD分解得到的各IMF分量分析時頻特性外，本文又對部分分量進行了功率譜分析計算，由于篇幅所限，在此僅選取了3個IMF（IMF2、IMF4和IMF6）的功率譜圖（圖3）進行說明。圖4為震中距14．2km的爐霍專業(yè)臺所獲取的加速度記錄NS向分量的原始記錄，通過FFT分析可以看出，原始信號的頻譜較豐富，且大部分分布在0～20Hz以內(nèi)。

圖3 三個IMF分量的功率譜圖Fig．3 The power spectrum diagrams of three IMFcomponents．

圖4 爐霍專業(yè)臺的加速度時程記錄Fig．4 Acceleration schedule records of Luhuo professional station．

2．2 Hilbert變換與譜分析

將EEMD分解后的IMF各分量經(jīng)Hilbert變換和文中1．2所列公式計算后，得到了圖5中的Hilbert邊際譜。與FFT譜對比可以發(fā)現(xiàn)，在低頻處FFT譜會低估地震動的幅值；隨著頻率的增加FFT譜又會放大地震動的幅值。同樣對其他兩個不同通道（EW向和UD向）以及所選取的不同震中距的臺站（爐霍和甘孜地辦臺）獲取的強震記錄進行傅里葉譜與邊際譜對比也可以發(fā)現(xiàn)這個現(xiàn)象。

3 能量計算與分析

3．1 能量計算

將震中距為14．2km的爐霍專業(yè)臺獲取到的加速度記錄NS向分量和震中距為24．0km的爐霍地辦臺獲取到的加速度記錄NS向分量分別記為DATA1NS和DATA2NS，從Hilbert能量譜（圖6）中可以看出此次地震信號能量在頻率上的集中程度。DATA1NS的能量分布在0～20Hz的主振頻域范圍以內(nèi)，DATA2NS的能量分布在0～10Hz以內(nèi)，其中主頻帶之內(nèi)又包括多個子振頻帶，具體計算量值列舉在表2和表3之中。

圖5 邊際譜與傅里葉譜對比Fig．5 The comparison between Marginal spectrum and Fourier spectral．

表2 DATA1NS信號能量集中區(qū)各頻段對應的能量

表3 DATA2NS信號能量集中區(qū)各頻段對應的能量

圖6 不同臺站記錄的NS向分量能量譜Fig．6 The energy spectrums of NS component from different station records．

從表2和表3中還可以看出，從不同震中距臺站所接收到的加速度記錄提取到的能量最大的優(yōu)勢頻段是不同的，DATA1NS的主要能量集中1～2 Hz，而DATA2NS主要能量集中4～5Hz。需要特別指出的是，經(jīng)過HHT得到的地震動能量譜圖并不是真正的能量值，它的量綱也不是國際單位制的焦耳，或工程單位制中的公斤力米或CGS制中的爾格，究其原因，是因為Hilbert變換本質(zhì)上是一個寬頻帶全通濾波，當?shù)卣饎蛹铀俣葐挝蝗al時，IE（t）的單位就是gal×gal×rad／s。所以 Hilbert能量譜只是間接地反映了輸入能量隨頻率的變化情況，至于能量隨著時間變化趨勢則用Hilbert瞬時能量譜來表示，后面的章節(jié)將會介紹。

3．2 基于能量概念的崩科結(jié)構(gòu)破壞機理初步分析

與國內(nèi)其他建筑物結(jié)構(gòu)形式不同的是，由于爐霍縣位于高原地區(qū)，建筑物大多屬于藏式民居，當?shù)胤Q之為“崩科”。“崩科”結(jié)構(gòu)的房屋使用了大量的木材，并用夯土墻或塊石墻打圍，內(nèi)部采用穿斗木架通柱、雙梁雙檁雙掛條，地嵌鎖腳井字形屋架（圖7）。雖然這種結(jié)構(gòu)的房屋具有良好的抗震性能，但在本次MS5．3地震中，各鄉(xiāng)鎮(zhèn)行政村的房屋還是出現(xiàn)了不同程度的破壞，具體表現(xiàn)為土墻或石墻出現(xiàn)裂縫、傾斜或局部垮塌，主體結(jié)構(gòu)破壞并不嚴重。

圖7 宜木鄉(xiāng)“崩科”結(jié)構(gòu)房屋Fig．7 “Bengke”building structure in Yimu village．

由于基于能量概念的分析能夠較好地反映振動強度、頻譜特性，特別是持時對建筑物的綜合影響，因此對它的研究成了研究抗震設計的一個重要發(fā)展方向［11－12］。對工程的抗震設計而言，我們在關心地震的總能量的基礎上，更應該關心的是地震動峰值能量的最大值，它表征了地震動對結(jié)構(gòu)破壞能力的大小，因此地震動峰值能量的最大值能夠較為全面地反映地震動三要素（幅值、頻譜和持時）的綜合作用效果。從表2和表3中可以看出，不同震中距的臺站所記錄到的地震動峰值能量的最大值是不同的，而且能量分布的優(yōu)勢頻段也是不同的，DATA1NS和DATA2NS的最大能量分別集中1～2 Hz和4～5Hz之間，這對我們通常的幅值越大，持時越長，累積能量和峰值能量越大的認識做了補充，而且對于不同頻帶而言，相鄰頻段上如果地震波都具有較高的能量，那么，它對于結(jié)構(gòu)的沖擊通常就會越大，結(jié)構(gòu)發(fā)生破壞的可能性也就越大。崩科結(jié)構(gòu)的房屋雖然有別于一般意義上的土木或磚木結(jié)構(gòu)的房屋，但也同樣適用于這一破壞機理。并且，從表3中所列舉的爐霍地辦臺加速度記錄的能量分布還可以發(fā)現(xiàn)，相鄰頻段都具有較高的能量造成的破壞會很大，這與震害調(diào)查中所確定該地（爐霍縣新都鎮(zhèn)）為宏觀震中相一致。所以，基于能量概念結(jié)合HHT計算，對建筑物的震害進行嘗試研究具有探索意義。

4 與汶川地震加速度記錄的對比

由于汶川MS8．0地震是龍門山構(gòu)造帶逆沖－右旋錯動的結(jié)果，而此次爐霍地震則是純走滑型的斷層錯動，這也導致了其能量的釋放過程中出現(xiàn)的差異。因此本文選取了汶川地震中距離斷層在30 km內(nèi)的德陽白馬臺所獲取的強震記錄NS向分量（簡記為BMNS）做了基于EEMD分解的HHT分析得到該分量的瞬時能量譜圖（圖9），圖8（a）、（b）、（c）分別為爐霍專業(yè)臺、爐霍地辦臺以及甘孜地辦臺（震中距93．7km）獲取到的加速度記錄NS向分量的瞬時能量譜圖。對比這三張圖可以發(fā)現(xiàn)，圖9所示的德陽白馬臺記錄NS向分量瞬時能量譜圖反映了該次地震能量的釋放是分段的，這與張勇等人［13］通過對震源時間函數(shù)分析得出汶川地震整個時間過程有5個主要的能量釋放階段基本一致；而從圖8則反映出爐霍地震的能量釋放大體上經(jīng)歷了一至兩次的釋放過程，主要能量釋放集中出現(xiàn)在5s左右。因此，Hilbert瞬時能量譜不但能夠清晰地表達能量隨時間的變化情況，而且還可以發(fā)現(xiàn)能量峰值出現(xiàn)的時刻，它在表達地震動信號的時頻特性方面有著較好的用途。圖10為爐霍專業(yè)臺獲取的SN向加速度時程信號進行基于sym4小波函數(shù)變換后得到的小波譜圖，從圖中雖然可以看出能量大部分集中在高尺度的低頻區(qū)，但它也只是定性描述，具體的時頻對應關系無法定量提取，所以小波譜無法表現(xiàn)出很好的局部化效果，更無法像Hilbert瞬時能量譜（圖9）和Hilbert邊際譜（圖5）一樣較為精確地表達 信號能量的時頻分布特性。

圖8 不同臺站記錄NS向分量瞬時能量譜Fig．8 Instantaneous energy spectrums of NS component from different station records．

圖9 德陽白馬臺記錄NS向分量瞬時能量譜Fig．9 Instantaneous energy spectrum of NS component from Baima station records in Deyang．

圖10 BMNS基于sym4小波基的小波譜Fig．10 Wavelet spectrum of MANS based on sym4 wavelet－basis．

5 結(jié)論與建議

（1）本文選取了2011年四川爐霍地震中不同震中距的臺站所獲取的三分量加速度時程記錄，采用基于EEMD分解的HHT進行了能量計算和時頻分析，不但有效抑制了以往EMD分解過程中所出現(xiàn)的模態(tài)混疊問題，而且還較好地提取了記錄的時頻特性和能量集中分布的時頻段，從而證明了該方法在處理非平穩(wěn)、非線性地震信號中的有效性和實用性，在地震工程領域有著較好的應用前景。然而，與Fourier變換和小波分析不同的是，HHT作為一種經(jīng)驗性的自適應方法目前仍存在一些不完善的地方，例如樣條差值的量化及邊界效應的處理等，在今后還需不斷完善其相關的數(shù)學理論，并研究分解過程中的快速算法，爭取將EEMD算法和一些可視化編程工具（Visual C＋＋、C＃或．NET）結(jié)合起來，開發(fā)出基于EEMD分解和頻譜分析的應用軟件，將會為以后的地震數(shù)據(jù)處理工作帶來很大的方便。

（2）針對藏區(qū)特有的崩科結(jié)構(gòu)的房屋，本文只是基于能量概念從相鄰頻段的能量集中分布會對建筑物造成的破壞情況進行了初步的定性－半定量的分析，如果需要詳細研究該崩科結(jié)構(gòu)詳細的振動特性和破壞情形，還必須結(jié)合各臺站具體的場地條件等資料，對加速度記錄數(shù)據(jù)進行基于HHT的瞬時輸入能量定量計算。除此之外，還需采用系統(tǒng)識別理論和結(jié)構(gòu)動力學等相關知識并量化分析，識別該種結(jié)構(gòu)的基本參數(shù)（質(zhì)量、阻尼、剛度、反應譜等），判別在結(jié)構(gòu)遭遇地震后房屋受到損傷后的阻尼、剛度和自振頻率變化情況，為抗震設計提供堅實的依據(jù)。

致謝：在本文的撰寫過程中，得到了成都理工大學地球探測與信息技術教育部重點實驗室李才明教授的悉心指導和幫助；在與李軍博士的探討過程中也聽取了很多建議；審稿專家對本文提出了寶貴的修改意見；臺站的加速度記錄數(shù)據(jù)由四川省地震局臺站管理中心的賴敏老師提供。在此一并向他們表示最衷心的感謝！

［1］ 盧大偉，姚凱，閔祥儀，等．蘭州觀象臺存放臺陣汶川MS8．0強震動記錄與分析［J］．西北地震學報，2011，33（2）：171－176．

［2］ 姚家駿，楊立明，馮建剛．常用時頻分析方法在數(shù)字地震波特征量分析中的應用［J］．西北地震學報，2011，33（2）：105－110．

［3］ 李大虎，何強，李才明，等．MATLAB與C＋＋混合編程實現(xiàn)航磁異常提取的小波分析方法研究［J］．地震研究，2011，25（2）：98－105．

［4］ 武安緒，吳培稚，蘭從欣，等．Hilbert－Huang變換與地震信號的時頻分析［J］．中國地震，2005，21（2）：207－215．

［5］ 公茂盛，謝禮立，連海寧，等．基于HHT的結(jié)構(gòu)強震記錄分析研究［J］．地震工程與工程振動，2007，27（6）：24－29．

［6］ Z Wu，N E Huang．A study of the characteristics of white noise using the empirical mode decomposition method［J］．Proc．R．Soc．1and．A，2004，460（2046）：1597－1611．

［7］ Wu Z H，Huang N E．Ensemble empirical mode decomposition：a noise assisted data analysis method［J］．Advances in A－daptive Data Analysis，2009，1：1－41．

［8］ 李大虎，賴敏，何強，等 ．基于聚類經(jīng)驗模態(tài)分解（EEMD）的汶川MS8．0強震動記錄時頻特性分析［J］．地震學報，2012，34（3）：350－362．

［9］ Rilling G，F(xiàn)landrin P，Gonc P，et al．On empirical mode decomposition and its algorithms［A］∥IEEE－EURASIP Workshop on Nonlinear Signal and Image Processing，NSIP－03［C］．Grado，Italy：IEEE，2003：1－5．

［10］ Flandrin P，Rilling G，et al．Empirical mode decomposition as a filter bank［J］．IEEE Signal Process，2004，11：112－114．

［11］ Akiyama H．Earthquake－resistant design for building［M］．Tokyo，Japan：University of Tokyo Press，1985．

［12］ Kamamnra H，Galambos，Theodore V．Earthquake loads structural for reliability［J］．J．STRUCT．ENG－ASCE．，1989，115（6）：1446－1462．

［13］ 張勇，馮萬鵬，許力生，等．2008年汶川大地震的時空破裂過程［J］．中國科學（D輯：地球科學），2008，38（10）：1186－1194．