逆行地球同步軌道特性探索

劉 贊，高益軍

（1.北京控制工程研究所，北京100190；2.空間智能控制技術國家級重點實驗室，北京 100190）

逆行地球同步軌道特性探索

劉 贊1，2，高益軍1

（1.北京控制工程研究所，北京100190；2.空間智能控制技術國家級重點實驗室，北京 100190）

為了探索逆行地球同步軌道在主要環(huán)境力作用下的軌道特性，給出了適合進行數(shù)值仿真的軌道動力學模型，并以先進的RKDP方法進行求解.對所得仿真數(shù)據利用求和取平均的方法去除攝動力產生的短周期效應，通過分析去短周期項后的數(shù)據揭示出了逆行地球同步軌道的演變特點.

逆行地球同步軌道；軌道動力學；RKDP方法；軌道特性

逆行地球同步軌道（RGSO，retrograde geosynchronous orbit）特指軌道周期等于地球自轉周期、軌道傾角等于180°的軌道.由于它一天之內可以在地球同步軌道高度繞地球飛行兩圈，故可以考慮用它對一系列不同經度位置的重要靜止軌道衛(wèi)星執(zhí)行一些有價值的任務，這是靜止軌道（GEO，geostationary orbit）或類靜止軌道衛(wèi)星所難以勝任的.

然而，國內外很少見到有討論逆行地球同步軌道的文獻，本文旨在對其軌道特性進行初步的探索，希望對相關工作的開展有所啟發(fā).

當?shù)厍蛐l(wèi)星軌道傾角等于180°時，各類常見攝動運動方程均會出現(xiàn)不同程度的奇異性［1-2］.為此，首先給出軌道動力學模型，然后采用先進的數(shù)值方法（RKDP法）求解該動力學方程，最終通過對數(shù)值仿真解的處理和分析揭示出逆行地球同步軌道的演變特性.

1 軌道動力學模型

地球衛(wèi)星的軌道動力學模型由軌道動力學方程和環(huán)境力模型組成.

軌道動力學方程描述了軌道運動狀態(tài)與航天器的質心所受的外力之間的動力學關系.軌道運動狀態(tài)由軌道根數(shù)定義，而選擇不同的軌道根數(shù)，往往方程的具體形式、適用性及方便性也不盡相同.綜合適用性和方便性的考慮，本文選用航天器在J2000地心平赤道坐標系下位置和速度矢量的笛卡爾分量作為動力學方程所用軌道根數(shù).這個做法有三個好處：動力學方程形式簡單并適用于任意類型航天器軌道；利用現(xiàn)有算法可方便地實現(xiàn)笛卡爾根數(shù)與其他各類常見根數(shù)之間的轉換［2］；日月引力及太陽光壓可直接表達在J2000地心平赤道坐標系中，而地球引力在地固坐標系中的分量又可利用現(xiàn)有的算法準確地轉換到J2000地心平赤道坐標系下［1］.

1.1 軌道動力學方程

衛(wèi)星運動一般方程的矢量形式為

式中，r、F0和Fε分別表示衛(wèi)星的地心距矢量、地球的球形中心引力和各種軌道攝動力的合力.

對于逆行地球同步軌道有攝動力合力

式中，F(xiàn)εEG、FεMG、FεSG和 FεSRP分別表示地球非球形引力攝動、月球引力攝動、太陽引力攝動和太陽光壓力攝動產生的攝動力.

在J2000地心平赤道坐標系下，矢量方程（1）可以改寫為如下標量形式：

方程（4）～（6）右端第一項和第二項分別對應方程（1）右端第一項和第二項在 J2000地心平赤道坐標系下的分量.r表示地心距，而 x、y和 z是 r在 J2000地心平赤道坐標系下的分量.該方程組適用于對任意類型的軌道進行數(shù)值仿真計算.計算所得笛卡爾根數(shù)可用文獻［2］中4.9節(jié)的算法轉換成經典軌道根數(shù)的形式；相應的，以經典軌道根數(shù)給出的初值可用該文獻4.8節(jié)的算法轉換成笛卡爾根數(shù)的形式.

1.2 環(huán)境力

1.2.1 地球非球形引力

在地固坐標系下，由地球引力模型EGM96所給地球引力勢函數(shù)可導出地球非球形引力在地心距方向、地心緯度圈切線方向和地心經度圈切線方向產生的攝動力分量為

