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

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

Info

Publication number
CN115292657A
CN115292657A CN202210776924.5A CN202210776924A CN115292657A CN 115292657 A CN115292657 A CN 115292657A CN 202210776924 A CN202210776924 A CN 202210776924A CN 115292657 A CN115292657 A CN 115292657A
Authority
CN
China
Prior art keywords
celestial body
small celestial
small
radius
atmosphere
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202210776924.5A
Other languages
English (en)
Other versions
CN115292657B (zh
Inventor
耿淑娟
周炳红
李明涛
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National Space Science Center of CAS
Original Assignee
National Space Science Center of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by National Space Science Center of CAS filed Critical National Space Science Center of CAS
Priority to CN202210776924.5A priority Critical patent/CN115292657B/zh
Publication of CN115292657A publication Critical patent/CN115292657A/zh
Application granted granted Critical
Publication of CN115292657B publication Critical patent/CN115292657B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/12—Simultaneous equations, e.g. systems of linear equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/12—Computing arrangements based on biological models using genetic models
    • G06N3/126—Evolutionary algorithms, e.g. genetic algorithms or genetic programming
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Biophysics (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Operations Research (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Evolutionary Computation (AREA)
  • Artificial Intelligence (AREA)
  • Genetics & Genomics (AREA)
  • Physiology (AREA)
  • Biomedical Technology (AREA)
  • Computational Linguistics (AREA)
  • Molecular Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Computing Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明提供了一种求解撞击地球小天体的物理特性的方法，所述方法为：根据小天体撞击数据，结合煎饼模型与遗传算法，考虑大气层对小天体各物理量以及轨道的影响，计算小天体的物理特性，利用已知小天体类型和物理特性对照表对计算结果进行筛选；其中，小天体的物理特性包括：质量、半径、密度和解体强度。本方法在判断小天体类型、计算小天体解体强度、反照率、直径等主要物理参数方面具有较大优势。

Description

一种求解撞击地球小天体的物理特性的方法
技术领域
本发明属于小行星防御与利用领域，具体涉及一种求解撞击地球小天体的物理特性的方法。
背景技术
小天体撞击地球的频率与尺寸呈负相关。小尺寸的小天体撞击地球频率高，但其较小的体积与较低的反照率使他们更加难以观测。同时，小尺寸小天体由于其撞击地球带来的危害效应相对轻微，往往并不是天文台或望远镜的观测目标。然而，1908年通古斯大爆炸与2013年车里雅宾斯克撞击事件表明，即使尺寸在几十米甚至几米的小天体，在撞击地球大气层过程中也能够产生剧烈的冲击波，从而可能导致严重的人员伤亡与经济损失。因此，中小尺寸(直径1至140米)的小天体正受到越来越广泛的关注与重视。
2022年我国航天局明确表示中国将建设近地小天体防御系统。为更好地了解小天体的撞击威胁，评估小天体撞击效应，有必要了解小天体(尤其是撞击频率高的中小尺寸小天体)的物理性质。小天体的物理特性主要指其结构、材质、密度、尺寸等性质。在研究中，结构通常由小天体的解体强度来反映。
获取撞击地球的小天体物理特性的方式主要有三种。
第一种为对小行星/彗星进行采样或收集陨石，直接获取小天体信息。目前人类已经实现了对“贝努”小行星和“龙宫”小行星的采样返回，获得了珍贵的材料样本。此外，若小天体的尺寸在米级以上，则很有可能在地球留下陨石。陨石或小天体岩石样本中包含小天体的结构信息，如小天体的材质、密度等，也能够提供太阳系早期演化的线索，是珍贵的研究资料。然而，由于小天体种类繁杂数量庞大，不同小天体之间结构差异也可能很大，而通过采样或陨石收集得到材料的小天体数量远不及日常撞击地球的小天体数量。因此，该方法虽然直接，但具有数量和统计上的局限性。
第二种方式为光学观测，主要指对小天体的天文学观测。通过这种方式，研究者能够获取或部分获取小天体的绝对星等、反照率、自转速度、近红外光谱等数据。这些数据可进一步用于推断小天体的尺寸、表面材质、内聚强度。然而，这类方法通常用于直径千米级的小天体物理特性的估算，很少应用于数量庞大且撞击地球频率更高的小尺寸小天体。
这首先是因为小尺寸小天体的观测数据少。截至2022年5月23日，在撞击地球前被人类观测到的小天体只有5个，他们分别为2008TC3，2014AA，2018LA，2019MO和2022 EB5。这意味着在可预见的未来，人类对可能撞击地球的小天体的观测仍有待完善。第二个原因是利用该方法估算小尺寸小天体性质，其相对误差比千米及小天体的误差要大。相较于小天体的材质与反照率，绝对星等更容易获得，此时通常会假设小天体的反照率为平均反照率0.147，借助绝对星等、反照率与尺寸的经验关系估算小天体尺寸。该方法用于估算千米级小天体的尺寸，相对误差较小，但估算小尺寸小天体，则会导致较大误差。
另外一种方法是通过观察小天体撞击地球过程中与大气层的相互作用确定小天体的物理特性，如小天体在大气进入过程中的解体现象、减速与烧蚀过程、发光现象等。全球很多地区已经建立了火流星观测网络。这些观测网通常为相机网络，能够在各自的视场观测并记录撞击地球的小天体在大气内发光飞行阶段的轨迹。对光辐射过程进行分析，能够还原小天体飞行的能量沉积过程。此外，NASA的近地小天体研究中心(CNEOS)自1988年起，通过卫星传感器记录并发布小天体在峰值亮度时刻的位置、速度及估算的总能量，且近几年已经实现了全球覆盖，已经成为确定小天体物理性质非常重要的数据来源。
然而，尽管小天体发光飞行阶段的轨迹已经被记录下来，但它的质量、尺寸及解体强度等信息仍然无法准确获知。前人的研究中，通常将小天体峰值亮度时刻的前缘动压作为小天体的解体强度，并假设小天体进入大气至达到峰值亮度期间，大气层对小天体的速度没有影响，再根据估算的总动能进一步计算小天体的质量。通过对小天体密度进行假设，给出小天体的名义直径。通过此种方法计算的小天体解体强度往往比小天体的真实解体强度大很多，小天体材质未知的情况下计算的名义直径也可能与真实的尺寸差别较大。
综合以上分析，能够发现，获取中小尺寸小天体的物理特性具有非常重要的意义；而对于撞击地球频率高的中小尺寸小天体，尚缺乏一种能够大规模地准确地推算小天体物理特性的方法。
发明内容
本发明的目的在于克服现有计算方法结果与真实值偏差大的缺陷。
为了实现上述目的，本发明提出了一种求解撞击地球小天体的物理特性的方法，所述方法为：根据小天体进入大气后的物理数据，结合煎饼模型与遗传算法，考虑大气层对小天体各物理量以及轨道的影响，计算小天体的物理特性，利用已知小天体类型和物理特性对照表对计算结果进行筛选，从而得到待求解小天体的物理特性；其中，小天体的物理特性包括：质量、半径、密度和解体强度。
作为上述方法的一种改进，所述方法具体包括：
步骤1：获取小天体进入大气信息；获取的信息包括峰值亮度时刻、峰值亮度时刻位置、峰值亮度时刻速度和撞击总能量；
步骤2：利用小天体峰值亮度时刻的物理量，通过四阶龙格-库塔方法对小天体飞行方程、烧蚀方程和半径变化方程构成的微分方程组进行逆向数值积分，得到小天体进入大气层边界处的速度矢量
Figure BDA0003727821180000031
位置矢量
Figure BDA0003727821180000032
半径R0、半径变化率dR0和质量m0；通过遗传算法计算小天体的解体强度S0、小天体在峰值亮度时刻的半径Rp、半径变化率dRp和质量mp；
步骤3：根据小天体类型与解体强度的对应关系，筛选出可能正确解。
作为上述方法的一种改进，所述飞行方程为：
Figure BDA0003727821180000033
其中，
Figure BDA0003727821180000034
为小天体相对地球的矢径；r为
Figure BDA0003727821180000035
的模，是小天体到地球的距离；m为小天体的质量；μ为地球引力常数，取3.986×1014m3/s2；t为时间；
Figure BDA0003727821180000036
表示小天体受到的大气阻力，满足下式：
Figure BDA0003727821180000037
其中，Cd为阻力系数，取值为1；Vr为小天体相对地球大气层的速度大小，亦为小天体在ECEF坐标系下的速度大小；S为迎风截面积；
Figure BDA0003727821180000038
为单位矢量，表示
Figure BDA0003727821180000039
的方向；ρa为大气密度。
作为上述方法的一种改进，所述烧蚀方程为：
Figure BDA00037278211800000310
其中，m为小天体的质量；t为时间；Ch为烧蚀系数；Qh为烧蚀热；ρa为大气密度；Vr为小天体相对地球大气层的速度大小，亦为小天体在ECEF坐标系下的速度大小；S为迎风截面积
作为上述方法的一种改进，所述半径变化方程为：
Figure BDA00037278211800000311
其中，R为小天体的半径；t为时间；ρa为大气密度；Vr为小天体相对地球大气层的速度大小，亦为小天体在ECEF坐标系下的速度大小；ρm为小天体的密度。
作为上述方法的一种改进，所述逆向数值积分计算公式为：
Figure BDA0003727821180000041
其中，Xb表示逆向数值积分结果；hp为小天体峰值亮度时刻的高度；hb为大气层边界高度hb＝100km；h为小天体高度；
Figure BDA0003727821180000042
t为时间；X代指小天体被积分的物理量，包括矢径
Figure BDA0003727821180000043
速度矢量
Figure BDA0003727821180000044
质量m，半径R，半径变化率
Figure BDA0003727821180000045
Figure BDA0003727821180000046
逆向数值积分过程由四阶龙格-库塔方法实现，最终得到小天体在大气层边界处ECI坐标系下的速度矢量
Figure BDA0003727821180000047
位置矢量
Figure BDA0003727821180000048
半径R0、半径变化率dR0和质量m0。
作为上述方法的一种改进，所述逆向数值积分时，未知变量包括：小天体峰值亮度时刻的质量mp、半径Rp及半径变化率dRp，小天体的密度ρm、解体强度S0、烧蚀热Qh，模型参数N；
其中，N取为4，当N取4无解时，N取2；
根据彗星、C类小行星、S类小行星、M类小行星四类小天体的类型与密度ρm、烧蚀热Qh的对应关系，取四组密度ρm、烧蚀热Qh分别求解；
彗星对应密度ρm为1000kg/m3，对应烧蚀热Qh为2.5×106J/kg；
C类小行星对应密度ρm为2000kg/m3，对应烧蚀热Qh为5×106J/kg；
S类小行星对应密度ρm为3500kg/m3，对应烧蚀热Qh为8×106J/kg；
M类小行星对应密度ρm为7900kg/m3，对应烧蚀热Qh为8×106J/kg；
另外四个未知量S0、mp、Rp及dRp，作为优化参数，利用遗传算法进行求解：
遗传算法的四个约束条件分别为：
(1)计算得到的大气层边界处的半径变化率为0，即dR0＝0；
(2)计算得到的大气层边界处的半径与质量满足球体假设，即R(m)-R0＝0，其中
Figure BDA0003727821180000051
(3)计算得到的大气层边界处的动能与获取的小天体进入大气信息中的撞击总能量相等，即E(V0,m0)-E0＝0，其中
Figure BDA0003727821180000052
(4)计算得到的大气层边界处半径应与峰值亮度时刻半径Rp满足煎饼模型的比例关系，即NR0-Rp＝0；
其中，E0为获取到的总撞击能量；m0为大气层边界处小天体质量；R0为大气层边界处小天体半径；ρap为小天体峰值亮度时刻大气密度；Vp为小天体峰值亮度时刻速度；
遗传算法中四个优化参数S0、mp、Rp及dRp取值范围为：
S0最大值为Pp，最小值为Pp/30，
Figure BDA0003727821180000053
mp最大值为2me，最小值为0.01me，
Figure BDA0003727821180000054
Rp最大值为NR0max，最小值为R0min，
Figure BDA0003727821180000055
Figure BDA0003727821180000056
dRp最大值为50N，最小值为0.1；
遗传算法中目标函数f由约束(1)至(4)加权组合而成，公式如下：
Figure BDA0003727821180000057
其中，a1、a2、a3、a4为四个约束的权重，初始值为10、100、0.01、50；
求得的解是否符合要求的判断标准为以下四项全部满足：
|R0-R(m0)|/R0＜1％
|NR0-Rp|/Rp＜1％
|m0-m(E0,V0)|/m0＜0.1％
|dR0|＜0.001
通过遗传算法最终计算得到小天体的解体强度S0、小天体在峰值亮度时刻的半径Rp、半径变化率dRp和质量mp。
作为上述方法的一种改进，根据计算出的小天体的物理特性和状态向量可以计算小天体日心轨道的轨道根数，包括：径向速度、比角动量、轨道倾角、节点线向量、升交点赤经、偏心率、近地点俯角、真近点角和轨道半长轴。
作为上述方法的一种改进，求出小天体轨道后，可根据小天体相对木星的蒂塞朗常数TJ，远日点距离Q、近日点距离q及木星半长轴aJ的距离之间的关系，对小天体的类型进行验证。
作为上述方法的一种改进，如果已知小天体的绝对星等，通过式(20)根据求得的直径，反推小天体的反照率，进一步根据不同类别小天体的反照率对小天体材质进行验证或约束：
Figure BDA0003727821180000061
其中，H为小天体的绝对星等，D为小天体的直径，pv为反照率。
与现有技术相比，本发明的优势在于：
1、本方法结合煎饼模型与遗传算法，在判断小天体类型、计算小天体解体强度、反照率、直径等主要物理参数方面具有较大优势。
2、本方法在飞行方程中考虑了大气层对小天体飞行的影响，计算的轨道根数比不考虑大气层的轨道根数更加准确。
附图说明
图1所示为求解撞击地球小天体的物理特性的方法流程图；
图2所示为煎饼模型示意图；
图3所示为2022 EB5反照率与直径关系图。
具体实施方式
本发明旨在利用小天体与地球大气层间的相互作用，对小天体的物理特性进行求解或约束，同时由此计算小天体的日心轨道根数。由于NASA近地小天体研究中心(CNEOS)能够记录大部分米级及以上小天体在峰值亮度时刻的速度、位置并根据经验关系估算撞击总能量，本发明以覆盖率广的CNEOS火球数据为基础，以使本方法具有更高的普遍性。本发明利用半解析解体模型，充分考虑大气层对小天体速度的影响，对小天体的物理特性进行推断。对于未知物理量，本发明将其视为优化参数，采用遗传算法对未知数进行求解。
下面结合附图对本发明的技术方案进行详细的说明。
本发明利用CNEOS提供的小天体撞击地球过程中峰值亮度时刻的速度、位置及撞击能量，推算小天体的物理特性。整体计算流程如图1所示。该过程中使用的数据与方法如下：
1.获取小天体进入大气信息
CNEOS公布的撞击地球的小天体信息列表如表1所示，其中Vx，Vy，Vz分别表示小天体的速度矢量在地心固连(ECEF)坐标系下x、y、z方向下的分量，Vr表示ECEF系下速度矢量的模，即速度的大小。
表1小天体峰值亮度时刻的速度、位置及撞击能量
Figure BDA0003727821180000071
表1中的数据将进一步转换为地心惯性(ECI)坐标系，得到ECI系下小天体峰值亮度时刻的状态向量，即位置矢量rp，速度矢量Vp，以及撞击总能量，也为小天体在地球大气层边界处的动能。
2.小天体进入大气过程飞行方程
小天体进入大气的速度理论上在11.2km/s至72km/s之间，因此小天体与大气层之间的撞击为超高速撞击。在大气进入阶段，小天体的速度将在大气阻力的作用下降低，即减速；同时与大气之间剧烈的相互作用使小天体表面温度升高，导致小天体的质量损失，又称为烧蚀。
减速由小天体的飞行方程控制，如式(1)所示。
Figure BDA0003727821180000072
其中，
Figure BDA0003727821180000073
为小天体相对地球的矢径，r为
Figure BDA0003727821180000074
的模，是小天体到地球的距离，m为小天体的质量；μ为地球引力常数，取3.986×1014m3/s2。
Figure BDA0003727821180000075
表示小天体受到的大气阻力，由式(2)决定：
Figure BDA0003727821180000081
上式中，Cd为阻力系数，在本计算中按球体假设取为1；Vr为小天体相对地球大气层的速度大小，亦为小天体在ECEF坐标系下的速度大小。S为迎风截面积，
Figure BDA0003727821180000082
为单位矢量，表示
Figure BDA0003727821180000083
的方向。ρa为大气密度，根据Wheeler等人对1976美国标准大气模型的拟合公式计算：
ρa＝-140.2e-0.000187h+141.4e-0.000186h (3)
其中，h为海拔高度，单位为m。
小天体的质量损失由烧蚀方程控制，如式(4)所示：
Figure BDA0003727821180000084
其中，Ch为烧蚀系数，取0.1；Qh为烧蚀热，其取值与小天体类型有关，后续将在介绍小天体分组时统一赋值。
3.煎饼模型
随着小天体进入地球大气层，小天体周围的大气密度也随高度的降低而近似指数增加，小天体的前缘所受的动压P＝ρaVr 2也随之增加。通常认为，当P超过小天体解体强度S0时，小天体会发生解体。
解体是小天体进入大气过程中非常重要的物理现象，也是剧烈的空爆现象的直接诱因。本发明使用煎饼模型对小天体的解体现象进行模拟，其基本原理示意如图2所示。
煎饼模型认为，小天体解体后，大大小小的碎片构成一个碎片云团。在解体初期，该云团仍然受到同一个激波的作用，因此在这一阶段依然可以将其视为一个整体。随后，各碎片逐渐分离，彼此的间距变大，整体表现为云团的迎风截面半径增加，该过程中，半径R的变化由式(5)求得：
Figure BDA0003727821180000085
其中，ρm为小天体的密度。
当云团的半径增加至初始半径的N倍时，各碎片发展出各自的激波，不能再作为一个整体进行计算。N为经验参数，是小天体半径变化程度的度量，在模型提出时取2，但在应用中发现，N通常取4-8才能够实现对解体过程的较好拟合。此时，小天体的减速与烧蚀作用迅速加剧，能量在短时间内大量释放，导致光辐射升高达到峰值，同时产生剧烈的冲击波，该现象即为空爆。从描述中可知，空爆与峰值亮度同步发生。
4.未知变量的求解
利用小天体峰值亮度时刻的物理量，通过对小天体飞行方程、烧蚀方程和半径变化方程组成的微分方程组进行逆向数值积分，可得到小天体进入大气层边界处的物理量，包括：，速度矢量
Figure BDA0003727821180000091
位置矢量
Figure BDA0003727821180000092
半径R0、半径变化率dR0和质量m0。
逆向数值积分过程从小天体峰值亮度时刻的高度hp向上积分至大气层边界处(大气层边界高度为hb＝100km)：
Figure BDA0003727821180000093
其中
Figure BDA0003727821180000094
X代指小天体被积分的物理量，包括矢径
Figure BDA0003727821180000095
速度矢量
Figure BDA0003727821180000096
(即
Figure BDA0003727821180000097
)，质量m，半径R，半径变化率
Figure BDA0003727821180000098
Figure BDA0003727821180000099
各物理量的微分方程表达如式(1)至式(5)所示。
逆向积分过程由四阶龙格-库塔方法实现。该方法为求解非线性常微分方程的常用方法，利用微分方程与初始条件，通过逐步迭代实现高精度的求解。
对飞行过程进行求解，会发现，仍有一些小天体的物理量未知：小天体峰值亮度时刻的质量mp、半径Rp及半径变化率dRp，小天体的密度ρm、解体强度S0、烧蚀热Qh，模型参数N。
N的取值范围为2-8，对于小尺寸的小天体而言，N的取值不宜过大，因此本方法中N取为4，当N取4无解时，N取2。密度ρm与烧蚀热Qh均与小天体类型有关，对于未捡到陨石的小天体而言，其类型往往未知，但它与小天体的解体强度存在定性的关联：通常情况下，彗星、C类小行星、S类小星星、M类小行星的密度递增，解体强度也递增。因此，本方法中，首先将计算按材质分为四组，每组的小天体密度均取该类型下典型的小天体密度值。四组均求得解后，再按照求得的解体强度对小天体的类型进行约束，并求出最终的小天体类型。按小天体类型的分组情况如表2所示：
表2小天体类型、密度、烧蚀热
Figure BDA0003727821180000101
小天体类型与材质存在对应关系，彗星主要材质为水冰，C类、S类、M类主要对应碳质、石质、铁质。
另外四个未知量S0、mp、Rp及dRp，将作为优化参数，利用遗传算法进行求解。
遗传算法是一种通过模拟自然选择以求问题最优解的优化方法，其工作流程为：
(1)在优化变量的取值范围内生成初始种群，优化变量包括小天体的解体强度、小天体峰值亮度时刻的质量、半径、半径变化率；
(2)计算目标函数值：利用四阶龙格库塔方法，对小天体飞行方程、半径变化方程、烧蚀方程进行逆向数值积分，计算小天体在大气层边界处的速度位置矢量、质量、半径、半径变化率。并根据边界处的物理量计算目标函数值。
(3)判断是否满足终止条件，若是，则输出优化结果；若否，则对初代种群进行遗传、交叉、变异，产生新一代个体，重复(2)和(3)。
在该问题中，存在四个独立的约束条件，未知量与条件数量相等，能够保证存在唯一解。四个约束条件分别为：
(1)计算得到的大气层边界处的半径变化率为0，即dR0＝0；
(2)计算得到的大气层边界处的半径与质量满足球体假设，即R(m)-R0＝0，其中
Figure BDA0003727821180000102
(3)计算得到的大气层边界处的动能与CNEOS估算的撞击总能量相等，即E(V0,m0)-E0＝0，其中
Figure BDA0003727821180000103
(4)计算得到的大气层边界处半径应与峰值亮度时刻半径Rp满足煎饼模型的比例关系，即NR0-Rp＝0。
上述描述中，除能量E0外，下标0表示通过对飞行方程逆向积分求得的大气层边界处的物理量。E0为CNEOS估算得到的总撞击能量。下标p表示峰值亮度时刻的物理量。
该求解过程由遗传算法实现。遗传算法是一种模拟自然选择以求解最小值问题的最优化算法。该算法首先从优化参数的取值范围中选取一定数目的个体，利用适应度函数或目标函数的取值作为个体的性状，并对个体的性状是否符合进化目标进行评估。适应性高的个体将有机会将基因传递给下一代，而性状不够优秀的个体将被淘汰。缺失的个体位将继续在参数的取值范围内随机选取。如此产生一代代个体，直至达到最大代数或达到进化要求。
本方法中，四个优化参数S0、mp、Rp及dRp取值范围的选取在充分考虑撞击事件个别性的基础上设定，其取值如表3所示。
表3优化变量取值范围
Figure BDA0003727821180000111
目标函数f由约束(1)至(4)加权组合而成，如式(6)所示：
Figure BDA0003727821180000112
其中，a1、a2、a3、a4为四个约束的权重，初始值为10、100、0.01、50以使四个约束值维持在同一量级，也能够根据撞击事件的实际情况进行调解。同时，求得的解是否符合要求的判断标准为：|R0-R(m0)|/R0＜1％，|NR0-Rp|/Rp＜1％，|m0-m(E0,V0)|/m0＜0.1％，|dR0|＜0.001。
通过(4)中方法求出未知量后，小天体大气层边界处的状态向量(位置矢量
Figure BDA0003727821180000113
和速度矢量
Figure BDA0003727821180000114
)也可求得。利用小天体的解体强度与材质之间的关系，可对小天体的材质进行判断。
5.轨道计算
经过对未知变量的求解，本发明中的核心计算已经完成，但为了使物理特性的推断更为全面，本方法将继续计算小天体日心轨道的轨道根数。对于有轨道观测数据的小天体，计算轨道可以验证本发明中方法的可靠性；对于缺乏轨道观测数据的小天体，计算轨道能够从轨道的统计规律角度对小天体类型(仅判断为彗星/小星星)进行验证。
获得小天体大气层边界处ECI坐标系下的状态向量
Figure BDA0003727821180000121
后，需要进一步对轨道方程进行逆向积分，求得小天体在地球影响球边界处的ECI系下的状态向量。进一步可将状态向量转换为小天体在日心黄道坐标系下的轨道根数。具体步骤如下：
(1)计算小天体在地球影响球边界处ECI坐标系下的状态向量：
该过程需要对轨道方程进行逆向积分，直至小天体与地球的距离达到地球影响球边界处。该阶段的轨道方程可以不考虑地球大气的影响，因为没有大气阻力项，如式(7)所示：
Figure BDA0003727821180000122
地球影响球半径
Figure BDA00037278211800001212
(2)将ECI坐标系下的状态向量转变为地心黄道坐标系下的状态向量：
Figure BDA0003727821180000123
其中，
Figure BDA0003727821180000124
与
Figure BDA0003727821180000125
为小天体在地球影响球边界处、地心黄道系下的状态向量，
Figure BDA0003727821180000126
Figure BDA0003727821180000127
为小天体在地球影响球边界处、ECI坐标系下的状态向量；ε为黄赤交角，取23.43929°。
(3)计算小天体在日心黄道坐标系下的状态向量：
日心黄道坐标系与地心黄道坐标系均为惯性系，且坐标平面一致，因此只需简单的计算即可求得：
Figure BDA0003727821180000128
其中，
Figure BDA0003727821180000129
与
Figure BDA00037278211800001210
为小天体日心黄道坐标系下的状态向量，
Figure BDA00037278211800001211
为地球在日心黄道坐标系下的状态向量。
(4)计算小天体日心黄道坐标系下的轨道根数
为方便表示，步骤(4)中日心黄道系下的
Figure BDA0003727821180000131
以
Figure BDA0003727821180000132
和
Figure BDA0003727821180000133
表示。
计算径向速度Vr：
Figure BDA0003727821180000134
计算比角动量
Figure BDA0003727821180000135
Figure BDA0003727821180000136
计算轨道倾角i：
Figure BDA0003727821180000137
计算节点线向量
Figure BDA0003727821180000138
Figure BDA0003727821180000139
计算升交点赤经Ω：
Figure BDA00037278211800001310
计算偏心率e：
Figure BDA00037278211800001311
计算近地点俯角ω：
Figure BDA00037278211800001312
计算真近点角θ：
Figure BDA00037278211800001313
Figure BDA0003727821180000141
计算轨道半长轴a：
Figure BDA0003727821180000142
6.小天体类型验证
当求出小天体轨道后，可根据小天体相对木星的蒂塞朗常数TJ，远日点距离Q、近日点距离q及木星半长轴aJ的距离之间的关系，对小天体的类型(仅能判断彗星/小行星)进行验证。类型判断标准如表4所示。
表4彗星、小行星的轨道统计规律
Figure BDA0003727821180000143
蒂塞朗常数计算方式如式(19)所示：
Figure BDA0003727821180000144
需要注意的是，该准则为统计规律，且仅能对小天体属于彗星或小行星做出大致判断，因此其分类无法作为判断小天体类型的绝对标准。
7.绝对星等
小天体的绝对星等H、直径D与反照率pv之间存在经验关系，如式(20)所示：
Figure BDA0003727821180000145
在利用遗传算法及煎饼模型对小天体大气进入过程进行求解后，最终得到四组小天体物理参数，通过对比解体强度与材质是否符合经验，可以对不合理的几组进行排除。然而，由于C类小行星与S类小行星的强度界限较为模糊，因此最终将剩余一或两组数据。
此时，若小天体的绝对星等已知，那么可以通过式(20)可以根据求得的直径，反推小天体的反照率，进一步根据不同类别小天体的反照率对小天体材质进行验证或约束。C类小行星的反照率范围通常为0.03-0.11，S类小行星的反照率范围约为0.13-0.29，M类小行星的反照率约为0.1-0.3.
2022年3月11日，匈牙利天文学家观测到一颗编号为2022 EB5的小天体将撞击地球。在被发现的两小时后，该小天体进入地球大气层，并在挪威海上空发生空爆，成为第五颗在撞击地球前被观测到的小天体。
为展示本申请的技术效果，同时验证本申请的正确性，现采用本申请中的方法对2022 EB5的物理特性进行分析，并对其日心轨道进行计算。
1.观测信息
1.1轨道观测信息
本方法在计算过程中，并不会使用到小天体的轨道观测信息，但已知轨道观测信息能够从侧面对本方法的正确性与合理性进行验证与说明，有利于技术效果的展示。
2022 EB5撞击地球前的轨道被多处天文台观测到，NASA的JPL实验室与ESA根据不同天文台的观测数据，分别计算了2022 EB5的日心轨道根数，如表5和表6所示。
表5 2022 EB5轨道根数(JPL解)
Figure BDA0003727821180000151
表6 2022 EB5轨道根数(ESA解)
Figure BDA0003727821180000152
其中，q为近日点距离，Q为远日点距离，M为平近点角。
可见，JPL与ESA的轨道解存在较大差异。在呈现本方法计算结果时，将利用本方法计算得到的轨道根数与以上两组轨道解进行对比，可以对两组解的可靠性进行判断。
1.2CNEOS峰值亮度观测信息
CNEOS的观测数据为小天体进入大气过程中峰值亮度时刻的信息，为本方法计算中使用的主要数据。
2022 EB5大气进入过程中，其峰值亮度时刻的速度、位置、时间及其撞击总能量如表7所示。
表7 2022 EB5的峰值亮度观测数据
Figure BDA0003727821180000161
表7中各参数含义与表1相同。
2.计算结果
2.1 2022 EB5物理特性分析
为更好地展示本方法的技术效果，本节首先将讨论在不使用本方法计算的情况下，利用已知信息能够对小天体的物理特性进行何种约束。
根据轨道根数，可以计算2022 EB5相对于木星的蒂瑟朗常数，T＝2.89。远日点为4.77AU，大于木星半长轴。根据蒂瑟朗常数判据，2022 EB5很可能为彗星小天体，但根据远日点距离、木星半长轴之间的关系判断，2022 EB5大概率为小行星。因此，根据轨道难以判断其类型。
对2022 EB5的观测给出了其绝对星等。根据式(20)可以得到2022 EB5直径与反照率之间的关系，如图3所示。然而，由于不确定2022 EB5的小天体类型，其反照率信息未知。此时，通常假设其反照率为平均反照率0.147，在此假设下，计算得到的2022 EB5直径为1.69m。
根据CNEOS对小天体峰值亮度时刻的观测数据，按照以往的估算方法，假设小天体进入大气时刻的速度与峰值亮度时刻的速度相等，即为17.24km/s，则根据初始动能，得到小天体质量约112640kg。按该估算方法，小天体解体强度约为峰值亮度时刻的总压，即为3.5MPa。
以上为不使用本方法，按以往常规处理方法估算得到的小天体物理特性。
使用本方法，考虑大气层对小天体速度、质量的影响，得到四组小天体大气层边界处速度、质量及小天体解体强度的解，如表8所示。其中，第四列为各自类型的小天体解体强度的经典范围。
表8 2022 EB5物理特性及大气层边界处物理参数表
Figure BDA0003727821180000171
从表8中第三列与第四列可以看到，根据考虑大气层得到的解体强度，2022 EB5不可能为彗星与M类小行星，其最有可能为C类小行星，然而，仅靠强度不能完全排除其为S类小行星的可能性。此外，考虑大气层的影响计算得到的解体强度(1.4-2.1MPa)，明显低于不考虑大气层(即不使用本方法)得到的解体强度3.5MPa，二者相对差高达40％-60％。
当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类小行星，且其反照率非常低。
此外，考虑到已知密度的C类小行星——“龙宫”与“贝奴”小行星，其密度均为～1200kg/m3，均小于表8中设定的2000kg/m3，因此本节按照以上两颗小行星的密度增加两组计算结果，如表9所示。
表9 2022 EB5物理特性及大气层边界处物理参数表
Figure BDA0003727821180000172
Figure BDA0003727821180000181
其中，C类1的密度设定与“龙宫”小天体一致，C类2的密度设定与“贝奴”小天体一致。
2.2 2022 EB5轨道计算结果
将大气层边界处的状态向量转换为日心轨道，其结果如表10所示。
表10使用本方法计算得到的2022 EB5日心轨道根数
Figure BDA0003727821180000182
表10中，组别中“JPL”与“ESA”分别表示JPL与ESA根据轨道观测数据计算的2022EB5的日心轨道参数(表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计算的解更合理。
此外，考虑大气层的三组C类小天体轨道根数比不考虑大气层(“无大气”组)的轨道根数更接近JPL的轨道解。使用本方法计算的轨道根数与JPL解的相对误差如表11所示。
表11 2022 EB5日心轨道根数相对误差(％)
Figure BDA0003727821180000183
其中，各相对误差的单位为％。
从表11能够更直观地看到，考虑大气层的三组轨道根数在半长轴、偏心率、轨道倾角三个方面比无大气组更加准确。而大气层对近地点俯角与升交点赤经的影响非常小，二者的偏差小于0.1％。
对小天体2022 EB5的计算表明，该方法在判断小天体类型、计算小天体解体强度、反照率、直径等主要物理参数方面具有较大优势。考虑大气层后的轨道根数比不考虑大气层的轨道根数更加准确。
3.总结
通过2022 EB5，对本申请的技术效果进行总结。
(1)本方法对小天体类型进行了较为准确的推断，仅使用轨道观测数据无法对小天体类型进行判断；
(2)本方法对小天体的解体强度、反照率、直径进行了更好的约束，与不考虑大气层相比，解体强度相差40％-60％；将反照率范围约束在0.01-0.035，该范围远小于0.147的平均反照率；直径范围5-6m，远大于使用平均反照率估算的1.69m。
(3)使用本方法计算得到2022 EB5的轨道半长轴、偏心率、轨道倾角，比不考虑大气层得到的轨道根数更接近JPL解。大气层对近地点俯角、升交点赤经的影响很小。但总体而言，考虑大气层对轨道计算非常重要，不可忽略。
本发明的创新点在于：
1、充分考虑地球大气层对小天体各物理量的影响，对小天体物理特性进行约束。
2、将遗传算法与煎饼模型相结合对小天体的未知变量进行求解。
3、根据小天体大气进入过程计算小天体的日心轨道根数，对小天体观测轨道的可靠性与小天体类型进行验证。
本发明还可提供的一种计算机设备，包括：至少一个处理器、存储器、至少一个网络接口和用户接口。该设备中的各个组件通过总线系统耦合在一起。可理解，总线系统用于实现这些组件之间的连接通信。总线系统除包括数据总线之外，还包括电源总线、控制总线和状态信号总线。
其中，用户接口可以包括显示器、键盘或者点击设备(例如，鼠标，轨迹球(trackball)、触感板或者触摸屏等。
可以理解，本申请公开实施例中的存储器可以是易失性存储器或非易失性存储器，或可包括易失性和非易失性存储器两者。其中，非易失性存储器可以是只读存储器(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 DataRate SDRAM，DDRSDRAM)、增强型同步动态随机存取存储器(Enhanced SDRAM，ESDRAM)、同步连接动态随机存取存储器(Synchlink DRAM，SLDRAM)和直接内存总线随机存取存储器(Direct Rambus RAM，DRRAM)。本文描述的存储器旨在包括但不限于这些和任意其它适合类型的存储器。
在一些实施方式中，存储器存储了如下的元素，可执行模块或者数据结构，或者他们的子集，或者他们的扩展集：操作系统和应用程序。
其中，操作系统，包含各种系统程序，例如框架层、核心库层、驱动层等，用于实现各种基础业务以及处理基于硬件的任务。应用程序，包含各种应用程序，例如媒体播放器(Media Player)、浏览器(Browser)等，用于实现各种应用业务。实现本公开实施例方法的程序可以包含在应用程序中。
在本上述的实施例中，还可通过调用存储器存储的程序或指令，具体的，可以是应用程序中存储的程序或指令，处理器用于：
执行上述方法的步骤。
上述方法可以应用于处理器中，或者由处理器实现。处理器可能是一种集成电路芯片，具有信号的处理能力。在实现过程中，上述方法的各步骤可以通过处理器中的硬件的集成逻辑电路或者软件形式的指令完成。上述的处理器可以是通用处理器、数字信号处理器(Digital Signal Processor，DSP)、专用集成电路(Application Specific IntegratedCircuit，ASIC)、现场可编程门阵列(Field Programmable Gate Array，FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件。可以实现或者执行上述公开的各方法、步骤及逻辑框图。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。结合上述公开的方法的步骤可以直接体现为硬件译码处理器执行完成，或者用译码处理器中的硬件及软件模块组合执行完成。软件模块可以位于随机存储器，闪存、只读存储器，可编程只读存储器或者电可擦写可编程存储器、寄存器等本领域成熟的存储介质中。该存储介质位于存储器，处理器读取存储器中的信息，结合其硬件完成上述方法的步骤。
可以理解的是，本发明描述的这些实施例可以用硬件、软件、固件、中间件、微码或其组合来实现。对于硬件实现，处理单元可以实现在一个或多个专用集成电路(Application Specific Integrated Circuits，ASIC)、数字信号处理器(Digital SignalProcessing，DSP)、数字信号处理设备(DSP Device，DSPD)、可编程逻辑设备(ProgrammableLogic Device，PLD)、现场可编程门阵列(Field-Programmable Gate Array，FPGA)、通用处理器、控制器、微控制器、微处理器、用于执行本申请所述功能的其它电子单元或其组合中。
对于软件实现，可通过执行本发明的功能模块(例如过程、函数等)来实现本发明技术。软件代码可存储在存储器中并通过处理器执行。存储器可以在处理器中或在处理器外部实现。
本发明还可提供一种非易失性存储介质，用于存储计算机程序。当该计算机程序被处理器执行时可以实现上述方法实施例中的各个步骤。
最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (10)

1.一种求解撞击地球小天体的物理特性的方法，所述方法为：根据小天体进入大气后的物理数据，结合煎饼模型与遗传算法，考虑大气层对小天体各物理量以及轨道的影响，计算小天体的物理特性，利用已知小天体类型和物理特性对照表对计算结果进行筛选，从而得到待求解小天体的物理特性；其中，小天体的物理特性包括：质量、半径、密度和解体强度。
2.根据权利要求1所述的求解撞击地球小天体的物理特性的方法，其特征在于，所述方法具体包括：
步骤1：获取小天体进入大气信息；获取的信息包括峰值亮度时刻、峰值亮度时刻位置、峰值亮度时刻速度和撞击总能量；
步骤2：利用小天体峰值亮度时刻的物理量，通过四阶龙格-库塔方法对小天体飞行方程、烧蚀方程和半径变化方程构成的微分方程组进行逆向数值积分，得到小天体进入大气层边界处的速度矢量
Figure FDA0003727821170000011
位置矢量
Figure FDA0003727821170000012
半径R0、半径变化率dR0和质量m0；通过遗传算法计算小天体的解体强度S0、小天体在峰值亮度时刻的半径Rp、半径变化率dRp和质量mp；
步骤3：根据小天体类型与解体强度的对应关系，筛选出可能正确解。
3.根据权利要求2所述的求解撞击地球小天体的物理特性的方法，其特征在于，所述小天体飞行方程为：
Figure FDA0003727821170000013
其中，
Figure FDA0003727821170000014
为小天体相对地球的矢径；r为
Figure FDA0003727821170000015
的模，是小天体到地球的距离；m为小天体的质量；μ为地球引力常数，取3.986×1014m3/s2；t为时间；
Figure FDA0003727821170000016
表示小天体受到的大气阻力，满足下式：
Figure FDA0003727821170000017
其中，Cd为阻力系数，取值为1；Vr为小天体相对地球大气层的速度大小，亦为小天体在ECEF坐标系下的速度大小；S为迎风截面积；
Figure FDA0003727821170000018
为单位矢量，表示
Figure FDA0003727821170000019
的方向；ρa为大气密度。
4.根据权利要求2所述的求解撞击地球小天体的物理特性的方法，其特征在于，所述烧蚀方程为：
Figure FDA0003727821170000021
其中，m为小天体的质量；t为时间；Ch为烧蚀系数；Qh为烧蚀热；ρa为大气密度；Vr为小天体相对地球大气层的速度大小，亦为小天体在ECEF坐标系下的速度大小；S为迎风截面积。
5.根据权利要求2所述的求解撞击地球小天体的物理特性的方法，其特征在于，所述半径变化方程为：
Figure FDA0003727821170000022
其中，R为小天体的半径；t为时间；ρa为大气密度；Vr为小天体相对地球大气层的速度大小，亦为小天体在ECEF坐标系下的速度大小；ρm为小天体的密度。
6.根据权利要求2所述的求解撞击地球小天体的物理特性的方法，其特征在于，所述逆向数值积分计算公式为：
Figure FDA0003727821170000023
其中，Xb表示逆向数值积分结果；hp为小天体峰值亮度时刻的高度；hb为大气层边界高度hb＝100km；h为小天体高度；
Figure FDA0003727821170000024
t为时间；X代指小天体被积分的物理量，包括矢径
Figure FDA0003727821170000025
速度矢量
Figure FDA0003727821170000026
Figure FDA0003727821170000027
质量m，半径R，半径变化率
Figure FDA0003727821170000028
Figure FDA0003727821170000031
逆向数值积分过程由四阶龙格-库塔方法实现，最终得到小天体在大气层边界处ECI坐标系下的速度矢量
Figure FDA0003727821170000032
位置矢量
Figure FDA0003727821170000033
半径R0、半径变化率dR0和质量m0。
7.根据权利要求6所述的求解撞击地球小天体的物理特性的方法，其特征在于，进行逆向数值积分时，未知变量包括：小天体峰值亮度时刻的质量mp、半径Rp及半径变化率dRp，小天体的密度ρm、解体强度S0、烧蚀热Qh，模型参数N；
其中，N取为4，当N取4无解时，N取2；
根据彗星、C类小行星、S类小行星、M类小行星四类小天体的类型与密度ρm、烧蚀热Qh的对应关系，取四组密度ρm、烧蚀热Qh分别求解；
彗星对应密度ρm为1000kg/m3，对应烧蚀热Qh为2.5×106J/kg；
C类小行星对应密度ρm为2000kg/m3，对应烧蚀热Qh为5×106J/kg；
S类小行星对应密度ρm为3500kg/m3，对应烧蚀热Qh为8×106J/kg；
M类小行星对应密度ρm为7900kg/m3，对应烧蚀热Qh为8×106J/kg；
另外四个未知量S0、mp、Rp及dRp，作为优化参数，利用遗传算法进行求解；
遗传算法的四个约束条件分别为：
(1)计算得到的大气层边界处的半径变化率为0，即dR0＝0；
(2)计算得到的大气层边界处的半径与质量满足球体假设，即R(m)-R0＝0，其中
Figure FDA0003727821170000034
(3)计算得到的大气层边界处的动能与获取的小天体进入大气信息中的撞击总能量相等，即E(V0,m0)-E0＝0，其中
Figure FDA0003727821170000041
(4)计算得到的大气层边界处半径应与峰值亮度时刻半径Rp满足煎饼模型的比例关系，即NR0-Rp＝0；
其中，E0为获取到的总撞击能量；m0为大气层边界处小天体质量；R0为大气层边界处小天体半径；ρap为小天体峰值亮度时刻大气密度；Vp为小天体峰值亮度时刻速度；
遗传算法中四个优化参数S0、mp、Rp及dRp取值范围分别为：
S0最大值为Pp，最小值为Pp/30，
Figure FDA0003727821170000042
mp最大值为2me，最小值为0.01me，
Figure FDA0003727821170000043
Rp最大值为NR0max，最小值为R0min，
Figure FDA0003727821170000044
Figure FDA0003727821170000045
dRp最大值为50N，最小值为0.1N；
遗传算法中目标函数f由约束(1)至(4)加权组合而成，公式如下：
Figure FDA0003727821170000046
其中，a1、a2、a3、a4为四个约束的权重，初始值分别为10、100、0.01、50；
求得的解是否符合要求的判断标准为以下四项全部满足：
R0-R(m0)/R0＜1％
NR0-Rp/Rp＜1％
m0-m(E0,V0)/m0＜0.1％
dR0＜0.001
通过遗传算法最终计算得到小天体的解体强度S0、小天体在峰值亮度时刻的半径Rp、半径变化率dRp和质量mp。
8.根据权利要求2所述的求解撞击地球小天体的物理特性的方法，其特征在于，所述方法还包括：根据计算出的小天体在大气层边界处的物理特性和状态向量计算小天体日心轨道的轨道根数，包括：径向速度、比角动量、轨道倾角、节点线向量、升交点赤经、偏心率、近地点俯角、真近点角和轨道半长轴。
9.根据权利要求8所述的求解撞击地球小天体的物理特性的方法，其特征在于，所述方法还包括：求出小天体轨道后，可根据小天体相对木星的蒂塞朗常数TJ，远日点距离Q、近日点距离q及木星半长轴aJ的距离之间的关系，对小天体的类型进行验证。
10.根据权利要求8所述的求解撞击地球小天体的物理特性的方法，其特征在于，所述方法还包括：如果已知小天体的绝对星等，根据求得的直径，反推小天体的反照率，进一步根据不同类别小天体的反照率对小天体材质进行验证或约束，公式如下：
Figure FDA0003727821170000051
其中，H为小天体的绝对星等，D为小天体的直径，pv为反照率。
CN202210776924.5A 2022-07-04 2022-07-04 一种求解撞击地球小天体的物理特性的方法 Active CN115292657B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210776924.5A CN115292657B (zh) 2022-07-04 2022-07-04 一种求解撞击地球小天体的物理特性的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210776924.5A CN115292657B (zh) 2022-07-04 2022-07-04 一种求解撞击地球小天体的物理特性的方法

Publications (2)

Publication Number Publication Date
CN115292657A true CN115292657A (zh) 2022-11-04
CN115292657B CN115292657B (zh) 2023-06-23

Family

ID=83821529

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210776924.5A Active CN115292657B (zh) 2022-07-04 2022-07-04 一种求解撞击地球小天体的物理特性的方法

Country Status (1)

Country Link
CN (1) CN115292657B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116306185A (zh) * 2023-02-14 2023-06-23 北京理工大学 基于不规则离散元的陨星进入地球大气层后运动仿真方法及系统

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
RU2012108021A (ru) * 2012-03-05 2013-09-10 Александр Викторович Багров Способ ударного воздействия на опасные космические объекты и устройство для его осуществления
CN110147598A (zh) * 2019-05-10 2019-08-20 上海理工大学 基于图像处理的超高速撞击碎片云建模及损伤评估方法
CN113596332A (zh) * 2021-07-29 2021-11-02 中国科学院国家天文台长春人造卫星观测站 一种火流星事件光电监测系统
US11312512B1 (en) * 2019-03-04 2022-04-26 United States Of America As Represented By The Secretary Of The Air Force Early warning reentry system comprising high efficiency module for determining spacecraft reentry time

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
RU2012108021A (ru) * 2012-03-05 2013-09-10 Александр Викторович Багров Способ ударного воздействия на опасные космические объекты и устройство для его осуществления
US11312512B1 (en) * 2019-03-04 2022-04-26 United States Of America As Represented By The Secretary Of The Air Force Early warning reentry system comprising high efficiency module for determining spacecraft reentry time
CN110147598A (zh) * 2019-05-10 2019-08-20 上海理工大学 基于图像处理的超高速撞击碎片云建模及损伤评估方法
CN113596332A (zh) * 2021-07-29 2021-11-02 中国科学院国家天文台长春人造卫星观测站 一种火流星事件光电监测系统

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
SHUJUAN GENG等: "On the capture of small stony asteroids into the Earth\'s orbit by atmospheric grazing" *
柳森;党雷宁;赵君尧;白智勇;黄洁;李毅;石义雷;: "小行星撞击地球的超高速问题" *
耿淑娟等: "小天体撞击地球大气层的空爆问题研究" *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116306185A (zh) * 2023-02-14 2023-06-23 北京理工大学 基于不规则离散元的陨星进入地球大气层后运动仿真方法及系统
CN116306185B (zh) * 2023-02-14 2023-11-03 北京理工大学 基于不规则离散元的陨星进入地球大气层后运动仿真方法及系统

Also Published As

Publication number Publication date
CN115292657B (zh) 2023-06-23

Similar Documents

Publication Publication Date Title
Michel et al. Science case for the asteroid impact mission (AIM): a component of the asteroid impact & deflection assessment (AIDA) mission
McNamara et al. Meteoroid Engineering Model (MEM): A meteoroid model for the inner solar system
Sánchez et al. Impact hazard protection efficiency by a small kinetic impactor
Korycansky et al. Modeling crater populations on Venus and Titan
Capannolo et al. Challenges in LICIA Cubesat trajectory design to support DART mission science
Newman et al. Impact erosion of planetary atmospheres: Some surprising results
CN115292657A (zh) 一种求解撞击地球小天体的物理特性的方法
CN107367941A (zh) 基于非线性增益的高超声速飞行器攻角观测方法
CN114936471A (zh) 一种基于并行计算的航天器碰撞预警分层快速筛选方法
Dvorak et al. Planetary systems and the formation of habitable planets
Kaplinger et al. GPU accelerated 3D modeling and simulation of a blended kinetic impact and nuclear subsurface explosion
Brykina et al. On the mass distribution of fragments of an asteroid disrupted in the Earth’s atmosphere
Dachwald et al. Solar sail trajectory optimization for intercepting, impacting, and deflecting near-earth asteroids
CN115848648A (zh) 一种动能撞击小行星的伪谱轨道优化方法和装置
ReVelle NEO fireball diversity: energetics-based entry modeling and analysis techniques
CN107480347A (zh) 一种分离体散布特性预示方法
Carrico et al. Investigating orbital debris events using numerical methods with full force model orbit propagation
Shuvalov et al. Numerical simulation of the LCROSS impact experiment
Beech et al. The chant meteor procession of 1913–Towards a descriptive model
Boone et al. Artificial debris collision risk following a catastrophic spacecraft mishap in lunar orbit
CN114997043A (zh) 一种撞击近地小行星的方法
Colombo Optimal trajectory design for interception and deflection of Near Earth Objects
Kaplinger et al. Preliminary results for high-fidelity modeling and simulation of orbital dispersion of asteroids disrupted by nuclear explosives
CN117826616B (zh) 一种基于序列凸优化的飞行器快速轨迹规划方法及装置
Prabhu et al. Physics-Based Modeling of Meteor Entry and Breakup

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant