基于精確星歷最小二乘擬合的軌道快速遞推方法*

車 征

（陜西星邑空間技術(shù)有限公司 西安 710043）

1 引言

高精度的軌道快速遞推模型可為衛(wèi)星敏感器、控制器等部件的在軌自主工作提供位置輸入信息。由于星上計(jì)算資源和運(yùn)算效率限制，無(wú)法使用高精度數(shù)值積分方式進(jìn)行軌道預(yù)報(bào)。傳統(tǒng)的攝動(dòng)分析方法采用一階或二階攝動(dòng)理論［1～2］，雖然能夠進(jìn)行較快速的平均根數(shù)與瞬時(shí)根數(shù)的轉(zhuǎn)換和軌道遞推，但需要計(jì)算大量的公式推導(dǎo)工作，且在攝動(dòng)模型項(xiàng)較多時(shí)運(yùn)算仍然較為復(fù)雜。GPS 星歷外推使用了一種參數(shù)擬合的方式，將攝動(dòng)規(guī)律轉(zhuǎn)為有限個(gè)參數(shù)描述的近似項(xiàng)［3～5］。參數(shù)個(gè)數(shù)從最初的16參數(shù)、發(fā)展出18 參數(shù)、22 參數(shù)等［6～9］。這種方法具有很高的擬合精度，但缺點(diǎn)是僅僅適用于短時(shí)間內(nèi)的位置解算，當(dāng)擬合時(shí)段變長(zhǎng)則誤差急劇增大，且沒(méi)有包含速度信息，不適合用于長(zhǎng)時(shí)間軌道遞推［10～13］。還有一種方法是多項(xiàng)式插值方式［14～15］，在擬合區(qū)間內(nèi)能夠?qū)崿F(xiàn)很高精度，但多項(xiàng)式階數(shù)較高時(shí)計(jì)算量較大，且用于預(yù)報(bào)存在發(fā)散趨勢(shì)。

本文從瞬根數(shù)變化規(guī)律中提取軌道根數(shù)各分量的頻率信息，使用含少量待定系數(shù)的三角函數(shù)和線性函數(shù)組合關(guān)系式描述軌道根數(shù)的變化趨勢(shì)，使用精確動(dòng)力學(xué)預(yù)報(bào)結(jié)果通過(guò)最小二乘擬合的方式確定系數(shù)值。避免了復(fù)雜的解析公式運(yùn)算，同時(shí)保證了與精確動(dòng)力學(xué)模型預(yù)報(bào)結(jié)果的吻合度。仿真算例表明擬合結(jié)果可以保證24h 平均預(yù)報(bào)誤差小于400m，且計(jì)算耗時(shí)小。

2 擬合模型

基本思想為根據(jù)瞬根周期規(guī)律進(jìn)行傅里葉分解，找出每個(gè)根數(shù)變化的短周期頻率，分別指定擬合函數(shù)，對(duì)每個(gè)軌道根數(shù)單獨(dú)擬合。以第一類無(wú)奇點(diǎn)根數(shù)為變量，定義各個(gè)分量的變化規(guī)律如下：1）a的周期為T/2，還包括隨時(shí)間線性變化項(xiàng)。2）ex和ey的周期為T/3和T。3）i的周期為T/2和1/2天。4）Ω的周期為T/2，還包括隨時(shí)間線性變化項(xiàng)。5）λ的周期為T和1/2天，還包括隨時(shí)間線性變化項(xiàng)。T為平均軌道周期。因此設(shè)置擬合函數(shù)如下式：

一共有30 個(gè)參數(shù)。對(duì)六根數(shù)分別擬合，設(shè)待擬合的參數(shù)符號(hào)為x（代表a，ex，ey，i，Ω，λ中的一個(gè)），待擬合參數(shù)的初值為x0到xn，則：

利用不同時(shí)間共m 組x 的精確動(dòng)力學(xué)模型預(yù)報(bào)數(shù)據(jù)，有

其中L=xpre-x為按照式（1）和時(shí)間序列預(yù)報(bào)的xpre與真實(shí)動(dòng)力學(xué)積分得到的x之間的偏差。

則根據(jù)L可以反算出x0到xn應(yīng)修正量和修正后的待擬合變量為

其中（a，ex，ey，i，Ω，λ）對(duì)各自待擬合變量的偏導(dǎo)數(shù)為

