近地小行星2016HO3表面温度建模研究

2020-01-19贾晓宇

深空探测学报 2019年5期

贾晓宇，杨晨，王彤，文毅

（1.北京空间飞行器总体设计部，北京 100094；2.中国空间技术研究院，北京 100094）

引言

2019年4月18日，中国国家航天局（China National Space Administration，CNSA）公布了小行星探测计划：其中探测器将开展对近地小行星2016HO3 的绕飞探测，随后择机附着在小行星表面并采集样品，随后携带样品返回地球，并通过返回器将样品送到地面。

小行星表面温度信息对于小行星探测载荷设计，特别是有近地表操作的原位探测、样品采集的任务十分必要。相比于采用一触即走（Touch ＆ Go）采样方式的“隼鸟号”（Hayabusa）、“隼鸟2号”（Hayabusa-2）[1]与“欧西里斯-雷克斯号”（OSIRIS-Rex）[2]，我国小行星探测将附着在小行星表面，采集样品，小行星表面热环境更为重要。

小行星表面热辐射所产生的净力与力矩会导致亚尔科夫斯基效应与YORP 效应，影响小行星轨道演化，对于足够小的天体，这两种效应更加显著[3-4]，准确分析亚尔科夫斯基效应与YORP效应需要小行星温度分布信息。温度及其在小行星整个生命周期中的演变可以改变小行星的表面组成和性质[5-6]。因此，小行星表面温度信息对于研究小行星轨道与表面物质演化具有重要意义。

本文将结合小行星2016HO3 的特征以及地面初步的观测数据，开展小行星2016HO3 温度分析，获得其上下限，并对其热特征进行初步总结。

1 目标小行星2016HO3基本信息

小行星2016HO3（编号469219，名称为komo’oalewa）是由位于夏威夷海勒卡拉山1.8 m 口径的“泛星1 号”小行星巡天望远镜2016年4月27日发现的。2016HO3 是一颗阿波罗型小行星。图1展示了2016HO3 相对地球运动轨迹，2016HO3 与地球以相同的周期共同绕太阳公转，以地球角度看，2016HO3也在绕地球转动，此类天体被称为地球准卫星（Earth Quasi-Satellite）。2016HO3 是目前发现的最稳定的地球准卫星，至少未来100年内不会飞离地球[7]。

目前地面观测未能直接得到2016HO3 的类型与大小，初步估计2016HO3可能是S/Q/L类小行星，其绝对星等Hv=-24.3±0.5[8-10]，反照率α=0.2[8]（但该数据并非观测所得，而是对S 类小行星的估计值），利用公式

可以计算出2016HO3 的直径为41 m[10]，而JPL的Agle 根据雷达探测数据估计2016HO3 的直径为40～100 m。本文在后续分析中将认为2016HO3的直径是40～100 m。

图1 2016HO3的轨道示意图［7］Fig.1 Orbit of 2016HO3［7］

表1 2016HO3的轨道参数［8］Tabel1 Orbital parameters of 2016HO3［8］

图2 2016HO3的光变曲线［11］Fig.2 Light curve of 2016HO3［11］

根据光变曲线（参见图2），目前可以得到2016HO3 较为准确的自转周期为0.467 ± 0.008 hr（28΄1.2" ± 28.8"）[10]。目前无法确定2016HO3 自转轴的方向，甚至无法判定2016HO3 是否存在稳定的自转轴。Burns （1973）曾经提出一个公式用于计算小天体自转轴趋于稳定所需时间尺度τdamp[12]

其中：Q代表品质因子，即每圈转动所损失的能量；μ代表物质硬度；ρ是小天体的密度；ω是小天体自转角速度；R是小天体半径；K32是一个与天体形状相关的无量纲量，对于球状小天体，K32接近于0.01，对于长椭球状小天体，K32接近于0.1。

Harris （1994）给出一个更简化的公式计算小天体自转轴趋于稳定所需时间τdamp[13]

