賈曉宇，楊 晨，王 彤，文 毅

（1.北京空間飛行器總體設(shè)計(jì)部，北京 100094；2.中國(guó)空間技術(shù)研究院，北京 100094）

引 言

2019年4月18日，中國(guó)國(guó)家航天局（China National Space Administration，CNSA）公布了小行星探測(cè)計(jì)劃：其中探測(cè)器將開(kāi)展對(duì)近地小行星2016HO3 的繞飛探測(cè)，隨后擇機(jī)附著在小行星表面并采集樣品，隨后攜帶樣品返回地球，并通過(guò)返回器將樣品送到地面。

小行星表面溫度信息對(duì)于小行星探測(cè)載荷設(shè)計(jì)，特別是有近地表操作的原位探測(cè)、樣品采集的任務(wù)十分必要。相比于采用一觸即走（Touch ＆ Go）采樣方式的“隼鳥(niǎo)號(hào)”（Hayabusa）、“隼鳥(niǎo)2號(hào)”（Hayabusa-2）[1]與“歐西里斯-雷克斯號(hào)”（OSIRIS-Rex）[2]，我國(guó)小行星探測(cè)將附著在小行星表面，采集樣品，小行星表面熱環(huán)境更為重要。

小行星表面熱輻射所產(chǎn)生的凈力與力矩會(huì)導(dǎo)致亞爾科夫斯基效應(yīng)與YORP 效應(yīng)，影響小行星軌道演化，對(duì)于足夠小的天體，這兩種效應(yīng)更加顯著[3-4]，準(zhǔn)確分析亞爾科夫斯基效應(yīng)與YORP效應(yīng)需要小行星溫度分布信息。溫度及其在小行星整個(gè)生命周期中的演變可以改變小行星的表面組成和性質(zhì)[5-6]。因此，小行星表面溫度信息對(duì)于研究小行星軌道與表面物質(zhì)演化具有重要意義。

本文將結(jié)合小行星2016HO3 的特征以及地面初步的觀測(cè)數(shù)據(jù)，開(kāi)展小行星2016HO3 溫度分析，獲得其上下限，并對(duì)其熱特征進(jìn)行初步總結(jié)。

1 目標(biāo)小行星2016HO3基本信息

小行星2016HO3（編號(hào)469219，名稱(chēng)為komo’oalewa）是由位于夏威夷海勒卡拉山1.8 m 口徑的“泛星1 號(hào)”小行星巡天望遠(yuǎn)鏡2016年4月27日發(fā)現(xiàn)的。2016HO3 是一顆阿波羅型小行星。圖1展示了2016HO3 相對(duì)地球運(yùn)動(dòng)軌跡，2016HO3 與地球以相同的周期共同繞太陽(yáng)公轉(zhuǎn)，以地球角度看，2016HO3也在繞地球轉(zhuǎn)動(dòng)，此類(lèi)天體被稱(chēng)為地球準(zhǔn)衛(wèi)星（Earth Quasi-Satellite）。2016HO3 是目前發(fā)現(xiàn)的最穩(wěn)定的地球準(zhǔn)衛(wèi)星，至少未來(lái)100年內(nèi)不會(huì)飛離地球[7]。

目前地面觀測(cè)未能直接得到2016HO3 的類(lèi)型與大小，初步估計(jì)2016HO3可能是S/Q/L類(lèi)小行星，其絕對(duì)星等Hv=-24.3±0.5[8-10]，反照率α=0.2[8]（但該數(shù)據(jù)并非觀測(cè)所得，而是對(duì)S 類(lèi)小行星的估計(jì)值），利用公式

可以計(jì)算出2016HO3 的直徑為41 m[10]，而JPL的Agle 根據(jù)雷達(dá)探測(cè)數(shù)據(jù)估計(jì)2016HO3 的直徑為40～100 m。本文在后續(xù)分析中將認(rèn)為2016HO3的直徑是40～100 m。

圖1 2016HO3的軌道示意圖［7］Fig.1 Orbit of 2016HO3［7］

表1 2016HO3的軌道參數(shù)［8］Tabel1 Orbital parameters of 2016HO3［8］

圖2 2016HO3的光變曲線(xiàn)［11］Fig.2 Light curve of 2016HO3［11］

根據(jù)光變曲線(xiàn)（參見(jiàn)圖2），目前可以得到2016HO3 較為準(zhǔn)確的自轉(zhuǎn)周期為0.467 ± 0.008 hr（28?1.2" ± 28.8"）[10]。目前無(wú)法確定2016HO3 自轉(zhuǎn)軸的方向，甚至無(wú)法判定2016HO3 是否存在穩(wěn)定的自轉(zhuǎn)軸。Burns （1973）曾經(jīng)提出一個(gè)公式用于計(jì)算小天體自轉(zhuǎn)軸趨于穩(wěn)定所需時(shí)間尺度τdamp[12]

