大角度機動下帶撓性附件航天器轉動慣量在軌辨識

何 驍，譚述君，吳志剛,2

(1.大連理工大學航空航天學院，大連 116024；2.大連理工大學工業裝備結構分析國家重點實驗室，大連 116024)

大角度機動下帶撓性附件航天器轉動慣量在軌辨識

何 驍1，譚述君1，吳志剛1,2

(1.大連理工大學航空航天學院，大連 116024；2.大連理工大學工業裝備結構分析國家重點實驗室，大連 116024)

針對大角度機動情況下帶撓性附件航天器轉動慣量在軌辨識的問題，提出一種將轉動慣量參數估計和撓性附件狀態估計相結合的并發遞推算法。該算法以大角度機動情況下帶撓性附件航天器的非線性動力學模型為基礎，分別利用廣義卡爾曼濾波做撓性附件振動模態的狀態估計，最小二乘法做轉動慣量的參數估計。最后通過并發遞推算法將二者結合，完成了帶撓性附件航天器的轉動慣量參數辨識。為了提高算法的效率，采用一步最小二乘、多步廣義卡爾曼濾波并發遞推的算法。仿真結果表明,該辨識方法兼具高精度、高效率，并且算法有一定的抗干擾能力。

撓性附件；轉動慣量；在軌參數辨識；廣義卡爾曼濾波；最小二乘

0 引 言

隨著航天技術的發展，許多大型航天器需要在軌展開和改變形狀，地面測量的質量特性參數很難真實地反映航天器在軌運行時的情況，因此需要在軌辨識。另外，航天材料向著輕質化、柔性化方向發展，航天器撓性附件振動對其質量特性參數辨識的影響不可忽略[1]。由于航天器的任務越來越復雜，存在大角度機動的工況，動力學模型需要考慮非線性項的影響。因此，研究大角度機動下帶撓性附件航天器轉動慣量在軌辨識具有很重要的價值。

目前，國內外學者針對航天器轉動慣量辨識問題的研究工作主要集中在剛體航天器方面。Bergman等[2]和Bergman等[3]提出一種高斯二階濾波辨識衛星質量特性參數的方法，通過在常規卡爾曼濾波算法的基礎上加入泰勒展開多項式二階項，以提高濾波精度。Tanygin等[4]采用最小二乘算法來辨識自旋穩定衛星的質量特性參數，包括轉動慣量和質心位置等參數。Wilson等[5-6]以SPHERE試驗衛星為背景提出一種多變量并發遞推最小二乘法在線辨識衛星的轉動慣量、質心位置、總質量等質量特性參數。在國內，王書廷等[7]將質心位置和慣量矩陣的辨識問題解耦為兩個最小二乘問題，提出一種遞推最小二乘法在線辨識衛星的轉動慣量和質心位置等參數。黃河等[8]針對小衛星轉動慣量辨識提出了一種基于最小二乘的閉環辨識方法，在航天器完成姿態機動任務的同時，能夠快速辨識出航天器的轉動慣量。許瑩等[9]研究帶太陽電池陣的衛星慣量辨識方法，在建立帶約束條件的優化辨識模型的基礎上，基于約束最小二乘算法精確求解本體慣量和太陽電池陣慣量在內的12個變量值。徐文福等[10]提出了基于參數解耦的最小二乘法和基于PSO的非線性優化兩種方法，辨識航天器的轉動慣量、質量和質心位置。然而這些研究成果均將航天器視為剛體，沒有考慮航天器撓性附件振動對航天器質量特性參數尤其是轉動慣量在軌辨識的影響。近年來已有學者開展帶撓性附件航天器轉動慣量在軌辨識方面的工作。蘭聰超等[11]針對小角度機動情況下帶撓性附件衛星質量特性參數在軌辨識方法的問題，利用最小二乘與卡爾曼濾波并發遞推的算法，辨識了帶撓性附件衛星的轉動慣量。朱東方等[12]針對復雜航天器的質量特性辨識問題，提出了基于擴展卡爾曼濾波的一種復雜航天器轉動慣量矩陣的辨識方法。但現有的針對帶撓性附件航天器的辨識算法對大角度機動工況下轉動慣量參數辨識精度較差，且收斂速度很慢，需進一步研究高精度、高效率的辨識方法。

針對以上問題以及工程發展的需求，本文以典型的航天器——帶撓性附件的衛星作為研究對象，控制力矩為輸入信號，陀螺儀測量的姿態角速度等為輸出信號，針對大角度機動情況，結合最小二乘法和廣義卡爾曼濾波算法，提出了一種可在軌辨識帶撓性附件航天器轉動慣量參數的并發遞推算法。

