拉格朗日/高斯無奇點衛(wèi)星運動方程推導(dǎo)與分析

蔣春華，徐天河，喬 晶，杜玉軍，王慶，許國昌

1. 山東大學(xué)空間科學(xué)研究院，山東 威海 264209； 2. 西安測繪研究所地理信息工程國家重點實驗室，陜西 西安 710054； 3. 東南大學(xué)儀器科學(xué)與工程學(xué)院，江蘇 南京 210096； 4. 香港理工大學(xué)，香港 999077

當(dāng)衛(wèi)星按圓軌道、赤道軌道以及圓赤道軌道運動時，某些衛(wèi)星開普勒(Kepler)軌道根數(shù)發(fā)生退化，用于表達衛(wèi)星軌道的拉格朗日運動方程和源于拉格朗日運動方程的高斯運動方程發(fā)生數(shù)學(xué)奇異，即為奇點問題，它一直是天體力學(xué)和N體問題研究中的核心問題之一[1-4]。衛(wèi)星軌道理論研究主要基于二體問題的攝動理論[5-7]，二階攝動理論解研究近年來取得了長足進步[8-12]，為奇點問題的進一步研究奠定了基礎(chǔ)。眾多學(xué)者在理論研究方面做了很多努力[13-16]，也取得了不少成果。第1類和第2類奇點問題為解決衛(wèi)星導(dǎo)航領(lǐng)域的奇點問題提供了可能，但其表達復(fù)雜且解決思路和最終表達很難統(tǒng)一[17-18]。文獻[19—20]提出參量代換和坐標系旋轉(zhuǎn)等方法，分別從圓軌道，赤道軌道和極軌道的情況，對其方法進行論證和闡述，為奇點問題的解決提供了新的思路。但其忽略了代換參量實際的物理和幾何意義，以及轉(zhuǎn)換過程對整個軌道系統(tǒng)求解的可逆性。2012年，文獻[21]以不定積分解的形式提出無奇點理論，并于2013年有了以微分方程描述的公式版本[22]。但其公式推導(dǎo)主要是邏輯推理，2014年給出了拉格朗日無奇點運動方程純數(shù)學(xué)的嚴格推導(dǎo)與證明[23]。2015年又給出了高斯無奇點運動方程純數(shù)學(xué)的嚴格推導(dǎo)與證明[24]。上述論文主要是從純數(shù)學(xué)角度進行推導(dǎo)證明，參數(shù)的物理意義和方程的連續(xù)性有待進一步分析。

為了解決上述問題，本文充分考慮參數(shù)的實際幾何物理意義、總體方程的統(tǒng)一性以及與原始方程的一致性，從拉格朗日無奇點運動方程出發(fā)，對奇點問題作了進一步的分析和探索。奇點情況時，某些Kepler軌道根數(shù)發(fā)生退化，本文基于原始拉格朗日和高斯運動方程及其物理意義，對衛(wèi)星無奇點方程的判據(jù)進行推導(dǎo)與分類，并分別對拉格朗日/高斯無奇點衛(wèi)星運動方程在圓軌道、赤道軌道和圓赤道軌道情況下存在奇異的軌道參數(shù)(升交點赤經(jīng)、近地點角距和平近點角)賦予新的實際物理和幾何意義。在拉格朗日-許與高斯-許衛(wèi)星運動方程的基礎(chǔ)上，根據(jù)參數(shù)定義對奇異項進行加和變換，從而消除奇異因子，最后給出一種新的拉格朗日/高斯無奇點衛(wèi)星運動方程的統(tǒng)一表達式。

1 衛(wèi)星無奇點方程的判據(jù)

衛(wèi)星運動方程的奇點問題源于某些特定情況時Kepler根數(shù)的物理意義的模糊，導(dǎo)致拉格朗日運動方程出現(xiàn)小分母或分母為零的問題。拉格朗日衛(wèi)星運動方程可表示為[7,14]

(1)

