小行星进入与撞击效应评估模型敏感性研究1)

党雷宁 柳森 白智勇 石义雷

(中国空气动力研究与发展中心超高速空气动力研究所,四川绵阳 621000)

引言

近地小行星对地球的撞击是人类生存和发展面临的潜在威胁之一.学术界已经基本认同:6500 万年前一颗直径10 km 的小行星撞击现墨西哥希克苏鲁伯地区,导致包括恐龙在内的全球物种大灭绝,史称K-T 事件[1].1908 年6 月30 日,在俄罗斯西伯利亚的Tunguska 发生小行星撞击事件,地面大范围的树木被破坏,释放出约千万吨TNT 当量的能量,受影响面积相当于美国华盛顿市[2].2013 年2 月15 日,在俄罗斯Chelyabinsk,一个尺寸近20 m 的流星在空中爆炸解体,释放出约50 万吨TNT 当量的能量,地面建筑物受损,数千名地面人员受到不同程度的伤害[3].

小行星撞击地球存在复杂的物理化学和力学过程,引起多种类型的进入与撞击效应及灾害[4].小行星在进入地球大气过程中发生空中爆炸,冲击波向地面传播,引起地面超压损伤;爆炸产生的热辐射传至地面,造成热辐射损伤.小行星穿透大气后,对陆地的直接撞击产生撞击坑,对海洋的直接撞击可能引发海啸,直接撞击还可能把部分物质抛向空中,产生反溅碎片云.对于千米量级的小行星,撞击产生的反溅碎片云将会引发全球气候变化.

计算分析小行星进入与撞击效应,有数值仿真和工程计算两种手段.数值仿真方法保真度高、耗时长,适用于具体撞击事件或灾害类型的精细分析.Boslough 和Crawford[5]以及Robertson[6]采用冲击动力学数值方法,对小行星进入过程的解体、空中爆炸以及撞击海洋海啸开展了研究.Aftosmis 等[7]以及Collins 等[8]基于可压缩空气动力学数值方法,研究了空中爆炸产生爆炸波的传播过程,评估了工程方法中常用的静态点爆炸模型的准确性.Johnston 等[9]采用高温气体动力学数值模拟手段,耦合高温烧蚀流场和辐射计算,研究了进入过程产生的热辐射.在国内,张庆明[10]以及赵君尧[11]利用冲击动力学数值软件,初步研究了撞击坑形成过程.李毅等[12]初步开展了空中爆炸与撞击海面的数值仿真.

与数值仿真相比,工程方法采用快速预测的模型,能够在短期内开展大量状态的灾害评估,并获得可接受准度的计算结果.1993 年,Chyba 等[2]提出并采用Pancake 解体模型,研究了Tunguska 小行星的空中爆炸,获得与数值仿真、地面试验一致的空爆高度.同一时期,Hills 和Coda[13]研究了解体、空中爆炸、超压、热辐射、撞击坑、海啸等各种类型的进入与撞击效应模型,特别是基于美国核试验数据给出地面超压半径的拟合公式.2005 年,Collins 等[4]广泛研究了各种类型的进入与撞击效应工程模型,推导出相应的解析公式,研制了网页版的小行星撞击灾害评估软件,放置在互联网上供开放使用.2017 年,Collins 等[8]又对该软件进行了改进.2002 和2017 年,美国NASA的科学定义组[14-15](Science Definition Team)利用当时工程模型发展的成果,考虑小行星物理特性和进入参数的不确定性以及撞击频率,给出小行星撞击年度受影响人口及其概率.美国NASA 还于2015—2017年,组织开发了小行星进入与撞击风险分析评估系统PAIR(Probabilistic Asteroid Impact Risk)[16-17],设定小行星撞击地球的场景,可定量分析进入和撞击过程对地球人口和基础设施的危害,并构成了美国政府、军方、学术界小行星撞击联合应急演习[18-19]的主要内容.在国内,柳森等[20]梳理归纳了小行星撞击地球的超高速问题,并采用文献中已有的模型,建立了小行星进入与撞击效应分析评估软件AICA,利用该软件对2019 年行星防御国际会议假想撞击事件开展分析评估[21].周炳红等[22]对小行星进入过程能量沉积开展了计算分析.

小行星进入与撞击效应评估模型是工程计算的基础与核心.随着人们对小行星进入与撞击过程认识的不断深入,随着与数值仿真最新成果的结合,这些工程模型不断迭代和发展,计算结果的准度不断提高,已经成功地应用于Tunguska,Chelyabinsk,Košice,Benešov 等实际撞击事件的计算分析[2,23-25],也在多次行星防御演习中获得应用[18-19].尽管如此,当前的工程模型输入参数不确定度大[15,17],给模型的使用带来困扰,也较大地影响了危害评估结果.首先,天文观测给出的潜在防御目标小行星的物理特性和运动参数不确定度大,即材料、尺寸、强度与进入参数通常在较大范围.其次,由于基础研究不充分,工程模型中的烧蚀系数、流明效率、扩散系数不确定度大.最后,工程模型特别是某些解体模型存在自由参数.

