基于改進HHT方法的密集模態(tài)結構參數識別

黃天立，邱發(fā)強，樓夢麟

(1. 中南大學 土木工程學院，湖南 長沙，410075；2. 同濟大學 土木工程防災國家重點實驗室，上海，200092)

Hilbert-Huang變換(Hilbert-Huang transform，HHT)[1-2]是一種全新的數據處理方法，由經驗模式分解(Empirical mode decomposition，EMD)方法及 Hilbert變換(HT) 2部分組成，其基本思想是將時間序列通過經驗模式分解成數個固有模式函數(Intrinsic mode function，IMF)，然后，利用 Hilbert變換構造解析信號，得出時間序列的瞬時頻率和振幅，進而得到Hilbert譜。EMD分解是依賴數據本身的時間尺度特征進行的，比傅里葉分析及小波分析更適合于處理非平穩(wěn)數據。目前，此方法已在振動工程、結構損傷診斷、系

統(tǒng)識別[3-7]等許多工程領域得到廣泛應用。Yang等[3]利用 HHT方法系統(tǒng)研究了線性多自由度結構系統(tǒng)的參數識別問題，提出了識別結構的模態(tài)參數(頻率、振型、阻尼)和物理參數(質量、剛度、阻尼矩陣)的方法；陳雋等[4-5]研究了 HHT方法在模態(tài)參數識別中的應用，并特別針對密頻結構的阻尼識別問題進行探討，發(fā)現基于 HHT的方法較傳統(tǒng)方法具有良好的識別密頻結構阻尼的性能；黃天立等[6]研究了HHT方法在非線性多自由度結構系統(tǒng)識別中的應用；Yan等[7]比較研究了小波變換和 HHT方法在結構模態(tài)參數識別中的相關問題。從這些基于 HHT方法的模態(tài)參數識別研究結果可以看出，HHT方法能否較好地識別出結構的模態(tài)參數在很大程度上取決于 2個因素：(1)EMD的模態(tài)分解能力或帶通濾波等前處理工具的模態(tài)分解能力；(2)EMD分解所得代表結構模態(tài)響應的IMF分量的質量。關于 EMD的模態(tài)分解能力問題，大量的實踐和理論分析結果[8]表明：EMD的模態(tài)分解能力與信號中不同模態(tài)的頻率之比密切相關。2個頻率成分之間的頻率之比越大，則 EMD完全分離這兩階模態(tài)的可能性越大，反之越小。為此，本文作者采用波組(Wave group)信號前處理工具[9]予以處理，通過對原始信號進行相應處理，使原始信號中各密集頻率成分之間的比值加大，從而能很好地適應 EMD有限的模態(tài)分解能力。關于EMD分解所得IMF分量的質量問題，文獻[10]的研究結果表明：EMD分解不能保證信號分解的能量守恒，即所得IMF分量之間的正交在一定程度上影響了模態(tài)參數識別結果的穩(wěn)定性和精確性。因此，本文采用正交化經驗模式分解(Orthogonal empirical mode decomposition，OEMD)方法[10]，以獲得完備且完全正交的IMF分量，使其更接近結構的真實模態(tài)響應，從而使得模態(tài)參數識別結果更加穩(wěn)定、精確。為敘述方便，將波組信號前處理結合正交化經驗模式分解的方法稱為改進的HHT方法。

1 正交化經驗模式分解

1.1 EMD分解

設X(t)為要分解的時程信號，EMD分解信號的步驟如下。

(1) 確定X(t)的所有極大、極小值點，用三次樣條函數曲線分別連接極大、極小值點形成上、下包絡線，求出上、下包絡線的均值線 m1(t)，定義 X(t)與 m1(t)的差為h1(t)=X(t)-m1(t)。若h1(t)滿足如下2個條件：① 在整個時程中，極值點數目和過零點數目相等或最多相差1個；② 在任意點，由局部極大值點和局部極小值點構成的2條包絡線平均值為0，則h1(t)為第1階IMF分量。一般來說，h1(t)并不滿足上述條件，則將h1(t)看作新的時程曲線，重復上述篩選過程。假設重復k次后，h1k(t)= h1(k-1)(t)- m1k(t)滿足上述條件，則令c1(t)= h1k(t)為X(t)的第1階IMF分量。

