楊彩紅 梁建文 張郁山

1）中國地震災害防御中心，北京 100029

2）天津大學土木工程系，天津 300072

引言

1998年，Huang等（1998）提出了處理非平穩(wěn)信號的HHT方法（Hilbert-Huang Transform，簡稱HHT）。該方法包括2個步驟：①任意信號首先經過經驗模態(tài)分解方法（Empirical Mode Decomposition，簡稱 EMD）被分解為一系列固有模態(tài)函數（Intrinsic Mode Function，簡稱IMF）；②對每個IMF進行Hilbert譜分析（Hilbert Spectral Analysis，簡稱HSA），得到相應的Hilbert譜，最后將所有IMF的Hilbert譜相加得到原始信號的Hilbert譜。Hilbert譜將原始時域信號表示在聯合的時間-頻率域上，是原始信號的一種時頻表示。HHT方法自提出后，已經在地震動的地震學與工程學特性研究（Huang等，2001；Zhang等，2003；Loh等，2001）與結構系統(tǒng)識別領域中得到了應用（Yang等，2003a；2003b；2003c；Salvino，2000；Pines等，2001；Zhang等，2005；石春香等，2005），而且一些學者針對HHT方法自身存在的問題也展開了研究工作（Huang等，1999；鄧擁軍等，2001；張郁山等，2003；黃大吉等，2003；羅奇峰等，2003）。

針對傳統(tǒng)的EMD方法分解長周期信號所面臨的問題，本文提出了基于折線包絡的經驗模態(tài)分解方法，通過與已有EMD方法的比較論證了該方法的分解效率，并利用數值算例驗證了該方法所得結果的意義。

1 基本方法

1.1 基于樣條包絡的“篩分”

在HHT方法中，IMF是滿足如下2個條件的信號：①在信號的整個持續(xù)時間內，零交點（波形曲線與零軸的交點）的數目與極值點的數目必須相等或至多相差一個；②在任意時刻，由極大值點定義的上包絡線和由極小值點定義的下包絡線之間的平均值為零。利用EMD方法，任意信號可以被分解為若干IMF分量。Huang等（1998）詳細描述了EMD方法分解信號的具體過程，此處不再贅述。EMD對信號分解的核心運算是其“篩分”處理（sifting process）。傳統(tǒng)的“篩分”過程處理信號x(t)的具體運算如下：

1）將x(t)波形中所有局部極大值點和局部極小值點識別出來；

2）用三次樣條曲線將所有局部極大值點連接起來構成原始波形的上包絡線u()e t~ ，同樣再用三次樣條曲線將所有局部極小值點連接起來構成原始波形的下包絡線d()e t~ ，上下包絡線應將原始波形包在中間；

上述過程即稱為一次“篩分”。本文的研究既是針對該篩分過程進行的改進，此處將該“篩分”處理定義為基于樣條包絡的“篩分”，記作

圖1 原始信號及其三次樣條包絡與折線包絡Fig. 1 Original signal and its cubic-spline envelopes and folding-line envelopes

圖2 原始信號及其兩種平均包絡Fig. 2 Original signal and its two kinds of mean envelopes

圖3 兩種篩分方法所得到的結果Fig. 3 Results of two sifting methods

傳統(tǒng)的EMD方法利用上述“篩分”處理分解長周期信號時，由于信號相鄰極值點之間的距離較長，連接信號極值點的三次樣條包絡線在信號的相鄰極值點之間將會出現較大的“抖動”，這將導致信號的三次樣條上、下包絡線不能很好地包絡原始信號，信號的平均樣條包絡線也不能很好地分割原始信號，從而會影響EMD方法的分解效率。考慮信號x(t)，其波形如圖1（a）中實線所示，基于樣條包絡EMD方法在對其進行“篩分”處理時，首先將其所有極大值點和所有極小值點分別用三次樣條曲線連接起來形成三次樣條上、下包絡線，分別如圖1（a）中的虛線與點線所示。從中可以看出，由于x(t)的兩個相鄰極大值點A、B和兩個相鄰極小值點C、D相距較遠，三次樣條上包絡線在A點與B點之間出現了比較大的“抖動”，同樣，三次樣條下包絡線在C點與D點之間也出現了比較大的“抖動”。這就導致了在C點與B點之間，三次樣條下包絡線反而位于上包絡線之上，如圖所示。此外，盡管下包絡線在三個相鄰極小值點E、F與G之間“抖動”不大，但在t=50之后依然出現了下包絡線位于上包絡線之上的現象。這種由于三次樣條曲線的“抖動”而導致的下包絡線在某些區(qū)間位于上包絡線之上的現象必然會使得三次樣條上、下包絡線不能很好地包絡原始信號，從而會影響EMD方法的分解效率?；谏鲜隹紤]，本文提出一種基于折線包絡的EMD方法。