因此,需要研究模型输入参数变化对输出参数的影响,即开展进入与撞击评估模型敏感性分析.Hills 和Coda[13]研究了不同小行星尺寸、进入速度和材料对超压损伤范围的影响,属于早期(1993 年)的工作,采用的工程模型与目前有所不同.Mathias[26]研究了输入参数对伤亡人口的影响,重点在于受影响人口,空爆模型也与目前不同.Wheeler 等[23]开展了FCM 解体模型的敏感性研究.这些研究还不够充分.着重分析输入参数对能量沉积的影响,缺乏对空爆高度和地面损伤范围的研究.

本文利用课题组建立的小行星进入与撞击效应分析评估软件AICA[20],研究输入参数的变化对能量沉积、空爆高度和地面损伤范围等输出参数的影响,期望得出一些有意义的结论,指导软件的使用和下一步的研究工作,同时提高对小行星进入与撞击效应的认识.

1 小行星进入与撞击效应分析评估软件AICA 介绍

1.1 AICA 软件中的工程模型

目前,AICA 软件包含了进入弹道方程、质量损失方程、解体判据、解体模型、空中爆炸模型,地面超压损伤范围模型、地面热辐射损伤范围模型等小行星进入与撞击效应工程模型,暂未包含撞击坑、海啸和全球气候变化模型.下面逐一阐述这些模型.

1.1.1 弹道方程

在AICA 软件中,考虑到后续研究精确求解弹道的需要,所用弹道方程包含圆球自转地球模型以及经纬度变量.若仅考虑飞行中的阻力,则单个飞行体的弹道方程[29]为

其中,h为飞行高度;θ 和Φ 分别是经度、纬度;V是飞行速度大小;γ 是飞行航迹角或进入角,定义为速度矢量与当地水平面的夹角,速度矢量在当地水平面上方为正;ψ 是航向角,定义为速度矢量在当地水平面的投影与当地纬线方向的夹角;Cd是阻力系数;A是横截面积;rt是小行星到地心的距离;ω 是地球自转角速度;采用美国1976 年标准大气计算大气密度ρair.

在本文的进入与撞击效应工程模型中,有两种类型的单个飞行体:以刚体运动的小行星和解体后的碎片,以变形体运动的解体后的碎片云(见1.1.4节).为了在后续分析中有所区别,本文分别记这两种类型单个飞行体阻力系数为Cd,f和Cd,c.

对于小行星和碎片,阻力系数与外形和飞行姿态相关,横截面积与外形相关.小行星和碎片是非规则外形,由于天文观测方面的基础研究不充分,大多数进入过程研究都假定圆球外形[4,8,13-18,23-28].Wheeler 等[24]在圆球外形的假设下利用改进的FCM 解体模型,获得了与多个小行星进入过程能量沉积观测数据一致的计算结果.Ceplecha 和Revelle[30]、Borovička 和Kalenda[31]、Borovička 等[32]在运动方程中利用参数ΓA′考虑外形和姿态对阻力系数的影响,成功地开展了多个小行星进入过程观测数据分析及弹道重构.其中,Γ 是Ceplecha 定义的阻力系数,等于0.5Cd,f;A′是形状因子,圆球的形状因子为1.21.小行星的气动升力与外形相关,受姿态影响较大.一方面,小行星在进入前往往存在自旋[33],进入大气后自旋有可能继续存在[34];另一方面非规则外形的姿态运动受气动力矩、烧蚀引射和解体爆炸影响,变化规律非常复杂,迄今为止还没有相关研究.当前的研究不考虑升力影响,从研究结果上看与观测数据基本一致[2,24,30,32].

对于以变形体运动的碎片云,根据Pancake 模型[2,13],其外形是压扁的“薄煎饼”形状,横截面积A=πr2(r为半径).由于碎片云不是刚体,在运动中不断变形,仅需要考虑阻力.

Wheeler 等[16-17]在PAIR 软件中以及Collins等[4,8]在其危害评估软件中对这两种类型的单个飞行体,都采用1.0 的阻力系数.白智勇等[35]采用高超声速化学非平衡流动数值模拟手段,获得了直径1 m的圆球小行星在速度10～20 km/s 的阻力系数,范围为0.87～0.92.Ceplecha 和Borovička[30-32]在多个进入过程观测数据分析及弹道重构中,采用的小行星及碎片ΓA′范围为0.7～1.14,若以圆球形状因子换算,则Cd,f(等于2Γ)的范围为1.15～1.88.

基于以上分析,本文假定小行星及碎片是圆球外形,阻力系数范围为0.8～2.0.通过考察阻力系数对进入与撞击效应的影响,间接研究外形和姿态的影响.对于以变形体飞行的碎片云,由于缺乏研究支撑,无法获得阻力系数取值范围,阻力系数取Wheeler 和Collins 研究中的1.0.

