基于引力锁孔面积计算的小行星撞击概率评估方法

1.本发明涉及一种小行星撞击概率评估方法，特别涉及一种基于计算引力锁孔面积的小行星撞击概率评估方法，属于航空航天技术领域。


背景技术：

2.潜在危险小行星(phas)是指那些最小轨道交叉距离(moid)小于等于0.05au且绝对星等小于等于22.0的小行星，根据iau的小行星数据库，目前有约2200颗小行星被列为phas。这些小行星受到地球和火星引力的摄动等，多次近距离飞越地球，具有撞击地球带来全球性灾难的可能性。因此，潜在危险小行星及行星防御受到广泛关注。量化评估phas撞击地球概率是行星防御中的第一步。
3.为了对phas可能的地球撞击做出预报，目前已发展了多种评估小行星撞击地球可能性的方法。在先技术【1】(参见chodas p.w.estimating the impact probability of a minor planet with the earth.[j]bull.am.astron.soc.,1993,25:1236.)评估小行星撞击地球概率研，航天器撞击地球概率计算方法被应用到小行星的地球撞击概率计算中，在小行星相对地球最近时刻的轨道状态进行采样，通过轨道积分得到下一次通过地球b平面的概率密度，对概率密度积分得到撞击概率。未考虑轨道误差长时间传播后分布规律，该方法概率密度函数的估计与实际情况相差较大，计算耗时长、精度低。
[0004]
在先技术【2】(参见chodas p.w.orbit uncertainties,keyholes,and collision probabilities.[j]bull.am.astron.soc.,1999,31:2804.)给出了地球共振小行星引力锁眼(keyhole)的定义，并预报了小行星1997xf11位于2028年的引力锁眼和导致小行星1999an10在2039、2044、2046年撞击地球的引力锁眼。该方法给出的引力锁孔对应弧段仅预报了小行星的近距离飞越位置与时间，定性分析了撞击风险，但未给出量化风险指标。
[0005]
在先技术【3】(参见chodas p.w.,yeomans d.k.could asteroid 1997xf11 collide with earth after 2028？[j]bull.am.astron.soc.,1999,31:703.)小行星1997xf11撞击地球风险研究，应用蒙特卡洛仿真在观测时刻对线性六维置信区域进行采样，然后使用非线性方程数值积分至共振轨道可能的撞击时刻，获得小行星在垂直速度平面上的狭长的位置误差椭圆，利用该椭圆的宽度计算最小飞越距离。该方法只能预报与地球形成共振轨道的小行星的撞击时刻并计算相应的撞击概率。
[0006]
在先技术【4】(参见milani a.,chesley s.r.,and valsecchi g.b.close approaches of asteroid 1999an10:resonant and non-resonant returns.[j]astron.astrophys.,1999,346:l65
–
l68.)小行星1999an
10
共振轨道与非共振轨道近距离飞越地球分析，同时考虑共振轨道与非共振轨道，分析近距离飞越时刻的非线性置信区域与lov的交点，预报了小行星1999an10的11个共振轨道返回并补充了14个非共振轨道返回，并进一步考虑moid的变化并预报小行星近距离返回时的最小飞越距离。该方法同样只能预报与地球形成共振轨道或与共振轨道十分相近的非共振轨道的小行星的撞击时刻并计算相应的撞击概率，不能对任意小行星轨道进行撞击概率计算。


技术实现要素：

