利用改進(jìn)TW-API方法在軌辨識撓性航天器時變模態(tài)參數(shù)

倪智宇，鄔樹楠，，吳志剛，，譚述君

(1.大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，大連116024;2.大連理工大學(xué)航空航天學(xué)院，大連116024)

0 引言

航天工程中許多問題需要考慮系統(tǒng)的時變特性，例如大型空間結(jié)構(gòu)裝配、太陽能帆板展開與旋轉(zhuǎn)等，這些動態(tài)過程造成系統(tǒng)模態(tài)參數(shù)發(fā)生改變。隨著航天技術(shù)的進(jìn)展，大型化與撓性化成為航天器結(jié)構(gòu)的一個重要發(fā)展趨勢。航天器撓性化導(dǎo)致固有頻率值很低，剛?cè)狁詈犀F(xiàn)象明顯。對于大型剛?cè)狁詈峡臻g結(jié)構(gòu)，通過有限元數(shù)值計算和地面動力學(xué)試驗辨識得到的參數(shù)必然存在著一定的誤差，甚至由于結(jié)構(gòu)過大過于復(fù)雜而無法進(jìn)行地面試驗。對于具有時變動力學(xué)特性的大型航天器，利用在軌飛行中的輸入輸出測量數(shù)據(jù)，采用在軌辨識算法獲取航天器模態(tài)參數(shù)隨時間的變化特性具有重要意義。

在航天器在軌辨識試驗中，往往利用飛輪或姿控推力器產(chǎn)生激勵信號，通過測量得到的姿態(tài)角和角速度，以及星體姿態(tài)運(yùn)動和大型附件振動的耦合關(guān)系來進(jìn)行航天器動力學(xué)參數(shù)辨識。為更精確地測量航天器的動態(tài)響應(yīng)，選擇在太陽能帆板上布設(shè)傳感器來獲得加速度信號用于參數(shù)辨識。但是如何產(chǎn)生有效的輸入信號以保證安全而有效地激發(fā)航天器的固有模態(tài)特性，以及如何在帆板上合理布置傳感器以獲得所需的輸出信號(實際在軌辨識中，傳感器數(shù)目有限)，即作動器的激振和傳感器配置優(yōu)化問題［1］，這些都是航天器在軌參數(shù)辨識研究中的重要課題和難點，需要進(jìn)一步研究。在實際的航天器在軌辨識試驗中，Juang等采用ERA系列方法［2］，使用控制力矩輸入信號和陀螺儀量測信號數(shù)據(jù)，對哈勃太空望遠(yuǎn)鏡進(jìn)行了系統(tǒng)頻率、阻尼等模態(tài)參數(shù)的在軌辨識。文獻(xiàn)［3］利用星上助推火箭產(chǎn)生的隨機(jī)激勵信號，測量帆板的加速度信號和星體的姿態(tài)信號，對ETS-VI衛(wèi)星進(jìn)行了系統(tǒng)的在軌辨識試驗和魯棒控制試驗。而Kasai等［4］對ETS-VIII衛(wèi)星的在軌模態(tài)參數(shù)進(jìn)行了一系列的辨識試驗。

在航天器模態(tài)參數(shù)辨識中，子空間系列方法得到了廣泛發(fā)展和應(yīng)用［5］。在各類子空間方法中，基于信號處理技術(shù)提出的遞推子空間方法由于具有較快的計算速度，在航天器動力學(xué)參數(shù)辨識方面具有一定優(yōu)勢。Yang［6］根據(jù)信號子空間投影原理，提出了投影估計子空間跟蹤(Projection Approximation Subspace Tracking，PAST)方法。龐世偉等［7］利用該方法辨識了移動質(zhì)量－簡支梁的時變模態(tài)參數(shù)。Badeau等［8］針對PAST方法信號子空間列向量不能嚴(yán)格正交的問題，提出了逼近冪迭代(Approximation Power Iteration，API)方法，對信號子空間投影乘以修正后的正交矩陣，使其投影子空間保持嚴(yán)格正交的同時可以跟蹤快變的主子空間。與利用整體數(shù)據(jù)的子空間辨識方法相比，遞推辨識方法不需要逐一時刻進(jìn)行SVD分解計算，而是通過信號子空間迭代以避免較高維數(shù)的Hankel矩陣求解問題，從而減小了計算量。

