引　言

近日点距离$q \leqslant 1.3\; \rm AU$的小天体被称为近地天体，其中大多数为日心轨道上的近地小行星。近地小行星与地球轨道发生相交或相切会给地球带来潜在风险，其撞击地球的过程可能会释放数十Mt（直径50 m）至数百万Mt（直径数km）的能量，对人类生命财产安全构成威胁。直径在140 m以上与地球最小轨道距离在750万 km以内的小行星被定义为潜在威胁小行星（Potentially Hazardous Asteroid，PHA），总数估计在几万颗以上，然而绝大多数还没有被发现。这些没有被发现的潜在威胁小行星随时可能突然出现并与地球发生碰撞，因此需要对其监测预警并建立防御系统，以应对近地小行星撞击的威胁。国际上早已开展关于近地小行星监测预警及防御的相关研究，并成立国际组织合作建立近地小行星的全球监测网，防范可能来自太空的威胁。

近地小行星撞击监测作为一个重要的研究领域，碰撞概率（Impact Probability）计算是其中的关键因素和前提，在20多年的时间里被逐步建立起来并不断改进和完善^[1]，以提高碰撞概率的计算精度，并在此基础上产生了多个评估系统。意大利比萨大学（University of Pisa）和空间动力学服务有限公司（SpaceDyS）研发了CLOMON和CLOMON2系统^[2]。CLOMON2在近地小行星动力学网站（Near Earth Objects-Dynamic Site，NEODys） ¹在线发布了潜在威胁小行星在前后100年与地球发生密近交会的时间、距离及概率。美国国家航空航天局（National Aeronautics and Space Administration，NASA）的喷气推进实验室（Jet Propulsion Laboratory，JPL）开发了SENTRY和SENTRY-Ⅱ系统，该系统在近地天体研究中心（Center for Near Earth Object Studies，CNEOS）网站 ²上给出了未来发生碰撞的概率和时间，以及后续发生密近交会的时间、距离，时间精度在1 min，距离精确到10^–7 AU。欧洲航天局（European Space Agency，ESA）AstOD系统的结果发表在近地天体协调中心（Near-Earth Objects Coordination Centre，NEOCC）网站 ³，给出了与地球发生密近交会的时间和距离，及发生碰撞的时间和概率，精度与SENTRY-Ⅱ相当。此外，对于临近小行星进行碰撞监测的系统也陆续被开发出来，如JPL的SCOUT^[3]、赫尔辛基大学（University of Helsinki）的NEORANGER^[4]和比萨大学/SpaceDyS公司的NEOScan^[5]。

近年来，中国也逐步展开了小行星监测预警、在轨处置等方面的研究。2021年航天日上提出了将论证建设近地小行星防御系统，2022年航天日上提出要完善建立地基天基对小行星的监测预警系统，组织编制近地小行星防御发展规划，开发近地小行星防御仿真推演软件并组织开展基本流程推演。

本文将对小行星碰撞概率计算方法研究的发展历程进行简要综述，同时提出对未来该研究方向的关键问题和未来发展趋势的思考。

1   碰撞概率理论基础

1.1   碰撞概率理论基础

直接通过小行星标称轨道来计算其与地球的碰撞是不可靠的。由于小行星的轨道根数存在误差，可通过小行星轨道根数的误差协方差矩阵确定未来小行星所处的置信区域，评估其内部是否包含会与地球发生撞击的小行星轨道簇，从而给出与地球的碰撞概率。Öpik等^[6-8]在20世纪50年代逐步建立了近地天体平均和长期碰撞概率的理论基础。Öpik^[9]提出了第1个关于行星相遇的数学理论：假设小行星的运动可被处理为不同两体问题解的组成，小行星相对于行星的相对速度定义了行星双曲轨道入射渐近线的方向和速度，这个方向和速度是小行星日心轨道的半长径 a 、偏心率 e 和倾角 i 的简单函数，其中忽略了交会距离的项。密近交会时刻被计算为速度矢量在双曲线轨道输出渐近线方向的瞬时偏转，忽略了太阳扰动以及小行星沿着连接2个渐近线的弯曲路径行进实际花费的时间。根据Öpik的理论，目标的轨道被认为是日心椭圆轨道，直到它进入由行星引力为主导的区域，此时轨迹变为双曲线上的一支，当其离开行星影响范围后进入一个新的日心椭圆轨道，其初始条件由双曲三体问题的解给出。Öpik的理论仅针对小行星和目标行星2个轨道发生实际接触的情况，因此理论仅在最小接近距离为0时才是精确的。此外，理论没有考虑到小行星随后与同一颗行星（或另一颗）的相遇并不独立于先前事件的发生。这一理论现在被称为共振回归，自20世纪70年代以来就被用于航天器导航，但直到20世纪90年代末才首次应用于小行星密近交会的研究^[10-11]。针对以上Öpik理论2个方面的不足，Valsecchi等^[12]进行了发展，使得该理论可用于非0近距离交会事件的计算。

1.2   目标平面和坐标系

