改進(jìn)的TW-API方法及其在時(shí)變模態(tài)參數(shù)辨識(shí)中的應(yīng)用

倪智 宇，譚述君，吳志剛，

（1.大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧大連116024；2.大連理工大學(xué)航空航天學(xué)院，遼寧 大連116024）

改進(jìn)的TW-API方法及其在時(shí)變模態(tài)參數(shù)辨識(shí)中的應(yīng)用

倪智 宇1，譚述君2，吳志剛1，2

（1.大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧大連116024；2.大連理工大學(xué)航空航天學(xué)院，遼寧 大連116024）

改進(jìn)了一種用于模態(tài)參數(shù)辨識(shí)的截?cái)啻氨平鼉绲涌臻g方法（Truncated Window Approximated Power Iteration，TW-API）。該算法采用輸入輸出信號(hào)的截?cái)啻盎ハ嚓P(guān)函數(shù)建立信號(hào)主子空間的遞推關(guān)系，應(yīng)用子空間投影來(lái)進(jìn)行辨識(shí)。跟傳統(tǒng)的TW-API方法相比，改進(jìn)后的TW-API方法基于信號(hào)子空間理論，簡(jiǎn)化了數(shù)據(jù)矩陣遞推的過(guò)程，顯著地減少了算法的整體計(jì)算量和計(jì)算時(shí)間，同時(shí)仍能較好地保持辨識(shí)精度。在數(shù)值算例中，選用二連桿機(jī)械臂和在軌航天器模型這兩個(gè)算例，分別進(jìn)行了時(shí)變模態(tài)參數(shù)的辨識(shí)，并將結(jié)果與經(jīng)典的投影估計(jì)子空間跟蹤方法（Projection Approximation Subspace Tracking，PAST）和逼近冪迭代（Approximated Power Iteration，API）遞推方法進(jìn)行比較，4種辨識(shí)方法的仿真結(jié)果表明本文提出的新算法能有效辨識(shí)系統(tǒng)的時(shí)變模態(tài)參數(shù)。

模態(tài)參數(shù)辨識(shí)；時(shí)變系統(tǒng)；遞推子空間方法；航天器

引 言

在現(xiàn)代工程應(yīng)用中，航空航天、機(jī)械制造、交通等領(lǐng)域的許多問(wèn)題都需要考慮到系統(tǒng)的時(shí)變特性，例如大型空間結(jié)構(gòu)的搬運(yùn)裝配、太陽(yáng)能帆板的展開(kāi)與旋轉(zhuǎn)、機(jī)械臂的運(yùn)動(dòng)以及高速列車(chē)的行駛等。由于線性時(shí)變系統(tǒng)的特性比定常系統(tǒng)更為復(fù)雜，而且定常系統(tǒng)的模態(tài)分析理論并不完全適用于時(shí)變系統(tǒng)，所以時(shí)變系統(tǒng)的理論到現(xiàn)在仍然處于建立的初步階段，還沒(méi)有在一個(gè)完整的框架下描述。從工程應(yīng)用的角度來(lái)說(shuō)，Liu等［1］提出了“偽模態(tài)參數(shù)”的概念，即按照定常系統(tǒng)中模態(tài)參數(shù)的概念，利用“時(shí)間凝固法”將其擴(kuò)展到時(shí)變系統(tǒng)中，即認(rèn)為在極小的時(shí)間間隔內(nèi)的模態(tài)參數(shù)是不變的，從而在各個(gè)采樣時(shí)刻分別計(jì)算出當(dāng)前時(shí)刻的模態(tài)參數(shù)，利用這些偽模態(tài)參數(shù)反映系統(tǒng)的時(shí)變特性。

