田偉

(上海歷元信息科技有限公司,上海 200241)

1 引言

近地小天體 (99942)Apophis于2004年6月19日首先被R.A.Tucker、D.J.Tholen和F.Bernardi發(fā)現(xiàn),因最早期的軌道分析顯示該顆小行星在2029年4月13日撞擊地球的概率超過了2%,得到了廣泛關(guān)注[1]。隨著對Apophis觀測資料的不斷積累,尤其是高精度的雷達(dá)數(shù)據(jù),研究人員發(fā)現(xiàn)該顆小行星撞擊地球的概率遠(yuǎn)小于早期預(yù)報(bào)[2-4]。在2021年3月Apophis再次飛掠地球時(shí),研究人員利用這次飛掠期間的觀測資料以及歷史資料已基本排除2068年之前Apophis撞擊地球的可能性[5]。

類似于Apophis,很多近地小天體會多次飛掠地球。由于受地月系統(tǒng)的影響,天體軌道和自轉(zhuǎn)狀態(tài)都會受到不同程度的改變。Brozovic等人[4]的研究結(jié)果顯示2029年Apophis飛掠地球期間,其日心半長徑將從0.92AU變化為1.10AU(相應(yīng)地,其軌道類型將從Aten型變化為Apollo型)。Souchay等人[6]的研究結(jié)果表明Apophis在2029年飛掠地球期間,在地月系統(tǒng)的引力作用下,其自轉(zhuǎn)狀態(tài)也將發(fā)生顯著改變。因此,在利用實(shí)測數(shù)據(jù)確定Apophis軌道之前,本文將定量分析2004-2022年期間幾次飛掠地球過程中太陽系歷表的誤差和相關(guān)動力學(xué)參數(shù)的不確定度對Apophis軌道的影響。圖1列出了觀測弧段期間Apophis飛掠地球、月球和金星的概況。

圖1 1994-2022期間,小行星 (99942)Apophis飛掠地球、月球和金星的概況Fig.1 Overview of encounters of asteroid(99942)Apophis with the Earth,Moon and Venus from 1994 to 2022

2 觀測資料

國際小行星中心 (MPC)收集了關(guān)于Apophis的地面光學(xué)觀測和無線電雷達(dá)測量數(shù)據(jù)。文中采用了截至2021年5月20日該網(wǎng)站公布的所有光學(xué)數(shù)據(jù),以及2005-2021年期間由美國Arecibo和Goldstone雷達(dá)站獲取的高精度測距和多普勒頻移觀測數(shù)據(jù)[4]。表1給出了部分觀測數(shù)據(jù)的統(tǒng)計(jì)信息。

表1 小行星 (99942)Apophis(部分)地面光學(xué) (赤經(jīng)/赤緯)和雷達(dá) (測距/多普勒頻移)觀測數(shù)據(jù)的統(tǒng)計(jì)Table 1 Statistics of partial ground-based optical observations(right ascension and declination)and radar observations(range and doppler frequency shift)of asteroid(99942)Apophis

3 動力學(xué)模型

文中考慮的天體包括太陽、八大行星、月球和冥王星,以及48顆大質(zhì)量的主帶小行星。動力學(xué)模型考慮的因素包括:(1)點(diǎn)質(zhì)量相互作用:EIH方程; (2)非球形引力場:太陽、地球和月球的低階引力場;(3)主帶小行星的攝動效應(yīng); (4)Yarkovsky效應(yīng)對Apophis的攝動影響等。

3.1 點(diǎn)質(zhì)量相互作用

因中心天體太陽與Apophis的相互作用中相對論效應(yīng)十分顯著,需要考慮點(diǎn)質(zhì)量相互作用中的一階后牛頓近似下的相對論效應(yīng)。太陽、月球、八大行星、冥王星對小行星 (99942)Apophis的引力作用由Einstein-Infeld-Hoffmann(EIH)方程給出[7]:

式中,下標(biāo)A表示Apophis,GMB和GMC分別是天體B和C的質(zhì)量與引力常數(shù)G的乘積,rAB=rB-rA和rAB=|rAB|分別是天體A和B之間的相對位置矢量和相對距離,rA(或rB)和vA(或vB)分別為對應(yīng)天體A(或B)的位置和速度矢量。相對位置矢量rAC和rBC的定義類似于rAB,而相對距離rAC和rBC的定義類似于rAB,以上各量均在太陽系質(zhì)心廣域參考系中給出。c=299792458m·s-1為光子在真空中的傳播速度。一階后牛頓參數(shù)β=γ=1。