式中，Re表示地球赤道半徑，λ和φ分別表示航天器地心經度和地心緯度.μe、和分別為地球引力常數(shù)和歸一化地球引力系數(shù)，其值可參考地球引力模型EMG96（sinφ）為歸一化締合勒讓德多項式.ˉPnm（sinφ）及其對sinφ的導數(shù)的計算可采用遞推算法［1，4］.采用歸一化的量和參數(shù)是為了提高計算的精度.對地球同步軌道而言，至多只需考慮到12階12次的地球引力球諧項［5］.

受歲差、章動的影響，航天器在地固坐標系中的坐標（x′，y′，z′）與其在 J2000地心平赤道坐標系中的坐標（x，y，z）之間存在如下關系：

式中（HG）表示從J2000地心平赤道坐標系到地固坐標系的坐標變換矩陣，其具體形式見參考文獻［1］的 1.4節(jié).

易得航天器的地心經緯度為

地球非球形引力在地固坐標系下的分量矩陣（Fr，F(xiàn)λ，F(xiàn)φ）與 J2000地心平赤道坐標系下的分量矩陣（FεEGx，F(xiàn)εEGy，F(xiàn)εEGz）之間有如下轉換關系：

式中（HG）－1是（HG）的逆矩陣，坐標軸旋轉矩陣 Tz（－λ）和 Ty（－φ）形式如下：

1.2.2 日月引力引起的攝動力

太陽引力產生的攝動力FεSG及月球引力產生的攝動力 FεMG分別為［1］

式中，rs和rm分別為太陽和月亮相對地球的位置矢量，rs和rm是它們的模；k為 J2000地心平赤道坐標系Z軸方向的單位矢量，zm則是rm在該方向上的分量；μs和 μm分別為太陽和月球的引力常數(shù).Δm和Δs分別為航天器相對太陽和月球的位置矢量，Δm和Δs是它們的模.太陽位置矢量的計算可參考文獻［6］，而高精度月球位置矢量可由DE405星歷文件解算得來.由于各位置矢量可直接表達在J2000地心平赤道坐標系下，故不需對 FεSG和 FεMG的分量矩陣做任何變換.

1.2.3 太陽光壓產生的攝動力

太陽光壓產生的攝動力為［1］

式中，k為與衛(wèi)星反射率相關的系數(shù)，Ae為衛(wèi)星等效截面積，m為衛(wèi)星質量，Δ0是日地平均距離，ρSRP0是Δ0處太陽光壓強度.該攝動力也可直接表達在J2000地心平赤道坐標系下，無需進行坐標變換.

1.3 求 解

用特殊攝動法求解衛(wèi)星軌道運動方程比較實用的算法有Gauss Jackson方法（簡記為GS法）、RKF7（8）方法［7］、Bulirsch Stoer外推法［8］（簡記為 BS法）及 RKDP方法［7-9］.文獻［10］顯示圓軌道情況下（定步長），GS法速度最快，而 BS法的速度則慢許多，RKDP方法的速度介于兩者之間.GS法和BS法的計算精度均能達到10－8，而 RKDP方法則可達到10－6.大橢圓情況下（變步長），RKDP方法速度最快，BS法次之，而GS法則不太適應變步長的積分.此時，RDKP方法的計算精度仍能達到10－6，而 BS法只能達到10－5.

實際上，經過多年的發(fā)展，從某些軟件的表現(xiàn)來看即便是近圓軌道，變步長的BS法和RKF7（8）方法也要比定步長的GS法快，尤其是變步長的RKF7（8）方法.RKF7（8）方法的優(yōu)異表現(xiàn)證明了高階RK方法加入步長控制機制后在計算速度上有很強的競爭力.這印證了前人所做的分析［7］，RK方法階次越高穩(wěn)定域越快，而穩(wěn)定域更寬則意味著可以選取更大的步長.

本文采用先進的八階 RKDP方法［7，9］，它采用PI控制器進行步長控制［11-12］，既能提高運算速度又能克服剛性問題.此外，它還具有密集輸出［7，9］（dense output）能力，可以快速地輸出積分軌跡.它最初由Dormand和 Prince在 70年代提出［13］，經過進一步的研究工作［11-12，14-17］才形成現(xiàn)在的八階算法［7，9］.為了使誤差系數(shù)最小，Dormand和 Prince對算法中公式的系數(shù)進行了優(yōu)化并采用嵌入式RK公式對進行內部誤差估計［14-17］.本文利用該方法對二體軌道進行了長達十年的積分，全局誤差的量級為微米級，完全能滿足軌道特性分析的要求.