其中：P是小天体自转周期，单位为h；D是小天体直径，单位为km，C是一个无量纲量，取17±2.5；τdamp的单位是 109年。则可以计算出2016HO3 的τdamp为2×106～1.3×107年，如果2016HO3的自转轴不稳定，则意味着其自转轴需要数百万年甚至更长时间才能稳定，目前尚不清楚2016HO3 演化历史，2016HO3的自转轴是否稳定无法判断。

小行星的热惯量是计算表面温度分布的关键参数，但是目前缺乏2016HO3 热红外观测数据对其温度、热惯量进行直接分析。

热惯量与小行星表面的颗粒大小有关，小行星表面物质颗粒越大，热惯量也越大，如月球表面等土壤颗粒非常细小，其热惯量仅为40～50 J·m-2·s-0.5·K-1，小行星数量巨大，表面形态不同，热惯量的跨度较大，根据目前公布的数据，最小的可以到1 J·m-2·-0.5（10199 Chariklo），大的可以到数百[14]。尽管根据“欧西里斯号”对贝努最新观测表明表面颗粒直径并非决定小行星表面热惯量唯一决定因素，但表面颗粒直径越大，小行星热惯量越大的结论仍然能够成立[15]。同时，小天体直径越小，其表面引力越微弱，则微小的风化层颗粒越难以保存在其表面上，因此小行星直径越大，其表面的热惯量也就越大。

图3 小行星的直径与热惯量关系[14]Fig.3 The relationship between the diameter and the thermal inertia of an asteroid[14]

Delbo 的研究给出了一个小行星直径（D）与热惯量（Γ）之间近似的拟合关系

其中：d0=300 ± 47 J·m-2·s-0.5·K-1，ξ=0.48±0.04[14]。

2016HO3的直径约在40～100 m之间，可以得到2016HO3 的热惯量为 700～1 850 J·m-2·s-0.5·K-1.。需要注意，该数据来源于公式（4）的外推，对于该数值应采取谨慎的态度。

图4 高分辨率贝努表面照片［18］Fig.4 High resolution photo for the surface of Bennu［18］

根据近期“欧西里斯号”对直径约为500 m小行星“贝努”（Bennu）[16]的探测结果可知，其表面无细微风化层的存在，几乎完全由石块构成。相比而言，2016HO3 直径仅有40～100 m，且自转速度高，经计算发现其赤道处自转速度显著大于逃逸速度，实际上2016HO3 表面重力加速度方向是指向空间，而非地心，如图5所示。这意味着即便2016HO3表面存在细小的风化层颗粒，这些小颗粒也可能由于温度交变、太阳风与紫外线作用产生的静电等因素而松动，并由于高速自转的离开小行星表面，进入空间。因此2016HO3 表面形态应当以较大石块为主，甚至不排除独石结构。可以推断，2016HO3 表面颗粒度大，热惯量高。

图5 考虑自转后2016HO3表面重力加速度方向Fig.5 The direction of the surface gravity acceleration with the rotation of 2016HO3

由于2016HO3可能为S/Q/L类小行星，且其直径显著小于同为S类的小行星“丝川”（Itokawa），可以推断2016HO3热惯量很可能大于“丝川”的热惯量。“丝川”的热惯量为700 J·m-2·s-0.5·K-1[17]左右，本文取 700 J·m-2·s-0.5·K-1作为 2016HO3 热惯量的下限。同时，热惯量越高，小行星表面温度变化幅度越低，对附着采样时的工程设计更友好，因此本文分析中将取 2016HO3 热惯量下限 700 J·m-2·s-0.5·K-1进行后续分析，确保分析结果可覆盖最恶劣情况。

2 小行星的热建模方法

