基于矩陣分式多項式時變結構模態(tài)參數(shù)最小二乘辨識

周思達， 劉 莉， 楊 武， 馬志賽

(北京理工大學 宇航學院，北京 100081)

由于傳統(tǒng)結構動力學研究中，通常假設被分析結構為時不變，故使結構動力學研究的諸多成果成功運用于工程結構分析、設計中。然而隨工程結構應用領域的不斷拓展，尤其航空航天、機械、橋梁等應用領域，越來越多結構隨時間變化特性不可忽略[1]。如鐵路橋與火車形成的組合結構系統(tǒng)，火車通過時的運動導致組合結構系統(tǒng)質量分布會發(fā)生改變[2-3]，甚至導致結構阻尼、剛度的改變[4]；制造業(yè)中機械臂在完成不同動作時幾何特征會發(fā)生變化，導致其結構動力學特性產生相應變化[5-6]；人與運動場形成的組合系統(tǒng)，運動場處于不同狀態(tài)時，無比賽時場地空置、比賽前觀眾入場、比賽進行等，人與運動場組合系統(tǒng)的結構特性隨時間變化，此在運動場結構健康監(jiān)測中亦不可忽略[7]。

模態(tài)分析中模態(tài)參數(shù)辨識已廣泛應用于結構動力學研究與應用。如結構健康監(jiān)測、故障診斷、模型更新、強迫振動分析等[8]。因此，進行線性時變結構的模態(tài)參數(shù)辨識研究十分必要。

在線性時變系統(tǒng)理論及系統(tǒng)辨識理論研究發(fā)展推動下，線性時變結構模態(tài)參數(shù)辨識研究亦在發(fā)展?，F(xiàn)有時變結構模態(tài)參數(shù)辨識方法主要分兩類，即參數(shù)化時域辨識方法與非參數(shù)化時頻域辨識方法。參數(shù)化時域辨識方法中亦分兩類：①基于時間相關自回歸滑動平均模型(TARMA)，如Poulimenos等[9]的泛函序列TARMA(FS-TARMA)方法；②基于狀態(tài)空間模型，如Liu[10]將傳統(tǒng)子空間方法(SSI)拓展到時變結構。吳日強等[11-12]提出改進遞推子空間方法；楊利芳等[13]對幾種用于時變結構模態(tài)參數(shù)辨識子空間方法進行比較研究。非參數(shù)化時頻域辨識方法為基于時頻分析的時變結構模態(tài)參數(shù)辨識方法，通過對時變結構響應信號進行時頻分析識別出目標結構相關模態(tài)參數(shù)，如Ghane等[14]提出基于小波變換的模態(tài)參數(shù)辨識方法；許鑫等[15]提出基于連續(xù)小波變化的線性時變結構瞬時頻率辨識方法；Roshan-Ghias等[16]提出基于平滑偽Wigner-Ville分布的辨識方法；續(xù)秀忠等[17-18]利用Gabor展開及Wigner-Ville分布辨識時變結構模態(tài)頻率；王學敏[19]提出基于Hilbert-Huang變換(HHT)的線性時變結構模態(tài)頻率辨識方法；Feldman[20]提出基于Hilbert振動分解(HVD)的時變結構模態(tài)參數(shù)辨識方法。

本文在現(xiàn)有針對時不變結構參數(shù)化最小二乘復頻域方法基礎上，將其拓展到時頻域，提出參數(shù)化時頻域模態(tài)參數(shù)辨識方法，其中參數(shù)化模型為時間相關的矩陣分式多項式模型，參數(shù)估計方法為最小二乘方法，并通過仿真算例驗證所提方法的有效性。本文時變結構模態(tài)參數(shù)最小二乘辨識方法由三部分組成：①建立時間相關矩陣分式多項式形式的時變結構動力學參數(shù)化模型，包括模型階數(shù)等，即時間與頻率多項式階數(shù)；②構建時頻域最小二乘參數(shù)估計方法，估計出時變結構參數(shù)化動力學模型中待估參數(shù)；③據(jù)所得待估參數(shù)，求解時變結構模態(tài)參數(shù)。

1 時間相關矩陣分式多項式模型

結構動力學系統(tǒng)頻率響應函數(shù)可表示為矩陣分式多項式形式：

(1)