[0007]
本发明公开的基于引力锁孔面积计算的小行星撞击概率评估方法要解决的技术问题是：提供一种适用于任意撞击时刻的撞击概率评估方法，能够建立小行星撞击条件与引力锁孔面积间的联系，通过计算引力锁孔面积与引力影响球截面的比值，得到对应时刻的撞击概率，实现小行星撞击概率评估。本发明具有适用于任意轨道周期的小行星任意时刻撞击地球概率的评估、扩展撞击概率评估的适用范围，节约计算资源、提高评估效率，提高概率评估可靠性等特点。
[0008]
本发明的目的是通过下述技术方案实现的：
[0009]
本发明公开的基于引力锁孔面积计算的小行星撞击概率评估方法，针对潜在危险小行星轨道范围特点，考虑中心天体引力、第三体引力摄动、太阳光压摄动、雅克夫斯基效应摄动，构建小行星高精度轨道动力学模型。在ti时刻测得小行星的轨道位置矢量和速度矢量，并添加观测误差，得到用于轨道预测的小行星轨道初值，基于构建的小行星高精度动力学模型进行蒙特卡洛打靶实现小行星轨道预测，并记录下小行星飞越地球位置速度矢量，根据位置速度矢量得到近心点半径和双曲超速参数；将所述近心点半径和双曲超速参数代入状态转移矩阵，利用状态转移矩阵推导保守撞击条件。判定小行星是否满足保守撞击条件，如果不满足保守撞击条件，将小行星撞击地球概率判定为0，排除蒙特卡洛打靶中不能够撞击地球的小行星轨道预测，节省小行星撞击概率评估运算资源，提高小行星撞击概率评估效率。如果满足保守撞击条件，利用变分法推导近心点半径偏差量和双曲超速方向角度偏差量的表达式，代入引力锁孔面积定义式得到引力锁孔面积计算公式；将已得到的小行星相对地球的近心点半径和双曲超速大小代入所述引力锁孔面积计算公式得到引力锁孔面积。根据蒙特卡洛打靶记录及保守撞击条件判定统计出满足保守撞击条件的预测的小行星轨道数目，根据计算的引力锁孔面积数值计算引力锁孔面积的累积分布函数，根据任务需求预设置信度选取引力锁孔面积估计值，根据满足保守撞击条件的小行星轨道预测数目占所有蒙特卡洛打靶的小行星轨道预测数目比例和引力锁孔面积估计值，计算出小行星在给定t
t
时刻撞击地球的概率，即基于引力锁孔面积计算实现小行星撞击概率评估。本发明评估的小行星撞击地球概率能够作为行星防御工程实施的判据，当小行星撞击地球概率超过阈值时，实施行星防御，避免小行星撞击地球。
[0010]
本发明公开的基于引力锁孔面积计算的小行星撞击概率评估方法，包括以下步骤：
[0011]
步骤一、针对潜在危险小行星轨道范围特点，考虑中心天体引力、第三体引力摄动、太阳光压摄动、雅克夫斯基效应摄动，构建小行星高精度轨道动力学模型。
[0012]
潜在危险小行星日心轨道主要位于地球至木星轨道范围，受到地球、火星和木星带来的第三体引力摄动。
[0013]
第三体引力摄动p
third
表达式为
[0014][0015]
其中，μ
pl
表示第三体的引力常数，r
pl
表示第三体相对中心天体位置矢量，δ表示
第三体相对小行星的位置矢量。
[0016]
太阳光压摄动p
srp
表达式为
[0017][0018]
其中，sa为小行星垂直太阳光横截面积，ε为太阳辐射通量，m为小行星质量，c为光速，n为太阳光方向单位矢量。
[0019]
雅克夫斯基效应是各向异性热辐射的辐射反冲导致小行星经历长期半长轴漂移，雅克夫斯基效应对小行星轨道也有较显著影响，表示为横向加速度
[0020][0021]
其中，a为反照率，φ(r)为标准辐射力因子，γ为倾角，θ为热参数且
[0022][0023]
考虑中心天体引力、如式(1)所述第三体引力摄动、如式(2)所述太阳光压摄动、如式(3)所述雅克夫斯基效应摄动，构建如式(5)所示的小行星高精度轨道动力学模型。
[0024]
小行星轨道动力学模型为
[0025][0026]
其中，μ为太阳引力常数，r为小行星相对太阳的位置矢量，为轨道面内与位置矢量垂直指向速度方向的横向单位矢量。
[0027]
步骤二、测得ti时刻小行星的轨道位置矢量和速度矢量，并添加观测误差，得到用于轨道预测的小行星轨道初值，基于步骤一构建的小行星高精度动力学模型进行蒙特卡洛打靶实现小行星轨道预测，并记录下小行星飞越地球位置速度矢量，根据位置速度矢量得到近心点半径和双曲超速参数；将所述近心点半径和双曲超速参数代入状态转移矩阵，利用状态转移矩阵推导小行星速度改变量与小行星-地球交会距离的关系，并构建保守撞击条件。判定小行星是否满足保守撞击条件，如果不满足保守撞击条件，将小行星撞击地球概率判定为0，排除蒙特卡洛打靶中不能够撞击地球的小行星，节省小行星撞击概率评估运算资源，提高小行星撞击概率评估效率；如果满足保守撞击条件，转移至步骤三。
[0028]
在ti时刻测得小行星的轨道位置矢量和速度矢量分别为rd和vd，将观测误差三维随机矢量ξ、三维随机矢量ζ、三维随机矢量η代入rd和vd，得到用于轨道预测的小行星轨道初值ri和vi分别为
[0029][0030]
其中，三维随机矢量ξ各分量分别为服从均值为0标准差为位置精度σr的随机数，三维随机矢量ζ各分量分别为服从均值为0标准差为速度精度σv的随机数，三维随机矢量η为导航误差导致的速度误差，其各分量分别为服从均值为0标准差为σg的随机数。
[0031]
以式(6)为初值带入模型(5)进行积分，当小行星在tk时刻前后近距离飞越地球，
即小行星轨道状态与地球轨道状态re,ve满足式(7)时，记录下该时刻和对应的小行星的轨道状态
[0032][0033]
此时，小行星相对地球的近心点半径和双曲超速大小根据式(8)(9)计算：
[0034][0035][0036]
其中，a、e分别为小行星在tk时刻前后近距离飞越地球时的轨道半长轴、偏心率，为μe为地球引力常数，||
·
||表示二范数。
[0037]
计算保守撞击条件涉及两条轨道，一条为标称撞击转移轨道，即在精确动力学模型下求解的从tk时刻地心位置到t
t
时刻地心位置的转移轨道；另一条轨道为实际小行星飞行轨道。其中，δr1,δv1分别为tk时刻实际小行星轨道与标称撞击轨道的位置、速度偏差，δr2,δv2分别为t
t
时刻实际小行星轨道与标称撞击轨道的位置、速度偏差。若实际小行星轨道与标称撞击转移轨道偏差δr1,δv1太大，则仅靠自然摄动力无法提供足够的轨道改变量使小行星轨道足够靠近t
t
时刻的地球导致撞击，只有当该小行星轨道与撞击轨道的偏差满足保守撞击条件时，该小行星才有撞击地球的可能性。
[0038]
已定义小行星轨道与撞击转移轨道的首、末偏差分别为δr1,δv1、δr2,δv2，由状态转移矩阵关系
[0039][0040]
得到
[0041][0042]
在日心系下，|δv1|即为小行星轨道与标称撞击轨道的双曲超速矢量之差|δv
∞
|，δr1等效为0，|δr2|为地球半径。因此式(11)具体为
[0043][0044]
由于式(12)为近似线性关系，出于保守估计考虑，当误差处于数倍e之内，仍将其考虑为可能的撞击情况，因此得保守撞击条件为
[0045][0046]
对于进入地球引力影响范围的小行星轨道，将式(8)、式(9)所计算数值，代入式(13)所示的保守撞击条件，只有当式(13)成立时才考虑该小行星在t
t
时刻有撞击地球的可
能性，将不满足条件的轨道判定在t
t
时刻撞击地球的可能性为0，排除蒙特卡洛打靶中不能够撞击地球的小行星，节省小行星撞击概率评估运算资源；如果满足保守撞击条件，转移至步骤三。
[0047]
步骤三、利用变分法推导近心点半径偏差量和双曲超速方向角度偏差量的表达式，代入引力锁孔面积定义式得到引力锁孔面积计算公式；将步骤二得到的小行星相对地球的近心点半径和双曲超速大小代入所述引力锁孔面积计算公式得到引力锁孔面积。
[0048]
小行星的引力锁孔面积定义式如式(14)所示。
[0049][0050]
式(14)中r即为小行星相对地球的双曲线轨道的近心点半径r
p
。利用变分法推导得到如式(15)所述近心点半径偏差量和式(16)所述双曲超速方向角度偏差量，
[0051][0052][0053]
其中|δv1|、δf
∞
和δv
∞
通过已经完成的实际小行星轨道数值结果和标称撞击轨道之差得到。
[0054]
将式(15)、(16)代入引力锁孔面积定义式(14)得到引力锁孔面积计算公式(17)；将步骤二得到的小行星相对地球的近心点半径和双曲超速大小带入所述引力锁孔面积计算公式得到引力锁孔面积。
[0055][0056]
步骤四、根据步骤二中蒙特卡洛打靶记录及保守撞击条件判定统计出满足保守撞击条件的预测的小行星轨道数目；然后根据步骤三中计算的引力锁孔面积数值计算引力锁孔面积的累积分布函数，并根据任务需求预设置信度选取引力锁孔面积估计值；最后根据满足保守撞击条件的小行星轨道预测数目占步骤二中所有蒙特卡洛打靶的小行星轨道预测数目比例和引力锁孔面积估计值，计算小行星在步骤二中给定t
t
时刻撞击地球的概率，即基于引力锁孔面积计算实现小行星撞击概率评估。
[0057]
步骤二蒙特卡洛打靶的k条小行星轨道预测中，共记录k条小行星轨道预测在tk时刻进入地球引力影响范围a
sphere
，所述k条轨道在步骤三中有n条轨道判定为满足式(13)所述的保守撞击条件。
[0058]
所述n条小行星轨道预测在步骤三中已分别计算出对应的引力锁孔面积。所计算引力锁孔面积略有不同，利用概率论基础理论计算所计算引力锁孔面积的累积分布函数(cdf)，并根据任务需求预设置信度，取置信度对应的引力锁孔面积数值作为该小行星在t
t
时刻撞击地球所对应的tk时刻引力锁孔面积估计值
[0059]
则在预设置信度下，小行星在t
t
时刻撞击地球概率计算式为
[0060][0061]
还包括步骤五：根据步骤四得到的小行星撞击地球概率能够作为行星防御工程实施的判据，当小行星撞击地球概率超过阈值时，实施行星防御，避免小行星撞击地球。
[0062]
有益效果：
[0063]
1、本发明公开的基于引力锁孔面积计算的小行星撞击概率评估方法，采用小行星高精度轨道动力学模型进行小行星轨道预测，而非共振轨道模型，因而不限于特定共振比的小行星轨道及共振返回时刻，能够实现任意轨道周期小行星在任意时刻撞击地球的概率评估，扩展小行星撞击概率评估的应用场景。
[0064]
2、本发明公开的基于引力锁孔面积计算的小行星撞击概率评估方法，由于采用状态转移矩阵推导出保守撞击条件，能够保守地判定出撞击地球概率为零的小行星轨道预测并不对所述撞击地球概率为零的小行星轨道预测进行概率计算，进而提高小行星撞击地球概率评估的效率。
[0065]
3、本发明公开的基于引力锁孔面积计算的小行星撞击概率评估方法，由于计算引力锁孔面积的累积分布函数并根据预设置信度取得引力锁孔面积估计值，而非简单的求取引力锁孔面积均值作为引力锁孔面积估计值，提高小行星撞击概率评估结果的可靠性。
[0066]
4、本发明公开的基于引力锁孔面积计算的小行星撞击概率评估方法，针对潜在危险小行星轨道范围特点，根据位置速度矢量得到近心点半径和双曲超速参数；将所述近心点半径和双曲超速参数代入状态转移矩阵，利用状态转移矩阵推导保守撞击条件。判定小行星是否满足保守撞击条件，如果不满足保守撞击条件，将小行星撞击地球概率判定为0，排除蒙特卡洛打靶中不能够撞击地球的小行星轨道预测，节省小行星撞击概率评估运算资源，提高小行星撞击概率评估效率。
[0067]
5、本发明公开的基于引力锁孔面积计算的小行星撞击概率评估方法，针对潜在危险小行星轨道范围特点，考虑中心天体引力、第三体引力摄动、太阳光压摄动、雅克夫斯基效应摄动，构建小行星高精度轨道动力学模型，能够提高小行星高精度轨道动力学模型预测精度，进而提高小行星撞击概率评估精度。
[0068]
6、本发明公开的基于引力锁孔面积计算的小行星撞击概率评估方法，评估的小行星撞击地球概率能够作为行星防御工程实施的判据，当小行星撞击地球概率超过阈值时，开启行星防御工程，避免小行星撞击地球。
附图说明
[0069]
图1基于引力锁孔面积计算的小行星撞击概率评估方法流程图。
[0070]
图2撞击偏差示意图。
[0071]
图3引力锁孔面积示意图。
[0072]
图4小行星a2轨道参数。
[0073]
图5小行星a2轨道图。
[0074]
图6撞击时刻4月1日引力锁孔面积累积分布函数。
具体实施方式
[0075]
为了更好的说明本发明的目的和优点，下面结合附图和实例对发明内容做进一步说明。
[0076]
实施例1：
[0077]
如图1所示，本实例公开的一种小行星撞击概率评估方法，具体实现步骤如下：
[0078]
步骤一、针对潜在危险小行星轨道范围特点，考虑中心天体引力、第三体引力摄动、太阳光压摄动、雅克夫斯基效应摄动，构建小行星高精度轨道动力学模型。
[0079]
本实例假设在2022年4月1日观测到小行星轨道参数和轨道图分别如图4和图5所示，可能导致下一次2026年4月1日小行星再次近距离飞越地球时撞击地球。潜在危险小行星日心轨道主要位于地球至木星轨道范围，主要考虑受到地球、火星和木星带来的第三体引力摄动。第三体引力摄动p
third
表达式为
[0080][0081]
其中，μ
pl
表示第三体的引力常数，r
pl
表示第三体相对中心天体位置矢量，δ表示第三体相对小行星的位置矢量。
[0082]
小行星轨道动力学方程为
[0083][0084]
其中，μ为太阳引力常数，为轨道面内与位置矢量垂直指向速度方向的横向单位矢量。
[0085]
步骤二、测得ti时刻(即2026年4月1日)小行星的轨道位置矢量和速度矢量，并添加观测误差，得到用于轨道预测的小行星轨道初值，基于步骤一构建的小行星高精度动力学模型进行蒙特卡洛打靶实现小行星轨道预测，并记录下小行星飞越地球位置速度矢量，根据位置速度矢量得到近心点半径和双曲超速参数；将所述近心点半径和双曲超速参数代入状态转移矩阵，利用状态转移矩阵推导小行星速度改变量与小行星-地球交会距离的关系，并构建保守撞击条件。判定小行星是否满足保守撞击条件，如果不满足保守撞击条件，将小行星撞击地球概率判定为0，排除蒙特卡洛打靶中不能够撞击地球的小行星，节省小行星撞击概率评估运算资源，提高小行星撞击概率评估效率；如果满足保守撞击条件，转移至步骤三。
[0086]
在ti时刻测得小行星的轨道位置矢量和速度矢量分别为rd和vd，将观测误差三维随机矢量ξ、三维随机矢量ζ、三维随机矢量η代入rd和vd，得到用于轨道预测的小行星轨道初值ri和vi分别为
[0087][0088]
其中，三维随机矢量ξ各分量分别为服从均值为0标准差为位置精度σr的随机数，三维随机矢量ζ各分量分别为服从均值为0标准差为速度精度σv的随机数，三维随机矢量η为导航误差导致的速度误差，其各分量分别为服从均值为0标准差为σg的随机数。其中，σr＝100km，σv＝1m/s且3σg＝20m/s。
[0089]
以式(21)为初值带入模型(20)进行积分，当小行星在tk时刻前后近距离飞越地球，即小行星轨道状态与地球轨道状态re,ve满足式(22)时，记录下该时刻和对应的小行星的轨道状态
[0090][0091]
此时，经过式(23)至式(25)的推导，小行星相对地球的近心点半径和双曲超速大小根据式(26)、(27)计算：
[0092][0093][0094][0095][0096][0097]
其中，为μe为地球引力常数，||
·
||表示二范数。
[0098]
计算保守撞击条件涉及2条轨道，如图2所示，一条为标称撞击转移轨道，即在精确动力学模型下求解的从tk时刻地心位置到t
t
时刻地心位置的转移轨道；另一条轨道为实际小行星飞行轨道。其中，δr1,δv1分别为tk时刻实际小行星轨道与标称撞击轨道的位置、速度偏差，δr2,δv2分别为t
t
时刻实际小行星轨道与标称撞击轨道的位置、速度偏差。若实际小行星轨道与标称撞击转移轨道偏差δr1,δv1太大，则仅靠自然摄动力无法提供足够的轨道改变量使小行星轨道足够靠近t
t
时刻的地球导致撞击，只有当该小行星轨道与撞击轨道的偏差满足保守撞击条件时，该小行星才有撞击地球的可能性。
[0099]
已定义小行星轨道与撞击转移轨道的首、末偏差分别为δr1,δv1、δr2,δv2，由状态转移矩阵关系
[0100][0101]
得到
[0102][0103]
在日心系下，|δv1|即为小行星轨道与标称撞击轨道的双曲超速矢量之差
[0104]
|δv
∞
|，δr1等效为0，|δr2|为地球半径。因此式(29)等效为
[0105][0106]
由于式(30)为近似线性关系，出于保守估计考虑，当误差处于数倍e之内，仍将其考虑为可能的撞击情况，因此得保守撞击条件为
[0107][0108]
对于进入地球引力影响范围的小行星轨道，将式(26)、式(27)所计算数值，代入式(31)所示的保守撞击条件，只有当式(13)成立时才考虑该小行星在t
t
时刻有撞击地球的可能性，将不满足条件的轨道判定在t
t
时刻撞击地球的可能性为0，排除蒙特卡洛打靶中不能够撞击地球的小行星，节省小行星撞击概率评估运算资源；如果满足保守撞击条件，转移至步骤三。
[0109]
步骤三、利用变分法推导近心点半径偏差量和双曲超速方向角度偏差量的表达式，代入引力锁孔面积定义式得到引力锁孔面积计算公式；将步骤二得到的小行星相对地球的近心点半径和双曲超速大小代入所述引力锁孔面积计算公式得到引力锁孔面积。
[0110]
小行星的引力锁孔面积(如3所示)定义式如式(32)所示。
[0111][0112]
式(32)中r即为小行星相对地球的双曲线轨道的近心点半径r
p
。利用变分法推导得到如式(33)所述近心点半径偏差量和式(34)所述双曲超速方向角度偏差量，
[0113][0114][0115]
其中|δv1|、δf
∞
和δv
∞
通过已经完成的实际小行星轨道数值结果和标称撞击轨道之差得到。
[0116]
将式(33)、式(34)代入引力锁孔面积定义式(32)得到引力锁孔面积计算公式(35)；将步骤二得到的小行星相对地球的近心点半径和双曲超速大小带入所述引
力锁孔面积计算公式得到引力锁孔面积。
[0117][0118]
步骤四、根据步骤二中蒙特卡洛打靶记录及保守撞击条件判定统计出满足保守撞击条件的预测的小行星轨道数目；然后根据步骤三中计算的引力锁孔面积数值计算引力锁孔面积的累积分布函数，并根据任务需求预设置信度选取引力锁孔面积估计值；最后根据满足保守撞击条件的小行星轨道预测数目占步骤二中所有蒙特卡洛打靶的小行星轨道预测数目比例和引力锁孔面积估计值，计算小行星在步骤二中给定t
t
时刻撞击地球的概率。
[0119]
步骤二蒙特卡洛打靶的k＝500条小行星轨道预测中，共记录k＝454条小行星轨道预测在tk时刻进入地球引力影响范围a
sphere
＝π(106)2km2，所述454条轨道在步骤三中有n＝11条轨道判定为满足式(13)所述的保守撞击条件。
[0120]
所述11条小行星轨道预测在步骤三中已分别计算出对应的引力锁孔面积。所计算引力锁孔面积略有不同，利用概率论基础理论计算所计算引力锁孔面积的累积分布函数(cdf)如图6，并根据任务需求预设置信度cdf＝99％，取置信度对应的引力锁孔面积数值作为该小行星在t
t
时刻撞击地球所对应的tk时刻引力锁孔面积估计值为2.4615
×
10
10
km2。则在99％的置信度下，小行星在2026年4月1日撞击地球概率计算为
[0121][0122]
表明，小行星a2在2026年4月1日撞击地球的概率有99％置信度不大于1.7237
×
10-4
。
[0123]
步骤五：根据步骤四得到的小行星撞击地球概率能够作为行星防御工程实施的判据，当小行星撞击地球概率超过阈值时，实施行星防御，避免小行星撞击地球。
[0124]
以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。