劉 林湯靖師侯錫云

(1南京大學(xué)天文與空間科學(xué)學(xué)院南京210046)

(2南京大學(xué)空間環(huán)境與航天動(dòng)力學(xué)研究所南京210093)

1 引言

關(guān)于定點(diǎn)在地-月系共線平動(dòng)點(diǎn)附近的探測器軌道外推(或稱預(yù)報(bào))問題,并不是一個(gè)新問題,該類型的飛行器,實(shí)際上就是一顆帶有一定軌道特征的遠(yuǎn)地衛(wèi)星.但是,盡管理論上繞地運(yùn)行周期與月球一致,但它受到的月球引力影響卻不是一個(gè)小擾動(dòng),而是一個(gè)幾乎與地球引力作用同等重要的外力源.因此,該問題也無法處理成一個(gè)簡單的受攝二體問題.如果仍然像構(gòu)造平動(dòng)點(diǎn)特殊軌道(如暈軌道)問題那樣,在地-月系中來處理該問題,也無法像衛(wèi)星運(yùn)動(dòng)那樣處理成一個(gè)簡單的受攝二體問題,其難點(diǎn)有如下兩個(gè)方面:

(1)即使基本模型采用簡單的圓型限制性三體問題[1?3],其基本解也較復(fù)雜.如果考慮月球軌道偏心率,處理成橢圓型限制性三體問題,不僅基本解更復(fù)雜,而且無濟(jì)于事.因?yàn)樘栆z動(dòng)影響與月球軌道偏心率的效應(yīng)相當(dāng),無論采用上述哪類基本模型,要在此基礎(chǔ)上構(gòu)造攝動(dòng)解,都不可能像一般受攝二體問題那樣簡單,無法滿足實(shí)際應(yīng)用的需求.

(2)即使形式上處理成一般的受攝二體問題,但由于探測器到月球的距離與探測器到地球的距離相當(dāng),這是處理第三體攝動(dòng)問題的一個(gè)難點(diǎn).

鑒于上述分析,對于地面測控而言,宜在J2000地心天球坐標(biāo)系中采用數(shù)值方法實(shí)現(xiàn)地-月系平動(dòng)點(diǎn)探測器的軌道外推,其數(shù)學(xué)模型就是一個(gè)表面形式上的受攝二體問題,狀態(tài)運(yùn)動(dòng)方程如下:

2 外力分析

對于航天器的軌道外推問題,實(shí)為一個(gè)定量問題,因涉及的各物理量有不同的量綱和大小,不便于問題的分析和表達(dá),擬采用無量綱單位(可簡稱為歸一化單位)來處理問題.以地球衛(wèi)星的軌道運(yùn)動(dòng)為背景,包括月球之類的遠(yuǎn)地衛(wèi)星(如地-月系的L1或L2點(diǎn)探測器)軌道,通常采用的歸一化單位系統(tǒng),即長度單位[L]、質(zhì)量單位[M]和時(shí)間單位[T]如下:

其中,GE是地心引力常數(shù),ae是地球參考橢球體的赤道半徑,時(shí)間單位[T]是導(dǎo)出單位,目的是使該計(jì)算單位系統(tǒng)中,引力常數(shù)G=1和中心天體引力常數(shù)μ=GE=1.如果地球引力模型采用當(dāng)今地固坐標(biāo)系所對應(yīng)的WGS-84系統(tǒng),則相應(yīng)的地球動(dòng)力學(xué)扁率為J2=1.082636022×10?3.

在采用上述歸一化單位后,狀態(tài)運(yùn)動(dòng)方程(1)式變?yōu)橄铝行问?

該式中的是月球、太陽等第三體的無量綱質(zhì)量,即

其中,GM和GS分別為月心和日心引力常數(shù).

其中,是航天器到第j個(gè)天體的位置矢量,是航天器到第j個(gè)天體的距離,是地心天球坐標(biāo)系中月球、太陽等第三體坐標(biāo)矢量.

