近地小行星基本物理性質(zhì)的測(cè)光研究

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

(1.中國科學(xué)院 云南天文臺(tái)，云南 昆明 650216；2.中國科學(xué)院大學(xué) 空間科學(xué)學(xué)院，北京 100049；3.中國科學(xué)院 天體結(jié)構(gòu)與演化重點(diǎn)實(shí)驗(yàn)室，云南 昆明 650216；4.云溪師范學(xué)院，云南 玉溪 653100)

1 引 言

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

近地小天體(近日點(diǎn)距離小于1.3天文單位的小行星或彗星)越來越受到人們的關(guān)注，源于這類天體對(duì)地球的碰撞威脅；同時(shí)，由于其距離地球最近，正成為人類未來太空采礦的目標(biāo)。為此，國內(nèi)外多個(gè)研究機(jī)構(gòu)開展了近地小行星巡天觀測(cè)研究項(xiàng)目。隨著觀測(cè)技術(shù)的發(fā)展和在天文觀測(cè)領(lǐng)域的應(yīng)用，自2000年以來近發(fā)小行星的發(fā)現(xiàn)數(shù)目迅速增長。迄今，已發(fā)現(xiàn)近地小天體22 253個(gè)，其中90%的近地小行星的直徑小于10 km。

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

本文旨在通過對(duì)所選小行星測(cè)光數(shù)據(jù)的分析，介紹團(tuán)隊(duì)所建立的小行星物理參數(shù)分級(jí)分析平臺(tái)。

2 小行星光度模型

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

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

眾所周知，小行星的亮度是其表面反射太陽光的結(jié)果?，F(xiàn)在也了解觀測(cè)得到的小行星亮度與小行星離太陽和觀測(cè)者的距離，小行星大小，觀測(cè)幾何、小行星表面散射太陽光的能力——反照率有關(guān)。當(dāng)以上任何一個(gè)因素發(fā)生變化時(shí)，小行星的亮度都會(huì)發(fā)生變化。對(duì)小行星光度及其變化的理解是一個(gè)漫長的歷史過過。

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

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

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

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

圖2 近地小行星星等分布

為了實(shí)現(xiàn)對(duì)近地小行星基本物理參數(shù)的擴(kuò)充與精確測(cè)定，云南天文臺(tái)太陽系小天體團(tuán)隊(duì)利用我國中、小型望遠(yuǎn)鏡開展對(duì)近地小行星的觀測(cè)，并建立了適用于近地小行星基本物理參數(shù)分級(jí)分析的工作平臺(tái)。平臺(tái)整合了近年來本研究團(tuán)隊(duì)自主研發(fā)、與赫爾辛基大學(xué)Muinonen研究團(tuán)隊(duì)合作研發(fā)以及國際上已有的光度模型及算法?；谒⒌墓ぷ髌脚_(tái)，開展了包括近地小行星在內(nèi)的小行星測(cè)光數(shù)據(jù)分級(jí)分析研究。本文從最基本的散射律出發(fā)，簡單介紹小行星面解析光度模型，詳細(xì)說明工作平臺(tái)中所涉及的積分光度模型。

2.1 散射函數(shù)

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

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

或者

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

μ0=cosi，μ=cose，

(1)

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

圖3 散射幾何示意圖

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

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

(2)

(3)

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

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

2.2 小行星盤解析光度模型

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

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

(4)

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

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

(5)

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

2.3 小行星積分光度模型

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

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)是小行星距太陽和觀測(cè)者均為1個(gè)天文單位，相位角為零時(shí)的積分亮度，該值與入射光的比值即為幾何反照率p。G稱為斜率因子，其值的大小與表面物質(zhì)的性質(zhì)相關(guān)。在這種假設(shè)下，公式(6)相當(dāng)于小行星光度隨相位角之間的關(guān)系，也稱之為小行星的相位函數(shù)。

實(shí)際上，大多數(shù)小行星的形狀都是非球形的。做為改進(jìn)，將三軸橢球體，甚至凸面體形狀引入到小行星積分光度模型中。如今，小行星的積分光度模型考慮了小行星形狀、空間姿態(tài)、轉(zhuǎn)動(dòng)及表面散射律等因素。將模型值與小行星的觀測(cè)值進(jìn)行比較，可以得到小行星的形狀、自轉(zhuǎn)周期、自轉(zhuǎn)軸在空間的指向及表面物質(zhì)性質(zhì)參數(shù)等信息。本文以實(shí)例展示3個(gè)小行星積分光度模型及其應(yīng)用結(jié)果。

3 小行星測(cè)光分析

本節(jié)以4個(gè)小行星測(cè)光數(shù)據(jù)的反演分析過程來說明本文整合的小行星測(cè)光反演平臺(tái)中包含的3種積分光度模型及應(yīng)用情況。所涉及的小行星見表1。

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

3.1 Lommel-Seeliger橢球模型

小行星的演化是一個(gè)碰撞演化的歷程。如今，大多數(shù)小行星是引力聚積而成的“碎石堆”結(jié)構(gòu)。按照流體靜力學(xué)轉(zhuǎn)動(dòng)平衡理論[23]，小行星的轉(zhuǎn)動(dòng)平衡形狀接近于三軸橢球。橢球體的軸比與小行星的密度及其自轉(zhuǎn)速率相關(guān)。

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