3.2 非球形引力場

除考慮主要天體 (作為點(diǎn)質(zhì)量)的引力作用外,文中也考慮了地球、月球和太陽的非球形引力場對Apophis軌道的影響。非球形引力場的引力作用取決于Apophis與非球形天體之間距離。通常近地小行星會 (近似)周期性地與地球或金星交會 (見圖1)。以地球?yàn)槔?其引力場高階項(xiàng) (J2,Cnm,Snm)對某外部天體的影響可由式 (2)表示[8]

式中,在ξ-η-ζ坐標(biāo)系中,ξ坐標(biāo)軸代表從地球到Apophis的指向,ξ-η平面包含了地球自轉(zhuǎn)軸,ζ軸與ξ和η軸形成右手坐標(biāo)系。GME和R分別代表了地球的質(zhì)量與引力常數(shù)乘積和平均半徑。表示n階m次的締合Legendre函數(shù)。n1和n2分別是Jn和 (Cnm,Snm)的最高階數(shù)。(r,φ,λ)是從地球到Apophis構(gòu)成的位置向量在地固坐標(biāo)系中的球面坐標(biāo)。需要說明的是,在式 (2)的基礎(chǔ)上計(jì)算太陽系質(zhì)心坐標(biāo)系下的相互作用力時(shí),需要實(shí)施一系列坐標(biāo)轉(zhuǎn)換。太陽、地球和月球的非球形引力場分別取J2項(xiàng)、3階以下項(xiàng)和2階以下項(xiàng)。限于篇幅,文中未考慮金星非球形引力場對Apophis的影響。

3.3 主帶小行星攝動模型

某一特定小行星對Apophis的攝動效應(yīng)取決于攝動小行星的質(zhì)量和相對距離。目前被精確測量質(zhì)量參數(shù)的小行星在已發(fā)現(xiàn)的小行星總數(shù)中占比較小,因而小行星攝動效應(yīng)是動力學(xué)模型中一個難以確定的因素。Giorgini等人[2]統(tǒng)計(jì)了約37.3萬顆已知小行星對Apophis軌道的攝動效應(yīng),并給出了4顆攝動影響最大的小行星:Ceres、Vesta、Pallas和Hygiea,其中Ceres、Pallas和Vesta是主帶小行星中質(zhì)量最大的3顆。文中未對單顆小行星的攝動效應(yīng)做逐個篩選,而是選取了48顆大質(zhì)量小行星作為攝動小天體。

3.4 Yarkovsky效應(yīng)

由于目前對Apophis的觀測主要依賴于地面設(shè)備,其形狀、自轉(zhuǎn)軸指向、平均密度、反照率、熱慣量等基本幾何和物理參數(shù)確定精度偏低,在計(jì)算Apophis受到的Yarkovsky效應(yīng)時(shí),文中采用了簡化的非保守力線性模型[3,9-10],

式中,標(biāo)準(zhǔn)參數(shù)A1、A2和A3分別代表與徑向、橫切線方向和軌道面法線方向加速度分量相關(guān)的參數(shù),r和r分別為小行星相對太陽的位置向量和距離,h=r×r˙。文中僅考慮了橫切線方向的作用力,即A1=A3=0,A2≠0。線性模型中的非零參數(shù)A2將與軌道參數(shù)一起解算。

4 部分太陽系歷表參數(shù)對Apophis軌道的影響

作為太陽系內(nèi)主要天體,太陽、地球和月球的質(zhì)量不確定度對Apophis軌道有著不同程度的影響。目前太陽、地球和月球的質(zhì)量參數(shù)可以通過不同技術(shù)手段精確確定[11-14],其數(shù)值和不確定度見表2。圖2給出了相應(yīng)天體質(zhì)量不確定度對Apophis軌道的攝動效應(yīng)。

表2 太陽、地球和月球的質(zhì)量參數(shù) (單位:km3/s2)Table 2 The mass parameters of the sun,Earth and Moon(Unit:km3/s2)

圖2 由太陽、地球和月球的質(zhì)量參數(shù)不確定度引起的小行星 (99942)Apophis軌道變化 (單位:km)Fig.2 Orbit changes of asteroid(99942)Apophis induced by uncertainties of mass parameters for the sun,Earth and Moon(Unit:km)

