劉林,楊志濤

(1.中國(guó)科學(xué)院國(guó)家天文臺(tái),北京 100101;2.南京大學(xué)天文與空間科學(xué)學(xué)院,南京 210093;3.中國(guó)科學(xué)院大學(xué),北京 100049)

1 引言

對(duì)太陽系近地天體 (主要是近地小行星)軌道運(yùn)動(dòng)的了解是防御的前提,而準(zhǔn)確的軌道確定和預(yù)報(bào)則取決于采用的動(dòng)力學(xué)模型。太陽系近地天體繞日運(yùn)動(dòng)的力學(xué)環(huán)境,與近地衛(wèi)星所處的力學(xué)環(huán)境有顯著差別,其外力源除中心天體太陽的質(zhì)點(diǎn)引力外,主要攝動(dòng)源來自各大行星的引力作用。至于表面力之類的攝動(dòng)源,不外乎太陽風(fēng)之類的耗散效應(yīng)和源于太陽輻射的各種力學(xué)效應(yīng)。對(duì)近地天體的軌道運(yùn)動(dòng)而言,前者無需討論,而后者歸結(jié)為輻射壓?jiǎn)栴},但太陽實(shí)為弱輻射恒星,與之相關(guān)的動(dòng)力學(xué)問題,特別是對(duì)太陽系內(nèi)各大小天體 (包括彗星)運(yùn)動(dòng)的影響,并不重要,多年前作者就作過分析[1],本文將進(jìn)一步給出該影響的定量體現(xiàn),為處理近地天體軌道運(yùn)動(dòng)問題,提供準(zhǔn)確的動(dòng)力學(xué)模型和相應(yīng)的計(jì)算方法。

2 日心系中小天體軌道運(yùn)動(dòng)的力學(xué)模型

2.1 小天體軌道運(yùn)動(dòng)的受攝運(yùn)動(dòng)方程

清晰的力學(xué)背景與相應(yīng)的數(shù)學(xué)模型,是準(zhǔn)確反映小天體軌道運(yùn)動(dòng)的重要前提?？紤]到小天體的質(zhì)量確實(shí)較小,可以引用限制性 (N+1)體模型,即小天體的質(zhì)量可以忽略 (m=0),N體包括中心天體太陽、各大行星、冥王星和月球。

通常采用的計(jì)算單位系統(tǒng)中,長(zhǎng)度單位[L],質(zhì)量單位 [M]和時(shí)間單位 [T]分別為

其中,AU=1.459 978 707 00×108(km)是天文單位,S=GS/G=1.9884×1030kg是太陽質(zhì)量,時(shí)間單位 [T]是導(dǎo)出單位。這種處理,一方面是為了各物理量的單位在一定程度上歸一化,便于對(duì)各有關(guān)物理量的量級(jí)進(jìn)行估計(jì)和比較,同時(shí)也可簡(jiǎn)化計(jì)算公式的表達(dá)。

在日心黃道坐標(biāo)系中,小天體的運(yùn)動(dòng)方程表達(dá)如下[2,3]：

在上述歸一化的無量綱單位系統(tǒng)中,有μ=G(S+m)=GS=1,各大行星、冥王星和月球的無量綱質(zhì)量以及日心距離分別記作和即

于是小天體軌道運(yùn)動(dòng)的受攝運(yùn)動(dòng)方程式 (2)可寫成下列形式：

該式中的攝動(dòng)項(xiàng)APN是源于太陽引力的后牛頓效應(yīng)產(chǎn)生的攝動(dòng)加速度,表達(dá)形式如下：

如果小天體離某大行星近到一定程度,還需要考慮相應(yīng)大行星的主要非球形項(xiàng) (如扁率J2)。至于與太陽輻射有關(guān)的表面力效應(yīng),對(duì)于已基本固化的小行星,其質(zhì)量密度與地球接近,該影響并不重要,對(duì)于小行星監(jiān)測(cè)預(yù)警和安全防御而言,不必夸大與其有關(guān)的輻射效應(yīng),后面2.4小節(jié)將會(huì)具體闡明并給出定量檢驗(yàn)。

2.2 各大天體的質(zhì)量與離日平均距離

這可參考下列表2列出的相關(guān)平均根數(shù)獲得。

關(guān)于月球離日的平均距離,在對(duì)近地天體受力估計(jì)中可取地球的值。

表1 質(zhì)量比Tab.1 Mass ratio

表1 質(zhì)量比Tab.1 Mass ratio

大行星 質(zhì)量比 大行星 質(zhì)量比 大行星 質(zhì)量比水星 1.660136795×10-7 火星 3.227156038E×10-7 天王星 4.366244043×10-5金星 2.447838340×10-6 木星 0.954791898E×10-3 海王星 5.151383713×10-5地球 3.003489616×10-6 土星 2.858856701E×10-4 冥王星 0.732246679×10-8月球 3.694303465×10-8

表2 各大行星平均軌道的半長(zhǎng)徑和偏心率i, i(i=1,2,…)Tab.2 Semi-major axis and eccentricity of mean orbit of planetsi, i(i=1,2,…)

