0 引言
近地天体（Near - Earth Objects，NEOs）是指轨道近日距小于1.3 au的小行星和彗星，分别称为近地小行星和近地彗星。截至目前，已确认的近地天体已达3万多颗（截至2024年11月11日，近地小行星有36485 颗，近地彗星有143颗，另还有90余颗彗星碎片。数据来自MPC）。从轨道范围的角度，近地小行星可分为4个类型，分别是：阿莫尔型（Amors）、阿波罗型（Apollos）、阿登型（Atens）和阿蒂娜型小行星（Atiras），其中阿蒂娜型又称阿波希利（Apohele）小行星，轨道完全在地球轨道内部（远日点Q < 0.983 au），因此也称为地内小行星，其受观测几何限制很难被观测到，目前只发现50颗左右。另外还有一类最新提出的亚类Vatiras，是从上述阿蒂娜型中分出来的亚类，其轨道完全在金星轨道以内（远日点Q < 0.718 au），已发现数量更少。
以上4个类型中，阿波罗型和阿登型小行星的轨道与地球轨道在黄道面上的投影有交叉，有相当的风险与地球相撞，两者又常被合并称为穿越地球轨道的小行星（Earth - Crossing Asteroids，ECAs）。阿莫尔型小行星在地球轨道的外侧，若Yarkovsky效应使其向内迁移，或者与火星等发生密近交会后轨道发生变化，便有与地球相撞的风险。阿蒂娜型小行星的轨道完全在地球轨道之内，若受到Yarkovsky效应等影响向外迁移，也有与地球相撞的风险，同时因其大部分时间与太阳的视角距很小，难被地面上的观测设备发现，所以很难提前预警。
为了增加行星防御研究的针对性，近地小行星中其轨道与地球靠近（最小轨道相交距离MOID < 0.05 au）且体积较大（直径D > 140 m或绝对星等H < 22）的一个亚类被定义为潜在威胁小行星（Potentially Hazardous Asteroids，PHAs）。满足PHAs条件的小天体虽然其瞬时轨道与地球轨道的距离很小，但由于对其平经度没有限制，一般不会与地球同时到达两者轨道的交会处，短期内不会产生真正的威胁。Zhou等[1]提出密近交会潜在威胁小行星（PHAs that truly undergo Close Approaches，CAPHAs）的概念，即与地球交会距离小于0.05 au的小行星，相比PHAs进一步聚焦到短期内对地球有威胁的小天体。为了进一步研究地球引力范围内小行星的行为特征，本文又进一步引入高威胁小行星的概念。高威胁小行星（Highly Hazardous Asteroids，HHAs）指进入地球希尔半径（约0.01 au）的小行星 2 ，从轨道类型上看是PHAs中的一个亚类（这里的HHAs没有尺寸限制）。小行星进入地球希尔半径后，若轨道能量足够大（相对地球的总能量为正值），则飞会直接掠地球引力场范围，其相对地球的轨道可近似为双曲线；若轨道能量较低，则会被短暂俘获，成为地球临时俘获的自然卫星（Natural Earth Satellites，NES），又称minimoons。截至2024年11月已发现五颗NES，包括1991 VG、2006 RH120、2020 CD3、2022 NX1和2024 PT5。
对于我国未来的行星防御监测的天基布局，为了更有效地监测从太阳方向靠近地球的小行星，Wang等[7]和Zhou等[8]分别提出了将设备放置在地球领航轨道（Earth-Leading Heliocentric Orbit，ELHO）上的监测方案和CROWN计划（Constellation of heteROgeneous Wide - field Nea surveyors）。为了分析监测效率，他们利用已知的2000 多颗PHAs作为观测样本，由此计算项目期间可以观测到的样本比例。实际上，由于观测选择效应的影响，已发现的PHAs并不能代表实际的PHAs群。除了尺寸上的选择效应，轨道分布的不同会影响观测效率的计算，从而导致不同轨道特征的PHAs被已有设备观测到的概率也不同。所以如果得到近地小行星或PHAs或高威胁小行星的“真实”轨道尺寸分布，就可以进一步更合理地评估待部署设备的监测效率。
很明显，行星防御的真正对象是地球撞击小行星，而这一群体与上面提到的诸多群体的轨道分布又各自有所差异。目前已知的地球撞击小行星的轨道样本极少，亟须以生成虚拟目标的方式来生成这类小行星的仿真样本数据库。文献[9]-[10]从近地小行星的轨道分布出发，首先用初始瞬时轨道筛选出MOID小于地球俘获横截半径（the Earth's Capture Cross Section，ECCS） 1 的小行星，然后计算小行星轨道处于ECCS内的部分占整个轨道的比例，选择这一比例大的作为下一步的候选体，最后再通过太阳系动力学模型推演来确定最终的撞击小行星样本。此筛选方法首先忽略了MOID短期内发生突变的小行星，其次会优先选择低偏心率低轨道倾角且半长径接近1 au的类地球轨道。由于文献[9]-[10]的工作是在“地球撞击体中类地球轨道占大多数”的假设前提下开展的，因此在根数分布上会存在选择性偏差。本文将通过计算不同轨道构型下与地球发生交会的占比，来分析这一假设的合理性。
1 近地小行星轨道/尺寸分布生成方法
本节基于近地小行星的形成演化理论，并结合当前的实测数据，给出近地小行星的轨道/尺寸分布的理论依据及实现过程。
主带小行星如果与大行星发生长期共振（两者的近点角距的进动率相近），或者发生轨道共振（周期比接近简单整数比，且共振角发生秤动），小行星的偏心率会逐渐被激发，从而导致其轨道近日距变小，小行星进入内太阳系。高偏心率的小行星会与内太阳系大行星发生频繁的密近交会，导致其轨道不稳定，从而散射出太阳系或与大天体发生碰撞而消亡；还有部分高偏心率小行星在与太阳的潮汐耗散作用下偏心率和半长径同时减少，其轨道变成内太阳系的近圆轨道；如果小行星在初始进入内太阳系时的偏心率足够大，其近日点离太阳足够近，还可能发生热裂解，从而分裂成更小的碎片，发生裂解的临界日心距与小行星的尺寸成反相关[16]。 基于上述的作用机制，处于共振中的小行星在约Myr时标内便被清空，若要维持约Gyr时标内对近地小行星群的持续补给，还需要Yarkovsky效应下小行星轨道的径向迁移机制。在Yarkovsky 效应下，处于共振附近的主带小行星会持续进入共振区域，从而为近地小行星提供一个长效的补给机制。综上，Yarkovsky效应下的径向迁移联合大行星长期共振、轨道共振等散射作用是目前普遍认为的近地小行星的补给机制，而近地小行星的轨道/尺寸分布也在这种变化中达到准静态平衡。当前所用的近地小行星样本即是在这一理论下模拟得到。
基于上述理论，文献[11]-[13]和[17]-[18]结合N体演化模拟结果和去偏观测结果，利用拟合分布参数的方法得出近地小行星的尺寸/轨道分布。首先，依据已发现的大尺寸主带小行星的轨道分布，以不同源区为单位，选取真实的或者生成虚拟的主带小行星及木星族彗星轨道，作为生成虚拟近地小行星样本的初值；其次，在太阳系动力学环境中演化若干百万年，记录已变为近地小行星的轨道，按照（a，e，i，H）四维空间中的区间来统计（变量含义见附录A）；再次，将上述模拟结果与去偏后的观测结果以区间为单位相对比，利用最大似然估计法计算预设分布函数中的未知参数。上述3个工作均用了类似的方法，不同点在于：1）源区的设定；2）先验的轨道/尺寸分布函数；3）实测样本数量。本文利用文献[11]-[13]中得出的分布及对直径10 m以上近地小行星数量的估计，随机生成700万颗虚拟近地小行星，作为研究工作的初始样本。关于初始样本数量的讨论见第4节。
2 计算模型和初值条件
本文采用自主开发的小天体轨道外推程序包（Small Body Orbit Propagation Package，SBOPP）进行小行星轨道计算，力模型的摄动项考虑了八大行星及月球的摄动影响、太阳的J2项、太阳广义相对论效应以及地球的J2项。通过对700万颗产生于NEOMOD3分布模型[13]的虚拟小行星进行轨道外推，从2025/01/01积分至2125/01/01（虚拟小行星的轨道根数初值历元设为MJD=60676.0 ），记录下进入地球Hill半径的小行星的运动状态和最近交会时刻。由于SBOPP程序包采用DE440历表提供的各大行星的运动状态作为小行星摄动力的计算输入，所以其中发生的小行星与地球的密近交会即是与真实地球轨道的交会，其计算结果可作为评估近地小行星监测预警能力的虚拟高威胁小行星样本。本文涉及名词及各变量含义见表1 。
表 1 本文涉及名词及各变量含义Table 1 Meaning of terms and variables in thispaper |
| 缩写/符号 | 全称/中文名称 | 含义 |
| NEOs | Near-Earth Objects, 近地小天体 | q < 1.3 au的小天体 |
| NEAs | Near-Earth Asteroids, 近地小行星 | q < 1.3 au的小行星，将q < 1.3 au的彗星排除在外 |
| ECAs | Earth-crossing Asteroids, 穿越地球轨道的小行星 | q < 1.017au 且 Q > 0.983 au；或者阿波罗型和阿登型 小行星的总和 |
| Apollos | 阿波罗型小行星 | a > 1.0 au 且q < 1.017 au |
| Amors | 阿莫尔型小行星 | 1.017 < q < 1.3 au |
| Atens | 阿登型小行星 | a < 1.0 AU，且Q > 0.983 au |
| Atiras/Apohele/IEOs | 阿蒂娜型小行星/阿波希利小行星/地内小行星 (interior-to-Earth objects) | Q < 0.983 au |
| Vatiras | — | Q < 0.718 au |
| PHAs | Potentially Hazardous Asteroids, 潜在威胁小行星 | H<22 mag, 且MOID<0.05 au; 或者D>140 m, 且MOID<0.05 au |
| PHA10 | 直径大于10 m且MOID<0.05 au的近地小行星 | MOID<0.05 au，D>10 m |
| CAPHAs | PHAs that truly undergo Close Approaches, 密近交会潜在威胁小行星 | |
| HHAs | Highly Hazardous Asteroids, 高威胁小行星 | |
| — | encounter popultations交会子类 | CAPHAs与HHAs的合称 |
| NES/minimoons | temporarily-captured Natural Earth satellites, 地球临时 俘获的小行星/迷你月亮/地球小卫星 | |
| — | Earth-impacting asteroids/Earth Impactor，撞击小行星 | |
| — | Earth-like orbits 类地球轨道 | a |
| — | 地球Hill半径/地球希尔半径 | 地球引力相对太阳引力起主导作用的地球周围球形 范围的半径，约为a(1-e)(m/3M)1/3 |
| a | 轨道半长径 | 六个轨道根数之一 |
| e | 轨道偏心率 | 六个轨道根数之一 |
| i | 轨道倾角 | 六个轨道根数之一 |
| q | 近日距 | a(1-e) |
| Q | 远日距 | a(1+e) |
| H | 绝对星等 | 假设小行星日心距是1 au，地心距是1 au，相位角为 0度情况下的亮度 |
| D | 小行星直径 | — |
| 小行星到地球质心的距离 | — | |
| R♁ | 地球半径 | — |
| MOID | minimum orbit intersection distances, 最小轨道相交距离 | 在两条轨道上各任取一点，两点之间的最小值 |
| MPC | minor planet center/国际小行星中心 | 网址为https://www.minorplanetcenter.net |
3 结果
3.1 各子类占比
本小节统计了各子类在近地小行星中的占比。其中PHAs的判定是基于小行星与地球在某时刻的瞬时轨道计算的，这也与MPC等现行数据库中的判定方法一致。而其他的子类则是考虑了100年内的轨道演化，凡是100年内出现满足某子类条件的小行星均归为其中。显而易见，这些子类的成员是随其轨道演化不断变化的，累积时段越长，子类的成员数量越多。此处选择100年作为一个累积时段，是综合当前历表精度及行星防御实施效率的一个折中考虑。CAPHAs与HHAs中同一天体多次回归按一次计算。
表 2 虚拟近地小行星中各子类的总数目以及两种占比Table 2 Total number of subclasses in virtual near-Earth asteroids and their proportions |
| 名称 | PHA10 | CAPHAs | HHAs | 迷你月亮 | 撞击小行星 |
| 数目 | 1 | 6 | |||
| 占比1/% | 24.7 | 12.2 | 0.8 | ||
| 占比2/% | 100 | 48 | 3.3 |
表 3 PHA10包含不同累积时间段内CAPHAs及HHAs的百分比Table 3 PHA10 contains the percentages of CAPHAs and HHAs in different accumulation periods |
| 累积时间段 | 10年 | 50年 | 100年 |
| CAPHAs | 99.8% | 99.0% | 97.7% |
| HHAs | 99.97% | 99.79% | 99.35% |
从表2 看出，迷你月亮发生的概率很低，甚至比撞击小行星的概率还低.这表明，在进入地球附近区域的小行星中，绝大多数（大于99.99%）是以双曲轨道擦肩而过的，仅有极少数会被地球短暂俘获。虽然迷你月亮和撞击小行星在试验中发生的次数很少，但两个子类间并没有交集。换句话说，成为迷你月亮的小行星演化趋势大概率是逃逸地球引力场，不会与地球发生碰撞，文献[19]中给出迷你月亮中发生撞击地球的小行星约占1%；而撞击小行星大多以双曲轨道撞向地球，撞击前不会经历迷你月亮的阶段，撞击小行星中迷你月亮的占比还需要更大的仿真样本来计算。初始PHA10小行星在之后的100年内有近一半会发生与地球小于0.05 au距离的交会。瞬时历元的PHA10占NEAs的24.7%，但却包含了100年累积时间段内的97.7%的CAPHAs和99.35%的高威胁小行星（见表3 ）。也就是说，NEAs中不属于PHA10的75.3%只包含了CAPHAs中的2.3%和HHAs中的0.65%，这说明MOID>0.05 au的轨道在短期内与地球发生近距离交会的概率很低，也说明将PHA10作为行星防御中的监测对象较为可靠。
也有极少数不包含在初始PHA10中的小行星，之后的100年演化中变为了CAPHAs和HHAs，即2.3%的CAPHAs和0.65%的HHAs，对应的数目分别为19381 颗和378颗，下一小节中我们进一步讨论这些小行星的轨道特征；
本文还计算了在D>140 m直径限制下的PHAs的“实际”数量，共4825 颗，与文献[20]中给出的4700 ±150的预测数目相近。相比目前已知的2446 颗（来自MPC，截至2024年11月11），约发现了其中的51%。
3.2 各子类轨道根数分布
本小节对比了近地小行星中各子类的轨道根数分布。虽然各子类都属于近地小行星群，保留了整体样本的轨道根数分布主要特征，但由于各自侧重不同，其根数分布还是有部分区别，从而也体现了其不同的几何构型或动力学特征。
（1）CAPHAs和HHAs的轨道半长径和偏心率分布相近，但轨道倾角分布有较大差异（见图1 （c））。4个子类相比，低倾角样本占比与地球威胁程度呈正相关，这也与文献[9]中的结论一致。
（2）与其余两个子类相比，CAPHAs和HHAs的轨道半长径分布明显向左偏移，且在1~2 au范围内的分布维持稳定（见图1 （a））；由于NEAs偏心率分布在0.5以上开始下降，从而导致大于2au的CAPHAs和HHAs数量减少；PHA10在1au附近的占比比NEAs略大，而CAPHAs和HHAs在1au附近的发生频率明显大于另外两个群的分布，这可能是由于1 au附近穿越地球轨道的小行星占比更多导致的。
（3）4个子类的偏心率分布差异不大，但值得注意的是，无论是PHA10、CAPHAs还是HHAs，都没有在低偏心率附近出现频率增大的倾向（图1 （b）），这一情况与文献[9]中给出的分布类似，但也说明，与地球发生交会或撞击的小行星轨道并不是类地球轨道占多数，高偏心率轨道占比更大。
交会子类的轨道倾角分布聚集在0°附近，但半长径和偏心率的分布范围很宽，下一小节将从交会占比角度对此问题做进一步的阐述。
另外，本文针对未包含在初始PHA10样本中的19381 颗CAPHAs和378颗HHAs的轨道根数分布特点作进一步分析。图2 给出了这两类小行星分别与CAPHAs、HHAs整体的轨道根数数密度分布对比，其中蓝色粗线表示CAPHAs中不包含在初始PHA10中的小行星；橙色粗线表示HHAs中不包含在初始PHA10中的小行星。可以看出，不包含在初始PHA10中的小行星的根数分布与CAPHAs、HHAs的整体分布差异明显：（1）轨道半长径在2.1 au、2.5 au、2.8 au附近聚集，对应木星的4∶1、3∶1、5∶2内平运动共振位置；（2）轨道偏心率和轨道倾角的分布均显示两子类的激发型轨道占比更大。这表明，未包含在初始PHA10中且100年内与地球密近交会的小天体大部分为木星轨道共振区附近的激发型轨道，无法通过计算其瞬时MOID来提前预测其危险性。这种轨道与木星发生周期性会合，受木星摄动影响更强，其轨道根数的大幅振荡导致其与地球的MOID会有短期内的快速变化。这种轨道类型的小行星需要格外引起注意，即使它目前的MOID比较大。
3.3 各轨道类型交会比例
（1）随偏心率的增大，这一比例变小的趋势略显平缓，到e~0.4时仅下降了不到50%，而且在e~0.6的时候又出现了第二个峰值。这说明低偏心率轨道在与地球发生交会的概率上并无明显优势。
（2）占比随轨道倾角的变化趋势表明，近零倾角轨道与地球交会的占比相比高倾角轨道是占绝对优势的，并且HHAs相比CAPHAs，近零倾角轨道的占比优势更大。这表明将与地球交会的小行星轨道或者与地球碰撞的小行星轨道假设成近零倾角的轨道是相对合理的。
（3）与地球交会比例在a~2 au和e~0.6处出现峰值，并不表示交会比例在a或e的对应范围内更大，而是由于初始近地小行星分布中在这两处低倾角轨道占比更多造成的。图4 统计了本文所用的近地小行星初始轨道中，不同轨道倾角范围内近地小行星占比随半长径和偏心率的变化。可以看出，在a~2 au附近和e∈（0.4，0.6）范围内，低倾角轨道（i∈（0，10˚））的小行星占比出现局部的峰值，而由于低倾角轨道在交会概率上有绝对的优势，从而导致交会概率随a，e增加时的第二个峰值。
3.4 各子类数据库的大样本生成
本小节将介绍各子类大样本数据库的生成方法。上文中已经提到，近地小行星中的CAPHAs和HHAs子类数量都是跟累积时间段的长短直接相关的，累积时间段越长，其间与地球发生交会的CAPHAs和HHAs也就越多。但对于负责小行星监测预警的空间任务来说，其寿命是有限的，3年、5年或10年不等。表4 给出了这些相应时间段内由初始700万颗近地小行星得到的CAPHAs和HHAs的数量，可以看出，对高威胁小行星，有限的时间段内的数量并不是太多。这主要是因为我们首先设置了直径大于10 m的前提条件，如果包含直径更小的小行星，这一数量肯定会成指数地增加。
表 4 相应时间段内由初始700万颗近地小行星得到的CAPHAs和HHAs的数量Table 4 The number of CAPHAs and HHAs obtained from the initial 7 million near-Earth asteroids in the corresponding time period |
| 累积时间段 | 3年 | 5年 | 10年 |
| CAPHAs | |||
| HHAs |
为了得到更大样本的高威胁小行星轨道，从而更有效地评估各种小行星监测方案的效率，我们利用现有样本中a、e、i的三维分布，从任取的大量与地球交会的轨道中选出满足分布的部分，从而得到在给定时间段内用户所需数量的高威胁小行星轨道样本（算法流程图如图5 所示）。再从这一大样本的高威胁小行星轨道出发，对每条轨道作短时间推演，从中筛选出撞击小行星，继而可得到给定时间段内的撞击小行星轨道样本。
4 结果讨论
Fang等[21]基于Bottke等[17]提供的近地小行星的轨道/尺寸分布，给出了与地球交会的小行星在无穷远处相对地球的速度分布，Heinze等[22]基于Granvik等[18]的近地小行星分布给出了交会距离小于0.05 au的小行星在交会时刻相对地球的速度分布。为了与其结果做对比，我们分别计算了CAPHAs小行星和高威胁小行星在无穷远处相对地球的速度的密度分布，以及CAPHAs小行星交会时刻相对地球的速度的密度分布（见图6 ），图中绿色实线是CAPHAS小行星在无穷远处相对地球的速度分布，红色实线是高威胁小行星在无穷远处相对地球的速度分布，绿色虚线是CAPHAs交会时刻相对地球的速度分布这里相对横坐标做了归一化。
通过比较可以看出，本文给出的速度分布在中等速度10~15 km/s区间有更大的占比，同时在较高速度20 km/s以上区间的占比相比其他两个工作的结果更低，而小于10 km/s的低速区间的占比则处于其他两个结果中间。高速度区间内占比更低的原因或许可以从所用近地小行星的初始分布中找到线索。图7 给出了文献[18]和[13]（本文所用）给出的近地小行星样本的轨道根数分布对比，纵坐标含义同图1 ，没有做归一化。可以看出，文献[13]所给的样本中，高偏心率和高倾角的轨道占比相对更少，而这部分轨道在与地球发生交会时会有更大的相对速度，从而导致本文给出的速度分布中高速度部分占比更低。理论上得出的相对速度分布是否符合实际情况，还需要与观测进行进一步的比对。
5 结论
本文利用文献[11-13]中提供的近地小行星样本，计算了700万颗直径大于10 m的虚拟近地小行星在2025−2125年的轨道演化，筛选出积CAPHAs、HHAs等多个子类，计算各子类在近地小行星中的发生概率： CAPHAs为12.2%、HHAs 0.8%、迷你月亮 0.00001 %、撞击小行星 0.00009 %，以及各个子类中属于PHA10的部分在PHA10中的占比：CAPHAs 48%、HHAs 3.3%。从中可看出：（1） 迷你月亮与撞击小行星发生的概率极低，且相互并没有包含关系；（2） PHA10中有48%在之后的100年内会变为CAPHAs，即发生与地球0.05au距离内的密近交会。进一步计算了PHA10中包含不同时段内CAPHA及HHAs的比例，结果显示PHA10包含100年内97.7%的CAPHAs和99.35%的HHAs，证实了初始瞬时PHA10群作为监测预警对象的可靠性。
通过CAPHAs及HHAs的轨道根数分布及随轨道根数变化的交会占比，说明了与地球交会轨道并不是类地球轨道占多数，因为轨道半长径和轨道偏心率分布更多取决于初始的近地小行星群的轨道分布。利用本文给出的CAPHAs和HHAs的根数分布，可以生成任意时间段内任意数量的CAPHAs和HHAs样本，从而为更精确地仿真行星防御中的监测预警布局和效率提供技术支持。
致谢
感谢两位审稿人的认真审阅及宝贵意见建议，增加了本文的逻辑性和可读性，并使内容更加完整。在此表示衷心的感谢！