2.1 各大天體的引力攝動(dòng)量級估計(jì)

地球的引力常數(shù)即GE=398600.4418 km3/s2,而月球、太陽、水星、金星、火星、木星和土星的引力常數(shù)相對地球的大小依次分別為0.0123000383、332946.050895、0.055273598、0.814998108、0.107446732、317.8942053和95.1574041.

攝動(dòng)量級的近似估計(jì)式為

該式中的即各大天體相對地球的質(zhì)量比,其值即上面給出的引力常數(shù)比,r和分別為探測器和攝動(dòng)天體到地球的距離,估計(jì)中可取的平均值,但因水星軌道的偏心率較大(e>0.2),對應(yīng)的r′值將分別按近日點(diǎn)和遠(yuǎn)日點(diǎn)兩種狀態(tài)設(shè)定平均距離.

考慮定點(diǎn)在地-月系共線平動(dòng)點(diǎn)L1和L2的探測器到地球的距離r分別為

其中是地-月平均距離.由(6)式給出上述各大天體對L1和L2點(diǎn)探測器軌道的攝動(dòng)量級分別如下:

2.2 地球非球形引力攝動(dòng)量級

主項(xiàng)(扁率J2)的攝動(dòng)量級估計(jì)

2.3 太陽光壓的攝動(dòng)量級

定點(diǎn)在地-月系共線平動(dòng)點(diǎn)L1和L2的一個(gè)通常尺度(包括質(zhì)量和承受光壓的等效截面)的探測器,太陽光壓攝動(dòng)的量級估計(jì)如下:

其中,κ=1.44,面質(zhì)比S/m=109,ρ⊙為1 au處的太陽輻射壓強(qiáng),ρ⊙=0.3169×10?17.

根據(jù)上述對外力因素的量級估計(jì),對于定點(diǎn)在地-月系共線平動(dòng)點(diǎn)L1和L2或其附近的探測器軌道,若考慮10?6以上的攝動(dòng)因素,相應(yīng)的力模型中只需要考慮如下攝動(dòng)源:月球、太陽和金星的質(zhì)點(diǎn)引力,地球非球形引力位的扁率J2和太陽光壓,其中最主要的是月球和太陽的質(zhì)點(diǎn)引力.

3 定點(diǎn)在地-月系L1和L2點(diǎn)的探測器軌道概貌

這里給出地-月系圓型限制性三體問題模型下兩個(gè)簡單的算例,初始?xì)v元為2016-09-30UTC0:00:0.0,對應(yīng)的TDT(地球力學(xué)時(shí))是57661.0007891667(MJD),探測器定點(diǎn)在地-月系L1和L2點(diǎn)的各一條軌道上.經(jīng)簡單坐標(biāo)轉(zhuǎn)換,即可獲得J2000.0地心天球坐標(biāo)系中相應(yīng)的兩條軌道的初始位置(x,y,z)、速度(˙x,˙y,˙z)和相應(yīng)的軌道根數(shù),分別列于表1和表2.

表1L 1點(diǎn)和L 2點(diǎn)的位置和速度Table 1 The positions and velocities of the points L1 and L2

表2L 1點(diǎn)和L 2點(diǎn)的軌道根數(shù)Table 2 The orbital elements of the points L1 and L2

兩條軌道的圖像見圖1–2,這表明在初始時(shí)刻瞬間,實(shí)際上都是一條偏心率較大的環(huán)繞地球的橢圓軌道,探測器均處于該軌道的遠(yuǎn)地點(diǎn)和近地點(diǎn)(對讀者而言,這一特點(diǎn)是容易理解的,無需做過多解釋),兩圖中的坐標(biāo)單位ae是地球參考橢球體的赤道半徑.就地-月+探測器系統(tǒng)而言,這都是初始瞬時(shí)軌道,而在月球的引力作用下,探測器與月球軌道“同步”做相同的圓軌道運(yùn)動(dòng).

4 平動(dòng)點(diǎn)軌道外推問題的相關(guān)分析

