楊彩紅 梁建文 張郁山

1）中國(guó)地震災(zāi)害防御中心，北京 100029

2）天津大學(xué)土木工程系，天津 300072

引言

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

針對(duì)傳統(tǒng)的EMD方法分解長(zhǎng)周期信號(hào)所面臨的問題，本文提出了基于折線包絡(luò)的經(jīng)驗(yàn)?zāi)B(tài)分解方法，通過與已有EMD方法的比較論證了該方法的分解效率，并利用數(shù)值算例驗(yàn)證了該方法所得結(jié)果的意義。

1 基本方法

1.1 基于樣條包絡(luò)的“篩分”

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

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

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

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

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

圖2 原始信號(hào)及其兩種平均包絡(luò)Fig. 2 Original signal and its two kinds of mean envelopes

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

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

1.2 基于折線包絡(luò)的“篩分”

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

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

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

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

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

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

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

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

1.3 基于折線包絡(luò)的EMD方法

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

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

1）記對(duì)原始信號(hào)x(t)進(jìn)行k1次與k1-1篩分所得結(jié)果為若它們之間的標(biāo)準(zhǔn)差：

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

2）令

若對(duì) r1( t)篩分k2次后的結(jié)果與k2-1次的結(jié)果之間的標(biāo)準(zhǔn)差（按照公式（4）計(jì)算）滿足Δs＜0.3，則為原始信號(hào)的第二個(gè)IMF分量：

3）將上述分解過程重復(fù)，即可得到原始信號(hào)的一系列IMF分量：

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

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

圖5 基于折線包絡(luò)的EMD方法分解圖4所示信號(hào)所得到的主要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 具有明確解析式的信號(hào)

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

與高頻正弦噪聲

共同合成。其中信號(hào)x1(t)的瞬時(shí)頻率fi1(t)定義為其相位的導(dǎo)數(shù)

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

圖6 基于折線包絡(luò)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)

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

圖7 原始信號(hào)的Hilbert譜Fig. 7 Hilbert spectrum of original signal

2.2 地震動(dòng)加速度記錄

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

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

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

圖9 基于折線包絡(luò)的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 結(jié)論

本文針對(duì)原始的經(jīng)驗(yàn)?zāi)B(tài)分解方法處理長(zhǎng)周期信號(hào)所面臨的問題，提出了一種基于折線包絡(luò)的經(jīng)驗(yàn)?zāi)B(tài)分解方法，通過對(duì)比兩者對(duì)信號(hào)“篩分”處理過程，論證了基于折線包絡(luò)EMD方法的分解效率，以一個(gè)具有明確解析表達(dá)式的信號(hào)為例驗(yàn)證了該方法所得結(jié)果的數(shù)學(xué)意義，并且分析了該方法分解 El Centro波所得 IMF分量的頻譜特性以及所得到的Hilbert譜的特征。

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

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

黃大吉，趙進(jìn)平，蘇紀(jì)蘭，2003. 希爾伯特-黃變換的端點(diǎn)延拓. 海洋學(xué)報(bào)，25（1）：3—6.

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

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

張郁山，梁建文，胡聿賢，2003. 應(yīng)用自回歸模型處理EMD方法中的邊界問題. 自然科學(xué)進(jìn)展，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.