Kizner^[13]在太阳系内星际航行中首次提出了小行星碰撞问题研究的主要工具—b平面，也称为目标平面。每当小行星和地球之间有近距离接触时，b平面即定义为以地球为原点且垂直于小行星关于地球相对速度矢量的平面，从而可通过轨道渐近线的位置b决定小行星与地球之间是否发生碰撞。实际应用时通常采用图1所示的近似：垂直于近地小行星吻切轨道地心双曲线的渐近线方向（又称无摄相对速度方向）的平面。在密近交会中当小行星的相对速度大于地球的逃逸速度时都可以采用此近似方法。以地心为原点在b平面建立右旋坐标系$(\xi ,\eta ,\zeta )$。图1中，η为无摄相对速度方向${v_{\mathrm{\infty }}}$，即渐近线方向，垂直于b平面； $\xi$为地球日心速度${v}_{\oplus}$在b平面上投影的反方向。令渐近线与b平面交点的坐标为$(\xi ,\zeta )$，其有明确的物理意义，其中$\xi$为近地小行星和地球的离迹交会距离，如图2所示 ，通过调整交会时间，可得到 $\xi$ 坐标在b平面的最小值，即小行星和地球吻切轨道的最小距离，称为地球最小的轨道距离（ Minimum Orbit Intersection Distance ，MOID）^[14]。$\zeta$为小行星和地球的沿迹交会距离，表示近地小行星提早或滞后于交会点的距离，间隔时间为

图  1  b平面示意图^[15]

Fig.  1  b-Plane^[15]

下载: 全尺寸图片 幻灯片

图  2  b平面坐标系^[14]

注：实线表示小行星轨道的渐近线，虚线表示实际受到引力影响的轨道；实线表示地球真实半径，虚线表示受引力影响，当小行星进入R_corss的范围后便会被吸引飞向地球；B表示小行星到b平面坐标系原点的距离。

Fig.  2  b-plane coordinates^[14]

下载: 全尺寸图片 幻灯片

其中：$\theta$为图1中${v}_{\infty }$和${v}_{\oplus}$ 的夹角。

其中：${r}_\oplus$为地球半径，当交点到地心的距离$b < {b}_{\oplus}$时表示撞击发生；${{\rm GM}}_{\oplus}$为地球引力常数。

2   碰撞概率计算方法

2.1   线性碰撞概率计算方法

1993年，Muinonen等^[16]发展了计算小行星碰撞概率的线性方法，将地球与小行星轨道根数误差椭球间的距离用于评估碰撞概率，于1994年提出了潜在威胁小行星的概念，把与地球最小轨道交会距离小于0.05 AU、绝对星等小于22 mag（magnitude，星等）的小行星定义为潜在威胁小行星，确定了主要研究目标^[17]。

JPL于90年代就开始研究，1994年Chodas等^[18]采用线性方法预报了小行星、彗星与地球的密近交会事件，计算了轨道的不确定度和碰撞概率。根据小行星初始历元状态矢量的误差协方差，通过线性过程求出目标时刻的误差协方差从而获得置信椭球体，置信椭球体与地球相交截面的概率密度即可用于刻画小行星的碰撞概率。

小行星轨道根数的不确定性分析采用最小二乘法来描述，在不确定性分析中，只关注误差协方差矩阵P。误差协方差矩阵取决于测量误差统计量和观测参数的偏导数，在进行误差协方差矩阵的状态转移计算P(t) =Φ(t_k)P(t_k)Φ(t_k)^T时，可采用序贯处理方法和伪序贯处理方法求取P(t_k)。其中

其中：e₁, e₂, ···, e₆是历时状态的分量；x₁, x₂, ···, x₆是观测时的状态分量。Chodas^[19]发现伪序贯处理方法比序贯处理方法更为稳健。在加入新的观测结果获得新的误差协方差矩阵后，通过状态转移将误差协方差矩阵映射到碰撞时刻。之后将由误差协方差矩阵得到的置信椭球体进行坐标转换，并投影到目标平面，得到置信椭圆。此时矩阵P(t)变换成C(t)，C(t)实际就是P(t)左上角的（2 × 2）子矩阵，使C(t)对角化可得

其中：$\sigma _{1}$、${\sigma }_{2}$分别为椭圆的半长轴和半短轴。椭圆与地球在目标平面截面内二维高斯密度函数的积分就是碰撞概率为

其中：$\left({\xi }_{\mathrm{I}}，{\zeta }_{\mathrm{I}}\right)$为椭圆中心的坐标。积分可通过计算机或采用Monte Carlo统计方法近似计算。大多数情况下椭圆区域远大于地球截面，截面面积事实上就是地球的横截面积。因此，碰撞概率可以近似考虑为以地球为中心的正态密度函数与地球截面面积的乘积，从而简化计算过程。

1998年，Marsden^[10]采用线性方法对近地小行星1997 XF11未来与地球发生碰撞的可能性进行了评估，认为将于2028年10月26日相撞，随后基于1990—1997年的累计观测数据进行重新计算排除了这一可能性。但研究发现，1997 XF11与地球在2028年发生密近交会后，轨道可能会受影响从而呈现非线性特征，使得线性的碰撞概率计算方法无法有效预测未来的碰撞概率。关于1997XF11碰撞概率的研究推动了撞击危险评估理论和算法的发展，出现了非线性碰撞概率的计算方法，同时提出了自动化碰撞监测预警系统的概念。

2.2   非线性碰撞概率计算方法

线性方法存在局限性，当初始状态矢量的置信椭球过于细长，或外推过程中出现非线性影响时，状态矢量映射到目标平面上的过程就会出现非线性化，导致初始历元状态矢量置信椭球体本身拉长或在目标平面上的投影椭圆非常狭长，从而引起线性方法的失效。另一方面，观测弧段过短、轨道误差较大的目标在计算时的过大误差，也会导致无法使用线性方法。由此产生了对非线性方法计算碰撞概率的研究。

