深空探測(cè)器在軌推力器參數(shù)辨識(shí)方法

信思博，趙訓(xùn)友，鄭藝裕，李綠萍

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

1 引言

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

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

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

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

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

(1)

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

(2)

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

噴氣卸載力矩為：

Tp=(rp-r0)×Fp

(3)

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

3 推力器布局

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

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

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

(4)

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

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

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

(5)

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

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

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

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

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

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

選取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)計(jì)算噴氣產(chǎn)生的角動(dòng)量，同時(shí)需扣除空間干擾力矩的影響。由于火星探測(cè)器在巡航段僅受光壓力矩這一種空間力矩影響，在數(shù)天時(shí)間內(nèi)光壓力矩可認(rèn)為是常值，因此可用一段時(shí)間的飛輪轉(zhuǎn)速曲線線性擬合得到，即式(2)中的Te，結(jié)果見表2中空間力矩這一項(xiàng)。

圖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 噴氣卸載前后干擾力矩計(jì)算結(jié)果

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

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

對(duì)比任務(wù)指標(biāo)，各推力器偏斜量均小于0.6°的任務(wù)指標(biāo)，滿足指標(biāo)需求。

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

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

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

5 結(jié)論

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

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

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

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