近地小行星基本物理性質的測光研究

王曉彬，黃佳寧，王 翱，劉怡思

(1.中國科學院 云南天文臺，云南 昆明 650216；2.中國科學院大學 空間科學學院，北京 100049；3.中國科學院 天體結構與演化重點實驗室，云南 昆明 650216；4.云溪師范學院，云南 玉溪 653100)

1 引 言

太陽星云假說[1]是被普遍接受的太陽系形成理論。太陽系小天體，特別是小行星帶的存在為這一理論提供了重要的佐證。主帶小行星(指火星與木星軌道之間的區域的小行星，這個區域也稱為小行星帶)碰撞形成理論認為：較早形成的木星對小行星帶的行星子的引力攝動，使得行星子之間的碰撞加劇，抑制了小行星帶內行星的形成。此后，隨著木星軌道向內和向外的遷移，小行星帶中的大量小天體被踢出該區域，被踢向太陽系內層的小行星有的撞向了地球、月球和其他類地行星，有些則被拋向太陽系外層。這也是現在小行星帶小天體總質量只有月球的4%的原因。小行星帶中幾個空隙(也稱，Kirkwood間隙)處于與木星和土星平均運動共振。進入共振帶區域小天體受到大行星的引力攝動作用而改變軌道，來到地球附近。

近地小天體(近日點距離小于1.3天文單位的小行星或彗星)越來越受到人們的關注，源于這類天體對地球的碰撞威脅；同時，由于其距離地球最近，正成為人類未來太空采礦的目標。為此，國內外多個研究機構開展了近地小行星巡天觀測研究項目。隨著觀測技術的發展和在天文觀測領域的應用，自2000年以來近發小行星的發現數目迅速增長。迄今，已發現近地小天體22 253個，其中90%的近地小行星的直徑小于10 km。

面對未來可能發生的近地小行星對地球的碰撞，人們正在考慮多種防災減災措施[2-5]。例如，利用核彈或火箭對入侵目標進行轟擊，通過粉碎目標或改變目標的原運行軌道來避免碰撞。還有的利用重力拖拽方法以及表面反射作用方式使近地小天體的軌道發生偏轉。具體選擇哪種方法以及如何實施都需要了解目標近地小行星的基本物理性質(大小、形狀、質量物質組成和內部結構)。然而，在已發現的近地小天體中僅有很少數的目標有基本物理參量信息。因此，近地小行星基本物理參數的測定是近地小行星碰撞預警工作中重要任務之一。

本文旨在通過對所選小行星測光數據的分析，介紹團隊所建立的小行星物理參數分級分析平臺。

2 小行星光度模型

按照軌道形狀，近地小天體被劃分為4類：阿莫型(Amor)、阿波羅型(Apollo)、阿坦型(Aten)及阿提拉型(Atira)(如圖1所示)。 近地小行星中，近60%的是阿波羅型。 阿波羅型及阿坦型都穿過地球軌道，存在著碰撞地球威脅，尤其半徑大于140 m的近地小天體。

圖1 近地小天體軌道及地球軌道

眾所周知，小行星的亮度是其表面反射太陽光的結果。現在也了解觀測得到的小行星亮度與小行星離太陽和觀測者的距離，小行星大小，觀測幾何、小行星表面散射太陽光的能力——反照率有關。當以上任何一個因素發生變化時，小行星的亮度都會發生變化。對小行星光度及其變化的理解是一個漫長的歷史過過。

對小行星亮度的最早觀測可以追溯到1887年[6]。當時人們觀測到小行星的光度會隨太陽相位角(小行星對光源和觀測者的張角)變化并且還具有幾個小時的周期性變化。早期研究認為小行星是球形的，小行星光度的變化是由觀測者和光源相對于小行星的幾何位置變化(或者說相位角變化)引起。后來，人們逐漸意識到小行星可能是非球形的，這一猜想在1993年被伽里略飛船飛近小行星(243)Ida后所證實。

對于一個旋轉的非球體小行星而言，其可視表面大小會隨自轉而變化，引起小行星視光度按其自轉周期的變化。本文中所說的小行星測光觀測研究就是借助小行星的光度模型來分析小行星的光度觀測值，估算出小行星形狀、自轉、空間姿態、表面散射性質等物理參數。