在飛掠大行星過程中,被飛掠天體的歷表位置誤差是Apophis軌道確定的一個誤差項(xiàng)。其影響程度取決于最近飛掠距離,圖3給出了2004-2022年期間,地球的歷表位置誤差 (假設(shè)位置誤差為1km)對Apophis軌道的攝動效應(yīng)。相比于與地球歷表誤差的攝動效應(yīng),月球的影響較小。

圖3 地球歷表誤差對Apophis軌道的影響 (假設(shè)位置誤差為1km)(單位:km)Fig.3 Influences on orbit of asteroid(99942)Apophis by ephemeris error of the Earth(the supposed position error is 1 km)(Unit:km)

文中考慮了主要天體非球形引力場的影響,其中太陽、地球和月球的J2項(xiàng)是主要影響項(xiàng)。與被飛掠天體的歷表位置誤差類似,其非球形引力場的引力攝動效應(yīng)與最近飛掠距離直接相關(guān)。在2004-2022年期間飛掠地球的最近距離約為0.1AU,地球和月球的J2項(xiàng)對Apophis軌道的影響可完全忽略 (其量級均小于10m)。太陽J2項(xiàng)的影響與飛掠事件無關(guān),在百年尺度上該效應(yīng)可以對Apophis軌道產(chǎn)生公里量級的攝動。

值得注意的是,當(dāng)Apophis的軌道穿越地月系統(tǒng)時(shí),比如2029年Apophis飛掠地球的最近距離約為5個地球半徑,地球和月球的歷表誤差和非球形引力場的影響將有量級上的變化,需要進(jìn)一步研究太陽系歷表誤差對Apophis飛掠前后運(yùn)動狀態(tài)的影響。

5 數(shù)據(jù)分析與定軌結(jié)果

在分析光學(xué)和雷達(dá)觀測資料時(shí),我們采用了經(jīng)典最小二乘原理下的批處理方法[15,16]。觀測方程和法方程的建立步驟參閱文獻(xiàn)[17]?？紤]到部分?jǐn)?shù)據(jù)存在由不同原因引起的系統(tǒng)性偏差,我們采用了如下的定權(quán)方法:首先將全部觀測數(shù)據(jù)根據(jù)臺站進(jìn)行分組,計(jì)算每組擬合后的均值a0和均方根σRMS;然后將作為該組每個觀測量的測量中誤差;考慮到觀測量N＜10的分組,我們將作為相應(yīng)觀測量的測量中誤差。雷達(dá)觀測資料的權(quán)重由觀測文件中的測量誤差信息給出?？紤]到觀測數(shù)據(jù)中由不同原因引起的粗差,本文采用了3σ0的原則對粗差進(jìn)行剔除,具體剔除數(shù)量見表1。

初始軌道在歷元T0=JD2454466.5處的狀態(tài)向量由JPL#199提供的軌道根數(shù)給出。由于初始軌道較為準(zhǔn)確,通過3次迭代便得到收斂數(shù)值解(驗(yàn)前和驗(yàn)后中誤差的比值趨于1)。解算參數(shù)包括Apophis的軌道參數(shù)和Yarkovsky效應(yīng)參數(shù)A2。我們采用2004-2021年期間的所有光學(xué)和雷達(dá)數(shù)據(jù)得到解 (Sol-A)。圖4和圖5分別給出了部分光學(xué)觀測 (在赤經(jīng)和赤緯方向上)和雷達(dá)測距數(shù)據(jù)的擬合后觀測殘差分布。部分?jǐn)M合后觀測殘差的加權(quán)均方根誤差 (WRMS)在表1中給出。圖6給出了2004-2022年期間Apophis軌道誤差的傳播情況 (9倍形式誤差或9σ)。

圖4 部分光學(xué)觀測 (赤經(jīng)和赤緯方向)的擬合后殘差分布Fig.4 Distribution of the post-fit residuals for partial optical observations(RA and Dec.)

圖5 無線電雷達(dá)測距的擬合后殘差分布Fig.5 Distributionof the post-fit residuals for the radio radar ranging measurements

圖6 2004-2022年期間小行星 (99942)Apophis的軌道誤差 (9σ)Fig.6 Orbit errors(9σ)during the period 2004-2022 for asteroid(99942)Apophis

