宋傳峰，秘金鐘，黨亞民，薛樹強

（中國測繪科學研究院，北京 100830）

1 引言

隨著GEO （Geostationary Earth orbit）與IGSO （Inclined Geo Synchronous orbit）衛(wèi)星在導航定位領域中的應用，對其軌道的精度要求越來越高［1］。利用數(shù)值算法積分衛(wèi)星軌道是人造衛(wèi)星軌道預報與精密定軌中的基礎工作［2］。由于二體問題可以直接得出衛(wèi)星位置、速度和加速度的解析解，因此本文避開復雜的攝動力計算，在二體意義下分別對GEO和IGSO衛(wèi)星進行了軌道數(shù)值積分，并與其理論值進行了分析比較，得出了有益結論。若考慮攝動力，則不便于結果比較。

目前，在MEO定軌中使用的積分方法有RKF與Admas預測－校正系統(tǒng)結合的方法 （如Gamit軟件），Collocation方法 （如Bernese軟件）等。本文驗證了RKF與Admas預測－校正系統(tǒng)結合的方法在GEO與IGSO衛(wèi)星軌道積分時的適用性。

2 數(shù)值積分方法

2.1 變步長RKF7（8）階單步積分法

Runge－Kutta－fehlberg方法 （RKF）是一種嵌套的RK方法，即同時給出n階和n＋1階兩組RK計算公式，用兩組公式計算出來的xi＋1之差來估計截斷誤差，根據(jù)截斷誤差的大小來控制步長［3］。以下給出7階和8階的RKF公式

截斷誤差

式中，h為積分步長，ck，βkλ為積分算法的系數(shù)。

2.2 定步長Adams線性多步積分法

在衛(wèi)星運動方程的數(shù)值積分中，單步法只是用于多步法的起步算法，當采用單步法推出足夠的步點后，就可采用高精度的多步法往前推算［3］。本文使用變步長的RKF7（8）階積分法積分10步后，再使用定步長的10階Adams預測－校正系統(tǒng)進行積分。以下給出顯式的Adams－Bashforth公式和隱式的Adams－Moulton公式。

顯式

隱式

3 GEO與IGSO衛(wèi)星二體問題

3.1 GEO衛(wèi)星二體問題

在二體意義下，GEO地球同步靜止軌道衛(wèi)星在地心地固坐標系中的坐標為一固定值，但軌道積分是在慣性坐標系下進行的，因此需要計算二體問題下地固坐標系與慣性坐標系間的相互轉換矩陣。

假設積分初始時刻t0（s），兩坐標系坐標軸完全重合，則任意t（s）時刻，地固系向慣性系轉換的矩陣為

在慣性系下，GEO衛(wèi)星繞Z軸旋轉的角速度與地球自轉角速度相等，因此易得下式

慣性系下GEO衛(wèi)星在X－Y平面內運動，即Z≡0，由此可得：

式中，X、Y、Z表示衛(wèi)星在慣性系下的坐標。

在二體意義下，GEO衛(wèi)星的軌跡為圓，且其在地固坐標系下的速度為0，因此給定初始坐標可得任意時刻GEO衛(wèi)星在地固和慣性坐標系下的位置、速度和加速度。

3.2 IGSO衛(wèi)星二體問題

IGSO衛(wèi)星軌道是特殊的地球同步圓軌道，軌道高度跟GEO一樣，均屬于地球同步軌道［4］。由于IGSO衛(wèi)星軌道在慣性系下也是圓軌道，因此可看做由GEO衛(wèi)星軌道繞X或Y軸或X－Y平面內一條直線旋轉i角，i為IGSO衛(wèi)星軌道傾角，本文中按繞X軸旋轉i角計算。

若t0時刻，IGSO衛(wèi)星在其軌道直角坐標系（慣性系繞X軸旋轉i角）下坐標為（，，0）T，速度為（－ωe，ωe，0）T，則t時刻其在軌道直角坐標系中坐標為（X＊，Y＊，0）T＝Rz（α），，0）T，速度為（－ωeY＊，ωeX＊，0）T。ωe，α與3.1節(jié)的意義相同。

t時刻IGSO衛(wèi)星在地心慣性坐標系下位置和速度：

式中，X、Y、Z、VX、VY、VZ為衛(wèi)星在慣性坐標系下位置和速度，i表示軌道傾角，X＊、Y＊為衛(wèi)星在軌道直角坐標系下平面坐標。

