協方差驅動隨機子空間的Toeplitz矩陣行數選擇方法

第一作者王燕女，碩士生，1990年生

通信作者韓曉林男，教授，1958年生

郵箱：xlhan@seu.edu.cn

協方差驅動隨機子空間的Toeplitz矩陣行數選擇方法

王燕1,2，杭曉晨1,2，姜東1,2，韓曉林1,2，費慶國1,2

(1.江蘇省工程力學分析重點實驗室,南京210096;2.東南大學土木工程學院,南京 210096)

摘要：隨機子空間識別算法是一種基于環境激勵的模態參數識別方法，僅需要響應時程便可識別模態參數。其中，協方差驅動隨機子空間方法中Toeplitz矩陣行數的選取直接影響識別精度。通過構造相關矩陣，研究了Toeplitz矩陣行數i對協方差驅動隨機子空間方法中奇異值分解去噪能力的影響。引入Toeplitz矩陣條件數，根據i與Toeplitz矩陣條件數的關系再次證明了i對識別精度的影響。研究了Toeplitz矩陣行數i的選擇方法。采用兩自由度彈簧振子系統和切尖三角翼模型兩個仿真算例研究了Toeplitz矩陣行數i的選擇方法。結果表明：在確定合適的系統階數的前提下，Toeplitz矩陣的條件數越小識別精度越高。

關鍵詞：隨機子空間方法；阻尼識別；Toeplitz矩陣；條件數

收稿日期：2014-01-28修改稿收到日期：2014-04-03

中圖分類號:TU317

文獻標志碼：A

DOI:10.13465/j.cnki.jvs.2015.07.011

Abstract:Stochastic Subspace Identification is a parameter identification method, which can effectively obtain modal parameters from the structural signal under ambient excitation. The choice of Toeplitz matrix row number directly influences the accuracy of identification. By constructing a correlation matrix, the influence of the dimension of Toeplitz matrix i on the denoising ability via SVD was derived. The concept of condition number was introduced in solving the system matrix. According to the relationship between i and condition number of Toeplitz matrix, it is proved once again that i has influence on identification accuracy. Then the selection method of Toeplitz matrix row number i was studied. Two examplic simulations in regard to a two-degree spring mass vibration system and a cropped delta wing model were presented to show the method in the selection of i. The results show that on the premise of determining a suitable system order, the smaller the Toeplitz matrix condition number is, the higher the identification accuracy is.

Selection method of Toeplitz matrix row number based on covariance driven stochastic subspace identification

WANGYan1, 2,HANGXiao-chen1, 2,JIANGDong1, 2,HANXiao-lin1, 2,FEIQing-guo1, 2(1. Jiangsu Key Laboratory of Engineering Mechanics, Nanjing 210096,China;2. Department of Engineering Mechanics, Southeast University, Nanjing 210096, China)

Key words:stochastic subspace identification; damping identification; Toeplitz matrix; condition number

隨機子空間方法可以在激勵未知、僅有響應數據的情況下識別系統的模態參數。根據識別過程中具體算法的不同，隨機子空間可以分為數據驅動隨機子空間和協方差驅動隨機子空間[1-2]。協方差驅動隨機子空間通過分解結構振動響應輸出數據協方差矩陣(Toeplitz矩陣)的奇異值得到系統矩陣，進而識別系統的模態參數[3]。

協方差驅動隨機子空間方法中Toeplitz矩陣行數i與數據驅動隨機子空間方法中的Hankel矩陣行數i概念相同，都是人為選擇的參數。一般認為，i的選取只會對計算速度產生影響。Overschee[4-6]最先提出了數據驅動隨機子空間方法，并對Hankel矩陣行數的取值范圍做了簡要介紹，他提出Hankel矩陣的行數需大于系統階數。陳愛林等[7]通過統計分析得到了識別精度最佳時Hankel矩陣行數與系統階數的關系，總結出實際應用中隨機子空間方法關于i的選取方法。劉東霞[8]研究了協方差驅動隨機子空間算法在橋梁中的模態參數識別，分析了Hankel矩陣行和列的關系。辛峻峰等[9]推導了Hankel矩陣維數與數據驅動隨機子空間方法去噪能力的理論關系，研究了數據驅動隨機子空間方法Hankel矩陣維數對識別精度帶來的影響。