其中：Q代表品質(zhì)因子，即每圈轉(zhuǎn)動(dòng)所損失的能量；μ代表物質(zhì)硬度；ρ是小天體的密度；ω是小天體自轉(zhuǎn)角速度；R是小天體半徑；K32是一個(gè)與天體形狀相關(guān)的無(wú)量綱量，對(duì)于球狀小天體，K32接近于0.01，對(duì)于長(zhǎng)橢球狀小天體，K32接近于0.1。

Harris （1994）給出一個(gè)更簡(jiǎn)化的公式計(jì)算小天體自轉(zhuǎn)軸趨于穩(wěn)定所需時(shí)間τdamp[13]

其中：P是小天體自轉(zhuǎn)周期，單位為h；D是小天體直徑，單位為km，C是一個(gè)無(wú)量綱量，取17±2.5；τdamp的單位是 109年。則可以計(jì)算出2016HO3 的τdamp為2×106～1.3×107年，如果2016HO3的自轉(zhuǎn)軸不穩(wěn)定，則意味著其自轉(zhuǎn)軸需要數(shù)百萬(wàn)年甚至更長(zhǎng)時(shí)間才能穩(wěn)定，目前尚不清楚2016HO3 演化歷史，2016HO3的自轉(zhuǎn)軸是否穩(wěn)定無(wú)法判斷。

小行星的熱慣量是計(jì)算表面溫度分布的關(guān)鍵參數(shù)，但是目前缺乏2016HO3 熱紅外觀測(cè)數(shù)據(jù)對(duì)其溫度、熱慣量進(jìn)行直接分析。

熱慣量與小行星表面的顆粒大小有關(guān)，小行星表面物質(zhì)顆粒越大，熱慣量也越大，如月球表面等土壤顆粒非常細(xì)小，其熱慣量?jī)H為40～50 J·m-2·s-0.5·K-1，小行星數(shù)量巨大，表面形態(tài)不同，熱慣量的跨度較大，根據(jù)目前公布的數(shù)據(jù)，最小的可以到1 J·m-2·-0.5（10199 Chariklo），大的可以到數(shù)百[14]。盡管根據(jù)“歐西里斯號(hào)”對(duì)貝努最新觀測(cè)表明表面顆粒直徑并非決定小行星表面熱慣量唯一決定因素，但表面顆粒直徑越大，小行星熱慣量越大的結(jié)論仍然能夠成立[15]。同時(shí)，小天體直徑越小，其表面引力越微弱，則微小的風(fēng)化層顆粒越難以保存在其表面上，因此小行星直徑越大，其表面的熱慣量也就越大。

圖3 小行星的直徑與熱慣量關(guān)系[14]Fig.3 The relationship between the diameter and the thermal inertia of an asteroid[14]

Delbo 的研究給出了一個(gè)小行星直徑（D）與熱慣量（Γ）之間近似的擬合關(guān)系

其 中 ：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.。需要注意，該數(shù)據(jù)來(lái)源于公式（4）的外推，對(duì)于該數(shù)值應(yīng)采取謹(jǐn)慎的態(tài)度。

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

根據(jù)近期“歐西里斯號(hào)”對(duì)直徑約為500 m小行星“貝努”（Bennu）[16]的探測(cè)結(jié)果可知，其表面無(wú)細(xì)微風(fēng)化層的存在，幾乎完全由石塊構(gòu)成。相比而言，2016HO3 直徑僅有40～100 m，且自轉(zhuǎn)速度高，經(jīng)計(jì)算發(fā)現(xiàn)其赤道處自轉(zhuǎn)速度顯著大于逃逸速度，實(shí)際上2016HO3 表面重力加速度方向是指向空間，而非地心，如圖5所示。這意味著即便2016HO3表面存在細(xì)小的風(fēng)化層顆粒，這些小顆粒也可能由于溫度交變、太陽(yáng)風(fēng)與紫外線(xiàn)作用產(chǎn)生的靜電等因素而松動(dòng)，并由于高速自轉(zhuǎn)的離開(kāi)小行星表面，進(jìn)入空間。因此2016HO3 表面形態(tài)應(yīng)當(dāng)以較大石塊為主，甚至不排除獨(dú)石結(jié)構(gòu)?？梢酝茢?，2016HO3 表面顆粒度大，熱慣量高。