傳統(tǒng)截斷窗逼近冪迭代(Truncated Window-Approximation Power Iteration，TW-API)遞推方法的計算量可以達(dá)到O(n3)，其中n為系統(tǒng)的階次。對于航天器在軌辨識而言，過長的計算時間必然會影響辨識效率。而本文簡化了數(shù)據(jù)遞推步驟中的計算，使得這種方法的計算量大幅減少，改進(jìn)后算法的總計算量只有 O(n2)，本文還對另外兩種經(jīng)典PAST和API遞推方法的計算量進(jìn)行了總結(jié)。為對比這幾種經(jīng)典遞推方法在航天器在軌辨識中的計算速度和計算精度，選取ETS-VIII衛(wèi)星作為數(shù)值仿真對象，辨識太陽能帆板轉(zhuǎn)動所造成的整星模態(tài)參數(shù)時變特性。結(jié)果表明本文提出的改進(jìn)TW-API方法與現(xiàn)有大多數(shù)方法相比，在計算時間上具有一定的優(yōu)勢，對于航天器在軌辨識有一定的實用價值。

1 航天器建模與輸入/輸出數(shù)據(jù)的前處理

1. 1 航天器約束模態(tài)姿態(tài)模型

常見的航天結(jié)構(gòu)如衛(wèi)星和空間站等，往往將其作為剛?cè)狁詈辖Y(jié)構(gòu)來建模。基于航天器剛體轉(zhuǎn)動方程和撓性附件振動方程，并且在剛?cè)狁詈舷禂?shù)中只考慮附件振動對剛體轉(zhuǎn)動的影響，建立約束模態(tài)剛?cè)狁詈蟿恿W(xué)方程［9］:

式中:θ為附件繞自身縱向中心軸旋轉(zhuǎn)的旋轉(zhuǎn)角;fm為航天器反作用飛輪產(chǎn)生的控制力矩。Ψ和J(θ)分別為整體結(jié)構(gòu)的姿態(tài)角矢量和轉(zhuǎn)動慣量;Pj(θ)為第j個附件振動和中心剛體轉(zhuǎn)動耦合的撓性系數(shù)耦合矩陣;μj為第j個撓性附件的模態(tài)坐標(biāo);Ω2j為中ωjp為撓性附件j的第p階模態(tài)頻率;ζj為對應(yīng)的附件j的阻尼比。

將動力學(xué)方程(1)和(2)改寫為如下形式:

其中，矩陣M、E和K的詳細(xì)形式見文獻(xiàn)［10］。方程中的狀態(tài)量q和系數(shù)矩陣L分別為

式中:矩陣L中的零矩陣維數(shù)為3pN，p為上文中提到的每個附件所選取的模態(tài)階數(shù)，而N為航天器結(jié)構(gòu)中所包含的附件數(shù)目。

從方程(1)和(2)可以看出，所構(gòu)建的動力學(xué)方程(3)里包含有結(jié)構(gòu)剛體運(yùn)動和彈性振動兩部分。姿態(tài)控制力矩通過剛?cè)狁詈舷禂?shù)矩陣對系統(tǒng)振動方程產(chǎn)生影響。因此本文利用飛輪產(chǎn)生的力矩作為激勵信號，系統(tǒng)的姿態(tài)角、姿態(tài)角速度以及帆板上的傳感器測得的附件振動響應(yīng)作為輸出信號進(jìn)行時變系統(tǒng)模態(tài)參數(shù)的辨識。利用飛輪力矩作為輸入信號的優(yōu)點是可以使用衛(wèi)星正常運(yùn)行的數(shù)據(jù)來進(jìn)行辨識試驗。例如ETS-VIII衛(wèi)星的在軌辨識試驗就將穩(wěn)態(tài)運(yùn)行階段的姿態(tài)機(jī)動信號作為辨識的輸入信號，從而有利于衛(wèi)星姿態(tài)保持，避免額外的燃料消耗。需要說明的是，本文只選用了飛輪的控制力矩作為輸入信號，暫未考慮其它激勵條件如熱致振動或沖擊輸入的影響。

本文在考慮附件旋轉(zhuǎn)角速度ωθ為常值的情況下，建立附件旋轉(zhuǎn)角θ和時間t之間的線性對應(yīng)關(guān)系θ=ωθt，從而將方程(3)中的線性變參數(shù)模型轉(zhuǎn)化為線性時變模型進(jìn)行分析。

1. 2 航天器輸入/輸出信號的前處理

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

