田偉

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

1 引言

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

類似于Apophis,很多近地小天體會多次飛掠地球。由于受地月系統(tǒng)的影響,天體軌道和自轉狀態(tài)都會受到不同程度的改變。Brozovic等人[4]的研究結果顯示2029年Apophis飛掠地球期間,其日心半長徑將從0.92AU變化為1.10AU(相應地,其軌道類型將從Aten型變化為Apollo型)。Souchay等人[6]的研究結果表明Apophis在2029年飛掠地球期間,在地月系統(tǒng)的引力作用下,其自轉狀態(tài)也將發(fā)生顯著改變。因此,在利用實測數(shù)據(jù)確定Apophis軌道之前,本文將定量分析2004-2022年期間幾次飛掠地球過程中太陽系歷表的誤差和相關動力學參數(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)收集了關于Apophis的地面光學觀測和無線電雷達測量數(shù)據(jù)。文中采用了截至2021年5月20日該網(wǎng)站公布的所有光學數(shù)據(jù),以及2005-2021年期間由美國Arecibo和Goldstone雷達站獲取的高精度測距和多普勒頻移觀測數(shù)據(jù)[4]。表1給出了部分觀測數(shù)據(jù)的統(tǒng)計信息。

表1 小行星 (99942)Apophis(部分)地面光學 (赤經(jīng)/赤緯)和雷達 (測距/多普勒頻移)觀測數(shù)據(jù)的統(tǒng)計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 動力學模型

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

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

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

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

3.2 非球形引力場

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

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

3.3 主帶小行星攝動模型

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

3.4 Yarkovsky效應

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

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

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

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

表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軌道確定的一個誤差項。其影響程度取決于最近飛掠距離,圖3給出了2004-2022年期間,地球的歷表位置誤差 (假設位置誤差為1km)對Apophis軌道的攝動效應。相比于與地球歷表誤差的攝動效應,月球的影響較小。

圖3 地球歷表誤差對Apophis軌道的影響 (假設位置誤差為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項是主要影響項。與被飛掠天體的歷表位置誤差類似,其非球形引力場的引力攝動效應與最近飛掠距離直接相關。在2004-2022年期間飛掠地球的最近距離約為0.1AU,地球和月球的J2項對Apophis軌道的影響可完全忽略 (其量級均小于10m)。太陽J2項的影響與飛掠事件無關,在百年尺度上該效應可以對Apophis軌道產(chǎn)生公里量級的攝動。

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

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

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

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

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

圖5 無線電雷達測距的擬合后殘差分布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給出了不同參考文獻的Yarkovsky效應參數(shù)A2估計值。與Pérez-Hernández等人[19]和JPL官方網(wǎng)站[20]給出的結果相比,Sol-A給出的Yarkovsky效應參數(shù)值A2更接近Yarkvosky效應的線性模型理論計算值。值得注意的是,解Sol-A 的形式誤差 (1σ) 與 Pérez-Hernández 等人[19]和JPL官方網(wǎng)站[20]給出的結果略有不同,這可能是因采用了不同于后兩者的定權策略而造成的。

表3 不同文獻對小行星 (99942)Apophis的Yarkovsky參數(shù)A2的估計 (單位: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效應測定的影響,我們利用2004-2020年期間的所有數(shù)據(jù)得到另一組解 (Sol-B)。比較表2中解Sol-A和Sol-B的Yarkovsky效應參數(shù)值不難發(fā)現(xiàn)2021年的觀測數(shù)據(jù)對Yarkovsky效應約束較強,被估計參數(shù)A2的信噪比有明顯提高。

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

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

6 總結與討論

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

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