基于夏氏最小二乘的軌道控制力系數辨識

王西京袁博孔大林卞燕山

1. 西北工業大學 航天學院, 西安 710072 2. 航天器在軌故障診斷與維修重點實驗室,西安 710043 3. 西安衛星測控中心, 西安 710043

航天器長期在軌運行，存在軌道攝動、初始狀態誤差、控制誤差等因素，使得實際軌道偏離標稱的設計軌道，不能完成預定任務[1-2]。為使衛星長期運行在理論設計的標稱軌道上，以及滿足衛星用戶新的應用需求，需要對衛星實施軌道機動。隨著衛星組網、編隊、交會對接、天地往返等多種主動變軌控制技術的日益成熟，以及空間目標碰撞規避次數的增加，還需要通過衛星軌道控制完成更復雜的空間任務[3]。

衛星軌道控制推力器推力是以推進劑貯箱壓力為自變量的多項式函數，其系數為裝訂常數，在衛星的整個壽命期保持不變[4-5]。隨著衛星軌控次數的增加，推進劑貯箱壓力不斷降低，以理論推力為基礎計算的衛星速度增量與實際速度增量出現誤差，導致軌道控制精度降低。傳統的軌道控制速度增量計算以前一次推力器標定系數作為計算依據[6]，并且推力系數維持不變[7]，沒有有效利用歷次軌道控制數據進行優化計算，導致理論速度增量與實際速度增量誤差不可測，給軌道控制策略制定帶來了一定困難。

本文統計分析在軌管理的典型航天器平臺及其發動機的軌道控制歷史數據，分析軌道控制理論和在軌控制數據擬合建立軌道控制經驗模型，用當前可測量的系統輸入和輸出預測系統輸出的未來演變，得到不同工作情況下實際軌道控制誤差與控制參數及其他主要影響因素之間關系的經驗公式，為軌道控制策略決策提供參考。選取軌道半長軸控制量300 m以上和300 m以下的兩類近地衛星，對其軌道控制歷史數據進行分析，經實際數據測試，采用夏氏改良法進行推力系數擬合后預測的速度變化量精度較高。

1 最小二乘夏氏改良法

最小二乘法的特點在于系統的輸入和輸出信號反復過濾。夏氏法是一種交替的廣義最小二乘法求解技術，由夏天長(T.C.Hsia)提出來的，它不需要數據反復過濾，因而計算效率較高。這種方法能夠克服最小二乘估計中的有偏估計問題，而且由這種方法所導出的計算方法也比較簡單，分為夏氏偏差修正法和夏氏改良法[8]。

當考慮系統噪聲影響時，系統的差分方程可表示為：

(1)

(2)

求式(2)中f的最小二乘估計:

(3)

于是有

(4)

這種算法稱之為夏氏改良法，夏氏改良法改進了收斂效率，工程上可用。

2 控制力系數辨識模型

根據動量與沖量之間的關系可得[9]：

(6)

(7)

式中：F為軌控前推力；Ff為軌控后推力。

(8)

式中：Di(i=1,2,3,4)為擬合多項式系數；P為軌控前儲箱壓力(Mpa)；Px為儲箱1、2當前推進劑密度(kg/cm3)；Pf為軌控后儲箱壓力，腳標x代表燃料儲箱編號(x=1，2)，推力器工作時儲箱1和儲箱2根據壓強變化互相切換[10]，

(9)

式中：VTx為儲箱容積；

ρx=1025.5-0.875×(Tx-273.15)

mx為燃料剩余質量，

ISx為比沖與壓力的關系，

IEx為本次機動推力器所需提供壓力所對應的總沖[11]，

其中衛星當前質量

msat=msat0-m01+m1-m02+m2

式中：i代表第i組測量值。

將式(10)展開后可得到：

(11)

根據式(11)，可得到測量矩陣Φ、待求參數向量θ、輸出向量y：

θ=[D0D1D2D3]T

y=[ΔV1xΔV2x… ΔVnx]T