1.1.2 质量损失方程

小行星由于烧蚀而不断损失质量,可以直接使用文献中广泛应用的质量损失方程[2,4,23-28]

其中,m是小行星瞬时质量,σ 是烧蚀系数,CH是热流系数,Q是烧蚀热.CH是小行星尺寸、飞行速度、高度的函数,Q与小行星的材料、结构相关.基于当前学术界的认识[2,17,36-37],烧蚀系数σ 取值范围为:3.5×10−10～7.0×10−8kg/J.

1.1.3 解体判据

在天文学中,广泛使用如下解体判据[2,4,17,36]

其中,P是小行星驻点压力,S是解体强度.即当驻点压力超过强度,小行星解体.

小行星与陨石的强度不同.陨石经历了严酷的气动力、热而降落至地面,代表了小行星结构中强度最大的部分.而小行星通常是包含裂纹、空穴的碎石堆结构,强度要小得多.根据Popova 等[38](2011)对大量观测结果的总结,0.1～106kg 质量范围的石质小行星的初次解体强度为0.1～10 MPa,而材料的平均拉伸强度为30 MPa,压缩和屈服强度更大.

Weibull 和Sweden[39]统计分析了地球上的岩石强度,认为强度和质量服从如下统计规律

其中,Sc是岩石质量mc的强度,Sp是同一岩石材料质量mp的强度;α 是强度指数,0 ＜α ＜1,不同材料有不同的取值.Cotto-Figueroa 等[40]通过对不同尺寸陨石开展强度试验,研究了两种陨石材料强度的统计散布,并与小行星进入观测得到的解体强度对比,初步认为小行星材料也服从Weibull 定律.

一些解体模型应用Weibull 定律计算解体后的碎片强度[8,23],此时下标“c”和“p”分别表示解体后的子碎片以及解体前的父碎片.Wheeler 等[23]在FCM解体模型的敏感性分析中,认为α 在0.1～0.3 范围内取值可以得到合理的、能够代表小行星结构变化对解体影响的能量沉积结果,在其利用PAIR 软件开展撞击风险分析中也取这个范围[17].NASA 科学定义组在2017 年的报告中[15],对含水和不含水的石质小行星分别采用0.1 及0.2 的取值,并认为铁质小行星不解体.

基于此,本文在敏感性分析中,认为强度指数α的取值范围为0.1～0.3.

1.1.4 解体模型

解体后,需要应用解体模型确定小行星的存在形式.目前存在多种解体模型,包括连续的Pancake模型[2,13]和Park 模型[41]、离散的“collective wake”模型[42]和“independent wake”模型[43],以及把离散和连续模型综合起来的耦合模型[44]及其改进后的FCM模型[23-24].

FCM 模型认为,小行星每次解体为一系列大尺度的碎片和一个由小碎片、液滴组成的碎片云.所有的碎片都是独立的飞行体,其运动服从方程组(1)和(2).碎片云的运动按照Pancake 模型计算,其整体运动在服从方程组(1)和(2)的同时,横截面积还在前后压力差的作用下逐渐增大.Hills 和Coda[13]根据能量守恒关系,推导出碎片云变形中横向扩散速度公式

其中,Vdisp是碎片云的扩散速度,V是飞行速度,ρair和ρ 分别是大气及小行星密度.Cdisp是扩散系数,Hills 和Coda[13]认为其量值为3.5,Passey 和Melosh略有不同,认为其量值为3.0.Wheeler 等[23]在FCM模型的敏感性分析中,认为Cdisp的变化对能量沉积影响很小.因此在本文中,不研究Cdisp的变化对模型输出的影响.

1.1.5 空中爆炸模型

在进入大气过程中,小行星可发生多次空中爆炸.目前的空中爆炸模型尚不能处理多次空爆情形,而认为空中爆炸是发生在某高度的点爆炸[2,4,7,13-17],称为静态点爆炸模型.点爆炸的高度称为空爆高度,点爆炸的能量称为空爆能量.

Chyba 等[2]以及Collins 等[4]基于Pancake 模型,认为当碎片云半径扩散到小行星初始半径的某倍数时,发生空中爆炸.NASA Ames 在建立PAIR 软件之初,认为空爆高度是瞬时动能减小到初始进入动能一半时的高度[26].NASA 科学定义组2003 和2017 年的报告[14-15]以及在PAIR 软件目前的版本[17]中,都认为空中爆炸发生在能量沉积峰值处.

Hills 和Coda[13]在1993 年认为空爆能量是全部动能,1998 年又认为空爆能量是全部动能的一部分[36],Toon 等[45]认为这个比例是50%.Boslough和Crawford[5]发现冲击动力学数值模拟预测的地面超压大于点爆炸模型预测结果,而Collins 等[8]却认为在距离较远处两种模型预测的地面超压基本一致.Wheeler 等[17]基于上述Boslough 的结果,认为把空爆能量定义为全部动能可以更准确地预测地面超压,并将其应用在PAIR 软件中.