式中，函數(shù)R稱為攝動位函數(shù)，包括了除中心引力位以外的全部項，6個開普勒根數(shù)a、e、ω、i、Ω、M，即半長軸、偏心率、近地點角距、傾角、升交點赤經(jīng)、平近點角。文獻[21]將其改寫成

(2)

式中，δ及其下標表示式(1)的右邊各項。當(dāng)e或者sini趨近于0時，即可能發(fā)生奇點問題。奇點問題解決的前提是需要對奇點域做具體的判斷。首先，定義徑向和赤道方向的誤差如下[22]

(3)

此處徑向誤差ae是軌道幾何中心O′與地心O的距離，如圖1所示，也是e引起的最大徑向誤差。赤道偏差asini≈ai表示衛(wèi)星由i引起的垂直于赤道的最大距離(傾角i是衛(wèi)星的最大緯度)。在積分時刻t0=0，e或者/和i接近于0，則在積分時刻t，偏差e和i就為Δe和Δi。式(3)可寫成

(4)

圖1 衛(wèi)星的扁心率和近點角Fig.1 The eccentricity and anomalies of a satellite

圓軌道、赤道軌道以及圓赤道軌道的軌道偏差定義為

(5)

式中，ε0表示在衛(wèi)星運動方向上由ΔΩ、Δω、ΔM引起的最大偏差。徑向偏差、赤道方向偏差和軌道方向偏差是軌道最大誤差的3個組成部分。3個偏差表示(e,ω,i,Ω,M)擾動對軌道產(chǎn)生的最大誤差域。實際應(yīng)用中，假定軌道3個方向偏差εe、εi和ε0分別為1 m，即這里定義的εe、εi和ε0可以作為判據(jù)標準。對于不同軌道，判據(jù)條件如表1所示，用于對奇點問題的判斷與進一步分析。

表1 奇點情況的判斷條件

2 拉格朗日無奇點衛(wèi)星運動方程推導(dǎo)分析

拉格朗日無奇點運動方程的數(shù)學(xué)推導(dǎo)已在文獻[21]中詳細給出。對于通常非奇點情況拉格朗日方程[15,22]可寫為

(6)

對于奇點情況還需進一步分析。由軌道參數(shù)的原始定義有：Ω表示在地球赤道平面上，升交點與春分點之間的地心夾角稱為升交點赤經(jīng)；ω表示在軌道平面上近地點與升交點之間的地心角距稱為近地點角距；f表示衛(wèi)星到近地點之間的地心角距，如圖2所示。為計算真近點角，二體問題中引入偏近點角和平近點角的定義，3種近點角(真近點角f、平近點角M、偏近點角E)均是時間的函數(shù)，用來表述衛(wèi)星在ECI坐標系中的位置隨著時間的變化。其中，平近點角M采用平均角速度n來描述衛(wèi)星繞地球的軌道運動，廣泛用于解析軌道的理論推導(dǎo)和應(yīng)用，本文的推導(dǎo)也是基于此。

圖2 軌道幾何圖Fig.2 Orbital geometry

奇點情況分為圓軌道、赤道軌道和圓赤道軌道3種情況。當(dāng)為圓軌道時，e=0近地點失去意義，無法確定。此時開普勒軌道根數(shù)發(fā)生退化，平近點角M和近地點角距ω意義發(fā)生改變。因此，從衛(wèi)星實際運動情況分析，可定義新的平近點M′表示衛(wèi)星與升交點的夾角，即

M′=M+ω

(7)

當(dāng)為赤道軌道時，sini=0，升交點失去意義，無法確定。此時開普勒軌道根數(shù)發(fā)生退化，升交點赤經(jīng)Ω和近地點角距ω意義發(fā)生改變。因此，可定義新的近地點角距ω′表示近地點與春分點的夾角，即

ω′=Ω+ω

(8)

當(dāng)為圓赤道軌道時，e=0且sini=0，近地點和升交點均失去意義，無法確定。此時開普勒軌道根數(shù)發(fā)生退化，平近點角M和近地點角距ω以及升交點赤經(jīng)Ω意義均發(fā)生變化。因此，定義一個新的平近點M″表示衛(wèi)星與春分點的夾角，即

