行星防御专栏

危地小行星撞击风险预报系统−短临预报场景(SHARP-IIS)设计与应用

  • 耿淑娟 , 1 ,
  • 赵洁凝 1, 2 ,
  • 王楷铎 1, 2 ,
  • 李明涛 , 1, 2, *
展开
  • 1 中国科学院国家空间科学中心, 北京100190
  • 2 中国科学院大学, 北京100190
(1982-),男,研究员,主要研究方向为小行星防御与利用、航天器轨道优化设计。通信地址:中国科学院国家空间科学中心(100190)电子邮箱:

(1996-),女,特别研究助理,主要研究方向为小行星防御与利用。通信地址:中国科学院国家空间科学中心(100190)电子邮箱:

网络出版日期: 2025-03-17

基金资助

科工局空间碎片与近地小行星防御科研专项-近地小行星时空协同监测与目标特性反演技术研究项目(KJSP202320105)

版权

版权所有 © 2024 空间科学与试验学报编辑部

Design and Application of Prediction and Analysis System for Imminent Impacts

  • Shujuan GENG , 1 ,
  • Jiening ZHAO 1, 2 ,
  • Kaiduo WANG 1, 2 ,
  • Mingtao LI , 1, 2, *
Expand
  • 1 National Space Science Center, CAS, Beijing 100190, China
  • 2 University of Chinese Academy of Sciences, Beijing 100190, China

Online published: 2025-03-17

Copyright

Copyright © 2024 Journal of Space Science and Experiment. All rights reserved.

摘要

近地小行星撞击是全人类共同面临的重大潜在威胁。通过监测预警发现撞击概率大的近地小行星,对其开展陨落预报,是防御小行星撞击灾害的重要前提。基于近地小行星短临预报任务场景的核心需求,中国科学院国家空间科学中心自主设计并开发了危地小行星撞击风险预报系统(简称“锐警”),其中的短临预报场景能够对进入地球一定范围内、存在撞击概率的近地小行星开展基于真实观测数据的轨道确定、轨道递推、陨落预报的全流程分析,为小行星撞击灾害防御提供参考。本文介绍了SHARP-IIS系统各功能模块的基本原理,并借助小行星2024 RW1、小行星2024 UQ两次撞击事件展示本系统在不同观测弧长下的表现。SHARP-IIS系统在观测弧段为7 h的2024 RW1场景中实现了较好的预测精度,与近地天体研究中心(CNEOS)记录的空爆信息相比,SHARP-IIS计算的空爆位置误差小于20 km,空爆时刻最小误差约1 s。

本文引用格式

耿淑娟 , 赵洁凝 , 王楷铎 , 李明涛 . 危地小行星撞击风险预报系统−短临预报场景(SHARP-IIS)设计与应用[J]. 空间科学与试验学报, 2024 , 1(4) : 20 -35 . DOI: 10.19963/j.cnki.2097-4302.2024.04.003

Abstract

Near-Earth asteroid impacts represent a significant potential threat faced by all of humanity. Monitoring and providing early warnings to identify near-Earth asteroids with a high impact probability, as well as conducting impact predictions, are essential prerequisites for planetary defense. According to the core requirements of short-term impact forecasting for near-Earth asteroid, National Space Science Center, CAS independently designed and developed the System for Hazardous Asteroids Risk Prediction (SHARP). The Imminent Impactor Scenario of SHARP (SHARP-IIS) performs orbit determination, orbit propagation, and impact prediction for near-Earth asteroids that reach a certain proximity to Earth based on observational data to provide critical references for planetary defense. This paper introduces the basic principles of each functional module of the SHARP-IIS system and demonstrates its performance through the impact events of asteroid 2024 RW1 and asteroid 2024 UQ with different observation arc lengths. In the scenario of asteroid 2024 RW1 with an observation arc of seven hours, the SHARP-IIS system achieved a relatively high-level prediction accuracy. Compared with the airburst data recorded by the Center for Near-Earth Object Studies (CNEOS), the calculated airburst location error by SHARP-IIS was less than 20 km, and the minimum airburst timing error was around 1 s.

0 引言

近地小行星撞击地球是全人类共同面临的威胁与挑战。直径10 km的小行星足以引发全球性灾害[1]。普遍认为,6500万年前的白垩纪物种大灭绝事件便是由一颗直径10 km左右的小行星或彗星撞击地球所致[2-3]。中小尺寸小行星撞击地球可能引发区域性灾害。1908年,一颗直径约50 m的小行星在俄罗斯通古斯卡河上空发生空爆或飞掠,产生的冲击波与热辐射效应导致超过2000 km2的森林损毁[4-5]。2013年,一颗直径约20 m的小行星在俄罗斯车里雅宾斯克地区发生空爆,致使3000多栋房屋受损,1500多人受伤[6-7]。由此可见,小行星撞击可能给人类带来严重灾害,开展小行星防御是保护人类文明的必然要求,也是保卫地球安全的重要举措。
对近地小行星开展监测预警是行星防御的重要环节。提前发现具有潜在威胁的小行星,有效确定其轨道并分析撞击概率,从而对小行星可能带来的灾害进行评估,在必要时采取灾害应对与减缓措施。根据NEOMOD 3近地天体尺寸分布模型[8],以及近地天体研究中心(CNEOS)的统计数据 1,截至2024年11月19日,直径1 km以上的近地小行星预计约830±60颗,已发现868颗,发现率近100%;直径100 m以上的近地小行星约30000±3000颗,已发现13587颗,发现率约45%;直径30 m以上的近地小行星约368000±60000颗,已发现25979颗,发现率约7%。由于小行星撞击地球的频率随尺寸的减小而增加,在可预见的未来,发现率较低的几十米及以下的近地小行星更有可能撞击地球。据估计,直径约1 m的小行星撞击事件约每两周发生一次,直径10 m的小行星撞击事件每10年发生一次,这部分小行星由于尺寸较小、反照率较低,不易在撞击地球前被观测到[9-10]。截至2024年12月25日,在撞击地球前被发现的小行星共11颗,分别为2008 TC3[11-12]、2014 AA[13]、2018 LA、2019 MO、2022 EB5[14]、2022 WJ1[15]、2023 CX1[16]、2024 BX1[17]、2024 RW1[18]、2024 UQ、2024 XA1。这些小行星被发现时距离地球较近,从发现至撞击的时间较短,最长约21 h,最短仅1 h。在短时间内预测小行星的撞击轨道及进入大气后的位置是短临预报场景的核心需求。
为了更好地预警小行星撞击威胁,欧空局(ESA)建立了近地天体扫描系统(NEOScan)[19-20],该系统每两分钟从国际小行星中心(MPC)的近地天体确认页(NEOCP)读取近地小行星观测数据,识别近地小行星类别,将NEOCP的数据进行筛选与管理,对即将撞击的小行星进行预警,计算撞击概率与风险等级,并开展跟踪观测。此外,ESA还开发了Meerkat小行星短临预警预报系统[21],针对发现的即将撞击地球的小行星开展快速定轨与威胁评估。美国航天局(NASA)开发了Sentry小行星预警软件系统[22]和Scout新发现小行星即时撞击预警系统[23]。Sentry系统能够自动读取MPC发布的小行星观测数据,并根据观测数据计算小行星轨道,评估撞击概率。Scout系统从MPC的NEOCP读取实时观测数据,并根据观测数据确定小行星的轨道、计算星历、估算撞击概率。在Scout系统建立前,NASA已经实现了对小行星2008 TC3陨落区域的预报。该小行星是第一颗人类成功预警的撞击地球的小行星,NASA计算了该小行星大气进入轨迹在地面的投影,估算了不同质量的陨石可能陨落的位置。尽管陨石的陨落位置在沿投影方向上存在偏差,但此次预报是人类对撞击地球小行星的首次预报,检验了小行星轨道确定与落点预测流程,提高了陨石搜集效率。Scout系统在建成后,又先后实现了对小行星2023 CX1、小行星2024 BX1的陨落预报,从发现到完成陨落预报的所需时间也从最初的20.5 h缩短为1.5 h。国内,中国科学院紫金山天文台具有丰富的小行星观测与定轨经验,建立了考虑太阳光压等多种效应的高精度轨道动力学模型,能够基于观测数据对小行星的轨道进行确定[24]。党雷宁等[25]开发了小行星进入与撞击效应分析评估软件AICA,在已知大气层边界处参数的情况下,能够计算小行星的能量沉积过程、空爆高度和地面损伤范围。然而,国内目前尚未形成针对短临预报场景的集轨道确定与陨落预报为一体的分析系统。
基于近地小行星短临预报工作场景的客观需求,中国科学院国家空间科学中心研发了危地小行星撞击风险预报系统(System for Hazardous Asteroids Risk Prediction,SHARP),简称“锐警”系统。根据不同的小行星撞击场景,“锐警”系统分短临预报场景(Imminent Impactor Scenario, SHARP-IIS)和威胁预警场景(Threat Warning Scenario, SHARP- TWS)。其中,锐警-短临预报场景针对距离地球较近的小行星,能够利用MPC公布的小行星观测数据,对撞击地球概率较大的近地小行星,开展轨道确定、轨道递推、陨落预报的全流程分析。本文旨在介绍SHARP-IIS系统的基本原理,计算流程,以及该系统在近两次撞击场景(小行星2024 RW1和2024 UQ)中的应用。

