辛峻峰，王樹青，劉福順

(中國海洋大學(xué) 工程學(xué)院，山東 青島 266100)

有效的模態(tài)識別方法能精確得到結(jié)構(gòu)的相關(guān)參數(shù)，從而可準確掌握大型結(jié)構(gòu)的健康狀況。Wang等［1-3］提出基于時域響應(yīng)的模態(tài)參數(shù)識別方法。如時間序列法、隨機減量法、自然激勵技術(shù)、隨機子空間法等。其中隨機子空間法是目前較先進的環(huán)境激勵下模態(tài)參數(shù)識別方法之一，能精確提取海洋平臺等大型結(jié)構(gòu)的模態(tài)參數(shù)。

DeMoor等［4］證明了數(shù)據(jù)驅(qū)動隨機子空間算法與協(xié)方差驅(qū)動隨機子空間算法在理論上的完全一致。Peters等［5］認為在應(yīng)用過程中，兩種方法的模態(tài)識別結(jié)果也一致。但在對導(dǎo)管架平臺模態(tài)分析過程中，兩種方法的識別結(jié)果并不一致。對此，本文研究探討兩種方法產(chǎn)生差異的原因，繼而進行對比分析。

1 隨機子空間法

簡單介紹隨機子空間法要點，討論詳見文獻［4-6］。在僅考慮隨機噪聲前提下，振動系統(tǒng)的離散狀態(tài)空間方程可表示為:

式中:xk∈Rn為離散時間狀態(tài)向量;A∈Rn×n為離散狀態(tài)矩陣;C∈Rl×n為輸出矩陣，描述內(nèi)部狀態(tài)轉(zhuǎn)化為外界測量值過程;wk∈Rn為由干擾、模型誤差造成的過程噪聲;vk∈Rl為由傳感器誤差等造成的測量噪聲。兩種噪聲均不可測量，在推導(dǎo)過程中常被假設(shè)為零均值，平穩(wěn)的白噪聲。兩種噪聲的協(xié)方差矩陣可表示為:

式中:E為數(shù)學(xué)期望算子;Q∈Rn×n，S∈Rn×l，X∈Rl×l;δpq表示 δpq=1(p=q)，δpq=0(p≠q)，且E［wp］=0，E［vp］=0。

兩種隨機子空間算法的本質(zhì)區(qū)別在于“過去”與“將來”相關(guān)性的計算方式不同:數(shù)據(jù)驅(qū)動隨機子空間法用投影計算:

式中:i為采樣時間間隔;()+表示偽逆(pseudoinverse)。而協(xié)方差驅(qū)動隨機子空間法用協(xié)方差矩陣計算:

由此數(shù)據(jù)驅(qū)動與協(xié)方差驅(qū)動隨機子空間法的關(guān)系主要為投影與協(xié)方差矩陣之間的關(guān)系。

2 數(shù)據(jù)驅(qū)動隨機子空間法與協(xié)方差驅(qū)動隨機子空間法之關(guān)系

推導(dǎo)［4］，投影Pi與協(xié)方差矩陣Ci為相似矩陣。據(jù)相似矩陣性質(zhì)知，二者奇異值相同，從而可得相同模態(tài)參數(shù)。與文獻［5］結(jié)果相同。但本文研究發(fā)現(xiàn)，被忽略的2個細節(jié)能使兩種方法在應(yīng)用中產(chǎn)生差異。

2.1 數(shù)據(jù)有限性

理論上式(4)、(5)成立的前提為:數(shù)據(jù)是無限的。而實際使用中由于計算時間、計算精度等因素影響，所提數(shù)據(jù)必然有限，故式(4)、(5)的右半部分(Yf/Yp與YfYTp)分別變成對投影Pi與協(xié)方差矩陣Ci估計，因數(shù)值計算方式不同，必會出現(xiàn)不同的估計誤差，故在應(yīng)用中投影Pi與協(xié)方差Ci不再相似，因而兩種方法的識別結(jié)果也會不一致。即有限的數(shù)據(jù)可能是造成差異的原因之一。

2.2 QR分解

除數(shù)據(jù)影響外，研究發(fā)現(xiàn)在數(shù)據(jù)驅(qū)動隨機子空間法的計算過程中，QR算法的使用也可能是差異產(chǎn)生的另一主要因素。由于常用的數(shù)據(jù)驅(qū)動隨機子空間法并非直接用式(4)計算投影，而在計算投影的過程中用QR分解算法:

因此數(shù)據(jù)驅(qū)動隨機子空間法應(yīng)準確為:基于QR分解的數(shù)據(jù)驅(qū)動隨機子空間法。本文由此推測QR算法的使用也會影響識別結(jié)果，為造成差異的另一原因。

3 數(shù)值算例

用一系列數(shù)值試驗證明差異的存在，繼而驗證差異原因，并對兩種隨機子空間法進行對比。建立五自由度質(zhì)量-彈簧-阻尼系統(tǒng)的數(shù)值模型［7］，見圖1。

圖1 自由度彈簧-質(zhì)量-阻尼系統(tǒng)Fig.1 A 5-DOF mass-spring-dashpot system