1 動力學模型

當衛星姿態角變化時，帶撓性附件衛星姿態動力學方程和撓性附件振動方程[13-15]為

(1)

(2)

2 帶撓性附件衛星轉動慣量參數在軌辨識的并發遞推算法

2.1轉動慣量參數的估計

將待辨識的轉動慣量參數表示成標稱值和殘差值相加的形式，

Jsat=Jnom+ΔJ

(3)

將式(3)代入式(1)，得到

(4)

式(4)等號左邊的處理為

(5)

(6)

式中：

基于式(5)和式(6)的描述形式，式(4)可以表示為

AJxJ=bJ

(7)

式中:

2.2撓性附件振動模態的估計

在帶撓性附件衛星轉動慣量的最小二乘描述形式中，當撓性附件振動模態已知時才可利用最小二乘法對衛星轉動慣量參數進行辨識。目前，針對航天器撓性模態參數辨識的工作很多，經典的如文獻[16-17]。對于撓性附件振動模態，本文將利用廣義卡爾曼濾波算法來估計辨識。

2.2.1 廣義卡爾曼濾波

非線性離散系統模型為

x(k+1)=f(x(k),k)+Γ(k)w(k)

(8)

y(k)=h(x(k),k)+v(k)

(9)

其中,測量噪聲v(k)和過程噪聲w(k)為互不相關的零均值白噪聲。測量噪聲和過程噪聲的方差矩陣分別為R和Q，且滿足

(10)

對上述這類非線性問題，可以通過局部線性化將卡爾曼濾波推廣到非線性系統[18]，得到廣義卡爾曼濾波的遞推方程如下，

(11)

P(k+ 1|k) =Φ(k+ 1|k)P(k)ΦT(k+ 1 |k) +Q

(12)

K(k+1)=P(k+1|k)HT(k+1)·

[H(k+1)P(k+1|k)HT(k+1)+R]-1

(13)

(14)

P(k+ 1) = [I-K(k+ 1)H(k+ 1)]P(k+ 1|k)

(15)

其中,狀態轉移矩陣Φ(k+1|k)和觀測矩陣H(k+1)分別是f和h的雅可比矩陣。濾波初值和濾波方差矩陣的初值分別為x(0)=E[x(0)],P(0)=var[x(0)]。

2.2.2 系統狀態的廣義卡爾曼濾波描述形式

將帶撓性附件衛星動力學方程式(1)和(2)轉化為一階微分方程的形式

(16)

式中:

設Aa(Jsat,xa)=D-1Axa,Ba(Jsat)=D-1B，式(16)可以表示為

(17)

引入觀測方程

ya=Caxa

(18)

其中,觀測值為姿態角和姿態角速度，即

采用差分近似將方程(17)離散化，差分的時間間隔為Ts

(19)

整理，得

x(k+1)= (TsAa(Jsat(k),x(k))+x(k))+

TsBa(Jsat(k))u(k)

(20)

同樣，觀測方程離散為

y(k)=Cax(k)

(21)

將式(20)和(21)與前面非線性離散系統模型的一般描述式(8)、(9)對比，有

f(x(k),k)= (TsAa(Jsat(k),x(k))+x(k))+

TsBa(Jsat(k))u(k)

(22)

h(x(k),k)=Cax(k)

(23)

將式(20)和式(21)代入廣義卡爾曼濾波遞推式(11)～(15)可實現狀態量的估計，其中線性化矩陣Φ(k+1|k)和H(k+1)如下所示

(24)

其中,

(25)

H(k+1)=Ca

(26)

2.3并發遞推算法

因為采用差分近似離散，對系統模型存在一定的近似，為了提高精度，廣義卡爾曼濾波估計的采樣周期應該盡量短。而最小二乘法做參數估計沒有近似，較長的采樣周期也能得到很好的辨識結果。綜合以上兩點，在保證算法精度的前提下，為了提高效率，本文采用q步廣義卡爾曼濾波與1步最小二乘并發遞推。其算法的具體流程如圖1所示。先用廣義卡爾曼濾波進行振動模態狀態估計q次，體現為圖1中內環；再用最小二乘法進行轉動慣量參數的辨識，體現為圖1中外環。反復迭代遞推，即可得到轉動慣量的辨識值。

圖1 辨識算法的流程圖Fig.1 Flow chart of the identification algorithm

3 仿真校驗

3.1算例模型

