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

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

1 引言

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

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

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

2 在軌衛(wèi)星姿態(tài)動(dòng)力學(xué)建模

2．1 坐標(biāo)系定義

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

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

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

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

本文所有參量均在星體坐標(biāo)系ObXbYbZb下定義。

2．2 衛(wèi)星姿態(tài)動(dòng)力學(xué)模型

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

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

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

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

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

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

即

式（5）中Ii、ω0和˙ω0量級(jí)均較小，因此其乘積項(xiàng)可以忽略，由式（5）簡(jiǎn)化可得

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

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

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

本文算法的核心就是將式（6）的姿態(tài)動(dòng)力學(xué)方程轉(zhuǎn)化成式（7）的標(biāo)準(zhǔn)形式，并通過(guò)式（8）的模式進(jìn)行求解。

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

3．1 質(zhì)心坐標(biāo)計(jì)算方法

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

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

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

令

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

3．2 轉(zhuǎn)動(dòng)慣量矩陣計(jì)算方法

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

接下來(lái)對(duì)I0˙ω0和?0I0ω0兩項(xiàng)進(jìn)行轉(zhuǎn)化：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖7 轉(zhuǎn)動(dòng)慣量Ⅰ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 不同初值選取情況下收斂特性對(duì)比Fig．11 Curve of different initial value

5 結(jié)束語(yǔ)

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

（References）

［1］王洪鑫，徐在峰，趙科，等．航天器質(zhì)量特性測(cè)試技術(shù)新進(jìn)展［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］徐文福，何勇，王學(xué)謙，等．航天器質(zhì)量特性參數(shù)的在軌辨識(shí)方法［J］．宇航學(xué)報(bào)，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ì)量特性的在線辨識(shí)算法研究［C］／／第25 屆中國(guó)控制會(huì)議．北京：中國(guó)自動(dòng)化學(xué)會(huì)，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)動(dòng)力學(xué)與控制［M］．北京：中國(guó)宇航出版社，1999

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

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

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