每個單元質(zhì)量、剛度、阻尼系數(shù)分別為:mn=50 kg，kn=2.9 ×107N/m，cn=1 000 N·s/m。位移記為xn，(n=1，…，5)。通過特征值分析，獲得五階模態(tài)頻率的理論值為:34.499 Hz，100.700 Hz，158.730 Hz，203.880 Hz，232.520 Hz;五階模態(tài)阻尼比的理論值為:0.003 737，0.010 909，0.017 197，0.022 092，0.025 198。

3.1 白噪聲激勵實驗

采用Matlab中l(wèi)sim函數(shù)［7］，以均值為零、標準差為1.006的高斯白噪聲為輸入，將其右向加載在靠固定端的首個質(zhì)量塊m1上(圖1)。采樣頻率500 Hz，從首個質(zhì)量塊采集9 000個數(shù)據(jù)點，約18 s，提取的響應(yīng)數(shù)據(jù)見圖2。

圖2 白噪聲激勵下首個質(zhì)量塊輸出信號Fig.2 Response signal of the first mass with white noise loading

本文構(gòu)建Hankel矩陣(i=100)，作為所有計算過程使用數(shù)據(jù)。

3.2 差異現(xiàn)象驗證

為驗證兩種方法差異的存在，分別用基于QR的數(shù)據(jù)驅(qū)動隨機子空間法與協(xié)方差驅(qū)動隨機子空間法分析測量數(shù)據(jù)。分析結(jié)果如圖3所示。

圖3 兩種方法穩(wěn)定圖(*為穩(wěn)定，○為不穩(wěn)定)Fig.3 Two methods of stabilization diagram

圖4 兩種方法穩(wěn)定圖(*為穩(wěn)定，○為不穩(wěn)定)Fig.4 Two methods of stabilization diagram

圖5 兩種方法穩(wěn)定圖(*為穩(wěn)定，○為不穩(wěn)定)Fig.5 Two methods of stabilization diagram

本文穩(wěn)定圖用信號的傅里葉變換為背景，顯示模態(tài)階數(shù)在0～30之間模態(tài)識別結(jié)果情況。對同一模態(tài)，若前后兩次識別結(jié)果同一模態(tài)頻率誤差為1%，同時響應(yīng)阻尼識別結(jié)果誤差為5%，即:

則此次識別結(jié)果標識為穩(wěn)定，否則為不穩(wěn)定［8］。

由圖3看出，兩種方法的識別結(jié)果并不一致，與基于QR分解的數(shù)據(jù)驅(qū)動隨機子空間法(圖3(b))相比，協(xié)方差驅(qū)動隨機子空間法(圖3(a))特點如下:

(1)在橫軸100 Hz、150 Hz左右，縱軸10～25之間形成穩(wěn)定的非常接近的虛假模態(tài)(放大圖);

(2)需高的階數(shù)才能識別較弱勢的第4個模態(tài)(區(qū)間在200 Hz左右);

(3)未能識別出較弱的第5個模態(tài)(區(qū)間在230 Hz左右)。

以上分析表明:差異現(xiàn)象是存在的。以下將展開對造成差異原因的驗證。

3.3 數(shù)據(jù)有限性影響驗證

本節(jié)將驗證數(shù)據(jù)有限性影響數(shù)據(jù)驅(qū)動隨機子空間法的識別結(jié)果。

未用QR分解的數(shù)據(jù)驅(qū)動隨機子空間法與協(xié)方差驅(qū)動隨機子空間方法唯一區(qū)別在于:有限數(shù)據(jù)造成的投影與協(xié)方差矩陣的估計誤差。若該兩種方法的分析結(jié)果不一致，則表示數(shù)據(jù)有限性造成了差異。穩(wěn)定圖的分析結(jié)果見圖4。

由圖4看出，兩種方法的模態(tài)識別結(jié)果不一致:圖4(a)中僅識別出第三階模態(tài)(區(qū)間在150左右)，而圖4(b)中卻識別出4個模態(tài)，由此表明:數(shù)據(jù)有限性導(dǎo)致了估計誤差的產(chǎn)生，是造成差異的原因之一。

3.4 QR分解算法影響驗證

本節(jié)將驗證QR分解算法影響了數(shù)據(jù)驅(qū)動隨機子空間法的識別結(jié)果。

基于QR分解與無QR分解的數(shù)據(jù)驅(qū)動隨機子空間法之間唯一區(qū)別在于:是否使用QR分解技術(shù)。若該兩種方法的模態(tài)識別結(jié)果不一致，即驗證QR分解算法會造成差異。

由圖5看出，兩種方法的識別結(jié)果差別較大:圖5(a)中僅識別出第三階模態(tài)(區(qū)間150左右)，圖5(b)中識別出全部五階模態(tài)。由此證明，QR分解提高了數(shù)據(jù)驅(qū)動隨機子空間法的計算精度與模態(tài)識別能力，是造成差異的另一個原因。

3.5 兩種子空間法對比分析