小行星热建模就是关于小行星表面和次表层的温度的计算，物理参数如反照率、热导率、热容量、发射率、密度和粗糙度，以及天体的形状，自转轴信息，乃至热历史都被考虑在内。粗略区分，热模型可以分成简化热模型和复杂热模型两类，其区别在于是否考虑小行星形状，是否忽略热传导（或简化其处理），是否处理表面粗糙度。小行星数据缺乏，简化热模型应用场景广泛，而随着红外观测数据以及小行星形状数据的增加，近10年来复杂热模型也得到广泛应用。本文将分别介绍两类热模型获得广泛应用的代表：近地小行星热模型与小行星热物理模型[20]。

2.1 近地小行星热模型

近地小行星热模型（Near Earth Asteroid Thermal Model，NEATM）[19]通常用于小行星形状和旋转的信息不足而无法使用小行星热物理模型的情况。通常，NEATM 允许对小行星直径和反照率进行初步而可靠的估计，但不提供热惯量或表面粗糙度的信息。典型的NEATM获得直径的精度为15%，反照率约为30%。其他简化热模型有标准热模型（STM）、快速旋转模型（FRM）。STM 和FRM 大部分情况下已经不再应用[20]。

NEATM 假设小行星为球形并且不直接考虑热惯量和表面粗糙度。表面温度由日照的瞬时热平衡给出，其与地表法线方向和太阳之间夹角的余弦成正比，而在夜间则取零。最高温度发生在日下点

其中：S⊙是1 AU处太阳辐射常量，取1 353 W/m；A为几何反照率；r是小行星到太阳的距离，单位为AU；参数η是考虑到粗糙表面引起的热发射增强效应，用于调整模型温度分布的参数，η也称为发射参数，对于高热惯量小行星，η值显著大于1 （例如，1.5～3，理论最大值约为3.5[21]），而对于低热惯量小行星，η=1，小行星表面粗糙度增加倾向于降低η的值，例如，某个主带小行星（MBA）的η约为0.8，表明该小行星热惯量低，表面粗糙度大（最小理论值为0.6～0.7[21]）；σ是斯忒藩-波尔兹曼常量，数值为5.670 51× 10-8Wm-2K-4；ε是发射率。

近地小行星热模型给出了太阳直射处最高温度。该模型认为太阳直射处的表面瞬间达到了热平衡，完全不考虑热容、热惯量、表面形状等因素，只通过一个发射参数η对多种因素进行调整，因此必然带来一定的不确定性。另外该模型只适合对温度上限进行估计，无法对永久阴影区，或者极夜区的温度进行分析。

表2给出了式（5）中关键参数范围判断结果。其中S⊙r-2代表了太阳能量的输入，这里分别分析了近地点和远地点之间的差异。由于2016HO3 几何反照率A与红外发射率ε无法获知，这里按照S 类小行星典型数值0.15～0.25 进行估算[22]。发射参数η受到粗糙度与热惯量的影响，而目前认为2016HO3 的热惯量大，粗糙度高，因此取η受范围1.5～3之间。

表2 2016HO3的近地小行星热模型（NEATM）参数选择Table 2 Parameter selection for the Near-Earth Asteroid Thermal Model（NEATM）of 2016HO3

综合以上结果，可以计算2016HO3 表面被太阳直射点的温度，见表3所示。而太阳直射点的温度代表了2016HO3 温度的上限值，非太阳直射点处的温度可以用

来计算，这里cosθ是太阳方向与所考虑平面夹角的余弦值。

表3 小行星2016HO3太阳直射点处温度分析Table 3 Temperature range analysis of asteroid 2016HO3 under direct sunlight K

2.2 小行星热物理模型

近地小行星热模型可以快速对温度上限进行估计，但是该模型过于简单，未考虑形状、发射面物理特性，未考虑热量的传输过程（纵向）；且只能用于估计温度上限，不能分析温度下限，不能用于计算瞬态温度。因此根据月球以及一些小行星的数据发展了小行星热物理模型（Thermal Physics Model）。该模型考虑了更多要素，可以将形状、表面粗糙度、热惯量等因素考虑在内。