關(guān)于時(shí)變系統(tǒng)的模態(tài)參數(shù)辨識(shí)，當(dāng)今主要的處理思想之一是利用子空間方法。子空間方法由于具有良好的魯棒性，所以得到了廣泛的應(yīng)用。子空間方法根據(jù)辨識(shí)形式的不同，可以分為整體數(shù)據(jù)辨識(shí)和遞推數(shù)據(jù)辨識(shí)方法。在利用整體數(shù)據(jù)辨識(shí)方面，K Liu［1-2］首先提出了一種利用輸入輸出整體數(shù)據(jù)辨識(shí)系統(tǒng)模型的時(shí)變子空間方法。于開(kāi)平［3］等利用該方法對(duì)移動(dòng)質(zhì)量簡(jiǎn)支梁系統(tǒng)進(jìn)行了辨識(shí)。李會(huì)娜［4-5］等在模態(tài)參數(shù)辨識(shí)的基礎(chǔ)上對(duì)時(shí)變系統(tǒng)物理參數(shù)進(jìn)行了辨識(shí)。在遞推辨識(shí)子空間方面，Tasker［6］和Bosse［7］提出了一種利用URV分解的遞推子空間方法。Yang［8-9］等根據(jù)信號(hào)處理的原理，提出了投影估計(jì)子空間跟蹤方法（Projection Approximation Subspace Tracking，PAST）。龐世偉［10］等在此基礎(chǔ)上提出了固定長(zhǎng)度平移窗改進(jìn)子空間方法，并利用PAST法辨識(shí)了移動(dòng)質(zhì)量-簡(jiǎn)支梁的時(shí)變模態(tài)參數(shù)［11］。Real［12］等提出了快速估計(jì)子空間方法（Fast Approximate Subspace Tracking，F(xiàn)AST）。HUA［13-14］等基于PAST法提出了自然冪迭代（Natural Power Iteration，NPI）方法。Badeau［15-16］等針對(duì)PAST法信號(hào)子空間列向量不能?chē)?yán)格正交的問(wèn)題，提出了逼近冪迭代（Approximated Power Iteration，API）方法，對(duì)信號(hào)子空間投影乘以修正后的正交矩陣，使其投影子空間保持嚴(yán)格正交的同時(shí)可以跟蹤快變的主子空間，并且在此基礎(chǔ)上提出了快速API方法（Fast API，F(xiàn)API）。與整體數(shù)據(jù)辨識(shí)方法相比，遞推辨識(shí)方法不需要逐一時(shí)刻進(jìn)行奇異值分解的計(jì)算，而是通過(guò)信號(hào)子空間的迭代以避免較高維數(shù)的Hankel矩陣求解問(wèn)題，從而減小了計(jì)算量。

本文基于API方法的理論，不同于一般的無(wú)限指數(shù)窗API方法，而是選用截?cái)啻白涌臻g方法（Truncated Window API，TW-API）進(jìn)行遞歸辨識(shí)。文［16］中提出的TW-API方法需要給出辨識(shí)時(shí)刻之前各個(gè)時(shí)刻的子空間輸入信號(hào)作為遞歸總體數(shù)據(jù)矩陣，給定初值以用來(lái)遞歸計(jì)算。與文［15-16］相比，本文簡(jiǎn)化了輸出信號(hào)的數(shù)據(jù)處理過(guò)程，不需要構(gòu)建遞歸總體數(shù)據(jù)矩陣；而且從信號(hào)子空間原理的角度出發(fā)，對(duì)遞推過(guò)程中需要較大計(jì)算量的步驟做了遞推的簡(jiǎn)化，使得方法在保持基本計(jì)算精度的前提下大幅減少了計(jì)算量；最后給出了改進(jìn)后的TWAPI的計(jì)算流程，并討論了截?cái)啻伴L(zhǎng)度對(duì)辨識(shí)的影響。數(shù)值仿真選用二連桿機(jī)械臂和技術(shù)試驗(yàn)衛(wèi)星-VIII（Engineering Test Satellite-VIII，ETS-VIII）模型這兩個(gè)算例，采用本文改進(jìn)后的遞推方法進(jìn)行時(shí)變模態(tài)參數(shù)的辨識(shí)。此外還將該方法與現(xiàn)有的TWAPI，PAST和指數(shù)窗API三種方法進(jìn)行了計(jì)算精度和計(jì)算效率的比較，從而驗(yàn)證了本文算法的有效性。

1 激勵(lì)與響應(yīng)數(shù)據(jù)的前處理

為了利用子空間跟蹤算法進(jìn)行參數(shù)的辨識(shí)，輸入輸出數(shù)據(jù)的前處理是必不可少的，無(wú)論是應(yīng)用FAST，PAST或是API遞推方法，均需要通過(guò)系統(tǒng)輸入輸出數(shù)據(jù)矩陣構(gòu)建對(duì)應(yīng)時(shí)刻的狀態(tài)量變化作為子空間跟蹤中的輸入信號(hào)參數(shù)，然后通過(guò)在子空間追蹤方法中不斷更新該狀態(tài)量（即不斷更新輸入輸出數(shù)據(jù)），來(lái)達(dá)到對(duì)系統(tǒng)的能觀性矩陣（信號(hào)子空間）的追蹤求解，進(jìn)而得到時(shí)變系統(tǒng)的模態(tài)參數(shù)。

1.1 信號(hào)子空間中輸入的迭代更新

在這一部分中，采用文［6-7］中給出的方法進(jìn)行輸入輸出數(shù)據(jù)的預(yù)處理?？紤]如下的離散時(shí)變系統(tǒng)其中方程中系統(tǒng)的階數(shù)為n，輸入的維數(shù)為r，輸出維數(shù)為m。則某一時(shí)刻的輸入矩陣U（k）和輸出矩陣Y（k）以及下一時(shí)刻的輸入數(shù)據(jù)ˉu（k+1）和輸出數(shù)據(jù)ˉy（k+1）的Hankel矩陣形式分別可以寫(xiě)為