4.1 平動(dòng)點(diǎn)軌道分析采用的數(shù)學(xué)模型

第1節(jié)引言中已指出,應(yīng)在J2000地心天球坐標(biāo)系中處理其軌道運(yùn)動(dòng)問題,并采用數(shù)值方法實(shí)現(xiàn)相應(yīng)的軌道外推.為了定量顯示這類軌道外推中誤差傳播狀態(tài)的主要特征,顯然應(yīng)選擇地-月-日+探測器的質(zhì)點(diǎn)引力系統(tǒng),相應(yīng)的狀態(tài)運(yùn)動(dòng)方程即

其中是地球引力加速度:而和分別為月球和太陽的無量綱質(zhì)量,見前面的(4)式.相應(yīng)的月球和太陽引力攝動(dòng)加速度、的具體形式分別為

圖1L1點(diǎn)初始軌道在J2000地心天球坐標(biāo)系中(赤道面內(nèi))的圖像Fig.1 The figure(in the equatorial plane)of point L1initial orbit in the J2000 reference system

圖2L 2點(diǎn)初始軌道在J2000地心天球坐標(biāo)系中(赤道面內(nèi))的圖像Fig.2 The figure(in the equatorial plane)of pointL 2initial orbit in the J2000 reference system

4.2 平動(dòng)點(diǎn)軌道類型的選擇

這里考查的平動(dòng)點(diǎn)軌道,包括如下3種類型:

(1)初始時(shí)刻定點(diǎn)在地-月系的L1點(diǎn)或L2點(diǎn)處的平動(dòng)點(diǎn)軌道,以下簡稱該類型軌道為L1點(diǎn)軌道或L2點(diǎn)軌道;

(2)初始時(shí)刻定點(diǎn)在地-月系的L1點(diǎn)或L2點(diǎn)附近的halo軌道上,以下簡稱該類型軌道為L1點(diǎn)暈軌道或L2點(diǎn)暈軌道;

(3)初始時(shí)刻定點(diǎn)在地-月系的L1點(diǎn)或L2點(diǎn)附近的Lissajous軌道上,以下簡稱該類型軌道為L1點(diǎn)Lissajous軌道或L2點(diǎn)Lissajous軌道.

經(jīng)初步設(shè)計(jì)(對應(yīng)所采用的質(zhì)點(diǎn)引力系統(tǒng))分別給出6條軌道,在J2000地心天球坐標(biāo)系中各對應(yīng)的軌道初值如下:所有初始時(shí)刻對應(yīng)歷元為2016-09-30UTC0:00:0.0(相應(yīng)TDT的MJD為57661.0007891667),位置、速度和相應(yīng)的軌道根數(shù)分別列于表3和表4.表中的軌道類型1、2、···、6依次為L1點(diǎn)軌道、L1點(diǎn)暈軌道、L1點(diǎn)Lissajous軌道、L2點(diǎn)軌道、L2點(diǎn)暈軌道、L2點(diǎn)Lissajous軌道,表4–8類同.

表3 6條軌道的位置和速度Table 3 The positions and velocities of 6 oribits

表4 6條軌道的軌道根數(shù)Table 4 The orbital elements of 6 oribits

下面首先對上述6類軌道作7d和27d的軌道外推,給出一個(gè)誤差傳播的定量輪廓,在此基礎(chǔ)上再作定性分析.

4.3 6條平動(dòng)點(diǎn)軌道外推7d的狀態(tài)

為簡單起見(也不失一般性),在考查誤差傳播中,將初始誤差全部集中在最重要的軌道半長徑上(根據(jù)目前定軌的實(shí)際狀況,軌道半長徑的精度為10 m量級),7d的軌道外推結(jié)果列于表5–6.

表5 平動(dòng)點(diǎn)軌道外推7d的軌道狀態(tài)Table 5 The states of the libration point orbits propagated for 7days

