基于空间光学成像的4179号小行星表面反射特性研究∗

赵东方 刘 鹏 赵 巍†黄长宁 张宏伟 唐降龙

(1哈尔滨工业大学模式识别与智能系统研究中心哈尔滨150001)(2中国空间技术研究院北京空间机电研究所北京100076)

基于空间光学成像的4179号小行星表面反射特性研究∗

赵东方1刘 鹏1赵 巍1†黄长宁2张宏伟2唐降龙1

(1哈尔滨工业大学模式识别与智能系统研究中心哈尔滨150001)
(2中国空间技术研究院北京空间机电研究所北京100076)

2012年12月13日,嫦娥二号在距地球约700万公里的深空,成功飞掠4179号小行星,获得了最高分辨率优于3 m的系列可见光图像.使用最小二乘拟合方法,对成像相机辐射定标数据进行处理,获取绝对定标系数和相对定标校正矩阵,校正原始小行星图像及定标数据,同时反演成像时刻小行星表面辐射亮度.根据Hapke对Nicodemus反射率定义标准在行星科学应用中的具体描述,获取小行星半球反照率.在R、G、B谱段上,小行星表面平均半球反照率分别为0.2083、0.1269和0.1346.小行星表面半球反照率为0.1566.首次使用基于空间光学成像的方法获得了4179号小行星表面光谱反照率,并获得了4179号小行星的表面反照率分布图,对行星科学研究有应用价值.

小行星:个别:4179 Toutatis,行星和卫星:探测,航天器,技术:图像处理

1 引言

行星表面反射特性是表壤(regolith)矿物类别判定、矿物含量反演、粗糙度反演和地形校正的重要依据.反照率是表征行星表面反射特性的重要参数,是建立行星表面热模型(thermal model of planetary surfaces)所需要的关键参数,反照率也与表壤颗粒种类和含量、行星质量等其它重要物理特征有关,通过观察行星表面反照率的变化可以对小行星表面物理结构进行分析.Li等[1]在对433号小行星Eros的研究中发现,表面反照率的变化趋势和表壤空间风化(space weathering)有很强的相关性.表面反射特性研究被认为是对小行星开展科学研究的重要环节之一.

4179号小行星Toutatis是一颗阿波罗型艾琳达族火星轨道穿越小行星.作为一颗近地小行星,其轨道远日点接近木星轨道,近日点在地球轨道附近.鉴于其公转时与地球在太空尺度距离太近(仅0.006 au,是地月距离的2.3倍),4179号小行星被美国航空航天局列为潜在危险小行星.

以往对4179号小行星表面反射特性的研究,其数据主要来源于地面光度测定和偏振测定,信息载体主要是可见光、紫外和红外波段的电磁波.Spencer等人在1992至1993年间,在世界25个地方对4179号小行星进行V谱段光度测定,获得了太阳相角范围在121◦～0.2◦的光变曲线[2].1998年,Hudson等人结合雷达数据,应用Hapke双向光谱反射模型,模拟4179号小行星光变曲线,将模拟结果和Spencer等人观测的光变曲线比较,计算两者之间的最小均方误差,以此寻找最优的Hapke模型参数[3].此外,Lupishko等人在1992年12月到1993年1月间对4179号小行星进行偏振测定,获得了UBVRI 5个谱段的偏振测定数据,并使用Zellner等人1977年研究的光谱斜率和反照率的经验关系[4],推导出相似相角时各谱段的反照率,首次通过地基偏振测定数据估计出4179号小行星表面平均反照率为0.13,R、G、B谱段反照率分别为0.15、0.13、0.13[5].相比于空间光学成像,使用地基光度测定和偏振测定的天文观测手段,受地球大气散射的影响.

2012年12月13日8时30分9秒(UTC),月球探测器嫦娥二号在距地球约700万公里的深空实施第2个拓展任务,以10.73 km/s的速度成功飞掠4179号小行星Toutatis,飞掠最近距离770 m,飞掠过程获得4.5 GB有效可见光图像数据,最高分辨率达2.8 m,是国际上首次近距离获取4179号小行星可见光图像[6–8].Huang等人最先对嫦娥二号传回的4179号小行星可见光图像进行科学分析,揭示了4179号小行星的基本物理特性、表面结构、内部结构以及可能的起源[6].Zhao等人根据嫦娥二号和4179号小行星的轨道数据分析成像时刻光线条件,对预期成像结果进行仿真,并结合地基雷达观测数据对4179号小行星的空间指向、自转特性、几何构型等进行研究[8–10].这些研究成果让人类多角度了解4179号小行星,对认识小行星的形成与演化都具有重要的科学价值.