基于目前的研究现状,可以认为空中爆炸是发生在能量沉积峰值处以全部初始动能为能量的点爆炸.

1.1.6 地面超压损伤模型

1993 年,Hills 和Coda[13]基于空中爆炸的静态点爆炸模型以及美国核试验数据,拟合出4 psi 地面超压损伤半径的经验公式,Wheeler 等[17]将其用于PAIR 软件中.该公式为

其中,rp为地面超压超过4 psi 的半径,单位是km;h是空爆高度,单位是km;E是空爆能量,单位为MT.MT 是能量单位,指的是106吨TNT 当量的能量,1 MT=4.184×1015J.4 psi 的超压可以使木框架房屋倒塌[4].

2005 年,Collins 等[4]对美国核试验数据进行拟合,推导出另一套地面超压的经验公式,并于2017 年进行了改进[8].Collins 的经验公式把自空爆点正下方(文献中称为ground zero) 冲击波的影响范围分为正则反射与马赫反射两个区域,并分别进行计算.在正则反射区域中

其中,p是地面超压,单位是Pa;r1和h1分别是1 KT爆炸当量时距空爆点正下方的水平距离及空爆高度,可通过缩尺率计算

式(9)中的1000 是把km 转换为m、MT 转换为kT所需的因子.在马赫反射区域中

正则反射区域与马赫反射区域的分界点为

Aftosmis 等[7]通过数值模拟,考察了在推导式(7)～式(11) 中使用的爆炸当量与距离的缩尺率关系,认为浮力破坏了缩尺率关系,在爆炸能量小于数十MT 时,可以使用缩尺率,而当爆炸能量大于200 MT 时,必须考虑浮力的影响.NASA 科学定义组在2017 年的报告中[15],把美国核试验数据与Aftosmis 的数值模拟结果综合起来,评估空中爆炸引起的地面超压.

在AICA 软件中,暂未包含数值模拟结果,主要使用Hills 和Coda 的经验公式(式(7))计算地面超压损伤范围;而在爆炸能量较低、超压小于1 psi 时,使用Collins 的经验公式(式(8)～式(11)) 计算超压损伤半径.图1 针对1KT 空爆能量下4 psi 超压损伤半径与空爆高度的关系,对比了Hills 和Coda 经验公式、Collins 经验公式的计算结果与美国核试验数据(图中标注为Glasstone&Dolan).可以看到,两种经验公式都能获得与美国核试验数据基本一致的结果.

图1 Hills 和Coda 的经验公式与Collins 经验公式计算结果的对比(1 kT 空爆能量,4 Psi 超压)Fig.1 Comparison of Hill&Coda’s and Collins’empirical formula for 1 kT airburst energy and 4 psi overpressure

1.1.7 地面热辐射损伤模型

Collins 等[4]基于美国核试验数据和理论分析,推导出在地面空爆时火球热辐射损伤范围的经验公式,Wheeler 等[17]将其修改到适用于在某高度的空爆情形

其中,rr是地面热辐射损伤半径,单位km;h是空爆高度,单位km;E是空爆能量,单位为MT.η 是流明效率,是辐射能占初始动能的比例.Collins[4]认为η不确定度较大,取值范围为:1.0×10−4～1.0×10−2.Φi是热暴露,即单位面积的传热.Zi是1MT TNT 当量的热暴露,对于3 级烧伤,Zi=0.42 MJ/(m2MT1/6).3级烧伤所需能量大于把草点燃的能量[4].AICA 软件也包含了1 级、2 级烧伤的输出.本文中若无特别说明,热辐射损伤半径指3 级烧伤的损伤半径.

1.2 AICA 软件的计算流程

AICA 软件的计算流程如图2 所示.该软件用弹道计算把各工程模型串联起来.通过数值求解弹道方程(1)和质量损失方程(2)获得单体运动参数和能量沉积.在弹道计算中通过式(4) 判断是否解体,若解体则通过FCM 模型得到解体后碎片及碎片云的数量和质量分数,通过式(5)计算解体后子碎片的强度,通过式(6) 计算碎片云的扩散.当全部碎片和碎片云都完全烧蚀或落地,则计算整个进入过程的能量沉积,并根据空中爆炸模型获得空爆高度和空爆能量.最后基于空爆高度和空爆能量,按照式(7)～式(13)计算地面超压损伤半径和热辐射损伤半径.

1.3 AICA 软件验证

作者2018 年的论文[20],通过Chelyabinsk 小行星进入和Tunguska 小行星进入两个算例,初步验证了AICA 软件.本文在此基础上补充危害评估结果,进一步对AICA 软件进行验证.

