一種基于推力器控制的衛星質心在軌估算方法研究

郭正勇，張增安，汪禮成，何益康，趙永德

(上海航天控制技術研究所,上海 201109)

一種基于推力器控制的衛星質心在軌估算方法研究

郭正勇，張增安，汪禮成，何益康，趙永德

(上海航天控制技術研究所,上海 201109)

為準確估計地球靜止軌道(GEO)衛星在軌道機動過程中因燃料消耗產生的質心位置的變化，提出了一種基于推力器連續噴氣的衛星質心在軌估算方法。采用推力器在固定時間內連續噴氣工作方式以形成恒定的推力器作用力和力矩，根據陀螺測量值用最小二乘法估算推力器產生的星體角加速度值，采用產生正負向相反控制力矩的兩個推力器同時工作，以減小對衛星姿態的擾動和三軸間的動力學耦合。給出了衛星質心計算公式，討論了質心估算中的推力器推力位置測量誤差、推力器推力矢量方向角度測量誤差、成對工作推力器推力大小偏差、陀螺組合測量噪聲、整星轉動慣量計算偏差，以及衛星姿態動力學耦合特性等主要誤差源對估算結果的影響。基于某GEO衛星的推力器數據，計算獲得了在軌質心的總測量偏差。仿真結果表明：理論計算值與仿真結果的誤差在允許范圍內，方法有效，可廣泛用于航天器的質心位置測量，方法有較大的工程應用價值。

地球靜止衛星； 衛星質心估算； 推力器； 陀螺； 連續噴氣； 最小二乘法； 誤差源； 測量偏差

0 引言

某地球靜止軌道衛星采用雙組元液體推進系統，4個燃料儲箱采用并聯布局方式，在軌道機動過程中約80%的燃料被消耗。由于儲箱排放不均衡及混合比偏差的影響，衛星的質心位置變化較大，且地面難以實現整個變軌期間的質量特性計算，有可能導致490 N變軌發動機因推力偏心產生的干擾力矩變大，增大姿控推力器消耗的燃料。因此，有必要對衛星的質心位置進行在軌估算，為儲箱的均衡排放控制和混合比調節提供參考[1]。目前，國內外對衛星在軌質量特性計算方法進行了研究，但多數采用對敏感器數據進行連續處理，數據量大、計算復雜，且未考慮敏感器測量誤差、推力器推力誤差等的影響[2-3]。文獻[4-6]用遞推最小二乘法對質心位置進行估算，以推力器作為執行機構，缺點是衛星在軌運行過程中，由于推力器噴氣為隨機的，導致辨識過程時間長；文獻[7]提出了基于角動量裝置的質心位置估算方法，不消耗星上燃料，但要求配置高精度的加速度測量裝置，增加了系統配置；文獻[8]提出了僅根據星敏感器測量信息，用卡爾曼濾波算法對衛星的三軸轉動慣量進行估計，但不能實現對衛星質心位置的辨識；文獻[9]將飛輪、推力器分別作為衛星轉動慣量、質心位置辨識的執行結構，用最小二乘法分步實現對衛星質量特性的辨識，但存在與文獻[3]相同的問題。

推力器產生的力矩與推力器的推力、衛星的質心位置有關，可采用推力器噴氣的方式，根據衛星姿態動力學方程，利用衛星角加速度、角速度信息，結合推力器的安裝可估算出衛星的質心位置[10-12]。因某GEO衛星未安裝角加速度計，需利用陀螺測量值較精確地估算星體角速度值，要求質心測量過程中衛星的角加速度為常值，即推力器產生的作用力及力矩恒定。為此，在不增加額外測量裝置的條件下，本文提出了一種基于推力器固定時間內連續噴氣的方法以實現對衛星質心位置估算，用最小二乘法對陀螺測量數據進行處理，估計衛星角加速度，進而計算衛星的質心位置，討論了不同誤差源的影響，分析了質心估算精度。

1 問題描述

1.1坐標系定義

衛星機械坐標系：為正交坐標系且與星體固定，是衛星安裝布局坐標系，星上單機的安裝布局及質心位置均在該坐標系中描述。坐標原點Oj為平臺承力筒分離面的幾何中心；OjZj軸垂直于星箭分離面，指向衛星本體的對地面板，此軸亦為衛星的縱軸；OjYj軸垂直于指向衛星本體的南面板，與太陽電池陣驅動軸線平行；OjXj軸由右手法則確定。

