信思博，趙訓友，鄭藝裕，李綠萍

上海衛(wèi)星工程研究所，上海 200240

1 引言

中國火星探測器天問一號于2020年7月23日發(fā)射，隨后進入地火轉(zhuǎn)移軌道。不同于常規(guī)地球高軌衛(wèi)星，天問一號探測器僅需不到一周時間完成變軌，在接下來的任務飛行中，需經(jīng)歷中途修正、近火制動捕獲、環(huán)火軌道調(diào)整和保持、器器分離等變軌流程[1-2]，涉及推進系統(tǒng)多次工作，其狀態(tài)的穩(wěn)定與健康與否直接決定任務的成敗。此外，地面需根據(jù)推力大小、方向偏差、燃料比沖等參數(shù)制定軌控策略，其實際工作參數(shù)的穩(wěn)定直接影響策略的制定和修正，因此通過在軌參數(shù)辨識，分析推力大小、推力方向、整器質(zhì)心等參數(shù)尤為必要[3-6]。例如日本的隼鳥號小行星探測器利用飛輪轉(zhuǎn)速變化識別出霍爾推力器故障，采取了相應的補救措施[7]，避免了航天器完全失效。

通常在軌推進系統(tǒng)參數(shù)標定主要依靠測定軌結(jié)果反算推力器產(chǎn)生的沖量。該方法主要針對主發(fā)動機的參數(shù)標定，可以根據(jù)點火時長遙測標定出推力大小和推力方向。而姿控推力器由于推力過小，難以對軌道產(chǎn)生較大影響，因而無法標定。此外，地面還常常根據(jù)軌控期間管路壓力變化，估算發(fā)動機推力，但該方法嚴重依賴發(fā)動機地面試車結(jié)果，當衛(wèi)星在軌的管路溫度、空間外熱流等環(huán)境發(fā)生變化時，通常該方法無法準確反映發(fā)動機推力特性。文獻[8-9]用遞推最小二乘法對航天器質(zhì)心參數(shù)進行標定，進而推算出推力器推力大小。缺點是需要建立在空間干擾力矩是偽隨機的條件下，保證多組推力器的噴氣時長盡可能接近，減少系統(tǒng)誤差。文獻[10]以連續(xù)小推力航天器為背景，提出了綜合考慮星載加速度計和推力器在軌標定的自主導航方案，但需要長期的軌道變化數(shù)據(jù)做反演迭代。此外還需建立高精度的引力場模型，在深空環(huán)境下各天體的引力攝動極為微弱，太陽系內(nèi)大天體的引力場建模難以達到迭代需求。文獻[11]描述了嫦娥一號探測器的在軌標定方法，在任務中使用了控前和控后的精密軌道數(shù)據(jù)以及加速度遙測，對沉底、軌控過程的速度增量進行了測量，同時對加速度計的刻度因數(shù)等參數(shù)進行了標定。另外將測量數(shù)據(jù)迭代進下一次軌控策略計算中，以提高控制精度，取得了良好的效果。該方法的缺陷是需要地面測定軌精度直接影響標定精度，在深空探測領(lǐng)域隨著航天器飛行距離越來越遠，測定軌精度逐步下降，發(fā)動機標定精度也隨之變差。文獻[12]針對一種海洋測高衛(wèi)星質(zhì)心測量的需求，提出了一種更符合工程實際情況的在軌估計衛(wèi)星質(zhì)心的改進算法，綜合考慮了實際推力器推力誤差及陀螺儀的測量誤差，采用總體最小二乘法求解，克服了以往算法不能適應實際工程中存在推力誤差的缺陷，具有一定的理論和工程實用價值。

針對深空任務的特殊性，探測器在飛行過程中通常需進行多次點火用于變軌、減速剎車等動作，在此過程中整器質(zhì)心發(fā)生劇烈變化，發(fā)動機推力參數(shù)嚴重影響變軌精度。且深空探測器飛行距離遠，軌道定位精度差，常規(guī)地球衛(wèi)星依靠的測定軌的標定方法在此無法適用。深空環(huán)境下溫度、光照處于時變狀態(tài)，無法為探測器提供穩(wěn)定的推進管路工作環(huán)境，因此也無法通過上一次的點火標定結(jié)果預測后續(xù)的推力器工作性能[13-15]。在深空環(huán)境下航天器的推進系統(tǒng)參數(shù)需依靠自身相關(guān)遙測進行辨識。

2 質(zhì)心和推力矢量估計算法