針對基于QR的數(shù)據(jù)驅(qū)動隨機子空間法與協(xié)方差驅(qū)動隨機子空間法的識別結(jié)果的差異，將采用蒙特卡洛模擬對比兩種方法的優(yōu)劣。

用Matlab程序中randn函數(shù)生成白噪聲激勵，進行100次蒙特卡洛激勵試驗，并統(tǒng)計模態(tài)階數(shù)為10的相關(guān)模態(tài)參數(shù)估計結(jié)果(模態(tài)頻率、模態(tài)阻尼)。如圖6、圖7所示。圖中行代表4個模態(tài);列代表協(xié)方差驅(qū)動隨機子空間法(SSI/cov)與基于QR的數(shù)據(jù)驅(qū)動隨機子空間法(SSI/data)。橫軸為激勵次數(shù)，縱軸為用真實值標準化的估計值，且所有估計值均用點表示，紅線代表每個模態(tài)的真實值，藍線代表每個模態(tài)100次估計結(jié)果的平均值。

圖6 100次仿真的模態(tài)頻率估計結(jié)果Fig.6 Modal eigenfrequency estimation results from 100 simulations

由圖6看出，基于QR分解的數(shù)據(jù)驅(qū)動隨機子空間法(SSI/data)識別出4個模態(tài)(4行)，而協(xié)方差驅(qū)動隨機子空間法(SSI/cov)僅識別出3個模態(tài)(3行)，說明基于QR分解的數(shù)據(jù)驅(qū)動隨機子空間法能更敏感地識別弱勢模態(tài)，估計結(jié)果更集中，計算精度更高。圖7與圖6類似。但模態(tài)阻尼估計結(jié)果計算精度稍差，相對較分散。

為進一步說明以上分析結(jié)果，分別統(tǒng)計模態(tài)頻率和模態(tài)阻尼估計結(jié)果均值及變異系數(shù)見表1、表2。由二表可見，數(shù)據(jù)驅(qū)動隨機子空間法的估計結(jié)果更接近真實值，同一模態(tài)變異系數(shù)更小。進一步表明基于QR分解的數(shù)據(jù)驅(qū)動隨機子空間法的優(yōu)勢。

圖7 100次仿真的模態(tài)阻尼估計結(jié)果Fig.7 Damping ratios estimation results from 100 simulations

表1 頻率估計值均值與變異系數(shù)(定階為10)Tab.1 Mean and variation coefficient of estimated frequencies

表2 阻尼估計值均值與變異系數(shù)(定階為10)Tab.2 Mean and variation coefficient of estimated damping ratios

由蒙特卡洛試驗結(jié)果知，無論在識別模態(tài)的數(shù)目上，或估計的精度上(尤其阻尼)，基于QR分解的數(shù)據(jù)驅(qū)動隨機子空間法明顯優(yōu)于協(xié)方差驅(qū)動隨機子空間法。

4 結(jié)論

本文分析并驗證了造成數(shù)據(jù)驅(qū)動與協(xié)方差驅(qū)動隨機子空間法差異的原因。

通過對數(shù)據(jù)有限性與QR分解算法的使用、研究發(fā)現(xiàn)，基于QR分解的數(shù)據(jù)驅(qū)動隨機子空間法無論在識別參數(shù)的精度上(尤其阻尼)，或在對較弱勢模態(tài)的識別能力上均明顯優(yōu)于協(xié)方差驅(qū)動隨機子空間法。因此，本文推薦用基于QR分解的數(shù)據(jù)驅(qū)動隨機子空間法作為海洋平臺等大型結(jié)構(gòu)的模態(tài)參數(shù)識別工具。

［1］ Wang S Q，Zhang Y T，F(xiàn)eng Y X.Comparative study of outputbased modal identification methods using measured signal from an offshore platform［C］.Proceedings of the ASME 29th International Conference on Ocean，Offshore and Arctic Engineering，2010，2:561-567.

［2］ Hu S，Li P，Vincent H，et al.Modal parameter estimation for jacket-type platforms using free-vibration data［J］.Journal of Waterway，Port，Coastal，and Ocean Engineering，2011，137(5):234-275.

［3］ Cunha A，Caetano E.Experimental modal analysis of civil engineering structures［M］.lvaro Cunha and Elsa Caetano:University of Porto，2006:12-20.

［4］ Overschee P，DeMoor B.Subspace identification for linear system:theory-implementation-applications［M］. Kluwer Academic Publishers Boston/London/Dordrecht，1996.

［5］ Peeters B，DeRoeck G.Stochastic system identification for operational modal analysis:a review［J］.Journal of Dynamic Systems，Measurement and Control，2001，123(4):659-667.

［6］Peeters B，DeRoeck G.Reference-based stochastic subspace identification for output-only modal analysis［J］.Mechanical Systems and Signal Processing，1999，13(6):855-878.

［7］ Hu S，Bao X，Li H.Model order determination and noise removal for modalparameterestimation［J］.Mechanical Systems and Signal Processing，2010，24(6):1605-1620.

［8］ Matlab.Control system toolbox user's guide［M］.Natic，Ma，2004.