考慮系統噪聲ξ影響時，定義噪聲擬合系數f=[f1f2f3f4]T，系統噪聲ξ的擬合公式為：

(12)

噪聲ε為不相關的隨機序列，故可用最小二乘法得到系數f的一致無偏估計[13]。若有n組測量數據，考慮到計算效率，采用最小二乘夏氏改良法[14]，則有

(13)

式中：e為殘差；Ω的具體表達式如下：

綜上，參照最小二乘夏氏改良法的迭代計算步驟即可求得估計系數[15]。

3 仿真驗證

軌道控制半長軸300 m對應的速度變化量約為0.15，在工程經驗上Δa=300 m為分界線，推力穩定性和保持環會有差別。軌道高度低的衛星軌跡保持環小，保持周期短，控制量大；軌道高度高的衛星軌跡保持環大，保持周期長，控制量小。所以分別選取軌道半長軸控制量300 m以上和300 m以下的兩類近地衛星，對其軌道控制歷史數據進行分析，得到的結論具有普遍性。

3.1 算例一

選取軌道高度300 km以上某低軌衛星共44批次控制數據計算推力擬合系數，其中前20批次軌控使用貯箱1推進劑，后24批次軌控使用貯箱2推進劑，相關仿真參數如表1所示。

綜上，可使用數據共44組，每次仿真去掉其中一組數據，用余下43組數據作為樣本，去掉的那組數據作為精度驗證標準，進行仿真，仿真結果如圖1所示。

圖1為理論速度變化量、預測速度變化量與標定速度變化量，其趨勢幾乎一致。對預測誤差百分比取絕對值，統計可得其誤差百分比分布。以標定速度變化量為近似真值，分別計算理論誤差和預測誤差，統計可得其相對誤差百分比分布如圖2所示，最大誤差為4.6%，平均誤差為1.39%，采用夏氏法進行推力系數擬合后預測的速度變化量精度較高。

表1 仿真參數

圖1 速度變化量分布Fig.1 Variation distribution of velocity

圖2 誤差百分比分布Fig.2 Percentage distribution of error

3.2 算例二

選取軌道高度300 km以下某低軌衛星共13個批次控制數據計算推力擬合系數，由于算例二所選衛星的實測軌控數據樣本偏少，且可直接提供軌控后的貯箱壓強，因此在估算模型中無需再重新計算軌控后貯箱壓強，直接代入即可。相關仿真參數如表2所示。

綜上，可使用的數據共有13組，每次仿真去掉其中一組數據，用余下的12組數據作為樣本，去掉的那組數據作為精度驗證標準，進行仿真，仿真結果如表3所示。

對上述預測誤差百分比取絕對值，統計可得其誤差百分比分布如圖3所示。

表2 仿真參數

表3 仿真結果

從表3和圖3可以看出，最大誤差為1.61%，平均誤差為0.85%，采用夏氏法進行推力系數擬合后預測的速度變化量精度較高。

參考第8～13次控制前的軌道參數(見表4)，將預估的速度增量代入軌道仿真模型中，進行仿真得到的結果與實際軌控后的半長軸對比如表5所示。

從表5的統計結果可以看出，所修正后的控制系數計算得到的速度增量代入軌道外推模型，得到的半長軸與實際半長軸對比都小于20 m，最大誤約19 m，最小誤差僅1.6 m。

圖3 誤差百分比分布Fig.3 Percentage distribution of error

表4 初始仿真參數

表5 仿真結果對比

4 結束語

通過分析歷次軌控理論速度變化量、實際速度變化量、推力器推力持續時間、軌控推力系數之間的相關性，采用最小二乘夏氏改良法估計表征推力的擬合系數，并運用夏氏法進行求解，建立一個能模仿推力器在軌工作情況的模型，用當前可測量的系統輸入和輸出預測系統輸出的未來演變，具有預測精度高的特點，且整個計算過程不涉及標定系數，可以降低主觀因素的影響，對軌道控制實施具有參考意義。