火星探測器在地火轉(zhuǎn)移階段正常飛行時保持對日姿態(tài)，由于此時已遠離地月系統(tǒng)，重力梯度和氣動力矩微乎其微，空間環(huán)境中僅有光壓力矩作用于航天器上。在對日姿態(tài)基準下探測器慣性姿態(tài)角一天僅變化不到0.2°，因此短期內(nèi)可認為探測器保持慣性姿態(tài)，且空間干擾力矩為常值。探測器通過周期性的噴氣卸載減少干擾力矩累積，每次卸載時滿足姿態(tài)動力學方程[16]。

(1)

式中：Tp為噴氣力矩；Te為空間干擾力矩；HB和Hw分別為整器和飛輪的角動量?；鹦翘綔y器在卸載前后的穩(wěn)定狀態(tài)下，整器慣性姿態(tài)角速度變化極小，忽略不計，主要由飛輪對角動量進行吸收，由此可根據(jù)飛輪轉(zhuǎn)速變化直接得到外力矩產(chǎn)生的角動量：

(2)

探測器保持對日固定姿態(tài)，空間干擾力矩短期內(nèi)近似常值，可通過穩(wěn)態(tài)時飛輪轉(zhuǎn)速變化曲線進行擬合，斜率即為干擾力矩幅值。由此解出每次噴氣卸載產(chǎn)生的角動量[17]。

噴氣卸載力矩為：

Tp=(rp-r0)×Fp

(3)

式中：rp為推力器喉部坐標；r0為探測器在軌質(zhì)心坐標；Fp為推力器推力矢量。推力器喉部坐標根據(jù)整器布局可事先確定，整器質(zhì)心和推力矢量為求解目標參數(shù)。根據(jù)各組推力器工作時角動量的不同變化，組成方程組聯(lián)合求解參數(shù)。

3 推力器布局

探測器的環(huán)繞器部分共安裝了1臺3 000 N發(fā)動機和6組標稱推力25 N的姿控推力器，用于飛行過程中軌道調(diào)整以及角動量卸載等操作。在衛(wèi)星底板-x面布置了4組姿控推力器，用于Y、Z軸控制，另外2組布置在環(huán)繞器+Z側(cè)控制X軸，推力器布局和局部坐標系定義如圖1所示。

圖1 推力器布局示意Fig.1 Layout diagram of thrusters

每組推力器包含2個噴嘴，安裝在同一個支架上，考慮到2個噴嘴距離較近且采用了共基準安裝，空間外部熱環(huán)境狀態(tài)接近，可認為二者的推力方向產(chǎn)生相同方向的偏差。推力器坐標系到衛(wèi)星本體系的標稱安裝矩陣用Cbp表示。由式(3)可知，推力器每次噴氣可得3組等式，質(zhì)心坐標為未知量，用3個參數(shù)表示坐標值。同樣的方式，推力矢量在空間中用推力大小|F0|、推力橫向角度偏差θ、推力縱向角度偏差φ三個參數(shù)表示。因此每個推力器噴氣用含6個未知參數(shù)的3個等式表示：

(4)

由于在軌飛行過程中各推力器每次噴氣受管路壓力、溫度、電磁閥閥芯運動特性等因素影響，推力大小存在一定波動，每次產(chǎn)生的噴氣沖量不完全相同，上式中參數(shù)推力大小|F0|為噴氣脈沖的平均推力。

顯然單個推力器噴氣得到的解算方程組，待求解參數(shù)數(shù)量大于線性方程數(shù)量，無法得到有效解析解。使用安裝在同一個支架上的2個推力器噴氣數(shù)據(jù)時，得到6個方程，待求解參數(shù)包括質(zhì)心坐標(3個參數(shù))、推力角度偏差(橫向和縱向偏差2個參數(shù))、推力大小(A和B推力器推力大小2個參數(shù))，共7個。以此類推，使用3組推力器數(shù)據(jù)時，可得到18個方程，而由于質(zhì)心參數(shù)相同，待求解參數(shù)只有15個，因此可得到有效的解析數(shù)值。

第m個推力器噴氣產(chǎn)生的姿態(tài)擾動動力學方程為：

(5)

當標稱推力方向與坐標軸平行時，在該坐標軸方向上的質(zhì)心坐標的解算方程易出現(xiàn)奇異，造成解析精度不高。為保證解析精度，選取的參與計算的3組推力器，至少有1組不與坐標軸重合。

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

采用中國首顆火星探測器天問一號在巡航段的噴氣卸載數(shù)據(jù)進行分析。此時探測器距地球超過3×106km，并且尚未到達火星影響球，空間外力矩僅有光壓力矩。姿態(tài)保持對日巡航姿態(tài)，姿態(tài)角速度變化極小，在數(shù)天的時間尺度內(nèi)可以視作慣性姿態(tài)，是研究空間光壓力矩的理想測試平臺。

表1給出了探測器6組推力器喉部的坐標以及標稱推力矢量，即式(4)(5)中的參數(shù)rp。

