基于J2 攝動(dòng)項(xiàng)的Runge-Kutta、Hermite 數(shù)值軌道預(yù)報(bào)方法?

張政超 朱 峰 郭偉鋒 周 磊

（63891部隊(duì) 洛陽(yáng) 471000）

1 引言

近年來(lái)，對(duì)衛(wèi)星軌道的研究方法層出不窮［1～5］，對(duì)衛(wèi)星軌道的描述方法和衛(wèi)星軌道的預(yù)測(cè)更是持續(xù)不斷的研究熱點(diǎn)［6～8］。一方面隨著人們對(duì)開普勒三定律的深入認(rèn)識(shí)，對(duì)衛(wèi)星軌道的描述使用慣性坐標(biāo)系下的模型普遍被接受，另一方面，由于星上計(jì)算資源的有限，利用當(dāng)前的位置、速度信息準(zhǔn)確、快速、高效地預(yù)測(cè)衛(wèi)星此后一段時(shí)間的軌道信息，為星上任務(wù)設(shè)備的姿態(tài)控制提供參考依據(jù)便變得十分有意義和重要。

雖然星上任務(wù)設(shè)備對(duì)衛(wèi)星的軌道信息的精度要求不是很高，但在無(wú)GNSS信息，僅靠上注信息的情況下，隨著時(shí)間的增長(zhǎng)，衛(wèi)星軌道的精度很快就會(huì)下降，甚至斷崖式下跌［9～12］，解決的辦法是選用合適的軌道預(yù)測(cè)方法、定時(shí)上注、多源信息融合等。本文僅對(duì)預(yù)測(cè)方法進(jìn)行分析和仿真，解決了一定精度下對(duì)衛(wèi)星軌道的預(yù)報(bào)問(wèn)題，其他不多贅述。

2 衛(wèi)星軌道的表示方法及坐標(biāo)轉(zhuǎn)換

2.1 表示方法

對(duì)衛(wèi)星軌道信息的描述主要分時(shí)間系統(tǒng)和坐標(biāo)系統(tǒng)。

時(shí)間系統(tǒng)包括UTC 時(shí)間和北京時(shí)間。坐標(biāo)系統(tǒng)包括WGS84坐標(biāo)系下的笛卡爾坐標(biāo)（3維位置、3維速度）、J2000.0 坐標(biāo)系下的軌道六根數(shù)（長(zhǎng)半軸、偏心率、傾角、升交點(diǎn)赤經(jīng)、近地點(diǎn)幅角、真近點(diǎn)角）或笛卡爾坐標(biāo)（3 維位置、3 維速度）。一般來(lái)說(shuō)，測(cè)控時(shí)選用的WGS84 坐標(biāo)系下的笛卡爾坐標(biāo)，研究軌道性質(zhì)選用的是J2000.0坐標(biāo)系下的軌道六根數(shù)或笛卡爾坐標(biāo)，三者之間可以互換。

2.2 WGS84坐標(biāo)系轉(zhuǎn)換到J2000.0坐標(biāo)系

J2000.0坐標(biāo)系與WGS84坐標(biāo)系主要在方面不同，WGS84 坐標(biāo)系（協(xié)議地球坐標(biāo)系）轉(zhuǎn)換到J2000.0 坐標(biāo)系（協(xié)議天球坐標(biāo)系）需進(jìn)行極移、格林尼治恒星時(shí)、章動(dòng)、歲差四方面的修正。其中三個(gè)過(guò)渡坐標(biāo)系依次為瞬時(shí)地球坐標(biāo)系、瞬時(shí)真天球坐標(biāo)系、瞬時(shí)平天球坐標(biāo)系。

經(jīng)上述過(guò)程，用Matlab 軟件編程后，可以實(shí)現(xiàn)從WGS84 坐標(biāo)系轉(zhuǎn)換到J2000.0 坐標(biāo)系下的坐標(biāo)變換。從J2000.0 坐標(biāo)系轉(zhuǎn)換到WGS84 坐標(biāo)系是上述過(guò)程的逆變換。經(jīng)與國(guó)內(nèi)“洞察者”衛(wèi)星軟件、衛(wèi)星導(dǎo)航系統(tǒng)坐標(biāo)轉(zhuǎn)換軟件、STK 軟件相比對(duì)，轉(zhuǎn)換的位置精度可達(dá)0.1m，速度精度可達(dá)0.001m/s。

3 衛(wèi)星軌道預(yù)報(bào)方法和流程

如前文所言，衛(wèi)星軌道預(yù)報(bào)的意義是為星上任務(wù)設(shè)備提供衛(wèi)星姿態(tài)參考，其內(nèi)容包括對(duì)此后一段時(shí)期的衛(wèi)星位置、速度的預(yù)報(bào)。具體選擇哪一種預(yù)報(bào)算法取決于星上的計(jì)算資源情況，需要在計(jì)算資源和預(yù)報(bào)精度之間取得平衡，從而得到合適的預(yù)報(bào)結(jié)果。

衛(wèi)星軌道預(yù)報(bào)的方法分為解析法和數(shù)值法等。文獻(xiàn)［3］詳細(xì)論述了二體運(yùn)動(dòng)的解析解，并針對(duì)第一類無(wú)奇點(diǎn)根數(shù)提出了簡(jiǎn)化的J2項(xiàng)攝動(dòng)的解析解，取得了不錯(cuò)的應(yīng)用效果，特別是在星上計(jì)算資源有限的情況下，能夠在要求的精度內(nèi)預(yù)測(cè)衛(wèi)星軌道。

