基于EMD的長周期地震動分量提取方法

譚潛，李英民，楊永斌

(1.中機中聯工程有限公司，重慶 400039；2. 重慶大學 土木工程學院,重慶 400045)

在實際地震記錄形成過程中，存在各種成分地震波的疊加(圖1)[1]，故地震記錄的主要特性由其主要成分特性來表達更具代表性。對地震動原始記錄的分析處理是以往研究地震動的主要方向，其忽視了地震動分量分析，而地震記錄中存在長周期分量的事實不能否認[2]。在長周期地震動界定中，徐龍軍等[3]以規準反應譜的離散性分析為基礎，指出地震動的反應譜特征與地震記錄的長周期分量是不可忽視的兩個方面。對地震記錄各類分量的提取和分析是確定地震動類型的可靠方法。地震動分量提取方法與信號分析方法息息相關。學者們在信號處理方法的基礎上提出了相關分量提取方法，Khanse等[4]利用傅里葉變換方法確定頻率，結合巴特沃斯濾波器對地震記錄進行分量分離，提取出長周期地震動脈沖分量，進而對相應的位移反應進行分析評估；陳紅等[5]采用分數階伽柏變換對地震信號進行時頻分析，其中一項重要的應用就是利用分數階伽柏變換來進行分數域譜分解，同時，也可以提取頻域中的不同尺度和不同方向上的特征，優勢明顯；Baker[6]以小波理論為基礎，提取地震記錄分量，從而對脈沖和非脈沖的地震記錄進行識別，具有較好的區分度，但是，不能對地震動分量的特性進一步的描述是此方法的不足；Mallat分解與重構是一種基于多尺度分析的小波變換分離分量的關鍵技術，趙國臣等[7]將其引入到地震記錄分量提取中，認為是一種有效的方法；Farid等[8]應用滑動平均濾波方法，結合截止頻率，將原地震記錄分解，提取脈沖型分量，剩下作為剩余分量。將提取的分量輸入結構，探索結構的彈性反應。商業軟件Seismo Signal中自帶數字濾波器，可以對地震動進行濾波，提取相應的分量，徐龍軍等[3]利用軟件的此功能提取脈沖分量，并對地震記錄反應譜及提取的分量反應譜的離散性進行分析和校驗。

在經驗模態分解(empirical mode decomposition， EMD)等信號處理與分量分離技術在地震動信號處理中應用的基礎上，采用合適評價長周期地震動的截止周期，實現了表征長周期地震動特性的長周期分量的提取，驗證了提取分量的合理性與有效性，提出可應用于長周期地震動識別及其規律特性研究的地震動分量提取方法。

圖1 長周期地震動的疊合Fig.1 Superposition of the long period ground motions

1 EMD理論

EMD作為一種新型自適應非平穩信號處理方法，不僅對線性平穩信號處理效果較好，而且對地震波等非線性、非平穩信號的分析處理有很好的適用性。本征模態函數(IMF)是EMD的前提。EMD方法認為信號都是由若干IMF組成，如果IMFs相互重組疊加，便可以重構為復合信號[9]。所以，EMD的重要目標就是為了得到IMFs，實現這個目標受到兩個條件的限制：首先，峰值點的個數與函數穿零數必須相等或者最多相差一個；其次，局部極大值構成的上包絡線和極小值構成的下包絡線的平均值為零，以IMFs為基礎可展開EMD，其過程如下：

1)基于給定的地震記錄信號x(t)，分別標記出其所有局部極大值xmax(t)和局部極小值xmin(t)；

2)將所有的局部極大值(xmax(t))點、局部極小值(xmin(t))點分別用3次樣條曲線連接起來，形成上、下包絡線；

3)對包絡線上的局部極大值、局部極小值取平均

(1)

4) 移除初始地震記錄信號中的平均趨勢

d(t)=x(t)-m(t)

(2)

5)移除平均趨勢的d(t)如果滿足上述兩個限定條件，令ci(t)=d(t)作為一個IMF，則可以算出殘量r(t)=x(t)-d(t)；反之，如d(t)不滿足兩個限制條件，則進入第下一步繼續篩選；

