基于加速度響應連續小波變換的線性時變結構瞬時頻率識別

許 鑫，史治宇，Wieslaw J.Staszewski，龍雙麗

(1.南京航空航天大學 飛行器結構力學與控制教育部重點實驗室，南京 210016;2.Dynamics Research Group，Department of Mechanical Engineering，University of Sheffield，Sheffield S1 3JD，England)

無論是航空航天、高速列車，還是土木橋梁、造船機械，力學結構在其服役期內總會由于載荷，環境，人為等因素的影響出現結構的力學參數隨時間變化的特性。這種時變特性所引起的力學問題近一二十年來備受關注。而結構參數識別是力學建模和力學分析的基礎，時變結構參數識別也由此成為近些年來參數識別領域的一個熱點問題。

國內外學者提出了不少的時變結構參數識別方法，它們都基于不同的信號處理技術，大致可以分為:經驗模態分解(EMD)和 HT(Hilbert Transform)變換[1－2]，狀 態 空 間 法 (State Space)[3－4]，小 波 變 換(Wavelet Transform)[5－7]，功能序列時變 ARMA 模型(FS － ARMA)[8－9]以及一些自適應算法[10－11]。

這些方法中，小波變換由于能夠自適應調節窗口的大小，具有多尺度分析的能力，可以對非穩態信號做精細的時頻局部處理，是目前非穩態信號分析中應用最廣泛的技術手段。連續小波變換最早是由Staszewski和 Cooper[12]應用到參數識別領域當中的。Le 等[13]根據小波脊與模態參數之間的對應關系進行了結構參數識別，并對比研究了三種分析小波函數的特性。對于時變問題，王超等[14]運用Morlet小波識別了時變結構的瞬時頻率。在離散小波方面，Ghanem等[6]運用Daubechies小波識別了時變系統的物理參數。Chen等[7]利用 Haar小波提取時變系統的參數。許鑫等[15－16]又將 Daubechies小波引入狀態空間模型，無需計算二階小波連接系數，達到了識別時變結構參數的目的。

本文首先推導了函數積分運算的連續小波變換方法，借助此法，僅利用時變結構自由振動的加速度響應信號，就可以計算出速度信號和位移信號的連續小波變換值，通過將一小段時刻點的線性方程組聯合構造最小二乘問題，識別了時變結構不同時刻點的阻尼和剛度，確定了結構的瞬時頻率。文中的5層剪切梁樓房模型和3自由度密集模態算例驗證了本方法的正確性，有效性和對噪聲的魯棒性。

本文方法與其他基于小波的識別方法[6－7，15]相比，優勢在于僅測量加速度信號和識別效率高。文獻[7]的方法需要測量時變結構的位移信號，而文獻[15]的方法則需要同時知道位移和速度兩種信號，這都給振動信號的測量帶來了不便。為了方便和實用，可以僅測量加速度信號，再將其數值積分，得到速度和位移信號，但這樣做:① 必須知道速度和位移信號的初始條件，即積分初值，這難度大，且精度沒保證;② 即使知道了積分初值(比如說零初值)，受數值積分的影響，積分誤差在所難免，并且還會累計，由此參數識別誤差更大。文獻[6]的方法需要同時計算Daubechies小波的一階和二階小波連接系數，識別過程復雜，計算時間耗費大，比文獻[15]的效率還要低[18]。

1 函數積分運算的連續小波變換

1.1 連續小波變換

函數y(t)的連續小波變換通常表示為:

其中a和b分別是尺度參數和平移參數，ψ(t)為小波分析函數(或稱母小波)，上式取其共軛，{Wψy}(a，b)符號表示利用小波分析函數ψ(t)通過小波變換將函數y(t)映射到小波平面。

本文假設小波分析函數具有兩階或兩階以上的消失矩，即滿足:

設小波分析函數ψ(t)存在一階和二階積分原函數Ψ1(t)和Ψ2(t)，并且小波分析函數和它的兩階積分原函數都在±∞處收斂為零，即:

1.2 函數積分運算的連續小波變換計算

設函數y(t)有一階積分原函數Y1(t)，則y(t)一階積分運算的連續小波變換可記為:

其中y0為積分常數，對上式用分部積分公式，

利用小波分析函數的特性(2)、(3)式，不難發現(5)式中等號右邊的第一項和第三項等于零，進而得出函數一階積分運算的連續小波變換計算方法:函數y(t)一階積分運算的連續小波變換可以通過小波分析函數ψ(t)的一階積分原函數Ψ1(t)對y(t)進行連續小波變換計算。公式表示為:

同理，函數y(t)二階積分運算的連續小波變換可以通過小波分析函數ψ(t)的二階積分原函數Ψ2(t)對y(t)進行連續小波變換計算。即:

