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

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

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

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

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

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

(1)

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

(2)

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

(3)

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

(4)

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

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

(5)

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

表1 奇點(diǎn)情況的判斷條件

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

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

(6)

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

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

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

M′=M+ω

(7)

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

ω′=Ω+ω

(8)

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

M″=M+ω+Ω

(9)

2.1 圓軌道情況

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

(10)

因此，圓軌道奇點(diǎn)情況的運(yùn)動(dòng)方程由式(6)中一、二、三、四式，以及下式組成

(11)

2.2 赤道軌道情況

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

(12)

因此，赤道軌道奇點(diǎn)情況的運(yùn)動(dòng)方程由式(6)中一、二、三、六式，以及下式組成

(13)

2.3 圓赤道軌道情況

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

(14)

因此，圓赤道軌道奇點(diǎn)情況的運(yùn)動(dòng)方程由式(6)中一、二、三式，以及下式組成

(15)

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

非奇點(diǎn)情況下高斯運(yùn)動(dòng)方程為[24]

(16)

式中，fa、fh、fr分別表示軌道坐標(biāo)系的3個(gè)方向的攝動(dòng)力。對(duì)于奇點(diǎn)情況下的高斯運(yùn)動(dòng)方程做進(jìn)一步推導(dǎo)。

3.1 圓軌道情況

ae<εe=1、asini≥εi=1，根據(jù)拉格朗日?qǐng)A軌道運(yùn)動(dòng)方程，得

(17)

因此，圓軌道奇點(diǎn)情況的運(yùn)動(dòng)方程由式(16)中一、二、三、五式，以及下式組成

(18)

3.2 赤道軌道情況

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

(19)

因此，赤道軌道奇點(diǎn)情況的運(yùn)動(dòng)方程由式(16)中一、二、三、六式，以及下式組成

(20)

3.3 圓赤道軌道情況

ae<εe=1,asini<εi=1，根據(jù)拉格朗日?qǐng)A赤道軌道運(yùn)動(dòng)方程，得

(21)

因此，圓赤道軌道奇點(diǎn)情況的運(yùn)動(dòng)方程由式(16)中一、二、三式，以及式(22)組成

(22)

4 無奇點(diǎn)衛(wèi)星運(yùn)動(dòng)方程的連續(xù)性分析

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

(23)

最終統(tǒng)一的高斯無奇點(diǎn)衛(wèi)星運(yùn)動(dòng)方程為

(24)

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

5 仿真與分析

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

(25)

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

(26)

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

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

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

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

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

6 結(jié) 論

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

參考文獻(xiàn)：

[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)星簡(jiǎn)化動(dòng)力學(xué)厘米級(jí)精密定軌[J]. 測(cè)繪學(xué)報(bào), 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]. 測(cè)繪學(xué)報(bào), 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] 韓星遠(yuǎn), 向開恒, 王海紅. 第一類無奇點(diǎn)變量的廣播星歷參數(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] 張中凱, 杜蘭, 旦增曲英, 等. 基于第二類無奇點(diǎn)根數(shù)的改進(jìn)根數(shù)[J]. 測(cè)繪科學(xué)技術(shù)學(xué)報(bào), 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] 許國(guó)昌, 陳武, 沈云中, 等. 高斯無奇點(diǎn)衛(wèi)星運(yùn)動(dòng)方程的數(shù)學(xué)推導(dǎo)——謹(jǐn)以本文恭賀師兄歐吉坤教授七十壽辰[J]. 導(dǎo)航定位學(xué)報(bào), 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è)計(jì)[D]. 武漢: 武漢大學(xué), 2008.

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