首先考慮將方程(3)改寫為如下離散狀態(tài)方程形式

式中:系統(tǒng)階數(shù)為n，輸入維數(shù)為r，輸出維數(shù)為m。任意k時刻的輸入矩陣U(k)和輸出矩陣Y(k)以及下一時刻輸入數(shù)據(jù)ˉu(k+1)和輸出數(shù)據(jù)ˉy(k+1)的Hankel矩陣形式分別寫為

式中:

其中，M為合適的Hankel矩陣維數(shù)。根據(jù)數(shù)據(jù)更新的 URV 遞歸方法［11］，令

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

式中:“+”表示偽逆。在式(7)和(9)中，下一個時刻k的(U(k)UT(k))－1和Y(k)U+(k)的更新方式為

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

2 用于參數(shù)辨識的兩種經(jīng)典子空間遞推方法

2. 1 PAST 方法

Yang基于信號子空間原理提出了PAST方法［6］。PAST方法總計算量為3mMn+O(n2)，計算量較小，這在航天工程中是比較適用的。但是該方法不能保證所得到的子空間列向量完全投影正交，所以PAST方法對于辨識快變系統(tǒng)的參數(shù)存在不足。為解決這種問題，文獻(xiàn)［8］提出API系列方法。

2. 2 API方法

相比較PAST方法，API方法(基于無限指數(shù)窗)在子空間遞推方面增加了正交項，以保證投影矩陣的正交性，同時還加強(qiáng)了對快變子空間的跟蹤性能，因此對于時變程度較快的系統(tǒng)仍能較好地跟蹤和辨識［8］。

API方法計算量為mM(n2+3n+1)+O(n3)，高于PAST方法，但是系統(tǒng)的正交性和對快變子空間的跟蹤性能更好。無論是PAST還是API方法，均是基于無限指數(shù)窗形式，對于系統(tǒng)信號發(fā)生突變的情況不能夠很好跟蹤，為此又在原有方法的基礎(chǔ)上，提出了基于滑移窗的截斷窗逼近冪迭代遞推方法(TW-API)。

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

本文只給出針對傳統(tǒng)TW-API方法所進(jìn)行的改進(jìn)，該方法的具體推導(dǎo)過程見文獻(xiàn)［8］。傳統(tǒng)TWAPI方法計算量為6mMn+mMn2+n2l+O(n3)，其中l(wèi)為選取的截斷窗的長度，可看出計算量明顯高于指數(shù)窗形式的API方法。在原TW-API方法中，復(fù)迭代的步驟。另外，傳統(tǒng)TW-API算法中，在更新矩陣Z(k)和正交矩陣Θ(k)這兩步驟中的計算量均為O(n3)，由于航天器系統(tǒng)的階次往往普遍較高，這對方法整體的計算量將會有極大的影響。為此，本文考慮減少原TW-API方法中有關(guān)矩陣Z(k)和Θ(k)的遞推計算量(Z(k)和Θ(k)的原表達(dá)式見文獻(xiàn)［8］的表2)。

首先從子空間輸出信號 s－(k)的自相關(guān)函數(shù)Css(k)的遞推關(guān)系出發(fā)，因為Z(k)=(k)，而自相關(guān)函數(shù)的遞推關(guān)系為

根據(jù)矩陣求逆引理，對于B=A+PJQ，存在有

那么利用式(14)的求逆引理對式(13)進(jìn)行求逆，并定義一個臨時變量h(k)(h(k)的具體形式見表1)可以得到

原TW-API方法中Z(k)的計算量為3n2+4n3，而式(15)給出的Z(k)新的遞推步驟的計算量為20+10n+5n2，從階次上可以看出明顯地減少了計算量。

接下來考慮Θ(k)的遞推步驟的改進(jìn)，定義ε－(k)為文獻(xiàn)［8］里原誤差矩陣 e－H(t)e－(t)的平方根

則有

利用文獻(xiàn)［8］中Θ(k)的表達(dá)式，可以得到

式中:In為單位矩陣，g(k)的具體形式見表1，將式(16)代入式(18)并應(yīng)用矩陣求逆引理，可得

其中，ρ－(k)為2×2的正定矩陣:

考慮式(19)，尋找一個特解σ－(k)使其滿足:

那么特解 σ－(k)= ρ－(k)+(ρ－H(k))1/2。

定義2×2的正定矩陣τ－(k)為

則