2.2.1   Monte Carlo方法

1996年，Chodas等^[20]利用非线性方法研究了Shoemaker-Levy彗星轨道的演化历史，是非线性方法的最早应用，在轨道元素空间生成1 000个随机点填充6维不确定性的椭球体，以获得与轨道解的协方差矩阵一致的初始条件集合，计算发现彗星可能在1929 ± 9年被木星俘获，同时确定该彗星碰撞的可能性高达95%，并给出了碎片撞击的时间。1999年Chodas等^[21]用Monte Carlo方法取样观测历元轨道根数的线性6维置信区间，然后采用非线性方法数值积分到碰撞时刻以计算碰撞概率，通过将初始历元置信椭球体划分为若干个子区域，对子区域进行线性化处理。Monte Carlo方法的主要思想是利用频次来估计概率，方法直接采用最小二乘法原理，对轨道元素空间中的概率分布进行采样，方法可对置信区域（Confidence Region）充分采样，并且计算思路简单，能对非线性方法进行完整的体现，但样本数量过大会带来计算成本的负担。

2.2.2   变化线（ Line Of Variation， LOV）方法

由于轨道的不确定性过大，需要考虑置信区域中一组残差在固定阈值范围内的轨道，通过抽样获得虚拟小行星（Virtual Asteroid）的轨道^[22]，判断它们在轨道外推时是否会与地球发生撞击。发生碰撞的虚拟小行星被称为虚拟碰撞体（Virtual Impactor），一旦确定发生碰撞的虚拟碰撞体，即获得了导致碰撞的初始条件，可以进一步计算碰撞概率。但虚拟碰撞体的碰撞概率与虚拟碰撞体在轨道元素空间中的体积成比例，如果碰撞概率较小，那么采样就会变得非常密集。因此需平衡搜索的完备性和计算成本。

意大利比萨大学Milani等^[11]采用多重解方法研究了小行星1999 AN10，首次证认了碰撞解和观测数据之间的相关性。同时还提出了一维抽样的变化线（ Line Of Variation， LOV）方法：在轨道根数空间中沿着一条可微的曲线进行采样{x(σ_i)}, i = –M, ···, M，对应于LOV的参数{σ_i}, i = –M, ···, M，在某些情况下该曲线可以代表置信区域的主轴，外推时轨道的不确定性沿着轨道延伸。沿着曲线计算对应于协方差矩阵最大特征值的单位特征向量λ₁ = (σ₁) ²，在特征空间中以h = σ/p步长前进，其中p为整数。由于是在曲线上采点，点x的下一步x + h × σ₁V₁可能不再属于该曲线，但相距不远，需要对其进行微分修正^[23]。在LOV变化线上进行采样的实质就是获得虚拟小行星集，也就是计算多个解，多解法在观测历元的轨道根数非线性置信区域的变量中心轴上取样，然后数值积分到碰撞时刻来估计碰撞概率^[23-24]。在对虚拟小行星外推时，目标是找到最小的交会距离。以交会距离的平方r²为LOV参数σ的函数。r²(σ) 是可微的，可以搜索使f(σ) = $\dfrac{{\rm{d}}\left({r}^{2}\right)}{{\rm{d}}\sigma }$为0的项。如果有闭合区间[σ₁, σ₂]，其中的σ对于所有初始条件都会发生密近交会，并且导数在2个极值的值为f(σ₁)<0，f(σ₂)>0，则至少有1个r²(σ)的最小值在区间内。Milani等^[25]对目标平面进行了修正，定义垂直于最接近时小行星的地心速度矢量的平面为目标平面。由于在目标平面上的分析是局部的，因此可在每个虚拟小行星的附近进行线性近似，通过对误差协方差矩阵的计算获得投影椭圆的半长轴和半短轴，即LOV上轨迹的拉伸长度和宽度。令$\bar{y}$ = f($\bar{x}$)是LOV轨道$\bar{x}$在目标平面上的轨迹。偏导数Df将置信椭球体Z_X 映射到目标平面上的置信椭圆Z_Y 上。利用协方差进行传播，定义椭圆Z_Y

其中

式（7）为目标平面上（2 × 2）的矩阵，其特征值的平方根是椭圆半长轴和半长轴的长度，也就是LOV轨道处拉伸的长度和宽度。

在检测每个虚拟小行星未来与地球的碰撞事件时，有些虚拟小行星完全不会发生碰撞事件，有些则可能发生多次。Milani等^[26]首先将密近交会的轨道按时间分类，然后每一类被进一步划分为连续的LOV段，通过分析LOV段搜索虚拟碰撞体。当给定的LOV段有多个点位于目标平面时，那么此LOV段的拉伸长度较短，可采用局部的线性方法进一步计算。相反，在强非线性情况下，LOV段拉伸很长，并且两点之间的变化很快，在这种情况下需要在每个虚拟小行星的邻域中进行局部分析。