基于小行星光度模型開展的小行星基本物理參數分析已有近三十多年的歷史。1989年Magnusson等[7]通過小行星的光度變化振幅實現了小行星自轉軸指向的估算，在這種方法中，僅考慮了小行星自轉時視截面大小的變化。1988年Drummond等[8]估算了小行星(21)Lutetia的三軸橢球體形狀。最近，Cellino等[9]應用三軸橢球體模型分析了Gaia巡天觀測中的小行星數據。基于Russell[10]在1906年提出的想法，2001年Kaasalainen等[11-12]實現了小行星凸面體反演。這一方法利用小行星光變曲線(光度隨時間的變化)反演出包裹小行星的凸面體的面元大小及頂點。所得到的小行星凸面體形狀與空間觀測得到的真實形狀很接近。但是這一方法所涉及到的待測參數較多，要求測光數據須覆蓋一定寬度的視界角(觀測者方向與小行星自轉軸的夾角)和相位角。

本文所涉及的研究目標近地小行星，95%的亮度暗于18等(如圖2)。超過60%的近地小行星亮度暗于22等。加上這些暗弱近地小行星的沖附近時運動速度快，使得近地小行星測光數據的獲得面臨較多困難，并阻礙近地小行星物理參數的獲得。

圖2 近地小行星星等分布

為了實現對近地小行星基本物理參數的擴充與精確測定，云南天文臺太陽系小天體團隊利用我國中、小型望遠鏡開展對近地小行星的觀測，并建立了適用于近地小行星基本物理參數分級分析的工作平臺。平臺整合了近年來本研究團隊自主研發、與赫爾辛基大學Muinonen研究團隊合作研發以及國際上已有的光度模型及算法。基于所建立的工作平臺，開展了包括近地小行星在內的小行星測光數據分級分析研究。本文從最基本的散射律出發，簡單介紹小行星面解析光度模型，詳細說明工作平臺中所涉及的積分光度模型。

2.1 散射函數

小行星的亮度是其表面散射太陽光的結果。對于小行星表面上面元(法線方向n)來說，觀測得到的亮度I可用散射律r乘入射流量J來表征。即：

I(i，e，α)=J×r(i，e，α)，

或者

I(μ0，μ，α)=J×r(μ0，μ，α)，

μ0=cosi，μ=cose，

(1)

其中：r(i，e，α)也稱為雙向散射函數/散射律。如圖3所示，i，e分別是入射光源及觀測者視線與面元法線的夾角，α是太陽相位角。

圖3 散射幾何示意圖

理論上，如果知道小行星表面的散射律，就可按光源入射角及出射角計算出小行星表面任一點上的亮度。而實際中，由于人們對小行星表面物質組成缺乏了解且性質信息缺乏，因此一直沒有找到適合于小行星表面的精確散射律。至今，人們多采用一些經驗的散射函數應用于小行星的光度模型中。常用的散射律有： Lambert散射律[13]、Lommel-Seeliger 散射律[14]、Minnaert 散射律[15]、 Lunar-Lambert 散射律[16]、 Hapker模型[17]及 Bowell-Lumme 模型[18-19]。在本文的分析平臺所包含的小行星光度模型中僅涉及Lambert散射律和Lommel-Seeliger 散射律。

Lambert散射律(式(2))適用于表面反照率較高的無大氣天體表面(例如月球表面)； Lommel-Seeliger 散射律(式(3))則能更好地描述了表面反照率低的無大氣天體表面：

(2)

(3)

其中：p(α)是單個粒子散射的相位函數，p(α)=1意味著表面顆粒各向同性散射性質。AL和ALS分別是Lambert和Lommel-Seeliger反照率。Minnaert散射律是Lambert的一種推廣形式，適合描述在有限相位角內各種性質行星表面。組合的Lunar-Lambert函數適用于較高反照率的小天體表面(例如S型、V型和E型小行星)。

Lumme-Bowell和Hapke散射律是光在粗糙的顆粒表面輻射轉移的近似解。 輻射轉移過程考慮了單次散射、多次散射、遮擋影應、沖效應以及單個粒子散射各向異性等因素。詳細公式可參見文獻[20]中式(15)和式(24)。

