協(xié)方差驅(qū)動(dòng)隨機(jī)子空間的Toeplitz矩陣行數(shù)選擇方法

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

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

郵箱：xlhan@seu.edu.cn

協(xié)方差驅(qū)動(dòng)隨機(jī)子空間的Toeplitz矩陣行數(shù)選擇方法

王燕1,2，杭曉晨1,2，姜東1,2，韓曉林1,2，費(fèi)慶國(guó)1,2

(1.江蘇省工程力學(xué)分析重點(diǎn)實(shí)驗(yàn)室,南京210096;2.東南大學(xué)土木工程學(xué)院,南京 210096)

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

關(guān)鍵詞：隨機(jī)子空間方法；阻尼識(shí)別；Toeplitz矩陣；條件數(shù)

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

中圖分類號(hào):TU317

文獻(xiàn)標(biāo)志碼：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

隨機(jī)子空間方法可以在激勵(lì)未知、僅有響應(yīng)數(shù)據(jù)的情況下識(shí)別系統(tǒng)的模態(tài)參數(shù)。根據(jù)識(shí)別過程中具體算法的不同，隨機(jī)子空間可以分為數(shù)據(jù)驅(qū)動(dòng)隨機(jī)子空間和協(xié)方差驅(qū)動(dòng)隨機(jī)子空間[1-2]。協(xié)方差驅(qū)動(dòng)隨機(jī)子空間通過分解結(jié)構(gòu)振動(dòng)響應(yīng)輸出數(shù)據(jù)協(xié)方差矩陣(Toeplitz矩陣)的奇異值得到系統(tǒng)矩陣，進(jìn)而識(shí)別系統(tǒng)的模態(tài)參數(shù)[3]。

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

本文推導(dǎo)了Toeplitz矩陣行數(shù)i與協(xié)方差驅(qū)動(dòng)隨機(jī)子空間方法奇異值分解的理論關(guān)系，并且從條件數(shù)著手，分析了i對(duì)求解系統(tǒng)矩陣的影響。證明了Toeplitz矩陣行數(shù)i不僅對(duì)計(jì)算速度產(chǎn)生影響，對(duì)協(xié)方差驅(qū)動(dòng)隨機(jī)子空間方法的識(shí)別精度也有著不可忽視的影響，繼而提出了協(xié)方差驅(qū)動(dòng)隨機(jī)子空間方法的Toeplitz矩陣行數(shù)i選擇方法。

1協(xié)方差驅(qū)動(dòng)隨機(jī)子空間的理論基礎(chǔ)

在協(xié)方差驅(qū)動(dòng)隨機(jī)子空間方法中，首先將響應(yīng)數(shù)據(jù)表達(dá)為Hankel矩陣，然后計(jì)算協(xié)方差序列組成Toeplitz矩陣，主要作用在保持信號(hào)原有信息的情況下縮減數(shù)據(jù)量，從而減少計(jì)算時(shí)間。再通過矩陣的奇異值分解和特征值分解求出系統(tǒng)矩陣，進(jìn)而識(shí)別結(jié)構(gòu)的模態(tài)參數(shù)。

Hankel矩陣定義為：

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

(2)

(3)

計(jì)算所得協(xié)方差序列組成Toeplitz矩陣，表示如下：

(4)

由狀態(tài)空間方程的定義可以推導(dǎo)：

T=ΓΘ

(5)

如果系統(tǒng)是可觀和可控的，則Toeplitz矩陣的秩為系統(tǒng)階次n。對(duì)Toeplitz矩陣進(jìn)行奇異值分解，秩反映在不為零的奇異值數(shù)量上：

T=USVT=

(6)

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

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

(7)

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

(8)

式中，F(xiàn)為n階非奇異矩陣，為簡(jiǎn)單起見，可取F為單位矩陣I，上式可以變?yōu)椋?/p>

(9)

同樣定義T2如下：

T2=ΓAΘ

(10)

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

(11)

根據(jù)Γ和Θ的定義可知，C等于Γ的前l(fā)行，其中l(wèi)為通道數(shù)。

對(duì)于離散系統(tǒng)，對(duì)狀態(tài)矩陣進(jìn)行特征值分解：

A=ΨΛΨ-1

(12)

最后，根據(jù)系統(tǒng)矩陣的特征值及特征向量識(shí)別出系統(tǒng)的模態(tài)參數(shù)頻率和阻尼。

2矩陣行數(shù)與奇異值分解的關(guān)系

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

(13)

式中S表示真實(shí)信號(hào)，C表示噪聲。

根據(jù)系統(tǒng)的階次n設(shè)定閥值，n的確定可以參見文獻(xiàn)[10-12]。對(duì)所計(jì)算的Toeplitz矩陣進(jìn)行奇異值分解得到：

(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)為對(duì)角矩陣。n≤i為系統(tǒng)的階數(shù)。

構(gòu)造矩陣：

(15)

由于信號(hào)與噪聲不相關(guān)，所以上式可以寫成：

(16)

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

HC=η2Ii

(17)

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

(18)

(19)

(20)

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

(21)

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

(22)

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

(23)

即:

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

(24)

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

3矩陣行數(shù)對(duì)Toeplitz矩陣條件數(shù)的影響

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

T2=TA

(25)

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