此外，仿真表明，對于RGSO而言，只需考慮到3階以下地球引力球諧項就可以使軌道預測精度達到10m級.這一結論對于簡化近似解析解的求取是非常有幫助的.

2 逆行地球同步軌道仿真分析

通過分析仿真數(shù)據，本節(jié)給出了逆行地球同步軌道各根數(shù)在三種主要攝動力作用下的演變規(guī)律，包括長期項和長周期項.長期項指非周期性的變化趨勢項，長周期項指周期比軌道周期長的周期性變化趨勢項.雖然偏置效應也是一種非周期性趨勢項，但由于它不隨時間變化，故單獨列出.

本節(jié)所給各圖的數(shù)據都經過了去短周期項處理，故反映的是長期項和長周期項效應.

2.1 半長軸的演變規(guī)律

地球非球形引力作用下，半長軸有厘米級的周期變化項和米級的負向偏置.月球引力作用下，半長軸有10 m級的周期變化項和百米級的正向偏置，周期項周期等于月球軌道周期.太陽引力作用下，半長軸有米級的周期變化項和百米級的負向偏置，周期項周期為一個季度.太陽光壓作用下，半長軸有分米級的周期變化項和10m級的正向偏置.

典型值見圖1.

圖1 半長軸的演變規(guī)律

2.2 偏心率的演變規(guī)律

地球非球形引力作用下，偏心率主要有10－5級的偏置.月球引力作用下，偏心率有10－5級的周期變化項，周期等于月球軌道周期.太陽引力作用下，偏心率有10－6級的周期變化項及10－5級的偏置.太陽光壓作用下，偏心率有10－4級的周期變化項及正向偏置，周期項周期為一年.日月引力及太陽光壓的攝動效應見圖2.

2.3 傾角的演變規(guī)律

地球非球形引力及太陽光壓對傾角的影響幾乎為零.月球引力作用下，傾角有10－3級的長期變化項和10－3級的周期變化項.太陽引力作用下，傾角有10－3級的長期變化項和10－2級的周期變化項，周期項周期為半年.

圖2 偏心率的演變規(guī)律

當?shù)厍蛲杰壍纼A角較大時，地球非球形引力對傾角的影響不容忽視.實際上，在它與日月引力的耦合作用下，地球同步軌道傾角變化沒有長期變化項，但考慮到RGSO衛(wèi)星的壽命期有限，尤其從它的定義出發(fā)，應該將周期為五十多年的長周期項視為長期變化項.典型值見圖3.

圖3 傾角的演變規(guī)律

2.4 升交點赤經的演變規(guī)律

地球非球形引力作用下，升交點赤經有10－2級的長期變化項.日月引力攝動對升交點赤經的影響與初始條件有關.然而，由升交點赤經和傾角定義的傾角矢量卻有著與初始條件“無關”的變化趨勢，且只有當軌道節(jié)線初始方向與該變化趨勢方向一致時才能獲得穩(wěn)定的升交點赤經.太陽光壓作用下，升交點赤經有10－5級的周期變化項.

圖4（b）中月球相關曲線的大斜率段對應圖3中月球相關曲線在180°附近的那段，這時升交點赤經快速變化并由一種平穩(wěn)狀態(tài)進入另一種平穩(wěn)狀態(tài).

圖4（b）表明，當初值為 270°時，升交點赤經將處于相對穩(wěn)定的狀態(tài).

圖4 升交點赤經的演變規(guī)律

2.5 近地點幅角的演變規(guī)律

地球非球形引力作用下，近地點幅角有10－2級的長期變化項.月球引力作用下，近地點幅角有10－2級的長期變化項和100級的周期項，周期項周期等于月球軌道周期.太陽引力作用下，近地點幅角有100級的偏置和100級的周期變化項，周期項周期為半年.太陽光壓作用下，近地點幅角有100級的偏置效應和101～102級的周期項，周期項的周期為一年.實際上，各攝動力對近地點幅角的周期項效應均與偏心率的初值有關，偏心率初值越小，周期項變化的幅值越大.

圖5（b）中月球相關曲線的大斜率段與圖4（b）中對應段相關，因為近地點幅角是從升交點起量的.

