近地小行星基本物理性质的测光研究

王晓彬，黄佳宁，王 翱，刘怡思

(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)小行星相位曲线的形状与其表面物质组成和性质有关。通过反演相位函数参数，可以推断小行星表面的物质组成及性质。但是，大太阳相位角跨度的小行星光度测量值往往是在不同可视期获得的，测光数据不可避免地含有非球形形状的影响，这种影响随视界变化而不同。首先进行小行星形状反演，依据得到的小行星自转轴指向及其形状，可以扣除由非球形形状引起的小行星光度变化，从而提高相位函数参数值的测定精度。