对4179号小行星成像的相机是设计用来监视嫦娥二号太阳翼展开状态的彩色CMOS相机.本文利用该相机的辐射定标数据来定量研究4179号小行星表面反射特性.通过相对定标数据处理获得校正矩阵,校正原始小行星图像及定标数据;又经绝对定标数据处理获得绝对定标系数;利用校正后的小行星图像和绝对定标系数反演小行星表面辐射亮度,计算出小行星R、G、B各谱段表面反照率和平均反照率.

2 嫦娥二号太阳翼监视相机辐射定标数据处理与分析

辐射定标是建立辐射量与探测器输出量的数值联系的过程.实验室辐射定标主要包括相对辐射定标和绝对辐射定标,其中,相对辐射定标数据分析用以校正相机传感器各像元之间响应不一致的情况,是图像数据信息复原的基础;绝对辐射定标数据是反演小行星表面辐射亮度和分析表面反射特性的基础.CMOS图像传感器的每个探测像元有其固有光学敏感特性.由于各个探测像元之间存在一定差异,相机对空间均匀景物的响应是不一致的.相对辐射定标就是为了使相机各个像元在相同谱段中相同输入辐射亮度下的数字化输出值相同,需要对输出结果进行校正,从而获得校正因子.将校正因子作用在4179号小行星的可见光图像上,得到均匀性校正后的小行星图像.绝对辐射定标是建立探测器像元的数字化输出与接收的辐射亮度之间关系的过程.确定绝对定标系数,将绝对定标系数作用于相对定标后的小行星图像,反演小行星表面辐射亮度.

辐射定标数据处理与表面反射特性分析过程包括预处理、相对辐射定标、绝对辐射定标、表面反射特性分析4个部分,如图1所示.

图1 表面反射特性分析流程Fig.1 Analysis process of surface reflection characteristics

2.1 辐射定标数据获取

使用1.6 m积分球(标准光源)、光谱辐射计(光谱辐射亮度测量)获取辐射定标数据.测试设备及其摆放布局如图2所示,测试时使用光谱辐射计测量积分球各级次下辐射开口处的光谱辐射亮度值,在积分球辐射开口前架设相机,使积分球辐射充满相机的有效通光口径,并保证积分球内置灯发光不会直接入射至相机,开启积分球至定标所需的最高级次,设定相机积分时间,并采集该级次下的定标图像;将积分球设置为下一级次,重复后两步,直至积分球所有辐射亮度级次的定标数据采集完毕.得到的辐射定标数据,由10个辐射亮度级别下共324张定标图像构成,每幅图像的灰度值均为特定辐射亮度下备份相机的数字化输出.

2.2 辐射定标数据预处理

为了减小辐射定标数据处理和反射特性分析工作的计算量,需对辐射定标测试中获得的定标数据进行预处理.

嫦娥二号太阳翼监视相机的焦面器件是Bayer格式的彩色CMOS图像传感器,图像尺寸为1024 pixel×1024 pixel.成像时刻,由于小行星成像只占焦面器件的局部区域,根据小行星成像区域裁定定标区域,在小行星BMP图像中,小行星成像位置在左上角(99,239),右下角(596,670),每幅图像的定标区域尺寸为498 pixel×432 pixel.后面对定标图像和小行星图像的数据处理工作均在该矩形区域上操作.

图2 定标测试设备Fig.2 Calibration equipments

辐射定标测试中,相同辐射亮度级别下被测相机获得多幅定标图像.为了减小随机误差,需要计算各辐射亮度级别下每个像素点的平均灰度值,计算得到了新的定标数据.计算方法如(1)式所示:

其中,DN(i,j,k)为由原始定标数据计算所得的(i,j)像元k辐射亮度级别的平均灰度值, DNraw(i,j,k,p)为定标数据(i,j)像元k辐射亮度级别第p幅图像的灰度值,nk为第k个辐射亮度级别的定标图像个数.实验得到每个像元10个辐射亮度级别下的平均灰度值,即新的定标数据,以498×432×10的矩阵表示.

2.3 相对辐射定标校正因子获取

嫦娥二号太阳翼监视相机的焦面器件为Bayer格式的传感器.按照谱段将每幅定标图像的定标区域划分为4个子区域,对各区域分别计算相对定标校正因子,再将每个像元的校正因子组合构成校正矩阵.具体做法如下：