M″=M+ω+Ω

(9)

2.1 圓軌道情況

圓軌道情況ae<εe=1、asini≥εi=1時，從物理意義分析，近地點無法確定，相關(guān)軌道參數(shù)對時間的偏導(dǎo)數(shù)也無法確定。從數(shù)學(xué)角度分析e除到式(2)右邊，分母存在零因子，即發(fā)生奇點。由式(7)，可得

(10)

因此，圓軌道奇點情況的運動方程由式(6)中一、二、三、四式，以及下式組成

(11)

2.2 赤道軌道情況

赤道軌道情況ae≥εe=1、asini<εi=1時，從物理意義分析，升交點無法確定，相關(guān)軌道參數(shù)對時間的偏導(dǎo)數(shù)也無法確定。從數(shù)學(xué)角度分析，sini除到式(2)右邊，分母存在零因子，即發(fā)生奇點。由式(8)可得

(12)

因此，赤道軌道奇點情況的運動方程由式(6)中一、二、三、六式，以及下式組成

(13)

2.3 圓赤道軌道情況

圓赤道軌道情況ae<εe=1、asini<εi=1時，從物理意義分析，升交點和近地點無法確定，相關(guān)軌道參數(shù)對時間的偏導(dǎo)數(shù)也無法確定。從數(shù)學(xué)角度分析，sini、e除到式(2)右邊，分母存在零因子，即發(fā)生奇點。由式(9)，可得

(14)

因此，圓赤道軌道奇點情況的運動方程由式(6)中一、二、三式，以及下式組成

(15)

3 高斯無奇點衛(wèi)星運動方程推導(dǎo)分析

非奇點情況下高斯運動方程為[24]

(16)

式中，fa、fh、fr分別表示軌道坐標系的3個方向的攝動力。對于奇點情況下的高斯運動方程做進一步推導(dǎo)。

3.1 圓軌道情況

ae<εe=1、asini≥εi=1，根據(jù)拉格朗日圓軌道運動方程，得

(17)

因此，圓軌道奇點情況的運動方程由式(16)中一、二、三、五式，以及下式組成

(18)

3.2 赤道軌道情況

ae≥εe=1、asini<εi=1，根據(jù)拉格朗日赤道軌道運動方程，得

(19)

因此，赤道軌道奇點情況的運動方程由式(16)中一、二、三、六式，以及下式組成

(20)

3.3 圓赤道軌道情況

ae<εe=1,asini<εi=1，根據(jù)拉格朗日圓赤道軌道運動方程，得

(21)

因此，圓赤道軌道奇點情況的運動方程由式(16)中一、二、三式，以及式(22)組成

(22)

4 無奇點衛(wèi)星運動方程的連續(xù)性分析

總結(jié)上述非奇點與奇點的各種情況，最終統(tǒng)一的拉格朗日無奇點衛(wèi)星運動方程為

(23)

最終統(tǒng)一的高斯無奇點衛(wèi)星運動方程為

(24)

衛(wèi)星的運動和力模型是連續(xù)的，進一步分析衛(wèi)星運動方程的連續(xù)性可知，式(23)、(24)對于任何一種奇點情況(圓軌道奇點、赤道軌道奇點、圓赤道軌道奇點)方程是非間斷的，奇點情況下運動方程中重新定義的組合量及非奇點情況下運動方程的軌道參數(shù)在其對應(yīng)的判據(jù)域內(nèi)均是連續(xù)的。因此，新的拉格朗日和高斯無奇點衛(wèi)星運動方程具有連續(xù)性。

5 仿真與分析

(1) 圓軌道情況。軌道坐標系下狀態(tài)矢量進行一次旋轉(zhuǎn)，分析角度物理意義為

(25)

(2) 赤道軌道情況。軌道坐標系下旋轉(zhuǎn)矩陣，分析角度物理意義可知

(26)