所有TPM 都将小行星形状表示为（三角形）小平面的网格，该小平面围绕给定的自转轴旋转。一般而言，小行星的形状模型来自雷达观测或者光变曲线的反演，如果没有相应模型，则只能假设小行星为球形。具体每个平面可以选择不同的粗糙度模型，如图6所示，（b）代表陨石坑截面[23]，（c）代表高斯表面[24]，（d）代表分形表面[23]，选择合适的粗糙度模型可反映辐射表面多次辐射反射所导致对阳光吸收的增强以及再现热红外光束效应的方向性[21]。

考虑自转因素后，小行星表面温度处于周期性交变状态，随着深度增加，温度变化幅度约来越小。体现温度变化所影响深度参数被称为热特征尺度或者热肤深度，定义为

其中：κ是热导率；C是表面材料比热容；ρ是表面材料密度；ω是自转角速度。对于高速自转的小行星，典型的热特征长度仅有几个厘米。

图6 小行星形状以及表面粗糙程度模型［21］Fig.6 Asteroid shape and surface roughness model［21］

通常小平面网格的尺度远远大于热特征尺度ls的大小，在该条件下，对于小平面中心点处，表面输入的热流将主要向地下传导而非在地表横向传导。因此TPM模型通常忽略横向传热，仅考虑深度方向z上的一维热传导方程，其形式为

边界条件为

其中：A为几何反照率；F⊙太阳在小行星处辐射常数，F⊙=S⊙r-2；siψi是太阳光照角度，这里ψi是太阳光线与平面i法向量的余弦值，si则代表太阳光是否照射到平面i，如果没有，si为0，否则为1；Fscat是来自其它表面太阳光散射项，对于一个凸多面体而言，Fscat=0；Fred是表面其它部位的红外辐射，对于一个凸多面体而言，Fred=0；Ath是热辐射反照率；ε是小行星表面的辐射系数；κ是表面物质的热导率；C是表面物质的比热容；ρ是表面物质的密度。

为了求解，对上述方程进行了一些简化[21]，即

其中：ls是特征长度；ω是小行星的自转角速度，为ω=2π/T；Te是特征温度，Te相当于NEATM 中未经过发射参数η修正的太阳直射点下的温度。

上述方程简化为

其中：太阳光照角度项为p1=siψi；散射光项p2=红外辐射项无量纲的系数热惯量

若使用热物理模型，则至少需要获取2016HO3的形状、几何反照率、热惯量、密度、比热容、发射率、红外辐射、红外反照率以及粗糙程度等参数。

实际上，目前2016HO3 先验知识相当缺乏，特别是，缺乏2016HO3 的形状信息，只能初步假设2016HO3 的形状是球体，由于球体的几何特性，球面上不同位置相互不可见，因此，Fscat与Fred为0，还需要知道小行星的几何反照率、热惯量、密度、比热容、发射率等参数才能进行初步分析。

2.3 两种模型的比较

对于观测信息缺乏的2016HO3 而言，仅需要少量输入参数NEATM模型便可提供对表面温度上限的分析，且计算简便，适合于工程上对2016HO3 表面温度上限进行初步分析，但是需要注意，NEATM 提供的并不是精确结果，其更多被应用于对近地小行星的直径和反照率进行初步而可靠的估计。另外NEATM 无法估计温度下限，特别是处于极夜的表面。

TPM 考虑形状、表面粗糙度、热惯量等因素，考虑了纵向热传导过程，适合于分析昼夜温度，但是TPM 模型是用于分析自转条件下热环境的，同样也无法直接用于分析极夜条件下的温度。TPM 需要了解小行星更多基本参数，这恰是2016HO3 目前所缺乏的，但是随着后续观测数据增加，TPM 模型将提供对更加详细的热环境信息。

3 TPM模型计算与结果分析

3.1 TPM模型参数分析