衛星本體坐標系：與衛星本體固連，由機械坐標系的坐標原點平移至衛星質心后形成。坐標原點Ob為衛星質心；ObXb、ObYb、ObZb軸為衛星的三個幾何軸，與機械坐標系的相應軸平行，分別稱為滾動軸、俯仰軸和偏航軸。

1.2推力器安裝布局

衛星姿控推力器的安裝布局如圖1所示。在標稱條件下，推力器在衛星本體坐標系中產生的力和力矩見表1。

推力器推力(本體系中)/N力矩(本體系中)/(N·m)FxFyFzMxMyMz6A-508660304599208439026557A5086603-04062-209468023458A0586603191304-04178024129A0-586603-1931750448302588

由表1可知：推力器6A/7A、8A/9A主要在俯仰、滾動向產生較大的控制力矩，同時會在其他方向產生附加力矩。

2 衛星質心位置估算

2.1推力器工作方式選擇

由于電磁閥動作時延的影響，推力器存在推力延時特性。在接收到噴氣脈沖時間指令時，推力器實際輸出的動態特性如圖2所示。

由圖1，因推力器自身存在的前、后時延等影響，推力器實際輸出的推力與標稱推力存在偏差，特別是當指令噴氣時間Tn較小時，推力器實際輸出的推力大小與標稱值存在偏差較大，且產生推力的大小變化較大。

為保證質心測量過程中推力器產生的推力近似為恒值，在質心測量過程中采用推力器在固定時間內連續噴氣的工作方式。另外，為防止衛星三軸角速度過大，減小對衛星姿態的擾動及三軸間的動力學耦合作用，選擇產生正、負方向相反控制力矩的兩個推力器同時工作。

2.2陀螺測量數據處理

2.3質心位置計算

衛星在軌運行過程中，環境干擾力矩(重力梯度力矩、地磁力矩、太陽光壓力矩等)遠小于推力器產生的力矩，可認為三軸慣性角速度變化由推力器引起。質心位置測量過程中，選擇產生正、負相反的控制力矩的兩個推力器同時噴氣(表1中推力器6A與7A、8A與9A)，對衛星產生的擾動力矩較小，因此在測量過程中衛星角速度較小，姿態動力學方程可簡化為

(1)

式中：Ixx，Iyy，Izz分別為衛星滾動、俯仰和偏航向主慣量；ωx，ωy，ωz分別為衛星三軸慣性角速度；Tx，Ty，Tz分別為衛星三個方向受到的外部力矩。

推力器6A、7A同時工作，噴氣時間長度為Δt1，陀螺組合測得衛星三軸慣性角速度增量分別為Δωx67,Δωy67,Δωz67，由姿態動力學方程及力矩計算公式可得

(2)

式中：py為衛星y向質心位置；F6a，F7a分別為推力器6A、7A產生的推力幅值；p6ay，α6，γ6分別為推力器6A在機械坐標系中的y向位置坐標、推力的滾動向余弦角和偏航向余弦角；p7ay，α7，γ7分別為推力器7A在機械坐標系中的y向位置坐標、推力滾動向余弦角和推力偏航向余弦角，均可在安裝時測量可得。

根據假設條件，質心測量過程中推力器產生的推力大小為恒定值，對式(2)分子、分母同時求導，則式(2)可變為

(3)

通過推力標定和選擇，成對工作的推力器產生的推力大小偏差≤2%，有F6a/F7a≈1，則式(3)可變為

(4)

(5)

(6)

式中：p8ax，β8，γ8分別為推力器8A在機械坐標系中的x向位置坐標、推力俯仰向余弦角和偏航向余弦角；p9ax，β9，γ9分別為推力器9A在機械坐標系中的x向位置坐標、推力俯仰向余弦角和偏航向余弦角，均可在安裝時測得。

3 測量誤差分析

由式(5)、(6)可知：質心測量的誤差源主要有推力器產生的推力位置測量誤差、推力器產生的推力矢量的方向角度測量誤差、成對工作的推力器的推力大小偏差、陀螺組合測量噪聲、整星轉動慣量計算偏差，以及衛星姿態動力學耦合特性等引起質心測量的偏差。

除整星轉動慣量計算偏差外，其他誤差源主要引起k1，k2值的變化。設由遙測數據算得的k1，k2的值相對標稱值的偏差分別為Δk1，Δk2，衛星質心位置的估算誤差Δp=[ΔpxΔpy]，則在小偏差條件下有

(7)

3.1陀螺測量噪聲導致的誤差

同理可得