图3 是Chelyabinsk 事件的计算结果,其中进入弹道(图3(a)) 和能量沉积结果(图3(b)) 来自作者2018 年论文.计算条件见作者2018 年的论文[20],其中小行星物理特性及进入条件来自观测结果,小行星强度、强度指数及碎片云质量分数来自Wheeler 等[17]的研究结果.可以看到进入轨迹和能量沉积曲线都与观测符合较好;能量沉积曲线能够复现 Chelyabinsk 小行星两次空爆的特征.在Chelyabinsk 小行星进入后,当地收到建筑物玻璃破坏的报告,图3(c) 列出玻璃无损、受损、严重受损以及收集到陨石的经纬度位置[3],根据Popova 等[3]的研究:500 Pa 超压的损伤半径可用来衡量玻璃破坏的最大范围,1000 Pa 超压的损伤半径可用来衡量玻璃严重破坏的范围.本文500 Pa 和1000 Pa 的超压损伤范围计算结果与观测结果基本一致.

图2 AICA 软件流程图Fig.2 Flow chart of AICA code

图3 Chelyabinsk 流星事件计算结果Fig.3 Computational results of Chelyabinsk event

图4 给出Tunguska 大爆炸事件的计算结果,其中能量沉积结果(图4(a))来自作者2018 年论文.计算条件见作者2018 年的论文[20],其中小行星种类和初始进入参数采用Chyba 等的研究成果,解体强度、烧蚀系数和流明效率在石质小行星的范围内取值,强度指数、每次解体后生成的子碎片(云)数及其质量分数通过假设给定.由图4(a)可见,空爆高度为9.2 km,与文献中的8.5 km[2]比较接近.4 psi 的超压可以使木框架房屋倒塌,Wheeler 等[25]也把4 psi 作为研究Tunguska 森林倒伏的超压阈值.由图4(b)可见,4 psi 超压损伤范围能够反映森林倒伏范围,比森林倒伏范围略小.当森林树木受到热辐射着火后,火势会蔓延,导致观测到的森林着火范围(图中标注为forest fire,observed)比最初的热辐射烧伤范围(图中标注为radiation burn,observed)大一些[9].2 级烧伤相当于点燃阔叶林[4],其损伤半径的计算结果与观测的热辐射烧伤范围相当.本文的计算与观测结果还存在一定偏差,原因如下:一是Tunguska 当地森林树木倒伏与超压的定量关系、当地森林树木着火与热辐射的定量关系还在继续研究中;二是地形变化可能对地面损伤有影响.

图4 Tunguska 大爆炸事件计算结果Fig.4 Computational results of Tunguska event

本节以Chelyabinsk 与Tunguska 两次小行星进入事件为算例,说明在给定小行星物理特性、进入参数以及仔细选择解体模型参数的情况下,AICA 软件的计算结果是可靠的.

2 工程模型的敏感性分析

2.1 输入参数范围

输入参数的范围和基准值如表1 所示.本文在基准值的基础上,通过改变一个或多个输入参数,研究进入与撞击效应工程模型的输出对输入参数的敏感性.其中重点研究的输出参数为空爆高度hb,4 psi超压损伤半径rp和3 级烧伤半径rr,也有对弹道和能量沉积的研究.

如表1 所示,小行星直径d0变化范围下限为60 m,旨在包含Tunguska 小行星直径(目前学术界认为是50～80 m[25]);而20 m 直径的Chelyabinsk 小行星不足以产生严重的地面超压和热辐射损伤,在表1 中没有包含.直径上限取300 m,不仅大于潜在威胁小行星的最小直径140 m,而且使得地面灾害以超压和热辐射为主,不用考虑海啸和全球气候变化等灾害[15].表1 中进入速度V0的范围取自Greenstreet等[46](2012 年)的数值模拟结果,忽略概率低的高于40 km/s 的值,基准值20 km/s 是行星防御研究的典型速度.进入角γ0的基准值−45◦是进入角范围中概率最大的取值[17].烧蚀系数σ、初次解体强度S0、强度因子α、流明效率η 的取值范围在1.1 节中已经叙述,基准值取自范围中间值或对数中间值.小行星及碎片的阻力系数Cd,f的取值范围在1.1 节中已经叙述,基准值为Wheeler 和Collins 研究中的取值1.0.

近地小行星绝大多数都是石质小行星,密度不仅与具体组成(含水、不含水)、陨落概率以及孔隙率有关,也和强度相关[15,17].本文不研究密度对进入与撞击评估模型输出的影响,而根据Mathias 等[17]的结果把基准值取为石质小行星中密度概率最大的值.在NASA 科学定义组2017 年的报告中,直接让碎片云质量分数等于50%来评估地面灾害.Wheeler 等在用PAIR 软件计算Chelyabinsk 流星事件的能量沉积中采用了大于80%的碎片云质量分数,在研究Tunguska 小行星可能物理特性时采用了80%的取值.在本文中,取80%作为碎片云质量分数的基准值,2.7 节中将会再次讨论这一取值的合理性.

表1 进入与撞击效应模型敏感性研究输入参数基准值及变化范围Table 1 Baseline value and variation range of AICA input parameters for sensitivity study