2016HO3 多数参数未知，但是可以进行合理估计。在初步模型中，由于缺乏小行星形状的信息，这里暂时不考虑散射光与其它面的红外辐射，即认为p2、p3为0。其相关参数信息见表4。

表4 2016HO3的TPM算例取值Table 4 Parameters for the TPM of 2016HO3

小行星的自转周期来源于光变曲线测量结果，是确切知道参数；2016HO3的反照率A目前没有测量结果，目前只能按照S 类小行星平均反照率进行估计；小行星的密度ρ同样来自于已知密度的S类小行星的平均值；热惯量分析详见第1节；比热容数值参考自小行星“贝努”；特征长度、特征温度与Φ均源于以上数据的推导。

3.2 TPM的数值计算方法

采用数值求解对TPM 进行计算。数值解中不存在∞，这里以x=10（10倍的特征长度）代替∞深处

其中：u是x与τ的函数，以ui,j表示数据点；i代表位置，j代表时间。

上述问题就变成了一个设定了边界条件的ui,j网格的求解问题。这里需要注意的是，p1是τ的函数，这里存在一个自然的边界，即u从τ从0变化到2π后，对应小行星已经自转完一圈，应该恢复τ=0的状态，利用该边界条件，不用再去设定温度的初始值了，上述差分方程的迭代本质上就是对一个4对角线矩阵的求解。但是考虑到该问题属于非齐次边界条件（热辐射正比于温度的四次方T4），这里只能采用迭代的方式进行计算。取后向差分格式，可以得到

这里需要注意Δτ与Δx取值必须满足a<1，否则迭代将发散，对于边界条件，有

3.3 计算结果分析

假设2016HO3 是球状，自转轴方向是影响小行星表面温度分布的关键参数。但是由于自转轴数据的缺乏，这里研究小行星位于1 AU 处，自转轴相对于公转轨道面不同倾角下的小行星2016HO表面的温度分析，具体分布见图7。这里小行星经度与纬度采用小行星本体坐标系，由于缺乏小行星的形状信息，假设小行星是球体，图7是太阳直射180°经度时小行星表面全球瞬间温度分布；由于球体的对称性，图7也可看作小行星不同纬度经历一个自转周期后的温度变化（需要把横坐标改为时间）。

图7 不同自转轴倾角下温度分布Fig.7 Temperature distribution with different inclination angles of the rotation axis

图7中设定太阳直射180°经度处，APM 模型计算出的表面温度最高点将出现在经度约210°处，即考虑热惯量后，小行星温度最高点相对太阳直射点有一定延迟。另外，小行星自转轴倾角越接近90°，则太阳直射点温度越高，最高甚至可以接近特征温度393 K，这较NEATM的分析结果更大，而小行星自转轴倾角越接近0°（或者180°），则最高温度将降低至310 K，因此自转轴对温度上线影响显著。同时需要注意，图7中极夜区温度数值并无使用价值，TPM模型无法计算极夜区温度分布。

另外，由于2016HO3 自转速度快，小行星同一点处昼夜温度变化不大，本算例中最大值约为30 K左右。同时，由于高速自转，温差仅存在于表层，当深度超过8 cm，小行星温度几乎是一个恒定的值。

图8 温度/温差与深度关系Fig.8 Relationship between Temperature/Temperature difference and depth

公转轨道位置同样对2016HO3 温度影响显著，由于2016HO3 的偏心率到达0.1，近日点和远日点太阳辐射能量差距达到了40%，这将显著影响小行星的最高温度，具体计算结果见表6所示。这里需要注意，最高温差出现时自转轴倾角为0°，而最高温度出现时自转轴倾角为90°，两者不会同时出现。

表5的比较发现，自转轴是影响小行星表面最高温度以及温差关键要素。表7的比较发现，在2016HO3 分析中，TPM 可以反应自转轴对温度的影响，而由于没有使用自转轴这个参数，NEATM计算结果则介于自转轴两种极端情况中间。