表2 各大行星平均軌道的半長(zhǎng)徑和偏心率i, i(i=1,2,…)Tab.2 Semi-major axis and eccentricity of mean orbit of planetsi, i(i=1,2,…)

大行星 ai(AU) e-i 大行星 ai(AU) e-i水星 0.38709831 0.20563175 土星 9.5549096 0.05550862金星 0.72332982 0.00677188 天王星 19.2184461 0.04629590地球 1.00000102 0.01670862 海王星 30.1103869 0.00898809火星 1.52367934 0.09340063 冥王星 39.4816868 0.24880766木星 5.20260319 0.04849485

2.3 各大行星引力攝動(dòng)加速度的量級(jí)估計(jì)

下面針對(duì)近地天體,對(duì)各大行星的引力攝動(dòng)加速度進(jìn)行量級(jí)估計(jì)。

(1)引力攝動(dòng)加速度的近似估計(jì)

對(duì)于離太陽和運(yùn)動(dòng)小天體較遠(yuǎn)的大行星 (土星、天王星等),其攝動(dòng)量級(jí)的估計(jì)與地球衛(wèi)星的日、月攝動(dòng)類似,同樣可采用下列估計(jì)式：

否則按下式估計(jì)：

上兩式中的m′即上述各大行星相對(duì)太陽的質(zhì)量比,r和r′各為小天體和攝動(dòng)天體到中心天體太陽的距離,攝動(dòng)天體到太陽的距離即可取上述表2中給出的平均軌道半長(zhǎng)徑值,對(duì)軌道偏心率較大的水星和冥王星,可采用其近日距或遠(yuǎn)日距。

對(duì)于近地天體,上述各大行星的攝動(dòng)量級(jí)為ε=O( 10-9～10-4)。

(2)后牛頓效應(yīng)的估計(jì)

這是引力攝動(dòng)的另一類效應(yīng),近似估計(jì)如下：

在太陽系中,對(duì)于水星繞日運(yùn)動(dòng)而言,這一攝動(dòng)量級(jí)為10-7,而對(duì)離太陽1～2個(gè)天文單位的近地天體,其攝動(dòng)量級(jí)為10-8。該效應(yīng)的特征是對(duì)小天體運(yùn)動(dòng)的軌道平面無影響,主要是對(duì)軌道近日點(diǎn)幅角ω有長(zhǎng)期效應(yīng),這正是廣義相對(duì)論的幾個(gè)論據(jù)之一,水星近日點(diǎn)進(jìn)動(dòng)的檢驗(yàn),即牛頓力學(xué)無法解釋的水星近日點(diǎn)進(jìn)動(dòng)的偏差。

2.4 小天體繞日運(yùn)動(dòng)中承受的表面力問題

對(duì)于本文論述的問題,太陽系中各天體運(yùn)動(dòng)所承受的表面力主要源于太陽輻射壓,相應(yīng)的量級(jí)估計(jì)式如下：

該估計(jì)式中的μ=G(S+m)=GS=1,是小天體所在位置 (日心距為r)的太陽輻射壓強(qiáng)度,ρs是地球附近rE=1AU處的太陽輻射壓強(qiáng)度,有

在日心系中的無量綱值 (本文采用的計(jì)算單位)為

面質(zhì)比 (S/M)的常用單位是m2/kg,注意,這里的m是長(zhǎng)度單位米的符號(hào)。在日心系統(tǒng)中面質(zhì)比換算為無量綱值的公式為

下面舉幾個(gè)小行星的例子：

(1)小行星65803(Didymos)

Didymos是Xk型近地雙小行星,其主星的直徑約為800m,平均密度為ρ=1.7 g/cm3。作為估計(jì),按球體形態(tài)考慮,相應(yīng)的圓球截面積和球體形態(tài)的質(zhì)量分別為

那么,面質(zhì)比 (S/M)為

相應(yīng)的太陽輻射壓攝動(dòng)量級(jí)為

(2)小行星10302(1989ML)

直徑約為600m,平均密度為2.0g/cm3,質(zhì)量可估為0.565×1011kg。相應(yīng)的太陽輻射壓攝動(dòng)量級(jí)為

(3)小行星175706(1996FG3)

直徑300m,質(zhì)量2.7×1010kg,相應(yīng)的太陽輻射壓攝動(dòng)量級(jí)為

表3 兩顆小行星的軌道根數(shù)Tab.3 Orbit elements of two asteroids

(4)小行星99942(Apophis)

直徑1.9km,質(zhì)量2.1×1012kg,密度1.4 g/cm3,相應(yīng)的太陽輻射壓攝動(dòng)量級(jí)為

(5)小行星433(Eros)

質(zhì)量6.69×1015kg,直徑16.8km,面質(zhì)比作為球形估計(jì),相應(yīng)的太陽輻射壓攝動(dòng)量級(jí)為

從上述幾個(gè)估計(jì)算例可以了解到,對(duì)于一般的近地小行星,作為表面力攝動(dòng)源的太陽輻射壓確實(shí)不重要。至于輻射壓的大小,這從上述估計(jì)式 (10)不難看出：在質(zhì)量密度相近的情況下,尺度越小的天體,承受的輻射壓越大。