1 原理与方法

小行星撞击短临预报与分析围绕三个核心过程展开:轨道确定、轨道递推、陨落预报。本节介绍SHARP-IIS采用的轨道确定方法、轨道动力学模型,以及大气进入过程计算方法。

1.1 轨道确定方法

轨道确定是指从观测数据确定某一时刻小行星轨道根数的过程。在SHARP-IIS中,轨道确定主要包括两个过程,即初定轨与轨道改进。

1.1.1 初定轨

初定轨的目的是根据观测数据快速求得初始解,为轨道改进过程提供合理的初值。
MPC提供的观测数据为不同时刻(记为$ {t}_{i} $)小行星相对观测站的经度$ {\alpha }_{i} $与纬度$ {\sigma }_{i} $,需要根据观测数据计算小行星相对地球的位置矢量。
首先,将小行星相对观测站的经纬度转化为小行星相对于每个观测站的位置矢量的单位向量$ {\boldsymbol{L}}_{i} $,该过程可通过式(1)实现。其中,$ \lambda $$ \mu $$ \nu $分别表示单位向量在站心赤道J2000坐标系下x、y、z轴方向的分量。
$ {\boldsymbol{L}}_{i}=\left(\genfrac{}{}{0pt}{}{\lambda }{\begin{array}{c}\mu \\ \nu \end{array}}\right)=\left(\genfrac{}{}{0pt}{}{\mathrm{cos}{\sigma }_{i}\mathrm{cos}{\alpha }_{i}}{\begin{array}{c}\mathrm{cos}{\sigma }_{i}\mathrm{sin}{\alpha }_{i}\\ \mathrm{sin}{\sigma }_{i}\end{array}}\right) $
在每个观测站相对于地心的位置矢量$ {\boldsymbol{R}}_{i} $已知的情况下,可得到小行星相对于地心的位置矢量$ {\boldsymbol{r}}_{i} $
$ {\boldsymbol{r}}_{i}={\rho }_{i}{\boldsymbol{L}}_{i}+{\boldsymbol{R}}_{i} $
式中,$ {\rho }_{i} $$ {t}_{i} $时刻小行星相对观测站的距离。
设首末时刻分别为$ {t}_{0} $$ {t}_{\mathrm{f}} $,若已知$ {\rho }_{0} $$ {\rho }_{\mathrm{f}} $以及观测站位置,则根据式(2)可求得对应时刻的$ {\boldsymbol{r}}_{0} $$ {\boldsymbol{r}}_{\mathrm{f}} $,并进一步通过求解兰伯特问题计算两个时刻小行星相对地球的速度矢量,进而求得小行星轨道,如式(3)所示。不同观测时刻下观测站的位置可从MPC提供的观测数据中获得。
$ \left({\boldsymbol{v}}_{\mathrm{f}},{\boldsymbol{v}}_{0}\right)=\mathrm{l}\mathrm{a}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r}\mathrm{t}\left({{\boldsymbol{r}}_{\mathrm{f}},\boldsymbol{r}}_{0},{t}_{\mathrm{f}}-{t}_{0}\right) $
由于初定轨阶段的目的是尽快求得轨道初始解,因此,为了在保障合理性的同时提高计算效率,此阶段所计算的小行星轨道需要满足:在不同观测时刻,由轨道计算得到的小行星相对观测站的位置矢量的单位向量与根据观测数据直接计算得到的位置矢量的单位向量之间的夹角之和最小。
因此,初定轨可以归纳为最优化问题,如式(4)所示。其中,$ {\rho }_{0} $$ {\rho }_{\mathrm{f}} $为优化变量,$ {{\boldsymbol{L}}_{i}}^{*} $$ {\rho }_{0} $$ {\rho }_{\mathrm{f}} $的函数,其含义为通过求解兰伯特问题得到的$ {t}_{i} $时刻的小行星相对于观测站的位置矢量的单位向量。
$ \underset{\left({\rho }_{0},{\rho }_{{{\mathrm{f}}}}\right)}{\mathrm{min}}\sum _{i=1}^{n}{\mathrm{cos}}^{-1}\left[{\boldsymbol{L}}_{i}\cdot {{\boldsymbol{L}}_{i}}^{*}\left({\rho }_{0},{\rho }_{\mathrm{f}}\right)\right] $
优化变量的初值及优化变量范围可通过视星等$ {V}_{i} $、绝对星等$ H $、小行星相对观测平台的位置矢量$ {\boldsymbol{\rho }}_{i}={\rho }_{i}{\boldsymbol{L}}_{i} $和小行星相对太阳的位置矢量$ {\boldsymbol{r}}_{{\mathrm{si}}} $之间的夹角$ \kappa $进行确定,如式(5)所示。假设小行星完全顺光观测,则有$ \kappa =0 $$ G $为斜率参数,此处取0.15。
$\begin{aligned}&H={V}_{i}-5\mathrm{log}\left({\boldsymbol{r}}_{{\mathrm{si}}}\cdot {\boldsymbol{\rho }}_{\boldsymbol{i}}\right)-\\&\quad\;\;\;\;2.5\mathrm{log}\left[\left(1-G\right){\varPhi}_{1}\left(\kappa \right)+G{\varPhi}_{2}\left(\kappa \right)\right]\\&{\varPhi}_{1}\left(\kappa \right)={\mathrm{e}}^{-3.339{\mathrm{tan}}^{0.63}\left(\kappa /2\right)} \\&{\varPhi}_{2}\left(\kappa \right)={\mathrm{e}}^{-1.87{\mathrm{tan}}^{1.22}\left(\kappa /2\right)} \end{aligned}$
利用遗传算法求解式(4)描述的优化问题,可得到初始轨道解。