其中：B(ω)為分子矩陣多項式；A(ω)為公分母多項式。

GXX(ω)=H(ω)GFF(ω)GH(ω)

(2)

其中：H(ω)為結構頻率響應函數(shù)矩陣；GFF(ω)為激勵(輸入)信號功率譜函數(shù)矩陣；GXX(ω)為結構響應(輸出)信號功率譜函數(shù)矩陣。

對僅輸出模態(tài)參數(shù)辨識而言，設作用于結構的激勵為白噪聲，結構響應功率譜函數(shù)輸入信號功率譜函數(shù)矩陣GFF(ω)為關于頻率的常數(shù)矩陣。則結構響應信號功率譜函數(shù)僅為結構參數(shù)(如剛度、阻尼系數(shù)、質量)的函數(shù)，也可寫為式(1)的矩陣分式形式。據(jù)Zadeh[21]關于線性時變系統(tǒng)頻域表示，慢時變系統(tǒng)仍可近似表示為分式多項式形式。因此，本文將時不變結構功率譜函數(shù)矩陣分式多項式模型拓展為時間相關的矩陣分式多項式模型，以表示時變結構時間相關功率譜函數(shù)(輸出o及參考點i)，定義為：

(3)

其中：下標τ=1,2,…,N為時間采樣；f=1,2,…,Nf為頻率采樣；k=1,2,…,NoNi(即k=(i-1)No+o)，No與Ni分別為輸出數(shù)及參考點數(shù)，分子多項式與公分母多項式分別定義為：

(4)

(5)

其中：φi,j(t,ω)為時頻基函數(shù)(時間多項式及頻率多項式分別為i，j階，i=0,1,m2,…,nt，j=0,1,2,…,nω)；多項式系數(shù)bk,i,j及ai,j可寫為向量形式：

(6)

(7)

為方便推導，將式(7)參數(shù)向量寫為：

(8)

其中：θ為參數(shù)化矩陣分式多項式模型待定參數(shù)向量。

2 時頻域最小二乘辨識方法

提出時頻域最小二乘方法，對式(3)線性時變結構參數(shù)化模型進行參數(shù)估計，獲得待定參數(shù)向量θ，并通過θ獲得結構時變模態(tài)參數(shù)。

2.1 最小二乘問題

據(jù)時頻域時變結構動力學參數(shù)化時間相關矩陣分式多項式模型(式(3))，時頻域線性最小二乘費用函數(shù)定義為：

(9)

其中：θ為式(8)的待定參數(shù)向量，方程誤差為：

(10)

由式(4)～式(8)，方程誤差可寫成矩陣形式，為：

(11)

其中：

據(jù)矩陣形式定義方程誤差(式(11))，時頻域線性最小二乘費用函數(shù)(式(9))可寫為：

(12)

考慮最小二乘中Jacobi矩陣分塊矩陣，式(12)可表示為：

?LS(θ)=θTRe(JHJ)θ

(13)

其中：

2.2 基于減縮正則方程最小二乘求解

為求解式(13)最小二乘問題，采用基于縮減正則方程方法求解。當式(9)費用函數(shù)取最小值時，其對待定參數(shù)導數(shù)為零，即：

(14)

(15)

將式(14)代入式(15)，得：

(16)

式(13)參數(shù)化模型具有參數(shù)冗余，給定常數(shù)γ，有：

(17)

式(17)的參數(shù)冗余會使縮減的最小二乘正則方程右端項為0，導致方程只有平凡解或無窮多解。因此，須對模型參數(shù)施加約束，選α中任意一分量設為常數(shù)1，則縮減的正則方程為：

D′α′=b′

(18)

其中：

或

對α中不同分量施加約束，會引起辨識結果包含不同情況的穩(wěn)定極點與不穩(wěn)定極點，文獻[22]給出關于時不變結構相關辨識的詳細分析與結論。本文將利用解析算例對該內容進行詳細討論。

3 求解模態(tài)參數(shù)

通過線性最小二乘辨識可得參數(shù)化模型(式(3))中待估的未知參數(shù)向量θ。若θ確定，則時間相關功率譜函數(shù)矩陣分式多項式模型確定。將矩陣分式多項式模型轉換為極點-留數(shù)模型可得結構模態(tài)參數(shù)：