圖5 考慮自轉(zhuǎn)后2016HO3表面重力加速度方向Fig.5 The direction of the surface gravity acceleration with the rotation of 2016HO3

由于2016HO3可能為S/Q/L類(lèi)小行星，且其直徑顯著小于同為S類(lèi)的小行星“絲川”（Itokawa），可以推斷2016HO3熱慣量很可能大于“絲川”的熱慣量?！敖z川”的熱慣量為700 J·m-2·s-0.5·K-1[17]左右，本文取 700 J·m-2·s-0.5·K-1作為 2016HO3 熱慣量的下限。同時(shí)，熱慣量越高，小行星表面溫度變化幅度越低，對(duì)附著采樣時(shí)的工程設(shè)計(jì)更友好，因此本文分析中將取 2016HO3 熱慣量下限 700 J·m-2·s-0.5·K-1進(jìn)行后續(xù)分析，確保分析結(jié)果可覆蓋最?lèi)毫忧闆r。

2 小行星的熱建模方法

小行星熱建模就是關(guān)于小行星表面和次表層的溫度的計(jì)算，物理參數(shù)如反照率、熱導(dǎo)率、熱容量、發(fā)射率、密度和粗糙度，以及天體的形狀，自轉(zhuǎn)軸信息，乃至熱歷史都被考慮在內(nèi)。粗略區(qū)分，熱模型可以分成簡(jiǎn)化熱模型和復(fù)雜熱模型兩類(lèi)，其區(qū)別在于是否考慮小行星形狀，是否忽略熱傳導(dǎo)（或簡(jiǎn)化其處理），是否處理表面粗糙度。小行星數(shù)據(jù)缺乏，簡(jiǎn)化熱模型應(yīng)用場(chǎng)景廣泛，而隨著紅外觀測(cè)數(shù)據(jù)以及小行星形狀數(shù)據(jù)的增加，近10年來(lái)復(fù)雜熱模型也得到廣泛應(yīng)用。本文將分別介紹兩類(lèi)熱模型獲得廣泛應(yīng)用的代表：近地小行星熱模型與小行星熱物理模型[20]。

2.1 近地小行星熱模型

近地小行星熱模型（Near Earth Asteroid Thermal Model，NEATM）[19]通常用于小行星形狀和旋轉(zhuǎn)的信息不足而無(wú)法使用小行星熱物理模型的情況。通常，NEATM 允許對(duì)小行星直徑和反照率進(jìn)行初步而可靠的估計(jì)，但不提供熱慣量或表面粗糙度的信息。典型的NEATM獲得直徑的精度為15%，反照率約為30%。其他簡(jiǎn)化熱模型有標(biāo)準(zhǔn)熱模型（STM）、快速旋轉(zhuǎn)模型（FRM）。STM 和FRM 大部分情況下已經(jīng)不再應(yīng)用[20]。

NEATM 假設(shè)小行星為球形并且不直接考慮熱慣量和表面粗糙度。表面溫度由日照的瞬時(shí)熱平衡給出，其與地表法線(xiàn)方向和太陽(yáng)之間夾角的余弦成正比，而在夜間則取零。最高溫度發(fā)生在日下點(diǎn)

其中：S⊙是1 AU處太陽(yáng)輻射常量，取1 353 W/m；A為幾何反照率；r是小行星到太陽(yáng)的距離，單位為AU；參數(shù)η是考慮到粗糙表面引起的熱發(fā)射增強(qiáng)效應(yīng)，用于調(diào)整模型溫度分布的參數(shù)，η也稱(chēng)為發(fā)射參數(shù)，對(duì)于高熱慣量小行星，η值顯著大于1 （例如，1.5～3，理論最大值約為3.5[21]），而對(duì)于低熱慣量小行星，η=1，小行星表面粗糙度增加傾向于降低η的值，例如，某個(gè)主帶小行星（MBA）的η約為0.8，表明該小行星熱慣量低，表面粗糙度大（最小理論值為0.6～0.7[21]）；σ是斯忒藩-波爾茲曼常量，數(shù)值為5.670 51× 10-8Wm-2K-4；ε是發(fā)射率。

近地小行星熱模型給出了太陽(yáng)直射處最高溫度。該模型認(rèn)為太陽(yáng)直射處的表面瞬間達(dá)到了熱平衡，完全不考慮熱容、熱慣量、表面形狀等因素，只通過(guò)一個(gè)發(fā)射參數(shù)η對(duì)多種因素進(jìn)行調(diào)整，因此必然帶來(lái)一定的不確定性。另外該模型只適合對(duì)溫度上限進(jìn)行估計(jì)，無(wú)法對(duì)永久陰影區(qū)，或者極夜區(qū)的溫度進(jìn)行分析。