1.1.2 轨道改进

轨道改进将在初定轨结果的基础上开展。与初定轨阶段不同,轨道改进要求小行星轨道满足所有时刻计算得到的观测方向与真实的观测方向之间残差的平方和最小。故轨道改进可归纳为以下最优化问题:
$ \underset{\left({\boldsymbol{r}}_{0},{\boldsymbol{v}}_{0}\right)}{\mathrm{min}}\sum _{i=1}^{n}{\left[{y}_{i}-f\left({{t}_{i},\boldsymbol{r}}_{0},{\boldsymbol{v}}_{0}\right)\right]}^{2} $
式中,$ \left({\boldsymbol{r}}_{0},{\boldsymbol{v}}_{0}\right) $是优化变量,为$ {t}_{0} $时刻小行星的位置与速度,其初值为初定轨阶段计算得到的轨道初始解;$ {y}_{i}=\left({\alpha }_{i},{\sigma }_{i}\right) $$ {t}_{i} $时刻的观测数据;$ f\left({{t}_{i},\boldsymbol{r}}_{0},{\boldsymbol{v}}_{0}\right)= \left({{\alpha }_{i}}^{*},{{\sigma }_{i}}^{*}\right) $为根据$ {t}_{0} $时刻的位置速度,考虑光行差后通过轨道递推计算得到的小行星相对观测站的经度、纬度。此阶段采用的轨道动力学模型将在第1.2节进行介绍。
采用最小二乘法对该问题进行求解。最终可得到轨道标称解$ \left({{\boldsymbol{r}}_{0}}^{*},{{\boldsymbol{v}}_{0}}^{*}\right) $。轨道的不确定性由位置、速度的协方差矩阵$ \mathbf{c}\mathbf{o}\mathbf{v}\left(\boldsymbol{\beta }\right) $进行描述,其计算方式如式(7)所示。
$\begin{aligned}{\boldsymbol{J}}_{i,j}=&\frac{\partial f\left({t}_{0},\boldsymbol{\beta }\right)}{\partial {\boldsymbol{\beta }}_{i}},i=\mathrm{1,2},\cdots ,6;j=\mathrm{1,2},\cdots ,6 \\&\quad\quad\;\;\;\mathbf{c}\mathbf{o}\mathbf{v}\left(\boldsymbol{\beta }\right)={\sigma }^{2}{\left({\boldsymbol{J}}^{\mathrm{T}}\boldsymbol{J}\right)}^{-1}\end{aligned}$
式中,$ \boldsymbol{J} $为雅可比矩阵;$ \boldsymbol{\beta }=\left({\boldsymbol{r}}_{0},{\boldsymbol{v}}_{0}\right) $$ {\sigma }^{2} $为残差的方差。
由于小行星观测条件与观测设备性能等因素的限制,MPC公布的观测数据质量水平可能存在差异,有些观测数据的存在可能会对轨道确定结果带来较大的误差,从而影响轨道确定的精度。因此,在初次获得定轨结果后,本系统设置了是否根据轨道解对观测数据进行筛选的选择按钮。若选择“是”,则在轨道改进阶段,系统将通过分析不同时刻的数据残差,剔除残差过大(>3$ \sigma $)的观测数据,而后利用未被剔除的数据重新开展轨道改进,求解式(6)。此时,第二次轨道改进的初始解为首次轨道改进得到的解,而非初定轨的解。若选择“否”,则系统将不会根据残差对观测数据进行筛选。

1.2 轨道动力学模型