2.2 小行星盤解析光度模型

Hapke和Lumme-Bowell散射律常被用于描述小行星盤解析光度。兩個模型包含表面物理參數：粒子單次反照率、表面粗糙度、沖效應的幅度和寬度、淺表體密度和不對稱因子。

考慮入射的太陽光為均勻光，Hapke面解析光度分布可寫成公式(4)：

(4)

其中：φ是光入射平面與觀測平面間的夾角；B(α)是沖效應函數；H(μ0/μ)是與多次散射有關的函數[21]；S(μe0，μe，φ)由表面粗糙度引起的遮擋函數；μe0，μe的具體計算公式可參見文獻[22]。

同樣考慮了單次散射和多次散Bowell-Lumme面解析光度模型簡寫成式(5)：

(5)

實際上大多數小行星的光度觀測數據是由地面望遠鏡得到的小行星積分光度，即望遠鏡終端上記錄的是小行星表面被照亮及可見部分散射太陽光的總和。

2.3 小行星積分光度模型

早期人們假設小行星是球形，直接對盤解析光度函數/或散射率進行積分得到小行星積分光度模型。例如，由Lumme-Bowell盤解析光度在球面積分，并加一些假設得到以下經驗的積分光度模型：

I(α)=I(0o)[(1-G)φ1(α)+Gφ2(α)]，

(6)

φi(α)=eAi(tan α/2)Bi，i=1，2，A1=3.33，B1=0.63，A2=1.87，B2=1.22，

其中：I(0o)是小行星距太陽和觀測者均為1個天文單位，相位角為零時的積分亮度，該值與入射光的比值即為幾何反照率p。G稱為斜率因子，其值的大小與表面物質的性質相關。在這種假設下，公式(6)相當于小行星光度隨相位角之間的關系，也稱之為小行星的相位函數。

實際上，大多數小行星的形狀都是非球形的。做為改進，將三軸橢球體，甚至凸面體形狀引入到小行星積分光度模型中。如今，小行星的積分光度模型考慮了小行星形狀、空間姿態、轉動及表面散射律等因素。將模型值與小行星的觀測值進行比較，可以得到小行星的形狀、自轉周期、自轉軸在空間的指向及表面物質性質參數等信息。本文以實例展示3個小行星積分光度模型及其應用結果。

3 小行星測光分析

本節以4個小行星測光數據的反演分析過程來說明本文整合的小行星測光反演平臺中包含的3種積分光度模型及應用情況。所涉及的小行星見表1。

表1 本文分析所包含小行星的軌道

3.1 Lommel-Seeliger橢球模型

小行星的演化是一個碰撞演化的歷程。如今，大多數小行星是引力聚積而成的“碎石堆”結構。按照流體靜力學轉動平衡理論[23]，小行星的轉動平衡形狀接近于三軸橢球。橢球體的軸比與小行星的密度及其自轉速率相關。

假設小行星的形狀為三軸橢球體，Muinonen 等人建立的Lommel-Seeliger 橢球光度模型[25]。對于橢球體表面上單位面元所散射的光度可寫成式(7)：

(7)

將式(7)按橢球體表面進行積分，可得到給定光源方向eΘ及觀測者方向e?小行星積分光度模型的解析式如公式(8)：

(8)

小行星Lommel-Seeliger橢球光度模型涉及12個待測參數：自轉軸指向(λ，β)、自轉周期Per、自轉的初相位角φ0、橢球的三個半長軸(a，b，c)、幾何照率p以及相位函數參數(H，G1，G2)。ΦHG1G2和ΦLS分別是三參數相位函數和Lommel-Seeliger積分相位函數[24]。公式(8)中輔助量S，S?，SΘ，λ′及α′與待測參數間的關系式可參看相關文獻[25]。

應用最小二乘法和馬可夫鏈蒙特卡洛方法(Markov Chain Monte Carlo，MCMC)進行測待參數求解都是基于小行星觀測數據與理論值之間卡方值的計算：

(9)

原則上，可以同時估算模型所涉及的12個參數。但在實際的分析過程中，由于小行星觀測的儀器星等向標準星等系統轉換的精度遠低于小行星相對測量(或者說較差測量)的精度，為了精確測定小行星的形狀和自轉參數，本文通常利用小行星的相對流量擬合小行星形狀(b/a，c/a)和自轉參數，然后再擬合相位函數參數。