(1)选取基准值

对于4个谱段的定标数据,分别做如(2)式的处理,获得某一辐射亮度的基准点.

其中,DNs(k)为由定标数据计算得到的k辐射亮度级别的基准值,DN(i,j,k)为定标数据(i,j)像元k辐射亮度级别的灰度值,nλ为中心波长在λ的谱段的像元个数,实验得到每个谱段10个辐射亮度级别下的基准值,以4×10的矩阵表示.

(2)获取相对定标校正矩阵

像元(i,j)的基准值DNs(i,j,k)与定标数据DN(i,j,k)的关系如(3)式所示：

其中,An(k),···,A1(k),A0(k)为待定系数,DNs(i,j,k)为基准值,DN(i,j,k)为定标数据.将定标数据和基准值代入(3)式,即可获得相机的相对定标校正因子,计算时将DNs(k)扩展为DNs(i,j,k),其中每个像元的基准值相同.本文使用1阶、2阶和3阶最小二乘拟合的方法处理定标数据和基准值,表1给出了R、G、B谱段3种方法计算得到的平均残差平方和.数据表明,1阶最小二乘拟合法偏差较大,3阶最小二乘拟合效果最好,但与2阶最小二乘法相比改善很小.因此,本文选用2阶最小二乘拟合法获取相对定标校正因子.用矩阵[A2,A1,A0]i,j表示,它是一个3×249×216的矩阵.实验得到4个3×249×216的校正矩阵,分别对应R、G和B谱段.

表1 不同阶数最小二乘拟合残差对比Table 1 The residual of least square fitting at di ff erent orders

2.4 校正小行星图像及定标数据

将相对定标数据处理获得的校正因子作用在原始小行星图像和定标数据,可以减少由于器件原因引入的图像误差.校正前和校正后小行星图像灰度值及定标数据灰度值之间的关系如(4)式所示:

其中,[A2(i,j),A1(i,j),A0(i,j)]为相对定标校正矩阵.DN(i,j)和DN′(i,j)根据不同的校正过程有不同的含义,当校正小行星图像时,DN(i,j)为原始小行星灰度值, DN′(i,j)为校正后小行星图像的灰度值;当校正定标数据时,DN(i,j)为10个辐射亮度等级的平均灰度,DN′(i,j)为校正后定标数据的灰度值.图3展示了4179号小行星BMP图像R、G、B谱段的相对辐射定标校正结果.

图3 相对辐射定标校正后的4179号小行星图像Fig.3 Relative radiometric calibrated images of Asteroid 4179

2.5 绝对辐射定标系数获取

绝对定标数据处理目的是获得绝对定标系数,探测器数字化输出与输入辐射亮度之间的函数关系如(5)式所示：

其中,DN(i,j)是相机探测像元的数字化输出,Le(i,j)是目标辐射亮度,A(i,j)是像元光电响应度,DN0(i,j)是探测像元输出的数字化暗电流信号.

[A(i,j),DN0(i,j)]是待定系数,也是像元(i,j)的绝对定标系数,将校正后的定标数据代入(5)式,即可获得相机的绝对定标系数,结果为498×432×2的矩阵,即每个像元均有一组绝对定标系数.图4展示了R、G、B谱段图像的中心点和4个顶点的辐射亮度-灰度值变化曲线,3个谱段共15条曲线,相同谱段的5条线基本符合相同的变化规律,因此图中曲线有重叠的情况,其中R谱段的像元在辐射亮度为58.28 W·m−2·sr−1·µm−1和65.19 W·m−2·sr−1·µm−1时有饱和的情况,饱和数据会影响绝对定标系数的精度.因此,实验中R谱段选用8个低辐射亮度级别的定标数据计算绝对定标系数.

图4 辐射亮度随灰度值变化曲线Fig.4 The curve of radiance varying with gray value

2.6 小行星表面辐射亮度反演

根据(5)式探测器数字化输出与输入辐射亮度之间的函数关系,可以反演成像时刻小行星表面真实辐射亮度:

其中,DN′(i,j)为校正后小行星图像的灰度值,DN0(i,j)为绝对定标系数.