3 結 論

不考慮短周期效應，在主要環(huán)境力作用下，逆行地球同步軌道的演變過程有如下特點：半長軸受日月引力影響，正向偏置0.3 km后基本保持不變，即軌道周期基本恒定；偏心率受太陽光壓影響有一個10－4級的周年變化，受月球引力影響有一個10－5級的周月變化；傾角受日月引力影響以0.85（°）／年的平均速度增加或減小，變化方向與升交點赤經初值有關且當傾角等于180°時改變變化方向；升交點赤經主要受日月引力影響，同時在地球非球形引力作用下緩慢東進；近地點幅角受太陽光壓影響做周年變化，其幅值與偏心率初值有關，而變化區(qū)間所在象限則與偏心率矢量的初值有關.

圖5 近地點幅角的演變規(guī)律

［1］劉林.航天器軌道理論［M］.北京：國防工業(yè)出版社，2000

［2］Vladimir A C.Oribtal mechanics［M］.3rd ed.Reston：American Institute of Aeronautics and Astronautics Inc.，2002

［3］章仁為.靜止衛(wèi)星的軌道和姿態(tài)控制［M］.北京：科學出版社，1987

［4］張強，劉林.關于勒讓德多項式的算法問題［J］.紫金山天文臺臺刊，1996，15（4）：259-283

［5］David A V.An analysis of state vector propagation using differing flight dynamics programs［C］.The AAS／AIAA Space Flight Mechanics，Copper Mountain，Colorado，USA，Jan 2005

［6］郗曉寧，王威.近地航天器軌道基礎［M］.長沙：國防科技大學出版社，2003

［7］Hairer E，Norsett S P，Wanner G.Solving ordinary differential equations I：non－stiff problems［M］.2nd ed.New York：Springer，1993

［8］Stoer J，Bulirsch R.Introduction to numerical analysis［M］.3rd ed.New York：Springer，2002

［9］Willam H P，Saul A T，William T V，Brian P F.Numerical recipes：the art of scientific computing edition［M］.3rd ed.New York：Cambridge University Press，2007

［10］Ken F.Numerical integration of the equations of motion of celestial mechanics［J］.Celestial Mechanics and Dynamical Astronomy，1984，33（2）：127-142

［11］Gustafsson K.Control theoretic techniques for stepsize selection in explicit Runge-Kutta methods［J］.ACM Transactions on Mathematical Software，1991，17（4）：533-554

［12］Soderlind G.Digital filters in adaptive time-stepping［J］.ACM Transactions on Mathematical Software，2003，29（1）：1-26

［13］Dormand J R，Prince P J.New Runge-Kutta-Nystrom algorithms for simulation in dynamical astronomy［J］.Celestial Mechanics and Dynamical Astronomy，1978，18（3）：223-232

［14］Dormand JR，Prince P J.A family of embedded Runge-Kutta formulae［J］.Journal of Computational and Applied Mathematics，1980，6（1）：19-26

［15］Dormand J R，Prince P J.Runge-Kutta triples［J］.Computers＆ Mathematics with Applications，1986，12A（9）：1007-1017

［16］Dormand J R，Prince P J.Runge-Kutta-Nystrom triples［J］.Computers＆ Mathematics with Applications，1987，13（12）：937-949

［17］Dormand J R，Prince P J.Practical Runge-Kutta processes［J］.SIAM Journal on Scientific and Statistical Computing，1989，10（5）：977-989

Study of Characteristics of Retrograde Geosynchronous Orbit

LIU Zan1，2，GAO Yijun1
（1.Beijing Institute of Control Engineering，Beijing 100190，China；2.National Laboratory of Space Intelligent Control，Beijing 100190，China）

To study evolving characteristics of a retrograde geosynchronous orbit influenced by major environmental forces，an orbit dynamics model suitable for numerical simulations is given and solved by the advanced RKDP method.Numerical solutions are periodically summed and averaged to get rid of shortperiodic-term effects.By analyzing polished data，its orbit characteristics are revealed to us.

retrograde geosynchronous orbit； orbital dynamics；RKDP method；orbit characteristics

V448

A

1674-1579（2009）04-0052-05

2008-11-04

劉 贊（1983—），男，湖南人，碩士研究生，研究方向為軌道動力學與控制（e-mail：liuzan＠yahoo.cn）.