轨道动力学模型是开展轨道递推过程的基础。在SHARP-IIS的分析中,主要存在两个轨道递推过程,第一个过程为轨道确定中轨道改进阶段的轨道递推,根据初始时刻小行星的位置、速度开展轨道递推,以获得不同时刻的位置、速度;第二个过程为获得小行星初始时刻的轨道解后,将小行星轨道递推至其近地点(不撞击地球的场景)或地球大气层边界处(撞击地球的场景),目前SHARP-IIS系统处理的场景以后者居多,此阶段的结果将作为大气进入阶段的输入。如无特殊说明,SHARP-IIS系统中的“轨道递推”默认为第二个过程。
在轨道改进以及远距离的轨道递推中,采用的轨道动力学模型考虑太阳的中心引力,八大行星、冥王星、月球的第三体引力摄动,相对论效应,如式(8)所示。
$ \frac{{\mathrm{d}}^{2}{\boldsymbol{r}}_{\mathbf{s}}}{\mathrm{d}{t}^{2}}=-\frac{{\mu }_{\mathrm{s}\mathrm{u}\mathrm{n}}}{{\left|{\boldsymbol{r}}_{\mathbf{s}}\right|}^{3}}{\boldsymbol{r}}_{\mathbf{s}}-\sum _{i=1}^{10}{\mu }_{i}\left(\frac{1}{{\left|{d}_{i}\right|}^{3}}{\boldsymbol{d}}_{\boldsymbol{i}}+\frac{1}{{\left|{{r}_{i}}'\right|}^{3}}{{\boldsymbol{r}}_{\boldsymbol{i}}}'\right)+{\boldsymbol{a}}_{\mathrm{R}\mathrm{E}} $
式中,$ {\boldsymbol{r}}_{\mathbf{s}} $为小行星相对太阳的位置矢量;$ {\mu }_{\mathrm{s}\mathrm{u}\mathrm{n}} $为太阳引力常数;$ {\boldsymbol{d}}_{\boldsymbol{i}} $为第$ i $颗第三体相对太阳的位置矢量,$ {{\boldsymbol{r}}_{\boldsymbol{i}}}' $为第$ i $颗第三体到小行星的位置矢量;$ i $从1到10分别代表水星、金星、地球、火星、木星、土星、天王星、海王星、冥王星、月球;$ {\mu }_{i} $为对应天体的引力常数;$ {\boldsymbol{a}}_{\mathrm{R}\mathrm{E}} $表示相对论效应。八大行星以及冥王星、月球星历采用JPL星历。
当小行星距离地球较近(如2008 TC3等11次成功预警的发现时距离撞击不足1天的事件)时,小行星主要受到地球、太阳、月球引力的作用,此时采用的轨道动力学模型如式(9)所示。
$ \frac{{\mathrm{d}}^{2}\boldsymbol{r}}{\mathrm{d}{t}^{2}}=-\frac{{\mu }_{\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{t}\mathrm{h}}}{{\left|\boldsymbol{r}\right|}^{3}}\boldsymbol{r}-\sum _{j=1}^{2}{\mu }_{j}\left(\frac{1}{{\left|{\mathrm{d}}_{j}\right|}^{3}}{\boldsymbol{d}}_{\boldsymbol{j}}+\frac{1}{{\left|{{r}_{j}}'\right|}^{3}}{{\boldsymbol{r}}_{\boldsymbol{j}}}'\right)+{\boldsymbol{a}}_{\mathrm{J}2} $
式中,j=1、j=2分别表示太阳与月球;$ {\boldsymbol{a}}_{\mathrm{J}2} $为地球非球形引力摄动。

1.3 大气进入模型

小行星在地球大气层内的行为主要包括减速、烧蚀、解体。在这一过程中,小行星速度、位置、质量随时间的变化可由式(10)求得。
$ \begin{array}{c}\dfrac{\mathrm{d}V}{\mathrm{d}t}=-\dfrac{{C}_{\mathrm{d}}\rho {AV}^{2}}{2m}-\dfrac{\mu }{{r}^{2}}\mathrm{sin}\theta\\\dfrac{\mathrm{d}\theta }{\mathrm{d}t}=\dfrac{{C}_{\mathrm{l}}\rho AV}{2m}-\dfrac{\mu }{{r}^{2}}\mathrm{cos}\theta +\dfrac{V}{r}\mathrm{cos}\theta \\\dfrac{\mathrm{d}h}{\mathrm{d}t}=V\mathrm{sin}\theta \\\dfrac{\mathrm{d}m}{\mathrm{d}t}=-\dfrac{{C}_{\mathrm{h}}\rho A{V}^{3}}{2Q} \\\dfrac{\mathrm{d}\lambda }{\mathrm{d}t}=\dfrac{V\mathrm{cos}\theta \mathrm{cos}\psi }{r\mathrm{cos}\varphi }\\\dfrac{\mathrm{d}\varphi }{\mathrm{d}t}=\dfrac{V\mathrm{cos}\theta \mathrm{sin}\psi }{r}\\\dfrac{\mathrm{d}\psi }{\mathrm{d}t}=-\dfrac{V\mathrm{cos}\theta \mathrm{cos}\psi \mathrm{tan}\varphi }{r} \end{array}$
式中,$ V $$ \theta $$ h $$ m $$ \lambda $$ \varphi $$ \psi $分别表示小行星的速度、进入角、高度、质量、经度、纬度、航向角;进入角为小行星速度矢量与当地水平线的夹角,航向角为速度在水平面的投影与东方向的夹角;$ {C}_{\mathrm{d}} $$ \rho $$ A $$ \mu $$ r $$ {C}_{\mathrm{l}} $$ Q $$ {C}_{\mathrm{h}} $分别为阻力系数、大气密度、小行星迎风截面积、地球引力常量、小行星地球轨道半径、升力系数、烧蚀热、热流系数。
当小行星前缘动压$ P=\rho {V}^{2} $超过其结构强度$ {S}_{0} $时,小行星开始解体。采用连续性解体模型模拟小行星的解体过程。该模型将解体后的小行星碎片组成的云团视为一个可以变形的整体,通过云团的横向形变模拟各碎片彼此分离的过程。在Chyba等[4]提出的模型中,当云团半径变为初始半径的$ N $倍时,认为小行星发生空爆。该云团的半径变化过程根据式(11)计算。
$ \frac{{\mathrm{d}}^{2}R}{\mathrm{d}{t}^{2}}=\frac{{C}_{\mathrm{d}}\rho {V}^{2}}{2{\rho }_{\mathrm{m}}R} $
式中,$ R $为小行星半径;$ {\rho }_{\mathrm{m}} $为小行星密度。

2 SHARP-IIS软件计算流程

小行星撞击短临预报与分析系统以软件的形式实现,该软件主要包括三个工作模块,分别为“轨道确定”模块、“轨道递推”模块、“大气进入”模块。其工作流程如图1所示。
图 1 SHARP-IIS计算流程

Fig.1 Calculation Process of SHARP-IIS

轨道确定模块主要功能是读取目标小行星的观测数据,并根据观测数据确定小行星轨道。其主要输入、输出参数如表1所示。
表 1 轨道确定模块的输入、输出参数

Table 1 The input and output parameters of orbit determination module

模块 输入参数 输出参数
轨道确定 1)观测数据文件路径 1)总观测数据条数
2)小行星编号 2)定轨使用的观测数据条数
3)绝对星等 3)剔除的观测数据条数
4)绝对星等上限 4)标称解轨道根数
5)绝对星等下限 5)标称解位置速度
6)观测数据组数 6)协方差矩阵
7)遗传算法种群数量 7)小行星日心轨道图
8)遗传算法最大迭代次数
9)蒙特卡罗仿真参数
10)是否对数据进行筛选
表1的输入参数中,观测数据文件由用户选取,当用户选取对应的观测数据文件后,软件将自动记录文件路径,并显示在对应的文本框中。小行星的绝对星等上限、下限用于初定轨阶段求解优化变量初值及对应的变量范围。观测数据组数为用户设置的定轨阶段使用的观测数据条数,当该条数超过观测数据文件中的数据总条数时,系统将自动替换为文件中的最大数据条数。蒙特卡罗仿真参数为蒙特卡罗仿真生成的场景数,用于轨道递推求解撞击概率、大气进入阶段计算空爆位置分布。在开展蒙特卡罗仿真时,将假设小行星位置速度分量服从联合正态分布,利用轨道标称解与协方差矩阵随机生成服从该分布的撞击场景,场景数量由蒙特卡罗仿真参数确定。若仅需要标称解的计算结果,则可将该值设定为1,在后续的轨道递推与大气进入计算过程中,系统将不开展蒙特卡罗仿真,仅更新与标称解对应的计算结果。用户可根据实际需求设置是否在轨道改进阶段对观测数据进行筛选。
表1的输出参数中,总观测数据条数为观测数据文件中所有观测数据的条数。剔除的观测数据条数为轨道改进阶段因数据残差过大而被剔除的观测数据条数。若用户在是否进行数据筛选的选项中选择“否”,则剔除的观测数据条数将显示为0。标称解轨道根数为小行星的日心轨道根数。标称解位置速度为小行星在地心J2000坐标系下的位置速度。协方差矩阵为地心J2000坐标系下的位置与速度分量对应的协方差矩阵。
轨道递推模块主要功能是根据轨道确定得到的标称解及协方差矩阵,确定小行星撞击地球的概率,并对存在撞击地球概率的场景计算其地球大气层边界处的位置速度。其主要输入、输出参数如表2所示。
表 2 轨道递推模块的输入、输出参数

Table 2 The input and output parameters of orbit propagation module