表2給出了式（5）中關(guān)鍵參數(shù)范圍判斷結(jié)果。其中S⊙r-2代表了太陽(yáng)能量的輸入，這里分別分析了近地點(diǎn)和遠(yuǎn)地點(diǎn)之間的差異。由于2016HO3 幾何反照率A與紅外發(fā)射率ε無(wú)法獲知，這里按照S 類(lèi)小行星典型數(shù)值0.15～0.25 進(jìn)行估算[22]。發(fā)射參數(shù)η受到粗糙度與熱慣量的影響，而目前認(rèn)為2016HO3 的熱慣量大，粗糙度高，因此取η受范圍1.5～3之間。

表2 2016HO3的近地小行星熱模型（NEATM）參數(shù)選擇Table 2 Parameter selection for the Near-Earth Asteroid Thermal Model（NEATM）of 2016HO3

綜合以上結(jié)果，可以計(jì)算2016HO3 表面被太陽(yáng)直射點(diǎn)的溫度，見(jiàn)表3所示。而太陽(yáng)直射點(diǎn)的溫度代表了2016HO3 溫度的上限值，非太陽(yáng)直射點(diǎn)處的溫度可以用

來(lái)計(jì)算，這里cosθ是太陽(yáng)方向與所考慮平面夾角的余弦值。

表3 小行星2016HO3太陽(yáng)直射點(diǎn)處溫度分析Table 3 Temperature range analysis of asteroid 2016HO3 under direct sunlight K

2.2 小行星熱物理模型

近地小行星熱模型可以快速對(duì)溫度上限進(jìn)行估計(jì)，但是該模型過(guò)于簡(jiǎn)單，未考慮形狀、發(fā)射面物理特性，未考慮熱量的傳輸過(guò)程（縱向）；且只能用于估計(jì)溫度上限，不能分析溫度下限，不能用于計(jì)算瞬態(tài)溫度。因此根據(jù)月球以及一些小行星的數(shù)據(jù)發(fā)展了小行星熱物理模型（Thermal Physics Model）。該模型考慮了更多要素，可以將形狀、表面粗糙度、熱慣量等因素考慮在內(nèi)。

所有TPM 都將小行星形狀表示為（三角形）小平面的網(wǎng)格，該小平面圍繞給定的自轉(zhuǎn)軸旋轉(zhuǎn)。一般而言，小行星的形狀模型來(lái)自雷達(dá)觀測(cè)或者光變曲線(xiàn)的反演，如果沒(méi)有相應(yīng)模型，則只能假設(shè)小行星為球形。具體每個(gè)平面可以選擇不同的粗糙度模型，如圖6所示，（b）代表隕石坑截面[23]，（c）代表高斯表面[24]，（d）代表分形表面[23]，選擇合適的粗糙度模型可反映輻射表面多次輻射反射所導(dǎo)致對(duì)陽(yáng)光吸收的增強(qiáng)以及再現(xiàn)熱紅外光束效應(yīng)的方向性[21]。

考慮自轉(zhuǎn)因素后，小行星表面溫度處于周期性交變狀態(tài)，隨著深度增加，溫度變化幅度約來(lái)越小。體現(xiàn)溫度變化所影響深度參數(shù)被稱(chēng)為熱特征尺度或者熱膚深度，定義為

其中：κ是熱導(dǎo)率；C是表面材料比熱容；ρ是表面材料密度；ω是自轉(zhuǎn)角速度。對(duì)于高速自轉(zhuǎn)的小行星，典型的熱特征長(zhǎng)度僅有幾個(gè)厘米。

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

通常小平面網(wǎng)格的尺度遠(yuǎn)遠(yuǎn)大于熱特征尺度ls的大小，在該條件下，對(duì)于小平面中心點(diǎn)處，表面輸入的熱流將主要向地下傳導(dǎo)而非在地表橫向傳導(dǎo)。因此TPM模型通常忽略橫向傳熱，僅考慮深度方向z上的一維熱傳導(dǎo)方程，其形式為

邊界條件為

其中：A為幾何反照率；F⊙太陽(yáng)在小行星處輻射常數(shù)，F(xiàn)⊙=S⊙r-2；siψi是太陽(yáng)光照角度，這里ψi是太陽(yáng)光線(xiàn)與平面i法向量的余弦值，si則代表太陽(yáng)光是否照射到平面i，如果沒(méi)有，si為0，否則為1；Fscat是來(lái)自其它表面太陽(yáng)光散射項(xiàng)，對(duì)于一個(gè)凸多面體而言，F(xiàn)scat=0；Fred是表面其它部位的紅外輻射，對(duì)于一個(gè)凸多面體而言，F(xiàn)red=0；Ath是熱輻射反照率；ε是小行星表面的輻射系數(shù)；κ是表面物質(zhì)的熱導(dǎo)率；C是表面物質(zhì)的比熱容；ρ是表面物質(zhì)的密度。