文獻(xiàn)［1］則從衛(wèi)星軌道的瞬時(shí)根數(shù)出發(fā)，推導(dǎo)出了一種更為通用的衛(wèi)星軌道預(yù)測(cè)方法，其精度也比文獻(xiàn)［3］顯著提高，但對(duì)于第一類無(wú)奇點(diǎn)根數(shù)則無(wú)法取得令人滿意的預(yù)測(cè)結(jié)果。

如此，可采取了多次上注、根據(jù)實(shí)際軌道精度需求選取合理的上注頻次。

4 基于J2 攝動(dòng)項(xiàng)的Runge-Kutta、Hermite數(shù)值軌道預(yù)報(bào)方法

4.1 基于J2攝動(dòng)項(xiàng)的Runge-Kutta 數(shù)值軌道預(yù)報(bào)方法

為得到更高精度的衛(wèi)星軌道預(yù)測(cè)結(jié)果，相比之下，數(shù)值方法能夠取得更為理想的結(jié)果，即間接引用泰勒展開，用積分區(qū)間上若干右函數(shù)值的線性組合確定相應(yīng)的系數(shù)。一般來(lái)說(shuō)，4 階Runge_Kutta就能夠取得較高的精度要求。現(xiàn)推導(dǎo)如下。

在J2000.0 坐標(biāo)系下的笛卡爾坐標(biāo)系中，設(shè)遞推量rn的表示方式為

各項(xiàng)分別表示衛(wèi)星當(dāng)前的三維位置和三維速度信息，其遞推公式為

其中h 為積分步長(zhǎng)，ki（i=1，2，3，4）為過(guò)程變量，其表達(dá)式為

式（3）中的f函數(shù)輸出為

其中，μ為地球引力常數(shù)，取3.986×1014m3/m2，。

其中J2為地球兩次帶諧系數(shù)，取1.082×10-3，Re為地球平均半徑，取6378137m。由此可遞推得到衛(wèi)星軌道的位置、速度信息。基于J2攝動(dòng)項(xiàng)的Runge-Kutta數(shù)值預(yù)報(bào)如圖1、圖2所示。

圖1 基于J2攝動(dòng)項(xiàng)的Runge-Kutta數(shù)值位置預(yù)報(bào)

圖2 基于J2攝動(dòng)項(xiàng)的Runge-Kutta數(shù)值速度預(yù)報(bào)

4.2 基于Hermite數(shù)值軌道預(yù)報(bào)方法

設(shè)t0+nh時(shí)刻（0

對(duì)n分別求一次導(dǎo)、二次導(dǎo)：

分別令n=0，n=1，則有

回代到式（7），則有

對(duì)式（12）中n求導(dǎo)，有

其中：

由此可遞推得到衛(wèi)星軌道的位置、速度信息。基于Hermite數(shù)值衛(wèi)星軌道預(yù)報(bào)如圖3、圖4所示。

圖3 基于Hermite數(shù)值位置預(yù)報(bào)

圖4 基于Hermite數(shù)值速度預(yù)報(bào)

5 軌道預(yù)報(bào)誤差分析

以STK 軟件的HPOP 算法為真值，對(duì)基于J2攝動(dòng)項(xiàng)的Runge-Kutta 數(shù)值軌道預(yù)報(bào)誤差分別如圖5、圖6 所示，基于Hermite 數(shù)值軌道預(yù)報(bào)誤差分別如圖7、圖8所示。

圖5 基于J2攝動(dòng)項(xiàng)的Runge-Kutta數(shù)值位置預(yù)報(bào)誤差

圖6 基于J2攝動(dòng)項(xiàng)的Runge-Kutta數(shù)值速度預(yù)報(bào)誤差

圖7 基于Hermite數(shù)值預(yù)報(bào)位置誤差

圖8 基于Hermite數(shù)值預(yù)報(bào)速度誤差

由圖5、圖6 可以看出，基于J2攝動(dòng)項(xiàng)的Runge-Kutta 數(shù)值位置預(yù)報(bào)各向誤差在12000s 內(nèi)均不超過(guò)350m，基于J2攝動(dòng)項(xiàng)的Runge-Kutta 數(shù)值速度預(yù)報(bào)各向誤差在12000s 內(nèi)均不超過(guò)0.4m/s。同樣以STK 軟件的HPOP 算法為真值由圖7、圖8 可以看出，基于Hermite 數(shù)值位置預(yù)報(bào)誤差在1800s 內(nèi)各向最大值約為6000m，基于Hermite 數(shù)值速度預(yù)報(bào)誤差在1800s 內(nèi)各向最大值約為12m/s，且呈“鐘形”發(fā)散。由此可見基于J2攝動(dòng)項(xiàng)的Runge-Kutta數(shù)值預(yù)報(bào)精度遠(yuǎn)高于基于Hermite 數(shù)值預(yù)報(bào)精度，實(shí)際使用時(shí)根據(jù)需要選擇合適的預(yù)報(bào)算法。

6 結(jié)語(yǔ)

本文分析了WGS84 坐標(biāo)與J2000.0 系下的坐標(biāo)互換，并重點(diǎn)研究了基于J2攝動(dòng)項(xiàng)的Runge-Kutta 數(shù)值預(yù)報(bào)、基于Hermite 數(shù)值預(yù)報(bào)的方法、精度等。一般來(lái)說(shuō)，衛(wèi)星軌道預(yù)報(bào)的數(shù)值計(jì)算方法的精度遠(yuǎn)高于解析方法，而且避開了第一類奇點(diǎn)導(dǎo)致的預(yù)報(bào)精度惡化現(xiàn)象［4～5］，但精度要求不高及計(jì)算資源有限時(shí)，解析方法亦不失為一種分析衛(wèi)星軌道特性的重要工具，比如計(jì)算軌道六根數(shù)的變化規(guī)律等，這也是下一步繼續(xù)努力的方向。