模块 输入参数 输出参数
轨道递推 1)轨道动力学模型 1)轨道标称解大气层边界处位置速度(ECEF坐标系、ICRF
坐标系、经纬高)
2)最长递推时间
3)最长积分步长
4)相对积分容差 2)撞击概率
5)绝对积分容差 3)标称解轨道递推投影初始
时刻速度、位置分布图
6)蒙特卡罗仿真参数
7)标称解位置速度 4)大气层边界处参数分布图
表2的输入参数中,最长递推时间、最长积分步长、相对积分容差、绝对积分容差为对轨道动力学方程积分时的参数设置。蒙特卡罗仿真参数、标称解位置速度由轨道确定模块传递而来。在表2的输出参数中,撞击概率为蒙特卡罗仿真中能够到达地球大气层边界的场景数与总场景数的比值。
大气进入模块主要功能是基于轨道递推阶段提供的小行星在地球大气层边界处的位置速度,计算小行星大气进入过程,开展空爆位置预报。其主要输入、输出参数如表3所示。
表 3 大气进入模块的输入、输出参数

Table 3 The input and output parameters of atmospheric entry module

模块 输入参数 输出参数
大气
进入
1)大气层边界处标称解位置速度(经纬高) 1)标称解下空爆时刻的位置
2)蒙特卡罗仿真参数 2)空爆位置与参考点之间的距离
3)撞击概率 3)标称解大气进入轨迹及
空爆位置投影图
4)反照率 4)小行星强度分布图
5)小行星直径 5)小行星空爆位置分布图
6)小行星密度 6)空爆位置与参考点的
距离分布图
7)小行星质量
8)小行星强度
9)小行星强度下限
10) 小行星强度上限
11) 形变因子
12) 参考点位置
表3的输入参数中,大气层边界处标称解位置速度、蒙特卡罗仿真参数、撞击概率由轨道递推模块传递而来。在小行星反照率未知的情况下,通常假设反照率为平均反照率0.15,利用绝对星等、直径、反照率之间的经验关系计算小行星直径。在小行星类型未知的情况下,考虑到石质小行星在近地小行星中占比最高,因此,可假设小行星密度为石质小行星的典型密度3500 kg/m3,也可根据需要假设为其他类型,SHARP-IIS将根据输入的密度自动判断类型,并设置对应的烧蚀热。小行星质量根据直径与密度,按照球体的质量公式计算。小行星强度为计算标称解时采用的强度。小行星强度通常未知,由用户凭借经验自行设定标称解的强度值。在蒙特卡罗仿真中,将设定小行星强度上下限,假设强度在该范围内服从均匀分布。形变因子为连续性解体模型中定义变形终止条件时使用的模型参数,通常取为4。参考点位置为观测到的小行星空爆位置,该值可以来源于CNEOS,可以来源于地基火流星观测网,也可以设定为任意用户期望比较的点位。若无参考点位置,系统将自动设置为0。
表3的输出参数中,若蒙特卡罗仿真参数为1,则小行星强度分布、距离分布、空爆位置分布图将不显示。若参考点位置为0,则小行星空爆与参考点的距离分布图将不显示,标称解下空爆时刻与参考点的距离将显示为0。

3 正确性验证

利用小行星2024 RW1验证SHARP-IIS的正确性,并展示SHARP-IIS软件界面。
小行星2024 RW1于UTC时间2024年9月4日被发现,为第9颗撞击地球前被发现的小行星。该小行星随后在菲律宾东北方空发生空爆,此次空爆被CNEOS记录在火球数据集中,具体参数如表4所示。表4中2024 RW1的空爆时刻及空爆位置将作为验证SHARP-IIS正确性的主要依据。
表 4 CNEOS记录的小行星2024 RW1的空爆信息

Table 4 Airburst information of 2024 RW1 recorded by CNEOS

参数 数值
空爆时间(UTC) 2024-09-04 16:39:32
空爆时刻经度 122.9°E
空爆时刻纬度 18.0°N
空爆时刻高度 25.0 km
空爆时刻速度及ECEF系下速度分量 19.7 km/s
Vx 3.9 km/s
Vy −19.1 km/s
Vz 2.6 km/s
总辐射能 6.2×1010 J
撞击能量 0.2 kt TNT
根据2024年9月4日从MPC下载的2024 RW1的观测数据(可用数据共40条,观测弧段约7 h)后,利用SHARP-IIS系统对2024 RW1场景进行计算。在轨道确定页面设置相应输入参数,点击“轨道确定”,SHARP-IIS即开始轨道确定过程。定轨完成后,轨道确定结果将显示在右侧的结果面板中,如图2图3所示。
图 2 2024 RW1轨道确定结果-轨道根数

Fig.2 Orbit determination results of 2024 RW1-orbital elements

图 3 2024 RW1轨道确定结果-小行星日心轨道

Fig.3 Orbit determination results of 2024 RW1-heliocentric orbit diagram

图2“观测数据使用情况”面板可以发现,选择的观测文件中观测数据共40条,实际用于定轨的数据为29条。“标称解-轨道根数/位置速度”面板中显示了标称解下小行星日心黄道轨道根数、地心J2000坐标系下的位置速度,以及对应的儒略日。地心J2000坐标系下位置和速度在x轴、y轴、z轴三个方向的分量对应的协方差矩阵显示在“协方差矩阵”面板中。该矩阵反映了定轨阶段不确定性的大小。图3展示了小行星的日心轨道,并将小行星轨道与水星、金星、地球、火星、木星轨道进行了对比。其中,黄色实心球表示太阳,小行星轨道以红色实线表示,水星、金星、地球、火星、木星轨道分别用蓝色实线、柿色实线、橙色实线、紫色实线、绿色实线表示。
SHARP-IIS系统计算的轨道根数与美国喷气式实验室(JPL)计算的轨道根数对比如表5所示。为方便对比,文中将SHARP-IIS的轨道根数递推至JPL的历元时刻。对比可以发现,二者半长轴差异约0.001 au;偏心率绝对差异在小数点后四位;轨道倾角绝对差异在小数点后五位;升交点经度绝对差异在小数点后三位;近日点幅角、平近点角的绝对差异在小数点后两位。从图2的协方差矩阵可以看出,2024 RW1的位置、速度方差较小,意味着即使在后续采用蒙特卡罗方法对误差范围内的轨道进行模拟,在较短的递推时间内,小行星的轨道也不会过于发散,小行星进入大气后的落点将相对集中。
表 5 2024 RW1日心轨道根数对比

Table 5 Comparison of heliocentric orbital elements of 2024 RW1

参数 SHARP-IIS JPL
历元时刻UTC 2024-09-03 23:58:51 2024-09-03 23:58:51
半长轴 $ a $/au 2.50807445497218 2.507107997333645
偏心率 $ e $ 0.706871858429814 0.7067745723554399
轨道倾角 $ i $/($ \text{°} $) 0.528031355989482 0.5280529674740915
升交点经度 $ \mathrm{\Omega } $/($ \text{°} $) 162.456926718385 162.4574634817969
近日点幅角 $ \omega $/($ \text{°} $) 249.613909121206 249.6223582485825
平近点角 $ M $/($ \text{°} $) 349.195675950465 349.1881578987681
在轨道确定页面中点击“轨道递推”,可进入轨道递推页面。在该页面设置积分参数后点击“轨道递推”,SHARP-IIS即开始计算。2024 RW1的轨道递推结果如图4图6所示。蒙特卡罗仿真显示,2024 RW1的撞击概率为100%。图4中的红色标记为标称解下小行星从轨道确定的时刻至进入地球大气层时刻的轨道在地图上的投影,黑色空心方框表示进入大气时刻的位置,即此段轨迹投影的终点。
图 4 2024 RW1轨道递推结果-轨迹投影

