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

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

引 言

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

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

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

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

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

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

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

可以計(jì)算出2016HO3 的直徑為41 m[10]，而JPL的Agle 根據(jù)雷達(dá)探測數(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的光變曲線［11］Fig.2 Light curve of 2016HO3［11］

根據(jù)光變曲線（參見圖2），目前可以得到2016HO3 較為準(zhǔn)確的自轉(zhuǎn)周期為0.467 ± 0.008 hr（28?1.2" ± 28.8"）[10]。目前無法確定2016HO3 自轉(zhuǎn)軸的方向，甚至無法判定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)的無量綱量，對于球狀小天體，K32接近于0.01，對于長橢球狀小天體，K32接近于0.1。

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

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

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

熱慣量與小行星表面的顆粒大小有關(guān)，小行星表面物質(zhì)顆粒越大，熱慣量也越大，如月球表面等土壤顆粒非常細(xì)小，其熱慣量僅為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)”對貝努最新觀測表明表面顆粒直徑并非決定小行星表面熱慣量唯一決定因素，但表面顆粒直徑越大，小行星熱慣量越大的結(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ù)來源于公式（4）的外推，對于該數(shù)值應(yīng)采取謹(jǐn)慎的態(tài)度。

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

根據(jù)近期“歐西里斯號(hào)”對直徑約為500 m小行星“貝努”（Bennu）[16]的探測結(jié)果可知，其表面無細(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)化層顆粒，這些小顆粒也可能由于溫度交變、太陽風(fēng)與紫外線作用產(chǎn)生的靜電等因素而松動(dòng)，并由于高速自轉(zhuǎn)的離開小行星表面，進(jìn)入空間。因此2016HO3 表面形態(tài)應(yīng)當(dāng)以較大石塊為主，甚至不排除獨(dú)石結(jié)構(gòu)?？梢酝茢啵?016HO3 表面顆粒度大，熱慣量高。

圖5 考慮自轉(zhuǎn)后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 熱慣量的下限。同時(shí)，熱慣量越高，小行星表面溫度變化幅度越低，對附著采樣時(shí)的工程設(shè)計(jì)更友好，因此本文分析中將取 2016HO3 熱慣量下限 700 J·m-2·s-0.5·K-1進(jìn)行后續(xù)分析，確保分析結(jié)果可覆蓋最惡劣情況。

2 小行星的熱建模方法

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

2.1 近地小行星熱模型

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

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

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

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

表2給出了式（5）中關(guān)鍵參數(shù)范圍判斷結(jié)果。其中S⊙r-2代表了太陽能量的輸入，這里分別分析了近地點(diǎn)和遠(yuǎn)地點(diǎn)之間的差異。由于2016HO3 幾何反照率A與紅外發(fā)射率ε無法獲知，這里按照S 類小行星典型數(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 表面被太陽直射點(diǎn)的溫度，見表3所示。而太陽直射點(diǎn)的溫度代表了2016HO3 溫度的上限值，非太陽直射點(diǎn)處的溫度可以用

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

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

2.2 小行星熱物理模型

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

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

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

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

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

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

邊界條件為

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

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

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

上述方程簡化為

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

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

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

2.3 兩種模型的比較

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

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

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

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

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

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

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

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

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

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

假設(shè)2016HO3 是球狀，自轉(zhuǎn)軸方向是影響小行星表面溫度分布的關(guān)鍵參數(shù)。但是由于自轉(zhuǎn)軸數(shù)據(jù)的缺乏，這里研究小行星位于1 AU 處，自轉(zhuǎn)軸相對于公轉(zhuǎn)軌道面不同傾角下的小行星2016HO表面的溫度分析，具體分布見圖7。這里小行星經(jīng)度與緯度采用小行星本體坐標(biāo)系，由于缺乏小行星的形狀信息，假設(shè)小行星是球體，圖7是太陽直射180°經(jīng)度時(shí)小行星表面全球瞬間溫度分布；由于球體的對稱性，圖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è)定太陽直射180°經(jīng)度處，APM 模型計(jì)算出的表面溫度最高點(diǎn)將出現(xiàn)在經(jīng)度約210°處，即考慮熱慣量后，小行星溫度最高點(diǎn)相對太陽直射點(diǎn)有一定延遲。另外，小行星自轉(zhuǎn)軸傾角越接近90°，則太陽直射點(diǎn)溫度越高，最高甚至可以接近特征溫度393 K，這較NEATM的分析結(jié)果更大，而小行星自轉(zhuǎn)軸傾角越接近0°（或者180°），則最高溫度將降低至310 K，因此自轉(zhuǎn)軸對溫度上線影響顯著。同時(shí)需要注意，圖7中極夜區(qū)溫度數(shù)值并無使用價(jià)值，TPM模型無法計(jì)算極夜區(qū)溫度分布。

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

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

公轉(zhuǎn)軌道位置同樣對2016HO3 溫度影響顯著，由于2016HO3 的偏心率到達(dá)0.1，近日點(diǎn)和遠(yuǎn)日點(diǎn)太陽輻射能量差距達(dá)到了40%，這將顯著影響小行星的最高溫度，具體計(jì)算結(jié)果見表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)軸對溫度的影響，而由于沒有使用自轉(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ù)η對自轉(zhuǎn)的效應(yīng)進(jìn)行了一定修正，而自轉(zhuǎn)軸傾角90°時(shí)極區(qū)全天受照，等同于無自轉(zhuǎn)，發(fā)射參數(shù)η在這種條件下又會(huì)降低對最高溫度的預(yù)期。

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

上述計(jì)算僅適用于有能量輸入的區(qū)域，即晝夜交替或者極晝的區(qū)域，而不適用于極夜區(qū)。不考慮橫向熱傳導(dǎo)時(shí)，極夜區(qū)無外部能量輸入，則達(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è)量級；由于自轉(zhuǎn)周期差別，熱特征長度上也相差2個(gè)數(shù)量級[見公式（7）]，因此自轉(zhuǎn)和公轉(zhuǎn)時(shí)熱過程分析在時(shí)間、空間尺度上均有較大的差異，如果放到一塊處理則由于時(shí)間、深度尺度不一，導(dǎo)致計(jì)算網(wǎng)格數(shù)量增加6個(gè)數(shù)量級，運(yùn)算量過大，給計(jì)算帶來困難。

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

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

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

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

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

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

使用該方法可以解決極夜區(qū)溫度的計(jì)算問題，但是同樣需要注意該方法也存在一定缺陷：對于自轉(zhuǎn)熱分析，熱特征深度僅有厘米量級，單元橫向尺度遠(yuǎn)大于熱特征深度，不考慮橫向熱傳導(dǎo)是合理的；對于極夜區(qū)處理，將日均光照平均后計(jì)算一年溫度變化，熱特征深度將達(dá)到1 m，在這種條件下單元橫向尺度與熱特征深度接近，對于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è)熱模型：近地小行星熱模型與小行星熱物理模型，對2016HO3 表面溫度進(jìn)行初步分析，得到主要結(jié)論：

1） 根據(jù)近地小行星熱模型（NEATM），2016HO3 的表面溫度上限最高不會(huì)超過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ì)超過31 ℃。

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

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

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