2.2 空爆高度、能量对地面损伤半径的影响

在本文的模型中,地面损伤范围与空爆高度、能量直接相关.为此,本文先讨论空爆高度、能量与地面损伤半径的关系,为下面的分析打下基础.

利用Hills 和Coda 的经验公式(式(7)),对表1 得到的空爆能量范围,计算4 psi 超压损伤半径随空爆高度的变化曲线,并如图5 所示,图中各条曲线表示不同的进入动能.正如Hills 和Coda[13]的研究并分析图5,超压损伤半径随空爆高度呈非线性变化;对于某空爆能量,存在使超压损伤范围最大的“最优高度”,小于此高度,空爆高度随超压损伤半径递增,大于此高度,则递减;若空爆高度不变,则超压损伤半径随空爆能量递增.

图5 不同空爆能量下4 psi 超压损伤半径随空爆高度变化曲线Fig.5 Radius of 4 psi overpressure as a function of airburst altitude for burst energy levels computed from Table 1

热辐射损伤半径随空爆高度、能量的变化则简单得多.分析式(12)～式(13),空爆高度的减小、空爆能量的增大,都使得热辐射损伤半径增大.

2.3 直径、进入速度和进入角的影响

图6 给出在小行星直径和进入速度都变化的情况下,空爆高度、4 psi 超压损伤范围、3 级烧伤范围的分布云图,图中的黑色等值线表示空爆能量,单位为MT.在图6 的计算中,除直径和进入速度外的其它输入参数取表1 中的基准值.根据图6(a),大尺寸、小进入速度的小行星具有较低的空爆高度.能量沉积的主要影响机理是解体爆炸.在其他参数都相同的情况下,小行星尺寸越大,解体后的碎片尺寸就越大,就越有能力在较低的高度存在飞行且继续解体的碎片,并对低高度的能量沉积产生贡献.进入速度越小,就需要在较低的高度获得足够的动压以解体,并在低高度沉积能量.

如前所述,地面损伤范围是空爆高度和能量的函数.在图6(a) 中,空爆高度总体上变化不大(0～20 km).而进入动能随直径的立方和速度的平方变化,对超压损伤半径的影响更大.因此,超压损伤半径总体上随小行星直径和进入速度,也即随着空爆能量的增大而增大(图6(b)).同理,3 级烧伤半径主要受空爆能量影响,总体上随着直径和进入速度的增大而增大(图6(c)).

图6 小行星直径和进入速度对空爆高度、4 psi 超压损伤半径以及3级烧伤半径的影响(图中黑色直线等值线为空爆能量,单位MT)Fig.6 Effects of pre-entry diameter and entry velocity on airburst altitude,radius of 4psi overpressure,and radius of 3rd degree burn.The black contour lines indicate burst energy levels with a unit of magaton(MT)

图7 小行星进入角和进入速度对空爆高度、4 psi 超压损伤半径以及3 级烧伤半径的影响Fig.7 Effects of entry angle and velocity on airburst altitude,radius of 4 psi overpressure and 3rd degree burn

图7 给出不同进入速度下进入角对空爆高度和地面损伤半径的影响,图中其他输入变量取表1 中的基准值.在小进入角下(大和小指绝对值,后文相同),小行星的水平速度可能大于第一宇宙速度,将会飞离地球而不会进入.图7(a) 中的灰色区域给出速度12～40 km/s、其他参数取基准值时飞离地球的进入角范围,约为−10◦～0◦.随着进入角的增大,小行星穿透大气层的时间更短,更有可能在较低的高度沉积动能,导致空爆高度降低.进入角对地面损伤半径的影响可从空爆高度的变化解释.在较大的进入角下,空爆高度较低且变化小(图7(a)),对损伤半径的影响小(图7(b)).当进入角较小时,空爆高度可能增大到“最优高度”以上,在进入角向飞离地球的临界值减小的过程中,空爆高度增大,使得地面超压损伤半径减小,也使得热辐射损伤半径减小.从图7(b)还可看出,相同计算条件下4 psi 地面超压半径要大于3 级烧伤半径,这与Wheeler 等[17]的结果相同.

2.4 烧蚀系数、初次解体强度和强度指数的影响

通过同时改变烧蚀系数、初次解体强度和强度指数,其他输入参数取表1 中基准值,计算得到如图8 所示的对空爆高度、4 psi 超压损伤半径、3 级烧伤损伤半径的影响曲线,图8(a) 和图8(b) 中每条曲线代表不同的烧蚀系数,图8(c) 中每条曲线代表烧蚀系数与地面损伤类型的组合.

图8 烧蚀系数、初次解体强度和强度指数对空爆高度、4 psi 超压损伤半径以及3 级烧伤半径的影响Fig.8 Effects of ablation coefficient,strength at 1st breakup and strength scaling exponent on airburst altitude,radius of 4 psi overpressure,and radius of 3rd degree burn