仿真算例選取某型號通信衛星模型，如圖2所示。衛星主要的撓性附件為對稱的兩個太陽能帆板，每個太陽能帆板長8.1 m，質量36.6 kg，衛星展開總跨度18.4 m，總質量2850.8 kg。該衛星是典型的帶有撓性附件的航天器。該衛星模型的動力學分析表明，附件振動對整星動力學特性的影響主要由左右兩側帆板的第一階振動模態決定，因此本文算例中只考慮第一階模態的影響，忽略高階模態。左右兩側帆板的一階模態頻率均為1.2754 Hz，模態阻尼比為0.005，轉動慣量真實值Jreal和帆板轉動剛柔耦合系數矩陣Prot為

圖2 含有兩個太陽帆板的帶撓性附件衛星Fig.2 A satellite with two flexible solar panels

在數值仿真中，姿態角、角速度的初始值為φ0=[0,0,0]T，ω0=[0,0,0]T。仿真的輸入力矩取為方波力矩信號，如圖3所示，仿真時長為100 s。在實際工程中，考慮對姿態角速度的測量誤差一般隨著角速度的增大而增大，因此在該數值仿真中，姿態角速度響應在2°/s以內，圖4為姿態角和姿態角速度響應曲線。

圖3 輸入力矩Fig.3 Input torque

3.2算法有效性校驗

為了校驗本文提出的廣義卡爾曼濾波與最小二乘并發遞推算法(EKL算法)的有效性，仿真算例還給出了不考慮非線性項的卡爾曼濾波與最小二乘并發遞推算法(KL算法)的辨識結果作為對比。

采用本文提出的q步廣義卡爾曼濾波與1步最小二乘相結合的并發遞推算法，其中q等于100，采樣間隔Ts為0.001 s。由于并發遞推算法中用到的姿態角加速度不能通過直接測量得到，本文在算法中利用測量的姿態角速度進行差分近似計算得到。遞推算法的初始值Jnom取為

圖4 力矩為M時的姿態角與角速度曲線Fig.4 The curve of attitude angle and angular velocity when the torque is M

當輸入力矩較小時(0.1M)，即衛星做小角度機動時，分別采用KL與EKL辨識算法進行轉動慣量辨識，得到的辨識結果如圖5所示。可以看出，兩種算法都得到了正確的辨識結果，這是因為當衛星做小角度機動時，衛星的角速度非常小，可做近似忽略，這也驗證了小角度機動時帶撓性附件衛星轉動慣量辨識算法忽略非線性項的合理性。當輸入力矩為M時，采用KL與EKL辨識算法得到的轉動慣量辨識結果如圖6所示。可以看出，當衛星做大角度機動時，KL算法對主慣量辨識結果較好，但對慣性積的辨識結果則很差，甚至得不到收斂的結果。而本文提出的EKL遞推算法可辨識出轉動慣量參數，且精度較高。

圖5 輸入力矩為0.1M時的辨識結果Fig.5 Identification results when the input torque is 0.1M

表1給出了輸入力矩為M時，分別采用KL和EKL算法得到的轉動慣量辨識結果。可以看出，無論是什么情況，本文提出的EKL算法對主慣量辨識結果的相對誤差在1%以內，慣性積相對誤差在10%以內。綜合上述對比分析表明，無論衛星做小角度機動，還是大角度機動，本文的EKL并發遞推算法都能高效、準確地獲得辨識結果，而KL算法由于忽略姿態角速度的非線性項影響，在大角度機動時已不能得到正確的辨識結果。因此，當衛星做大角度機動時，應該考慮非線性項的影響，而本文提出的多步廣義卡爾曼濾波和最小二乘法相結合的并發遞推算法則是處理帶撓性附件衛星做大角度機動轉動慣量參數辨識的有效方法。

圖6 輸入力矩為M的辨識結果Fig.6 Identification results when the input torque is M

辨識參數JxJyJzJxyJxzJyz真實值3035.43691800.28923934.274449.017-23.46-27.892KL方法辨識結果2993.981780.893740.0448.80-160.91-49.97相對誤差1.37%1.08%4.94%0.43%585.89%79.17%EKL方法辨識結果3030.621784.453929.2945.46-21.50-30.24相對誤差0.16%0.88%0.13%7.26%8.37%8.43%

3.3算法抗干擾性校驗

上面的仿真算例中，由于只是對算法的有效性進行校驗，并沒有考慮噪聲干擾對算法的影響，而在實際情況中不可能沒有干擾，下面對并發遞推算法的抗干擾能力進行研究。在其他條件都相同的情況下，在測量的姿態角、角速度的信號中加入信噪比(Signal to noise ratio, SNR)為140 dB的高斯白噪聲。信噪比是指系統中信號與噪聲的比例，計量單位為dB。在實際工程中，輸入力矩的測量也會存在一定的誤差，在仿真中考慮10%的誤差，所以在輸入力矩中加入信噪比30 dB的噪聲。

