數據驅動與協方差驅動隨機子空間法差異化分析

辛峻峰，王樹青，劉福順

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

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

DeMoor等［4］證明了數據驅動隨機子空間算法與協方差驅動隨機子空間算法在理論上的完全一致。Peters等［5］認為在應用過程中，兩種方法的模態識別結果也一致。但在對導管架平臺模態分析過程中，兩種方法的識別結果并不一致。對此，本文研究探討兩種方法產生差異的原因，繼而進行對比分析。

1 隨機子空間法

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

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

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

兩種隨機子空間算法的本質區別在于“過去”與“將來”相關性的計算方式不同:數據驅動隨機子空間法用投影計算:

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

由此數據驅動與協方差驅動隨機子空間法的關系主要為投影與協方差矩陣之間的關系。

2 數據驅動隨機子空間法與協方差驅動隨機子空間法之關系

推導［4］，投影Pi與協方差矩陣Ci為相似矩陣。據相似矩陣性質知，二者奇異值相同，從而可得相同模態參數。與文獻［5］結果相同。但本文研究發現，被忽略的2個細節能使兩種方法在應用中產生差異。

2.1 數據有限性

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

2.2 QR分解

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

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

3 數值算例

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

圖1 自由度彈簧-質量-阻尼系統Fig.1 A 5-DOF mass-spring-dashpot system

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

3.1 白噪聲激勵實驗

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

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

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

3.2 差異現象驗證

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

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

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

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

本文穩定圖用信號的傅里葉變換為背景，顯示模態階數在0～30之間模態識別結果情況。對同一模態，若前后兩次識別結果同一模態頻率誤差為1%，同時響應阻尼識別結果誤差為5%，即:

則此次識別結果標識為穩定，否則為不穩定［8］。

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

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

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

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

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

3.3 數據有限性影響驗證

本節將驗證數據有限性影響數據驅動隨機子空間法的識別結果。

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

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

3.4 QR分解算法影響驗證

本節將驗證QR分解算法影響了數據驅動隨機子空間法的識別結果。

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

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

3.5 兩種子空間法對比分析

針對基于QR的數據驅動隨機子空間法與協方差驅動隨機子空間法的識別結果的差異，將采用蒙特卡洛模擬對比兩種方法的優劣。

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

圖6 100次仿真的模態頻率估計結果Fig.6 Modal eigenfrequency estimation results from 100 simulations

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

為進一步說明以上分析結果，分別統計模態頻率和模態阻尼估計結果均值及變異系數見表1、表2。由二表可見，數據驅動隨機子空間法的估計結果更接近真實值，同一模態變異系數更小。進一步表明基于QR分解的數據驅動隨機子空間法的優勢。

圖7 100次仿真的模態阻尼估計結果Fig.7 Damping ratios estimation results from 100 simulations

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

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

由蒙特卡洛試驗結果知，無論在識別模態的數目上，或估計的精度上(尤其阻尼)，基于QR分解的數據驅動隨機子空間法明顯優于協方差驅動隨機子空間法。

4 結論

本文分析并驗證了造成數據驅動與協方差驅動隨機子空間法差異的原因。

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

［1］ Wang S Q，Zhang Y T，Feng 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.