基于Hilbert-Huang變換與理想帶通濾波器的系統識別*

王其昂， 吳子燕， 劉 露

(西北工業大學力學與土木建筑學院 西安，710129)

?

基于Hilbert-Huang變換與理想帶通濾波器的系統識別*

王其昂， 吳子燕， 劉 露

(西北工業大學力學與土木建筑學院 西安，710129)

針對經驗模態分解存在模態混疊現象，提出基于Hilbert-Huang變換與理想帶通濾波器的系統識別方法。該方法利用傅里葉變換得到結構加速度響應頻響函數，粗略估計固有頻率范圍，通過半功率帶寬法設計理想帶通濾波器，定量化確定通帶帶寬，使信號在經過濾波器后頻域內零相移，同時不改變其幅值譜。結構響應通過指定頻帶的理想帶通濾波器產生若干窄帶信號，利用經驗模態分解獲取結構模態響應，經Hilbert變換構造模態響應解析信號，并通過線性最小二乘擬合提取結構模態參數與物理參數。結果表明：半功率帶寬法可實現帶通濾波器頻帶的定量化設計,理想帶通濾波器的零相移特點較好契合Hilbert-Huang變換用于系統識別的要求，兩者結合可有效地解決模態混疊現象，減少虛假模態，大大提高結構系統識別精度。

系統識別; Hilbert-Huang變換; 理想帶通濾波器; 經驗模態分解; 半功率帶寬法; 模態響應

引 言

工程結構系統識別與損傷診斷是結構健康檢測的核心技術，處于該學科研究前沿。其首要任務是從測得數據中確定結構模態參數，從整體上描述系統動力特性，進一步辨識結構物理參數，為故障診斷、可靠性評估提供理論依據。

系統識別和損傷檢測Benchmark模型在哥倫比亞大學UBC實驗室建立[1]，并發展了卡爾曼濾波、小波變換、子空間識別等一系列方法[2-4]。小波變換[5-6]具有良好的時頻分辨能力，避免傅里葉變換的Gibbs效應，利用時域脈沖響應提取結構模態參數，但其本質上為線性變換，不能有效處理非線性問題。基于經驗模態分解(empirical mode decomposition, 簡稱EMD)的希爾伯特黃變換(Hilbert-Huang transform, 簡稱HHT)[7-8]，具有完全自適應性，不受Heisenberg測不準原理制約等優點，能處理非線性、非平穩信號。當結構兩階固有頻率相似，基于Hilbert-Huang變換的模態參數識別較小波變換更為準確。但EMD方法存在模態混疊，限制了它在實際中的應用，其他學者提出改進的Hilbert-Huang方法[9-10]，利用傅里葉變換，估計固有頻率范圍，然后讓信號通過指定頻帶的帶通濾波器，改善模態混疊，但未涉及濾波器帶寬定量化設計。 Yang等[11]利用Hilbert-Huang變換識別多自由度體系模態參數，并指出該方法用于系統識別中，帶通濾波器相移應盡可能小。

由此，筆者提出基于Hilbert-Huang變換與理想帶通濾波器的系統識別方法，利用傅里葉變換得到結構加速度頻響函數，估計固有頻率范圍。為避免模態混疊，減少虛假模態響應，在使用EMD方法進行模態分解前，使信號通過指定頻帶的理想帶通濾波器，將寬頻信號分解為若干窄帶信號。設計理想帶通濾波器，使結構響應經濾波器后頻域內零相移，同時不改變其幅值譜，并使用半功率帶寬法定量化設計頻帶范圍，減少模態混疊，提高結構系統識別精度。

1 脈沖激勵下模態響應

多自由度系統外力作用下的動力學方程為

(1)

其中：X(t)=[x1,x2, …,xn]T為系統位移向量；M，C，K分別為質量、阻尼和剛度矩陣；p(t)為外部激勵。

通過模態變換，位移、加速度響應分解為n階實模態，考慮模態正交性可得

(2)

其中：ψj為第j階模態向量；qj(t)為模態坐標；ωj，ζj和mj分別為第j階模態頻率、模態阻尼和模態質量。

若有以一激勵作用在第r個自由度上，即當j=r時，pr(t)=p0δ(t)；當j≠r時，pj(t)=0。pj(t)為激勵向量第j個元素。由此可知，第j個模態坐標下加速度響應為

(3)

因此，結構在第s(s= 1, 2, …,n)自由度上加速度脈沖響應可以表示為

(4)

(5)

(6)

其中：φsj,r為第j階模態向量第s，r個元素相位差。