6)反復運行1)～5)步，讓d(t)滿足限制條件要求，最終得到多個IMF和殘量rn(t)，rn(t)為單調函數或者為常數

x(t)-c1(t)=r1(t)

r1(t)-c2(t)=r2(t)

…

rn-1(t)-cn(t)=rn(t)

(3)

反過來，可以通過n個IMF分量(ci(t))與殘量(rn(t))的線性組合得到原始信號

(4)

2 截止周期確定

數字信號處理領域常用數字濾波技術調制頻率分量，文獻[4]采用頻率1.67 Hz作為識別脈沖特征的界限。Rathje等[10-11]通過對地震記錄的頻譜周期參數對比后發現，加速度反應譜平均周期Tavg(式(5))在反映長周期地震動的低頻特性方面具有無可比擬的優勢，且采用頻譜參數Tavg能將地震記錄的高、低頻分量分布情況更準確地區分，低頻分量更容易被識別。杜東升等[12]研究表明，Tavg的變化與震級的變化具有一致性，Tavg是敏感性更好的頻譜參數，在評價長周期地震動時更具優勢。由Tavg計算公式可以知道其物理意義表示的是周期關于對應譜值平方的一個加權平均值，反映了地震動頻譜在整個計算頻域內的分布情況[11]。

0.05 s≤Ti≤10 s

(5)

式中：Ti為加速度譜對應的離散等間隔周期點；Sa(Ti)為周期點的加速度反應譜；PGA為地震動加速度時程對應的峰值。

為了更好地體現Tavg對地震動特性區分的實用性，以Tavg在0.4～5.1 s范圍內的10條汶川地震記錄為數據基礎，結合地震動的速度時程、歸一化傅里葉幅值譜和速度譜(圖2)，分析地震動頻譜特性隨Tavg變化而變化的規律：1)從速度時程看，地震動的速度時程曲線隨Tavg的增大而越來越稀疏，穿零點越來越少，表示其長周期成分不斷豐富，表明Tavg的變化能在一定程度上反映地震動周期的變化；2)從傅里葉譜(0.1～20 Hz)看,灰色背景的主頻區域隨Tavg的增大而向低頻方向平移，低頻成分占比逐漸增加，顯示出地震記錄主要頻譜成分隨Tavg的變化而變化的趨勢，以Tavg=2 s為界，1 Hz以下的低頻成分占比迅速增加，1 Hz以上的高頻部分占比迅速減小；3)從速度譜看，譜幅值相對較大的區段則隨著Tavg的增大向長周期段平移，同樣以Tavg=2 s為界，越大長周期成分越豐富，速度反應譜峰值2 s后基本不斷向右移動。可見，Tavg=2 s界限區分明顯，可作為長周期分量提取的頻譜指標。

圖2 地震記錄速度時程、歸一化傅里葉譜和速度譜伴隨Tavg的變化特征Fig.2 Variation characteristics of velocity time-history, normalized Fourier spectrum and velocity spectrum with vary Tavg

3 分量提取方法步驟

1)地震記錄的提取。以汶川地震中的長周期地震動chnua370505為例，其加速度、速度、位移時程曲線見圖3(a)，采用EMD方法的步驟對原始地震記錄(Original)加速度時程分解，提取地震記錄中的IMFs和殘余分量(R)，提取的IMFs以I1、I2、I3…I12表示，見圖3(b)，通過EMD的步驟可以知道，最終提取的IMFs和R均需符合EMD中的兩個限制條件。

圖3 chnua37050時程及IMF分量Fig.3 Time history and IMF components of chnua370505

2)地震動長、短周期分量的重構。基于地震動chnua370505的IMFs求解對應的Tavg值，如表1所示；由于地震記錄周期10 s以后不確定性的存在，導致在限定周期范圍時僅取值10以內。以Tavg=2 s為界，將提取的2 s及以上的IMFs與R疊加重構為長周期分量(long period components, LPC)，小于2 s的IMFs重新組合為短周期分量(short period components, SPC)。表1顯示：地震動的前4個IMFs的Tavg小于2 s，疊加組合為初始短周期分量，后面7個IMFs與R的Tavg大于2 s，重新組合為初始長周期分量，組合后的長、短周期分量時程見圖4。根據EMD中地震信號重構公式(4)可知，重構后的加速度LPC與SPC應符合式(6)、式(7)的范式，形成的長、短周期分量與IMFs疊加的誤差應在可接受范圍內，式中k表示Tavg≤2 s的IMFs最大分量的下標。進一步對加速度分量積分，得到速度分量(圖4)，為后續校正做準備。

