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

王西京袁博孔大林卞燕山

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

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

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

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

1 最小二乘夏氏改良法

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

當(dāng)考慮系統(tǒng)噪聲影響時，系統(tǒng)的差分方程可表示為：

(1)

(2)

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

(3)

于是有

(4)

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

2 控制力系數(shù)辨識模型

根據(jù)動量與沖量之間的關(guān)系可得[9]：

(6)

(7)

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

(8)

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

(9)

式中：VTx為儲箱容積；

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

mx為燃料剩余質(zhì)量，

ISx為比沖與壓力的關(guān)系，

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

其中衛(wèi)星當(dāng)前質(zhì)量

msat=msat0-m01+m1-m02+m2

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

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

(11)

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

θ=[D0D1D2D3]T

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

考慮系統(tǒng)噪聲ξ影響時，定義噪聲擬合系數(shù)f=[f1f2f3f4]T，系統(tǒng)噪聲ξ的擬合公式為：

(12)

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

(13)

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

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

3 仿真驗證

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

3.1 算例一

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

綜上，可使用數(shù)據(jù)共44組，每次仿真去掉其中一組數(shù)據(jù)，用余下43組數(shù)據(jù)作為樣本，去掉的那組數(shù)據(jù)作為精度驗證標(biāo)準(zhǔn)，進行仿真，仿真結(jié)果如圖1所示。

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

表1 仿真參數(shù)

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

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

3.2 算例二

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

綜上，可使用的數(shù)據(jù)共有13組，每次仿真去掉其中一組數(shù)據(jù)，用余下的12組數(shù)據(jù)作為樣本，去掉的那組數(shù)據(jù)作為精度驗證標(biāo)準(zhǔn)，進行仿真，仿真結(jié)果如表3所示。

對上述預(yù)測誤差百分比取絕對值，統(tǒng)計可得其誤差百分比分布如圖3所示。

表2 仿真參數(shù)

表3 仿真結(jié)果

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

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

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

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

表4 初始仿真參數(shù)

表5 仿真結(jié)果對比

4 結(jié)束語

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