其中

其中k為對(duì)應(yīng)的時(shí)刻，M為選取的合適的Hankel矩陣維數(shù)。根據(jù)數(shù)據(jù)更新的遞歸方法［6-7］，令

則通過(guò)輸入輸出數(shù)據(jù)生成的狀態(tài)量z（k）的更新方式為

在式（4）和（6）中，下一個(gè)時(shí)刻k的［U（k）UT（k）］-1和［Y（k）U?（k）］的更新方式為

［U（k）UT（k）］-1=［U（k-1）UT（k-1）］-1-

通過(guò)式（4）～（8）的遞推，利用各個(gè)時(shí)刻的ˉu（k）和ˉy（k）構(gòu)建出每個(gè)時(shí)刻的狀態(tài)量z（k）（即在信號(hào)子空間中的每個(gè)時(shí)刻的輸入信號(hào)），以滿足子空間迭代求解中的數(shù)據(jù)更新條件。

1.2 初始信號(hào)主子空間的構(gòu)建

選取時(shí)變系統(tǒng)的任意k時(shí)刻作為初始時(shí)刻，構(gòu)建輸入矩陣U（k）的投影正交矩陣U⊥（k），使得U（k）U⊥（k）=0，則

式中U?（k）為U（k）的偽逆，I為單位矩陣。

通過(guò)等式右乘矩陣的投影正交矩陣，消去第2項(xiàng)，則Y（k）U⊥（k）=Γ（k）X（k）U⊥（k），對(duì)Y（k）U⊥（k）進(jìn)行奇異值分解，得到

式中R（k）為mM×mM正交矩陣，則R（k）前n列組成的矩陣即為該時(shí)刻的能觀性矩陣Γ（k）。

在實(shí)際計(jì)算中，通過(guò)輸入輸出數(shù)據(jù)構(gòu)建出每個(gè)時(shí)刻的Z（k）作為子空間追蹤的輸入信號(hào)。并利用奇異值分解求取某一時(shí)刻Y（k）U⊥（k）對(duì)應(yīng)的能觀性矩陣Γ（k）。將z（k）和Γ（k）分別作為子空間追蹤的輸入量x（k）和初始信號(hào)子空間W（0），這樣就完成了對(duì)系統(tǒng)輸入輸出數(shù)據(jù)的前處理過(guò)程。

2 改進(jìn)的TW-API方法

基于篇幅所限，本文只給出針對(duì)傳統(tǒng)TW-API方法所進(jìn)行修改的部分，詳細(xì)的TW-API方法推導(dǎo)過(guò)程請(qǐng)參考文［16］。該文中給出的傳統(tǒng)TW-API方法流程如表1所示。

從表1中可以看出，原文構(gòu)建了遞歸總體數(shù)據(jù)矩陣X（k）?［x（k-l+1），x（k-l+2），…，x（k）］和?V（k），通過(guò)X（k）和?V（k）的遞歸運(yùn)算來(lái)計(jì)算v（k），為此還另外定義了一個(gè)?v（k）作為狀態(tài)量v（k）的近似值，而在本文中，認(rèn)為遞歸運(yùn)算中v（k）=?v（k），即直接建立（k）與v（k）的關(guān)系，從而省略文［16］中構(gòu)建輸入/輸出總體數(shù)據(jù)矩陣X（k）和?V（k）的步驟，直接利用（k），v（k）和（k）建立子空間追蹤的關(guān)系。這樣也就不需要給定總體數(shù)據(jù)矩陣的初值X（0）和?V（0），有效地減少了計(jì)算量和計(jì)算時(shí)間。

此外，從表1中給出的每一步驟的計(jì)算量可以看出，在傳統(tǒng)TW-API算法的步驟（10）和（12）中，這兩步的計(jì)算量均為O（n3），當(dāng)系統(tǒng)的階次較高時(shí)對(duì)算法整體的計(jì)算量將有極大的影響。為此，本文考慮減少步驟（10）和（12）中有關(guān)Z（k）和Θ（k）的遞推計(jì)算量。

表1 傳統(tǒng)TW-API方法遞推步驟［16］Tab.1 The recursive step of traditional TW-API method

首先從子空間輸出信號(hào)s（k）的自相關(guān)函數(shù)Css（k）的遞推關(guān)系出發(fā)，因?yàn)閆（k）=Css（k）-1，而自相關(guān)函數(shù)的遞推關(guān)系為

利用矩陣求逆引理，即對(duì)于B=A+PJQ，有

那么利用式（13）的求逆引理對(duì)式（12）進(jìn)行求逆，可以得到