將慣性系下位置左乘矩陣RT，R同式 （6），即可得地固系下衛(wèi)星位置，再轉成大地坐標系 （B，L，H）即可得IGSO衛(wèi)星的星下點軌跡，圖1為經度為45°，傾角為10°和30°的IGSO衛(wèi)星星下點軌跡 （二體意義下）。

4 算例分析與比較

對不同經度的GEO衛(wèi)星、不同經度和傾角的IGSO衛(wèi)星，先使用變步長的RKF7（8）階積分法積分10步，截斷誤差設定為10－12，再使用定步長75s的10階Adams預測－校正系統(tǒng)進行積分48h（2圈），每隔15min輸出一組結果。不同經度的GEO衛(wèi)星在地心地固坐標系下的積分初值見表1，其速度均為0；不同經度和傾角的IGSO衛(wèi)星在慣性坐標系下的積分初值見表2。

圖1 IGSO衛(wèi)星星下點軌跡 （二體意義）

表1 GEO衛(wèi)星積分初值

表2 IGSO衛(wèi)星積分初值

由于篇幅原因，僅給出了Y方向差異的統(tǒng)計信息和經度為30°的GEO衛(wèi)星差異變化圖。圖2是經度為30°的GEO衛(wèi)星軌道數(shù)值積分結果與二體意義下結果在地心地固坐標系下的差異，歷元間隔為15分鐘，表3是GEO衛(wèi)星軌道積分結果與真值 （二體意義）在地固系下Y方向差異的統(tǒng)計信息。

圖2 GEO衛(wèi)星軌道數(shù)值積分差異

表3 GEO衛(wèi)星軌道數(shù)值積分Y方向差異統(tǒng)計信息

由于篇幅原因，僅給出了X方向差異的統(tǒng)計信息和經度為45°傾角為30°的IGSO衛(wèi)星差異變化圖。圖3是經度為45°傾角為30°的IGSO衛(wèi)星軌道數(shù)值積分結果與二體意義下結果在地心地固坐標系下的差異，歷元間隔為15min，表4是IGSO衛(wèi)星軌道積分結果與真值 （二體意義）在地固系下X方向差異的統(tǒng)計信息。

表4 IGSO衛(wèi)星軌道數(shù)值積分X方向差異統(tǒng)計信息

圖3 IGSO衛(wèi)星軌道數(shù)值積分差異

5 結論

由以上圖表可以得出，二體意義下GEO、IGSO衛(wèi)星軌道數(shù)值積分結果與真值相差甚小，均在亞毫米級或以內，與現(xiàn)有精密定軌精度相比可以忽略，適當減小積分步長可獲得更小的理論值差異，但這將明顯增加計算耗時。結果表明對MEO（如GPS）衛(wèi)星適用的軌道數(shù)值積分方法對GEO和IGSO衛(wèi)星也同樣適用。

本文僅在二體意義下驗證了數(shù)值積分方法的適用性，現(xiàn)實中導航衛(wèi)星攝動力復雜多變，因此建立合適的力學模型再結合數(shù)值積分方法進行軌道積分是統(tǒng)計定軌中的重要步驟。

［1］劉吉華，歐吉坤，鐘世民，等.GEO與IGSO衛(wèi)星星間差分的精密定軌［C］//中國衛(wèi)星導航學術年會組委會.第一屆中國衛(wèi)星導航學術年會文集.北京：出版社不詳，2010：1－7.

［2］李得海，袁運斌，歐吉坤，等.12階Runge－kutta2次算法的衛(wèi)星軌道積分研究［J］.武漢大學學報：信息科學版，2010，35（11）：1335－1338.

［3］趙齊樂.GPS導航星座及低軌衛(wèi)星的精密定軌理論和軟件研究［D］.武漢：武漢大學，2004.

［4］楊 洋，范 麗，董緒榮.IGSO星座分析與優(yōu)化設計［J］.兵工自動化，2008，27（11）：53－55.

［5］ES－HAGH M.Step Variable Numerical Orbit Integration of a Low Earth Orbiting Satellite［J］.Journal of the Earth ＆Space Physics，2005，31（1）：1－12.

［6］XU Guo－chang.GPS Theory Algorithms and Applications［M］.Berlin：Springer－Verlag，2003.

［7］IERS Convention Centre.IERS Conventions：2010［M］.Frankfurt am Main：Verlag des Bundesamts für Kartographie und Geod?sie，2010.

［8］LEICK A.GPS Satellite Surveying［M］.Hoboken：John Wiley ＆Sons Inc.，2004.