由于虚拟小行星是一组对平滑曲线采样获得的点，Milani等^[26]和Tommei^[27]假设2个连续的虚拟小行星在目标平面上的轨迹连线穿越地球截面，如果轨迹的几何形状简单，则可采用插值方法获得位于地球截面内的点，即假设有2个连续的虚拟小行星 x_i和x_i+1在目标平面上的轨迹点y_i和y_i+1的连线穿越地球截面，则插值方法可提供 x_i+δ，0 < δ < 1，使y_i+δ在地球撞击截面内，而x_i+δ附近有一个虚拟碰撞体，通过计算以该点为中心的概率密度函数可得到碰撞概率。

Milani等^[24-26]根据分析，计算2004 FU4在2010年的碰撞概率为4 × 10^–8；1997 AE12与地球的密近交会距离0.09 AU，随着进一步观测之后排除了对地球的威胁；同时计算了1997 XF11在2028年的碰撞概率，发现可能性极小。Tommei^[27]讨论了小行星Apophis在2029年和2036年的碰撞概率，其中2029年的碰撞概率高达2%。

LOV变化线方法较为复杂，但相对Monte Carlo方法减少了计算量。在采样方式上，LOV变化线上的概率密度是由高斯方法定义，但LOV曲线中间的概率密度相比LOV两端要高得多，采用固定步长进行采样可能导致虚拟小行星的遗漏，因此采取与概率密度成反比的步长进行采样是更优的选择^[28]，但该采样方式也会导致计算成本的增加。当置信区形状不稳定，出现多个方向时，一维采样方式可能会出现遗漏虚拟小行星的情况。同时小行星在与地球密近交会后可能导致置信椭球体映射过程中的强非线性，从而导致置信区域拉伸较大且变化迅速或LOV采样点处的拉伸非常大，使得LOV变化线变得复杂，带来虚拟碰撞体搜索失败的问题。此时需适当地增加采样密度，使得给定的LOV段在目标平面上有多个点，从而进行连续虚拟小行星的插值计算^[29]。

3   碰撞概率计算方法的改进

随着对小行星探测技术的提升和观测数据的增加，发现了越来越多对地球具有潜在危险的小行星。而小行星与地球发生碰撞事件的场景常常是较为复杂的，仅仅依靠经典的线性和非线性碰撞概率计算方法，无法满足对碰撞概率高精度快速预报的要求，例如小行星飞掠地球的速度过快或者过慢，或已经发生过密近交会的小行星再次飞掠地球，都会使经典碰撞概率计算方法的结果产生偏差。因此根据实际需求，对碰撞概率计算方法进行改善的研究不断涌现。

3.1   变化流形（ Manifold Of Variations ，MOV）方法

巡天观测的发展使得小行星的观测数据大大增加，但同时巡天观测获得的弧段过短，给小行星轨道确定和危险评估带来了很大的问题^[30]。多解法是在置信区域内进行采样，当观测到的弧段足够长时，置信区域的形状是拉长状的，因此可采用LOV方式以一维曲线作为该区域的主轴进行近似处理。但当观测到的弧段很短时，置信区域变成平坦的盘面，在单一方向上不再具有明显优势，而是在2个方向上都有一定的宽度，导致一维LOV变化线无法完全代表置信区域的分布，从而无法获得有代表性的采样点。

对于由于弧段太短而无法采用传统方法确定完整轨道的观测弧段，称之为极短弧（Too Short Arc，TSA）^[31]。Milani等^[31]通过线性回归和其它拟合程序来获得极短弧的直线弧段，将极短弧由可归属项（Attributable）来表示，可归属项由参考时间（观测时间的平均值）下的2个平均角坐标和2个相应的角速率构成。令 r 和 q 分别表示目标和观测站为时刻t在地球上的日心位置向量。利用少量的观测数据计算赤经α、赤纬δ，以及其导数$\dot{\alpha }$和$\dot{\delta }$，通过多项式拟合可得到可归属项（α，δ，$\dot{\alpha }$，$\dot{\delta }$）^[32]。可归属项中包含的信息无法求得目标相对太阳的距离r和径向速度$\dot{r}$，由此引入容许区域 （Admissible Region）的概念，在受约束的范围内随机生成距离r和径向速度$\dot{r}$，与可归属项共同计算小行星的轨道根数。约束条件包括日心两体能量、地心两体能量、地球影响范围半径和地球半径。可表示为

其中：A = [α, δ, $\dot{\alpha }$, $\dot{\delta }$]，B = [r, $\dot{r}$]。由于与通常情况不同，（6 × 6）的协方差矩阵不可用，不确定性需要描述为

其中：$\mu > 0$是一个参数；A₀为可归属项角坐标的标称值；C_A0是对应的正规矩阵。此集合不是笛卡尔乘积，尽管在许多情况下它可以近似为A空间中置信椭球与由A₀计算得到的容许区域$\mathcal{D} (A)$的乘积

由于容许区域是紧致的，可以用有限数量的点对其进行采样^[5]。如果存在标称解，则在标称解的附近进行蛛网采样^[27]，否则采用（r, $\dot{r}$）上的矩形网格在容许区域进行扫描^[33]。

蛛网采样利用最小化观测残差RMS目标函数二次近似的水平曲线来获得。目标函数的水平曲线是6维轨道根数空间中的同心5维椭球体，因此（r, $\dot{r}$）空间上的水平曲线由5维的边缘椭球表示，然后对该曲线进行采样。

矩形网格的采样分为两步，首先根据日心能量和容许区域中连接部分的数量划分采样区域并进行采样；然后根据已有的网格应用双约束微分校正，获得完整的轨道样本，并找到（r, $\dot{r}$）中的最小值和最大值，以此为约束在新的网格中进一步采样。