對(duì)比表1中的步驟（10）的計(jì)算量O（n3）（具體的計(jì)算量為3n2+4n3），式（14）給出的Z（k）新的遞推步驟的計(jì)算量為O（n2）（具體計(jì)算量為20+10n+5n2），從計(jì)算量的階次上可以看出很明顯地減少了計(jì)算量。

接下來(lái)考慮Θ（k）的遞推步驟的改進(jìn)，定義ε（k）為表1中誤差矩陣的平方根

從表1的步驟（12）中可以得到Θ（k）和它的轉(zhuǎn)置矩陣Θ（k）H的關(guān)系為

將式（15）代入式（17）并應(yīng)用矩陣求逆引理，可得

定義2×2的正定矩陣

這樣就得到了Θ（k）的新的遞推關(guān)系。和表1中遞推Θ（k）的步驟對(duì)比計(jì)算量，表1中原方法計(jì)算Θ（k）的計(jì)算量為4m M+O（n3），而新方法中有關(guān)Θ（k）遞推的總計(jì)算量為O（n2），從而簡(jiǎn)化了計(jì)算量。新方法與原方法相比，計(jì)算時(shí)間的優(yōu)化將在下一部分內(nèi)容中通過(guò)數(shù)值仿真來(lái)證明。改進(jìn)后的TW-API方法步驟見(jiàn)表2。

表2 改進(jìn)后的TW-API方法遞推步驟Tab.2 The recursive step of improved TW-API method

利用改進(jìn)的TW-API法計(jì)算得到各時(shí)刻的信號(hào)子空間W（k）后，離散的系統(tǒng)矩陣（k）可以通過(guò)下式得到

式中W1（k）和W2（k）分別為矩陣W（k）的前（M-1）×m行和后（M-1）×m行元素所組成的矩陣。這里（k）的上標(biāo)“?”表示辨識(shí)值，符號(hào)“?”表示偽逆。對(duì)辨識(shí)得到的?A（k）矩陣進(jìn)行特征值分解可得

式中Λ（k）=diag（λ1（k），λ2（k），…，λn（k））為對(duì)角特征根矩陣，λi（k）是時(shí)變共軛復(fù)特征根，ψ（k）為時(shí) 變特征向 量。第i階 偽特 征根λi（k）=。其中ζi（k）為系統(tǒng)的第i階偽阻尼比，ωdi（k）為系統(tǒng)的第i階含阻尼的偽固有頻率，Δt為采樣時(shí)間，圖1給出了時(shí)變模態(tài)參數(shù)辨識(shí)算法的簡(jiǎn)要流程圖。

圖1 時(shí)變模態(tài)參數(shù)辨識(shí)算法流程圖Fig.1 Flow charts of identification algorithm of time-varying modal parameters

3 仿真算例

3.1 算例1——二連桿機(jī)械臂模型

第1個(gè)數(shù)值算例選用一個(gè)二連桿機(jī)械臂的例子，結(jié)構(gòu)的示意圖如圖2所示。機(jī)械臂為均勻剛性桿，質(zhì)量均為m，長(zhǎng)度均為l。u1和u2為在關(guān)節(jié)上施加的作用力矩。關(guān)節(jié)1和關(guān)節(jié)2的阻尼為d1和d2，轉(zhuǎn)動(dòng)剛度為k1和k2。一個(gè)時(shí)變的作用力f（t）施加在連桿的自由端，f（t）與水平方向x軸的夾角為φ3。角度φ1和φ2為連桿相對(duì)于x軸的位置。當(dāng)連桿受到擾動(dòng)時(shí)，連桿將會(huì)在它們的平衡位置附近振動(dòng)，則實(shí)際角度為φ1=φ10+φ11和φ2=φ20+φ21，其中φ10和φ20為連桿與水平方向的初始夾角。假定該系統(tǒng)做微幅振動(dòng)，則系統(tǒng)的線性模型可以寫(xiě)為

圖2 二連桿機(jī)械臂模型Fig.2 Model of two-link robotic manipulator

其中，a1=

仿真中使用如下的參數(shù)：l=1 m，m=1 kg，d1=d2=0.8N·m·s/rad，k1=k2=80 N·m/rad，初始角度φ10=0°，φ20=90°，φ3=90°。作用力f（t）按以下方式變化：

其中f0=20 N，Δf=10 N。輸入u（t）為單位方差的高斯噪聲信號(hào)。狀態(tài)空間方程中的輸出矩陣直接給定為單位陣，即量測(cè)方程為y=x，采樣頻率100 Hz。定義輸出信號(hào)的信噪比（Signal Noise Ratio，SNR）為

式中σr表示原輸出信號(hào)的標(biāo)準(zhǔn)差，σnr表示含有噪聲的輸出信號(hào)標(biāo)準(zhǔn)差。在本算例中，輸入輸出Hankel矩陣中M=20，遺忘因子β=0.98。圖3給出了在信噪比SNR=50情況下的頻率辨識(shí)值。由圖3可以看出，對(duì)于這樣一個(gè)同時(shí)包含突變和周期變化特點(diǎn)的系統(tǒng)，時(shí)變模態(tài)參數(shù)可以被有效地跟蹤得到。

