基于矩陣分式多項式時變結構模態參數最小二乘辨識

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

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

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

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

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

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

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

結構動力學系統頻率響應函數可表示為矩陣分式多項式形式：

(1)

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

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

(2)

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

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

(3)

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

(4)

(5)

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

(6)

(7)

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

(8)

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

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

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

2.1 最小二乘問題

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

(9)

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

(10)

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

(11)

其中：

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

(12)

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

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

(13)

其中：

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

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

(14)

(15)

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

(16)

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

(17)

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

D′α′=b′

(18)

其中：

或

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

3 求解模態參數

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

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

(19)

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

(2) 據所得系統極點與參數化功率譜函數獲得系統留數：

(20)

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

Rr=U∑V

(21)

4 算例一

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

4.1 算例設置

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

(22)

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

圖1 三自由度彈簧-阻尼-質量系統

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

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

4.2 辨識結果

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

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

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

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

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

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

5 算例二

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

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

(23)

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

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

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

圖10 辨識所得時變模態參數

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

6 結 論

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

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

參 考 文 獻

[1]鄒經湘, 于開平, 楊炳淵. 時變結構的參數識別方法[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]吳日強, 于開平, 鄒經湘. 改進的子空間方法及其在時變結構參數辨識中的應用[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]龐世偉, 于開平, 鄒經湘. 用于時變結構模態參數識別的投影估計遞推子空間方法[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]楊利芳, 于開平, 龐世偉,等. 用于線性時變結構系統辨識的子空間方法比較研究[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,等.基于加速度響應連續小波變換的線性時變結構瞬時頻率識別[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]續秀忠,張志誼,華宏星,等. 結構時變模態參數辨識的時頻分析方法[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]續秀忠, 華宏星, 張志誼,等. 應用時頻表示進行結構時變模態頻率辨識[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變換的橋梁監測信號分析與處理和時變模態參數識別[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.