式中:Y2(t)是函數y(t)的二階原函數，y1和y2是積分常值。

以上函數積分運算的連續小波變換方法，是由Sone等[17]提出的，他們借助這個想法識別了時不變結構的力學參數，本文則要把這個方法推廣運用到時變結構。

2 線性時變結構的瞬時頻率識別

2.1 線性時變結構自由振動微分方程

一個p自由度線性時變動力學結構，其自由振動可以用下列方程式表示:

其中:M(t)、E(t)、K(t)分別為p×p的時變質量矩陣、阻尼矩陣、剛度矩陣，x(t)為p×1的位移向量。速度向量和位移向量x(t)均可由加速度向量積分表示:

2.2 時變結構物理參數確定

將等式(9)、式(10)代入振動微分方程(8)，選取滿足條件(2)、條件(3)的小波分析函數ψ(t)，對式(8)兩邊進行連續小波變換，進而利用式(6)、式(7)的推導結論:

改寫上式，

假設時變結構中各單元的質量參數為時不變常數，或者其隨時間變化的規律已被經驗掌握，這一假設在工程實踐中不失一般性和合理性，

從式(13)不難發現，線性時變結構的時變阻尼和時變剛度可以通過求解不同時刻線性代數方程組(13)確定。

為了增加識別方法的穩定性和抗噪聲能力，我們將一小段時刻點的線性代數方程組聯立，構成最小二乘問題求解，這大大增加了識別的準確性與抗噪聲能力。這將在算例中有明確的說明。

2.3 線性時變結構瞬時頻率識別

在時變結構的阻尼和剛度確定之后，時變結構的瞬時頻率識別就轉化為普通的特征值問題。

3 仿真算例

3.1 多自由度時變結構

圖1所示為剪切梁樓房模型，各樓層的質量為m1=m2=m3=m4=m5=5 kg，各層之間的剪切剛度為K1=K2=K3=K4=K5=80 000 N/m，各層之間的阻尼系數E1=E2=E3=E4=E5=3 N·s/m。5階固有頻率分別是 5.7 Hz、16.7 Hz、26.4 Hz、33.9 Hz、38.6 Hz。5個自由度的初始位移和初始速度均為零，初始加速度均為5 000 m/s2。結構的自由響應使用Newmark－β法計算，加速度信號的采樣周期Δt=1×10－3s(采樣頻率f=1 000 Hz)，加速度信號的采樣時間總長度為T=2 s。

圖1 五層剪切梁樓房模型Fig.1 A 5 － story shear-beam building model

小波分析函數選取滿足(2)、(3)式的墨西哥草帽小波(Mexican hat wavelet)，其數學表達式為:ψ(t)=(－t4+6t2－ 3)e－t2/2。設其存在區間為 t0＜t＜t1，參數 t0=－8;t1=8。算例中進行小波分析時，可設尺度參數α=αj，平移參數 b= αj(β0k+β1)。其中 α ＞1，β0＞0，j和k均為整數。β0取 t1－t0，它是一個描述連續小波分析長度的參數;β1取 －t0;k=0，1，…，[T/αjβ0]－ 1，符 號[x]表示不大于 x的整數。本算例中 α =20.01，j= －680。在構造最小二乘問題時，1 000個相鄰時刻點的線性代數方程聯立，也就是說1 000個連續時刻點的線性方程(如等式13)聯立求最小二乘解。

為了說明識別方法的適應性，算例中考慮三種不同的時變情況(突變，線性變化和周期變化)。識別過程中物理量的單位:質量單位為kg，剛度單位為N/m，阻尼系數單位為N·s/m，時間為s。

(1)參數突變情況

結構的力學參數隨時間變化關系為:

其他力學參數不隨時間變化。時變結構的五階瞬時頻率識別結果如圖2所示:

圖2 (a) 參數突變結構頻率的識別結果與理論值對比(a)Comparison of true value and the identified instantaneous frequencies for the abruptly time-varying structure

圖2 (b) 參數突變結構頻率的識別結果與理論值對比Fig.2 (b)Comparison of true value and the identified instantaneous frequencies for the abruptly time－varying structure

從圖2(a)和圖2(b)可以看出，在參數突變的瞬間，識別值與理論值之間存在較大的識別誤差。這說明算法在跟蹤參數突變的瞬間，識別能力有限，但算法很快又再次準確地捕捉到了時變結構的瞬時頻率，整體識別效果很好。

(2)參數線性變化情況

結構的力學參數隨時間變化關系為:

其他力學參數不隨時間變化。時變結構的五階瞬時頻率識別結果如圖3所示。

(3)參數周期變化情況