(1) 通過求解含系數(shù)α的各時刻公分母多項式A(t,ω,α)的根獲得系統(tǒng)極點，再利用系統(tǒng)極點獲得結構模態(tài)頻率及模態(tài)阻尼比：

(19)

其中：fr(t)為t時刻結構模態(tài)頻率；ξr(t)為結構模態(tài)阻尼比。

(2) 據(jù)所得系統(tǒng)極點與參數(shù)化功率譜函數(shù)獲得系統(tǒng)留數(shù)：

(20)

(3) 通過留數(shù)矩陣計算工作模態(tài)參數(shù)向量(無縮放的結構模態(tài)振型)[23]：對留數(shù)矩陣進行奇異值分解(SVD)獲得酉矩陣U第一列即工作模態(tài)參數(shù)向量：

Rr=U∑V

(21)

4 算例一

為考察本文所提線性時變結構模態(tài)參數(shù)辨識方法性能與待估參數(shù)約束對辨識結果影響，給出解析表達功率譜函數(shù)算例。

4.1 算例設置

三自由度彈簧-阻尼-質量系統(tǒng)見圖1。其中質量隨時間變化，質量矩陣為：

(22)

其中：tend=1 s為終止時刻；pc=0.5為比例系數(shù)；M0為初始時刻質量矩陣。

圖1 三自由度彈簧-阻尼-質量系統(tǒng)

設激勵為理想白噪聲，系統(tǒng)響應功率譜函數(shù)僅為結構參數(shù)函數(shù)。據(jù)給定系統(tǒng)參數(shù)可直接獲得此系統(tǒng)解析的無噪聲時間相關功率譜函數(shù)，第一、二自由度響應互功率譜函數(shù)見圖2(限于篇幅，僅給出G12)。

圖2 解析的無噪聲時間相關功率譜函數(shù)(G12)

4.2 辨識結果

以圖1的解析時間相關功率譜函數(shù)(包含G11，G12，G13)為測量值，采用所提線性時變結構模態(tài)參數(shù)最小二乘辨識方法進行模態(tài)參數(shù)辨識，獲得辨識的時間相關功率譜函數(shù)，其中第一、二自由度響應互功率譜函數(shù)見圖3。辨識所得時變模態(tài)頻率及模態(tài)阻尼比見圖4(圓圈為理論值，實線為辨識結果)。

由圖3、圖4看出，所提線性時變結構模態(tài)參數(shù)最小二乘辨識方法能十分準確的辨識出無噪聲功率譜情況下時變結構模態(tài)參數(shù)，且辨識結果不依賴待估參數(shù)的約束。圖1的解析時間相關無噪聲功率譜函數(shù)上疊加高斯白噪聲，獲得信噪比10 dB有噪聲情況下時間相關功率譜。0.25 s時刻有噪聲時間相關功率譜見圖5(其中實線為無噪聲，點為有噪聲)。

圖5 解析的有噪聲時間相關功率譜函數(shù)(有噪聲，G12)

圖6 辨識所得時間相關功率譜函數(shù)(有噪聲，G^12)

以上所述有噪聲時間相關功率譜為測量值，采用所提辨識方法進行模態(tài)參數(shù)辨識。分別施加式(18)兩種待估參數(shù)約束，辨識獲得時間相關功率譜函數(shù)見圖6。辨識所得模態(tài)頻率及阻尼比見圖7。

由以上辨識結果可知，有噪聲時α的不同約束方式對辨識結果影響較大。原因為對α的不同約束對辨識所得穩(wěn)定極點、不穩(wěn)定極點的分布產生影響，即在約束時不穩(wěn)定極點被視作穩(wěn)定極點。而無噪聲時，穩(wěn)定極點與不穩(wěn)定極點成對出現(xiàn)，即將不穩(wěn)定極點視為穩(wěn)定極點亦不影響結果。所提方法具有的此性質與用于時不變結構模態(tài)參數(shù)辨識的最小二乘復頻域法相關性質[22]類似。若給定正確約束(式(18)的約束二)，所提線性時變結構模態(tài)參數(shù)最小二乘辨識方法能較好辨識出有噪聲時線性時變結構模態(tài)頻率及模態(tài)阻尼比。

5 算例二