Fig.4 Orbit propagation results of 2024 RW1-trajectory projection

图 6 2024 RW1轨道递推结果-大气进入参数分布

Fig.6 Orbit propagation results of 2024 RW1-distribution of atmospheric entry parameters

图5绘制了轨道递推初始时刻小行星地心J2000坐标系下位置与速度分量的分布,可看到大致呈正态分布,与蒙特卡罗仿真中生成初始虚拟轨道阶段所假设的分布一致。从该分布图中也可以看到,轨道的位置分布范围在百公里级,速度的分布范围约十米每秒级,意味着轨道确定过程收敛性较好,不确定性较低。图6中,蒙特卡罗仿真后大气层边界处的参数分布也同样较为集中,进入角、方向角的分布范围均在0.02°内,速度的分布范围小于0.01 km/s。纬度分布范围在0.02°内,经度分布范围在0.03°内,位置投影分布呈现出狭长的椭圆形。
图 5 2024 RW1轨道递推结果-位置速度分布

Fig.5 Orbit propagation results of 2024 RW1-distribution of position and velocity

在轨道递推界面点击“大气进入计算”,SHARP-IIS将弹出大气进入界面。在该界面可开展大气进入过程的计算。由于小行星2024 RW1的反照率与材质未知,因此,在计算大气进入过程时,假设其反照率为平均反照率0.15,根据反照率、绝对星等、直径之间的经验关系计算小行星直径。在蒙特卡罗仿真中,假设小行星为石质小行星,密度为3500 kg/m3,小行星强度服从0.1~20 MPa的均匀分布。在计算标称解时,尝试根据不同的小行星材质与强度进行分组计算,分组情况如表6所示,对于石质小行星,假设其强度分别为1 MPa、2 MPa、5 MPa、10 MPa、15 MPa;对于碳质小行星,假设其强度为0.1 MPa;对于铁质小行星,假设其强度为20 MPa。分组情况与每组对应的空爆位置与空爆时刻计算结果如表6所示。其中,石质小行星密度为3500 kg/m3,碳质小行星密度为2000 kg/m3,铁质小行星密度为7900 kg/m3。为方便对比,表6中增加了CNEOS记录的2024 RW1的空爆位置与时刻,其数值与表4一致。位置差为计算得到的空爆位置与CNEOS记录下的空爆位置在大地上投影的距离。
表 6 不同强度及材质下标称解计算得到的2024 RW1空爆信息

Table 6 Airburst information of 2024 RW1 calculated from nominal solution with various strengths and materials

参数 空爆时刻(UTC) 空爆位置(经纬高) 位置差/km
石质,强度1 MPa 2024-9-4 16:39:31 17.90$ \text{°} $N,122.84$ \text{°} $E,40.23 km 12.91
石质,强度2 MPa 2024-9-4 16:39:31 17.91$ \text{°} $N,122.85$ \text{°} $E,36.83 km 10.74
石质,强度5 MPa 2024-9-4 16:39:31 17.94$ \text{°} $N,122.87$ \text{°} $E,31.88 km 7.633
石质,强度10 MPa 2024-9-4 16:39:31 17.95$ \text{°} $N,122.89$ \text{°} $E,27.78 km 5.164
石质,强度15 MPa 2024-9-4 16:39:31 17.97$ \text{°} $N,122.90$ \text{°} $E,25.16 km 3.834
碳质,强度0.1 MPa 2024-9-4 16:39:30 17.85$ \text{°} $N,122.79$ \text{°} $E,50.52 km 19.51
铁质,强度20 MPa 2024-9-4 16:39:32 17.97$ \text{°} $N,122.91$ \text{°} $E,23.37 km 2.946
CNEOS 2024-9-4 16:39:32 18$ \text{°} $N,122.9$ \text{°} $E,25 km
对比表6可发现,与CNEOS记录的空爆位置相比,各场景下标称解得到的位置差均在20 km内,高度差在26 km内。其中,“石质、强度15 MPa”的空爆高度与CNEOS记录下2024 RW1的空爆高度最为接近,相差0.16 km,该组的经纬度与CNEOS记录的空爆时刻经纬度之间的位置差为3.834 km;“铁质、强度20 MPa”的空爆位置差最小,为2.946 km,但其高度差较“石质、强度15 MPa”更大,为1.63 km。考虑到近地小行星中石质小行星最为普遍,且该组的高度差最小,因此接下来以“石质、强度15 MPa”组的计算结果进行展示,如图7所示。此外,图6结果面板中“UTC”时刻表明,标称解下小行星进入大气时刻(即约100 km高度处)UTC时间2024年9月4日16时39分27秒,对比表6中各组的空爆时刻可知,不同材质、不同强度场景下,从大气进入至空爆时刻所经历的时间为3~5 s,与CNEOS记录下的空爆时刻相差1~2 s。因此可认为在2024 RW1的场景中,空爆时刻误差相对较小。需要注意的是,表6中分组的目的为展示不同材质、不同强度下标称解计算得到的空爆位置、时刻与CNEOS记录值之间的差异,选择以“石质、强度15 MPa”组的计算结果进行后续展示不代表2024 RW1必然为强度15 MPa的石质小行星。
图 7 2024 RW1大气进入计算结果

Fig.7 Atmospheric entry results of 2024 RW1

图7中,左侧“大气进入参数”面板中“参考点位置”即表4中CNEOS记录的2024 RW1的空爆位置。“标称解(空爆时刻)”面板中为标称解下空爆时刻的经度、纬度、高度、UTC时刻以及与参考点的距离。从标称解空爆位置投影图可发现,空爆位置位于菲律宾东北方向。从强度与距离分布图可知,在蒙特卡罗仿真中,虚拟小行星的强度呈均匀分布,所有虚拟轨道对应的空爆位置与参考点位置之间的距离均在20 km内,绝大多数场景的距离差在10 km内,意味着陨落预报结果具有较高的可靠性。

4 应用案例

以2024 UQ为例,展示SHARP-IIS系统在更短观测弧段的撞击场景中的应用。小行星2024 UQ于2024年10月22日被发现,在被发现后的1 h内撞击地球,是人类成功预警的第10颗近地小行星。与2024 RW1不同,MPC公布的该小行星的观测数据共9条(观测弧段1.53 h)。该小行星进入大气后仍旧被CNEOS监测并记录,空爆时刻的信息如表7所示,CNEOS并未推算出空爆时刻的速度大小及方向。
表 7 CNEOS记录的小行星2024 UQ的空爆信息

Table 7 Airburst information of 2024 UQ recorded by CNEOS