由式(4)可知，Toeplitz矩陣為i×i的方陣(通道數(shù)為l時(shí))，i為人為選擇的參數(shù)。當(dāng)i取不同值時(shí)，Toeplitz矩陣的條件數(shù)也會(huì)發(fā)生變化，從而導(dǎo)致輸出數(shù)據(jù)的擾動(dòng)對(duì)系統(tǒng)矩陣的求解產(chǎn)生不同的影響。當(dāng)Toeplitz矩陣條件數(shù)越小，響應(yīng)數(shù)據(jù)的擾動(dòng)對(duì)系統(tǒng)矩陣的擾動(dòng)越小，從而導(dǎo)致根據(jù)系統(tǒng)矩陣求解的模態(tài)參數(shù)誤差更小。所以，Toeplitz矩陣條件數(shù)越小，系統(tǒng)的識(shí)別精度越高。

4Toeplitz矩陣行數(shù)i的選擇方法

設(shè)響應(yīng)數(shù)據(jù)總長(zhǎng)度為s，Hankel矩陣行數(shù)和列數(shù)滿足關(guān)系：

s=2i+j-1

(26)

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

理論上j→∞，根據(jù)經(jīng)驗(yàn)值j>20i，所以i需滿足

(27)

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

i>n

(28)

式中n為系統(tǒng)的階數(shù)。

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

5仿真算例

5.1兩自由度彈簧振子系統(tǒng)

圖1　彈簧振子系統(tǒng) Fig.1 Spring mass system

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

圖2為該系統(tǒng)的穩(wěn)定圖，根據(jù)穩(wěn)定圖確定系統(tǒng)階數(shù)為4。由于頻率的識(shí)別結(jié)果精度比較高，這里不做比較。表1為Toeplitz矩陣行數(shù)i按公差為2的等差數(shù)列選取的阻尼識(shí)別結(jié)果以及相應(yīng)的Toeplitz矩陣譜條件數(shù)。

圖2　彈簧振子系統(tǒng)的穩(wěn)定圖 Fig.2 Stabilization diagram of spring mass system

觀察表1可以發(fā)現(xiàn)：

(1)協(xié)方差驅(qū)動(dòng)隨機(jī)子空間算法能比較精確的識(shí)別出該彈簧振子系統(tǒng)的阻尼參數(shù)；

(2)當(dāng)i=12時(shí)，第一階模態(tài)的識(shí)別誤差最小，達(dá)到-0.19%。當(dāng)i=6時(shí)，第二階模態(tài)的識(shí)別誤差最小。

(3)當(dāng)i=10時(shí)，Toeplitz矩陣的譜條件數(shù)最小，為7.22e15，此時(shí)兩階識(shí)別結(jié)果的平均絕對(duì)誤差最小，為1.07%，說明在條件數(shù)最小時(shí)獲得最優(yōu)的識(shí)別結(jié)果。

表1　彈簧振子系統(tǒng)的阻尼識(shí)別結(jié)果

5.2切尖三角翼的阻尼識(shí)別

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

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

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

同樣地，由于頻率的識(shí)別結(jié)果精度比較高，這里只比較阻尼的識(shí)別結(jié)果。采用仿真一的方案選取i進(jìn)行識(shí)別計(jì)算。根據(jù)如圖四所示的穩(wěn)定圖，確定系統(tǒng)的階數(shù)為10階。阻尼識(shí)別結(jié)果和相應(yīng)的譜條件數(shù)如表2所示。

由表2可以觀察到：

(1)協(xié)方差驅(qū)動(dòng)隨機(jī)子空間算法能夠比較精確地識(shí)別出系統(tǒng)的前四階模態(tài)參數(shù)。但當(dāng)i取值較大時(shí)，識(shí)別結(jié)果出現(xiàn)遺漏模態(tài)。

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

表2　切尖三角翼的阻尼識(shí)別結(jié)果及Toeplitz矩陣的譜條件數(shù)

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

6結(jié)論

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

參考文獻(xiàn)

[1]辛峻峰, 王樹青, 劉福順. 數(shù)據(jù)驅(qū)動(dòng)與協(xié)方差驅(qū)動(dòng)隨機(jī)子空間法差異化分析[J]. 振動(dòng)與沖擊, 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]李永軍, 馬立元, 王天輝, 等. 協(xié)方差驅(qū)動(dòng)子空間模態(tài)參數(shù)辨識(shí)方法改進(jìn)分析[J]. 中國(guó)機(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ī)子空間方法的模態(tài)參數(shù)識(shí)別研究[J].中國(guó)航空學(xué)報(bào)，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]劉東霞. 基于隨機(jī)子空間法的梁橋模態(tài)參數(shù)識(shí)別 [D]. 成都: 西南交通大學(xué), 2008.

[9]辛峻峰, 盛進(jìn)路, 張永波. 數(shù)據(jù)驅(qū)動(dòng)隨機(jī)子空間法矩陣維數(shù)選擇與噪聲問題研究[J]. 振動(dòng)與沖擊, 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]樊江玲. 基于輸出響應(yīng)的模態(tài)參數(shù)辨識(shí)方法研究 [D]. 上海：上海交通大學(xué), 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ī)顫振試驗(yàn)數(shù)據(jù)預(yù)處理中的應(yīng)用[J]. 振動(dòng)與沖擊, 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]易大義, 陳道琦. 數(shù)值分析引論[M]. 杭州：浙江大學(xué)出版社, 1998, 9.