本文推導了Toeplitz矩陣行數i與協方差驅動隨機子空間方法奇異值分解的理論關系，并且從條件數著手，分析了i對求解系統矩陣的影響。證明了Toeplitz矩陣行數i不僅對計算速度產生影響，對協方差驅動隨機子空間方法的識別精度也有著不可忽視的影響，繼而提出了協方差驅動隨機子空間方法的Toeplitz矩陣行數i選擇方法。

1協方差驅動隨機子空間的理論基礎

在協方差驅動隨機子空間方法中，首先將響應數據表達為Hankel矩陣，然后計算協方差序列組成Toeplitz矩陣，主要作用在保持信號原有信息的情況下縮減數據量，從而減少計算時間。再通過矩陣的奇異值分解和特征值分解求出系統矩陣，進而識別結構的模態參數。

Hankel矩陣定義為：

Hankel矩陣Y0/2i-1是包括“2i行塊”、j列塊的矩陣。每一塊行由l行組成，l是輸出通道數。理論上假設j→∞，但實際上j不可能是無窮大的。Hankel矩陣分塊為“過去”和“將來”兩部分，下標0/2i-1表示Hankel矩陣的第一列的第一塊行和最后一塊行的下標，下標p和f分別表示“過去”和“將來”。矩陣分成Yp和Yf將Y0/2i-1分成相等的兩部分，每一部分都有i塊行。

(2)

(3)

計算所得協方差序列組成Toeplitz矩陣，表示如下：

(4)

由狀態空間方程的定義可以推導：

T=ΓΘ

(5)

如果系統是可觀和可控的，則Toeplitz矩陣的秩為系統階次n。對Toeplitz矩陣進行奇異值分解，秩反映在不為零的奇異值數量上：

T=USVT=

(6)

式中，U，V都是正交矩陣，S是由正奇異矩陣組成的對角陣，U1∈Rli×n,S1∈Rn×n,V1∈Rli×n，奇異值按降序排列：

σ1≥σ2≥…≥σn≥0

(7)

比較式(4)和式(5)，可觀矩陣和可控矩陣可以表示如下：

(8)

式中，F為n階非奇異矩陣，為簡單起見，可取F為單位矩陣I，上式可以變為：

(9)

同樣定義T2如下：

T2=ΓAΘ

(10)

T2和T具有相同的結構，T2包含的協方差從2到i+1。根據式(9)和式(10)可以解得A:

(11)

根據Γ和Θ的定義可知，C等于Γ的前l行，其中l為通道數。

對于離散系統，對狀態矩陣進行特征值分解：

A=ΨΛΨ-1

(12)

最后，根據系統矩陣的特征值及特征向量識別出系統的模態參數頻率和阻尼。

2矩陣行數與奇異值分解的關系

由式(3)和(4)可知Toeplitz矩陣由以下兩部分組成：

(13)

式中S表示真實信號，C表示噪聲。

根據系統的階次n設定閥值，n的確定可以參見文獻[10-12]。對所計算的Toeplitz矩陣進行奇異值分解得到：

(14)

式中，正交矩陣U∈Ri×i、U1∈Ri×n、U2∈Ri×(i-n)，且V∈Ri×i、V1∈Ri×n、V2∈Ri×(i-n)也為正交矩陣。S∈Ri×i、S1∈Rn×n、S2∈R(i-n)×(i-n)為對角矩陣。n≤i為系統的階數。

構造矩陣：

(15)

由于信號與噪聲不相關，所以上式可以寫成：

(16)

由上式可知，HC為噪聲C的協方差矩陣，設已知噪聲的方差為η2，則噪聲C的協方差矩陣為[13]

HC=η2Ii

(17)

式中，Ii為i階單位陣。根據式(6)有

(18)

(19)

(20)

由式(17)和式(20)可以得到：

(21)

則結合式(18)~(21)知：

(22)

為了便于分析，將S2表示為如下形式：

(23)

即:

S1=(S2n-iη2In)1/2

(24)

式(24)表明，在數據總量不變時，S1取決于Sn，Toeplitz矩陣行數i和系統的階數。而Sn取決于系統的階數n。所以，在數據總量固定時Toeplitz矩陣行數i和噪聲方差η都會影響SVD的去噪能力，從而影響算法的識別精度。

3矩陣行數對Toeplitz矩陣條件數的影響

由式(10)對T2的定義可知：

T2=TA

(25)

為求系統矩陣A則需求解式(25)所示的線性方程組，根據系統矩陣A可以解得系統的模態參數。因此，系統矩陣A的精度決定了模態參數識別結果的識別精度。式(25)可能是病態的，此時T和T2的微小擾動將引起系統矩陣A較大的改變。由此可以根據Toeplitz矩陣即T矩陣的條件數來判斷該線性方程組是否是病態方程組。