其他力學參數不隨時間變化。時變結構的五階瞬時頻率識別結果如圖4所示:

圖4 (a) 參數周期變化結構頻率的識別結果與理論值對比Fig.4 (a)Comparison of true value and the identified instantaneous frequencies for the periodically time-varying structure

圖4 (b) 參數周期變化結構頻率的識別結果與理論值對比Fig.4 (b)Comparison of true value and the identified instantaneous frequencies for the periodically time-varying structure

從圖3(a)和圖3(b)的識別結果可以看出，時變結構的參數線性變化時，算法可以準確地識別出結構的瞬時頻率，誤差很小。

圖4(a)和圖4(b)中的對比結果顯示，結構參數周期變化時，算法準確地識別出了5階瞬時頻率，但是在識別的起始時段和終止時段的誤差要明顯大于中間時段，這是由于連續小波變換的邊界效應引起的。但這種邊界效應對識別算法的整體識別效果影響不大。

為考察算法的抗噪聲能力，我們給加速度信號中添加了高斯白噪聲，其統計參數設置為零均值，方差為1。定義信噪比(The Signal-To-Noise，SNR)為信號的均方差與噪聲的均方差之比:

考慮在不同信噪比，時變結構參數線性變化和周期變化兩種情況，5階瞬時頻率識別結果可由MAPE值來衡量，見表1和表2。

表1 不同信噪比參數線性變化頻率識別結果的MAPE值Tab.1 Errors of identification(MAPE)using noisy data(different SNR values)with smoothly varying parameters

表2 不同信噪比參數周期變化頻率識別結果的MAPE值Tab.2 Errors of identification(MAPE)using noisy data(different SNR values)with periodically varying parameters

從表1和表2可以看出，隨著高斯白噪聲添加量的增大，即信噪比減小，無論參數線性變化還是周期變化，5階時變頻率的識別精度都較高，且對噪聲不敏感，MAPE都在2%左右，即使在信噪比非常惡劣的情況下(SNR=30)，MAPE仍然不超過4.5%，可見本文識別方法具有良好的抗噪聲能力。

3.2 模態密集時變結構

構造一個3自由度彈簧－質量塊模態密集時變結構，對其進行瞬時頻率識別，考察算法對模態密集結構的適用性。圖5中質量塊m1=1 000 kg，m2=0.2 kg，m3=0.02 kg，剛度 K1=6 ×106N/m，K2=1.2 ×104N/m，K3=120 N/m，阻尼系數 E1=E2=E3=0.01 N·s/m，三階固有頻率12.29、12.34、39.20 Hz，前兩階頻率相差0.05 Hz。結構的初始速度、初始位移均為零，三個自由度上的初始加速度為5 000 m/s。系統的響應使用Newmark－β法求得，其中采樣周期 Δt=0.001 s。小波分析函數及其分析參數的選取同算例一。

為了能夠激發學生學習數學的興趣，在教學中，教師應該根據教學內容合理地創設合作學習的情境，為學生進行小組合作學習提供平臺與條件，也能夠調動學生的積極性，促使學生主動參與小組合作學習的教學。由于小學生具有爭強好勝的心理特點，比賽、競爭類的游戲對學生有相當大的吸引力。在課堂教學中可以適當地引入競賽、比賽等方法來刺激學生的學習興趣。

圖5 三自由度密集模態結構Fig.5 Three degree-of-freedom structure with closely spaced modes

識別過程中力學參數的單位:質量單位為kg，剛度單位為N/m，阻尼系數單位為N·s/m，時間為s。

(1)參數突變情況

結構的力學參數隨時間變化關系為:

其他力學參數不隨時間變化。密集模態時變結構的三階瞬時頻率識別結果如圖6所示:

圖6 參數突變密集模態結構頻率的識別結果與理論值對比Fig.6 Comparison of true value and the identified instantaneous frequencies for the abruptly time-varying structure with closely spaced modes

圖7 參數線性變化密集模態結構頻率的識別結果與理論值對比Fig.7 Comparison of true value and the identified instantaneous frequencies for the smoothly time-varying structure with closely spaced modes

(2)參數線性變化情況

結構的力學參數隨時間變化關系為:

其他力學參數不隨時間變化。密集模態時變結構的三階瞬時頻率識別結果如圖7所示。

從圖6和圖7可以看出，識別算法能夠很有效地識別密集模態結構的瞬時頻率，在前兩階頻率相差不大于0.1 Hz的情況下，無論是結構參數突變還是線性變化，識別算法都能準確地識別出時變結構的瞬時頻率，頻率分辨率很高。

4 結論

(1)本文推導了函數積分運算的連續小波變換計算方法。