(6)

(7)

圖4 疊加后的加速度與速度分量Fig.4 Acceleration and velocity components after superposition

表1 chnua370505的IMFs對應的Tavg值Table 1 Tavg values corresponding to IMFs of chnua370505 s

3)分量基線漂移的消除。由圖4的加速度與速度時程曲線可看出，重新組合的加速度分量時程無基線漂移，速度時程偏移明顯。與原始地震動進行基線校正同樣的目的，為了消除基線偏移對信號處理結果造成的干擾，保證分析結果的有效性，對重構后的加速度分量時程進行基線校正與濾波，形成校正后最終的加速度長周期分量與短周期分量，積分得到速度分量時程，消除基線漂移后的加速度、速度時程如圖5所示。可以看出，基線校正對基線漂移的控制效果明顯，速度分量時程的偏向大幅度減小，校正后的地震動分量能更好地代表實際地震動，可將地震動中的干擾因素最大程度地排除。

圖5 消除基線漂移的分量時程Fig.5 Component time history by eliminating baseline drift

4 長周期分量提取方法的有效性驗證

4.1 地震波分量提取

以Chi-Chi地震中的地震記錄ILA056-NS與TCU052-EW驗證本文提出的提取方法的有效性。基于地震動記錄的波形特點，王博等[13]將長周期地震動分為遠場長周期地震動與近斷層長周期地震動，從而可將ILA056-NS歸類于前者，而TCU052-EW歸類于后者；圖6展示了兩條地震記錄的加速度與速度時程與其分量。可以看出：在原始記錄中，基線漂移均未在加速度時程和速度時程中出現；觀察加速度、速度分量的時程，基線偏移在加速度分量中沒有出現，而速度分量偏移明顯，是對上述分量提取方法二次校正的必要性的進一步確認。

圖6 ILA056-NS 、TCU052-EW原地震記錄及分量加速度和速度時程Fig.6 Acceleration and velocity time history of original seismic records (ILA056-NS 、TCU052-EW)and their components

4.2 原地震記錄與其長周期分量相關性

驗證提取分量的代表性和有效性，可以分析地震記錄分量與原地震動的相關性。圖4與圖6顯示，未校正的地震記錄分量速度時程存在明顯基線偏移；對比圖5中校正后的分量速度時程，漂移消除效果較好。分析校正后的地震動分量速度時程與原地震動的相關性更可靠。從圖7可看出相關性趨勢：從地震記錄的波形判斷，兩條地震記錄校正后的LPC與原始地震動吻合較好，SPC與原地震動吻合較差，說明長周期地震動的LPC在很大程度上能夠代表長周期地震動，而SPC則不能代表，表明上述提取分量的代表性與有效性，同時，表2中的相關性系數計算值也更進一步說明基于速度的LPC與長周期地震動速度時程的相關性更好。

圖7 原地震記錄與重構校正后分量的相關性Fig.7 Correlation between original seismic records and their reconstructed and baseline drift corrected components

4.3 地震地面運動反應譜的離散性分析

分量提取方法的適用性已經通過相關性驗證做出說明，進一步可以通過地震記錄分解前后反應譜的變化來驗證。在設計反應譜準確估計的影響原因中，基于統計分析的地震反應譜的離散性不可忽視，且加速度反應譜長周期段的離散性更加突出[3]。圖8對比了兩條原始地震記錄與其長、短周期分量的反應譜的吻合程度，ILA056-NS反應譜與長周期分量反應譜在長周期段整體吻合較好，TCU052-EW吻合稍差，短周期分量與兩條原地震記錄在短周期段吻合，長周期段偏差顯著。此種差異，是由遠場長周期地震動與近斷層長周期地震動的長周期分量更豐富的分布特征決定的。徐龍軍等[3]以離散性系數k來評價地震反應譜的離散性，其計算公式見式(8)。借用參數k來評估LPC反應譜與原始地震記錄反應譜的離散性。離散系數k與離散性是正相關關系。基于兩條地震記錄與其分量反應譜，計算得到k，見表3。從表中k的值可判斷，長周期分量與原長周期地震動的離散性較小，短周期分量與原長周期地震動的離散性較大，同圖8結果一致。