由式(4)可知，Toeplitz矩陣為i×i的方陣(通道數為l時)，i為人為選擇的參數。當i取不同值時，Toeplitz矩陣的條件數也會發生變化，從而導致輸出數據的擾動對系統矩陣的求解產生不同的影響。當Toeplitz矩陣條件數越小，響應數據的擾動對系統矩陣的擾動越小，從而導致根據系統矩陣求解的模態參數誤差更小。所以，Toeplitz矩陣條件數越小，系統的識別精度越高。

4Toeplitz矩陣行數i的選擇方法

設響應數據總長度為s，Hankel矩陣行數和列數滿足關系：

s=2i+j-1

(26)

式中i表示Hankel矩陣有2i個行塊，同時i與Toeplitz矩陣的行塊或列塊相等。j表示Hankel矩陣有j個列塊。

理論上j→∞，根據經驗值j>20i，所以i需滿足

(27)

另一方面，Toeplitz矩陣的行數i需滿足[4]：

i>n

(28)

式中n為系統的階數。

綜上所述，提出一種通用的Toeplitz矩陣行數i取值方法：在滿足(27)和(28)的前提條件下，確定合適的系統階數，計算不同i值對應的Toeplitz矩陣的條件數，選取條件數較小時對應的i值，進行識別計算。

5仿真算例

5.1兩自由度彈簧振子系統

圖1　彈簧振子系統 Fig.1 Spring mass system

如圖1所示，k1=10 000N/m，m1=0.5 kg，k2=5 000 N/m，m2=0.2 kg。設模態阻尼比為3%，該系統的第一階和第二階頻率分別為17.062 Hz和33.135 Hz。在質量塊m2施加白噪聲激勵，響應數據采樣頻率為160 Hz (截止頻率為50 Hz)，采樣時間為100 s，數據總長為16 000。

圖2為該系統的穩定圖，根據穩定圖確定系統階數為4。由于頻率的識別結果精度比較高，這里不做比較。表1為Toeplitz矩陣行數i按公差為2的等差數列選取的阻尼識別結果以及相應的Toeplitz矩陣譜條件數。

圖2　彈簧振子系統的穩定圖 Fig.2 Stabilization diagram of spring mass system

觀察表1可以發現：

(1)協方差驅動隨機子空間算法能比較精確的識別出該彈簧振子系統的阻尼參數；

(2)當i=12時，第一階模態的識別誤差最小，達到-0.19%。當i=6時，第二階模態的識別誤差最小。

(3)當i=10時，Toeplitz矩陣的譜條件數最小，為7.22e15，此時兩階識別結果的平均絕對誤差最小，為1.07%，說明在條件數最小時獲得最優的識別結果。

表1　彈簧振子系統的阻尼識別結果

5.2切尖三角翼的阻尼識別

圖三為切尖三角翼的幾何尺寸，翼板厚度為0.016 m。彈性模量為71.161 2 GPa，剪切模量為27.100 4 GPa，泊松比為0.324 8。選取翼板翼根為固定端，建立有限元模型。在翼板的上表面施加白噪聲面壓。響應數據采樣頻率為1 800 Hz，采樣時間為10 s，數據總長度為18 000。設阻尼為3%的模態阻尼，對應各階頻率均為3%，提取如圖4所示的六個點的位移響應。

圖3　切尖三角翼的幾何尺寸 Fig.3 The geometric size of the cropped delta wing

圖4　切尖三角翼的有限元模型 Fig.4 The finite element model of the cropped delta wing

同樣地，由于頻率的識別結果精度比較高，這里只比較阻尼的識別結果。采用仿真一的方案選取i進行識別計算。根據如圖四所示的穩定圖，確定系統的階數為10階。阻尼識別結果和相應的譜條件數如表2所示。

由表2可以觀察到：

(1)協方差驅動隨機子空間算法能夠比較精確地識別出系統的前四階模態參數。但當i取值較大時，識別結果出現遺漏模態。

(2)當i=12~28時，各階模態的阻尼識別的最小誤差不盡相同。當i=12時，Toeplitz矩陣的譜條件數最小，為1.73e17，此時前四階模態阻尼的平均絕對誤差最小，為1.23%。可以得到與兩自由度算例相同的結論，即：條件數最小時獲得最優的識別結果。

表2　切尖三角翼的阻尼識別結果及Toeplitz矩陣的譜條件數