重復(fù)式（3）～（6）直到計(jì)算的x0到xn應(yīng)修正量接近0，即表示待擬合變量收斂。由于最小二乘方法需要一個(gè)較為準(zhǔn)確的初值才能穩(wěn)定收斂，下節(jié)給出初值的計(jì)算方法。

3 擬合初值計(jì)算方法

3.1 a的系數(shù)初值計(jì)算

對(duì)于給定的數(shù)值積分獲得的一系列時(shí)刻的a（ti），i=1 到N，首先根據(jù)大氣衰減較小的假設(shè)，取a3=0。然后根據(jù)式（1），可知

當(dāng)tf遠(yuǎn)大于T時(shí)，可以認(rèn)為

其中Δt為數(shù)據(jù)時(shí)間間隔，且(N-1)ΔtΔn是一個(gè)小量，則

T作為常值用于隨后各個(gè)軌道根數(shù)系數(shù)的擬合過(guò)程中。根據(jù)式（1）的性質(zhì)，可知：

則：

3.2 ex，ey 的系數(shù)初值計(jì)算

則

求解方程組得：

同理利用

可計(jì)算：

ey的初值計(jì)算與式（15）到式（18）完全相同。

3.3 i的系數(shù)初值計(jì)算

同樣地，i0初值直接取i 的平均值，。根據(jù)和i(T)的特性，有：

3.4 Ω的系數(shù)初值計(jì)算

根據(jù)：

可得：

3.5 λ 的系數(shù)初值

首先λ3初值取為式（10）中的軌道角速度n。

根據(jù)下式：

得：

再根據(jù)下式：

得：

4 測(cè)試結(jié)果

初始瞬根為［7130km，0，98°，0°，0°，360°］。則擬合一天的參數(shù)為

a=7120.9688+9.0790 sin(2πt/T+1.5751)-3.3803e-08t

ex=-0.00046978+0.0007439 sin(6πt/T+1.5752)+0.0002941 sin(2πt/T-1.5749)

ey=2.2291e-5+0.0007435 sin(6πt/T+0.004545)+0.0009327 sin(2πt/T-3.1340)

i=1.7105+8.9871e-05 sin(4πt/T-1.5839)+2.7833e-05 sin(2πt/43200-0.6091)Ω=2.0126e-05+9.0666e-05 sin(4πt/T-3.1393)+1.9014e-07t

λ=0.0001174+0.0009410 sin(2πt/T-3.0928)+0.001049t+9.7748e-05 sin(2πt/43200-2.1751)+9.4583e-05 sin(4πt/T+0.008591)

對(duì)應(yīng)的殘差換算到三維位置偏差見(jiàn)圖1，平均總誤差小于400m。設(shè)去除3 倍周期以內(nèi)的短周期項(xiàng)后的剩余長(zhǎng)期項(xiàng)（常數(shù)項(xiàng)、長(zhǎng)周期項(xiàng)和時(shí)間線性項(xiàng)）為平根數(shù)，與瞬根的對(duì)比如圖2～圖7。

圖1 位置殘差

圖2 半長(zhǎng)軸長(zhǎng)期項(xiàng)提取

圖3 ex長(zhǎng)期項(xiàng)提取

圖4 ey長(zhǎng)期項(xiàng)提取

圖5 i長(zhǎng)期項(xiàng)提取圖

圖6 Ω長(zhǎng)期項(xiàng)提取

圖7 λ 長(zhǎng)期項(xiàng)提取

可見(jiàn)擬合獲得的平根能夠很好地反映軌道根數(shù)長(zhǎng)期變化。

5 結(jié)語(yǔ)

本文使用三角函數(shù)和時(shí)間線性函數(shù)來(lái)逼近低軌近圓軌道的軌道根數(shù)變化趨勢(shì)，給出了最小二乘擬合函數(shù)系數(shù)的算法和流程，并設(shè)計(jì)了快速獲得擬合系數(shù)初值的計(jì)算方法。仿真實(shí)例表明，擬合方法能夠快速獲得較高精度的軌道根數(shù)遞推函數(shù)，24h數(shù)據(jù)擬合精度達(dá)到400m，且計(jì)算量很小。算法可用于星上快速軌道遞推和瞬平根數(shù)轉(zhuǎn)換。