利用Lommel-Seeliger橢球光度模型，本文分析了小行星(155140)2005UD的測光數據。小行星(155140)是一個阿波羅型近地小行星。其絕對星等(距離太陽和觀測均為1天文單位、相位角為零時的星等)約17.4等。2018年10月，利用位于西藏阿里的一臺望遠鏡，對該小行星進行了11個夜晚的測光觀測。

圖4 小行星(155140)的光變曲線

圖5 小行星(155140)橢球軸比及自轉參數分布

首先，利用Nelder-Mead downhill算法，得到該小行星的自轉軸指向的黃道坐標為(285.41°，-23.36°)，橢球體的軸比為b/a=0.76和c/a=0.40。圖4中紅線是Lommel-Seeliger 橢球模理論值。應用馬可夫鏈方法蒙特卡洛模擬，得到待測參數的后驗分布(如圖5)。并給出小行星形狀反演的最佳解及誤差情況(彩圖見期刊電子版)。

第二步，利用小行星歸算后的光度值，進行相位函數參數的反演。在小行星(155140)的相位曲線分析中，除了團隊自已觀測的數據外，還包括了ZTF (Zwichky Transient Facility)的數據[26]。同樣地，利用MCMC方法得到小行星相位函數參數的最優解：H=17.19，G1=0.573，和G2=0.004。圖6中紅線即為小行星相位函數對觀測值的最佳擬合結果。圖7給出相位函數參數的后驗分布(彩圖見期刊電子版)。

圖6 小行星(155140)的相位曲線

圖7 小行星(155140)相位函數參數的后驗分布

3.2 Cellinoid模型

當觀測的小行星自轉光變曲線的形狀呈現不對稱的極值時，三軸橢球體模形型很難擬合好這類光變曲線。對于這種情形，本文引入Cellinoid 橢球——由8個八分之一的橢球拼接而來(如圖8(a)所示)。

圖8 (a)Cellinoid 橢球示意，(b)面元分割方式

(10)

分析用到(106)Dione在1981年至2015年之間獲得的26條光變曲線(圖9)。應用MCMC方法得到一對自轉軸的解(58.0°，+21.1°) 和(242.3°，+21.3°)；相應的形狀參數為：b/a=0.9，c/a=0.63，a1/a=1.65，b1/b=0.86，c1/c=1.64。從圖9中可以看出，光度模型很好地擬合了觀測的光變曲線。

圖9 小行星(106)的光變曲線及最佳模型解

然后，取每條光變曲線的平均值(歸算后的星等)再進行H-G1G2相位函數擬合。經過MCMC模擬得到該小行星相位函數的參數的最佳解：H=7.66，G1=0.68，G2=0.08。圖10給出小行星(106)的相位曲線及最佳相位函數擬合值。

圖10 小行星(106)Dione的相位曲線

3.3 凸面體模型

從以上兩個例子可以看到，三軸橢球體模型及Cellinoid橢球體模型對兩個小行星的光度曲線能很好地擬合，這種情形下所求解的自轉參數，尤其是自轉軸指向的精度是有保證的。如果小行星的光變曲呈現更復雜的不規則變化，使用以上兩種小行星光度模型很難很可靠地估算小行星物理參量。為此，引入能夠更精確描述小行星形狀的光度模型。

以(103)Hera的測光數據分析為例，下面給出應用Lommel-Seeliger橢球體模型與凸面體模型的所得自轉參數結果的差異[29]。

小行星凸面體模型是假設小行星的形狀為凸面體，這里同樣用Lommel-Seeliger與Lambert律按權重線性組合的散射率。類似地，凸面體積分光度的也是按對單位球面上三角分割所對應的多面體面元反射太陽光的總和來計算。這里所用的凸面體模型及程序是由Kaasalainen所建立的[30-31]。具體形式如公式(11)：

(11)