圖5　切尖三角翼模型的穩定圖 Fig.5 Stabilization diagram of the cropped delta wing

6結論

本文利用奇異值分解的理論推導和線性方程組中條件數這一數學工具，分析了協方差驅動隨機子空間中Toeplitz矩陣行數i的選擇對識別精度的影響。提出了一種協方差驅動隨機子空間算法Toeplitz矩陣行數i的選擇方法：在確定合適的系統階數后，選取不同的i值，分別計算不同i值對應的Toeplitz矩陣的條件數，選取條件數較小時對應的i值，進行識別計算。通過兩個算例仿真表明：i越大，求解系統矩陣A的線性方程組病態越嚴重；Toeplitz矩陣行數i對識別精度的影響至關重要，不同系統為達到最好識別精度，i的取值不是固定的；當條件數較小時，識別精度比較理想，由此證明了本文提出的Toeplitz矩陣行數i選擇方法的可行性。根據本文的結論，更加方便地得到i的選取方法，并且該方法確定的i能更加精確的識別出系統的模態參數。

參考文獻

[1]辛峻峰, 王樹青, 劉福順. 數據驅動與協方差驅動隨機子空間法差異化分析[J]. 振動與沖擊, 2013,32(9):1-5.

XIN Jun-feng, WANG Shu-qing, LIU Fu-shun. Performance comparison for data-driven and covariance-driven stochastic subspace identification method[J]. Journal of Vibration and Shock, 2013,32(9):1-5.

[2]Boonyapinyo V, Janesupasaeree T. Data-driven stochastic subspace identification of flutter derivatives of bridge decks[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2010, 98(12): 784-799.

[3]李永軍, 馬立元, 王天輝, 等. 協方差驅動子空間模態參數辨識方法改進分析[J]. 中國機械工程, 2012, 23(13): 1533-1536.

LI Yong-jun, MA Li-yuan, WANG Tian-hui, et al. An improvement on subspace modal parameter identification algorithm driven by covariance[J]. China Mechanical Engineering, 2012, 23(13):1533-1536.

[4]Van Overschee P, De Moor B. Subspace identification for linear systems: theory, implementation, applications. dordrecht, Netherlands: Kluwer Academic Publishers, 1996.

[5]Van Overschee P, De Moor B. Subspace algorithms for the stochastic identification problem[J]. Proceedings of the 30thIEEE Conference on Decision and Control, pages 1321-1326,Brighton,UK,1991.

[6]Van Overschee P, De Moor B. Subspace algorithms for the stochastic identification problem[J]. Automatica, 1993, 29(3): 649-660.

[7]陳愛林,張令彌.基于隨機子空間方法的模態參數識別研究[J].中國航空學報，2001(4):222-228.

CHEN Ai-lin, ZHANG Ling-mi. Study of stochastic subspace system identification method[J].Chinese Journal of Aeronautics,2001(4):222-228.

[8]劉東霞. 基于隨機子空間法的梁橋模態參數識別 [D]. 成都: 西南交通大學, 2008.

[9]辛峻峰, 盛進路, 張永波. 數據驅動隨機子空間法矩陣維數選擇與噪聲問題研究[J]. 振動與沖擊, 2013, 32(16): 152-157.

XIN Jun-feng, SHENG Jin-lu, ZHANG Yong-bo. Relation between noise and matrix dimension of data-driven stochastic subspace identification method[J]. Journal of Vibration and Shock, 2013, 32(16):152-157.

[10]樊江玲. 基于輸出響應的模態參數辨識方法研究 [D]. 上海：上海交通大學, 2007.

[11]Hong A L, Ubertini F, Betti R. New stochastic subspace approach for system identification and its application to long-span bridges[J]. Journal of Engineering Mechanics, 2012, 139(6): 724-736.

[12]Loh C H, Liu Y C, Ni Y Q. SSA-based stochastic subspace identification of structures from output-only vibration measurements[J]. Smart Structures and Systems, 2012, 10(4-5): 331-351.

[13]張波, 李健君. 基于 Hankel 矩陣與奇異值分解 (SVD) 的濾波方法以及在飛機顫振試驗數據預處理中的應用[J]. 振動與沖擊, 2009, 28(2): 162-166.

ZHANG Bo, LI Jian-jun. Denoising method based on Hankel matrix and SVD and its application in flight flutter testing data preprocessing[J]. Journal of Vibration and Shock, 2009, 28(2):162-166.

[14]易大義, 陳道琦. 數值分析引論[M]. 杭州：浙江大學出版社, 1998, 9.