圖7給出了考慮噪聲干擾情況下當輸入力矩為M時的轉動慣量辨識結果。與圖6的結果相比，由于考慮了噪聲干擾，辨識值收斂速度較慢，精度有所降低，但還是能獲得收斂的結果。這說明本文提出的并發遞推算法具有一定的抗干擾能力，而當測量信號的噪聲繼續增加時，則會進一步影響到辨識結果的精度與收斂速度。這就對辨識過程中用到的測量儀器的精度提出了一定的要求。

圖7 考慮噪聲的結果圖Fig.7 Identification results in consideration of noise

4 結 論

本文針對大角度機動情況下帶撓性附件航天器轉動慣量參數辨識的問題，以帶撓性附件衛星的動力學模型為基礎導出了轉動慣量參數的最小二乘辨識和撓性附件振動模態狀態估計相結合的并發遞推算法。本文考慮了姿態方程中的非線性項，得到了解決帶撓性附件航天器轉動慣量參數辨識問題更通用的方法。為了提高遞推算法的精度和效率，本文采用了多步廣義卡爾曼濾波與一步最小二乘法并發遞推的算法。該算法只需要借助衛星原有的測量裝置和執行機構，不需要額外的設備，工程可行性強。仿真算例的結果表明：在大角度機動的情況下，忽略動力學方程中非線性項的影響會減低辨識精度、甚至得到錯誤的辨識結果。而本文提出的多步廣義卡爾曼濾波與一步最小二乘法的并發遞推算法則可以實現帶撓性附件航天器轉動慣量參數的高精度、高效率辨識。

[1] 倪智宇, 鄔樹楠, 吳志剛, 等. 利用改進 TW-API 方法在軌辨識撓性航天器時變模態參數[J]. 宇航學報, 2015, 36(7): 769-776. [Ni Zhi-yu, Wu Shu-nan, Wu Zhi-gang, et al. On-orbit identification of time-varying modal parameters of flexible spacecraft by an improved TW-API method [J]. Journal of Astronautics, 2015, 36(7): 769-776.]

[2] Bergmann E, Dzielski J. Spacecraft mass property identification with torque-generating control[J]. Journal of Guidance, Control, and Dynamics, 1990, 13(1): 99-103.

[3] Bergmann E V, Walker B K, Levy D R. Mass property estimation for control of asymmetrical satellites[J]. Journal of Guidance, Control, and Dynamics, 1987, 10(5): 483-491.

[4] Tanygin S, Williams T. Mass property estimation using coasting maneuvers[J]. Journal of Guidance, Control, and Dynamics, 1997, 20(4): 625-632.

[5] Wilson E, Sutter D W, Mah R W. MCRLS for on-line spacecraft mass-and thruster-property identification[J]. Parameters, 2004, 446(155): 324.

[6] Wilson E, Sutter D W, Mah R W. Multiple concurrent recursive least squares identification [C]. Proc IASTED International Conference on Intelligent Systems and Control, Honolulu,USA,August 23-25,2004.

[7] 王書廷, 曹喜濱. 衛星質量特性的在線辨識算法研究[C]. 第25屆中國控制會議論文集, 哈爾濱, 中國, 2006年8月7-11日. [Wang Shu-ting, Cao Xi-bin. On-line mass-property identifica-tion algorithm research for satellite [C]. The 25th Chinese Control Conference, Harbin, China, August 7-11, 2006.]

[8] 黃河, 周軍, 劉瑩瑩. 航天器轉動慣量在線辨識[J]. 系統仿真學報, 2010,22(5): 1117-1120. [Huang He, Zhou Jun, Liu Ying-ying. On-line identification of spacecraft moment of inertia [J]. Journal of System Simulation, 2010, 22(5): 1117-1120.]

[9] 許瑩, 呂旺, 李云端, 等. 多剛體衛星轉動慣量在軌辨識[J]. 空間控制技術與應用, 2015, 41(6): 31-36. [Xu Ying, Lv Wang, Li Yun-duan, et al. Identification of rotary inertia for multi-body satellite [J]. Aerospace Control and Application, 2015, 41(6): 31-36.]