极短弧误差协方差矩阵的特征值往往存在有2个特征值比其它特征值大得多的情况，表明不确定性存在于2个方向上，因此置信区域具有二维结构^[1]。可以在容许区域上定义二维流形，即变化流形MOV，进行二维采样，并采用条件概率的方法将残差空间的概率密度函数非线性传播到采样空间。采用蛛网采样的方式，利用最小化观测残差函数曲线来构建覆盖平面（r, $\dot{r}$）内置信区域的蛛网。在每条水平曲线上选择与固定方向相对应的点作为虚拟小行星，在椭圆极坐标系（R, θ）中建立正则网格，$0 \leqslant \theta < 2 \text{π}$，$0 \leqslant R \leqslant M_{\rm RMS}$ （$M_{\rm RMS}$是定义RMS最大值的参数）。根据初始轨道的协方差矩阵和轨道本身，对网格的每个点进行变换。

1）将（R, θ）空间的网格映射到距离速度平面

其中：λ₁和λ₂是变量（r, $\dot{r}$）的（2 × 2）协方差矩阵的特征值；V₁是对应于较大特征值λ₁的特征向量。

2）由于标称解在距离速度平面的原点，需要利用（r_norm, $\dot{r}$_norm）向量对所有点进行移动，表示用最小二乘拟合计算出的初始轨道的真实位置

随后通过微分修正的方法进行后续校正。

（1）由4维向量A₀和（r, $\dot{r}$）平面上的（r₀, $\dot{r}$₀）构成6维向量；

（2）考虑固定方向上的下一个点（r₁, $\dot{r}$₁），构成（A₀, r₁, $\dot{r}$₁）；

（3）对A₀进行微分修正得到${{\boldsymbol{A}}}_{1}^{h}$，从而得到（${{\boldsymbol{A}}}_{1}^{h}$, r₁, $\dot{r}$₁）；

（4）重复前面的步骤，分析沿着h方向的所有点。

Vigna^[34]对2014 AA进行了计算，根据有3个观测数据的短弧发现其有3%的碰撞概率，结合后续7个观测数据，其碰撞概率提升为100%。

3.2   轨道测距方法 （Ranging Methods）

对于在即将撞击前几天或几小时才检测到的小行星，经典碰撞概率算法不再适用。除了可采用容许区域的方法，还可利用轨道测距方法来计算，测距方法主要分为2类，轨道统计测距（Statistical Methods）和轨道系统测距（Systematic Methods）。

Muinonen等^[34]采用轨道统计测距方法严格映射轨道根数的概率密度，预测了小行星1998 OX4的碰撞概率。在Virtanen等^[35]利用轨道统计测距技术计算轨道根数概率密度的基础上，Muinonen等不再从整个观测集中随机选择两对观测值，而是使用一对观测值，对赤经和赤纬随机采样生成样本，同时样本内观测值和计算值之间的残差在定义阈值内，并且改进了迭代过程，减少了计算成本。Muinonen等^[34]首先生成大量无偏的轨道根数，将其分布作为映射轨道根数的概率密度；其次在每个轨道根数集附近，沿着最小二乘协方差矩阵的主特征向量计算一维区间；然后在每个一维区间中计算在给定时间内导致碰撞的根数的子区间，由此可以根据导致碰撞的子区间轨道根数的概率密度计算小行星的碰撞概率。则小行星的碰撞概率为

其中：N₁为初始轨道区间${\Delta }_{i}$的数目；N₂为每个轨道区间的样本轨道数；${\Delta }_{i}^{c}$表示发生碰撞的轨道区间i；$\,\boldsymbol{\varLambda }$为观测误差的协方差矩阵；$\Delta \boldsymbol{\psi }\left(P\right)$为天空平面的残差；${P}_{ij}$ 和${P}_{ij}^{c}$分别为轨道区间和碰撞轨道区间中的样本轨道j。通常碰撞概率可以简化计算为

Muinonen等^[34]据此给出了1998 OX4在2014、2038、2044、2046年的碰撞概率分别为5 × 10^–7、2 × 10^–5、9 × 10^–6、4 × 10^–5。

