一種在軌衛(wèi)星質(zhì)量特性計算方法

張洪波 武向軍 劉天雄 叢飛 周耀華

（北京空間飛行器總體設(shè)計部，北京 100094）

1 引言

衛(wèi)星需要通過高精度軌道和姿態(tài)控制，來實現(xiàn)星體及有效載荷的高精度指向，以完成復(fù)雜的空間任務(wù)。為了實現(xiàn)衛(wèi)星姿態(tài)高精度控制，需要已知星體準確的質(zhì)量特性參數(shù)（包括質(zhì)心、轉(zhuǎn)動慣量和慣性積）。在國內(nèi)外實際航天工程應(yīng)用中，衛(wèi)星質(zhì)量特性參數(shù)多源于地面測量，但地面質(zhì)量特性測量雖耗費大量人力、財力和研制周期，卻無法實現(xiàn)高精度、全生命周期和狀態(tài)的測量［1］，這是由于：第一，重力和測試裝置本身的影響無法完全消除；第二，太陽翼、天線等部件展開狀態(tài)質(zhì)量特性無法測量；第三，衛(wèi)星慣性積量級較小，與一些系統(tǒng)測量誤差量級相當，導(dǎo)致慣性積測量結(jié)果不可用。為此，可借助星上敏感器和執(zhí)行機構(gòu)計算衛(wèi)星在軌質(zhì)量特性，這樣不僅可以實時反映質(zhì)量特性在軌變化情況，控制系統(tǒng)可自動適應(yīng)變化來實現(xiàn)衛(wèi)星高效、高精度的姿態(tài)控制，而且計算分析結(jié)果對于衛(wèi)星總體構(gòu)型布局優(yōu)化、衛(wèi)星表面展開部件設(shè)計優(yōu)化、地面測量系統(tǒng)誤差修正也具有指導(dǎo)意義。

截至目前，國外已有一些學者開展了衛(wèi)星在軌質(zhì)量特性計算方法研究。Bergmann等人提出了一種基于高斯二階濾波的計算方法［2-3］，此算法邏輯復(fù)雜、計算量大，并對動力學模型進行了簡化，只適用于衛(wèi)星轉(zhuǎn)動角速度足夠小的情況；Wilson和Rock等人提出了一種基于指數(shù)加權(quán)遞歸的最小二乘法［4-5］，此算法需要衛(wèi)星配置線加速度計，且計算過程中通過動力學方程的線性近似，人為割裂了質(zhì)量、慣量和質(zhì)心的內(nèi)部耦合關(guān)系，雖然采用“遞歸”方法可以逐步減小或忽略耦合引起的誤差，但在衛(wèi)星角速度較大的情況下誤差也較大；Tanygin 和Williams等人提出了一種使用最小二乘法在自旋衛(wèi)星執(zhí)行在軌機動時辨識其質(zhì)量特性的方法［6］，不適用于三軸穩(wěn)定衛(wèi)星；Kim 等人提出利用遞歸最小二乘法辨識三軸航天器的慣量和質(zhì)心位置［7-8］，存在與Wilson等人算法類似的問題。國內(nèi)在該方面的研究還不多［9］，王書廷等人提出了一種基于遞歸最小二乘法的衛(wèi)星慣量和質(zhì)心位置計算方法［10］，思路和問題與文獻［5］類似；徐文福等人提出了基于POS優(yōu)化的計算方法［9］。對于國內(nèi)外的上述計算方法，有以下4點因素限制了其精度和應(yīng)用范圍：①需要衛(wèi)星配置線加速度計；②割裂參數(shù)內(nèi)部耦合關(guān)系，需要假設(shè)衛(wèi)星旋轉(zhuǎn)足夠慢；③未考慮星體與飛輪的動量交換，并假設(shè)星體采用零動量控制方式來簡化耦合關(guān)系；④一些不進行在軌姿態(tài)大角度機動（連續(xù)轉(zhuǎn)動，且繞三軸均有角速度和角加速度分量）的衛(wèi)星，不適用遞歸算法。