1.2 基于折線包絡的“篩分”

基于折線包絡的EMD方法依然利用“篩分”處理來分解原始信號。針對原始的基于樣條包絡的EMD方法所面臨的問題，基于折線包絡的EMD方法對信號x(t)的“篩分”處理過程如下：

1）將原始信號x(t)的所有極大值點用直線段連接起來形成原始信號的折線上包絡e?u( t)；將原始信號的所有極小值點用直線段連接起來形成原始信號的折線下包絡e?d( t)；

3）將折線平均包絡m?0( t)的所有折點用三次樣條曲線連接起來即形成原始信號的均值包絡 m?( t)；

4）原始信號x(t)減去上述均值包絡 m?( t)即得到一個新的信號y(t)。

上述過程即為基于折線包絡的EMD方法的一次“篩分”，此處稱之為基于折線包絡的“篩分”，該過程可用符號?S表示，即：

下面比較這兩種“篩分”的處理過程?？紤]到圖1所示信號，在傳統(tǒng)的基于樣條包絡的“篩分”處理中：首先將原始信號x(t)的所有極大值點和所有極小值點分別用三次樣條曲線連接起來形成上、下包絡線，如圖1（a）中虛線與點線所示；這兩條包絡的均值即為基于樣條包絡的平均包絡 ()m t~ ，如圖2（a）中點線所示；原始信號x(t)減去此平均包絡0()m t~ 即得到一個新的信號 ()y t~ ，如圖3（a）所示。

基于折線包絡的“篩分”處理，首先將原始信號 x(t)的所有極大值點和所有極小值點分別用直線段連接起來形成折線上、下包絡，它們的均值即為原始信號的折線平均包絡，如圖1（b）所示；該折線平均包絡的折點所在位置與原始信號極值點所在位置相同，如圖所示。將該折線平均包絡的所有折點用三次樣條曲線連接起來即形成基于折線包絡的平均包絡如圖2（b）所示；相比基于樣條包絡的平均包絡，基于折線包絡的平均包絡 m~ ( t)更好地分割了原始信號的波形，如圖2所示。將原始信號x(t)減去此平均包絡 ()m t~ 即得到一個新的信號?()y t，如圖3（b）所示。相比圖3（a）所示基于樣條包絡“篩分”所得到的 ()y t~ ，?()y t的波形更接近本征模態(tài)函數（IMF）的要求：①相比 ()y t~ ，?()y t的波形關于零軸更加對稱；②?()y t的所有極大值點均位于零軸以上，所有極小值點均位于零軸以下，而在 ()y t~ 的波形中依然存在著位于零軸以下的極大值點，如圖3（a）所標記的A、B兩點。

通過以上比較可以看出，基于折線包絡的平均包絡線更好地分割了原始信號，利用它進行一次“篩分”處理所得結果的性態(tài)更接近IMF的要求，因此相比基于樣條包絡的“篩分”方法，基于折線包絡的“篩分”方法具有更高的分解效率。

1.3 基于折線包絡的EMD方法

由于信號的波形中存在著隆起及其他高階導數不連續(xù)的奇點，上述一次“篩分”并不能得到滿足要求的IMF分量，通常需要多次“篩分”過程進行處理，最終所得信號才能滿足IMF分量的要求。對原始信號進行k次“篩分”處理可以記為：

仿照傳統(tǒng)EMD方法的處理過程，本文提出的基于折線包絡EMD方法的處理流程如下：

1）記對原始信號x(t)進行k1次與k1-1篩分所得結果為若它們之間的標準差：

滿足Δs＜0.3，則為信號x(t)的第一個IMF分量c1(t)：

2）令

若對 r1( t)篩分k2次后的結果與k2-1次的結果之間的標準差（按照公式（4）計算）滿足Δs＜0.3，則為原始信號的第二個IMF分量：

3）將上述分解過程重復，即可得到原始信號的一系列IMF分量：

4）如果最后得到的IMF分量cn(t)或殘余信號rn(t)的值非常低，低于預先設定好的值；或者最后的殘余信號rn(t)為時間的單調函數，在其波形中不存在極值點，那么上述分解過程結束。