(2)借助此法，假設時變結構的質量系數為時不變常數，或者其隨時間變化的規律已被經驗掌握(這一假設在工程實踐中不失一般性和合理性)，僅利用時變結構自由振動的加速度響應信號，計算出速度信號和位移信號的連續小波變換值，進而將結構的時變剛度和時變阻尼的識別問題轉化為求解不同時刻線性代數方程組聯立之后的最小二乘問題。

(3)五層剪切梁樓房時變結構和三自由度模態密集結構的仿真研究，充分驗證了識別算法的正確性、有效性和抗噪聲能力。對于實際工程有較好的應用前景。

(4)本文方法需要時變結構的加速度自由響應數據，在工程實用中有局限性，有待后續的研究改進。

(5)文中最小二乘求解，聯立方程個數的選取是個非線性優化問題，通常情況下，方程聯立越多，識別方法的抗噪聲能力越強，識別精度能夠提高，但這是以犧牲識別效率為代價的，計算速度會變慢。當聯立方程多到一定程度，就會出現過度擬合，此時識別結果反而不好，精度不高。綜上，方程個數的選取需要折中考慮，具體的個數選取與操作者的經驗有關。

致 謝

文中的部分工作是博士生許鑫和龍雙麗在受到中國國家留學基金委資助公派英國后完成的，在此對國家留學基金委表示感謝。

非常感謝University of Patras機械工程與航空系的Prof.S.D.Fassois。和他進行交流與討論是一件令人非常高興的事情。

[1]Feldman M.Non-linear system vibration analysis using Hilbert transform-I:free vibration analysis method FREEV B[J].Mechanical Systems and Signal Processing，1994，8(2):119－127.

[2]Shi Z Y，Law S S，Xu X.Identification of linear time-varying mdof dynamical systems from forced excitation using Hilbert transform and EMD method[J].Journal of Sound and Vibration，2009，321:572－589.

[3]Liu K.Extension of modal analysis to linear time-varying systems[J]. Journal of Sound and Vibration，1999，226(1):149－167.

[4]Shi Z Y，Law S S，Li H N.Subspace-based identification of linear time-varying system[J].AIAA Journal，2007，45:2042－2050.

[5]Staszewski W J.Identification of non-linear systems using multi-scale ridges and skeletons of the wavelet transform[J].Journal of Sound and Vibration，1998，214:639－658.

[6]Ghanem R，Romeo F.A wavelet-based approach for the identification of linear time-varying dynamical systems[J].Journal of Sound and Vibration，2000，234:555－576.

[7]Chen S L，Lai H C，Ho K C.Identification of linear time varying systems by Haar wavelet[J].International Journal of System Science，2006，37:619－628.

[8]Poulimenos A G，Fassois S D.Parametric time-domain methods for non-stationary random vibration modelling and analysis a critical survey and comparison[J].Mechanical Systems and Signal Processing，2006，20:763－816.

[9]Spiridonakos M D，Fassois S D.Parametric identification of a time-varying structure based on vector vibration response measurements[J]. Mechanical Systems and Signal Processing，2009，23:2029－2048.

[10]Smyth A W，Masri S F，Kosmatopoulos E B，et al.Development of adaptive modeling techniques for non-linear hysteretic systems[J].International Journal of Non-Linear Mechanics，2002，37:1435 －1451.

[11]Yang J N，Lin S.Identification of parametric variations of structures based on least squares estimation and adaptive tracking technique[J].Journal of Engineering Mechanics，2005，131:290－298.

[12]Staszewski W J，Cooper J E.Flutter data analysis using the wavelet transform[C].Proceedings of international Congress on MV2:New AdvancesinModalSynthesisofLarge Structures， Non-linear， Damped and Non-Deterministic Cases，Lyon，France，October 1995:549 －561.

[13]Le T P，Argoul P.Continuous wavelet transform for modal identification using free decay response[J].Journal of Sound and Vibration，2004，277:73 －100.

[14]王 超，任偉新，黃天立.基于復小波變換的結構瞬時頻率識別[J].振動工程學報，2009，22(5):492－496.

[15]許 鑫，史治宇.狀態空間下基于小波變換的時變系統參數識別[J].振動工程學報，2010，23(4):415－419.

[16]許 鑫，史治宇.用于時變系統參數識別的狀態空間小波方法[J].工程力學，2011，28(3):23－28.

[17]Sone A，Hata H，Masuda A.Identification of structural parametersusing the wavelettransform ofacceleration measurements[J].Journal of Pressure Vessel Technology，2004，126:128－133.

[18]Xu X，Shi Z Y，You Q.Identification of linear time-varying systems using a wavelet-based state-space method [J].Mechanical Systems and Signal Processing， Sep.2011，published online.