可知當推力器工作時間長度越長，陀螺測量噪聲的影響相對較小。以推力器工作時間10 s，陀螺組合隨機測量噪聲δω=0.06 (°)/h(3σ)，可算得陀螺組合測量噪聲引起的x、y向質心測量誤差分別為Δpx=0.002 mm，Δpy=0.003 mm。

3.2推力器推力方向偏差導致的誤差

因推力器本身的性能導致推力器的推力分量與標稱值存在偏差，由k1，k2的表達式可得推力器推力矢量角度偏差引起的k1，k2偏差分別近似為

(8)

(9)

由式(8)、(9)可得推力器的推力變化對k1，k2的影響，由此可確定k1，k2值影響的最壞工況條件。

考慮推力器的推力方向角度偏差0.1°，最壞工況條件下可得推力器推力方向偏差引起的質心估算誤差為Δpx=0.045 mm，Δpy=0.041 mm。

3.3推力器推力位置偏差導致的誤差

因推力器產生的推力實際位置與噴管法蘭中心存在偏差，由k1，k2的表達式可得推力器推力矢量角度偏差引起的k1，k2的偏差分別近似為

(10)

分析可得：推力器6A、7A的y向位置坐標存在同向偏差時引起k1的誤差最大；推力器8A、9A的x向位置坐標存在同向偏差時引起k2的誤差最大。則推力器的推力位置偏差為1.5 mm時，引起的衛星x、y向質心位置坐標的最大估算誤差為1.5 mm。

3.4成對工作推力器推力大小偏差導致的誤差

考慮成對工作的兩個推力器產生的推力大小偏差對質心估算精度的影響。令推力器6A、7A產生的推力大小之比kF67=F6a/F7a，推力器8A、9A產生的推力大小之比kF89=F8a/F9a，則由式(3)可得

(11)

(12)

同時工作的兩推力器產生的推力大小存在偏差，在小偏差條件下近似有

(13)

考慮成對工作的推力器產生的推力大小偏差小于2%，最壞工況條件下由式(7)、(13)可算得Δpx=0.5 mm，Δpy=0.498 mm。

3.5衛星轉動慣量計算偏差導致的誤差

實際工程中，衛星在軌的轉動慣量理論計算值與實際值存在誤差。考慮衛星的實際轉動慣量與理論值間的誤差為10%，考慮最壞工況組合條件

0.18k1non

(14)

0.18k2non

(15)

式中：下標non表示標稱值。

代入式(5)，可得最壞工況條件下衛星質量特性計算偏差引起衛星x向質心最大測量誤差約0.352 mm，y向約0.62 mm。通過在軌對衛星轉動慣量進行估算后，由轉動慣量計算偏差導致的質心測量誤差可忽略。

3.6三軸耦合運動導致的誤差

由于慣量積、陀螺力矩的存在，導致三軸間的耦合會對質心測量精度產生一定的影響。三軸間耦合運動引起的k1，k2的變化可分別表示為

將上述計算值代入式(5)，可算得三軸耦合運動引起的x、y向質心測量誤差分別為Δpx=0.01 mm，Δpy=0.434 mm。

3.7太陽帆板振動的影響

在質心測量過程中，推力器工作會引起太陽帆板撓性振動，導致衛星三軸慣性角速度存在周期性變化，其變化為高頻分量。數據處理過程中通過數據擬合的方式可基本消除帆板振動的影響，因此可忽略帆板振動對質心測量的影響。

3.8總測量誤差

影響測量誤差的不同因素綜合結果見表2。由表2可知：三軸耦合運動導致衛星y向質心位置的測量誤差相對較大，最大測量偏差約3.1 mm。

4 仿真驗證

通過控制推力器6A與7A、8A與9A同時噴氣，將噴氣期間的陀螺測得衛星三軸慣性角速度數據遙測下傳，通過對陀螺測量數據進行處理，可得推力器產生的衛星慣性角加速度，衛星三軸慣性角速度變化如圖3所示。代入式(3)、(4)后可算得衛星質心位置，質心估算結果見表3。

由仿真結果可知：

a)不考慮拉偏時，衛星的質心測量精度優于0.5 mm；考慮成對工作的推力器產生的推力偏差，質心測量誤差約0.4 mm；考慮衛星轉動慣量誤差，質心測量誤差約1.5 mm。

b)數學仿真結果與理論分析結果的誤差在可允許范圍內，數學仿真結果表明理論分析的準確性。

5 結束語