表5 倾角与温差、最高温度的关系Table 5 Relationship between inclination and temperature difference/highest temperature

表6 公转轨道不同位置的最高温度与最大温差Table 6 The maximum temperature difference/temperatureat different positions of the orbit

表7 TPM模型与NEATM模型计算最高温度比较Table 7 The highest temperature between TPM and NEATM

如果取自转轴倾角0°，TPM 模型计算出的最高温度显著低于NEATM模型，这是由于2016HO3高热惯量、高转速有效平均了表面温度、降低了最高温度；NEATM 模型则认为受晒表面已经达到瞬间热平衡，这个条件只有在低转速或者低热惯量条件下才能实现。而自转轴倾角90°时，TPM 结果又显著高于NEATM，这是由于NEATM采用发射参数η对自转的效应进行了一定修正，而自转轴倾角90°时极区全天受照，等同于无自转，发射参数η在这种条件下又会降低对最高温度的预期。

4 极夜区的处理方式

上述计算仅适用于有能量输入的区域，即昼夜交替或者极昼的区域，而不适用于极夜区。不考虑横向热传导时，极夜区无外部能量输入，则达到平衡状态下表面温度是0 K，这显然不符合实际情况。实际上极区一个公转周期内经历的光照条件可以分为4个部分的循环：极昼→昼夜交替→极夜→昼夜交替→极昼，考虑热惯量后，便可计算出一个公转周期内极夜条件下最低温度。

2016HO3公转周期为一年，自转周期仅有0.467 h，两者相差了4个量级；由于自转周期差别，热特征长度上也相差2个数量级[见公式（7）]，因此自转和公转时热过程分析在时间、空间尺度上均有较大的差异，如果放到一块处理则由于时间、深度尺度不一，导致计算网格数量增加6个数量级，运算量过大，给计算带来困难。

自转和公转时热过程分析在时间空间尺度均有数量级的差别，两者可以分开进行处理。在处理公转问题时，可以将每个自转周期太阳光照提供的能量平均后作为公转分析时的输入，而在处理自转问题时，可以将公转计算得到的准稳态的地下的温度作为边界条件，叠加每日光照变化再计算表面温度，此时边界条件由原来的第二类边界条件式（8）变为第一类边界条件

阳光输入的平均值就是对昼夜效应进行平均，具体算法为

其中：小行星自转轴太阳光的夹角为β；所研究平面与自转轴夹角为γ。

图9给出了一个算例，即太阳光与自转轴夹角为60°，小行星的轨道的夏至点与近日点重合，小行星纬度为70°区域的光照变化情况：蓝色线代表每时刻光照情况，而红线代表了日均的光照，利用日均光照问题就可以计算出公转轨道上小行星日均温度随时间变化。

具体计算过程类似于3.2节中，参数可参见表4，但是需要注意自转周期由0.467个小时变为了366天，假设太阳光与自转轴夹角为60°，小行星纬度为70°区域一年日平均温度变化如下：

可以得到极区纬度70°处最低温度在一年最冷一天为137.9 K。

使用该方法可以解决极夜区温度的计算问题，但是同样需要注意该方法也存在一定缺陷：对于自转热分析，热特征深度仅有厘米量级，单元横向尺度远大于热特征深度，不考虑横向热传导是合理的；对于极夜区处理，将日均光照平均后计算一年温度变化，热特征深度将达到1 m，在这种条件下单元横向尺度与热特征深度接近，对于2016HO3 这样的小目标，忽略横向传热不再合理，横向热传导起到的作用是降低小行星表面同一时刻的温差，因此该方法分析出极区温度最低值将比实际情况偏低，工程上可以用作最低温度下限值使用。

未知自转轴朝向的条件下，分析了不同纬度下一年内小行星表面日均温度的变化，如图11所示，极夜条件下日均温度就是其瞬时温度，因此使用该方法计算出2016HO3温度下限值为115 K。