基于上述問題，本文提出一種基于批量最小二乘法的質(zhì)量特性計算方法，僅需要衛(wèi)星配置陀螺，針對偏置動量控制方式建模并考慮星體與飛輪動量交換，通過交互迭代的方式逐步消除耦合誤差，計算基于若干離散噴氣運動狀態(tài)，不需衛(wèi)星連續(xù)噴氣轉(zhuǎn)動。此方法對于衛(wèi)星配置和運動狀態(tài)要求低，充分考慮了參數(shù)耦合關(guān)系，適用范圍更廣。

2 在軌衛(wèi)星姿態(tài)動力學建模

2．1 坐標系定義

衛(wèi)星坐標系定義［11］如下：

（1）星體坐標系ObXbYbZb：原點Ob位于星箭分離面中心，對地模式下ObXb軸指向飛行方向，ObZb軸指向地心，ObYb軸指向根據(jù)右手定則確定。星體坐標系是星上幾何參考基準，推力器安裝方位、衛(wèi)星質(zhì)心等均以此為參考；

（2）質(zhì)心坐標系OmXmYmZm：原點Om位于星體質(zhì)心，三軸指向與星體坐標系相同；

（3）軌道坐標系OoXoYoZo：原點Oo位于星體質(zhì)心，OoZo軸由星體質(zhì)心指向地心，OoXo位于軌道平面內(nèi)，垂直于OoZo軸且沿向飛行方向，OoYo軸指向根據(jù)右手定則確定。

本文所有參量均在星體坐標系ObXbYbZb下定義。

2．2 衛(wèi)星姿態(tài)動力學模型

假設(shè)衛(wèi)星采用偏置動量控制方式，衛(wèi)星配置n個推力器、p個飛輪和q個磁力矩器，衛(wèi)星采用陀螺測量三軸角速度。衛(wèi)星角動量［11］可以表示為

式中：I0為衛(wèi)星轉(zhuǎn)動慣量矩陣（包括飛輪）；ω0為星體角速度；Ii為第i個飛輪的繞自身轉(zhuǎn)軸的轉(zhuǎn)動慣量矩陣；ωi為第i個飛輪相對于星體的角速度。

根據(jù)角動量定理，可得

式中：T為除飛輪控制力矩以外衛(wèi)星受到的其他所有外力矩之和，本文模型中指推力器和磁力矩器引起的力矩，則有

式中：ri為第i個推力器噴口中心點坐標，rm為衛(wèi)星質(zhì)心坐標，F(xiàn)i為第i個推力器的推力矢量，τi為第i個磁力矩器產(chǎn)生的卸載力矩，?0為角速度ω0的叉乘操作數(shù)，對于ω0＝有?0＝

將式（1）代入式（3）可得

即

式（5）中Ii、ω0和˙ω0量級均較小，因此其乘積項可以忽略，由式（5）簡化可得

3 質(zhì)量特性計算方法

式（6）中共含有9個未知數(shù)：Ⅰ0xx、Ⅰ0yy、Ⅰ0zz、Ⅰ0xy、Ⅰ0xz、Ⅰ0yz、rmx、rmy、rmz，耦合關(guān)系強，計算復(fù)雜；文獻［2-3］等假設(shè)星體運動足夠慢，從而忽略?0、I0、ω0項來簡化模型，但同時也降低了計算精度。本文提出一種基于交互式迭代和批量最小二乘法的計算方法，先給定衛(wèi)星轉(zhuǎn)動慣量初值，根據(jù)I00計算r1m，然后根據(jù)計算，再根據(jù)I10計算r2m…… 以此類推，最終算法收斂可以得出I0和rm；每一步計算都使用多元線性批量最小二乘回歸算法，其標準形式［12］為

式中：A為由已知量組成的參數(shù)矩陣，x為待求解未知向量，b為無噪聲測量已知量，ξ為測量噪聲矢量。最小二乘算法的目標是求解，使得-b平方和達到最小。對于式（7）所示的標準形式可以直接使用批量最小二乘算法求解：

本文算法的核心就是將式（6）的姿態(tài)動力學方程轉(zhuǎn)化成式（7）的標準形式，并通過式（8）的模式進行求解。

