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

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

無論是航空航天、高速列車，還是土木橋梁、造船機(jī)械，力學(xué)結(jié)構(gòu)在其服役期內(nèi)總會由于載荷，環(huán)境，人為等因素的影響出現(xiàn)結(jié)構(gòu)的力學(xué)參數(shù)隨時間變化的特性。這種時變特性所引起的力學(xué)問題近一二十年來備受關(guān)注。而結(jié)構(gòu)參數(shù)識別是力學(xué)建模和力學(xué)分析的基礎(chǔ)，時變結(jié)構(gòu)參數(shù)識別也由此成為近些年來參數(shù)識別領(lǐng)域的一個熱點問題。

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

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

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

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

1 函數(shù)積分運(yùn)算的連續(xù)小波變換

1.1 連續(xù)小波變換

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

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

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

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

1.2 函數(shù)積分運(yùn)算的連續(xù)小波變換計算

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

其中y0為積分常數(shù)，對上式用分部積分公式，

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

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

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

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

2 線性時變結(jié)構(gòu)的瞬時頻率識別

2.1 線性時變結(jié)構(gòu)自由振動微分方程

一個p自由度線性時變動力學(xué)結(jié)構(gòu)，其自由振動可以用下列方程式表示:

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

2.2 時變結(jié)構(gòu)物理參數(shù)確定

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

改寫上式，

假設(shè)時變結(jié)構(gòu)中各單元的質(zhì)量參數(shù)為時不變常數(shù)，或者其隨時間變化的規(guī)律已被經(jīng)驗掌握，這一假設(shè)在工程實踐中不失一般性和合理性，

從式(13)不難發(fā)現(xiàn)，線性時變結(jié)構(gòu)的時變阻尼和時變剛度可以通過求解不同時刻線性代數(shù)方程組(13)確定。

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

2.3 線性時變結(jié)構(gòu)瞬時頻率識別

在時變結(jié)構(gòu)的阻尼和剛度確定之后，時變結(jié)構(gòu)的瞬時頻率識別就轉(zhuǎn)化為普通的特征值問題。

3 仿真算例

3.1 多自由度時變結(jié)構(gòu)

圖1所示為剪切梁樓房模型，各樓層的質(zhì)量為m1=m2=m3=m4=m5=5 kg，各層之間的剪切剛度為K1=K2=K3=K4=K5=80 000 N/m，各層之間的阻尼系數(shù)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。結(jié)構(gòu)的自由響應(yīng)使用Newmark－β法計算，加速度信號的采樣周期Δt=1×10－3s(采樣頻率f=1 000 Hz)，加速度信號的采樣時間總長度為T=2 s。

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

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

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

(1)參數(shù)突變情況

結(jié)構(gòu)的力學(xué)參數(shù)隨時間變化關(guān)系為:

其他力學(xué)參數(shù)不隨時間變化。時變結(jié)構(gòu)的五階瞬時頻率識別結(jié)果如圖2所示:

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

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

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

(2)參數(shù)線性變化情況

結(jié)構(gòu)的力學(xué)參數(shù)隨時間變化關(guān)系為:

其他力學(xué)參數(shù)不隨時間變化。時變結(jié)構(gòu)的五階瞬時頻率識別結(jié)果如圖3所示。

(3)參數(shù)周期變化情況

其他力學(xué)參數(shù)不隨時間變化。時變結(jié)構(gòu)的五階瞬時頻率識別結(jié)果如圖4所示:

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

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

從圖3(a)和圖3(b)的識別結(jié)果可以看出，時變結(jié)構(gòu)的參數(shù)線性變化時，算法可以準(zhǔn)確地識別出結(jié)構(gòu)的瞬時頻率，誤差很小。

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

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

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

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

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

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

3.2 模態(tài)密集時變結(jié)構(gòu)

構(gòu)造一個3自由度彈簧－質(zhì)量塊模態(tài)密集時變結(jié)構(gòu)，對其進(jìn)行瞬時頻率識別，考察算法對模態(tài)密集結(jié)構(gòu)的適用性。圖5中質(zhì)量塊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，阻尼系數(shù) E1=E2=E3=0.01 N·s/m，三階固有頻率12.29、12.34、39.20 Hz，前兩階頻率相差0.05 Hz。結(jié)構(gòu)的初始速度、初始位移均為零，三個自由度上的初始加速度為5 000 m/s。系統(tǒng)的響應(yīng)使用Newmark－β法求得，其中采樣周期 Δt=0.001 s。小波分析函數(shù)及其分析參數(shù)的選取同算例一。

為了能夠激發(fā)學(xué)生學(xué)習(xí)數(shù)學(xué)的興趣，在教學(xué)中，教師應(yīng)該根據(jù)教學(xué)內(nèi)容合理地創(chuàng)設(shè)合作學(xué)習(xí)的情境，為學(xué)生進(jìn)行小組合作學(xué)習(xí)提供平臺與條件，也能夠調(diào)動學(xué)生的積極性，促使學(xué)生主動參與小組合作學(xué)習(xí)的教學(xué)。由于小學(xué)生具有爭強(qiáng)好勝的心理特點，比賽、競爭類的游戲?qū)W(xué)生有相當(dāng)大的吸引力。在課堂教學(xué)中可以適當(dāng)?shù)匾敫傎?、比賽等方法來刺激學(xué)生的學(xué)習(xí)興趣。

圖5 三自由度密集模態(tài)結(jié)構(gòu)Fig.5 Three degree-of-freedom structure with closely spaced modes

識別過程中力學(xué)參數(shù)的單位:質(zhì)量單位為kg，剛度單位為N/m，阻尼系數(shù)單位為N·s/m，時間為s。

(1)參數(shù)突變情況

結(jié)構(gòu)的力學(xué)參數(shù)隨時間變化關(guān)系為:

其他力學(xué)參數(shù)不隨時間變化。密集模態(tài)時變結(jié)構(gòu)的三階瞬時頻率識別結(jié)果如圖6所示:

圖6 參數(shù)突變密集模態(tài)結(jié)構(gòu)頻率的識別結(jié)果與理論值對比Fig.6 Comparison of true value and the identified instantaneous frequencies for the abruptly time-varying structure with closely spaced modes

圖7 參數(shù)線性變化密集模態(tài)結(jié)構(gòu)頻率的識別結(jié)果與理論值對比Fig.7 Comparison of true value and the identified instantaneous frequencies for the smoothly time-varying structure with closely spaced modes

(2)參數(shù)線性變化情況

結(jié)構(gòu)的力學(xué)參數(shù)隨時間變化關(guān)系為:

其他力學(xué)參數(shù)不隨時間變化。密集模態(tài)時變結(jié)構(gòu)的三階瞬時頻率識別結(jié)果如圖7所示。

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

4 結(jié)論

(1)本文推導(dǎo)了函數(shù)積分運(yùn)算的連續(xù)小波變換計算方法。

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

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

(4)本文方法需要時變結(jié)構(gòu)的加速度自由響應(yīng)數(shù)據(jù)，在工程實用中有局限性，有待后續(xù)的研究改進(jìn)。

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

致 謝

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

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

[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，F(xiàn)assois 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，F(xiàn)assois 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，F(xiàn)rance，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]王 超，任偉新，黃天立.基于復(fù)小波變換的結(jié)構(gòu)瞬時頻率識別[J].振動工程學(xué)報，2009，22(5):492－496.

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

[16]許 鑫，史治宇.用于時變系統(tǒng)參數(shù)識別的狀態(tài)空間小波方法[J].工程力學(xué)，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.