由图8(a) 和图8(b) 可见,烧蚀系数和初次解体强度对空爆高度影响较大,强度指数对空爆高度的影响较小.初次解体强度越大,小行星就越趋向于在较低的高度获得足够的动压、初次解体并沉积能量,使得空爆高度减小.烧蚀系数越大,在下降相同高度的过程中碎片和碎片云就能沉积更多动能,使得在以后的飞行中需要沉积的动能减小,导致能量沉积峰值高度即空爆高度上移.

在图8 的计算条件下,空爆高度总体上小于“最优高度”,与地面超压损伤半径呈递增关系(图5),与热辐射损伤半径呈递减关系,使得随着初次解体强度的增大以及烧蚀系数的减小,4 psi 超压损伤半径减小,3 级烧伤半径增大(图8(c)).总体上看,在图8的计算条件下,地面损伤范围变化较小.

2.5 碎片云质量分数和解体后子碎片数的影响

碎片云质量分数是影响能量沉积、空爆高度和地面损伤范围的关键参数之一.

图9(a)给出不同碎片云质量分数下进入过程的能量沉积情况.可以看到,能量沉积曲线不仅在整体上呈现出大的“鼓包”,而且在局部有小峰值.曲线的这种形态主要由小行星解体行为决定.若仅有一次解体,如Pancake 模型计算结果[2,13](相当于Ccloud=100%),或者当碎片云质量分数较大时(如70%,80%),峰值由初次解体决定,只存在一个大的“鼓包”.若碎片云质量分数较小(如30%～50%),则大质量碎片的再次解体在能量沉积中占主导作用,会呈现数个局部峰值.

由于本文把空爆高度取为能量沉积峰值所在的高度,碎片云质量分数较低时的局部峰值现象,将使空爆高度结果呈现较大波动,并影响地面损伤范围.图9(b)、图9(c) 分别展示了空爆高度、地面损伤范围随碎片云质量分数的变化,图中也给出解体后子碎片数n取2,4,6 的对比.可以看出,空爆高度总体上随着碎片云质量分数的增大而减小.当碎片云质量分数小于50%,空爆高度和地面损伤范围曲线存在明显波动,且空爆高度量值较大;而当碎片云质量分数大于50%时,空爆高度和地面损伤范围曲线则平滑变化,空爆高度量值小得多.

图9 碎片云质量分数和解体后子碎片数对能量沉积、空爆高度和地面损伤范围的影响(4 Psi 超压半径,3 级烧伤半径)Fig.9 Effects of mass fraction of debris cloud and number of son fragments at each breakup on airburst altitude,radius of 4-psi overpressure,and radius of 3rd degree burn

当解体后的子碎片数n从2 增大到6,图9(b)中的各条空爆高度曲线逐渐靠拢,n取4 和6 时的曲线已经非常接近,说明本文采用n=4 的基准值可以得到n取更大值时计算结果.还可看到,当碎片云质量分数大于70%时,解体后的子碎片数对空爆高度和地面损伤范围影响较小.这是因为,在碎片云质量分数较大时,碎片云的减速和烧蚀主导了能量沉积过程.

2.6 流明效率和阻力系数的影响

在本文的进入与撞击效应评估模型中,流明效率η 仅用于计算地面热辐射损伤半径(式(12)),因此仅影响热辐射损伤范围,而不影响超压损伤范围.

在前文的分析中,发现相同的计算条件下4 psi超压损伤半径总大于3 级烧伤半径.那么在流明效率具有表1 所示的不确定度情况下,这个结论是否仍然成立？为此,计算了流明效率取表1 中最小、最大值时的3 级烧伤半径,并与相同条件下的4 psi 超压损伤半径进行比较,如图10 所示.图10 的横轴为进入速度,旨在通过速度变化说明不同空爆能量的影响.可见,3 级烧伤半径在大多数情况下要小于4 psi超压损伤半径,而在进入速度或进入能量较大时,由流明效率最大值得到的3 级烧伤半径略大.Aftosmis等[7]指出,当进入能量大于数百MT 时,公式(7)预测的地面超压范围小于精确的数值模拟结果.因此,可以认为在考虑流明效率不确定度的条件下,4 psi超压损伤半径大于3 级烧伤半径.Tunguska 流星事件中树木倒伏范围大于着火范围(见图4(b) 或文献[47]),也是旁证.

图10 流明效率对3 级烧伤半径的影响Fig.10 Effect of luminous efficiency on radius of 3rd degree burn

图11 给出不同小行星及碎片阻力系数下的能量沉积曲线,图中其他输入参数来自表1 中的基准值.在高度100 km 到初次解体位置43 km 之间,小行星的气动阻力是产生能量沉积的重要原因,阻力系数越大,能量沉积越大,但影响总体较小.小行星解体后,在当前较大的碎片云质量分数基准值下,对能量沉积起主要作用是碎片云的减速和烧蚀,碎片的阻力系数对能量沉积影响较小,也对空爆高度和地面损伤范围影响较小.当然,在小的碎片云质量分数情况下,碎片的阻力系数可能对能量沉积有较大影响,但是正如2.7 节的讨论,在本文研究的小行星直径范围内,80%的碎片云质量分数基准值有较大的合理性.此外,在图11 的计算条件下,碎片云的减速对能量沉积贡献较大.由于缺乏碎片云阻力系数基础研究支撑,无法获得合理的阻力系数变化范围,本文没有研究工程模型输出对碎片云阻力系数的敏感性.