(7)

2 結構響應Hilbert-Huang變換

Hilbert-Huang變換[7]由EMD方法和Hilbert變換組成，其核心思想是將時間序列通過EMD分解成數個本征模函數(intrinsic mode function, 簡稱IMF)，之后利用Hilbert變換構造解析信號，計算得到信號瞬時頻率和振幅。

2.1 EMD概述

EMD將一個復雜的信號分解為若干個IMF之和，每一個IMF均具有相同極值點與過零點數目，相鄰零點之間僅一個極值點，且上下包絡線關于時間軸對稱。對于給定的信號x(t)(如加速度響應)，EMD基本步驟[7]如下：a.確定信號所有局部極值點;b.用三次樣條曲線構造信號上(極大值點)、下(極小值點)包絡線，分別為u(t)和l(t)；c.計算上下包絡線的平均值，記為m1(t)，求出x(t)-m1(t)=h1(t)，若h1(t)滿足IMF特性，即為第1個IMF分量;d. 若h1(t)不滿足IMF的條件，則將其視為原始數據,重復步驟a～c，得到上下包絡線平均值m11(t)，判斷h11(t)=h1(t)-m11(t)是否滿足IMF的條件，如不滿足，則循環k次，得到h1k(t)=h1(k-1)(t)-m1k(t)，使得h1k(t)滿足IMF的條件,記c1(t)=h1k(t)，為信號x(t)的第1個IMF分量;e.將c1(t)從x(t)中分離，得到r1(t)=x(t)-c1(t)，將r1(t)作為原始數據重復以上過程，直至分離出全部IMF，將x(t)分解為n個IMF與殘余分量之和

(8)

2.2 半功率帶寬法設計理想帶通濾波器

為減少虛假模態，在使用EMD方法進行模態分解前，首先使信號通過不同頻段的帶通濾波器，將寬頻信號分解為若干窄帶信號。筆者設計了理想帶通濾波器，使結構響應經濾波器后頻域內零相移，同時不改變其幅值譜，提高模態參數識別精度。

理想帶通濾波器具有完全平坦的通帶，通帶內沒有放大或衰減，通帶外頻率被完全衰減，且通帶外的轉換在極小頻率范圍完成。具體設計思路如下：將加速度響應通過傅里葉變換得頻響函數，通帶范圍內頻響函數數據不變，通帶外變為零，再通過逆傅里葉變換轉化為時域數據，獲得窄帶加速度響應。根據半功率帶寬法，保留信號峰值處多數能量，理想帶通濾波器上、下限截止頻率設計為幅值譜半功率點處，其帶寬的功率譜密度比峰值低3 dB。

2.3 模態響應Hilbert變換及其解析信號

(9)

對于cj(t)的解析信號zj(t)為

(10)

其中：Aj(t)和θj(t)分別為IMF的瞬時振幅和瞬時相位。

(11)

(12)

由瞬時相位計算信號瞬時頻率

(13)

信號瞬時振幅和瞬時相位均為時間函數。測量數據x(t)首先應用EMD分解成單一頻率或窄帶信號IMF，經Hilbert變換可得到有意義的結果，在任意時刻t僅有單一頻率，使得瞬時頻率具有實際物理意義。對于常規信號，瞬時頻率隨時間改變而非單一頻率。

加速度模態響應由式(5)獲得，其Hilbert變換可由Bedrosian定理計算

(14)

其中

(15)

(16)

由式(10)得模態響應解析信號為

(17)

若阻尼較小、固有頻率較大，模態響應Hilbert變換式(14)簡化為

(18)

瞬時幅值、瞬時相位公式簡化為

(19)

(20)

3 系統識別

3.1 固有頻率與阻尼比的識別

3.2 模態振型及物理參數識別

理論上固有頻率、阻尼比的識別，僅需獲取某一自由度上加速度響應,但對于模態振型、質量、剛度和阻尼矩陣的識別則需測得所有自由度上的響應。由式(6)可得到模態元素絕對值之比

(21)

兩模態元素相位差通過式(22)可得到

(22)

從式(21)和式(22)可以看出，將t0作為所取時間段的中間某點，得到θsj和θrj在t0處的擬合值，進而可以確定ψsj相對于ψrj的相位。根據式(7)，確定ψsj相對于ψrj符號(正或負)。由此，在j階模態向量中的所有元素相對于某一特定元素的絕對值及符號都可確定。

得到模態參數后，結構質量、剛度和阻尼矩陣便可得到。將式(19)、式(20)帶入式(6)中，第j階模態質量為