表6 平動(dòng)點(diǎn)軌道外推7d的空間位置狀態(tài)Table 6 The positions of the libration point orbits prop a gated for 7days

4.4 6條軌道外推27d的狀態(tài)

盡管探測器的定點(diǎn)只是近似的,實(shí)際運(yùn)行過程中必須通過不斷軌控才能保持這類軌道,外推弧段增長至27d,只是為了進(jìn)一步了解這類特殊軌道的動(dòng)力學(xué)特征及其相應(yīng)的誤差傳播規(guī)律.初始誤差仍全部集中在軌道半長徑上,27d的軌道外推結(jié)果列于表7–8,其中,第2和3兩類軌道(即L1點(diǎn)暈軌道和L1點(diǎn)Lissajous軌道)只外推了22 d,其原因?qū)⒃?.5小節(jié)中具體說明.

表7 平動(dòng)點(diǎn)軌道外推27d的軌道狀態(tài)Table 7 The states of the libration point orbits propagated for 27days

對于探測器的軌道運(yùn)動(dòng)而言,通常所說的長期位置預(yù)報(bào)和短期位置預(yù)報(bào)中的長期或短期并不是簡單的時(shí)間間隔,而是運(yùn)行弧段的長短.因此,為了比較上述各條軌道之間外推誤差傳播的定量大小,需要了解它們的軌道運(yùn)行周期,這6條特殊軌道的初始運(yùn)行周期TS值依次為

L1類:TS=13.708768 d,12.221658 d,13.113512 d;

L2類:TS=126.58784 d,58.847196 d,93.089267d.

由此便于了解軌道外推7d和27d對上述6條軌道所對應(yīng)的弧段長短.

在已知6條軌道自身的動(dòng)力學(xué)特性和軌道外推弧段長短的前提下,不難看出表5所列出的外推位置誤差所反映的一些動(dòng)力學(xué)規(guī)律,基本上可歸納如下:

(1)軌道外推7d均為短弧,位置誤差都在1 km以內(nèi);

(2)對于L1點(diǎn)軌道和L2點(diǎn)軌道,外推7d或27d,位置誤差的累積仍不嚴(yán)重,其誤差傳播的特征就是一般Kepler運(yùn)動(dòng)特征的反映.由于初值誤差(?a0=10 m)較小,實(shí)為小擾動(dòng),既不會(huì)激發(fā)其初值不穩(wěn)定的固有特征,又不會(huì)明顯改變短弧段誤差累積的效果,反而周期性的效果比長期累積效應(yīng)更明顯,見表5中軌道1(即L1點(diǎn)軌道)的誤差定量狀態(tài).

(3)對于L1點(diǎn)和L2點(diǎn)的halo軌道和Lissajous軌道而言,由于其位置已經(jīng)“遠(yuǎn)離”不穩(wěn)定平動(dòng)點(diǎn)L1和L2,而嚴(yán)格的halo軌道和Lissajous軌道設(shè)計(jì)又無法實(shí)現(xiàn),探測器的定點(diǎn)只是近似的,在同樣是?a0=10 m的初值誤差情況下,已不能再簡單地只看作對halo軌道和Lissajous軌道的小擾動(dòng),而更重要的是對平動(dòng)點(diǎn)的大擾動(dòng)起作用,在不太長的外推弧段內(nèi),其平動(dòng)點(diǎn)本身的不穩(wěn)定特征即顯現(xiàn)無遺.見表5中的軌道2、軌道3、軌道5和軌道6,特別是軌道2和3,相對而言,27d的弧段顯得更長,外推弧段超過22 d后,其軌道偏心率很快就會(huì)達(dá)到e≈1.0的狀態(tài).

表8 平動(dòng)點(diǎn)軌道外推27d的空間位置狀態(tài)Table 8 The positions of the libration point orbits propagated for 27days

4.5 平動(dòng)點(diǎn)軌道運(yùn)行中誤差傳播狀態(tài)的基本特征