图11 小行星和碎片的阻力系数对能量沉积的影响Fig.11 Effect of drag coefficient of asteroid and fragment on energy deposition

2.7 讨论

本章通过同时改变1～3 个表1 中的输入参数,研究进入与撞击效应评估模型对输入参数敏感性.从研究方法上看,没有覆盖更多变量的同时变化.然而,正如前文对表1 的说明,输入参数的基准值代表了取值范围中概率最大或文献研究的推荐值,因此分析结果仍然具有一定的意义.

由2.3～2.6 节的分析可见,除碎片云质量分数之外的其他输入参数,对空爆高度和地面损伤范围的影响比较“平滑”,没有出现“跳跃”的变化.可以简单推论,在这些输入参数基准值下仅同时改变少数变量得到的定性分析结论,与更多输入变量同时变化的情形相比,不会有质的变化.

在碎片云质量分数小于50%时,空爆高度和地面损伤范围随碎片云质量分数的变化趋势出现了波动.本文没有研究碎片云质量分数小于50%时,小行星直径、进入速度、进入角、初次解体强度和强度指数的影响.Wheeler 等[24]在利用改进的FCM 模型研究Košice,Benešov,Tagish Lake 等小能量进入事件能量沉积的过程中,碎片云质量分数总体大于50%;在进入能量更大的Chelyabinsk 事件能量沉积计算中,碎片云质量分数大于80%[23].而对于近代以来进入能量最大的Tunguska 事件,Chyba 等[2]在计算时采用了Pancake 模型,意味着碎片云质量分数为100%.Wheeler 等和Chyba 等的能量沉积计算结果都与观测或数值模拟符合较好.因此,在本文研究的小行星尺寸条件下(下限是Tunguska 小行星),有理由认为采用80%碎片云质量分数基准值获得的结果,能够代表大多数进入情况.

3 结束语

小行星进入与撞击效应评估模型输入参数不确定度大,给模型的使用带来困扰,也较大地影响了危害评估结果.本文首先阐述了模型的由来和发展现状.然后分析了输入参数的取值范围,并以范围中概率最大或文献的推荐值作为敏感性研究的基准值.最后利用小行星进入与撞击效应分析评估软件AICA,在输入参数基准值的基础上,通过改变一个或多个输入参数,研究能量沉积、空爆高度、4 psi 超压损伤半径和3 级烧伤半径等模型输出参数对输入参数的敏感性.研究的模型输入参数包括:小行星直径,进入速度,进入角,烧蚀系数,初次解体强度,强度指数,碎片云质量分数,解体后的子碎片数,流明效率,小行星和碎片阻力系数.

对于直径60～300 m、进入速度12～40 km/s 的石质小行星,本文的研究指出:

(1)大尺寸、小进入速度的小行星具有较低的空爆高度.超压损伤半径和热辐射损伤半径总体上随着小行星尺寸和进入速度的增大而增大.

(2)在小进入角下,小行星的水平速度可能大于第一宇宙速度,将会飞离地球而不会进入.

(3) 随着初次解体强度的增大或烧蚀系数的减小,空爆高度减小.小行星及碎片阻力系数、强度指数对空爆高度和地面损伤范围影响较小.

(4)碎片云质量分数是影响能量沉积、空爆高度和地面损伤半径的关键参数之一.若碎片云质量分数小于50%,能量沉积曲线出现数个局部峰值,使得空爆高度和地面损伤半径评估结果出现波动.空爆高度总体上随着碎片云质量分数的增大而减小.通过分析Wheeler,Chyba 等对若干小行星进入事件能量沉积的计算结果,可以认为在本文研究的小行星直径范围内,较大的碎片云质量分数基准值(如80%)能够代表大多数进入情况.

(5)流明效率仅影响热辐射损伤范围,不影响超压损伤范围.在本文的计算条件下,3 级烧伤半径小于4 psi 超压损伤半径,小行星进入与撞击危害以空中爆炸产生的地面超压为主.

小行星撞击地球危害评估是行星防御研究的重要方向之一,当前国内的研究处于起步阶段,国外的研究中还有大量问题没有解决.通过本文的敏感性分析,可以认为在后续的工作中,一方面需要在天文和进入过程观测中,建立小行星尺寸、材料、进入速度、进入角、强度等物理和力学特性的高保真度模型;另一方面,需要在进入过程研究中,着重研究烧蚀、空中爆炸、冲击波效应和热辐射问题,持续改进工程模型并减小相关参数的不确定度.