(2) 令r1(t)=X(t)-c1(t)并將其視作新的時程信號，重復上述過程，即可得到第2階IMF分量c2(t)。當第n階IMF分量cn(t)滿足：① cn(t)或rn(t)小于預先設定的數值；② 殘差rn(t)成為單調函數，不可能再從中提取固有模式函數，或滿足這2個條件之一時，則整個篩選過程完成。

由此，X(t)被分解為n個IMF分量及余量rn(t)之和。

式(1)中，信號皆表示時間的連續(xù)形式?？紤]到采集的信號都是離散的，是在N個時間點上的不同數值，這里給出信號的離散向量表達形式：

1.2 OEMD分解

通過對EMD分解得到的各階IMF分量進行正交化處理，可得到完全正交的各階IMF分量，其步驟如下。

(1) 由EMD分解得到時程信號X(t)的第1階IMF分量，記作c1(t)。令c1( t) = c1(t)，稱為時程信號X(t)的第1階正交化IMF分量。

(2) r1(t)=X(t)-c1(t)，將r1(t)當作一新的時程信號進行 EMD分解，可得到第 2階初始 IMF分量，記作c2(t)。從形成c2(t)的過程中可以看出：并不能保證c2(t)與c1(t)對時間的分布具有正交的性質。為了得到時程信號X(t)的第2階正交化IMF分量，應從c2(t)中消除所含的c1(t)分量，即

式中：c2(t)稱為時程信號X(t)的第2階正交化IMF分量；β21稱為c2(t)與 c1(t)之間的正交化系數。為了得到β21，可將式(2)的兩邊同乘以c1(t)并對時間t進行積分，利用c2(t)與c1(t)的正交性，得到：

將β21表示成離散形式：

(3) 采用與上述相同的方法，從EMD分解得到時程信號X(t)的第j+1階初始IMF分量中消除所有前j階正交化IMF分量，則可得到時程信號X(t)的第j+1階正交化IMF分量cj+1(t)(j=2, …, n)。具體的計算過程如下。

(a) 計算：

將rj(t)當作一新的時程信號進行EMD分解，可得到第j+1階初始IMF分量，記作。正交化計算公式為：

(b) 為了得到1,jiβ+，將式(6)兩邊同乘以ck(t)(k≤j)并對時間t進行積分，利用ck(t)與ci(t)(i≠k)以及cj+1(t)間的正交性，得到：

當i=k時，即可得到：

(c) 將1,jiβ+表示成離散形式：

當滿足：①cn(t)或 rn(t)小于預先設定的數值和②殘差rn(t)成為單調函數這2個條件之一時，則分解終止，整個篩選過程完成。

經過上述步驟以后，X(t)被分解成如下形式：

