小行星(99942)Apophis的轨道改进与Yarkovsky效应测定

田伟

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

1 引言

近地小天体 (99942)Apophis于2004年6月19日首先被R.A.Tucker、D.J.Tholen和F.Bernardi发现,因最早期的轨道分析显示该颗小行星在2029年4月13日撞击地球的概率超过了2%,得到了广泛关注[1]。随着对Apophis观测资料的不断积累,尤其是高精度的雷达数据,研究人员发现该颗小行星撞击地球的概率远小于早期预报[2-4]。在2021年3月Apophis再次飞掠地球时,研究人员利用这次飞掠期间的观测资料以及历史资料已基本排除2068年之前Apophis撞击地球的可能性[5]。

类似于Apophis,很多近地小天体会多次飞掠地球。由于受地月系统的影响,天体轨道和自转状态都会受到不同程度的改变。Brozovic等人[4]的研究结果显示2029年Apophis飞掠地球期间,其日心半长径将从0.92AU变化为1.10AU(相应地,其轨道类型将从Aten型变化为Apollo型)。Souchay等人[6]的研究结果表明Apophis在2029年飞掠地球期间,在地月系统的引力作用下,其自转状态也将发生显著改变。因此,在利用实测数据确定Apophis轨道之前,本文将定量分析2004-2022年期间几次飞掠地球过程中太阳系历表的误差和相关动力学参数的不确定度对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的地面光学观测和无线电雷达测量数据。文中采用了截至2021年5月20日该网站公布的所有光学数据,以及2005-2021年期间由美国Arecibo和Goldstone雷达站获取的高精度测距和多普勒频移观测数据[4]。表1给出了部分观测数据的统计信息。

表1 小行星 (99942)Apophis(部分)地面光学 (赤经/赤纬)和雷达 (测距/多普勒频移)观测数据的统计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颗大质量的主带小行星。动力学模型考虑的因素包括:(1)点质量相互作用:EIH方程; (2)非球形引力场:太阳、地球和月球的低阶引力场;(3)主带小行星的摄动效应; (4)Yarkovsky效应对Apophis的摄动影响等。

3.1 点质量相互作用

因中心天体太阳与Apophis的相互作用中相对论效应十分显著,需要考虑点质量相互作用中的一阶后牛顿近似下的相对论效应。太阳、月球、八大行星、冥王星对小行星 (99942)Apophis的引力作用由Einstein-Infeld-Hoffmann(EIH)方程给出[7]:

式中,下标A表示Apophis,GMB和GMC分别是天体B和C的质量与引力常数G的乘积,rAB=rB-rA和rAB=|rAB|分别是天体A和B之间的相对位置矢量和相对距离,rA(或rB)和vA(或vB)分别为对应天体A(或B)的位置和速度矢量。相对位置矢量rAC和rBC的定义类似于rAB,而相对距离rAC和rBC的定义类似于rAB,以上各量均在太阳系质心广域参考系中给出。c=299792458m·s-1为光子在真空中的传播速度。一阶后牛顿参数β=γ=1。

3.2 非球形引力场

除考虑主要天体 (作为点质量)的引力作用外,文中也考虑了地球、月球和太阳的非球形引力场对Apophis轨道的影响。非球形引力场的引力作用取决于Apophis与非球形天体之间距离。通常近地小行星会 (近似)周期性地与地球或金星交会 (见图1)。以地球为例,其引力场高阶项 (J2,Cnm,Snm)对某外部天体的影响可由式 (2)表示[8]

式中,在ξ-η-ζ坐标系中,ξ坐标轴代表从地球到Apophis的指向,ξ-η平面包含了地球自转轴,ζ轴与ξ和η轴形成右手坐标系。GME和R分别代表了地球的质量与引力常数乘积和平均半径。表示n阶m次的缔合Legendre函数。n1和n2分别是Jn和 (Cnm,Snm)的最高阶数。(r,φ,λ)是从地球到Apophis构成的位置向量在地固坐标系中的球面坐标。需要说明的是,在式 (2)的基础上计算太阳系质心坐标系下的相互作用力时,需要实施一系列坐标转换。太阳、地球和月球的非球形引力场分别取J2项、3阶以下项和2阶以下项。限于篇幅,文中未考虑金星非球形引力场对Apophis的影响。

3.3 主带小行星摄动模型

某一特定小行星对Apophis的摄动效应取决于摄动小行星的质量和相对距离。目前被精确测量质量参数的小行星在已发现的小行星总数中占比较小,因而小行星摄动效应是动力学模型中一个难以确定的因素。Giorgini等人[2]统计了约37.3万颗已知小行星对Apophis轨道的摄动效应,并给出了4颗摄动影响最大的小行星:Ceres、Vesta、Pallas和Hygiea,其中Ceres、Pallas和Vesta是主带小行星中质量最大的3颗。文中未对单颗小行星的摄动效应做逐个筛选,而是选取了48颗大质量小行星作为摄动小天体。

3.4 Yarkovsky效应

由于目前对Apophis的观测主要依赖于地面设备,其形状、自转轴指向、平均密度、反照率、热惯量等基本几何和物理参数确定精度偏低,在计算Apophis受到的Yarkovsky效应时,文中采用了简化的非保守力线性模型[3,9-10],

式中,标准参数A1、A2和A3分别代表与径向、横切线方向和轨道面法线方向加速度分量相关的参数,r和r分别为小行星相对太阳的位置向量和距离,h=r×r˙。文中仅考虑了横切线方向的作用力,即A1=A3=0,A2≠0。线性模型中的非零参数A2将与轨道参数一起解算。