這樣就得到了Θ(k)的新的遞推關(guān)系。原TW-API方法計算Θ(k)的計算量為4mM+O(n3)，而新方法中Θ(k)遞推步驟的總計算量為O(n2)，同樣簡化了計算量。改進(jìn)后的TW-API方法的計算量為mM(n2+6n+4)+O(n2)，改進(jìn)后的TW-API方法的遞推步驟見表1，其中遞推計算中的初始條件選為W(0)=過程中，由第1.2節(jié)中生成的各個時刻的z(k)作為各個時刻的輸入值x(k)。表2則給出了以上提到的四種子空間遞推方法的計算量和適用范圍的總結(jié)。

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

表2 四種方法的復(fù)雜性的比較Table 2 Comparison of the four algorithms'complexities

利用子空間遞推方法計算得到各時刻的信號子空間W(k)后，離散的系統(tǒng)矩陣A^(k)可以通過下式得到

式中:W1(k)和 W2(k)分別為矩陣 W(k)的前(M－1)×m行和后(M－1)×m行元素所組成的矩陣。矩陣上標(biāo)“^”表示辨識值。

根據(jù)時變系統(tǒng)的偽模態(tài)理論，即認(rèn)為在很小的采樣時間內(nèi)，系統(tǒng)的特征根是不變的，對應(yīng)于該時刻的特征根稱其為“偽”特征根，將這種利用定常系統(tǒng)的理論來描述得到時變系統(tǒng)模態(tài)參數(shù)的方法稱為“偽”模態(tài)分析方法。從航天工程的角度來說，只要采樣間隔足夠小，那么偽特征根對于系統(tǒng)的時變情況的反映是近似而合理的。

對辨識得到的A^(k)矩陣進(jìn)行特征值分解可得

式中:Λ(k)=diag(λ1(k)，…，λn(k))為對角特征根矩陣，ψ(k)為時變特征向量，λi(k)是時變共軛復(fù)特征根。第i階偽特征根λi(k)=exp(－ζi(k)Δt± iωdi(k)Δt)，ζi(k)為系統(tǒng)第 i階偽阻尼比，ωdi(k)為系統(tǒng)的第i階含阻尼的偽固有頻率，Δt為采樣時間。

4 數(shù)值實例

4. 1 衛(wèi)星模型的建立與簡化

本文在數(shù)值算例中，以日本的ETS-VIII衛(wèi)星為研究對象，它的基本構(gòu)型以衛(wèi)星本體為中心，一對太陽能帆板(以下簡稱帆板n和s)以及一對巨大的拋物面通訊天線(以下簡稱天線a和b)通過硬質(zhì)連桿鉸接到衛(wèi)星本體上［10］。利用前文提到的基于約束模態(tài)模型的航天器剛?cè)狁詈蟿恿W(xué)方程，建立該衛(wèi)星的狀態(tài)空間方程。

本文在對衛(wèi)星進(jìn)行有限元建模時將天線的殼型薄板簡化為矩形平板，且天線和帆板處在同一平面上，該平面經(jīng)過衛(wèi)星本體的質(zhì)心;衛(wèi)星本體簡化為長方體，而附件坐標(biāo)系均建立在各自的質(zhì)心上，總體坐標(biāo)系的原點取在整星的質(zhì)心上。簡化后的衛(wèi)星結(jié)構(gòu)如圖1所示。

由于ETS-VIII為地球同步衛(wèi)星，帆板在電機(jī)的驅(qū)動下繞衛(wèi)星的俯仰軸進(jìn)行轉(zhuǎn)動以保證始終朝向太陽以獲得足夠的電力，其旋轉(zhuǎn)角θ是周期變化的，從而導(dǎo)致衛(wèi)星的模態(tài)參數(shù)是時變的，令帆板旋轉(zhuǎn)角速度為0.6(°)/s。該系統(tǒng)中前三階為剛體運(yùn)動對應(yīng)的模態(tài)參數(shù)，從第四階開始對應(yīng)彈性振動。

圖1 簡化后的衛(wèi)星示意圖Fig.1 Simplified image of satellite

4. 2 仿真辨識結(jié)果