從上述6條軌道的外推計(jì)算結(jié)果已能看出拉格朗日點(diǎn)軌道外推中位置誤差傳播的基本特征,首先將由?a0=10 m導(dǎo)致的位置誤差集中列于表9.

表9 拉格朗日點(diǎn)軌道外推位置誤差的定量狀態(tài)(單位:km)Table 9 The quantitative state of position error s of the Lagrange point orbit propagation(unit:km)

綜上幾小節(jié),就平動(dòng)點(diǎn)探測器軌道運(yùn)行誤差傳播狀態(tài)的簡單計(jì)算和分析,可以表明:在地-月系L1和L2平動(dòng)點(diǎn)軌道設(shè)計(jì)中,確實(shí)很難實(shí)現(xiàn)較長弧段的無動(dòng)力控制運(yùn)行軌道,而不是設(shè)計(jì)者本身的問題.那么,在這樣較短的弧段內(nèi),就軌道預(yù)報(bào)而言,要達(dá)到較高的精度是容易實(shí)現(xiàn)的,本節(jié)最后一小節(jié)將會(huì)給出具體算例,并有實(shí)測結(jié)果的檢驗(yàn).

4.6L1點(diǎn)暈軌道的運(yùn)行狀態(tài)

為了進(jìn)一步揭示平動(dòng)點(diǎn)飛行器軌道外推中誤差傳播的動(dòng)力學(xué)機(jī)制,這里再給出一條與上述L1點(diǎn)暈軌道相近的環(huán)繞地球運(yùn)行的大橢圓逆行軌道,同樣對應(yīng)一遠(yuǎn)地衛(wèi)星.初始時(shí)刻仍為2016-09-30UTC0:00:0.0,初始位置和速度及相應(yīng)的軌道根數(shù)分別列于表10和表11.該軌道的初始運(yùn)行周期為12.208831 d.對此軌道同作27d的外推,計(jì)算結(jié)果列于表12–13.

表10L1點(diǎn)暈軌道的位置和速度Table 10 The position and velocity of the halo orbit at pointL 1

表11L 1點(diǎn)暈軌道的軌道根數(shù)Table 11 The orbital elements of the halo or bit at pointL 1

表12L 1點(diǎn)暈軌道外推27d的軌道狀態(tài)Table 12 The states of the halo orbit prop a gated for 27days at pointL 1

表13L 1點(diǎn)暈軌道外推27d的空間位置狀態(tài)Table 13 The positions of the halo orbit prop a gated for 27days at pointL 1

這樣一條有別于L1點(diǎn)暈軌道繞地運(yùn)行的大橢圓逆行軌道,本質(zhì)上就是一條普通的Kepler軌道,其固有的初值不穩(wěn)定性在繞地運(yùn)行27d(僅2圈)的短弧內(nèi)不會(huì)有明顯體現(xiàn),這是一個(gè)常規(guī)問題,本文引進(jìn)這一算例的目的,是從另一角度體現(xiàn)地-月系平動(dòng)點(diǎn)軌道誤差快速傳播的固有不穩(wěn)定性特征.至于逆行軌道自身的運(yùn)動(dòng)規(guī)律及其動(dòng)力學(xué)特征,已超出本文論述的范疇,不再介紹,如有需要,可見文獻(xiàn)[4]及其有關(guān)作者的研究工作.

4.7 平動(dòng)點(diǎn)軌道的定軌和短弧預(yù)報(bào)精度的檢驗(yàn)

既然嚴(yán)格的拉格朗日點(diǎn)軌道設(shè)計(jì)無法實(shí)現(xiàn),探測器的定點(diǎn)只是近似的,運(yùn)行過程中必須通過不斷的軌控才能保持,那么,對于地面測控和星上控制,只能從短弧角度來考慮問題.

