一种求解撞击地球小天体的物理特性的方法

1.本发明属于小行星防御与利用领域，具体涉及一种求解撞击地球小天体的物理特性的方法。


背景技术：

2.小天体撞击地球的频率与尺寸呈负相关。小尺寸的小天体撞击地球频率高，但其较小的体积与较低的反照率使他们更加难以观测。同时，小尺寸小天体由于其撞击地球带来的危害效应相对轻微，往往并不是天文台或望远镜的观测目标。然而，1908年通古斯大爆炸与2013年车里雅宾斯克撞击事件表明，即使尺寸在几十米甚至几米的小天体，在撞击地球大气层过程中也能够产生剧烈的冲击波，从而可能导致严重的人员伤亡与经济损失。因此，中小尺寸(直径1至140米)的小天体正受到越来越广泛的关注与重视。
3.2022年我国航天局明确表示中国将建设近地小天体防御系统。为更好地了解小天体的撞击威胁，评估小天体撞击效应，有必要了解小天体(尤其是撞击频率高的中小尺寸小天体)的物理性质。小天体的物理特性主要指其结构、材质、密度、尺寸等性质。在研究中，结构通常由小天体的解体强度来反映。
4.获取撞击地球的小天体物理特性的方式主要有三种。
5.第一种为对小行星/彗星进行采样或收集陨石，直接获取小天体信息。目前人类已经实现了对“贝努”小行星和“龙宫”小行星的采样返回，获得了珍贵的材料样本。此外，若小天体的尺寸在米级以上，则很有可能在地球留下陨石。陨石或小天体岩石样本中包含小天体的结构信息，如小天体的材质、密度等，也能够提供太阳系早期演化的线索，是珍贵的研究资料。然而，由于小天体种类繁杂数量庞大，不同小天体之间结构差异也可能很大，而通过采样或陨石收集得到材料的小天体数量远不及日常撞击地球的小天体数量。因此，该方法虽然直接，但具有数量和统计上的局限性。
6.第二种方式为光学观测，主要指对小天体的天文学观测。通过这种方式，研究者能够获取或部分获取小天体的绝对星等、反照率、自转速度、近红外光谱等数据。这些数据可进一步用于推断小天体的尺寸、表面材质、内聚强度。然而，这类方法通常用于直径千米级的小天体物理特性的估算，很少应用于数量庞大且撞击地球频率更高的小尺寸小天体。
7.这首先是因为小尺寸小天体的观测数据少。截至2022年5月23日，在撞击地球前被人类观测到的小天体只有5个，他们分别为2008tc3，2014aa，2018la，2019mo和2022 eb5。这意味着在可预见的未来，人类对可能撞击地球的小天体的观测仍有待完善。第二个原因是利用该方法估算小尺寸小天体性质，其相对误差比千米及小天体的误差要大。相较于小天体的材质与反照率，绝对星等更容易获得，此时通常会假设小天体的反照率为平均反照率0.147，借助绝对星等、反照率与尺寸的经验关系估算小天体尺寸。该方法用于估算千米级小天体的尺寸，相对误差较小，但估算小尺寸小天体，则会导致较大误差。
8.另外一种方法是通过观察小天体撞击地球过程中与大气层的相互作用确定小天体的物理特性，如小天体在大气进入过程中的解体现象、减速与烧蚀过程、发光现象等。全
球很多地区已经建立了火流星观测网络。这些观测网通常为相机网络，能够在各自的视场观测并记录撞击地球的小天体在大气内发光飞行阶段的轨迹。对光辐射过程进行分析，能够还原小天体飞行的能量沉积过程。此外，nasa的近地小天体研究中心(cneos)自1988年起，通过卫星传感器记录并发布小天体在峰值亮度时刻的位置、速度及估算的总能量，且近几年已经实现了全球覆盖，已经成为确定小天体物理性质非常重要的数据来源。
9.然而，尽管小天体发光飞行阶段的轨迹已经被记录下来，但它的质量、尺寸及解体强度等信息仍然无法准确获知。前人的研究中，通常将小天体峰值亮度时刻的前缘动压作为小天体的解体强度，并假设小天体进入大气至达到峰值亮度期间，大气层对小天体的速度没有影响，再根据估算的总动能进一步计算小天体的质量。通过对小天体密度进行假设，给出小天体的名义直径。通过此种方法计算的小天体解体强度往往比小天体的真实解体强度大很多，小天体材质未知的情况下计算的名义直径也可能与真实的尺寸差别较大。
10.综合以上分析，能够发现，获取中小尺寸小天体的物理特性具有非常重要的意义；而对于撞击地球频率高的中小尺寸小天体，尚缺乏一种能够大规模地准确地推算小天体物理特性的方法。


技术实现要素：