参数数值
空爆时间(UTC)2024-10-22 10:54:48
空爆时刻经度136.0°W
空爆时刻纬度30.0°N
空爆时刻高度38.2 km
总辐射能4.5×1010 J
撞击能量0.15 kt TNT
通过读取MPC记录的观测数据,利用SHARP-IIS系统对2024 UQ的空爆位置进行计算。轨道确定结果如图8图9所示。轨道改进阶段实际采用的数据条数为8条,从图8中小行星位置、速度的协方差矩阵可看出,观测数据的减少使得标称解的不确定性增加,与2024 RW1相比,2024 UQ的协方差矩阵中,位置方差高1个数量级,速度方差高1~3个数量级,这意味着在蒙特卡罗仿真中,小行星的轨道将更加分散,最终得到的空爆位置也将更加分散。
图 8 2024 UQ轨道确定结果-轨道根数

Fig.8 Orbit determination results of 2024 UQ-orbital elements

图 9 2024 UQ轨道确定结果-小行星日心轨道

Fig.9 Orbit determination results of 2024 UQ-heliocentric orbit diagram

表8展示了SHARP-IIS计算得到的轨道根数与ESA计算结果的对比情况。同样地,为方便对比,文中将SHARP-IIS的轨道根数递推至ESA解的历元时刻。可发现,二者的半长轴绝对差约0.001 au,偏心率绝对差异在小数点后四位,轨道倾角与升交点经度的绝对差异在小数点后五位,近日点幅角与平近点角的绝对差异在小数点后两位。两组解的整体差异较小。此外,由于JPL在2024 UQ的轨道确定过程中将CNEOS记录的空爆位置也作为了参考,其计算结果无法代表轨道预报场景的定轨情况,因此本文未与JPL计算的2024 UQ的轨道根数进行对比。
表 8 2024 UQ日心轨道根数对比

Table 8 Comparison of heliocentric orbital elements of 2024 UQ

参数 SHARP-IIS ESA
历元时刻UTC 2024-10-22 8:41:12 2024-10-22 8:41:12
半长轴 $ a $/au 2.19653295665941 2.1954829373825215
偏心率 $ e $ 0.729968829187115 0.72985807323496321
轨道倾角 $ i $/($ \text{°} $) 1.72114897841168 1.7211369336442
升交点经度 $ \mathrm{\Omega } $/($ \text{°} $) 209.138619808433 209.1386348272669
近日点幅角 $ \omega $/($ \text{°} $) 267.609121468536 267.6194525259180
平近点角 $ M $/( $ \text{°} $) 346.198341305274 346.18693332422004
小行星2024 UQ的轨道递推结果如图10~图12所示。从图10可知,标称解下SHARP-IIS系统预测的2024 UQ进入大气时刻为2024年10月22日10时53分23秒,撞击概率为100%。其进入大气时刻速度相对较高,为23.45 km/s;进入大气后预计从西南向东北方向飞行。在蒙特卡罗仿真中,2024 UQ进入大气时刻的经纬度分布范围更广,其经度跨度约0.7$ \text{°} $,纬度跨度约0.1$ \text{°} $,进入角与方向角分布范围约0.8$ \text{°} $,速度分布范围约0.1 km/s,约比2024 RW1对应范围大1个数量级。
图 10 2024 UQ轨道递推结果-轨迹投影

Fig.10 Orbit propagation results of 2024 UQ-trajectory projection

图 11 2024 UQ轨道递推结果-位置速度分布

Fig.11 Orbit propagation results of 2024 UQ-distribution of position and velocity

图 12 2024 UQ轨道递推结果-大气进入参数分布

Fig.12 Orbital propagation results of 2024 UQ-distribution of atmospheric entry parameters

在确定2024 UQ的撞击概率及进入大气时刻的位置、速度后,开展大气进入计算与空爆预报。与2024 RW1场景类似,为确定材质与强度对空爆时刻及位置的影响,对不同材质、不同强度的场景下标称解计算的空爆位置与时刻进行对比,其结果如表9所示。
表 9 不同标称解强度及材质计算得到的2024 UQ空爆信息

Table 9 Airburst information of 2024 UQ calculated from nominal solution with various strengths and materials

参数 空爆时刻(UTC) 空爆位置(经纬高) 位置差/km
石质,强度1 MPa 2024-10-22 10:53:26 29.91$ \text{°} $N,136.16$ \text{°} $W,41.78 km 18.85
石质,强度2 MPa 2024-10-22 10:53:27 29.91$ \text{°} $N,136.14$ \text{°} $W,38.78 km 17.03
石质,强度5 MPa 2024-10-22 10:53:27 29.93$ \text{°} $N,136.11$ \text{°} $W,33.48 km 13.83
石质,强度10 MPa 2024-10-22 10:53:27 29.93$ \text{°} $N,136.09$ \text{°} $W,29.43 km 11.42
石质,强度15 MPa 2024-10-22 10:53:27 29.94$ \text{°} $N,136.08$ \text{°} $W,26.88 km 9.929
碳质,强度0.1 MPa 2024-10-22 10:53:26 29.89$ \text{°} $N,136.23$ \text{°} $W,52.03 km 25.10
铁质,强度20 MPa 2024-10-22 10:53:27 29.94$ \text{°} $N,136.06$ \text{°} $W,25.05 km 8.881
CNEOS 2024-10-22 10:53:48 30$ \text{°} $N,136$ \text{°} $W,38.2 km
表9可知,空爆高度随强度的增加而递减,位置差随强度的增加而减小,这意味着CNEOS记录的参考点位于标称解大气进入轨迹的下游,强度大的小行星其空爆位置更靠近轨迹下游,因此与CNEOS记录的位置更为接近。然而,高强度场景下的空爆高度较低,与CNEOS组的空爆高度差更高,说明标称解下,大气进入的位置与实际情况仍存在误差。表9中除CNEOS外所有场景的时间差为1~2 s,进入大气时刻至空爆时刻所经历的时间段为3~4 s,与CNEOS组的空爆时刻相差21~22 s,说明时间差主要由轨道确定结果的误差引起。ESA预测的2024 UQ撞击时刻为2024-10-22 10:53:29,与本文标称解计算的空爆时刻接近,这说明轨道确定结果误差很可能由观测数据条数过少、观测弧段较短所致。与CNEOS记录的2024 UQ空爆时刻与高度相比,在所设置的场景中,“石质、强度2 MPa”对应的空爆高度与38.2 km最为接近,高度差为0.58 km,此时空爆位置差为17.03 km,时间差为21 s。该组的大气进入结果如图13所示。
图 13 2024 UQ大气进入计算结果

Fig.13 Atmospheric Entry Results of 2024 UQ

图13中,参考点位置为表5所示的CNEOS记录下的2024 UQ的空爆位置。对比计算结果与表8可发现,SHARP-IIS预测的2024 UQ标称解空爆位置误差在20 km内,蒙特卡罗仿真得到的位置误差在40 km内,超过半数场景的位置误差在20 km内。这意味着在观测弧长较短、观测数据较少的撞击场景下,SHARP-IIS的计算结果依然具有较好的参考意义。