将2.5节获得的绝对定标系数和2.4节校正后的小行星图像代入(6)式,得到成像时刻相机接收小行星反射的辐射亮度,每幅校正后的小行星图像的每个像元都有对应的辐射亮度.成像时刻小行星几何中心与嫦娥二号交会距离约1.32 km[10],查阅欧洲空间局近地天体动态网站,成像时刻4179号小行星距太阳距离是1.0136 au1http://newton.dm.unipi.it/neodys/index.php?pc=0,因此可以忽略小行星反射太阳光在传播过程中的损失,即可以认为此处反演的辐射亮度等同于小行星表面反射太阳光的辐射亮度值.对小行星表面反射特性的分析就是在反演的表面辐射亮度的基础上进行的.

3 4179号小行星表面反射特性分析

1977年,Nicodemus系统地给出了反射率定义的术语标准[11].该标准以双向反射分布函数为基础,推导出反射率(reflectance)和反射率因子(reflectance factor).反射率是反射通量( flux)和入射通量之比,其取值范围在0到1之间;反射率因子则是物体表面反射通量与同样光照、同样几何观测条件下被理想朗伯体反射的辐射通量之比.

Nicodemus将入射辐射分为方向型(directional)入射、锥形(conical)入射和半球型(hemispherical)入射;同样地,反射辐射也可分为方向型反射、锥形反射和半球型反射,如图5所示.当入射或反射辐射的立体角Ω无限小,即在dΩ内时,被称为方向型入射或反射;当入射或反射辐射的立体角Ω满足dΩ＜Ω＜2π时,被称为锥形入射或反射;当入射或反射辐射散布在整个半球空间即Ω=2π时,被称为半球入射或反射.

图5 Nicodemus模型Fig.5 Nicodemus model

当入射辐射为方向型、出射辐射为半球型时,求得的反射率即为方向半球反射率(Directional-Hemispherical Re flectance),也称为半球反照率(Hemispherical Albedo),它是建立行星表面太阳热模型所需的重要参数.4179号小行星表面入射辐射和出射辐射特征符合此模型.

1981年,Hapke在他的双向反射光谱理论[12]中,对方向半球反射率在行星科学中的应用进行了具体描述,Hapke将它表述为(7)式：

其中,Leλ(i,j)是辐射定标反演的小行星表面辐射亮度,Nλ是小行星成像点数,λ是中心波长,红光、绿光、蓝光谱段的中心波长分别为700 nm、550 nm和440 nm.

2013年,黄江川等[7]分析了小行星成像时刻光线条件,如图6所示.小行星表面反射为朗伯体反射,将成像半球假设为一个整体平面,相机光轴与成像半球平面垂直.此时,成像平面太阳光线的入射角即为太阳与嫦娥二号相对小行星的夹角(Solar-Asteroid-Probe,SAP),黄江川等[7]给出小行星成像时刻SAP角为38◦,观测时刻是2012年12月13日08:30(UTC).由于高清成像时间内小行星姿态未发生改变,用于分析表面反射特性的12幅校正后的小行星图像入射角均为38◦[13].

图6 4179号小行星成像时刻光线条件示意图[7]Fig.6 Illustration on light conditions during Asteroid 4179 imaging[7]

为了降低小行星局部区域反照率对整体反射特性的影响,本文选用12幅校正后的高清晰度可见光小行星完整图像作为分析主体.4179号小行星的自转周期和进动周期分别为5.41 d和7.35 d[14],嫦娥二号相对4179号小行星观测时间为100 s,其中高清晰度观测时间小于20 s,观测时间远小于自转周期和进动周期,可以认为在嫦娥二号飞掠4179号小行星的过程中小行星姿态未发生改变.

成像时刻4179号小行星距太阳距离是1.013 6 au.查阅太阳光学特性表格[15]得知,大气团airmass=0(不考虑大气因素)时,1 au距离下的太阳光谱辐射照度.根据辐射照度的平方反比定律,可以计算成像时刻4179号小行星表面接收的太阳光谱辐照度,即距太阳1.013 6 au时的太阳光谱辐照度.表2展示了距太阳1 au和1.013 6 au距离的太阳光谱辐照度.

表2 太阳光谱辐照度(单位:W·m−2·µm−1)Table 2 The solar spectral irradiance(unit:W·m−2·µm−1)

观测图像的饱和区域未参与反照率的计算,选取的12幅图像中有效像元数列于表3.

表3 实验图像有效像元数Table 3 Number of effective pixels

图7显示了基于空间可见光图像计算得到的小行星表面反照率分布图.其中,灰度值越高的区域,反照率越高,饱和区域未参与计算,反照率用0表示.本文所获得的基于空间可见光图像的小行星表面R、G、B谱段反照率分布图,具有数据精度高、数据量充足的特点,可以为进一步分析小行星表面的局部光学特性以及形态学特征等提供有力支持.本文基于空间光学成像方法以及Lupishko等[5]基于光变曲线得到的R、G、B谱段和可见光谱段平均方向半球反照率列于表4.本文方法与文献[5]得到的各谱段反照率的关系是一致的,均表现为4179号小行星对红光的反射能力更强.