顯然，分量cj(t)(j=1, 2, …, n)之間是完全正交的，將各 cj(t)分量作線性變換不會改變各分量之間的正交性，因此，)(t (j=1, 2, …, n)之間是完全正交的。這樣，信號X(t)被分解成為n個正交IMF分量)(t (j=1,2, …, n)及余量rn(t)之和。

2 波組信號前處理方法

EMD的模態(tài)分解能力與信號中不同模態(tài)頻率之比密切相關，2個頻率成分之間的頻率之比越大，則EMD完全分離這兩階模態(tài)的可能性越大，反之便越小。為了增加信號不同頻率之間的比值，Wang[9]提出了一個簡單的針對原始信號進行前處理的方法，其基本思想如下。

考慮 1條由 2個不同頻率余弦信號構成的信號x(t)，即所謂波組(Wave Group，WG)信號，

定義頻率比α=ω2/ω1，當2個頻率成分比較靠近時，頻率之比較小。若將頻率比的分子和分母同時減去 1個合適的瞬時頻率ω0，則頻率比值α=(ω2-ω0)/(ω1-ω0)被顯著放大。對于實值信號，可通過如下步驟來實現這一基本思想。

(1) 求取實值信號 x(t)的 Hilbert變換 H[x(t)]，構造相應的復解析信號X(t)為：

(2) 將復解析信號X(t)乘以指數函數 e-iω0t構造新的復解析信號Z(t)：

這樣就實現了信號頻率的下移(Downshifting)。

(3) 通過選擇合適的0ω，使得相應的頻率比α增大，則Zr(t)和Zj(t)就能由EMD分解為：

式中：Crk(t)和 Cjk(t)分別為 Zr(t)和 Zj(t)的 IMF分量；rnr和rnj為相應的余量。

(4) 將復解析信號Z(t)進一步表示為：

式中：Ck(t)和Rn為復數。

(5) 由式(12)可知：

因此，原始信號x(t)可表示為：

3 基于改進HHT方法的結構模態(tài)參數識別

從波組信號前處理方法中看出，該方法的基本思想是對從原始信號構造出的對應解析信號的實部和虛部分別進行 EMD分解，而避免了直接對原始信號進行分解所帶來的 EMD頻率分離能力不足問題。由于要進行2次EMD分解，故EMD分解所帶來的能量不守恒問題將更加嚴重，因而，必須采用正交化經驗模式分解方法對EMD分解所產生的IMF分量進行正交化處理。將正交化經驗模式分解方法融入上述波組信號前處理方法之中，對EMD分解所得IMF分量進行正交化處理。最后，從分解所得Ck(t)中選取代表相應模態(tài)響應的IMF分量，按基于Hilbert變換(HT)的模態(tài)參數識別方法[6]即可識別出相應的模態(tài)參數。

4 數值算例

考慮圖1所示三自由度密集模態(tài)線性系統(tǒng)，系統(tǒng)物理參數如下：m1=m2=m3=1 t，k1=k4=30 kN/m，k2=k3=15 kN/m，c1=c4=120 N·s/m，c2=c3=45 N·s/m。計算得到該系統(tǒng)的各階頻率f1分別為0.672 9，1.180 3和1.303 9 Hz；模態(tài)阻尼比1ξ分別為0.64%，1.11%和1.23%。

圖2所示為在m3上作用脈沖荷載時的加速度響應時程曲線(信號采樣頻率為10 Hz，時間長度為100 s)及其Fourier幅值譜。由圖2可見：結構的第2和第3階模態(tài)比較密集。

圖1 三自由度密集模態(tài)系統(tǒng)Fig.1 3-DOF system with closely spaced modes

圖2 加速度響應及Fourier幅值譜Fig.2 Acceleration response and Fourier amplitude spectrum

對該加速度響應時程進行EMD分解所得前2階IMF分量及其Fourier幅值譜如圖3所示。由第1階IMF分量的Fourier幅值譜可知：第1階IMF分量中包含結構的第2和第3階模態(tài)，這表明EMD不具備將該信號的第2和第3階密集模態(tài)分離的能力。

本文分別采用3種方法研究該密集模態(tài)結構的參數識別問題。

4.1 帶通濾波+EMD+HT

為增強 EMD的分解能力，采用帶通濾波作為前處理工具，而后采用基于 EMD+HT的方法進行模態(tài)參數識別。針對圖 2所示加速度響應信號，選擇0.50～0.80，1.10～1.20 和 1.25～1.35 Hz 3 個帶通濾波范圍，分別將該信號進行帶通濾波，得到3條窄帶信號，然后，將這3條窄帶信號進行EMD分解得到與之相關系數最大的 IMF分量代表結構的模態(tài)響應，基于HT的模態(tài)參數識別方法即可識別出模態(tài)參數，識別結果如表1所示。

從表1可以看出：第1階模態(tài)的頻率和阻尼識別結果較好，而第2和第3階模態(tài)的頻率和阻尼識別效果均不理想。這主要是結構第2和第3階模態(tài)頻率比較靠近，使得帶通濾波器的參數取值困難且設計出的帶通濾波器效果較差而造成的，因而需要進一步改進。

圖3 經驗模式分解Fig.3 Empirical mode decomposition

表1 密集模態(tài)結構參數識別結果Table 1 Results of modal parameters of structure with closely spaced modes

4.2 濾波+WG+EMD/OEMD+HT

為了進一步增強 EMD的分解能力，在濾波的基礎上增加波組信號前處理(WG)，而后采用EMD/OEMD+HT的方法進行模態(tài)參數識別。

(1) 針對信號中的第 1階模態(tài)響應，選用 10階Chebyshev Ⅱ型低通濾波器，設定截止頻率為0.8 Hz，得到結構的第1階模態(tài)響應，對其進行EMD分解，選擇相應的可代表此模態(tài)響應分量的IMF分量，按前述基于HT的模態(tài)參數識別方法進行識別，識別結果見表1。結果表明：第1階模態(tài)的頻率和阻尼識別結果很好。

(2) 針對信號中的第2和第3階模態(tài)響應，選用10階Chebyshev Ⅱ型高通濾波器，截止頻率為0.8 Hz得到結構的響應 x(t)及其 Fourier幅值譜，如圖 4所示。從圖4可看出：結構第2和第3階模態(tài)皆包含在此分量中。

(3) 根據前述波組信號前處理的基本原理，選擇ω0= 1.0 Hz，使得此兩階密集頻率比由 α =ω3/ω2≈1.1提高到 α =(ω3-ω0)/(ω2-ω0)≈1.6，將圖4所示高通濾波信號x(t)進行Hilbert變換，構造出相應的解析信號X(t)，并乘以指數函數 e-iω0t構造出新的復解析信號 Z (t) = X(t) exp(- iω0t)。Z(t)的實部和虛部如圖 5所示。

(4) 對Z(t)的實部和虛部分別進行EMD分解，得到Z(t)實部和虛部前2階IMF分量；采用OEMD處理后得到實部和虛部前2階正交IMF分量，如圖6所示。由圖6可見：二者在形狀上差別不大，但幅值上有些差別。這是EMD分解中存在能量泄漏造成的。

圖4 高通濾波后結構響應x(t)及其Fourier幅值譜Fig.4 x(t) response after high-pass filtering and Fourier amplitude spectrum

圖5 Z(t)的實部和虛部Fig.5 Real and imaginary part of Z(t)

圖6 Z(t)實部和虛部IMF及其正交IMF分量Fig.6 IMF/OIMF components of real and imaginary part of Z(t)

(5) 將按IMF分量或正交IMF分量展開的Z(t)乘以指數函數 eiω0t即得到以IMF分量或正交IMF分量表示的原始信號，即

圖7和圖8所示分別為x(t)的前2階IMF/正交IMF分量及其Fourier幅值譜。由圖7和圖8可見：由于采用了波組信號前處理方法，第2和第3階模態(tài)響應得到了較好分離。對比IMF和正交IMF分量可見：正交IMF分量c1比IMF分量c1有所提高，表現在部分消除了c1分量Fourier幅值譜中代表第2階模態(tài)頻率的小峰值；正交IMF分量c2與IMF分量c2相比沒有太大變化。

(6) 對此3階IMF/正交IMF分量按基于HT的模態(tài)參數識別方法即可識別出相應的模態(tài)參數，相應的IMF/正交 IMF分量對數幅值和相位曲線如圖 9和圖10所示(圖中虛線為實際曲線，實線為最小二乘擬合曲線)，識別結果見表1。從表1可以看出：

圖7 x(t)前2階IMF分量及其Fourier幅值譜Fig.7 The first two IMF component of x(t) and their Fourier amplitude spectrum

圖8 x(t)前2階正交IMF分量及其Fourier幅值譜Fig.8 The first two OIMF component of x(t) and their Fourier amplitude spectrum

圖9 IMF分量對數幅值和相位曲線Fig.9 Log amplitude and phase curve of IMF component

圖10 正交IMF分量對數幅值和相位曲線Fig.10 Log amplitude and phase curve of OIMF component

(1) 由于采用了波組信號前處理使得密集模態(tài)得到了很好分離，因而基于“濾波+WG+EMD/OEMD+HT”的模態(tài)參數識別方法優(yōu)于基于“帶通濾波+EMD+HT”的模態(tài)參數識別方法，且識別結果更接近理論值。

(2) 由于基于濾波+WG+OEMD獲得的正交 IMF分量在能量上很好地修正了基于濾波+WG+EMD獲得的IMF分量，因此，基于“濾波+WG+OEMD+HT”的阻尼識別結果優(yōu)于基于“濾波+WG+EMD+HT”的阻尼識別結果，且識別結果更接近理論值。這表明對IMF分量進行正交化處理有利于阻尼的精確識別。

(3) 在頻率識別方面，基于“濾波+WG+OEMD+HT”的識別方法與基于“濾波+WG+EMD+HT”的識別方法沒有明顯改進。

5 結論

(1) 針對HHT方法在識別結構模態(tài)參數中存在的EMD模態(tài)分解能力不足以及IMF分量之間不正交這2個問題，提出了采用波組信號前處理和正交化經驗模式分解方法的改進措施即改進的 Hilbert-Huang變換方法，并將其應用于具有低頻密集模態(tài)的結構參數識別。

(2) 由于波組信號前處理對 EMD模態(tài)分解能力增強，在低頻密集模態(tài)結構的參數識別中，基于改進HHT方法的模態(tài)參數識別方法優(yōu)于基于HHT的模態(tài)參數識別方法。

(3) 由于正交化經驗模式分解對IMF分量在能量上進行修正，基于“濾波+WG+OEMD+HT”的阻尼識別結果優(yōu)于基于“濾波+WG+EMD+HT”的阻尼識別結果，且識別結果更接近理論值。

[1] Huang N E, Shen Z, Long S R, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J]. Proc R Soc Lond A, 1998,454: 903-995.