表2 原地震記錄與分量速度時程的相關性系數Table 2 Correlation coefficients of velocity time history of original seismic records and their components

圖8 地震記錄(ILA056-NS、TCU052-EW)及其長、短周期分量反應譜Fig.8 Seismic records(ILA056-NS、TCU052-EW)and their long , short periodic components response spectra

(8)

式中：S1T為原地震動在周期T時的反應譜值；S2T為各分量在周期T時的反應譜值；ΔT為反應譜的計算周期間距，取為0.02；Ttol為反應譜的計算最大周期，取15 s。

表3 長、短周期分量與原地震記錄反應譜的離散性系數Table 3 Discretization coefficients of response spectra of original ground motions and their long, short period components

5 不同分量提取方法的結果對比

Baker[6]對脈沖型地震動進行分量分離時采用了基于小波分析的分量提取方法，小波分析過程中，母波選擇多貝西四階小波，將地震記錄分解為15個分量，并計算每個分量的Tavg值,其值大于2.0 s的重構為長周期分量，小于2 s的組合為短周期分量；Farid Ghahari等[8]采用滑動平均濾波方法提取脈沖型分量。其步驟：1)用短時傅里葉變換計算速度時程的卓越周期Tp；2)采用式(9)結合卓越周期Tp得到參數m的值，式(9)中參數α采用文獻[8]中的經驗取值0.25，dt為地震記錄采樣點的時間間隔；3)以m值代入式(10)計算截止頻率fc，fc取倒數得到截止周期Tc(式(11))；4)利用Tc，結合滑動平均濾波來識別長、短周期分量。將上述兩種方法提取的分量反應譜與以EMD方法提取的分量轉換反應譜一起與原地震動的反應譜進行吻合性對比，由圖9可以看出，基于EMD的長周期分量提取方法提取的長周期分量反應譜在長周期段與原地震動反應譜整體上吻合較好；與基于小波分析的方法相比，在ILA056-NS中，長周期段具有相近的效果，在TCU052-EW中，本文的方法提取效果更好；而與基于滑動平均濾波的方法比，本文方法提取的長周期分量均能更好地預測地震記錄和反應譜的趨勢。

(9)

(10)

(11)

圖9 EMD、MA、Wavelet分量提取結果對比Fig.9 Comparison of component extraction results by EMD、MA、Wavelet

6 結論

研究了基于EMD的長周期地震動分量提取方法，對方法的有效性進行了分析，并與以往的分量提取方法進行對比，主要結論包括：

1)Tavg大于2 s的長周期成分在長周期地震動中起著至關重要的作用。由于對于大多數長周期地震記錄，周期高于2 s的速度譜振幅通常比周期低于2 s的速度譜振幅更具優勢，2 s以下基本不再出現峰值點。歸一化傅里葉譜的頻率成分偏向于低頻，通過地震動特性分析確定了分量分離的周期參數(Tavg)及取值，該結果被建議作為長周期地震動識別程序中的關鍵參數。

2)采用經驗模態分解方法，將地震動波分解成IMFs；結合截止周期Tavg=2 s，將相應的IMFs分別重組為長周期分量、短周期分量并校正，進而建立地震動長周期分量提取方法，并提出了分量提取實現的程序。

3)基于原地震動與長、短周期分量的相關性，地震動反應譜與長、短周期分量的反應譜的離散性對分量提取方法的有效性進行了驗證，并將基于EMD的分量提取方法與基于小波、滑動平均濾波的方法進行了比較，驗證了方法的合理性與有效性。提取的長周期分量是長周期地震動識別模型化的數據基礎與參數指標來源。