圖3 頻率值辨識(shí)結(jié)果的比較（SNR=50）Fig.3 Comparison of the identified results of pseudo natural frequencies values（SNR=50）

為了定量地反映本文方法的計(jì)算精度，定義辨識(shí)結(jié)果的平均絕對(duì)百分誤差（Mean Absolute Percentage Error，MAPE）eMAPE如下

式中ωk和分別為k時(shí)刻的頻率的理論值和辨識(shí)值，L為采樣點(diǎn)數(shù)。當(dāng)考慮不同的信噪比時(shí)（這里分別選取SNR=30，50和100時(shí)的情況），得到的系統(tǒng)各階頻率的MAPE誤差比較如表3所示。從表中可以看出，改進(jìn)后的TW-API方法在不同信噪比下均能滿足頻率辨識(shí)的計(jì)算精度要求。值得說(shuō)明的是，目前時(shí)變模態(tài)參數(shù)的辨識(shí)重點(diǎn)主要關(guān)注于頻率的辨識(shí)，而對(duì)于時(shí)變阻尼和振型的辨識(shí)還比較少見(jiàn)?，F(xiàn)有的子空間方法往往在其辨識(shí)中存在著較大的誤差（例如對(duì)時(shí)變阻尼的辨識(shí)，有時(shí)在噪聲的影響下可以達(dá)到30%；而某些時(shí)變系統(tǒng)還可能會(huì)有某幾階振型錯(cuò)位的現(xiàn)象），故對(duì)于阻尼和振型的辨識(shí)研究還需要做進(jìn)一步的研究工作。

表3 利用傳統(tǒng)和改進(jìn)后的TW-API方法得到的頻率值的MAPE的對(duì)比Tab.3 The MAPE comparison of frequencies by traditional and improved TW-API method

3.2 算例2——ETS-VIII衛(wèi)星模型

在第2個(gè)數(shù)值算例中，以日本的ETS-VIII衛(wèi)星為研究對(duì)象，它的基本構(gòu)型以衛(wèi)星本體為中心，一對(duì)太陽(yáng)能帆板（以下簡(jiǎn)稱(chēng)帆板n和s）以及一對(duì)巨大的拋物面通訊天線（以下簡(jiǎn)稱(chēng)天線a和b）通過(guò)硬質(zhì)連桿鉸接到衛(wèi)星本體上，具體的結(jié)構(gòu)描述請(qǐng)參見(jiàn)文［17］。

由于ETS-VIII為地球同步衛(wèi)星，為了保證太陽(yáng)能帆板朝向太陽(yáng)以獲得足夠的電力，帆板在電機(jī)的驅(qū)動(dòng)下繞衛(wèi)星的俯仰軸進(jìn)行轉(zhuǎn)動(dòng)，其旋轉(zhuǎn)角θ是周期變化的，從而導(dǎo)致衛(wèi)星的模態(tài)參數(shù)是時(shí)變的。對(duì)于這樣的航天器，一般將衛(wèi)星的中心星體作為剛體、星體所連接的附件作為柔性體來(lái)進(jìn)行考慮。則對(duì)于這樣一個(gè)呈現(xiàn)出剛?cè)狁詈咸匦缘慕Y(jié)構(gòu)，考慮到剛體的轉(zhuǎn)動(dòng)和附件的振動(dòng)，構(gòu)建對(duì)應(yīng)的約束模態(tài)模型，利用衛(wèi)星的剛體轉(zhuǎn)動(dòng)方程和柔性附件的振動(dòng)方程，建立模型的剛?cè)狁詈蟿?dòng)力學(xué)方程［17］。

式中fm為衛(wèi)星反作用飛輪產(chǎn)生的控制力矩。Ψ和J（θ）分別為整星的姿態(tài)角矢量和轉(zhuǎn)動(dòng)慣量；P（θ）為附件振動(dòng)和衛(wèi)星轉(zhuǎn)動(dòng)耦合的柔性系數(shù)耦合矩陣；μj為每個(gè)柔性附件的模態(tài)坐標(biāo)為附件的模態(tài)剛度矩陣，=diag［ω2j1，…，ω2jp］，其中ωjp為柔性附件j的第p階模態(tài)頻率；ζj為對(duì)應(yīng)的附件j的阻尼比。

將方程（28）和（29）改寫(xiě)為一般的動(dòng)力學(xué)方程形式

其中詳細(xì)的質(zhì)量矩陣等參數(shù)的形式參見(jiàn)文［17］中的表述，而方程中的狀態(tài)量q為