為了求解，對(duì)上述方程進(jìn)行了一些簡(jiǎn)化[21]，即

其中：ls是特征長(zhǎng)度；ω是小行星的自轉(zhuǎn)角速度，為ω=2π/T；Te是特征溫度，Te相當(dāng)于NEATM 中未經(jīng)過(guò)發(fā)射參數(shù)η修正的太陽(yáng)直射點(diǎn)下的溫度。

上述方程簡(jiǎn)化為

其中：太陽(yáng)光照角度項(xiàng)為p1=siψi；散射光項(xiàng)p2=紅外輻射項(xiàng)無(wú)量綱的系數(shù)熱慣量

若使用熱物理模型，則至少需要獲取2016HO3的形狀、幾何反照率、熱慣量、密度、比熱容、發(fā)射率、紅外輻射、紅外反照率以及粗糙程度等參數(shù)。

實(shí)際上，目前2016HO3 先驗(yàn)知識(shí)相當(dāng)缺乏，特別是，缺乏2016HO3 的形狀信息，只能初步假設(shè)2016HO3 的形狀是球體，由于球體的幾何特性，球面上不同位置相互不可見(jiàn)，因此，F(xiàn)scat與Fred為0，還需要知道小行星的幾何反照率、熱慣量、密度、比熱容、發(fā)射率等參數(shù)才能進(jìn)行初步分析。

2.3 兩種模型的比較

對(duì)于觀測(cè)信息缺乏的2016HO3 而言，僅需要少量輸入?yún)?shù)NEATM模型便可提供對(duì)表面溫度上限的分析，且計(jì)算簡(jiǎn)便，適合于工程上對(duì)2016HO3 表面溫度上限進(jìn)行初步分析，但是需要注意，NEATM 提供的并不是精確結(jié)果，其更多被應(yīng)用于對(duì)近地小行星的直徑和反照率進(jìn)行初步而可靠的估計(jì)。另外NEATM 無(wú)法估計(jì)溫度下限，特別是處于極夜的表面。

TPM 考慮形狀、表面粗糙度、熱慣量等因素，考慮了縱向熱傳導(dǎo)過(guò)程，適合于分析晝夜溫度，但是TPM 模型是用于分析自轉(zhuǎn)條件下熱環(huán)境的，同樣也無(wú)法直接用于分析極夜條件下的溫度。TPM 需要了解小行星更多基本參數(shù)，這恰是2016HO3 目前所缺乏的，但是隨著后續(xù)觀測(cè)數(shù)據(jù)增加，TPM 模型將提供對(duì)更加詳細(xì)的熱環(huán)境信息。

3 TPM模型計(jì)算與結(jié)果分析

3.1 TPM模型參數(shù)分析

2016HO3 多數(shù)參數(shù)未知，但是可以進(jìn)行合理估計(jì)。在初步模型中，由于缺乏小行星形狀的信息，這里暫時(shí)不考慮散射光與其它面的紅外輻射，即認(rèn)為p2、p3為0。其相關(guān)參數(shù)信息見(jiàn)表4。

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

小行星的自轉(zhuǎn)周期來(lái)源于光變曲線(xiàn)測(cè)量結(jié)果，是確切知道參數(shù)；2016HO3的反照率A目前沒(méi)有測(cè)量結(jié)果，目前只能按照S 類(lèi)小行星平均反照率進(jìn)行估計(jì)；小行星的密度ρ同樣來(lái)自于已知密度的S類(lèi)小行星的平均值；熱慣量分析詳見(jiàn)第1節(jié)；比熱容數(shù)值參考自小行星“貝努”；特征長(zhǎng)度、特征溫度與Φ均源于以上數(shù)據(jù)的推導(dǎo)。

3.2 TPM的數(shù)值計(jì)算方法

采用數(shù)值求解對(duì)TPM 進(jìn)行計(jì)算。數(shù)值解中不存在∞，這里以x=10（10倍的特征長(zhǎng)度）代替∞深處

其中：u是x與τ的函數(shù)，以u(píng)i,j表示數(shù)據(jù)點(diǎn)；i代表位置，j代表時(shí)間。