(7)

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

(8)

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

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

(9)

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

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

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

圖5 小行星(155140)橢球軸比及自轉(zhuǎn)參數(shù)分布

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

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

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

圖7 小行星(155140)相位函數(shù)參數(shù)的后驗(yàn)分布

3.2 Cellinoid模型

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

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

(10)

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

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

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

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

3.3 凸面體模型

從以上兩個(gè)例子可以看到，三軸橢球體模型及Cellinoid橢球體模型對(duì)兩個(gè)小行星的光度曲線能很好地?cái)M合，這種情形下所求解的自轉(zhuǎn)參數(shù)，尤其是自轉(zhuǎn)軸指向的精度是有保證的。如果小行星的光變曲呈現(xiàn)更復(fù)雜的不規(guī)則變化，使用以上兩種小行星光度模型很難很可靠地估算小行星物理參量。為此，引入能夠更精確描述小行星形狀的光度模型。

以(103)Hera的測(cè)光數(shù)據(jù)分析為例，下面給出應(yīng)用Lommel-Seeliger橢球體模型與凸面體模型的所得自轉(zhuǎn)參數(shù)結(jié)果的差異[29]。

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

(11)

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

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

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

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

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

表2 小行星(103)Hera的自轉(zhuǎn)參數(shù)

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

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

在利用觀測(cè)量進(jìn)行相關(guān)物理參數(shù)測(cè)定時(shí)都需要回答測(cè)定精度問題。這對(duì)具有潛存威脅的近地小行星的物理參數(shù)的測(cè)定尤為重是。

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

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

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

圖12 自轉(zhuǎn)周期值掃描搜尋結(jié)果

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

圖13 (346)的觀測(cè)值(黑色點(diǎn))及理論值(紅色線)

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

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

Kaasalainen[31]利用設(shè)置不同散射律對(duì)得到自轉(zhuǎn)軸指向的不確定性進(jìn)行評(píng)估，他們認(rèn)為凸面體模型測(cè)定的自轉(zhuǎn)軸精度約在±10°。始終沒有一個(gè)方法對(duì)于所得形狀的不確定性進(jìn)行評(píng)估。

本文利用蒙特卡洛方法建立了一個(gè)新的估算方法。新的方法可以對(duì)由于觀測(cè)誤差引起的自轉(zhuǎn)參數(shù)及形狀不確定性的評(píng)估。

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

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

圖15 小行星(346)的自轉(zhuǎn)參數(shù)聯(lián)合分布

(12)

(a)最佳解

(b)與最佳解比較的解

4 結(jié) 論

由小行星測(cè)光數(shù)據(jù)進(jìn)行的形狀、自轉(zhuǎn)參數(shù)及表面散射參數(shù)的反演已逐漸成熟。關(guān)鍵的問題是如何將這一反演技術(shù)應(yīng)用于近地小行星的基本性質(zhì)的分析上。不同于主帶小行星，暗弱的近地小行星在沖附近時(shí)具有運(yùn)動(dòng)速度快，加上大多數(shù)近地小行星是2000年以后發(fā)現(xiàn)的，測(cè)光數(shù)據(jù)的積累遠(yuǎn)不如主帶小行星多。大部分近地小行星缺乏基本物理參數(shù)信息。

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

(1)應(yīng)用Lommel-Seeliger橢球模型，可以在近地小行星發(fā)現(xiàn)后的短時(shí)間內(nèi)對(duì)其基本物理參數(shù)進(jìn)行初步估算。如果目標(biāo)小行星的形狀或者說光變曲線較規(guī)則，Lommel-Seeliger橢球模型可以給出可靠的形狀和自轉(zhuǎn)參數(shù)的解。

(2)對(duì)于光變曲線不規(guī)則(例如，光變曲線的峰值不對(duì)稱)的目標(biāo)，可用Cellinoid 橢球模型改進(jìn)由三軸橢球體估算的近地小行星形狀和自轉(zhuǎn)參數(shù)的測(cè)定精度。

(3)對(duì)于那些光變曲線形狀極不規(guī)則的近地小行星，需要用凸面體光度模型來分析其基本物理參數(shù)。該光度模型要求測(cè)光數(shù)據(jù)的視界角及相位角的跨度足夠大。利用凸面體模型進(jìn)行小行星光變曲線反演計(jì)算包括兩個(gè)求解參數(shù)過程：由測(cè)光數(shù)據(jù)求解光度模型中的待測(cè)參數(shù)(自轉(zhuǎn)參數(shù)、形狀參數(shù)及散射性質(zhì)參數(shù))；由形狀參數(shù)重構(gòu)凸面體形狀(構(gòu)成閉合凸面體的頂點(diǎn))。復(fù)雜且費(fèi)時(shí)形狀重構(gòu)過程制約了對(duì)小行星凸面體反演解不確定性的估計(jì)。利用本文建立的統(tǒng)計(jì)學(xué)方法，實(shí)現(xiàn)凸面體反演中參數(shù)，包括形狀不確定性的評(píng)估。

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