該衛(wèi)星系統(tǒng)中前三階模態(tài)參數(shù)表示為剛體運(yùn)動(dòng)，從第4階開(kāi)始是彈性振動(dòng)對(duì)應(yīng)的模態(tài)參數(shù)。即該動(dòng)力學(xué)方程里包含有剛體運(yùn)動(dòng)和彈性振動(dòng)兩部分，在這里只考慮辨識(shí)除剛體運(yùn)動(dòng)外的前四階彈性振動(dòng)對(duì)應(yīng)的頻率值。

在本文中考慮帆板做勻速轉(zhuǎn)動(dòng)、即旋轉(zhuǎn)角速度ωθ為常值的情況，建立帆板旋轉(zhuǎn)角θ和時(shí)間t之間的線性對(duì)應(yīng)關(guān)系θ=ωθt，那么就可以將線性變參數(shù)模型轉(zhuǎn)化為一般的線性時(shí)變模型來(lái)進(jìn)行分析。

本文在對(duì)衛(wèi)星進(jìn)行有限元建模時(shí)將原衛(wèi)星天線的殼型薄板簡(jiǎn)化為矩形平板，且天線和帆板處在同一平面上，該平面經(jīng)過(guò)衛(wèi)星本體的質(zhì)心；衛(wèi)星本體是一個(gè)規(guī)則的長(zhǎng)方體，天線和帆板各自的鉸接點(diǎn)均在衛(wèi)星本體上。為了計(jì)算得到方程（28），（29）中的整星轉(zhuǎn)動(dòng)慣量和附件的剛?cè)狁詈舷禂?shù)等方程參數(shù)，本文將總體坐標(biāo)系的原點(diǎn)取在整星的質(zhì)心，而附件坐標(biāo)系均建立在各自的質(zhì)心上。簡(jiǎn)化后的衛(wèi)星結(jié)構(gòu)如圖4所示。

圖4 簡(jiǎn)化后的衛(wèi)星示意圖Fig.4 Image of simplified satellite

在數(shù)值仿真中，輸入選取為方波激勵(lì)信號(hào)，利用Newmark法計(jì)算系統(tǒng)響應(yīng)，并在輸出中加入2%的量測(cè)白噪聲干擾。信號(hào)的采樣頻率為10 Hz，Hankel矩陣中M=15，遞推算法中的遺忘因子β= 0.98。分別利用改進(jìn)的TW-API和傳統(tǒng)TW-API方法辨識(shí)系統(tǒng)除剛體運(yùn)動(dòng)外的前四階振動(dòng)模態(tài)頻率，其中關(guān)于截?cái)啻暗拈L(zhǎng)度l的選取，通常情況下令0.5m M≤l≤2m M［16］，其中m M即為信號(hào)子空間W（k）的行數(shù)，窗的長(zhǎng)度太小的話將直接影響計(jì)算的精度。在本文中l(wèi)的長(zhǎng)度取為信號(hào)子空間W（k）的行數(shù)的3/2，即l=1.5m M。頻率值的辨識(shí)結(jié)果如圖5～9所示。

圖5 第1階頻率辨識(shí)結(jié)果Fig.5 Identification result of the first-order frequency

圖6 第2階頻率辨識(shí)結(jié)果Fig.6 Identification result of the second-order frequency

圖7 第3階頻率辨識(shí)結(jié)果Fig.7 Identification result of the third-order frequency

從圖5～9可以看出，改進(jìn)后的TW-API與傳統(tǒng)的TW-API方法相比，仍能很好地辨識(shí)系統(tǒng)的模態(tài)參數(shù)。接下來(lái)從計(jì)算時(shí)間上和原TW-API方法進(jìn)行比較，進(jìn)行10次重復(fù)試驗(yàn)，分別統(tǒng)計(jì)采樣點(diǎn)數(shù)為1×104，5×104和10×104個(gè)時(shí)的仿真計(jì)算時(shí)間。仿真的系統(tǒng)硬件條件為：Windows 8 64位系統(tǒng)、CPU處理器Intel酷睿i7 4700MQ、內(nèi)存容量8G，數(shù)值仿真使用的軟件為Matlab，計(jì)算機(jī)仿真所需要的時(shí)間對(duì)比結(jié)果見(jiàn)表4。

圖8 第4階頻率辨識(shí)結(jié)果Fig.8 Identification result of the fourth-order frequency

圖9 傳統(tǒng)TW-API方法頻率辨識(shí)結(jié)果Fig.9 Identification result of frequencies by the traditional TW-API method

從表4的仿真時(shí)間的結(jié)果可以看出，與傳統(tǒng)的TW-API方法相比，改進(jìn)后的TW-API方法可以提高約50%的計(jì)算效率，證明該方法的遞推步驟的計(jì)算量簡(jiǎn)化是有效的。再利用文［8］和［15］中分別提出的經(jīng)典的PAST和指數(shù)窗API方法，在相同條件下對(duì)該算例進(jìn)行辨識(shí)，各階頻率的辨識(shí)結(jié)果見(jiàn)圖10 和11。上述4種方法得到的系統(tǒng)各階頻率的MAPE誤差比較如表5所示。