上述問(wèn)題就變成了一個(gè)設(shè)定了邊界條件的ui,j網(wǎng)格的求解問(wèn)題。這里需要注意的是，p1是τ的函數(shù)，這里存在一個(gè)自然的邊界，即u從τ從0變化到2π后，對(duì)應(yīng)小行星已經(jīng)自轉(zhuǎn)完一圈，應(yīng)該恢復(fù)τ=0的狀態(tài)，利用該邊界條件，不用再去設(shè)定溫度的初始值了，上述差分方程的迭代本質(zhì)上就是對(duì)一個(gè)4對(duì)角線(xiàn)矩陣的求解。但是考慮到該問(wèn)題屬于非齊次邊界條件（熱輻射正比于溫度的四次方T4），這里只能采用迭代的方式進(jìn)行計(jì)算。取后向差分格式，可以得到

這里需要注意Δτ與Δx取值必須滿(mǎn)足a<1，否則迭代將發(fā)散，對(duì)于邊界條件，有

3.3 計(jì)算結(jié)果分析

假設(shè)2016HO3 是球狀，自轉(zhuǎn)軸方向是影響小行星表面溫度分布的關(guān)鍵參數(shù)。但是由于自轉(zhuǎn)軸數(shù)據(jù)的缺乏，這里研究小行星位于1 AU 處，自轉(zhuǎn)軸相對(duì)于公轉(zhuǎn)軌道面不同傾角下的小行星2016HO表面的溫度分析，具體分布見(jiàn)圖7。這里小行星經(jīng)度與緯度采用小行星本體坐標(biāo)系，由于缺乏小行星的形狀信息，假設(shè)小行星是球體，圖7是太陽(yáng)直射180°經(jīng)度時(shí)小行星表面全球瞬間溫度分布；由于球體的對(duì)稱(chēng)性，圖7也可看作小行星不同緯度經(jīng)歷一個(gè)自轉(zhuǎn)周期后的溫度變化（需要把橫坐標(biāo)改為時(shí)間）。

圖7 不同自轉(zhuǎn)軸傾角下溫度分布Fig.7 Temperature distribution with different inclination angles of the rotation axis

圖7中設(shè)定太陽(yáng)直射180°經(jīng)度處，APM 模型計(jì)算出的表面溫度最高點(diǎn)將出現(xiàn)在經(jīng)度約210°處，即考慮熱慣量后，小行星溫度最高點(diǎn)相對(duì)太陽(yáng)直射點(diǎn)有一定延遲。另外，小行星自轉(zhuǎn)軸傾角越接近90°，則太陽(yáng)直射點(diǎn)溫度越高，最高甚至可以接近特征溫度393 K，這較NEATM的分析結(jié)果更大，而小行星自轉(zhuǎn)軸傾角越接近0°（或者180°），則最高溫度將降低至310 K，因此自轉(zhuǎn)軸對(duì)溫度上線(xiàn)影響顯著。同時(shí)需要注意，圖7中極夜區(qū)溫度數(shù)值并無(wú)使用價(jià)值，TPM模型無(wú)法計(jì)算極夜區(qū)溫度分布。

另外，由于2016HO3 自轉(zhuǎn)速度快，小行星同一點(diǎn)處晝夜溫度變化不大，本算例中最大值約為30 K左右。同時(shí)，由于高速自轉(zhuǎn)，溫差僅存在于表層，當(dāng)深度超過(guò)8 cm，小行星溫度幾乎是一個(gè)恒定的值。

圖8 溫度/溫差與深度關(guān)系Fig.8 Relationship between Temperature/Temperature difference and depth

公轉(zhuǎn)軌道位置同樣對(duì)2016HO3 溫度影響顯著，由于2016HO3 的偏心率到達(dá)0.1，近日點(diǎn)和遠(yuǎn)日點(diǎn)太陽(yáng)輻射能量差距達(dá)到了40%，這將顯著影響小行星的最高溫度，具體計(jì)算結(jié)果見(jiàn)表6所示。這里需要注意，最高溫差出現(xiàn)時(shí)自轉(zhuǎn)軸傾角為0°，而最高溫度出現(xiàn)時(shí)自轉(zhuǎn)軸傾角為90°，兩者不會(huì)同時(shí)出現(xiàn)。

表5的比較發(fā)現(xiàn)，自轉(zhuǎn)軸是影響小行星表面最高溫度以及溫差關(guān)鍵要素。表7的比較發(fā)現(xiàn)，在2016HO3 分析中，TPM 可以反應(yīng)自轉(zhuǎn)軸對(duì)溫度的影響，而由于沒(méi)有使用自轉(zhuǎn)軸這個(gè)參數(shù)，NEATM計(jì)算結(jié)果則介于自轉(zhuǎn)軸兩種極端情況中間。