Oszkiewicz等^[36]利用Markov链根据 Bayesian反演理论得出的概率密度分布生成无偏轨道样本序列，代替随机采样的方式，改进了统计测距方法，该方法被称为MCMC（Markov-Chain Monte Carlo ）方法。从观测值中选择2对（α_A, δ_A）、（α_B, δ_B），通过对2个球面位置Q = {ρ_A, α_A, δ_A；ρ_B, α_B, δ_B}引入高斯候选密度（Gaussian proposal densitiy） p_t（${Q'}$；${Q}_{{\rm{t}}}$）进行轨道根数空间的MCMC映射，分别用${Q'}$和${Q}_{{\rm{t}}}$表示候选和最后一个被接受的球面位置，利用该密度来产生连续的轨道元素链（P⁰, P¹, ···, P^N）。由此产生候选距离$\rho' _{{\rm{A}}}、\rho' _{{\rm{B}}}$和位置（$\alpha' _{{\rm{A}}}$, $\delta' _{{\rm{A}}}$）、（$\alpha '_{{\rm{B}}}$, $\delta' _{{\rm{B}}}$），根据坐标可以得到一条轨道${P'}$，检测所有观测值的拟合，判断${P'}$是否被接受，将高斯候选密度转移至该轨道对应的球坐标，重复此过程直到密度平稳。新的候选轨道${P'}$根据候选轨道根数的概率密度、最后被接受的样本p_p(P_t)的概率密度以及轨道根数的密度被接受或拒绝。可以计算$a_{\rm{r}}=\dfrac{{p}_{{\rm{p}}}\left({P'}\right){J}_{{\rm{t}}}}{{p}_{{\rm{p}}}\left({P}_{{\rm{t}}}\right){J'}}$，其中${J'}$和${J}_{{\rm{t}}}$是候选样本和最后接受样本的轨道根数的雅可比行列式。当$a_{\rm{r}} \geqslant 1$时接受所候选根数。随后Muinonen等^[37]研究了随机游走测距的方法，其中轨道空间被均匀采样，同样采用了MCMC方法。

Chesley^[38]和Farnocchia等^[39]引入了轨道系统测距，Farnocchia等^[39]采用系统测距的方法在弱约束的地心距离和速率的空间中进行扫描，而位置和速率与观测值直接相关。根据观测数据可以获得小行星在天球上位置和速度的4个分量，即可归属项（α, δ, $\dot{\alpha }$, $\dot{\delta }$），再利用扫描的方式在均匀分布的网格中采样获取地心距离和速率两个分量（r, $\dot{r}$），即可得到小行星在极坐标的完整坐标。对每个网格固定r = r_i，$\dot{r}$ = $\dot{r}$_j，找到使函数最小的可归属项A_ij最佳拟合

其中：v为天体测量残差的向量；W为权重矩阵。Q的最小值通过微分校正得到

起始的A为(α₁, δ₁，(α_N–α₁)/(t_N–t₁)，(δ_N–δ₁)/(t_N–t₁))，N表示最后一个观测值。受约束的最佳拟合解可以转换为轨道，然后传播该轨道以寻找未来的密近交会。由于每个网格点都可找到可归属项的最佳拟合解，并且包含了协方差矩阵，因此可以进行碰撞概率的计算。Farnocchia等^[39]给出了2008 TC3、2014AA等多颗小行星的碰撞概率，根据观测点数量的不同，撞击概率均约为10^–3～1.0。

3.3   “锁眼”（Keyholes）

1998年关于小行星1997 XF11是否会在未来与地球相撞的讨论，推动了对撞击危险评估的重大研究，并导致了“锁眼”概念的确定。若地球轨道周期和密近交会后小行星的轨道周期构成简单公约，则未来再次在相同位置发生的密近交会称为共振回归。如果周期的比率很接近，则可能会发生后续的交会，但时间将比上一次早或晚^[40]。导致共振碰撞的b平面坐标位于可预测的Valsecchi圆上，轨道不确定性区域和Valsecchi圆之间的交点称为“锁眼”，由Chodas^[41]最早提出这一概念，它表示特定近距离相遇的b平面上的小区域，如果小行星经过其中的一个“锁眼”，那么它将在随后的回归中撞击地球。“锁眼”也可以用于指示b平面上的某个区域，该区域不一定导致之后的碰撞，但可以导致非常近的密近交会。“锁眼”与在给定日期发生的下一次碰撞的小行星的轨道半长径相关联。

Chodas^[41]等对 1997 XF11在 2028 年经过不确定性区域内的“锁眼”后引起的碰撞进行了计算，发现最高的碰撞概率为10^–5。Valsecchi等^[40]给出了求解半长径的方法，并讨论了对应于给定半长径的b平面上预测圆的结构。研究假设b平面上的水平距离ξ在2次相遇之间基本保持不变，而沿ζ的垂直距离在2次相遇传播的过程中在改变。选取ξ₁，ζ₁和ξ₂，ζ₂计算$\xi ''_{1}$，$\zeta ''_{1}$，$\xi ''_{2}$，$\zeta'' _{2}$；选择ξ₁和ξ₂，使得ξ₂=ξ₁，这里$\xi ''_{1}$与$\xi'' _{2}$相差较小。检查是否有$\zeta''_{1}\zeta ''_{2}$<0，若大于0则换一对$\zeta$直到式子成立。平行于$\zeta$轴的线段需要横跨一个共振回归，在该点上，找到使${\zeta'' }$=0的点，此点设为（ξ₀，ζ₀），如果$| \xi''_{0}|$<${b}_{\oplus}$（${b}_{\oplus}$为考虑引力作用下的地球半径），

其中$：{r}_{\oplus}$为地球半径。则($\xi''_0,\sqrt{b^2_{\oplus}-\xi''^{2}_0}$)的预图像为($\xi_0,\sqrt{\dfrac{b^2_{\oplus}-\xi''^{2}_0}{\partial\zeta''/\partial\zeta}}$)，预图像位于第2次相遇的b平面上，当小行星初次密近交会导致轨道发生改变后，在下一次交会前会经过该点。同时Valsecchii等^[40]还预测了1999 AN10在2037年和2032年的回归。

“锁眼”可用于危险小行星的轨道偏转策略优化，如果“锁眼”的尺寸较小，则在碰撞发生前使小行星发生小幅度的偏转并移出“锁眼”相比其它的偏转方案会更加容易。

3.4   Yarkovsky效应

对于大多数小行星来说，影响碰撞预测的不确定性是由计算轨道的不确定性引起的。但对于某些小行星，不确定性的主要来源是非引力摄动，尤其是Yarkovsky效应^[42]，主要表现在半长径的长期变化。当轨道受到很好的约束，并且当小行星碰撞概率需要考虑长期传播的影响时，Yarkovsky效应也会成为撞击危险评估问题的关键考虑因素。Milani等^[43]发现在模拟Bennu的轨迹和评估其撞击危险时，需要考虑的关键因素是Yarkovsky效应，这种由各向异性发射的热辐射引起的反冲加速度会导致半长轴的变化^[44]。Yarkovsky效应引起的半长径变化为

其中$：\rho$为小行星的密度；R为小行星的半径；$\gamma$为小行星赤道平面相对于其轨道平面的倾角。当自转方向发生变化时，半长径变化的符号会改变。对于1999RQ36，计算得到半长径的变化约为–12.5 ± 5 × 10^–4 AU/10⁶ a。Chesley等^[45]发现了2175—2196年的几种可能的密近交会和撞击，累积的撞击概率为3.7 × 10^–4。

Farnocchia等^[46]通过统计分析的方法评估了小行星（99942）Apophis撞击地球的风险，将轨道解和Yarkovsky效应的不确定性考虑在内。采用Monte Carlo方法模拟了Yarkovsky效应，充分考虑了物理特性的不确定性，尤其是自转方向的不确定性。之后将不确定性信息映射到2029年的b平面上并识别了与后续撞击相对应的“锁眼”，评估了未来碰撞的危险等级。之后Farnocchia等^[47]又对Bennu小行星在考虑Yarkovsky效应的情况下进行了碰撞危害的评估。Vokrouhlický^[48]结合Yarkovsky效应将不确定性映射到2029年的密近交会，并计算了导致共振回归的已知“锁眼”的权重，修改了21世纪下半叶Apophis的撞击概率。Spoto等^[49]在加入Yarkovsky效应后将轨道空间从6维扩展到7维，并采用LOV方法和Monte Carlo方法计算了2009 FD的碰撞概率。Tardioli等^[50]在Farnocchia的基础上利用“锁眼”和Yarkovsky效应研究了小行星碰撞概率的上限。

3.5   碰撞概率计算方法效率和采样完备性的改进

3.5.1   Monte Carlo方法的效率提升

Romano等^[51]采用2种基于Monte Carlo方法的采样方式，线采样^[52]和子集模拟^[53-55]，提高了碰撞概率经典的Monte Carlo方法的性能。两者以不同的方式对初始不确定性区域进行采样。线采样方法通过使用线而不是点来搜索不确定性区域，然后沿着线的方向进行积分从而计算碰撞概率，提高了碰撞概率计算的准确性。在确定采样方向后，在采样区域内随机产生${\theta }^{k}$（k = 1 ,2, ···, N）以生成线进行搜索，以函数Y(c)的值来评估沿着线搜索不确定性区域时线与碰撞区域的交点c

对于每条线${\theta }^{k}$，和撞击区域之间最多可以找到2个交点使得$Y({c}_{1}^{k})=Y({c}_{2}^{k})=0$。则随机数${\theta }^{k}$对应的碰撞概率为

其中：$\phi$表示单位高斯概率，总的碰撞概率为

子集模拟方法将碰撞概率计算过程变为条件概率的乘积，通过逐步识别向碰撞事件趋近的条件，减少采样所需的样本总数。给定碰撞事件I，令I₁$\supset$I₂$\supset$···$\supset$I_n=I为一系列的中间事件，则I_k=$\bigcap _{i=1}^{k}{I}_{i}$，碰撞概率可表示为$P$(${I}_{n}$)=$P$(${I}_{1}$)$\prod _{i=1}^{n-1}P\left({I}_{i+1}\left|{I}_{i}\right.\right)$。以计算最小地心距离的时间区间为区分不同事件的标准，一旦单个事件及其历元时刻被确认，可以得到计算最小地心距离的样本的时间区间，计算最小地心距离后可得到事件${I}_{1}$，采用MCMC方法在新区域生成样本，从而得到事件${I}_{2}$，重复此过程则得到最终的碰撞概率。Romano等^[51]对2010 RF12、2017 RH16和Apophis小行星的碰撞概率进行了计算，在计算结果精度相当的情况下提高了计算效率。

3.5.2   碰撞伪观测值法（Impact Pseudo-Observation，IOBS）

Roa等^[56]提出了在轨道确定的过程中，将碰撞条件作为观测值寻找参数空间中导致碰撞区域的方法。撞击伪观测值的残差是近距离接近时 b 平面的坐标，不确定性设为撞击位置的地球横截面弦的1/2，若与观测值兼容，则进行轨道外推最终收敛到发生碰撞的解。观测数据通常约束轨道不确定性分布的强方向，而撞击伪观测值约束弱方向。该方法外推轨道时沿着弱方向迭代，以最小化MOID的值飞向地球，在不简化动力学假设的情况下，探索了轨道不确定性的多维分布空间。

搜索虚拟碰撞体的步骤分为2步：首先确定会发生密近交会的轨道，之后依据密近交会的时间对轨道进行分类寻找会发生碰撞的轨道。IOBS方法在对密近交会的轨迹进行分类后，对在 b 平面上的交点加入残差作为伪观测值，约束不确定性的弱方向

其中：${b}_{{\rm{s}}}\left({t}^{*}\right)$表示历元${t}^{*}$时刻相对于地球的b平面上的坐标，$\|{b}_{{\rm{s}}}\left({t}^{*}\right)\| < {R}_{\oplus}$时即发生碰撞；${R}_{\oplus}$为考虑到大气的地球半径。沿着穿过地球中心的撞击轨道残差${e}_{{\rm{b}}}$为0。加入一个极小的固定的偏移量$\epsilon$，${e}_{{\rm{b}}}=\epsilon-{b}_{{\rm{s}}}\left({t}^{*}\right)$ 。撞击伪观测量的不确定性为

其中：$\alpha$为b平面上不确定性椭球的缩放因子，通常$\alpha$=1。在加入伪观测值后，外推过程经过不确定性区域向碰撞事件演变，而不是在不确定区域中寻找会发生碰撞的轨道。

最后根据轨道不确定性的初始分布利用Monte Carlo方法采样可计算碰撞概率。当虚拟碰撞体在参数空间中的分布被概率密度函数完整描述时，则可以采用重点抽样等方差缩减方法对虚拟碰撞体进行采样后计算碰撞概率。引入碰撞概率函数p_i(x)，则碰撞概率可以表示为

其中：F_i为参数空间中在某一日期撞击地球的轨道所属的区域；$q\left(x\right)$为观测值置信椭球的概率密度分布。计算的难点在于找到适当的${{p}_{i}\left(x\right)}$对虚拟碰撞体的分布进行建模，而引入伪观测值则可以解决这一问题，因为碰撞解的协方差对虚拟碰撞体分布进行了描述。

Roa等^[56]采用此方法计算了2008 TS10、2017 LD等小行星的碰撞概率，与SENTRY系统相比，在10^–7碰撞概率的情况下，二者计算结果相当，当碰撞概率减小到4 × 10^–8后，该方法更加稳健。

4   碰撞概率计算方法总结和展望

4.1   碰撞概率计算方法框架

综合以上对于碰撞概率计算方法的梳理，目前主要的碰撞概率计算流程的框架图如图3所示，在置信椭球体较为简单时可直接采用线性方法快速计算碰撞概率。当线性方法失效时，可根据获得的轨道根数是否为可靠的最小二乘解，选择在置信区域或容许区域内进行采样，依据计算成本和精度要求选择不同采样方式获得虚拟小行星，对虚拟小行星进行外推获得可与地球发生碰撞的虚拟碰撞体，最终得到小行星的碰撞概率。同时，在计算过程中可以加入IOBS等方法对计算的效率和稳定性进行改善。通常碰撞概率的计算软件采用其中的一种方法进行计算，也会同时列出线性方法和非线性方法进行自主选择。

图  3  碰撞概率计算流程图

Fig.  3  Flow chart for impact probability

下载: 全尺寸图片 幻灯片

第1个被开发出来的碰撞监测系统CLOMON，及后续改进的Sentry 和 CLOMON2 系统均采用了置信区域下的LOV方法。针对临近小行星进行碰撞监测的系统SCOUT，NEORANGER和NEOScan系统则采用了容许区域下的MOV方法。2021年7月Sentry-Ⅱ采用了新的IOBS方法进行碰撞概率的计算，在计算速度上Sentry-Ⅱ相比Sentry系统要慢，但搜索虚拟碰撞体的稳定性和准确性得到了大幅提升。

4.2   关键问题和未来展望

目前小行星碰撞概率计算方法在国外已较为成熟，后续改进主要集中于提高计算效率和提升搜索的完备性。

大量巡天观测计划的实施使得近地小行星观测数据海量剧增，其中短弧数据越来越多，而基于短弧定轨对碰撞概率的快速评估仍然是危险评估的难点，因此有效利用短弧数据对近地小行星进行危险评估也成为未来研究的重点。另一方面，对于信息量缺乏的短弧数据，也可通过轨道关联方法将多段短弧连接为足以进行精确计算的弧长，进而计算准确的碰撞概率。我国可以通过建立可靠的短弧段关联的计算方法极大地促进小行星碰撞概率计算准确性的提高。

实际问题中小行星的轨道较为复杂，采用单一评估方法计算负荷高且可能导致撞击概率的计算遗漏，后续可对不同观测情况的小行星分类分阶段处理，对拥有长期观测弧段的小行星进行高精度的碰撞概率计算，对即将发生碰撞的小行星则进行快速预测，从而建立可应对多种场景的中国独立自主的小行星监测预警系统。

小行星碰撞概率计算的准确性依赖于轨道精度，高精度的观测资料可以很好地约束置信区域，减少计算成本，简化计算方法。因此，采用多源数据观测可提高观测精度，从根本上提升小行星碰撞概率计算的准确性。中国目前正在建立的“中国复眼”，以及规划中的大口径高精度天体测量望远镜，对小行星监测预警系统的建立有着重要意义，准确的多源观测资料将会有力支持小行星防御系统的建设。

5   结束语

小行星防御需要对小行星进行危险评估和潜在威胁小行星的编目管理，而碰撞概率的计算是评估近地小行星威胁程度最关键的因素。本文梳理了现有计算小行星碰撞概率的基本方法和框架流程，并对中国未来的发展方向进行了展望。鉴于中国目前在碰撞预警方面的能力与国外有较大差距，中国应找准突破口，从自身特色出发建立适用于中国的监测预警系统。