圖10 PAST方法的頻率辨識(shí)結(jié)果Fig.10 Identification result of frequencies by the PAST method

圖11 API方法頻率辨識(shí)結(jié)果Fig.11 Identification result of frequencies by the API method

表4 傳統(tǒng)和改進(jìn)后的TW-API方法的計(jì)算時(shí)間對(duì)比Tab.4 The computation time comparison of traditional and improved TW-API method

從表5可以看出，本文提出的改進(jìn)TW-API方法由于采用了截?cái)啻靶问剑耘c指數(shù)窗形式的PAST和API方法相比，誤差稍有增加，但誤差的精度量級(jí)仍能很好地滿足辨識(shí)要求；而與傳統(tǒng)的TW-API方法相比，保留了原有的計(jì)算精度。將表4和5的計(jì)算結(jié)果結(jié)合起來(lái)分析，可以看出改進(jìn)的TW-API方法在保證了原有計(jì)算精度的前提下大幅減少了計(jì)算量，證明本文提出的改進(jìn)TW-API方法在辨識(shí)系統(tǒng)的動(dòng)力學(xué)參數(shù)方面是有效的。

4 結(jié) 論

本文從截?cái)啻氨平鼉绲涌臻g方法出發(fā)，對(duì)現(xiàn)有的TW-API方法進(jìn)行了簡(jiǎn)化和改進(jìn)，使得在保持原有逼近冪迭代方法的投影子空間正交性的同時(shí)，不再需要構(gòu)建輸入輸出總體數(shù)據(jù)矩陣用于迭代運(yùn)算；同時(shí)從信號(hào)子空間遞推理論出發(fā)，通過(guò)建立新的數(shù)據(jù)矩陣的遞推關(guān)系，成功地減少了矩陣遞推運(yùn)算中需要的計(jì)算量，從而有效地減少了運(yùn)算時(shí)間。通過(guò)二連桿機(jī)械臂和ETS-VIII衛(wèi)星模型的例子，采用本文改進(jìn)的TW-API方法進(jìn)行時(shí)變模態(tài)參數(shù)的辨識(shí)。從辨識(shí)結(jié)果可以看出改進(jìn)的TW-API方法與原來(lái)的方法相比，計(jì)算時(shí)間減少了約50%，而且仍然能夠保證計(jì)算的精度。在仿真算例中還與常用的PAST和指數(shù)窗API遞推辨識(shí)方法進(jìn)行了比較，結(jié)果證明本文提出的算法能夠有效地辨識(shí)時(shí)變系統(tǒng)的模態(tài)參數(shù)。

［1］ Liu K.Identification of linear time-varying systems ［J］.Journal of Sound and Vibration，1997，206（4）：487—505.

［2］ Liu K.Extension of modal analysis to linear time-varying systems［J］.Journal of Sound and Vibration，1999，226（1）：149—167.

［3］ 于開(kāi)平，謝禮立，樊久銘，等.移動(dòng)質(zhì)量簡(jiǎn)支梁系統(tǒng)的參數(shù)辨識(shí)［J］.地震工程與工程振動(dòng)，2002，22（5）：14—17. Yu Kaiping，Xie Lili，F(xiàn)an Jiuming，et al.A parameter identification of simply supported beams system carrying a moving mass［J］.Earthquake Engineering and Engineering Vibration，2002，22（5）：14—17.

［4］ 李會(huì)娜，史治宇.基于自由響應(yīng)數(shù)據(jù)的時(shí)變系統(tǒng)物理參數(shù)識(shí)別［J］.振動(dòng)工程學(xué)報(bào)，2007，20（4）：348—352. Li Huina，Shi Zhiyu.Physical parameter identification of time-varying system based on free response data ［J］.Journal of Vibration Engineering，2007，20（4）：348—352.

［5］ 李會(huì)娜，史治宇.基于隨機(jī)激勵(lì)響應(yīng)的時(shí)變系統(tǒng)物理參數(shù)子空間識(shí)別方法研究［J］.工程力學(xué)，2008，25 （9）：46—51. Li Huina，Shi Zhiyu.Subspace-based physical parameter identification of time-varying system using random excitation response［J］.Engineering Mechanics，2008，25（9）：46—51.

［6］ Tasker F，Bosse A，F(xiàn)isher S.Real-time modal parameter estimation using subspace methods：theory［J］. Mechanical Systems and Signal Processing，1998，12 （6）：797—808.