表5 傾角與溫差、最高溫度的關(guān)系Table 5 Relationship between inclination and temperature difference/highest temperature

表6 公轉(zhuǎn)軌道不同位置的最高溫度與最大溫差Table 6 The maximum temperature difference/temperatureat different positions of the orbit

表7 TPM模型與NEATM模型計(jì)算最高溫度比較Table 7 The highest temperature between TPM and NEATM

如果取自轉(zhuǎn)軸傾角0°，TPM 模型計(jì)算出的最高溫度顯著低于NEATM模型，這是由于2016HO3高熱慣量、高轉(zhuǎn)速有效平均了表面溫度、降低了最高溫度；NEATM 模型則認(rèn)為受曬表面已經(jīng)達(dá)到瞬間熱平衡，這個(gè)條件只有在低轉(zhuǎn)速或者低熱慣量條件下才能實(shí)現(xiàn)。而自轉(zhuǎn)軸傾角90°時(shí)，TPM 結(jié)果又顯著高于NEATM，這是由于NEATM采用發(fā)射參數(shù)η對(duì)自轉(zhuǎn)的效應(yīng)進(jìn)行了一定修正，而自轉(zhuǎn)軸傾角90°時(shí)極區(qū)全天受照，等同于無(wú)自轉(zhuǎn)，發(fā)射參數(shù)η在這種條件下又會(huì)降低對(duì)最高溫度的預(yù)期。

4 極夜區(qū)的處理方式

上述計(jì)算僅適用于有能量輸入的區(qū)域，即晝夜交替或者極晝的區(qū)域，而不適用于極夜區(qū)。不考慮橫向熱傳導(dǎo)時(shí)，極夜區(qū)無(wú)外部能量輸入，則達(dá)到平衡狀態(tài)下表面溫度是0 K，這顯然不符合實(shí)際情況。實(shí)際上極區(qū)一個(gè)公轉(zhuǎn)周期內(nèi)經(jīng)歷的光照條件可以分為4個(gè)部分的循環(huán)：極晝→晝夜交替→極夜→晝夜交替→極晝，考慮熱慣量后，便可計(jì)算出一個(gè)公轉(zhuǎn)周期內(nèi)極夜條件下最低溫度。

2016HO3公轉(zhuǎn)周期為一年，自轉(zhuǎn)周期僅有0.467 h，兩者相差了4個(gè)量級(jí)；由于自轉(zhuǎn)周期差別，熱特征長(zhǎng)度上也相差2個(gè)數(shù)量級(jí)[見(jiàn)公式（7）]，因此自轉(zhuǎn)和公轉(zhuǎn)時(shí)熱過(guò)程分析在時(shí)間、空間尺度上均有較大的差異，如果放到一塊處理則由于時(shí)間、深度尺度不一，導(dǎo)致計(jì)算網(wǎng)格數(shù)量增加6個(gè)數(shù)量級(jí)，運(yùn)算量過(guò)大，給計(jì)算帶來(lái)困難。

自轉(zhuǎn)和公轉(zhuǎn)時(shí)熱過(guò)程分析在時(shí)間空間尺度均有數(shù)量級(jí)的差別，兩者可以分開(kāi)進(jìn)行處理。在處理公轉(zhuǎn)問(wèn)題時(shí)，可以將每個(gè)自轉(zhuǎn)周期太陽(yáng)光照提供的能量平均后作為公轉(zhuǎn)分析時(shí)的輸入，而在處理自轉(zhuǎn)問(wèn)題時(shí)，可以將公轉(zhuǎn)計(jì)算得到的準(zhǔn)穩(wěn)態(tài)的地下的溫度作為邊界條件，疊加每日光照變化再計(jì)算表面溫度，此時(shí)邊界條件由原來(lái)的第二類(lèi)邊界條件式（8）變?yōu)榈谝活?lèi)邊界條件

陽(yáng)光輸入的平均值就是對(duì)晝夜效應(yīng)進(jìn)行平均，具體算法為

其中：小行星自轉(zhuǎn)軸太陽(yáng)光的夾角為β；所研究平面與自轉(zhuǎn)軸夾角為γ。

圖9給出了一個(gè)算例，即太陽(yáng)光與自轉(zhuǎn)軸夾角為60°，小行星的軌道的夏至點(diǎn)與近日點(diǎn)重合，小行星緯度為70°區(qū)域的光照變化情況：藍(lán)色線(xiàn)代表每時(shí)刻光照情況，而紅線(xiàn)代表了日均的光照，利用日均光照問(wèn)題就可以計(jì)算出公轉(zhuǎn)軌道上小行星日均溫度隨時(shí)間變化。