圖4 含高頻正弦噪聲的衰減余弦調頻波Fig. 4 Attenuated frequency-modulated cosine wave containing sinusoidal noise

圖5 基于折線包絡的EMD方法分解圖4所示信號所得到的主要IMF分量與殘余分量Fig. 5 Main IMF components and residue obtained by decomposing the signal shown in Fig. 4 using the folding-line-envelope-based EMD method

2 算例

2.1 具有明確解析式的信號

考慮如圖4所示含噪聲的衰減余弦調頻波x(t)=x1(t)+x2(t)，該信號由幅值衰減的余弦調頻波

與高頻正弦噪聲

共同合成。其中信號x1(t)的瞬時頻率fi1(t)定義為其相位的導數

可以看出，其瞬時頻率隨時間作簡諧波動，其幅值隨時間呈指數衰減，衰減因子為0.2。應用基于折線包絡的EMD方法對信號x(t)進行分解將會得到12個IMF分量，其中前兩個分量c1(t)和c2(t)幅值明顯高于其余分量，具有明確的數學意義；其余分量的幅值較小，它們的出現主要是由數值誤差引起的，它們的和為r2(t)。c1(t)、c2(t)與r2(t)的波形如圖5所示，可以看出：r2(t)的幅值在初始時刻較大，但仍未超過0.1，它主要是由于EMD方法中數據邊界誤差引起的；在其余大部分時段內，r2(t)的幅值在一個非常小的范圍內波動。IMF分量c1(t)與 c2(t)分別對應于原始信號的合成分量x2(t)與x1(t)，它們之間的比較在圖6中給出，從中可以看出基于折線包絡的EMD方法分解所得到的IMF分量與原始信號合成分量的符合程度是非常好的。其中，c1(t)與x2(t)之間的比較是在區(qū)間[1.0,2.0]中給出的，其他區(qū)間二者的吻合程度與該區(qū)間類似。

圖6 基于折線包絡EMD方法分解x(t)所得主要IMF分量與x(t)的合成分量之間的比較Fig. 6 Comparison between the main IMF components obtained by using folding-line-envelope-based EMD decompose x(t) and the synthesizing components of x(t)

對圖5所示IMF分量進行Hilbert譜分析之后，可得到原始信號的Hilbert譜，如圖7所示。其中，頻率保持在15.0Hz的高頻分量為原始信號中的高頻噪聲部分x2(t)，其譜值非常低；頻率在1.0Hz波動的分量為原始信號中的余弦調頻波，其譜值隨時間衰減。針對該圖有兩點需要說明：①由于離散Hilbert變換與數值微分所帶來的誤差，低頻分量瞬時頻率只是大體上遵循式（11）所示的余弦變化；②由于在進行Hilbert譜分析時沿時間軸進行了平滑處理，以及數值計算所導致的瞬時頻率微小的波動，使得這兩個分量在其瞬時頻率附近均出現了較窄的頻帶。

圖7 原始信號的Hilbert譜Fig. 7 Hilbert spectrum of original signal

2.2 地震動加速度記錄

圖8 給出的波形為1940年美國El Centro強震加速度記錄。應用基于折線包絡EMD方法分解該記錄所得的IMF分量及每個IMF分量所對應的Fourier幅值譜，如圖9所示。從圖9（a）、圖9（c）中可以看出：①所有分量的波形都對稱于零軸；②每個分量波形的所有極大值點均位于零軸以上，所有極小值點均位于零軸之下。因此該方法分解El Centro波所得到的分量滿足IMF的要求，此外每個IMF分量的波形還具有類似于窄帶波包的特性。從圖9所示IMF分量的波形與Fourier幅值譜中可以看出：每個IMF分量均代表著一種振動模態(tài)，它們的幅值與頻譜含量各不相同，自第一個分量起，IMF分量的高頻成分逐步減少。

圖8 El Centro地震加速度記錄Fig. 8 Acceleration records obtained at El Centro station

對圖9所示IMF分量進行Hilbert譜分析，可得到原始El Centro波的Hilbert譜，如圖10所示。從中可以看出，地震動的Hilbert譜揭示了原始地震動的能量在時間-頻率域上的分布特征，與時間-頻率反應譜有相同的作用，對于研究地震動的工程特性具有重要意義（羅奇峰等，2003）。而且，相比小波譜，與原始EMD方法所得結果相似，本文提出基于折線包絡的EMD所得地震動Hilbert譜同樣具有更高的時頻分辨率（Huang等，1998）。