式中，R1、R3為x和z中的旋轉(zhuǎn)矩陣[16]，圓赤道軌道兩點兼有。

下面利用本文推導(dǎo)的拉格朗日無奇點衛(wèi)星運動方程描述衛(wèi)星運動，軌道初值如表2所示(其中后3個可表示新定義的根數(shù))。首先基于式(23)求出各時刻的衛(wèi)星的軌道根數(shù)，然后利用式(25)、(26)將其轉(zhuǎn)換成衛(wèi)星狀態(tài)矢量，最后利用自主編寫的衛(wèi)星軌道三維可視化仿真軟件分別繪制二體運動圓軌道、赤道軌道以及圓赤道軌道軌跡[25]，如圖3所示。將計算軌道與積分軌道對比，軌道差在厘米級。高斯方程同理。

表2 圓軌道、赤道軌道及圓赤道軌道的軌道初值

圖3 圓軌道、赤道軌道以及圓赤道軌道衛(wèi)星軌跡Fig.3 The circular/equatorial/circular and equatorial satellite orbit

從上述過程與結(jié)果可知，本文推導(dǎo)的無奇點衛(wèi)星運動方程，可描述二體運動中的圓軌道、赤道軌道以及圓赤道軌道衛(wèi)星運動，且計算精度與非奇點二體運動相當(dāng)，可用于解決二體運動中的奇點問題，從而在一定程度上證明本文推導(dǎo)公式的正確性和可用性。

6 結(jié) 論

本文推導(dǎo)了衛(wèi)星拉格朗日/高斯無奇點運動方程，無任何近似，與原始拉格朗日/高斯運動方程具有邏輯一致性；奇點情況下，對參數(shù)含義重新定義，所得方程在奇點情況下仍具有幾何和物理意義；圓軌道、赤道軌道和圓赤道軌道下的衛(wèi)星運動方程的表達更加簡潔，便于理解和計算；軌道方程右端完全消除了零因子，從根本上解決衛(wèi)星運動方程的奇點問題，并且該衛(wèi)星運動方程具有連續(xù)性。

參考文獻：

[1] XIA Zhihong. The Existence of Noncollision Singularities in Newtonian Systems[J]. Annals of Mathematics, 1992, 135(3): 411-468.

[2] DIACU F. The Solution of The N-body Problem[J]. The Mathematical Intelligencer, 1996, 18(3): 66-70.

[3] DIACU F, HOLMES P. Celestial Encounters: the Origins of Chaos and Stability[M]. Princeton, NJ: Princeton University Press, 1996.

[4] WANG Qiudong. The Global Solution of the N-body Problem[J]. Celestial Mechanics and Dynamical Astronomy, 1990, 50(1): 73-88.

[5] XU Guochang. Sciences of Geodesy-I: Advances and Future Directions[M]. Berlin: Springer, 2010: 105-154.

[6] BROUWER D, CLEMENCE G M. Methods of Celestial Mechanics[M]. Burlington, MA: Elsevier, 2013.

[7] CHOBOTOV V A. Orbital Mechanics[M]. 3rd ed. Washington DC: AIAA, 2002.

[8] XU Guochang, XU Tianhe, YEH T K, et al. Analytical Solution of a Satellite Orbit Disturbed by Lunar and Solar Gravitation[J]. Monthly Notices of the Royal Astronomical Society, 2011, 410(1): 645-653.

[9] XU Yan, YANG Yuanxi, ZHANG Qin, et al. Solar Oblateness and Mercury’s Perihelion Precession[J]. Monthly Notices of the Royal Astronomical Society, 2011, 415(4): 3335-3343.

[10] XU Guochang, XU Jia. On Orbital Disturbing Effects of the Solar Radiation[J]. Monthly Notices of the Royal Astronomical Society, 2013, 432(1): 584-588.

[11] 張兵兵, 聶琳娟, 吳湯婷, 等. SWARM衛(wèi)星簡化動力學(xué)厘米級精密定軌[J]. 測繪學(xué)報, 2016, 45(11): 1278-1284. DOI: 10.11947/j.AGCS.2016.20160284.