具體計(jì)算過(guò)程類(lèi)似于3.2節(jié)中，參數(shù)可參見(jiàn)表4，但是需要注意自轉(zhuǎn)周期由0.467個(gè)小時(shí)變?yōu)榱?66天，假設(shè)太陽(yáng)光與自轉(zhuǎn)軸夾角為60°，小行星緯度為70°區(qū)域一年日平均溫度變化如下：

可以得到極區(qū)緯度70°處最低溫度在一年最冷一天為137.9 K。

使用該方法可以解決極夜區(qū)溫度的計(jì)算問(wèn)題，但是同樣需要注意該方法也存在一定缺陷：對(duì)于自轉(zhuǎn)熱分析，熱特征深度僅有厘米量級(jí)，單元橫向尺度遠(yuǎn)大于熱特征深度，不考慮橫向熱傳導(dǎo)是合理的；對(duì)于極夜區(qū)處理，將日均光照平均后計(jì)算一年溫度變化，熱特征深度將達(dá)到1 m，在這種條件下單元橫向尺度與熱特征深度接近，對(duì)于2016HO3 這樣的小目標(biāo)，忽略橫向傳熱不再合理，橫向熱傳導(dǎo)起到的作用是降低小行星表面同一時(shí)刻的溫差，因此該方法分析出極區(qū)溫度最低值將比實(shí)際情況偏低，工程上可以用作最低溫度下限值使用。

未知自轉(zhuǎn)軸朝向的條件下，分析了不同緯度下一年內(nèi)小行星表面日均溫度的變化，如圖11所示，極夜條件下日均溫度就是其瞬時(shí)溫度，因此使用該方法計(jì)算出2016HO3溫度下限值為115 K。

圖9 公轉(zhuǎn)后的平均光照因子Fig.9 Average illumination factor of the orbit

圖10 極夜區(qū)一年日平均溫度分布（緯度為70°）Fig.10 Average daily temperature for Polar zone（latitude 70°）

圖11 不同自轉(zhuǎn)軸傾角下不同緯年日均溫度Fig.11 Daily average temperature at different latitudes with different inclination angles

5 結(jié) 論

本文介紹了小行星的兩個(gè)熱模型：近地小行星熱模型與小行星熱物理模型，對(duì)2016HO3 表面溫度進(jìn)行初步分析，得到主要結(jié)論：

1） 根據(jù)近地小行星熱模型（NEATM），2016HO3 的表面溫度上限最高不會(huì)超過(guò)380 K（約100 ℃）；根據(jù)小行星熱物理模型（TPM），2016HO3的表面溫度上限為412 K，但是僅發(fā)生在其自轉(zhuǎn)軸傾角為90°且夏至點(diǎn)與近日點(diǎn)重合的條件；

根據(jù)小行星熱物理模型計(jì)算，發(fā)現(xiàn)2016HO3 表面不同位置溫差雖大，由于高速自轉(zhuǎn)效果，小行星表面同一點(diǎn)處的晝夜溫差不會(huì)超過(guò)31 ℃。

鑒于兩個(gè)模型均無(wú)法處理極夜區(qū)的溫度，本文在TPM 模型基礎(chǔ)上提出以年平均思路計(jì)算極夜區(qū)溫度下限為115 K，該方法存在一定缺陷，計(jì)算結(jié)果較實(shí)際偏低，計(jì)算結(jié)果作為下限值可以為工程提供參考。

小行星表面溫度分布以及溫度范圍的決定性因素之一是自轉(zhuǎn)軸方向，而2016HO3 的自轉(zhuǎn)軸傾角目前尚未獲知，導(dǎo)致目前分析2016HO3溫度包絡(luò)為115～412 K。需等待地面觀測(cè)或者空間觀測(cè)結(jié)果進(jìn)行進(jìn)一步分析，獲得更加精確的結(jié)果。

由于兩個(gè)熱模型尚且無(wú)法分析極夜區(qū)溫度，同時(shí)對(duì)于“極區(qū)”需要將公轉(zhuǎn)要素考慮在內(nèi)，最后本文在TPM 模型基礎(chǔ)上提出以年平均思路計(jì)算極夜極區(qū)溫度的分析方法，但是對(duì)于2016HO3 尺寸僅有40～100 m的小行星仍有不足，后續(xù)需考慮橫向傳熱因素影響，改進(jìn)計(jì)算方法，獲得更準(zhǔn)確的溫度下限估計(jì)。