5 结语

危地小行星撞击风险预报系统-短临预报场景(SHARP-IIS)实现了从近地小行星观测数据到轨道确定、轨道递推、陨落预报的全流程分析。对于观测弧长7 h、观测数据密度较高的小行星2024 RW1撞击场景,SHARP-IIS对空爆场景的预测实现了较好的精度,在强度15 MPa石质小行星的场景下,标称解计算得到的空爆位置误差在4 km内,时间误差约1 s;在蒙特卡罗仿真中,空爆位置误差在20 km内,验证了本系统的正确性与可行性,在预警的小行星撞击事件中为陨落区域预报、撞击灾害分析提供有力参考。在观测弧长1.5 h、观测数据较少的小行星2024 UQ撞击场景中,SHARP-IIS预测的标称解下石质小行星2 MPa场景下的空爆位置误差在20 km内,蒙特卡罗仿真中的空爆位置误差在40 km内,对于小行星撞击的灾害应对、风险评估等具有良好的参考价值。
目前,除短临预报场景外,“锐警”威胁预警场景模块已开发完毕。在未来,“锐警”系统将进一步在功能上进行完善:系统将在保留手动下载并读取观测数据手段的基础上,增加实时监测功能,从MPC实时获取更新的观测数据并进行分析;同时,系统将通过引入混合解体模型实现多次解体过程的模拟,增加用户选择解体模型选项,计算小行星碎片的落点分布,实现更加全面的陨落预报分析;此外,系统还将进一步优化计算方法,在保证准确率的同时提高计算效率。

1数据获取方式:https://cneos.jpl.nasa.gov/stats/size.html

1
赵坚, 张如生, 李明涛, 等. 近地小行星监测预警六度分析框架[J]. 科学通报, 2023, 68 (8): 981- 992.

2
OH J M, VENTERS C C, DI C, et al. U1 snRNP regulates cancer cell migration and invasion in vitro[J]. Nature Communications, 2020, 11 (1): 1.

DOI

3
CHIARENZA A A, FARNSWORTH A, MANNION P D, et al. Asteroid impact, not volcanism, caused the end-Cretaceous dinosaur extinction[J]. Proceedings of the National Academy of Sciences, 2020, 117 (29): 17084- 17093.

DOI

4
CHYBA C F, THOMAS P J, ZAHNLE K J. The 1908 Tunguska explosion: atmospheric disruption of a stony asteroid[J]. Nature, 1993, 361 (6407): 40- 44.

DOI

5
BRONSHTEN V A. Nature and destruction of the Tunguska cosmical body[J]. Planetary and Space Science, 2000, 48 (9): 855- 870.

DOI

6
POPOVA O P, JENNISKENS P, EMEL'YANENKO V, et al. Chelyabinsk airburst, damage assessment, meteorite recovery, and characterization[J]. Science, 2013, 342 (6162): 1069- 1073.

DOI

7
BOROVIČKA J, SPURNÝ P, BROWN P, et al. The trajectory, structure and origin of the Chelyabinsk asteroidal impactor[J]. Nature, 2013, 503 (7475): 235- 237.

DOI

8
NESVORNÝ D, VOKROUHLICKÝ D, SHELLY F, et al. NEOMOD 3: The debiased size distribution of Near-Earth Objects[J]. Icarus, 2024, 417: 116110.

DOI

9
BROWN P, WIEGERT P, CLARK D, et al. Orbital and physical characteristics of meter-scale impactors from airburst observations[J]. Icarus, 2016, 266: 96- 111.

DOI

10
DEVILLEPOIX H A R, BLAND P A, SANSOM E K, et al. Observation of metre-scale impactors by the Desert Fireball Network[J]. Monthly Notices of the Royal Astronomical Society, 2019, 483 (4): 5166- 5178.

DOI

11
JENNISKENS P, SHADDAD M H, NUMAN D, et al. The impact and recovery of asteroid 2008 TC3[J]. Nature, 2009, 458 (7237): 485- 488.

DOI

12
FARNOCCHIA D, JENNISKENS P, ROBERTSON D K, et al. The impact trajectory of asteroid 2008 TC3[J]. Icarus, 2017, 294: 218- 226.

DOI

13
FARNOCCHIA D, CHESLEY S R, BROWN P G, et al. The trajectory and atmospheric impact of asteroid 2014 AA[J]. Icarus, 2016, 274: 327- 333.

DOI

14
GENG S, ZHOU B, LI M. Near-Earth object 2022 EB5: From atmospheric entry to physical properties and orbit[J]. Astronomy & Astrophysics, 2023, 670: A27.

15
KARETA T, VIDA D, MICHELI M, et al. Telescope-to-Fireball Characterization of Earth Impactor 2022 WJ1[J]. The Planetary Science Journal, 2024, 5 (11): 253.

DOI

16
BISCHOFF A, PATZEK M, DI ROCCO T, et al. Saint-Pierre-le-Viger (L5-6) from asteroid 2023 CX1 recovered in the Normandy, France—220 years after the historic fall of L'Aigle (L6 breccia) in the neighborhood[J]. Meteoritics & Planetary Science, 2023, 58 (10): 1385- 1398.

17
SPURNÝ P, BOROVIČKA J, SHRBENÝ L, et al. Atmospheric entry and fragmentation of the small asteroid 2024 BX1: Bolide trajectory, orbit, dynamics, light curve, and spectrum[J]. Astronomy & Astrophysics, 2024, 686: A67.

18
SEKIGUCHI T. Meteor shower originating from asteroid 2024 RW1 (CAQTDL2) by SonotaCo Network in Japan[J]. eMeteorNews, 2024, 9 (6): 403- 405.

19
TOMMEI G. A new way of thinking about impact monitoring of near-earthobjects[C]//Proc. 1st NEO and Debris Detection Conference,Darmstadt,Germany,22-24 January 2019. 2019.

20
SPOTO F, DEL VIGNA A, MILANI A, et al. Short arc orbit determination and imminent impactors in the Gaia era[J]. Astronomy & Astrophysics, 2018, 614: A27.

21
FRÜHAUF M,MICHELI M,OLIVIERO D,et al. Meerkat Asteroid Guard imminent impactor warning service of the European Space Agency[C]//7th IAA Planetary Defense Conference,2021.

22
ROA J, FARNOCCHIA D, CHESLEY S R. A novel approach to asteroid impact monitoring[J]. The Astronomical Journal, 2021, 162 (6): 277.

DOI

23
FARNOCCHIA D, CHESLEY S R, MICHELI M. Systematic ranging and late warning asteroid impacts[J]. Icarus, 2015, 258: 18- 27.

DOI

24
HU S, LI B, JIANG H, et al. Peculiar orbital characteristics of Earth quasi-satellite 469219 Kamo, oalewa: Implications for the Yarkovsky detection and orbital uncertainty propagation[J]. The Astronomical Journal, 2023, 166 (4): 178.

DOI

25
党雷宁, 柳森, 白智勇, 等. 小行星进入与撞击效应评估模型敏感性研究[J]. 力学学报, 2021, 53 (1): 278- 292.

文章导航

/