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

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

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

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

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

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

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

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

引 言

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

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

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

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

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

1.1 信號子空間中輸入的迭代更新

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

其中

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

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

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

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

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

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

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

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

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

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

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

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

基于篇幅所限，本文只給出針對傳統(tǒng)TW-API方法所進(jìn)行修改的部分，詳細(xì)的TW-API方法推導(dǎo)過程請參考文［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），通過X（k）和?V（k）的遞歸運(yùn)算來計算v（k），為此還另外定義了一個?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），有效地減少了計算量和計算時間。

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

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

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

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

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

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

接下來考慮Θ（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）的步驟對比計算量，表1中原方法計算Θ（k）的計算量為4m M+O（n3），而新方法中有關(guān)Θ（k）遞推的總計算量為O（n2），從而簡化了計算量。新方法與原方法相比，計算時間的優(yōu)化將在下一部分內(nèi)容中通過數(shù)值仿真來證明。改進(jìn)后的TW-API方法步驟見表2。

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

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

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

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

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

3 仿真算例

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

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

圖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）為單位方差的高斯噪聲信號。狀態(tài)空間方程中的輸出矩陣直接給定為單位陣，即量測方程為y=x，采樣頻率100 Hz。定義輸出信號的信噪比（Signal Noise Ratio，SNR）為

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

本文從截斷窗逼近冪迭代子空間方法出發(fā)，對現(xiàn)有的TW-API方法進(jìn)行了簡化和改進(jìn)，使得在保持原有逼近冪迭代方法的投影子空間正交性的同時，不再需要構(gòu)建輸入輸出總體數(shù)據(jù)矩陣用于迭代運(yùn)算；同時從信號子空間遞推理論出發(fā)，通過建立新的數(shù)據(jù)矩陣的遞推關(guān)系，成功地減少了矩陣遞推運(yùn)算中需要的計算量，從而有效地減少了運(yùn)算時間。通過二連桿機(jī)械臂和ETS-VIII衛(wèi)星模型的例子，采用本文改進(jìn)的TW-API方法進(jìn)行時變模態(tài)參數(shù)的辨識。從辨識結(jié)果可以看出改進(jìn)的TW-API方法與原來的方法相比，計算時間減少了約50%，而且仍然能夠保證計算的精度。在仿真算例中還與常用的PAST和指數(shù)窗API遞推辨識方法進(jìn)行了比較，結(jié)果證明本文提出的算法能夠有效地辨識時變系統(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］ 于開平，謝禮立，樊久銘，等.移動質(zhì)量簡支梁系統(tǒng)的參數(shù)辨識［J］.地震工程與工程振動，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］ 李會娜，史治宇.基于自由響應(yīng)數(shù)據(jù)的時變系統(tǒng)物理參數(shù)識別［J］.振動工程學(xué)報，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］ 李會娜，史治宇.基于隨機(jī)激勵響應(yīng)的時變系統(tǒng)物理參數(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］龐世偉，于開平，鄒經(jīng)湘.用于線性時變系統(tǒng)辨識的固定長度平移窗投影估計遞推子空間方法［J］.機(jī)械工程學(xué)報，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］龐世偉，于開平，鄒經(jīng)湘.用于時變結(jié)構(gòu)模態(tài)參數(shù)識別的投影估計遞推子空間方法 ［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

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