表3給出了不同參考文獻(xiàn)的Yarkovsky效應(yīng)參數(shù)A2估計(jì)值。與Pérez-Hernández等人[19]和JPL官方網(wǎng)站[20]給出的結(jié)果相比,Sol-A給出的Yarkovsky效應(yīng)參數(shù)值A(chǔ)2更接近Yarkvosky效應(yīng)的線性模型理論計(jì)算值。值得注意的是,解Sol-A 的形式誤差 (1σ) 與 Pérez-Hernández 等人[19]和JPL官方網(wǎng)站[20]給出的結(jié)果略有不同,這可能是因采用了不同于后兩者的定權(quán)策略而造成的。

表3 不同文獻(xiàn)對小行星 (99942)Apophis的Yarkovsky參數(shù)A2的估計(jì) (單位:10-14AU/d2)Table 3 Estimations of Yarkovsky parameter A2from different references for asteroid(99942)Apophis(Unit:10-14AU/d2)

為了考察2021年期間觀測數(shù)據(jù)對Apophis軌道確定和Yarkovsky效應(yīng)測定的影響,我們利用2004-2020年期間的所有數(shù)據(jù)得到另一組解 (Sol-B)。比較表2中解Sol-A和Sol-B的Yarkovsky效應(yīng)參數(shù)值不難發(fā)現(xiàn)2021年的觀測數(shù)據(jù)對Yarkovsky效應(yīng)約束較強(qiáng),被估計(jì)參數(shù)A2的信噪比有明顯提高。

為了檢驗(yàn)Yarkovsky效應(yīng)參數(shù)A2對數(shù)據(jù)的依賴性。我們將所有光學(xué)和雷達(dá)數(shù)據(jù)分割為四組,每組觀測數(shù)據(jù)的數(shù)量和類型近似等同。利用每組數(shù)據(jù)進(jìn)行單獨(dú)解算Apophis軌道參數(shù)和Yarkovsky效應(yīng)參數(shù)值A(chǔ)2,其中Yarkovsky效應(yīng)參數(shù)值A(chǔ)2的結(jié)果見表4。表4中四組解的加權(quán)平均值與解Sol-A保持一致。與解 Sol-A中給出的形式誤差 (1σ)相比,表4中加權(quán)平均值的WRMS(約9σ)更能反映Yarkovsky效應(yīng)參數(shù)A2的真實(shí)不確定度。相應(yīng)地,我們在圖6中取9σ來評估2004-2022年期間的Apophis軌道誤差。

表4 數(shù)據(jù)分組后,四組不同解中的Yarkovsky參數(shù)A2估計(jì)值 (單位:10-14AU/d2)Table 4 Results of Yarkovsky parameter A2from four solutions with sub-grouped observations(Unit:10-14AU/d2)

6 總結(jié)與討論

文中初步分析了大行星歷表誤差,以及部分歷表動力學(xué)參數(shù)誤差對小行星Apophis軌道的影響,發(fā)現(xiàn)觀測期間2004-2021年的幾次飛掠,大行星歷表誤差對Apophis軌道的影響較小,不會對軌道參數(shù)和Yarkovsky效應(yīng)參數(shù)A2的確定產(chǎn)生顯著影響。大行星歷表誤差的影響程度取決于小行星與大行星的最近飛掠距離。2029年Apophis將以5倍地球半徑的距離飛掠地球,期間大行星歷表誤差的影響將會顯著增強(qiáng),需進(jìn)一步研究。

Yarkovsky效應(yīng)是Apophis軌道確定和預(yù)報(bào)的主要誤差來源。本文利用2004-2021年期間的數(shù)據(jù)同時(shí)對Apophis的軌道和Yarkovsky效應(yīng)參數(shù)A2進(jìn)行確定。通過比較Sol-A和Sol-B兩組解,不難發(fā)現(xiàn)2021年的觀測數(shù)據(jù)顯著地提高了Apophis軌道的確定精度和Yarkovsky效應(yīng)的測量信噪比。考慮到解算結(jié)果對觀測數(shù)據(jù)的依賴性,我們給出了在2004-2022年期間Apophis軌道誤差小于10km(9σ)和Yarkovsky效應(yīng)參數(shù)值A(chǔ)2=-2.73±0.12(對應(yīng)信噪比約為22.8)的結(jié)果。