4 部分太阳系历表参数对Apophis轨道的影响

作为太阳系内主要天体,太阳、地球和月球的质量不确定度对Apophis轨道有着不同程度的影响。目前太阳、地球和月球的质量参数可以通过不同技术手段精确确定[11-14],其数值和不确定度见表2。图2给出了相应天体质量不确定度对Apophis轨道的摄动效应。

表2 太阳、地球和月球的质量参数 (单位:km3/s2)Table 2 The mass parameters of the sun,Earth and Moon(Unit:km3/s2)

图2 由太阳、地球和月球的质量参数不确定度引起的小行星 (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轨道产生公里量级的摄动。

值得注意的是,当Apophis的轨道穿越地月系统时,比如2029年Apophis飞掠地球的最近距离约为5个地球半径,地球和月球的历表误差和非球形引力场的影响将有量级上的变化,需要进一步研究太阳系历表误差对Apophis飞掠前后运动状态的影响。

5 数据分析与定轨结果

在分析光学和雷达观测资料时,我们采用了经典最小二乘原理下的批处理方法[15,16]。观测方程和法方程的建立步骤参阅文献[17]。考虑到部分数据存在由不同原因引起的系统性偏差,我们采用了如下的定权方法:首先将全部观测数据根据台站进行分组,计算每组拟合后的均值a0和均方根σRMS;然后将作为该组每个观测量的测量中误差;考虑到观测量N＜10的分组,我们将作为相应观测量的测量中误差。雷达观测资料的权重由观测文件中的测量误差信息给出。考虑到观测数据中由不同原因引起的粗差,本文采用了3σ0的原则对粗差进行剔除,具体剔除数量见表1。

初始轨道在历元T0=JD2454466.5处的状态向量由JPL#199提供的轨道根数给出。由于初始轨道较为准确,通过3次迭代便得到收敛数值解(验前和验后中误差的比值趋于1)。解算参数包括Apophis的轨道参数和Yarkovsky效应参数A2。我们采用2004-2021年期间的所有光学和雷达数据得到解 (Sol-A)。图4和图5分别给出了部分光学观测 (在赤经和赤纬方向上)和雷达测距数据的拟合后观测残差分布。部分拟合后观测残差的加权均方根误差 (WRMS)在表1中给出。图6给出了2004-2022年期间Apophis轨道误差的传播情况 (9倍形式误差或9σ)。

图4 部分光学观测 (赤经和赤纬方向)的拟合后残差分布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效应参数A2估计值。与Pérez-Hernández等人[19]和JPL官方网站[20]给出的结果相比,Sol-A给出的Yarkovsky效应参数值A2更接近Yarkvosky效应的线性模型理论计算值。值得注意的是,解Sol-A 的形式误差 (1σ) 与 Pérez-Hernández 等人[19]和JPL官方网站[20]给出的结果略有不同,这可能是因采用了不同于后两者的定权策略而造成的。

表3 不同文献对小行星 (99942)Apophis的Yarkovsky参数A2的估计 (单位:10-14AU/d2)Table 3 Estimations of Yarkovsky parameter A2from different references for asteroid(99942)Apophis(Unit:10-14AU/d2)

为了考察2021年期间观测数据对Apophis轨道确定和Yarkovsky效应测定的影响,我们利用2004-2020年期间的所有数据得到另一组解 (Sol-B)。比较表2中解Sol-A和Sol-B的Yarkovsky效应参数值不难发现2021年的观测数据对Yarkovsky效应约束较强,被估计参数A2的信噪比有明显提高。

为了检验Yarkovsky效应参数A2对数据的依赖性。我们将所有光学和雷达数据分割为四组,每组观测数据的数量和类型近似等同。利用每组数据进行单独解算Apophis轨道参数和Yarkovsky效应参数值A2,其中Yarkovsky效应参数值A2的结果见表4。表4中四组解的加权平均值与解Sol-A保持一致。与解 Sol-A中给出的形式误差 (1σ)相比,表4中加权平均值的WRMS(约9σ)更能反映Yarkovsky效应参数A2的真实不确定度。相应地,我们在图6中取9σ来评估2004-2022年期间的Apophis轨道误差。

表4 数据分组后,四组不同解中的Yarkovsky参数A2估计值 (单位:10-14AU/d2)Table 4 Results of Yarkovsky parameter A2from four solutions with sub-grouped observations(Unit:10-14AU/d2)

6 总结与讨论

文中初步分析了大行星历表误差,以及部分历表动力学参数误差对小行星Apophis轨道的影响,发现观测期间2004-2021年的几次飞掠,大行星历表误差对Apophis轨道的影响较小,不会对轨道参数和Yarkovsky效应参数A2的确定产生显著影响。大行星历表误差的影响程度取决于小行星与大行星的最近飞掠距离。2029年Apophis将以5倍地球半径的距离飞掠地球,期间大行星历表误差的影响将会显著增强,需进一步研究。

Yarkovsky效应是Apophis轨道确定和预报的主要误差来源。本文利用2004-2021年期间的数据同时对Apophis的轨道和Yarkovsky效应参数A2进行确定。通过比较Sol-A和Sol-B两组解,不难发现2021年的观测数据显著地提高了Apophis轨道的确定精度和Yarkovsky效应的测量信噪比。考虑到解算结果对观测数据的依赖性,我们给出了在2004-2022年期间Apophis轨道误差小于10km(9σ)和Yarkovsky效应参数值A2=-2.73±0.12(对应信噪比约为22.8)的结果。