2.5 小天體繞日運(yùn)動(dòng)中承受表面力作用大小的數(shù)值檢驗(yàn)

這里考慮兩顆小行星,一顆近地小行星99942和一顆主帶小行星261936,看其在軌道演化中太陽輻射壓是否有明顯作用。兩顆星的初始?xì)v元時(shí)刻分別對(duì)應(yīng)UTC：2000年1月01日12時(shí)0分0.0秒,2016年7月31日0時(shí)0分0.0秒,相應(yīng)時(shí)刻的軌道根數(shù)列于表3。

表3中半長(zhǎng)徑a的單位是AU,其它4個(gè)角度的單位是 (°)。對(duì)上述兩條軌道作一個(gè)繞日運(yùn)行周期的軌道外推,給出相應(yīng)的軌道變化結(jié)果和各攝動(dòng)因素影響的定量狀態(tài)。

(1)近地小行星99942

面質(zhì)比(S/m)=10-6(m2/kg),運(yùn)行一個(gè)周期 (395.238592天),結(jié)果列于表4和表5。

(2)主帶小行星261936

面質(zhì)比(S/m)=10-6(m2/kg),運(yùn)行一個(gè)周期 (1487.960956天),結(jié)果列于表6和表7。

上列各表中的力模型序號(hào) (1)-(10),依次為

表4 瞬時(shí)根數(shù)狀態(tài)Tab.4 Instantaneous elements states

表5 空間位置、速度狀態(tài)Tab.5 Space position and speed states

表6 瞬時(shí)根數(shù)狀態(tài)Tab.6 Instantaneous elements states

表7 空間位置、速度狀態(tài)Tab.7 Space position and speed states

序號(hào)中的1或0,分別表示考慮或不考慮對(duì)應(yīng)的攝動(dòng)源,相應(yīng)的12個(gè)攝動(dòng)源,依次對(duì)應(yīng)8個(gè)大行星 (水星,金星,…,海王星)、冥王星、月球、后牛頓效應(yīng)和太陽輻射壓。

從上述結(jié)果不難看出：前面對(duì)各攝動(dòng)影響的定量分析是準(zhǔn)確的,對(duì)于基本固化的小天體,面質(zhì)比較小,就軌道運(yùn)動(dòng)而言,表面力之類的外力攝動(dòng)影響確實(shí)無需考慮,特別在空間防御領(lǐng)域,不應(yīng)受其“干擾”。作者曾采用上述力模型對(duì)應(yīng)的受攝運(yùn)動(dòng)方程 (2),處理過“近地小行星與地球交會(huì)問題”的研究[4],所獲結(jié)果與國(guó)際小行星中心 (MPC)公布的結(jié)果相符。

3 小天體軌道運(yùn)動(dòng)解擬采用的數(shù)學(xué)方法[5,6]

與人造地球衛(wèi)星的運(yùn)動(dòng)狀態(tài)不同,對(duì)于近地天體的軌道運(yùn)動(dòng),若要構(gòu)造攝動(dòng)分析解 (即小參數(shù)冪級(jí)數(shù)解),將會(huì)涉及第三體攝動(dòng)函數(shù)展開中出現(xiàn)的外攝 (r/r′＜1)和內(nèi)攝 (r/r′＞1)情況,特別是兩種狀態(tài)出現(xiàn)在同一攝動(dòng)天體的狀態(tài),這給攝動(dòng)分析解的具體構(gòu)造帶來困難。

鑒于近地天體的運(yùn)動(dòng)不同于人造地球衛(wèi)星,幾乎一年左右才繞日運(yùn)行一圈,軌道預(yù)報(bào)的緊迫性不同于近地空間目標(biāo)監(jiān)測(cè)防御的要求,擬直接采用數(shù)值方法實(shí)現(xiàn)軌道預(yù)報(bào)計(jì)算。對(duì)于近地天體軌道運(yùn)動(dòng)方程的數(shù)值解,其狀態(tài)量的選擇可考慮采用相應(yīng)的軌道根數(shù)σ,其運(yùn)動(dòng)方程表達(dá)如下：

考慮到算法的普適性,擬采用同時(shí)消除小偏心率和小傾角的無奇點(diǎn)根數(shù)σ(即適用于0≤e＜1.0,0°≤i＜180°.0) 作為上述運(yùn)動(dòng)方程 (19)的積分狀態(tài)量。所有攝動(dòng)源對(duì)應(yīng)的攝動(dòng)加速度分量S,T,W可由Fε形成,即

采用上述處理方法,導(dǎo)致整個(gè)計(jì)算過程與直接采用坐標(biāo)、速度(r,)的流程沒有明顯差別,但卻避開了原運(yùn)動(dòng)方程 (2)或 (4)右端無攝項(xiàng)的計(jì)算,該項(xiàng)既是數(shù)值積分過程中制約積分步長(zhǎng)增大的關(guān)鍵項(xiàng),又是積分截?cái)嗾`差傳播最嚴(yán)重的誤差源。有關(guān)具體細(xì)節(jié)不再介紹,僅供參考。