其中：eΘ，e?分別是光源和觀測者方向單位矢量。以有限項球諧函數表征的高斯面密度G(?，φ)是用于計算凸面體上法線方向(?，φ)的面元的大小。相位函數f(α)采用了線性函數及指數函數的組合形狀。凸面體模型中所含的待測參數包括：自轉參數(λ，β，per)、形狀參數(球諧函數的系數alm)、相位函數參數(A0，D，κ，c)。形狀參數alm的個數取決于球諧函數展式的截斷級數l和階數m的值。針對現有地面測光觀測的精度，設置l=m=8足夠描述小行星光變曲線中的信號。

由于形狀參數的增加，凸面體光度模型中所含待測參數遠多于前2種模型。求解參數所需要的觀測數據相應也多(需要至少3～4個不同可視期的覆蓋小行星自轉周期的測光觀測)。具體地，Kaasalainen的凸面體反演程序應用了Levenberg-Marquardt最小二乘法來求解所涉參數(詳細可參看相關文獻[30])。

利用Kaasalainen的凸面體模型及Muinonen的Lommel-Seeliger橢球模型，本文分別分析了小行星(103)Hera 的測光數據。數據包含分布于5個可視期的40條光變曲線，數據的相位角α的跨度為1.6°～22.7°，視界角θ的跨度為65°～125°。圖11給出利用兩種方法擬合觀測數據情況，可以看出凸面體模型(橙色實線)可以很好地擬合觀測的光變曲線(黑色散點)，而Lommel-Seeliger橢球模型(黑色虛線)則偏差較大(彩圖見期刊電子版)。

圖11 小行星(103)Hera的光變曲線擬合

比較兩種方法所得到自轉軸結果[27](如表2)，可以看出，自轉周期值符合較好，自轉軸指向在經度方向有9°差異。表2在第3列還例舉了Hanus等人利用凸面體模型分析(103)的結果[32]，與本文分析結果在自轉軸經度上是基本一致，但緯度相差5°。相比于Hanus等人的研究，本文分析使用了更多的測光數據。

表2 小行星(103)Hera的自轉參數

雖然，本文可以比較不同形狀模型得到的自轉軸參數的差異，甚至相同的形狀模形，但使用不同的觀測數據所得到的自轉參數的解的差異。對形狀的不確定性的評估一直是未能很好解決的問題。為此，本文建立了一個新的基于蒙特卡洛方法的評估的方法。

3.4 小行星凸面體反演的精度估算

在利用觀測量進行相關物理參數測定時都需要回答測定精度問題。這對具有潛存威脅的近地小行星的物理參數的測定尤為重是。

以主帶小行星(346)Hermentaria為例，詳細說明測光反演及基本物理參數不確性估算流程。小行星(346)Hermentaria是直徑110 km石質類小行星。早期的測光觀測數據給出幾個不同的自轉周期：28.33 h[30-31]，19.4 h[33]及9.7 h[34]。這主要由長周期光變小行星觀測缺限引起。對于一個站點的望遠鏡來講，每次觀測長周期或者周期值是24 h的整倍數、1/2倍數等小行星時，得到的多是光變曲線片段，以這樣的數據進行周期分析時就會出現以上的結果。小行星光變曲線反演則能精確測定這類小行星自轉周期。

在小行星(346)的凸面體反演分析中總共用到分布在4個可視期：1980，2001，2002及2015的23條光變曲線(如圖13)。數據的相位角跨度為2.4～°20.5°，視界角跨度為36°～155°。

第一步：利用Kaasalainen的凸面體反演程序求解反演的最小二乘解。為找到Levenberg-Marquardt最小二乘全局最優解，對自轉參數空間進行掃描。首先，對自轉周期范圍16～24 h進行掃描，掃描步長為per2/2ΔT， 其中ΔT是測光數據最大的時間跨度。最小的χ2出現在17.79 h附近，次極小則在20.3 h附近(如圖12)。

圖12 自轉周期值掃描搜尋結果

然后，對自轉軸指向在整個球面進行粗的網格掃描(步長90°：90°)。找到最可能的自轉軸指向后，以掃描得到的周期值及軸指向為初值，再進行所有參數的求解。最終得到一對RMS相近的軸指向解(134.4°，16.8°) 和(322.4°，14.7°)，相應的自轉周期值為17.790 00 h和17.790 04 h。圖13是第一軸指向解的光度理論值與觀測的比較(彩圖見期刊電子版)。

圖13 (346)的觀測值(黑色點)及理論值(紅色線)