選取采樣時間為Δt=0.1 s，進(jìn)行10次隨機(jī)試驗，Hankel矩陣中M=20，仿真時在輸出中加入2%的量測白噪聲干擾。分別利用改進(jìn)TW-API、傳統(tǒng)TW-API、PAST和指數(shù)窗API四種辨識方法進(jìn)行衛(wèi)星系統(tǒng)振動模態(tài)頻率的辨識。關(guān)于截斷窗長度l的選取，通常情況下令0.5mM ≤ l≤2mM［8］，其中mM即為矩陣W(k)的行數(shù)，窗的長度太小的話將直接影響計算的精度。在本文中l(wèi)的長度取為l=1.5mM。這里均以系統(tǒng)第4階頻率值的辨識結(jié)果為例，辨識結(jié)果如圖2～圖5所示。

圖2 改進(jìn)TW-API方法的頻率辨識結(jié)果Fig.2 Identification results of frequencies by the improved TW-API method

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

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

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

從圖2～圖5可以看出，這幾種方法均能很好地辨識衛(wèi)星的模態(tài)參數(shù)。接下來從計算時間上將四種方法進(jìn)行比較，為此分別進(jìn)行10次蒙特卡羅試驗，其中統(tǒng)計采樣點數(shù)為10，000個，計算所需要的仿真時間。仿真的系統(tǒng)硬件條件為:Windows 8 64位系統(tǒng)、CPU處理器Intel酷睿i7 4700MQ、內(nèi)存容量8G，數(shù)值仿真使用的軟件為Matlab，各辨識方法在計算機(jī)仿真中所需要的平均計算時間見表3。

表3 四種遞推方法的計算時間Table 3 The computation time of the four recursive methods

從表3可以看出，四種方法所需要的計算時間:PAST＜改進(jìn)TW-API＜API＜傳統(tǒng)TW-API。這與之前計算量分析的結(jié)論相吻合。與傳統(tǒng)的TW-API方法相比，改進(jìn)后的TW-API方法可以提高近50%的計算效率，雖然在計算時間上仍高于PAST方法，但是低于指數(shù)窗形式的API方法，證明該方法有效簡化了計算量。

為定量比較這幾種遞推方法在航天器辨識中的計算精度，定義辨識結(jié)果的平均絕對百分誤差(Mean Absolute Percentage Error，MAPE)如下

從表4可以看出，本文提出的改進(jìn)TW-API方法由于采用截斷窗形式，所以與采用指數(shù)窗形式的PAST和API方法相比，誤差稍有增加，但誤差的精度量級仍能很好地滿足辨識要求。綜合分析表3和表4給出的數(shù)據(jù)結(jié)果，可以看出改進(jìn)的TW-API方法在保證計算精度的前提下，大幅減少了遞推過程中所需的計算量，證明本文提出的方法在航天器系統(tǒng)的模態(tài)參數(shù)辨識方面是有效的。

表4 頻率的MAPE的比較Table 4 Comparison of the frequencies'MAPE

在數(shù)值仿真中，如果令衛(wèi)星的帆板轉(zhuǎn)速增加為6(°)/s，即令系統(tǒng)成為一個快速變化的系統(tǒng)，比較這幾種方法對于快變系統(tǒng)的辨識效果。圖6給出了第四階頻率前60 s的辨識結(jié)果，從圖中可以看出，當(dāng)時變系統(tǒng)變化較快時，API/TW-API系列方法比較準(zhǔn)確地跟蹤了系統(tǒng)的頻率值，而PAST方法對于快變系統(tǒng)的跟蹤性能有所下降。特別是在曲線的拐點處，該方法得到的頻率值與API/TW-API方法相比，存在著比較大的誤差。

圖6 四種方法的頻率辨識結(jié)果Fig.6 Identification result of frequency by the four methods

在該算例的最后，我們考慮在航天器時變模態(tài)參數(shù)在軌辨識時，飛輪實際施加的力矩輸入信號值與理論值存在誤差的情況。在實際工程應(yīng)用中，飛輪通過控制轉(zhuǎn)速來得到相應(yīng)的力矩，但是由于飛輪容易存在輸出擾振的情況，從而導(dǎo)致獲得的實際輸入信號值與想要的力矩值存在一定的計算誤差。這里將這種誤差作為干擾信號加入到原輸入信號中，分別考慮具有如下兩種輸入噪聲干擾的情況:1)原方波輸入信號u(k)中含有1%白噪聲干擾;2)方波輸入信號u(k)中含有大小為0.1×u(k)的測量誤差。

圖7 輸入信號含有1%白噪聲干擾時的第四階頻率結(jié)果Fig.7 The result of the fourth-order frequency with the 1%white noise in the input signal