[2] Huang N E, Shen Z, Long S R. A new view of nonlinear water waves: the Hilbert Spectrum[J]. Annual Review of Fluid Mechanics, 1999, 31: 417.

[3] Yang J N, Lei Y, Pan S, et al. System identification of linear structures based on Hilbert Huang spectral analysis. Part 1:Normal modes[J]. Earthquake Engineering and Structural Dynamics, 2003, 32: 1443.

[4] 陳雋, 徐幼麟. HHT方法在結構模態(tài)參數識別中的應用[J].振動工程學報, 2003, 16(3): 383-388.CHEN Jun, XU You-lin. Application of HHT for modal parameter identification to civil structures[J]. Journal of Vibration Engineering, 2003, 16(3): 383-388.

[5] 陳雋, 徐幼麟, 李杰. Hilbert-Huang變換在密頻結構阻尼識別中的應用[J]. 地震工程與工程振動, 2003, 23(4): 34-42.CHEN Jun, XU You-lin, LI Jie. Hilbert-Huang transform for damping ratio identification of structures with closely spaced modes of vibration[J]. Earthquake Engineering & Engineering Vibration: Chinese Edition, 2003, 23(4): 34-42.

[6] 黃天立, 樓夢麟. 基于 HHT的非線性結構系統(tǒng)識別研究[J].地震工程與工程振動, 2006, 26(3): 48-52.HUANG Tian-li, LOU Meng-lin. System identification of nonlinear structures based on HHT[J]. Earthquake Engineering& Engineering vibration: Chinese Edition, 2006, 26(3): 48-52.

[7] Yan B F, Miyamoto A. A comparative study of modal parameter identification based on wavelet and Hilbert-Huang transforms[J].Computer-Aided Civil and Infrastructure Engineering, 2006, 21:9-23.

[8] 沈國際, 陶利民, 陳仲生. 多頻信號經驗模態(tài)分解的理論研究及應用[J]. 振動工程學報, 2005, 18(1): 91-94.SHEN Guo-ji, TAO Li-min, CHEN Zhong-sheng. Theoretical research on empirical mode decomposition of multi-frequency signal and its application[J]. Journal of Vibration Engineering,2005, 18(1): 91-94.

[9] Wang W. Decomposition of wave groups with EMD method[C]//Huang N E. The Hilbert-Huang transform in engineering. Taylor& Francis Group. London: CRC Press, 2005: 267-280.

[10] 樓夢麟, 黃天立. 正交化經驗模式分解方法[J]. 同濟大學學報: 自然科學版, 2007, 35(3): 293-298.LOU Meng-lin, HUANG Tian-li. The orthogonal empirical mode decomposition[J]. Journal of Tongji University: Natural Science, 2007, 35(3): 293-298.