表1 推力器布局參數(shù)

根據(jù)布局關(guān)系，Y軸、Z軸推力器裝在底面，標稱安裝矩陣Cbp為單位陣。X軸推力器裝在側(cè)面，根據(jù)布局關(guān)系可知正負兩組推力器安裝矩陣為：

選取2020年8月的6次不同推力器的噴氣卸載數(shù)據(jù)，噴氣前后三軸飛輪轉(zhuǎn)速的變化情況如圖2～圖7所示。

圖2 X軸(A2推力器)噴氣前后飛輪轉(zhuǎn)速變化Fig.2 Flywheel rotation changes of X axial before and after ignition (A2 thruster)

圖3 X軸(B2推力器)噴氣前后飛輪轉(zhuǎn)速變化Fig.3 Flywheel rotation changes of X axial before and after ignition (B2 thruster)

圖4 Y軸(A4推力器)噴氣前后飛輪轉(zhuǎn)速變化Fig.4 Flywheel rotation changes of Y axial before and after ignition (A4 thruster)

圖5 Y軸(B4推力器)噴氣前后飛輪轉(zhuǎn)速變化Fig.5 Flywheel rotation changes of Y axial before and after ignition (B4 thruster)

根據(jù)每次推力器噴氣前后的飛輪轉(zhuǎn)速變化，利用式(2)計算噴氣產(chǎn)生的角動量，同時需扣除空間干擾力矩的影響。由于火星探測器在巡航段僅受光壓力矩這一種空間力矩影響，在數(shù)天時間內(nèi)光壓力矩可認為是常值，因此可用一段時間的飛輪轉(zhuǎn)速曲線線性擬合得到，即式(2)中的Te，結(jié)果見表2中空間力矩這一項。

圖6 Z軸(A10推力器)噴氣前后飛輪轉(zhuǎn)速變化Fig.6 Flywheel rotation changes of Z axial before and after ignition (A10 thruster)

圖7 Z軸(B10推力器)噴氣前后飛輪轉(zhuǎn)速變化Fig.7 Flywheel rotation changes of Z axial before and after ignition (B10 thruster)

表2 噴氣卸載前后干擾力矩計算結(jié)果

從上表看出，不同時段空間力矩有一定差異。在排除空間力矩影響后，每次噴氣卸載引起的整器角動量變化也有較大不同，即使同一組推力器的兩個分支，由于推力作用點，產(chǎn)生的角動量方向有一定區(qū)別。

根據(jù)式(5)，僅有質(zhì)心r0和推力矢量Fp未知(包括推力大小|F0|、推力橫向角度偏差θ、推力縱向角度偏差φ三個參數(shù))，利用6次噴氣卸載數(shù)據(jù)得到18個等式，而由于質(zhì)心參數(shù)相同，待求解參數(shù)只有15個，最小二乘違逆求解推力器和質(zhì)心偏差方向，結(jié)果如表3、表4所示。

對比任務指標，各推力器偏斜量均小于0.6°的任務指標，滿足指標需求。

和理論質(zhì)心位置相比，由于燃料裝填量因素不可控，以及器上電纜等附件質(zhì)量難以精確計算，質(zhì)心在X軸上有一定偏差，其余兩軸質(zhì)心得到較好的控制。

表3 推力器偏斜量計算結(jié)果

表4 質(zhì)心偏差計算結(jié)果

5 結(jié)論

1)本文針對深空探測器飛行距離遠、無法通過測定軌精確標定推力器參數(shù)的難點，利用探測器飛行過程中空間干擾力矩穩(wěn)定的特點，提出了一種利用噴氣卸載前后飛輪轉(zhuǎn)速變化辨識整器質(zhì)心和卸載推力器推力矢量的方案。和已有方法相比，脫離了額外的器上單機、地面測定軌的約束，任務面影響小。

2)在使用天問一號火星探測器的遙測數(shù)據(jù)進行推進系統(tǒng)參數(shù)辨識的過程中，由于A4、B4、A10、B10推力器方向與X軸平行，質(zhì)心解析產(chǎn)生一定耦合，YZ軸質(zhì)心具有較高的辨識精度，X軸質(zhì)心辨識精度稍差。若要進一步提高精度，需在各軸上均進行噴氣，在6個姿態(tài)維度上進行解析。

3)質(zhì)心估計和推力器方向偏斜辨識可為后續(xù)軌控任務的點火方向制定、燃料預算提供輸入依據(jù)，提高軌控精度，降低軌控風險。

4)后續(xù)主要針對推力軸與質(zhì)心測量軸耦合問題、復雜空間干擾力矩影響問題開發(fā)解耦的辨識方法，進一步擴展本方法在航天器推進系統(tǒng)辨識上的應用。