實際工程應用中結構響應時間相關功率譜函數(shù)無法直接測量或通過解析方式獲得。因此通過數(shù)值積分求解響應仿真算例進一步驗證所提線性時變結構模態(tài)參數(shù)最小二乘辨識方法。由于激勵是隨機的且響應由數(shù)值積分直接求得，因此本算例考慮輸入噪聲。所用時變系統(tǒng)與算例一類似，仍為三自由度彈簧-阻尼-質量時變系統(tǒng)，質量隨時間按線性關系減少，即：

M(t)=M0(1-pct/tend)

(23)

其中：tend=1 s為終止時刻；pc=0.5為比例系數(shù)。

給定激勵為作用在系統(tǒng)第一自由度的仿真高斯白噪聲。系統(tǒng)響應采用Newmark-β數(shù)值積分方法[24]計算獲得，其中積分步長為1/4 096 s。響應信號采樣頻率1 024 Hz，采樣時間1 s。采用SPWVD非參數(shù)化估計時間相關功率譜函數(shù)，平均次數(shù)40次。由平滑偽Wigner-Ville分布(SPWVD)[25]非參數(shù)化估計的第一、二、三自由度響應時間相關互功率譜函數(shù)見圖8。

圖8 SPWVD非參數(shù)化估計時間相關功率譜函數(shù)(G12，G13)

圖10 辨識所得時變模態(tài)參數(shù)

以上仿真算例表明，所提線性時變結構模態(tài)參數(shù)最小二乘辨識方法能在僅輸出情況下較好辨識出線性時變結構的模態(tài)頻率，但由于阻尼辨識本身困難及最小二乘辨識未考慮測量值的不確定性，阻尼辨識結果不及模態(tài)頻率好，尤其第一、三階模態(tài)阻尼比與理論值有一定偏差。

6 結 論

(1) 本文所提基于矩陣分式多項式模型的線性時變結構模態(tài)參數(shù)最小二乘辨識方法，給出線性時變結構時頻域參數(shù)化模型-時間相關矩陣分式多項式模型；所提模態(tài)參數(shù)最小二乘辨識方法，其中基于縮減正則方程的最小二乘求解方法大大降低對計算資源要求。

(2) 通過解析算例驗證所提線性時變結構模態(tài)參數(shù)最小二乘辨識方法中待估參數(shù)約束對辨識結果影響，指出正確約束方式。并在僅輸出條件下較好辨識出線性時變結構模態(tài)頻率，方法的有效性獲得驗證。

參 考 文 獻

[1]鄒經湘, 于開平, 楊炳淵. 時變結構的參數(shù)識別方法[J]. 力學進展, 2000, 30(3): 370-377.

ZOU Jing-xiang,YU Kai-ping,YANG Bing-yuan. Methods of time-varying structural parameter identification[J]. Advances In Mechanics, 2000, 30(3): 370-377.

[2]Au F T K, Jiang R J, Cheung Y K. Parameter identification of vehicles moving on continuous bridges[J]. Journal of Sound and Vibration, 2004, 269(1-2): 91-111.

[3]Marchesiello S, Bedaoui S, Garibaldi L, et al. Time-dependent identification of a bridge-like structure with crossing loads[J]. Mechanical Systems and Signal Processing,2009,23(6):2019- 2028.

[4]Jiang R J, Au F T K, Cheung Y K. Identification of vehicles moving on continuous bridges with rough surface [J]. Journal of Sound and Vibration, 2004, 274(3-5): 1045-1063.

[5]Kopmaz O, Anderson K S. On the eigenfrequencies of a flexible arm driven by a flexible shaft[J]. Journal of Sound and Vibration, 2001, 240(4): 679-704.

[6]Chang P H, Park H S. Time-varying input shaping technique applied to vibration reduction of an industrial robot[J]. Control Engineering Practice, 2005, 13(1): 121-130.

[7]Peeters B, Van Der Auweraer H, Vanhollebeke F, et al. Operational modal analysis for estimating the dynamic properties of a stadium structure during a football game[J]. Shock and Vibration, 2007, 14(4): 283-303.

[8]Heylen W, Lammens S, Sas P. Modal analysis theory and testing [M]. Leuven: Katholieke Universiteit Leuven, 2007.

