一種基于推力器控制的衛(wèi)星質(zhì)心在軌估算方法研究

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

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

一種基于推力器控制的衛(wèi)星質(zhì)心在軌估算方法研究

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

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

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

地球靜止衛(wèi)星； 衛(wèi)星質(zhì)心估算； 推力器； 陀螺； 連續(xù)噴氣； 最小二乘法； 誤差源； 測量偏差

0 引言

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

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

1 問題描述

1.1坐標(biāo)系定義

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

衛(wèi)星本體坐標(biāo)系：與衛(wèi)星本體固連，由機(jī)械坐標(biāo)系的坐標(biāo)原點(diǎn)平移至衛(wèi)星質(zhì)心后形成。坐標(biāo)原點(diǎn)Ob為衛(wèi)星質(zhì)心；ObXb、ObYb、ObZb軸為衛(wèi)星的三個幾何軸，與機(jī)械坐標(biāo)系的相應(yīng)軸平行，分別稱為滾動軸、俯仰軸和偏航軸。

1.2推力器安裝布局

衛(wèi)星姿控推力器的安裝布局如圖1所示。在標(biāo)稱條件下，推力器在衛(wèi)星本體坐標(biāo)系中產(chǎn)生的力和力矩見表1。

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

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

2 衛(wèi)星質(zhì)心位置估算

2.1推力器工作方式選擇

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

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

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

2.2陀螺測量數(shù)據(jù)處理

2.3質(zhì)心位置計算

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

(1)

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

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

(2)

式中：py為衛(wèi)星y向質(zhì)心位置；F6a，F(xiàn)7a分別為推力器6A、7A產(chǎn)生的推力幅值；p6ay，α6，γ6分別為推力器6A在機(jī)械坐標(biāo)系中的y向位置坐標(biāo)、推力的滾動向余弦角和偏航向余弦角；p7ay，α7，γ7分別為推力器7A在機(jī)械坐標(biāo)系中的y向位置坐標(biāo)、推力滾動向余弦角和推力偏航向余弦角，均可在安裝時測量可得。

根據(jù)假設(shè)條件，質(zhì)心測量過程中推力器產(chǎn)生的推力大小為恒定值，對式(2)分子、分母同時求導(dǎo)，則式(2)可變?yōu)?/p>

(3)

通過推力標(biāo)定和選擇，成對工作的推力器產(chǎn)生的推力大小偏差≤2%，有F6a/F7a≈1，則式(3)可變?yōu)?/p>

(4)

(5)

(6)

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

3 測量誤差分析

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

除整星轉(zhuǎn)動慣量計算偏差外，其他誤差源主要引起k1，k2值的變化。設(shè)由遙測數(shù)據(jù)算得的k1，k2的值相對標(biāo)稱值的偏差分別為Δk1，Δk2，衛(wèi)星質(zhì)心位置的估算誤差Δp=[ΔpxΔpy]，則在小偏差條件下有

(7)

3.1陀螺測量噪聲導(dǎo)致的誤差

同理可得

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

3.2推力器推力方向偏差導(dǎo)致的誤差

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

(8)

(9)

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

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

3.3推力器推力位置偏差導(dǎo)致的誤差

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

(10)

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

3.4成對工作推力器推力大小偏差導(dǎo)致的誤差

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

(11)

(12)

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

(13)

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

3.5衛(wèi)星轉(zhuǎn)動慣量計算偏差導(dǎo)致的誤差

實(shí)際工程中，衛(wèi)星在軌的轉(zhuǎn)動慣量理論計算值與實(shí)際值存在誤差?？紤]衛(wèi)星的實(shí)際轉(zhuǎn)動慣量與理論值間的誤差為10%，考慮最壞工況組合條件

0.18k1non

(14)

0.18k2non

(15)

式中：下標(biāo)non表示標(biāo)稱值。

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

3.6三軸耦合運(yùn)動導(dǎo)致的誤差

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

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

3.7太陽帆板振動的影響

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

3.8總測量誤差

影響測量誤差的不同因素綜合結(jié)果見表2。由表2可知：三軸耦合運(yùn)動導(dǎo)致衛(wèi)星y向質(zhì)心位置的測量誤差相對較大，最大測量偏差約3.1 mm。

4 仿真驗證

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

由仿真結(jié)果可知：

a)不考慮拉偏時，衛(wèi)星的質(zhì)心測量精度優(yōu)于0.5 mm；考慮成對工作的推力器產(chǎn)生的推力偏差，質(zhì)心測量誤差約0.4 mm；考慮衛(wèi)星轉(zhuǎn)動慣量誤差，質(zhì)心測量誤差約1.5 mm。

b)數(shù)學(xué)仿真結(jié)果與理論分析結(jié)果的誤差在可允許范圍內(nèi)，數(shù)學(xué)仿真結(jié)果表明理論分析的準(zhǔn)確性。

5 結(jié)束語

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

該方法在計算過程中用到了衛(wèi)星的轉(zhuǎn)動慣量，由于衛(wèi)星實(shí)際的轉(zhuǎn)動慣量與計算值存在偏差，從而帶來衛(wèi)星質(zhì)心位置的測量誤差。在下一步研究中，應(yīng)考慮對衛(wèi)星轉(zhuǎn)動慣量進(jìn)行在軌估算，減少轉(zhuǎn)動慣量計算值與實(shí)際值之間的偏差，以提高質(zhì)心測量的精度。

表3 質(zhì)心測量數(shù)學(xué)仿真結(jié)果

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

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

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

[4] 王書廷, 曹喜濱. 衛(wèi)星質(zhì)量特性的在線辨識算法研究[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] 徐文福, 何勇, 王學(xué)謙, 等. 航天器質(zhì)量特性參數(shù)的在軌辨識方法[J]. 宇航學(xué)報, 2010, 31(8): 1905-1914.

[10] 劉偉霞, 熊智, 郁豐, 等. 組合航天器轉(zhuǎn)動慣量在軌兩步辨識標(biāo)定[J]. 中國空間科學(xué)技術(shù), 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—)，男，工程師，主要從事衛(wèi)星姿態(tài)動力學(xué)建模及高精度姿態(tài)控制等研究。

V412.4

A

10.19328/j.cnki.1006-1630.2017.05.012