ZHANG Bingbing, NIE Linjuan, WU Tangting, et al. Centimeter Precise Orbit Determination for SWARM Satellite via Reduced-dynamic Method[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(11): 1278-1284. DOI: 10.11947/j.AGCS.2016.20160284.

[12] 鄒賢才, 李建成, 姜衛(wèi)平, 等. 衛(wèi)星重力資料分析的同解法研究及其仿真[J]. 測繪學(xué)報, 2010, 39(4): 344-348.

ZOU Xiancai, LI Jiancheng, JIANG Weiping, et al. Research on the Simultaneous Solution Method for Satellite Gravity Data Analysis and Its Simulation[J]. Acta Geodaetica et Cartographica Sinica, 2010, 39(4): 344-348.

[13] HAVEL K. N-body Gravitational Problem: Unrestricted Solution[M]. Brampton, ON, Canada: Grevyt Press, 2008.

[14] BATTIN R H. An Introduction to the Mathematics and Methods of Astrodynamics[M]. Reston, VA: AIAA, 1999.

[15] KAULA W M. Theory of Satellite Geodesy: Applications of Satellites to Geodesy[M]. Mineola, NY: Dover Publications Inc, 2000.

[16] XU Guochang. Orbits[M]. Berlin: Springer, 2008.

[17] 韓星遠, 向開恒, 王海紅. 第一類無奇點變量的廣播星歷參數(shù)擬合算法[J]. 航天器工程, 2011, 20(4): 54-59.

HAN Xingyuan, XIANG Kaiheng, WANG Haihong. Research on Broadcast Ephemeris Parameters Fitting Algorithm Based on the First Class of No Singularity Variables[J]. Spacecraft Engineering, 2011, 20(4): 54-59.

[18] 張中凱, 杜蘭, 旦增曲英, 等. 基于第二類無奇點根數(shù)的改進根數(shù)[J]. 測繪科學(xué)技術(shù)學(xué)報, 2012, 29(4): 257-261.

ZHANG Zhongkai, DU Lan, DAN Zengquying, et al. Improved Elements Based on Second Class of No-singularity Variables[J]. Journal of Geomatics Science and Technology, 2012, 29(4): 257-261.

[19] BROUCKE R A, CEFOLA P J. On the Equinoctial Orbit Elements[J]. Celestial Mechanics, 1972, 5(3): 303-310.

[20] BATTIN R H. An Introduction to the Mathematics and Methods of Astrodynamics[M]. New York: American Institute of Aeronautics and Astronautics. 1987.

[21] XU Guochang, XU Jia. On the Singularity Problem in Orbital Mechanics[J]. Monthly Notices of the Royal Astronomical Society, 2013, 429(2):1139-1148.

[22] XU Guochang, XU Jia. Orbits: 2nd Order Singularity-free Solutions[M]. Berlin: Springer, 2013.

[23] XU G, LV Z P, SHEN Y Z, et al. A Mathematical Derivation of Singularity-free Lagrange Equations of Planetary Motion, Special Issue for Celebration 80th Birthday of Academician Houze Xu[J]. Journal of Surveying and Mapping, 2014.

[24] 許國昌, 陳武, 沈云中, 等. 高斯無奇點衛(wèi)星運動方程的數(shù)學(xué)推導(dǎo)——謹以本文恭賀師兄歐吉坤教授七十壽辰[J]. 導(dǎo)航定位學(xué)報, 2015, 3(3): 5-12.

XU Guochang, CHEN Wu, SHEN Yunzhong, et al. A Mathematical Derivation of Singularity-free Gaussian Equations of Planetary Motion[J]. Journal of Navigation and Positioning, 2015, 3(3): 5-12.

[25] 杜玉軍. 衛(wèi)星軌道的三維可視化程序設(shè)計[D]. 武漢: 武漢大學(xué), 2008.

DU Yujun. 3D Visualization Programming of Satellite Orbiting[D]. Wuhan: Wuhan University, 2008.