圖9 基于折線包絡的EMD方法分解El Centro波所得到的IMF分量及其Fourier幅值譜Fig. 9 The IMF components and the Fourier spectra of El Centro wave decomposed by the folding-line-envelope-based EMD method

3 結論

本文針對原始的經驗模態(tài)分解方法處理長周期信號所面臨的問題，提出了一種基于折線包絡的經驗模態(tài)分解方法，通過對比兩者對信號“篩分”處理過程，論證了基于折線包絡EMD方法的分解效率，以一個具有明確解析表達式的信號為例驗證了該方法所得結果的數學意義，并且分析了該方法分解 El Centro波所得 IMF分量的頻譜特性以及所得到的Hilbert譜的特征。

圖10 El Centro波的Hilbert譜Fig. 10 The Hilbert spectrum of El Centro wave

鄧擁軍，王偉，錢成春等，2001. EMD方法及Hilbert變換中邊界問題的處理. 科學通報，46（3）：257—263.

黃大吉，趙進平，蘇紀蘭，2003. 希爾伯特-黃變換的端點延拓. 海洋學報，25（1）：3—6.

羅奇峰，石春香，2003. Hilbert-Huang變換理論及其計算中的問題. 同濟大學學報，31（6）：637—640.

石春香，羅奇峰，施衛(wèi)星，2005. 基于Hilbert-Huang變換方法的結構損傷診斷. 同濟大學學報（自然科學版），33（1）：17—20.

張郁山，梁建文，胡聿賢，2003. 應用自回歸模型處理EMD方法中的邊界問題. 自然科學進展，13（10）：1054—1059.

Huang N.E., Shen Z., Long S.R. et al., 1998. The Empirical Mode Decomposition and Hilbert Spectrum for Nonlinear and Nonstationary Time Series Analysis. Proc. Roy. Soc. Lond., A (454)∶ 903—995.

Huang N.E., Shen Z. and Long R.S., 1999. A New View of Nonlinear Water Waves—Hilbert Spectrum. Ann. Rev.Fluid Mech., 31∶ 417—457.

Huang N.E. et al., 2001. A New Spectral Representation of Earthquake Data∶ Hilbert Spectral Analysis of Station TCU129, Chi-Chi, Taiwan, 21, September 1999. Bull. Soc. Seism. Am., 91 (5)∶ 1310—1338.

Loh C.H., Wu T.C. and Huang N.E., 2001. Application of the Empirical Mode Decomposition—Hilbert Spectrum Method to Identify Near-fault Ground Motion Characteristics and Structural Responses. Bull. Soc. Seism. Am.,91 (5)∶ 1339—1357.

Pines and Salvino L.W., 2001. Health Monitoring of Structures Using Empirical Mode Decomposition and Phase Dereverberation. See∶ Proc the 2001 Mech and Materials Summer Conf. UCSD∶ San Diego, CA∶ 228—229.

Salvino L.W., 2000. Empirical Mode Analysis of Structural Response and Damping. See∶ Proc the 18th Inter Modal Analysis Conf, San Antonio, TX∶ 503—509.

Yang J.N., Lei Y., Pan S.W. et al., 2003a. System Identification of Linear Structures Based on Hilbert- Huang Spectral Analysis. Part 2∶ Complex Modes. Earthq. Engrg. Struct. Dyn., 32∶ 1533—1554.

Yang J.N., Lei Y., Pan S.W. et al., 2003b. System Identification of Linear Structures Based on Hilbert- Huang Spectral Analysis. Part 1∶ Normal Modes. Earthq. Engrg. Struct. Dyn., 32∶ 1443—1467.

Yang J.N., Lei Y., Lin S. et al., 2003c. Hilbert-Huang Based Approach for Structural Damage Detection. J. Engrg.Mech., ASCE, 130 (1)∶ 85—95.

Zhang R.R., Ma S., Safak E. et al., 2003. HHT Analysis of Earthquake Recordings. J. Engrg. Mech., ASCE, 129 (8)∶861—875.

Zhang Y.S., Liang J.W. and Hu Y.X., 2005. Hilbert Spectrum and Intrinsic Oscillation Mode of Dynamic Response of a Bilinear SDOF System∶ Influence of Harmonic Excitation Amplitude. Earthq. Engrg. Engrg. Vibr., (1)∶17—26.