(23)

其中:p0為可測沖擊力荷載。

模態剛度ki、阻尼ci為

(24)

根據模態正交性得結構質量、剛度及阻尼陣為

(25)

其中:Φ為模態向量矩陣。

4 算例分析

圖1 結構系統識別流程圖Fig.1 Flow chart for the structural system identification

圖2 4自由度剪切梁模型Fig.2 Four degrees of freedom shear-beam structure

筆者以某4層剪切梁模型為算例，進行實驗模態分析，并與理論分析結果作對比，以驗證上述系統識別方法的有效性及準確性。具體流程見圖 1。4層剪切梁模型見圖 2，樓層質量均為4.989 5 kg，樓層剛度k1，k2，k3與k4分別為1 401，1 576，1 226與1 051 N/m，阻尼類型為瑞利阻尼。由質量矩陣和剛度矩陣按比例組合構造而成，質量比例系數Alpha為-0.001 6，剛度比例系數Beta為3.918 8×10-4。在第4樓層施加單位脈沖激勵，構造系統狀態空間函數公式。利用Matlab系統仿真獲取各樓層加速度，而實際上任何實測數據都不可避免受到噪聲信號的影響,因此本研究在響應信號中添加相同長度的高斯白噪聲，信噪比為26.02 dB。

該系統4階自振頻率及對應阻尼比均可由任一樓層加速度脈沖響應獲取。圖 3為頻域內加速度脈沖響應幅值譜，響應中主要包含4階模態的振動成分，固有頻率大致在0.925,2.475,3.85和4.95 Hz左右。

圖3 脈沖響應頻域幅值譜Fig.3 Frequency spectra of impulse responses

圖4 帶通濾波器相移圖Fig.4 Phase shift of different band-pass filters

圖5 4階模態響應Fig.5 Four modal responses

圖6 解析信號相位圖Fig.6 Phase plot of the analytical signal

圖7 解析信號幅值圖Fig.7 Amplitude plot of the analytical signal

由固有頻率和阻尼比識別步驟3得第1階模態頻率、阻尼。同時，對于其余3層響應可得第1階模態響應，通過Hilbert-Huang變換理論計算可得第1階模態參數，取均值得第1階固有頻率、阻尼比。重復同樣步驟得其余3階模態頻率、阻尼比，并與Butterworth濾波器識別結果相比較，結果如表 1所示，由式(7)、式(21)確定模態向量。

表1 結構模態參數識別結果

模態頻率識別，兩種濾波器估算結果一致，最大誤差均小于0.2%。模態阻尼識別對數值誤差較敏感，尚缺乏較為理想的識別方法，理想帶通濾波器(最大誤差23%)略優于Butterworth濾波器(最大誤差28%)。模態向量識別結果與理論值相比得模態置信度(modal assurance criterion，簡稱MAC)值(見表 1)。4階MAC值均大于0.99(實踐中，兩模態向量MAC值大于0.9時為強相關)，基于理想帶通濾波器的識別結果與理論振型基本重合，模態向量識別精度高。

與常規Butterworth帶通濾波器相比，理想帶通濾波略提高模態參數識別精度。由于誤差累積效應，這一提高對于結構物理參數的識別至關重要。1～4樓層質量識別結果分別為4.941，5.047，4.892和4.999kg，最大誤差為1.95%，而Butterworth濾波器最大誤差為11.83%。1～4樓層剛度識別結果分別為1 323.1，1 654.0，1 166.4和1 037.1，最大誤差為5.56%，Butterworth濾波器最大誤差為19.31%。阻尼矩陣主對角元素識別結果分別為Cd(1,1)=1.08，Cd(2,2)=0.98，Cd(3,3)=0.84，Cd(4,4)=0.41，最大誤差為10%，常規濾波器識別結果誤差為43.1%。為驗證該方法識別靈敏度，在原有算例基礎上分別修改質量、剛度矩陣。樓層質量修改為5.2 kg，以剛度識別結果為例，理想帶通濾波器最大誤差為3.7%，而Butterworth濾波器識別結果最大誤差為11.4%。修改1～4樓層剛度，理想帶通濾波器識別結果同樣優于常規濾波器。由此，對于結構物理參數識別，理想帶通濾波器處理后有明顯優勢，大大提高了識別精度。

5 結 論

1) 筆者提出基于改進的Hilbert-Huang變換的系統識別方法，利用傅里葉變換得到結構加速度頻響函數，粗略估計固有頻率范圍。為減少虛假模態，使用EMD模態分解前，使信號通過指定頻帶的理想帶通濾波器，將寬頻信號分解為若干窄帶信號，構造模態響應解析信號，提取模態參數與物理參數。