本文提出的批量最小二乘法與文獻［10］提出的遞歸最小二乘法相比，可以適用于在軌不進行連續(xù)大角度姿態(tài)機動的衛(wèi)星。這類衛(wèi)星只在地球捕獲等工況進行間歇噴氣姿態(tài)調(diào)整，推力器單次噴氣時間約為幾十至幾百毫秒，而噴氣間隔時間可達幾秒、十幾秒甚至幾十秒，如使用遞歸算法，往往會由于相鄰幀數(shù)據(jù)中含有不噴氣的無效數(shù)據(jù)而導(dǎo)致ATA病態(tài)，無法求逆。

3．1 質(zhì)心坐標計算方法

質(zhì)心坐標計算時假設(shè)I0為已知，由式（6）可得

根據(jù)矩陣叉乘交換律有

其中~Fi為Fi的叉乘操作數(shù)，對于

令

則根據(jù)（8）式可以求解質(zhì)心坐標rm。

3．2 轉(zhuǎn)動慣量矩陣計算方法

轉(zhuǎn)動慣量計算時假設(shè)質(zhì)心rm為已知，根據(jù)式（6）有

接下來對I0˙ω0和?0I0ω0兩項進行轉(zhuǎn)化：

綜合式（15）、式（18）和式（19）即可利用批量最小二乘法求出x，即求衛(wèi)星轉(zhuǎn)動慣量矩陣中的6個未知參數(shù)。

4 在軌數(shù)據(jù)分析與驗證

利用某衛(wèi)星在軌實際飛行數(shù)據(jù)對本文算法進行驗證，通過Matlab仿真軟件進行數(shù)值仿真。驗證狀態(tài)選取為星箭分離后衛(wèi)星進行地球捕獲的工況，此工況下衛(wèi)星繞三軸均有角速度分量，姿控推力器雖為間歇噴氣，但會在三軸上均產(chǎn)生角加速度。

衛(wèi)星姿控推力器噴口中心點位置和安裝角度見表1，其中安裝角度指推力器軸線（沿羽流方向）與星體坐標系三軸的夾角；推力指向如圖1所示。衛(wèi)星陀螺角速度和推力器噴氣累計時間見表2。

表1 衛(wèi)星姿控推力器指向及噴口中心點坐標Table 1 Directions and coordinates of attitude control thrusters

圖1 衛(wèi)星姿控推力器推力指向示意圖Fig．1 Directions of attitude control thrusters

表2 陀螺角速度和推力器噴氣累計時間Table 2 Angular velocity of gyro and accumulative working-time of attitude control thrusters

根據(jù)地面質(zhì)量特性測試結(jié)果推算和根據(jù)本文件算法計算的在軌質(zhì)量特性參數(shù)對比見表3。

表3 地面與在軌質(zhì)量特性計算結(jié)果對比Table 3 Results of measurement and calculation

在軌計算結(jié)果與地面計算結(jié)果略有差異，但非常相近，綜合考慮兩種計算過程的誤差（包括地面測量系統(tǒng)誤差、太陽翼展開模型和推進劑微重力模型與真實狀態(tài)差異、衛(wèi)星在軌姿態(tài)機動引起推進劑晃動等），可以認為本文提出的算法正確并具備較高精度。

當初值取為I0＝時，質(zhì)心坐標收斂曲線分別如圖2～圖4所示，轉(zhuǎn)動慣量和慣性積收斂曲線如圖5～圖10所示，本文算法具有較快的收斂特性，在迭代500次左右時就已經(jīng)完全收斂。

此外本文算法對于初值選取有一定的適應(yīng)性。

圖2 質(zhì)心坐標X 收斂曲線Fig．2 Curve of centroid in Xdirection

圖3 質(zhì)心坐標Y 收斂曲線Fig．3 Curve of centroid in Y direction

圖4 質(zhì)心坐標Z 收斂曲線Fig．4 Curve of centroid in Zdirection

圖5 轉(zhuǎn)動慣量Ⅰxx收斂曲線Fig．5 Curve of Ⅰxx

圖6 轉(zhuǎn)動慣量Ⅰyy收斂曲線Fig．6 Curve of Ⅰyy

圖7 轉(zhuǎn)動慣量Ⅰzz收斂曲線Fig．7 Curve of Ⅰzz

圖8 慣性積Ⅰxy收斂曲線Fig．8 Curve of Ⅰxy

圖9 慣性積Ⅰzx收斂曲線Fig．9 Curve of Ⅰzx