[10] 徐文福, 何勇, 王學謙, 等. 航天器質量特性參數的在軌辨識方法[J]. 宇航學報, 2010,31(8): 1906-1914. [Xu Wen-fu, He yong, Wang Xue-qian, et al. On orbit identification of mass characteristic parameters for spacecraft [J]. Journal of Astronautics, 2010, 31(8): 1906-1914.]

[11] 蘭聰超, 譚述君, 吳志剛, 等. 帶撓性附件衛星轉動慣量的在軌辨識[J]. 振動與沖擊, 2017, 36(8): 16-21.[Lan Cong-chao, Tan Shu-jun, Wu Zhi-gang, et al. On-orbit identification of moments of inertia for satellites with flexible appendages [J]. Journal of Vibration and Shock, 2017, 36(8): 16-21.]

[12] 朱東方, 王衛華, 宋婷, 等. 復雜撓性航天器轉動慣量在線辨識算法研究[J]. 上海航天, 2015, 32(5): 1-8, 14.[Zhu Dong-fang, Wang Wei-hua, Song Ting, et al. On-line identification of flexible spacecraft moment of inertia[J]. Aeros-pace Shanghai, 2015, 32(5):1-8,14.]

[13] 游斌弟, 溫建民, 張廣玉, 等. 航天器薄殼柔性附件展開耦合行為特性研究[J]. 宇航學報, 2015, 36(6): 640-647. [You Bin-di, Wen Jian-min, Zhang Guang-yu, et al. Study on coupling behavior of spacecraft deployment with flexible appendages of shell structure [J]. Journal of Astronautics, 2015, 36(6): 640-647.]

[14] 呂旺, 向明江, 葉文郁, 等. 撓性衛星在軌非約束模態計算研究[J]. 宇航學報, 2014, 35(4): 404-409. [Lv Wang, Xiang Ming-jiang, Ye Wen-yu, et al. Research on calculation of on-orbit unconstrained modal of flexible satellite [J]. Journal of Astronautics, 2014, 35(4): 404-409.]

[15] 張嘉鐘, 魏英杰, 曹偉. 飛行器動力學與控制[M]. 哈爾濱:哈爾濱工業大學出版社, 2011.

[16] 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.

[17] Juang J N. An eigensystem realization algorithm using data correlations (ERA/DC) for modal parameter identification[J]. Journal of Control Theory and Advanced Technology, 1988, 4(1): 5.

[18] 蔡金獅. 飛行器氣動參數辨識學[M]. 北京:國防工業出版社, 2003.

On-OrbitIdentificationoftheMomentofInertiaforaSpacecraftwithFlexibleAppendagesduringaLarge-AngleManeuver

HE Xiao1, TAN Shu-jun1, WU Zhi-gang1,2

(1. School of Aeronautics and Astronautics, Dalian University of Technology, Dalian 116024，China;2. State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian 116024，China)

For the problem of the on-orbit identification of the moment of inertia for a spacecraft with flexible appendages during a large-angle maneuver, a concurrent recursive algorithm is proposed in this paper, which combines the estimation of the inertia parameters with the state estimation of the flexible appendages’ vibration modes. The algorithm is based on the nonlinear dynamic model of a spacecraft with flexible appendages during a large-angle maneuver. The extended Kalman filter is used to estimate the states of the flexible appendages’ vibration modes, and the least square method is used to estimate the parameters of the moment of inertia. Finally, the inertia parameters of the spacecraft are obtained by the concurrent recursive algorithm. In order to improve the efficiency of the algorithm, a recursive algorithm based on one-step least square method and multi-step extended Kalman filter is proposed. Simulation results show that the proposed algorithm not only improves the computational accuracy and efficiency, but also has the capability of anti-interference.

Flexible appendages；Moment of inertia；On-orbit parameter identification；Extended Kalman filter (EKF)；Least square

V557+.1

A

1000 -1328(2017)09- 0927- 09

10.3873/j.issn.1000-1328.2017.09.005

2017- 01-19;

2017- 06-22

國家自然科學基金(11572069，11502040，11432010)；中央高校基本科研業務費專項資金(DUT16ZD225)

何驍(1992-)，男，碩士生，主要從事飛行器動力學系統辨識研究。

通信地址：遼寧省大連理工大學綜合1號實驗樓406(116023)

電話：(0411)84706213

E-mail: jsychx@mail.dlut.edu.cn

譚述君(1979-)，男，博士，副教授，主要從事航天器動力學與控制、在軌參數辨識、保結構數值方法、魯棒與最優控制研究。本文通信作者。

通信地址：遼寧省大連理工大學綜合1號實驗樓409A(116023)

電話：(0411)84706213

E-mail: tansj@dlut.edu.cn