图7 小行星表面反照率分布Fig.7 Asteroid surface albedo distribution of Toutatis

表4 4179号小行星反照率Table 4 Albedo of Asteroid 4179

4 结论

用1.6 m积分球和光谱辐射计对嫦娥二号太阳翼监视相机备份产品进行辐射定标,获得了相机定标数据.采用最小二乘拟合的方法,获得相对定标校正矩阵,校正了原始小行星图像及定标数据.在绝对定标数据处理中,使用校正后的定标数据,获得了绝对定标系数.对嫦娥二号飞掠4179号小行星任务获取的可见光观测图像进行处理和分析,反演了小行星表面辐射亮度.最后,根据小行星成像时刻太阳光谱辐照度,利用Nicodemus对半球反照率的定义和Hapke对它在行星科学应用中的具体推导,对4179号小行星表面反射特性进行分析,获得了R、G、B谱段小行星表面反照率分布图,以及表面平均半球反照率.结果表明:4179号小行星表面对红光反射能力较强.本文研究结果为进一步分析4179号小行星表面局部光学特性、表面形态学特征以及建立行星表面热模型等研究提供科学依据.

[1]Li J Y,Micheal F A,Lucy A M.Icar,2004,172:415

[2]Spencer J R,Leonid A A,Claudia A.Icar,1995,117:71

[3]Hudson R S,Ostro S J.Icar,1998,135:451

[4]Zellner B,Leake M,Lebertre T,et al.Lunar and Planetary Science Conference Proceedings.New York:Pergamon Press,1977,1:1091

[5]Lupishko D F,Valilyev S V,E fimov J S,et al.Icar,1995,113:200

[6]Huang J C,Ji J H,Ye P J,et al.NatSR,2013,3:3411

[7]黄江川,王晓磊,孟林智,等.中国科学E辑,2013,43:596

[8]赵玉晖,王素,胡寿村,等.天文学报,2013,54:447

[9]Zhao Y H,Wang S,Hu S C,et al.ChA&A,2014,38:163

[10]Zhao Y H,Ji J H,Huang J C,et al.MNRAS,2015,450:3620

[11]Nicodemus F E,Richmond J C,Hsia J J,et al.Geometrical Considerations and Nomenclature for Re flectance.National Bureau of Standards Monograph,1977,9:12

[12]Hapke B.JGR,1981,3039:3054

[13]Bu Y L,Tang G S,Di K C,et al.AJ,2015,149:1

[14]Hudson R S,Ostro S J.Science,1995,270:84

[15]姚连兴.目标和环境的光学特性.北京:中国宇航出版社,1995:61-63

Re flectance of Asteroid 4179 Toutatis Based on Space Optical Image

ZHAO Dong-fang1LIU Peng1ZHAO Wei1HUANG Chang-ning2ZHANG Hong-wei2TANG Xiang-long1
(1 Pattern Recognition and Intelligence System Research Center,Harbin Institute of Technology,Harbin 150001)
(2 Beijing Institute of Space Mechanics Electricity,China Academy of Space Technology, Beijing 100076)

On 2012 December 13,Chang’e-2 probe made a success flyby for Asteroid 4179 Toutatis in deep space of about 7 million kilometers away from the earth,and acquired a series of optical images with high resolution better than 3 m.In this paper,we process the radiation calibration data of imaging camera by least square fitting method, to obtain the absolute calibration coefficient and relative calibration correction matrix, and to recover original intensity of asteroid and its real surface radiance.According to the Nicodemus’reflectance de finition proposed by Hapke,the directional-hemispherical reflectance of Toutatis is obtained.The average surface albedo in R,G,and B spectrum bands are 0.2083,0.1269,and 0.1346,respectively,and the asteroid’s surface albedo is 0.1566.Data indicate that,Toutatis is,somewhat,a red body in visible spectrum.

minor planets,asteroids:individual:4179 Toutatis,planets and satellites: detection,space vehicles,techniques:image processing

P185;

:A

10.15940/j.cnki.0001-5245.2016.01.005

2015-04-14收到原稿,2015-09-02收到修改稿

∗航天科技集团公司重点创新基金资助

†zhaowei@hit.edu.cn