辛峻峰，王樹青，劉福順

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

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

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

1 隨機(jī)子空間法

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

式中:xk∈Rn為離散時(shí)間狀態(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。

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

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

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

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

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

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

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

2.2 QR分解

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

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

3 數(shù)值算例

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

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

每個(gè)單元質(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 白噪聲激勵(lì)實(shí)驗(yàn)

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

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

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

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

為驗(yàn)證兩種方法差異的存在，分別用基于QR的數(shù)據(jù)驅(qū)動隨機(jī)子空間法與協(xié)方差驅(qū)動隨機(jī)子空間法分析測量數(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%，同時(shí)響應(yīng)阻尼識別結(jié)果誤差為5%，即:

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

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

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

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

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

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

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

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

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

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

3.4 QR分解算法影響驗(yàn)證

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

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

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

3.5 兩種子空間法對比分析

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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