[9]Poulimenos A G, Fassois S D. Output-only stochastic identification of a time-varying structure via functional series TARMA models[J]. Mechanical Systems and Signal Processing, 2009, 23(4): 1180-1204.

[10]LIU K. Identification of linear time-varying systems[J]. Journal of Sound and Vibration, 1997, 206(4): 487-505.

[11]吳日強, 于開平, 鄒經湘. 改進的子空間方法及其在時變結構參數(shù)辨識中的應用[J]. 工程力學,2002,19(4): 67-70,89.

WU Ri-qiang, YU Kai-ping, ZOU Jing-xiang. An improved subspace method and its application to parameter identification of time-varying structures[J]. Engineering Mechanics, 2002, 19(4): 67-70,89.

[12]龐世偉, 于開平, 鄒經湘. 用于時變結構模態(tài)參數(shù)識別的投影估計遞推子空間方法[J].工程力學,2005,22(5):115-119.

PANG Shi-wei,YU Kai-ping,ZOU Jing-xiang. A projection approximation recursive subspace method for identification of modal parameters of time-varying structures[J]. Engineering Mechanics, 2005, 22(5): 115-119.

[13]楊利芳, 于開平, 龐世偉,等. 用于線性時變結構系統(tǒng)辨識的子空間方法比較研究[J]. 振動與沖擊,2007,26(3):8-12, 153.

YANG Li-fang,YU Kai-ping, PANG Shi-wei, et al. Comparison study on identification methods applied to linear time-varying structures[J]. Journal of Vibration and Shock, 2007, 26(3): 8-12,153.

[14]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(4): 555-576.

[15]許 鑫,史治宇,Wieslaw J. Staszewski,等.基于加速度響應連續(xù)小波變換的線性時變結構瞬時頻率識別[J].振動與沖擊, 2012,31(20):166-171.

XU Xin,SHI Zhi-yu, Staszewski W J, et al.Instantaneous frequnencies identification of a linear time-varying structure using continuous wavelet transformation of free decay acceleration response[J].Journal of Vibration and Shock, 2012, 31(20):166-171.

[16]Roshan-ghias A, Shamsollahi M B, Mobed M, et al. Estimation of modal parameters using bilinear joint time-frequency distributions[J]. Mechanical Systems and Signal Processing, 2007, 21(5): 2125-2136.

[17]續(xù)秀忠,張志誼,華宏星,等. 結構時變模態(tài)參數(shù)辨識的時頻分析方法[J]. 上海交通大學學報, 2003, 37(1): 122-126.

XU Xiu-zhong, ZHANG Zhi-yi, HUA Hong-xing, et al. Time-varying modal parameter identification with time-frequency analysis methods[J]. Journal of Shanghai Jiaotong University, 2003, 37(1): 122-126.

[18]續(xù)秀忠, 華宏星, 張志誼,等. 應用時頻表示進行結構時變模態(tài)頻率辨識[J]. 振動與沖擊, 2002, 21(2): 36-40,44.

XU Xiu-zhong, HUA Hong-xing, ZHANG Zhi-yi, et al. Time-varying modal frequency identification by using time-frequency representation [J]. Journal of Vibration and Shock, 2002, 21(2): 36-40,44.

[19]王學敏. 基于Hilbert-Huang變換的橋梁監(jiān)測信號分析與處理和時變模態(tài)參數(shù)識別[D]. 長沙:中南大學, 2008.

[20]Feldman M. Time-varying vibration decomposition and analysis based on the Hilbert transform [J]. Journal of Sound and Vibration, 2006, 295(3-5): 518-530.

[21]Zadeh L A. Frequency analysis of variable networks[J]. Proceedings of the IRE, 1950, 38(3): 291-299.

[22]Cauberghe B, Guillaume P, Verboven P, et al. On the influence of the parameter constraint on the stability of the poles and the discrimination capabilities of the stabilisation diagrams[J]. Mechanical Systems and Signal Processing,2005,19(5): 989-1014.

[23]Verboven P. Frequency-domain system identification for modal analysis [D]. Brussels:Vrije Universiteit Brussel, 2002.

[24]Newmark N M. A method of computation for structural dynamics [J]. Journal of Engineering Mechanics, ASCE, 1959, 85(7): 67-94.

[25]Qian S. Introduction to time-frequency and wavelet transforms [M]. Prentice Hall, 2002.