本文提出了一種基于推力器、陀螺組合的衛星質心位置估計的方法。該方法不需配置角加速度測量裝置，通過將陀螺測量數據進行遙測下傳，經由地面數據處理估算出衛星的質心位置。該方法簡化了數據處理過程，避免了以往利用卡爾曼濾波估計等算法時計算量大的缺點；通過采用一對產生方向相反力矩的推力器在固定時間內連續同時工作的方式，避免了推力器的推力大小變化引起的角加速度變化偏差，提高了質心測量精度。該方法可廣泛用于航天器的質心位置測量，為航天器儲箱混合比偏差、不均衡排放的調整提供依據，具有較大的工程應用價值。

該方法在計算過程中用到了衛星的轉動慣量，由于衛星實際的轉動慣量與計算值存在偏差，從而帶來衛星質心位置的測量誤差。在下一步研究中，應考慮對衛星轉動慣量進行在軌估算，減少轉動慣量計算值與實際值之間的偏差，以提高質心測量的精度。

表3 質心測量數學仿真結果

[1] 屠善澄. 衛星姿態動力學與控制(2)[M]. 北京: 中國宇航出版社, 2009.

[2] 李春華. 航天器動態質量特性分析[J]. 航天器工程, 2005, 14(3): 15-20.

[3] 張洪波, 武向軍, 劉天雄, 等. 一種在軌衛星質量特性計算方法[J]. 航天器工程, 2013, 22(6): 30-36.

[4] 王書廷, 曹喜濱. 衛星質量特性的在線辨識算法研究[C]∥第25屆中國控制會議. 哈爾濱： [s. n.]， 2006: 524-528.

[5] WILSON E, SUTTER D W. Motion-based mass and thruster-property identification for thruster-controlled spacecraft[R]. AIAA, 2005-6907, 2005.

[6] WILSON E, LAGES C, MAH R. On-line, gyro-based, mass-property identification for thruster-controlled spacecraft using recursive least squares[C]∥The 45th Midwest Symposium on Circuits and Systems. Tuls: [s. n.], 2002: 334-337.

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

[8] SU J J, BBAN H. Mass Property estimation of gyroless spacecraft[C]∥ 11thInternational Conference on Control, Automation and Systems. KINTEX, Gyeonggi-do： [s. n.], 2011: 1884-1887.

[9] 徐文福, 何勇, 王學謙, 等. 航天器質量特性參數的在軌辨識方法[J]. 宇航學報, 2010, 31(8): 1905-1914.

[10] 劉偉霞, 熊智, 郁豐, 等. 組合航天器轉動慣量在軌兩步辨識標定[J]. 中國空間科學技術, 2013, 33(2): 32-39.

[11] 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: 483-491.

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

StudyofOn-OrbitEstimationMethodofSatellite’CentroidBasedonThrustControl

GUO Zheng-yong， ZHANG Zeng-an， WANG Li-cheng， HE Yi-kang， ZHAO Yong-de

(Shanghai Institute of Spaceflight Control Technology, Shanghai 201109, China)

To identify the centroid change of on-orbit geostationary satellite caused by propellant consumption, an on-orbit centroid estimation method with thruster’s continuous jet was put forward in this paper. Constant thrust and torque of the thruster were performed through continuous jet during fixed time period. Angular acceleration of satellite produced by thruster was estimated by the least squares-based algorithms according to the gyro measurements. The two thrusters which produced control torques with opposite directions were operating at the same time, which would reduce disturbance of satellite attitude and dynamic coupling among three axes. The computation equations for centroid were given. The effects of main factors on centroid estimation error were discussed, which were thrust position measurement error of thruster, angular measurement error of thrust vector of thruster, thrust value difference between the two thrusters operated, measurement noise of gyros unit, calculation error of rotational inertia of the whole satellite, and dynamic coupled characteristics of satellite attitude. Based on the trust of some GEO satellite, the total on-orbit estimation error of centroid was obtained. The simulation results showed that the difference between the computation and simulation was allowed. The method proposed is effective and has its value in engineering application, which can be used in the measurement of spaceraft centrold.

geostationary orbit satellite; satellite centroid estimation; thruster; gyro; continuous jet; the least squares-based algorithms; error source; estimation error

1006-1630(2017)05-0076-07

2016-12-28；

2017-03-14

國家重大航天工程

郭正勇(1986—)，男，工程師，主要從事衛星姿態動力學建模及高精度姿態控制等研究。

V412.4

A

10.19328/j.cnki.1006-1630.2017.05.012