2) 理想帶通濾波器的零相移特點較好契合Hilbert-Huang變換用于系統識別的要求，兩者結合有效地解決了EMD中模態混疊問題。

3) 獲取加速度響應幅值譜，由半功率帶寬法，確定理想帶通濾波器上、下限截止頻率為幅值譜半功率點處，實現頻帶定量化設計。

4) 與常規Butterworth帶通濾波器相比，理想帶通濾波略微提高模態參數識別精度。但由于誤差累積效應，這一提高對于結構物理參數的識別至關重要，大大提高了結構質量、剛度及阻尼矩陣的識別精度。

[1] Black C J, Ventura C E. Blind test on damage detection of a steel frame structure[C]∥Proceedings of SPIE. Santa Barbara CA: Society for Experimental Mechanics, 1998: 623-629.

[2] 任宜春，易偉建. 結構物理參數識別的多尺度參數卡爾曼濾波方法[J]. 工程力學，2008，25(5)：1-5.

Ren Yichun, Yi Weijian. Identification of physical parameters by multi-scale parameter Kalman filter[J]. Engineering Mechanics, 2008, 25(5): 1-5. (in Chinese)

[3] 史治宇，沈林. 基于小波方法的時變動力系統參數識別[J]. 振動、測試與診斷，2008，28(2): 108-112.

Shi Zhiyu, Shen Lin. Parameter identification of linear time-varying dynamical system based on wavelet method[J]. Journal of Vibration, Measurement & Diagnosis, 2008, 28(2): 108-112. (in Chinese)

[4] 徐良，江見鯨，過靜珺. 隨機子空間識別在懸索橋實驗模態分析中的應用[J]. 工程力學，2002，19(4): 46-49.

Xu Liang, Jiang Jianjing, Guo Jingjun. Application of stochastic subspace method to experimental modal analysis of suspension bridges[J]. Engineering Mechanics, 2002, 19(4): 46-49. (in Chinese)

[5] Vafaei M, Alih S C, Abd Rahman A B, et al. A wavelet-based technique for damage quantification via mode shape decomposition[J]. Structure and Infrastructure Engineering, 2015, 11(7): 869-883.

[6] Dziedziech K, Staszewski W J, Uhl T. Wavelet-based modal analysis for time-variant systems[J]. Mechanical Systems and Signal Processing, 2015, 50-51: 323-337.

[7] 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]. Oyal Society of London Proceedings, 1998, 454(1971): 903-995.

[8] Lee J, Hussain S H, Wang S, et al. Enhanced method to reconstruct mode shapes of continuous scanning measurements using the Hilbert Huang transform and the modal analysis method[J]. Review of Scientific Instruments, 2014, 85(9): 1-11.

[9] 付春，姜紹飛. 基于改進EMD-ICA的結構模態參數識別研究[J]. 工程力學，2013，30(10)：199-204.

Fu Chun, Jiang Shaofei. Study on the modal parameter identification based on improved EMD and independent component analysis[J]. Engineering Mechanics, 2013, 30(10): 199-204. (in Chinese)

[10]秦毅，秦樹人，毛永芳. 改進的Hilbert-Huang變換在信號瞬態特征提取中的應用[J]. 振動與沖擊， 2008，27(11)：129-133.

Qin Yi, Qin Shuren, Mao Yongfang. Application of improved Hilbert-Huang transform in transient feature extraction[J]. Journal of Vibration and Shock, 2008, 27(11): 129-133. (in Chinese)

[11]Yang J N, Lei Ying, Pan Shuwen, 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(9): 1443-1467.

[12]Shi Z Y, Law S S, Xua X. Identification of linear time-varying mdof dynamic systems from forced excitation using Hilbert transform and EMD method[J]. Journal of Sound and Vibration, 2009, 321(3-5): 572-589.

10.16450/j.cnki.issn.1004-6801.2016.06.005

*國家自然科學基金資助項目(51278420)；博士創新基金資助項目(CX201408)

2015-06-08；

2015-07-13

TB122; TB123

王其昂,男，1986年8月生，博士生。主要研究方向為結構健康檢測、系統識別及可靠性評估。曾發表《Seismic fragility analysis of highway bridges considering multi-dimensional performance limit state》(《Earthquake Engineering and Engineering Vibration》2012,Vol.11,No.2)等論文。E-mail:qawang@mail.nwpu.edu.cn