圖8 輸入信號含有0.1×u(k)噪聲干擾時第四階頻率結(jié)果Fig.8 The result of the fourth-order frequency with the noise 0.1 × u(k)in the input signal

在這兩種情況下利用本文的改進(jìn)TW-API方法對航天器時變模態(tài)參數(shù)(以系統(tǒng)第4階頻率為例)的辨識結(jié)果如圖7～圖8所示。從圖7～圖8中的結(jié)果可以看出，當(dāng)輸入信號存在一定的白噪聲干擾或是測量上的誤差時，撓性航天器的時變模態(tài)參數(shù)仍能夠有效地得到辨識。

5 結(jié)論

對于大型航天器模態(tài)參數(shù)的在軌辨識，所選方法的計算量和計算速度是工程技術(shù)人員需要考慮的一個重點。本文提出的改進(jìn)的TW-API方法，通過遞推步驟的改進(jìn)，大大提高了這種遞推方法的計算效率。同時還將改進(jìn)的TW-API方法與傳統(tǒng)的TWAPI方法以及經(jīng)典的PAST和API遞推方法進(jìn)行了計算效率的比較，詳細(xì)分析了這四種方法的計算效率和適用條件，結(jié)果證明本文提出的方法與其他三種方法相比，具有很好的辨識計算速度和精度，數(shù)值仿真結(jié)果也證明了這一點。

［1］ 徐敏強(qiáng)，宋其江，王日新.基于可觀測性和可靠性的傳感器分布優(yōu)化設(shè)計［J］.宇航學(xué)報，2010，31(11):2618-2622.［Xu Min-qiang，Song Qi-jiang，Wang Ri-xin.Optimization design of sensors location based on fault observability and reliability［J］.Journal of Astronautics，2010，31(11):2618 －2622.］

［2］ Juang JN，Phan M，Horta L G，et al.Identification of observer/Kalman filter Markov parameters-theory and experiments［J］.Journal of Guidance，Control，and Dynamics，1993，16(2):320－329.

［3］ Adachi S，Yamaguchi I，Kida T，et al.On-orbit system identification experiments on engineering test satellite-VI［J］.Control Engineering Practice，1999，7(7):831 －841.

［4］ Kasai T，Yamaguchi I，Igawa H，et al.On-orbit system identification experiments of the engineering test satellite-VIII［J］.Transactions of Space Technology Japan，2009，7(26):79－84.

［5］ 孔憲仁，張紅亮，張也弛.用于統(tǒng)計能量分析參數(shù)辨識的子空間方法研究［J］.宇航學(xué)報，2011，32(12):2471－2477.［Kong Xian-ren，Zhang Hong-liang，Zhang Ye-chi.Study on subspace method applied to SEA parameter identification［J］.Journal of Astronautics，2011，32(12):2471 －2477.］

［6］ Yang B.Projection approximation subspace tracking［J］.IEEE Transactions on Signal Processing，1995，43(1):95－107.

［7］ 龐世偉，于開平，鄒經(jīng)湘.用于時變結(jié)構(gòu)模態(tài)參數(shù)識別的投影估計遞推子空間方法［J］.工程力學(xué)，2005，22(5):115－119.［Pang Shi-wei，Yu Kai-ping，Zou Jing-xiang.A projection approximation recursive subspace method for identification of modal parameters of time-varying structures［J］.Engineering Mechanics，2005，22(5):115 －119.］

［8］ Badeau R，David B，Richard G.Fast approximated power iteration subspace tracking［J］.IEEE Transactions on Signal Processing，2005，53(8):2931-2941.

［9］ 史紀(jì)鑫，曲廣吉.可變構(gòu)型復(fù)合柔性結(jié)構(gòu)航天器動力學(xué)建模研究［J］.宇航學(xué)報，2007，28(1):130-135.［Shi Ji-xin，Qu Guang-ji.Mathematical modeling of a class of variable structure spacecraft with flexible multibody appendages［J］.Journal of Astronautics，2007，28(1):130 －135.］

［10］ Yoshiro H，Takashi O，Takashi K，et al.Synthesis of a linearly interpolated gain scheduling controller for large flexible spacecraft ETS-VIII［J］.Control Engineering Practice，2011(19):611 －625.

［11］ Tasker F，Bosse A， Fisher S. Real-time modal parameter estimation using subspace methods:theory［J］.Mechanical Systems and Signal Processing，1998，12(6):797－808.