基于J2 攝動項的Runge-Kutta、Hermite 數值軌道預報方法?

張政超 朱 峰 郭偉鋒 周 磊

（63891部隊 洛陽 471000）

1 引言

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

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

2 衛星軌道的表示方法及坐標轉換

2.1 表示方法

對衛星軌道信息的描述主要分時間系統和坐標系統。

時間系統包括UTC 時間和北京時間。坐標系統包括WGS84坐標系下的笛卡爾坐標（3維位置、3維速度）、J2000.0 坐標系下的軌道六根數（長半軸、偏心率、傾角、升交點赤經、近地點幅角、真近點角）或笛卡爾坐標（3 維位置、3 維速度）。一般來說，測控時選用的WGS84 坐標系下的笛卡爾坐標，研究軌道性質選用的是J2000.0坐標系下的軌道六根數或笛卡爾坐標，三者之間可以互換。

2.2 WGS84坐標系轉換到J2000.0坐標系

J2000.0坐標系與WGS84坐標系主要在方面不同，WGS84 坐標系（協議地球坐標系）轉換到J2000.0 坐標系（協議天球坐標系）需進行極移、格林尼治恒星時、章動、歲差四方面的修正。其中三個過渡坐標系依次為瞬時地球坐標系、瞬時真天球坐標系、瞬時平天球坐標系。

經上述過程，用Matlab 軟件編程后，可以實現從WGS84 坐標系轉換到J2000.0 坐標系下的坐標變換。從J2000.0 坐標系轉換到WGS84 坐標系是上述過程的逆變換。經與國內“洞察者”衛星軟件、衛星導航系統坐標轉換軟件、STK 軟件相比對，轉換的位置精度可達0.1m，速度精度可達0.001m/s。

3 衛星軌道預報方法和流程

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

衛星軌道預報的方法分為解析法和數值法等。文獻［3］詳細論述了二體運動的解析解，并針對第一類無奇點根數提出了簡化的J2項攝動的解析解，取得了不錯的應用效果，特別是在星上計算資源有限的情況下，能夠在要求的精度內預測衛星軌道。

文獻［1］則從衛星軌道的瞬時根數出發，推導出了一種更為通用的衛星軌道預測方法，其精度也比文獻［3］顯著提高，但對于第一類無奇點根數則無法取得令人滿意的預測結果。

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

4 基于J2 攝動項的Runge-Kutta、Hermite數值軌道預報方法

4.1 基于J2攝動項的Runge-Kutta 數值軌道預報方法

為得到更高精度的衛星軌道預測結果，相比之下，數值方法能夠取得更為理想的結果，即間接引用泰勒展開，用積分區間上若干右函數值的線性組合確定相應的系數。一般來說，4 階Runge_Kutta就能夠取得較高的精度要求。現推導如下。

在J2000.0 坐標系下的笛卡爾坐標系中，設遞推量rn的表示方式為

各項分別表示衛星當前的三維位置和三維速度信息，其遞推公式為

其中h 為積分步長，ki（i=1，2，3，4）為過程變量，其表達式為

式（3）中的f函數輸出為

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

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

圖1 基于J2攝動項的Runge-Kutta數值位置預報

圖2 基于J2攝動項的Runge-Kutta數值速度預報

4.2 基于Hermite數值軌道預報方法

設t0+nh時刻（0

對n分別求一次導、二次導：

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

回代到式（7），則有

對式（12）中n求導，有

其中：

由此可遞推得到衛星軌道的位置、速度信息。基于Hermite數值衛星軌道預報如圖3、圖4所示。

圖3 基于Hermite數值位置預報

圖4 基于Hermite數值速度預報

5 軌道預報誤差分析

以STK 軟件的HPOP 算法為真值，對基于J2攝動項的Runge-Kutta 數值軌道預報誤差分別如圖5、圖6 所示，基于Hermite 數值軌道預報誤差分別如圖7、圖8所示。

圖5 基于J2攝動項的Runge-Kutta數值位置預報誤差

圖6 基于J2攝動項的Runge-Kutta數值速度預報誤差

圖7 基于Hermite數值預報位置誤差

圖8 基于Hermite數值預報速度誤差

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

6 結語

本文分析了WGS84 坐標與J2000.0 系下的坐標互換，并重點研究了基于J2攝動項的Runge-Kutta 數值預報、基于Hermite 數值預報的方法、精度等。一般來說，衛星軌道預報的數值計算方法的精度遠高于解析方法，而且避開了第一類奇點導致的預報精度惡化現象［4～5］，但精度要求不高及計算資源有限時，解析方法亦不失為一種分析衛星軌道特性的重要工具，比如計算軌道六根數的變化規律等，這也是下一步繼續努力的方向。