對于平動(dòng)點(diǎn)軌道的短弧定軌和軌道預(yù)報(bào)而言,相對地球低軌或高軌衛(wèi)星的同類問題,實(shí)無任何特殊困難和特別需要處理的難題.采用南京大學(xué)空間環(huán)境與航天動(dòng)力學(xué)研究所自主編寫的定軌軟件和利用國內(nèi)的USB(Unified S-band)測量數(shù)據(jù),在沒有任何其他輔助信息的前提下,對嫦娥3號的相關(guān)任務(wù)探測器進(jìn)行了定軌,并與北京航天飛行控制中心的事后定軌結(jié)果作了對比.在此定軌的基礎(chǔ)上,采用非常簡單的數(shù)值外推方法(只考慮地、月、日三體的質(zhì)點(diǎn)引力和簡單的光壓模型,外推中的6個(gè)軌道初值采用相關(guān)任務(wù)的定軌結(jié)果)進(jìn)行了軌道預(yù)報(bào),毫無困難地達(dá)到了較高精度.略去不必要的細(xì)節(jié)說明,將有關(guān)結(jié)果一并列于表14–15.

表14L 2點(diǎn)暈軌道短弧外推3 d與事后精密定軌結(jié)果的比較Table 14 The comparison of two precise orbit determination method s for the pointL 2orbit after propagating 3 days

表15L 2點(diǎn)暈軌道短弧外推7d與事后精密定軌結(jié)果的比較Table 15 The comparison of two precise orbit determination methods for the pointL 2orbit after prop a gating 7days

表14–15中的精密定軌A和B分別對應(yīng)北京航天飛行控制中心和南京大學(xué)空間環(huán)境與航天動(dòng)力學(xué)研究所的結(jié)果,表中的結(jié)果基本上已能說明問題,但為了讓讀者對這類探測器的定軌和外推精度有更清晰的了解,下面進(jìn)一步作些必要的說明.

(1)關(guān)于光壓模型,在不了解探測器的具體細(xì)節(jié)情況下,作者們根據(jù)獨(dú)立定軌中獲得的有關(guān)估計(jì)值,獲得了包括衛(wèi)星表面熱性能在內(nèi)的等效面質(zhì)比,從而給出了相應(yīng)的經(jīng)驗(yàn)?zāi)Ｐ?一個(gè)等效的平面模型.

(2)盡管沒有具體給出兩個(gè)單位的定軌(包括測量數(shù)據(jù))細(xì)節(jié),但表14–15所給出的計(jì)算結(jié)果,已能說明本文要體現(xiàn)的這類特殊軌道的定軌和外推精度了.因?yàn)楸碇械慕Y(jié)果是外推3 d和7d與事后精密定軌結(jié)果的比較,且兩個(gè)單位的定軌結(jié)果之差基本上在500 m之內(nèi),這樣的比較更能體現(xiàn)兩個(gè)單位定軌結(jié)果的真實(shí)性以及本文所采用的力模型的合理性.

上述計(jì)算結(jié)果和兩點(diǎn)補(bǔ)充說明充分表明:盡管這類探測器的軌道特殊,初值誤差的傳播程度遠(yuǎn)比一般的環(huán)繞型探測器的軌道顯著,但相應(yīng)的短弧定軌和高精度軌道預(yù)報(bào)并無特殊困難.

4.8 關(guān)于兩種坐標(biāo)系中平動(dòng)點(diǎn)探測器軌道形式的轉(zhuǎn)換問題

就定軌和預(yù)報(bào)的需求,顯然是在J2000地心天球坐標(biāo)系中進(jìn)行相關(guān)問題的處理,而對這類具有特殊性質(zhì)的軌道,軌控又必須通過相應(yīng)的地-月系旋轉(zhuǎn)坐標(biāo)系來處理.這就涉及到兩種坐標(biāo)系之間的轉(zhuǎn)換問題,其本身是容易實(shí)現(xiàn)的,而在具體的航天任務(wù)中,各有關(guān)部門根據(jù)實(shí)際需求,對相應(yīng)的地-月系旋轉(zhuǎn)坐標(biāo)系實(shí)有不同取法,故這里不再做相應(yīng)討論.