［7］ Bosse A，Tasker F，F(xiàn)isher S.Real-time modal parameter estimation using subspace methods：applications［J］.Mechanical Systems and Signal Processing，1998，12（6）：809—823.

［8］ Yang B.Projection approximation subspace tracking ［J］.IEEE Transactions on Signal Processing，1995，43（1）：95—107.

［9］ Yang B.Asymptotic convergence analysis of the projection approximation subspace tracking algorithms ［J］.Signal Processing，1996，50：123—136.

［10］龐世偉，于開(kāi)平，鄒經(jīng)湘.用于線性時(shí)變系統(tǒng)辨識(shí)的固定長(zhǎng)度平移窗投影估計(jì)遞推子空間方法［J］.機(jī)械工程學(xué)報(bào)，2005，41（10）：117—122. Pang Shiwei，Yu Kaiping，Zou Jingxiang.Using recursive subspace method based on projection approximation with moving window to estimate parameters of linear time-varying structure［J］.Chinese Journal of Mechanical Engineering，2005，41（10）：117—122.

［11］龐世偉，于開(kāi)平，鄒經(jīng)湘.用于時(shí)變結(jié)構(gòu)模態(tài)參數(shù)識(shí)別的投影估計(jì)遞推子空間方法 ［J］.工程力學(xué)，2005，22（5）：115—119. Pang ShiWei，Yu Kaiping，Zou Jingxiang.A projection approximation recursive subspace method for identification of modal parameters of time-varying structures［J］.Engineering Mechanics，2005，22（5）：115—119.

［12］Real E C，Tufts D W，Cooley J W.Two algorithms for fast approximate subspace tracking［J］.IEEE Transactions on Signal Processing，1999，47（7）：1 936—1 945.

［13］Hua Y B，Xiang Y，Chen T，et al.A new look at the power method for fast subspace tracking［J］.Digital Signal Processing，1999，9（4）：297—314.

［14］Hua Y B，Xiang Y，Chen T，et al.Natural power method for fast subspace tracking［A］.Neural Networks for Signal Processing IX，1999.Proceedings of the 1999 IEEE Signal Processing Society Workshop. IEEE［C］.Madison，USA.1999：176—185.

［15］Badeau R，Richard G，David B，et al.Approximated power iterations for fast subspace tracking［A］.Signal Processing and Its Applications，2003.Proceedings. Seventh International Symposium on.IEEE［C］.Paris，F(xiàn)rance.2003，2：583—586.

［16］Badeau R，David B，Richard G.Fast approximated power iteration subspace tracking［J］.IEEE Transactions on Signal Processing，2005，53（8）：2 931—2 941.

［17］Hamada Y，Ohtani T，Kida T，et al.Synthesis of a linearly interpolated gain scheduling controller for large flexible spacecraft ETS-VIII［J］.Control Engineering Practice，2011，19：611—625.

（1.State Key Laboratory of Structural Analysis for Industrial Equipment，Dalian University of Technology，Dalian 116024，China；2.School of Aeronautics and Astronautics，Dalian University of Technology，Dalian 116024，China）

An improved TW-API recursive method and its application on time-varying modal parameters identification

NI Zhi-yu1，TAN Shu-jun2，WU Zhi-gang1，2

An improved method based on Truncated Window Approximated Power Iteration（TW-API）Tracking for modal parameter identification is proposed in this paper.The truncation window cross-correlation function of input and output signals is adopted to establish the recursion relation of signal dominant subspace，and subspace projection is then used to identify the modal parameters.Compared with the traditional TW-API method，the improved one，which is based on the theory of signal subspace，simplifies the recursion procedure of data matrix，significantly reduces the computing cost and computation time，and meanwhile guarantees the calculation accuracy.In numerical simulation，a structure of two-link robotic manipulator and a model of on-orbit spacecraft are used to identify the time-varying modal parameters.Moreover，the computation results of these two methods are also compared with the classical Projection Approximation Subspace Tracking（PAST）and Approximated Power Iteration（API）recursive methods.The results of the four identification methodologies demonstrate that the proposed new algorithm in this paper can effectively identify the time-varying modal parameters of the system.

modal parameter identification；time-varying system；recursive subspace method；spacecraft

TB123；N945.14

A

1004-4523（2015）05-0721-09

10.16385/j.cnki.issn.1004-4523.2015.05.006

倪智宇（1985—），男，博士研究生。電話：（0411）84706772；E-mail：nizhiyu-8597@mail.dlut.edu.cn

吳志剛（1971—），男，博士，教授。電話：（0411）84706772；E-mail：wuzhg@dlut.edu.cn

2014-04-03

2015-05-11

國(guó)家自然科學(xué)基金資助項(xiàng)目（11072044，11372056，11002032）；高等學(xué)校博士學(xué)科點(diǎn)專(zhuān)項(xiàng)科研基金資助項(xiàng)目（20110041130001）