圖10 慣性積Ⅰyz收斂曲線Fig．10 Curve of Ⅰyz

圖11 不同初值選取情況下收斂特性對比Fig．11 Curve of different initial value

5 結(jié)束語

準確獲取質(zhì)量特性參數(shù)，是衛(wèi)星實現(xiàn)高精度姿態(tài)控制進而完成復(fù)雜空間任務(wù)的基礎(chǔ)。基于敏感器和執(zhí)行機構(gòu)的在軌質(zhì)量特性計算方法，能夠?qū)崿F(xiàn)全真實狀態(tài)下、全生命周期的質(zhì)量特性實時計算。本文提出的衛(wèi)星在軌質(zhì)量特性計算方法，對衛(wèi)星姿態(tài)傳感器配置要求低，對衛(wèi)星控制模型適應(yīng)性強，對動力學模型參數(shù)耦合關(guān)系考慮充分，對衛(wèi)星在軌運動狀態(tài)要求低，適用范圍更廣泛。通過某衛(wèi)星實際在軌遙測數(shù)據(jù)的計算，結(jié)果表明，此方法計算精度高，收斂快速、穩(wěn)定，可為衛(wèi)星在軌高精度軌道和姿態(tài)控制提供參考。

（References）

［1］王洪鑫，徐在峰，趙科，等．航天器質(zhì)量特性測試技術(shù)新進展［J］．航天器環(huán)境工程，2011，28（2）：171-174

Wang Hongxin，Xu Zaifeng，Zhao Ke，et al．Recent advances of mass property measuring technology for spacecraft［J］．Spacecraft Environment Engineering，2011，28（2）：171-174（in Chinese）

［2］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（2）：483-492

［3］Bergmann E V，Dzielski J．Spacecraft mass property identification with torque-generating control［J］．Journal of Guidance，Control，and Dynamics，1990，13（2）：99-103

［4］ Wilson E，Lages C，Mah R．On-line，gyro-based，mass-property identification for thruster-controlled spacecraft using recursive least squares［C］／／The 45thMidwest Symposium on Circuits and Systems．New York：IEEE，2002

［5］Wilson E，Sutter D W． Motion-based mass and thruster-property identification for thruster-controlled spacecraft，AIAA 2005-6907［C］／／AIAA Infotech Aerospace Conference．Washington D．C．：AIAA，2005

［6］Tanygin S，Williams T．Mass property estimation using coasting maneuvers［J］．Journal of Guidance，Control，and Dynamics，1997，20：625-632

［7］Kim J，Agrawal B N．System identification and automatic mass balancing of ground-based three-axis spacecraft simulator，AIAA 2006-6595［C］／／ AIAA Guidance，Navigation，and Control Conference and Exhibit．Washington D．C．：AIAA，2006

［8］Kim J，Agrawal B N．Automatic mass balancing of airbearing-based three-axis rotational spacecraft simulator［J］．Journal of Guidance，Control，and Dynamics，2009，32（3）：1005-1017

［9］徐文福，何勇，王學謙，等．航天器質(zhì)量特性參數(shù)的在軌辨識方法［J］．宇航學報，2010，31（8）：1906-1914

Xu Wenfu，He Yong，Wang Xueqian，et al．On orbit identification of mass characteristic parameters for spacecraft［J］．Journal of Astronautics，2010，31（8）：1906-1914（in Chinese）

［10］王書廷，曹喜濱．衛(wèi)星質(zhì)量特性的在線辨識算法研究［C］／／第25 屆中國控制會議．北京：中國自動化學會，2006：519-524

Wang Shuting，Cao Xibin．On-line mass-property identification algorithm research for satellite ［C］／／Proceedings of the 25thChinese Control Conference．Beijing：Chinese Association of Automation，2006：519-524（in Chinese）

［11］屠善澄．衛(wèi)星姿態(tài)動力學與控制［M］．北京：中國宇航出版社，1999

Tu Shancheng．Satellite attitude dynamics and control［M］．Beijing：China Astronautics Press，1999 （in Chinese）

［12］張賢達．矩陣分析與應(yīng)用［M］．北京：清華大學出版社，2004

Zhang Xianda．Matrix analysis and application［M］．Beijing：Tsinghua University Press，2004（in Chinese）