11.本发明的目的在于克服现有计算方法结果与真实值偏差大的缺陷。
12.为了实现上述目的，本发明提出了一种求解撞击地球小天体的物理特性的方法，所述方法为：根据小天体进入大气后的物理数据，结合煎饼模型与遗传算法，考虑大气层对小天体各物理量以及轨道的影响，计算小天体的物理特性，利用已知小天体类型和物理特性对照表对计算结果进行筛选，从而得到待求解小天体的物理特性；其中，小天体的物理特性包括：质量、半径、密度和解体强度。
13.作为上述方法的一种改进，所述方法具体包括：
14.步骤1：获取小天体进入大气信息；获取的信息包括峰值亮度时刻、峰值亮度时刻位置、峰值亮度时刻速度和撞击总能量；
15.步骤2：利用小天体峰值亮度时刻的物理量，通过四阶龙格-库塔方法对小天体飞行方程、烧蚀方程和半径变化方程构成的微分方程组进行逆向数值积分，得到小天体进入大气层边界处的速度矢量位置矢量半径r0、半径变化率dr0和质量m0；通过遗传算法计算小天体的解体强度s0、小天体在峰值亮度时刻的半径r
p
、半径变化率dr
p
和质量m
p
；
16.步骤3：根据小天体类型与解体强度的对应关系，筛选出可能正确解。
17.作为上述方法的一种改进，所述飞行方程为：
[0018][0019]
其中，为小天体相对地球的矢径；r为的模，是小天体到地球的距离；m为小天体的质量；μ为地球引力常数，取3.986
×
10
14
m3/s2；t为时间；表示小天体受到的大气阻力，满足下式：
[0020]
[0021]
其中，cd为阻力系数，取值为1；vr为小天体相对地球大气层的速度大小，亦为小天体在ecef坐标系下的速度大小；s为迎风截面积；为单位矢量，表示的方向；ρa为大气密度。
[0022]
作为上述方法的一种改进，所述烧蚀方程为：
[0023][0024]
其中，m为小天体的质量；t为时间；ch为烧蚀系数；qh为烧蚀热；ρa为大气密度；vr为小天体相对地球大气层的速度大小，亦为小天体在ecef坐标系下的速度大小；s为迎风截面积
[0025]
作为上述方法的一种改进，所述半径变化方程为：
[0026][0027]
其中，r为小天体的半径；t为时间；ρa为大气密度；vr为小天体相对地球大气层的速度大小，亦为小天体在ecef坐标系下的速度大小；ρm为小天体的密度。
[0028]
作为上述方法的一种改进，所述逆向数值积分计算公式为：
[0029][0030]
其中，xb表示逆向数值积分结果；h
p
为小天体峰值亮度时刻的高度；hb为大气层边界高度hb＝100km；h为小天体高度；t为时间；x代指小天体被积分的物理量，包括矢径速度矢量质量m，半径r，半径变化率
[0031][0032]
逆向数值积分过程由四阶龙格-库塔方法实现，最终得到小天体在大气层边界处eci坐标系下的速度矢量位置矢量半径r0、半径变化率dr0和质量m0。
[0033]
作为上述方法的一种改进，所述逆向数值积分时，未知变量包括：小天体峰值亮度时刻的质量m
p
、半径r
p
及半径变化率dr
p
，小天体的密度ρm、解体强度s0、烧蚀热qh，模型参数n；
[0034]
其中，n取为4，当n取4无解时，n取2；
[0035]
根据彗星、c类小行星、s类小行星、m类小行星四类小天体的类型与密度ρm、烧蚀热qh的对应关系，取四组密度ρm、烧蚀热qh分别求解；
[0036]
彗星对应密度ρm为1000kg/m3，对应烧蚀热qh为2.5
×
106j/kg；
[0037]
c类小行星对应密度ρm为2000kg/m3，对应烧蚀热qh为5
×
106j/kg；
[0038]
s类小行星对应密度ρm为3500kg/m3，对应烧蚀热qh为8
×
106j/kg；
[0039]
m类小行星对应密度ρm为7900kg/m3，对应烧蚀热qh为8
×
106j/kg；
[0040]
另外四个未知量s0、m
p
、r
p
及dr
p
，作为优化参数，利用遗传算法进行求解：
[0041]
遗传算法的四个约束条件分别为：
[0042]
(1)计算得到的大气层边界处的半径变化率为0，即dr0＝0；
[0043]
(2)计算得到的大气层边界处的半径与质量满足球体假设，即r(m)-r0＝0，其中
[0044]
(3)计算得到的大气层边界处的动能与获取的小天体进入大气信息中的撞击总能量相等，即e(v0,m0)-e0＝0，其中
[0045]
(4)计算得到的大气层边界处半径应与峰值亮度时刻半径r
p
满足煎饼模型的比例关系，即nr
0-r
p
＝0；
[0046]
其中，e0为获取到的总撞击能量；m0为大气层边界处小天体质量；r0为大气层边界处小天体半径；ρ
ap
为小天体峰值亮度时刻大气密度；v
p
为小天体峰值亮度时刻速度；
[0047]
遗传算法中四个优化参数s0、m
p
、r
p
及dr
p
取值范围为：
[0048]
s0最大值为p
p
，最小值为p
p
/30，
[0049]mp
最大值为2me，最小值为0.01me，
[0050]rp
最大值为nr
0max
，最小值为r
0min
，，
[0051]
dr
p
最大值为50n，最小值为0.1；
[0052]
遗传算法中目标函数f由约束(1)至(4)加权组合而成，公式如下：
[0053][0054]
其中，a1、a2、a3、a4为四个约束的权重，初始值为10、100、0.01、50；
[0055]
求得的解是否符合要求的判断标准为以下四项全部满足：
[0056]
|r
0-r(m0)|/r0＜1％
[0057]
|nr
0-r
p
|/r
p
＜1％
[0058]
|m
0-m(e0,v0)|/m0＜0.1％
[0059]
|dr0|＜0.001
[0060]
通过遗传算法最终计算得到小天体的解体强度s0、小天体在峰值亮度时刻的半径rp
、半径变化率dr
p
和质量m
p
。
[0061]
作为上述方法的一种改进，根据计算出的小天体的物理特性和状态向量可以计算小天体日心轨道的轨道根数，包括：径向速度、比角动量、轨道倾角、节点线向量、升交点赤经、偏心率、近地点俯角、真近点角和轨道半长轴。
[0062]
作为上述方法的一种改进，求出小天体轨道后，可根据小天体相对木星的蒂塞朗常数tj，远日点距离q、近日点距离q及木星半长轴aj的距离之间的关系，对小天体的类型进行验证。
[0063]
作为上述方法的一种改进，如果已知小天体的绝对星等，通过式(20)根据求得的直径，反推小天体的反照率，进一步根据不同类别小天体的反照率对小天体材质进行验证或约束：
[0064][0065]
其中，h为小天体的绝对星等，d为小天体的直径，pv为反照率。
[0066]
与现有技术相比，本发明的优势在于：
[0067]
1、本方法结合煎饼模型与遗传算法，在判断小天体类型、计算小天体解体强度、反照率、直径等主要物理参数方面具有较大优势。
[0068]
2、本方法在飞行方程中考虑了大气层对小天体飞行的影响，计算的轨道根数比不考虑大气层的轨道根数更加准确。
附图说明
[0069]
图1所示为求解撞击地球小天体的物理特性的方法流程图；
[0070]
图2所示为煎饼模型示意图；
[0071]
图3所示为2022 eb5反照率与直径关系图。
具体实施方式
[0072]
本发明旨在利用小天体与地球大气层间的相互作用，对小天体的物理特性进行求解或约束，同时由此计算小天体的日心轨道根数。由于nasa近地小天体研究中心(cneos)能够记录大部分米级及以上小天体在峰值亮度时刻的速度、位置并根据经验关系估算撞击总能量，本发明以覆盖率广的cneos火球数据为基础，以使本方法具有更高的普遍性。本发明利用半解析解体模型，充分考虑大气层对小天体速度的影响，对小天体的物理特性进行推断。对于未知物理量，本发明将其视为优化参数，采用遗传算法对未知数进行求解。
[0073]
下面结合附图对本发明的技术方案进行详细的说明。
[0074]
本发明利用cneos提供的小天体撞击地球过程中峰值亮度时刻的速度、位置及撞击能量，推算小天体的物理特性。整体计算流程如图1所示。该过程中使用的数据与方法如下：
[0075]
1.获取小天体进入大气信息
[0076]
cneos公布的撞击地球的小天体信息列表如表1所示，其中v
x
，vy，vz分别表示小天体的速度矢量在地心固连(ecef)坐标系下x、y、z方向下的分量，vr表示ecef系下速度矢量的模，即速度的大小。
[0077]
表1小天体峰值亮度时刻的速度、位置及撞击能量
[0078][0079]
表1中的数据将进一步转换为地心惯性(eci)坐标系，得到eci系下小天体峰值亮度时刻的状态向量，即位置矢量r
p
，速度矢量v
p
，以及撞击总能量，也为小天体在地球大气层边界处的动能。
[0080]
2.小天体进入大气过程飞行方程
[0081]
小天体进入大气的速度理论上在11.2km/s至72km/s之间，因此小天体与大气层之间的撞击为超高速撞击。在大气进入阶段，小天体的速度将在大气阻力的作用下降低，即减速；同时与大气之间剧烈的相互作用使小天体表面温度升高，导致小天体的质量损失，又称为烧蚀。
[0082]
减速由小天体的飞行方程控制，如式(1)所示。
[0083][0084]
其中，为小天体相对地球的矢径，r为的模，是小天体到地球的距离，m为小天体的质量；μ为地球引力常数，取3.986
×
10
14
m3/s2。表示小天体受到的大气阻力，由式(2)决定：
[0085][0086]
上式中，cd为阻力系数，在本计算中按球体假设取为1；vr为小天体相对地球大气层的速度大小，亦为小天体在ecef坐标系下的速度大小。s为迎风截面积，为单位矢量，表示的方向。ρa为大气密度，根据wheeler等人对1976美国标准大气模型的拟合公式计算：
[0087]
ρa＝-140.2e-0.000187h
+141.4e-0.000186h
ꢀꢀꢀ
(3)
[0088]
其中，h为海拔高度，单位为m。
[0089]
小天体的质量损失由烧蚀方程控制，如式(4)所示：
[0090][0091]
其中，ch为烧蚀系数，取0.1；qh为烧蚀热，其取值与小天体类型有关，后续将在介绍
小天体分组时统一赋值。
[0092]
3.煎饼模型
[0093]
随着小天体进入地球大气层，小天体周围的大气密度也随高度的降低而近似指数增加，小天体的前缘所受的动压p＝ρ
avr2
也随之增加。通常认为，当p超过小天体解体强度s0时，小天体会发生解体。
[0094]
解体是小天体进入大气过程中非常重要的物理现象，也是剧烈的空爆现象的直接诱因。本发明使用煎饼模型对小天体的解体现象进行模拟，其基本原理示意如图2所示。
[0095]
煎饼模型认为，小天体解体后，大大小小的碎片构成一个碎片云团。在解体初期，该云团仍然受到同一个激波的作用，因此在这一阶段依然可以将其视为一个整体。随后，各碎片逐渐分离，彼此的间距变大，整体表现为云团的迎风截面半径增加，该过程中，半径r的变化由式(5)求得：
[0096][0097]
其中，ρm为小天体的密度。
[0098]
当云团的半径增加至初始半径的n倍时，各碎片发展出各自的激波，不能再作为一个整体进行计算。n为经验参数，是小天体半径变化程度的度量，在模型提出时取2，但在应用中发现，n通常取4-8才能够实现对解体过程的较好拟合。此时，小天体的减速与烧蚀作用迅速加剧，能量在短时间内大量释放，导致光辐射升高达到峰值，同时产生剧烈的冲击波，该现象即为空爆。从描述中可知，空爆与峰值亮度同步发生。
[0099]
4.未知变量的求解
[0100]
利用小天体峰值亮度时刻的物理量，通过对小天体飞行方程、烧蚀方程和半径变化方程组成的微分方程组进行逆向数值积分，可得到小天体进入大气层边界处的物理量，包括：，速度矢量位置矢量半径r0、半径变化率dr0和质量m0。
[0101]
逆向数值积分过程从小天体峰值亮度时刻的高度h
p
向上积分至大气层边界处(大气层边界高度为hb＝100km)：
[0102]
其中x代指小天体被积分的物理量，包括矢径速度矢量(即)，质量m，半径r，半径变化率
[0103]
各物理量的微分方程表达如式(1)至式(5)所示。
[0104]
逆向积分过程由四阶龙格-库塔方法实现。该方法为求解非线性常微分方程的常用方法，利用微分方程与初始条件，通过逐步迭代实现高精度的求解。
[0105]
对飞行过程进行求解，会发现，仍有一些小天体的物理量未知：小天体峰值亮度时刻的质量m
p
、半径r
p
及半径变化率dr
p
，小天体的密度ρm、解体强度s0、烧蚀热qh，模型参数n。
[0106]
n的取值范围为2-8，对于小尺寸的小天体而言，n的取值不宜过大，因此本方法中n取为4，当n取4无解时，n取2。密度ρm与烧蚀热qh均与小天体类型有关，对于未捡到陨石的小天体而言，其类型往往未知，但它与小天体的解体强度存在定性的关联：通常情况下，彗星、c类小行星、s类小星星、m类小行星的密度递增，解体强度也递增。因此，本方法中，首先将计算按材质分为四组，每组的小天体密度均取该类型下典型的小天体密度值。四组均求得解后，再按照求得的解体强度对小天体的类型进行约束，并求出最终的小天体类型。按小天体类型的分组情况如表2所示：
[0107]
表2小天体类型、密度、烧蚀热
[0108][0109]
小天体类型与材质存在对应关系，彗星主要材质为水冰，c类、s类、m类主要对应碳质、石质、铁质。
[0110]
另外四个未知量s0、m
p
、r
p
及dr
p
，将作为优化参数，利用遗传算法进行求解。
[0111]
遗传算法是一种通过模拟自然选择以求问题最优解的优化方法，其工作流程为：
[0112]
(1)在优化变量的取值范围内生成初始种群，优化变量包括小天体的解体强度、小天体峰值亮度时刻的质量、半径、半径变化率；
[0113]
(2)计算目标函数值：利用四阶龙格库塔方法，对小天体飞行方程、半径变化方程、烧蚀方程进行逆向数值积分，计算小天体在大气层边界处的速度位置矢量、质量、半径、半径变化率。并根据边界处的物理量计算目标函数值。
[0114]
(3)判断是否满足终止条件，若是，则输出优化结果；若否，则对初代种群进行遗传、交叉、变异，产生新一代个体，重复(2)和(3)。
[0115]
在该问题中，存在四个独立的约束条件，未知量与条件数量相等，能够保证存在唯一解。四个约束条件分别为：
[0116]
(1)计算得到的大气层边界处的半径变化率为0，即dr0＝0；
[0117]
(2)计算得到的大气层边界处的半径与质量满足球体假设，即r(m)-r0＝0，其中
[0118]
(3)计算得到的大气层边界处的动能与cneos估算的撞击总能量相等，即e(v0,m0)-e0＝0，其中
[0119]
(4)计算得到的大气层边界处半径应与峰值亮度时刻半径r
p
满足煎饼模型的比例关系，即nr
0-r
p
＝0。
[0120]
上述描述中，除能量e0外，下标0表示通过对飞行方程逆向积分求得的大气层边界处的物理量。e0为cneos估算得到的总撞击能量。下标p表示峰值亮度时刻的物理量。
[0121]
该求解过程由遗传算法实现。遗传算法是一种模拟自然选择以求解最小值问题的最优化算法。该算法首先从优化参数的取值范围中选取一定数目的个体，利用适应度函数或目标函数的取值作为个体的性状，并对个体的性状是否符合进化目标进行评估。适应性高的个体将有机会将基因传递给下一代，而性状不够优秀的个体将被淘汰。缺失的个体位将继续在参数的取值范围内随机选取。如此产生一代代个体，直至达到最大代数或达到进化要求。
[0122]
本方法中，四个优化参数s0、m
p
、r
p
及dr
p
取值范围的选取在充分考虑撞击事件个别性的基础上设定，其取值如表3所示。
[0123]
表3优化变量取值范围
[0124][0125]
目标函数f由约束(1)至(4)加权组合而成，如式(6)所示：
[0126]
[0127]
其中，a1、a2、a3、a4为四个约束的权重，初始值为10、100、0.01、50以使四个约束值维持在同一量级，也能够根据撞击事件的实际情况进行调解。同时，求得的解是否符合要求的判断标准为：|r
0-r(m0)|/r0＜1％，|nr
0-r
p
|/r
p
＜1％，|m
0-m(e0,v0)|/m0＜0.1％，|dr0|＜0.001。
[0128]
通过(4)中方法求出未知量后，小天体大气层边界处的状态向量(位置矢量和速度矢量)也可求得。利用小天体的解体强度与材质之间的关系，可对小天体的材质进行判断。
[0129]
5.轨道计算
[0130]
经过对未知变量的求解，本发明中的核心计算已经完成，但为了使物理特性的推断更为全面，本方法将继续计算小天体日心轨道的轨道根数。对于有轨道观测数据的小天体，计算轨道可以验证本发明中方法的可靠性；对于缺乏轨道观测数据的小天体，计算轨道能够从轨道的统计规律角度对小天体类型(仅判断为彗星/小星星)进行验证。
[0131]
获得小天体大气层边界处eci坐标系下的状态向量后，需要进一步对轨道方程进行逆向积分，求得小天体在地球影响球边界处的eci系下的状态向量。进一步可将状态向量转换为小天体在日心黄道坐标系下的轨道根数。具体步骤如下：
[0132]
(1)计算小天体在地球影响球边界处eci坐标系下的状态向量：
[0133]
该过程需要对轨道方程进行逆向积分，直至小天体与地球的距离达到地球影响球边界处。该阶段的轨道方程可以不考虑地球大气的影响，因为没有大气阻力项，如式(7)所示：
[0134][0135]
地球影响球半径
[0136]
(2)将eci坐标系下的状态向量转变为地心黄道坐标系下的状态向量：
[0137][0138]
其中，与为小天体在地球影响球边界处、地心黄道系下的状态向量，为小天体在地球影响球边界处、地心黄道系下的状态向量，为小天体在地球影响球边界处、eci坐标系下的状态向量；ε为黄赤交角，取23.43929
°
。
[0139]
(3)计算小天体在日心黄道坐标系下的状态向量：
[0140]
日心黄道坐标系与地心黄道坐标系均为惯性系，且坐标平面一致，因此只需简单的计算即可求得：
[0141][0142]
其中，与为小天体日心黄道坐标系下的状态向量，为地球在日心黄道坐标系下的状态向量。
[0143]
(4)计算小天体日心黄道坐标系下的轨道根数
[0144]
为方便表示，步骤(4)中日心黄道系下的以和表示。
[0145]
计算径向速度vr：
[0146][0147]
计算比角动量
[0148][0149]
计算轨道倾角i：
[0150][0151]
计算节点线向量
[0152][0153]
计算升交点赤经ω：
[0154][0155]
计算偏心率e：
[0156][0157]
计算近地点俯角ω：
[0158][0159]
计算真近点角θ：
[0160]
[0161][0162]
计算轨道半长轴a：
[0163][0164]
6.小天体类型验证
[0165]
当求出小天体轨道后，可根据小天体相对木星的蒂塞朗常数tj，远日点距离q、近日点距离q及木星半长轴aj的距离之间的关系，对小天体的类型(仅能判断彗星/小行星)进行验证。类型判断标准如表4所示。
[0166]
表4彗星、小行星的轨道统计规律
[0167][0168]
蒂塞朗常数计算方式如式(19)所示：
[0169][0170]
需要注意的是，该准则为统计规律，且仅能对小天体属于彗星或小行星做出大致判断，因此其分类无法作为判断小天体类型的绝对标准。
[0171]
7.绝对星等
[0172]
小天体的绝对星等h、直径d与反照率pv之间存在经验关系，如式(20)所示：
[0173][0174]
在利用遗传算法及煎饼模型对小天体大气进入过程进行求解后，最终得到四组小天体物理参数，通过对比解体强度与材质是否符合经验，可以对不合理的几组进行排除。然而，由于c类小行星与s类小行星的强度界限较为模糊，因此最终将剩余一或两组数据。
[0175]
此时，若小天体的绝对星等已知，那么可以通过式(20)可以根据求得的直径，反推小天体的反照率，进一步根据不同类别小天体的反照率对小天体材质进行验证或约束。c类小行星的反照率范围通常为0.03-0.11，s类小行星的反照率范围约为0.13-0.29，m类小行星的反照率约为0.1-0.3.
[0176]
2022年3月11日，匈牙利天文学家观测到一颗编号为2022 eb5的小天体将撞击地球。在被发现的两小时后，该小天体进入地球大气层，并在挪威海上空发生空爆，成为第五颗在撞击地球前被观测到的小天体。
[0177]
为展示本技术的技术效果，同时验证本技术的正确性，现采用本技术中的方法对2022 eb5的物理特性进行分析，并对其日心轨道进行计算。
[0178]
1.观测信息
[0179]
1.1轨道观测信息
[0180]
本方法在计算过程中，并不会使用到小天体的轨道观测信息，但已知轨道观测信
息能够从侧面对本方法的正确性与合理性进行验证与说明，有利于技术效果的展示。
[0181]
2022 eb5撞击地球前的轨道被多处天文台观测到，nasa的jpl实验室与esa根据不同天文台的观测数据，分别计算了2022 eb5的日心轨道根数，如表5和表6所示。
[0182]
表5 2022 eb5轨道根数(jpl解)
[0183][0184]
表6 2022 eb5轨道根数(esa解)
[0185][0186]
其中，q为近日点距离，q为远日点距离，m为平近点角。
[0187]
可见，jpl与esa的轨道解存在较大差异。在呈现本方法计算结果时，将利用本方法计算得到的轨道根数与以上两组轨道解进行对比，可以对两组解的可靠性进行判断。
[0188]
1.2cneos峰值亮度观测信息
[0189]
cneos的观测数据为小天体进入大气过程中峰值亮度时刻的信息，为本方法计算中使用的主要数据。
[0190]
2022 eb5大气进入过程中，其峰值亮度时刻的速度、位置、时间及其撞击总能量如表7所示。
[0191]
表7 2022 eb5的峰值亮度观测数据
[0192][0193]
表7中各参数含义与表1相同。
[0194]
2.计算结果
[0195]
2.1 2022 eb5物理特性分析
[0196]
为更好地展示本方法的技术效果，本节首先将讨论在不使用本方法计算的情况下，利用已知信息能够对小天体的物理特性进行何种约束。
[0197]
根据轨道根数，可以计算2022 eb5相对于木星的蒂瑟朗常数，t＝2.89。远日点为4.77au，大于木星半长轴。根据蒂瑟朗常数判据，2022 eb5很可能为彗星小天体，但根据远日点距离、木星半长轴之间的关系判断，2022 eb5大概率为小行星。因此，根据轨道难以判断其类型。
[0198]
对2022 eb5的观测给出了其绝对星等。根据式(20)可以得到2022 eb5直径与反照率之间的关系，如图3所示。然而，由于不确定2022 eb5的小天体类型，其反照率信息未知。此时，通常假设其反照率为平均反照率0.147，在此假设下，计算得到的2022 eb5直径为1.69m。
[0199]
根据cneos对小天体峰值亮度时刻的观测数据，按照以往的估算方法，假设小天体进入大气时刻的速度与峰值亮度时刻的速度相等，即为17.24km/s，则根据初始动能，得到小天体质量约112640kg。按该估算方法，小天体解体强度约为峰值亮度时刻的总压，即为3.5mpa。
[0200]
以上为不使用本方法，按以往常规处理方法估算得到的小天体物理特性。
[0201]
使用本方法，考虑大气层对小天体速度、质量的影响，得到四组小天体大气层边界处速度、质量及小天体解体强度的解，如表8所示。其中，第四列为各自类型的小天体解体强度的经典范围。
[0202]
表8 2022 eb5物理特性及大气层边界处物理参数表
[0203][0204]
从表8中第三列与第四列可以看到，根据考虑大气层得到的解体强度，2022 eb5不可能为彗星与m类小行星，其最有可能为c类小行星，然而，仅靠强度不能完全排除其为s类小行星的可能性。此外，考虑大气层的影响计算得到的解体强度(1.4-2.1mpa)，明显低于不考虑大气层(即不使用本方法)得到的解体强度3.5mpa，二者相对差高达40％-60％。
[0205]
当2022 eb5为c类或s类时，其直径范围为3.91-4.68米，与不使用本方法计算得到的直径大1.3-1.76倍。利用绝对星等数据与图3，可以得到2022 eb5的反照率约为0.01-0.035。该范围与c类小行星的反照率最为符合。(c类、s类、m类反照率范围分别为0.03-0.11，0.13-0.39，0.1-0.3)。因此，2022 eb5极有可能为一颗c类小行星，且其反照率非常低。
[0206]
此外，考虑到已知密度的c类小行星——“龙宫”与“贝奴”小行星，其密度均为～1200kg/m3，均小于表8中设定的2000kg/m3，因此本节按照以上两颗小行星的密度增加两组计算结果，如表9所示。
[0207]
表9 2022 eb5物理特性及大气层边界处物理参数表
[0208][0209][0210]
其中，c类1的密度设定与“龙宫”小天体一致，c类2的密度设定与“贝奴”小天体一致。
[0211]
2.2 2022 eb5轨道计算结果
[0212]
将大气层边界处的状态向量转换为日心轨道，其结果如表10所示。
[0213]
表10使用本方法计算得到的2022 eb5日心轨道根数
[0214][0215]
表10中，组别中“jpl”与“esa”分别表示jpl与esa根据轨道观测数据计算的2022 eb5的日心轨道参数(表5与表6)，为1组；c类1、c类2、c类三组为本方法计算得到的日心轨道参数，为2组；无大气为不考虑大气层的影响直接计算得到的日心轨道参数，为3组。比较三个组别的半长轴a可以发现，2组计算得到的半长轴更接近jpl计算的半长轴，二者最低相差～0.31au，而esa计算的半长轴与2组最低相差～0.74au。对比偏心率e与轨道倾角i，可以发现同样的规律。这说明，对2022 eb5而言，jpl计算的解比esa计算的解更合理。
[0216]
此外，考虑大气层的三组c类小天体轨道根数比不考虑大气层(“无大气”组)的轨道根数更接近jpl的轨道解。使用本方法计算的轨道根数与jpl解的相对误差如表11所示。
[0217]
表11 2022 eb5日心轨道根数相对误差(％)
[0218][0219]
其中，各相对误差的单位为％。
[0220]
从表11能够更直观地看到，考虑大气层的三组轨道根数在半长轴、偏心率、轨道倾角三个方面比无大气组更加准确。而大气层对近地点俯角与升交点赤经的影响非常小，二者的偏差小于0.1％。
[0221]
对小天体2022 eb5的计算表明，该方法在判断小天体类型、计算小天体解体强度、反照率、直径等主要物理参数方面具有较大优势。考虑大气层后的轨道根数比不考虑大气层的轨道根数更加准确。
[0222]
3.总结
[0223]
通过2022 eb5，对本技术的技术效果进行总结。
[0224]
(1)本方法对小天体类型进行了较为准确的推断，仅使用轨道观测数据无法对小天体类型进行判断；
[0225]
(2)本方法对小天体的解体强度、反照率、直径进行了更好的约束，与不考虑大气层相比，解体强度相差40％-60％；将反照率范围约束在0.01-0.035，该范围远小于0.147的平均反照率；直径范围5-6m，远大于使用平均反照率估算的1.69m。
[0226]
(3)使用本方法计算得到2022 eb5的轨道半长轴、偏心率、轨道倾角，比不考虑大气层得到的轨道根数更接近jpl解。大气层对近地点俯角、升交点赤经的影响很小。但总体
而言，考虑大气层对轨道计算非常重要，不可忽略。
[0227]
本发明的创新点在于：
[0228]
1、充分考虑地球大气层对小天体各物理量的影响，对小天体物理特性进行约束。
[0229]
2、将遗传算法与煎饼模型相结合对小天体的未知变量进行求解。
[0230]
3、根据小天体大气进入过程计算小天体的日心轨道根数，对小天体观测轨道的可靠性与小天体类型进行验证。
[0231]
本发明还可提供的一种计算机设备，包括：至少一个处理器、存储器、至少一个网络接口和用户接口。该设备中的各个组件通过总线系统耦合在一起。可理解，总线系统用于实现这些组件之间的连接通信。总线系统除包括数据总线之外，还包括电源总线、控制总线和状态信号总线。
[0232]
其中，用户接口可以包括显示器、键盘或者点击设备(例如，鼠标，轨迹球(track ball)、触感板或者触摸屏等。
[0233]
可以理解，本技术公开实施例中的存储器可以是易失性存储器或非易失性存储器，或可包括易失性和非易失性存储器两者。其中，非易失性存储器可以是只读存储器(read-only memory，rom)、可编程只读存储器(programmable rom，prom)、可擦除可编程只读存储器(erasable prom，eprom)、电可擦除可编程只读存储器(electrically eprom，eeprom)或闪存。易失性存储器可以是随机存取存储器(random access memory，ram)，其用作外部高速缓存。通过示例性但不是限制性说明，许多形式的ram可用，例如静态随机存取存储器(static ram，sram)、动态随机存取存储器(dynamic ram，dram)、同步动态随机存取存储器(synchronous dram，sdram)、双倍数据速率同步动态随机存取存储器(double data rate sdram，ddrsdram)、增强型同步动态随机存取存储器(enhanced sdram，esdram)、同步连接动态随机存取存储器(synchlink dram，sldram)和直接内存总线随机存取存储器(direct rambus ram，drram)。本文描述的存储器旨在包括但不限于这些和任意其它适合类型的存储器。
[0234]
在一些实施方式中，存储器存储了如下的元素，可执行模块或者数据结构，或者他们的子集，或者他们的扩展集：操作系统和应用程序。
[0235]
其中，操作系统，包含各种系统程序，例如框架层、核心库层、驱动层等，用于实现各种基础业务以及处理基于硬件的任务。应用程序，包含各种应用程序，例如媒体播放器(media player)、浏览器(browser)等，用于实现各种应用业务。实现本公开实施例方法的程序可以包含在应用程序中。
[0236]
在本上述的实施例中，还可通过调用存储器存储的程序或指令，具体的，可以是应用程序中存储的程序或指令，处理器用于：
[0237]
执行上述方法的步骤。
[0238]
上述方法可以应用于处理器中，或者由处理器实现。处理器可能是一种集成电路芯片，具有信号的处理能力。在实现过程中，上述方法的各步骤可以通过处理器中的硬件的集成逻辑电路或者软件形式的指令完成。上述的处理器可以是通用处理器、数字信号处理器(digital signal processor，dsp)、专用集成电路(application specific integrated circuit，asic)、现场可编程门阵列(field programmable gate array，fpga)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件。可以实现或者执行上述公开的
各方法、步骤及逻辑框图。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。结合上述公开的方法的步骤可以直接体现为硬件译码处理器执行完成，或者用译码处理器中的硬件及软件模块组合执行完成。软件模块可以位于随机存储器，闪存、只读存储器，可编程只读存储器或者电可擦写可编程存储器、寄存器等本领域成熟的存储介质中。该存储介质位于存储器，处理器读取存储器中的信息，结合其硬件完成上述方法的步骤。
[0239]
可以理解的是，本发明描述的这些实施例可以用硬件、软件、固件、中间件、微码或其组合来实现。对于硬件实现，处理单元可以实现在一个或多个专用集成电路(application specific integrated circuits，asic)、数字信号处理器(digital signal processing，dsp)、数字信号处理设备(dsp device，dspd)、可编程逻辑设备(programmable logic device，pld)、现场可编程门阵列(field-programmable gate array，fpga)、通用处理器、控制器、微控制器、微处理器、用于执行本技术所述功能的其它电子单元或其组合中。
[0240]
对于软件实现，可通过执行本发明的功能模块(例如过程、函数等)来实现本发明技术。软件代码可存储在存储器中并通过处理器执行。存储器可以在处理器中或在处理器外部实现。
[0241]
本发明还可提供一种非易失性存储介质，用于存储计算机程序。当该计算机程序被处理器执行时可以实现上述方法实施例中的各个步骤。
[0242]
最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。