第二步：利用得到的球諧函數系數(alm) 重構凸面體。這里應用了Minkowski 理論[37]，依據所給的三角分割方式和球諧函數系數(alm)來找到包裹小行星的凸面體的頂點。求解的過程是另一個最小二乘求解過程(詳細可參看文獻[24]的附錄)。圖14是對應于第一軸解重構的(346)的凸面體形狀，紅色箭頭為X軸，綠色箭頭為Y軸，藍色箭頭為Z軸，也是自轉軸方向(彩圖見期刊電子版)。

從小行星光變曲線反演模型參數是一個較快的計算過程，而從球諧函數系數重構凸面體則是一個費時的計算過程(重構(346)一個解大約需要10 min計算機時間)。

Kaasalainen[31]利用設置不同散射律對得到自轉軸指向的不確定性進行評估，他們認為凸面體模型測定的自轉軸精度約在±10°。始終沒有一個方法對于所得形狀的不確定性進行評估。

本文利用蒙特卡洛方法建立了一個新的估算方法。新的方法可以對由于觀測誤差引起的自轉參數及形狀不確定性的評估。

圖14 小行星(346)的凸面體

第三步，在前面兩步計算基礎上，執行蒙特卡洛模擬。具體地，第一步產生一系列虛擬測光數據(在觀測數據上加入高斯噪聲)；第二步求解基于虛擬光變曲線的最小二乘解，得到一系列凸面體反演虛擬解。通常以上過程會重復超過10 000次。所有的虛擬解構成了自轉參數聯合分布可以用于估算自轉參數的最佳解如圖15中紅色線所示，以及誤差范圍如籃色線表標的1-σ的范圍(彩圖見期刊電子版)。

圖15 小行星(346)的自轉參數聯合分布

(12)

(a)最佳解

(b)與最佳解比較的解

4 結 論

由小行星測光數據進行的形狀、自轉參數及表面散射參數的反演已逐漸成熟。關鍵的問題是如何將這一反演技術應用于近地小行星的基本性質的分析上。不同于主帶小行星，暗弱的近地小行星在沖附近時具有運動速度快，加上大多數近地小行星是2000年以后發現的，測光數據的積累遠不如主帶小行星多。大部分近地小行星缺乏基本物理參數信息。

從近地小行星碰撞防御及空間探測的需求上看，更需要了解包括物質組成在內的近地小行星的基本物理參數。據此，本文針對近地小行星的特點及其觀測數據的積累情況，開展近地小行星測光數據分級反演研究。未來可應用不同的光度模型開展近地小行星的基本物理性質的反演研究。

(1)應用Lommel-Seeliger橢球模型，可以在近地小行星發現后的短時間內對其基本物理參數進行初步估算。如果目標小行星的形狀或者說光變曲線較規則，Lommel-Seeliger橢球模型可以給出可靠的形狀和自轉參數的解。

(2)對于光變曲線不規則(例如，光變曲線的峰值不對稱)的目標，可用Cellinoid 橢球模型改進由三軸橢球體估算的近地小行星形狀和自轉參數的測定精度。

(3)對于那些光變曲線形狀極不規則的近地小行星，需要用凸面體光度模型來分析其基本物理參數。該光度模型要求測光數據的視界角及相位角的跨度足夠大。利用凸面體模型進行小行星光變曲線反演計算包括兩個求解參數過程：由測光數據求解光度模型中的待測參數(自轉參數、形狀參數及散射性質參數)；由形狀參數重構凸面體形狀(構成閉合凸面體的頂點)。復雜且費時形狀重構過程制約了對小行星凸面體反演解不確定性的估計。利用本文建立的統計學方法，實現凸面體反演中參數，包括形狀不確定性的評估。

(4)小行星相位曲線的形狀與其表面物質組成和性質有關。通過反演相位函數參數，可以推斷小行星表面的物質組成及性質。但是，大太陽相位角跨度的小行星光度測量值往往是在不同可視期獲得的，測光數據不可避免地含有非球形形狀的影響，這種影響隨視界變化而不同。首先進行小行星形狀反演，依據得到的小行星自轉軸指向及其形狀，可以扣除由非球形形狀引起的小行星光度變化，從而提高相位函數參數值的測定精度。