基于序列影像的小天体三维形状重建方法研究

蓝朝桢,耿迅,徐青,崔平远

(1.信息工程大学地理空间信息学院,郑州450052;2.北京理工大学深空探测技术研究所,北京100081)

基于序列影像的小天体三维形状重建方法研究

蓝朝桢1,2,耿迅1,徐青1,崔平远2

(1.信息工程大学地理空间信息学院,郑州450052;2.北京理工大学深空探测技术研究所,北京100081)

提出了一种基于序列影像的小天体三维形状重建方法。利用小天体序列影像间的内在几何约束关系,在无精确位置、姿态数据的条件下,基于稀疏光束法平差求解小天体序列影像的相对位置、姿态。综合利用核线几何约束、半全局匹配、多视最小二乘匹配等策略进行小天体序列影像的密集匹配,重建小天体三维形状信息。由于小天体序列影像尺度变化大、纹理贫乏,导致影像匹配的误匹配率较高,在匹配过程中采用随机KD树搜索策略与拟合透视变换模型的RANSAC算法剔除粗差。利用“黎明号”探测器的绕飞段实测影像,对灶神星(VESTA)的三维形状重建进行了实验,结果证明了该方法的有效性。

小天体;三维重建;深空探测;密集匹配;序列影像

0 引言

小天体探测是人类深空探测最复杂的任务之一,也是目前国际上航天活动的热点。小天体三维形状的准确获取是实现探测器着陆过程导航制导、安全着陆的前提,也是小天体质量、体积、转动惯量、惯性主轴以及引力场模型精确确定的依据。因此,获取小天体的三维形状是各个小天体探测任务的重要目标。由于成像传感器具有范围大、功耗低、分辨率高等特点,被广泛应用于深空探测任务中,利用图像数据重建小天体三维形状是目前普遍采用的方法。但是,由于摄影条件不严格,影像尺度变化大、重叠率较为随意,影像基本处于无序的状态,且影像间的光照条件不统一,辐射变化大。此外,小天体特殊的表面形貌造成影像上纹理稀疏,信息量小。这些因素给小天体三维形状的重建造成很大困难。因此,利用小天体序列影像重建小天体三维形状一直是深空探测领域的一个研究难点。

国内外学者对序列影像的小天体形状重建技术开展了卓有成效的研究和实践。2000年,“NEAR号”探测任务的图像处理小组通过识别估计Eros小天体表面的弹坑分布来恢复小天体的三维模型[1]。国外学者对“黎明号”(DAWN)探测器获取的灶神星(VESTA)的三维形貌提取开展了研究[23],但所采用的方法依赖于高精度的位置和姿态测量数据。国内学者近年来也在此方面开展了相关研究,邵巍(2009)[4]和崔平远等(2010)[5]提出了一种绕飞过程中小天体三维模型重构方法,采用PCA—SIFT进行特征点降维处理和弱透视模型解算,这样虽然提高了速度,但会损失地形重建的精度。

本文提出了一种基于序列影像的小天体三维形状自动重建方法。在连接点提取方面,利用SIFT (Scale-invariant feature transform,尺度不变特征转换)匹配同名点,结合随机KD树搜索策略与RANSAC(Random sample consensus,随机抽样一致算法)剔除粗差,实现连接点的可靠自动匹配;在无影像精确位置、姿态数据的条件下,采用稀疏光束法平差(Sparse bundle adjustment)方法计算序列影像的相对位置、姿态;最后综合利用核线几何约束、多视最小二乘匹配、半全局匹配等策略进行密集匹配,重建小天体三维形状信息。

1 小天体影像连接点自动提取

序列影像间的同名连接点的自动提取与可靠匹配,是探测器相对位置和姿态估计的前提,是利用二维影像重建小天体三维形状的第一步。

1.1 尺度不变影像特征提取

由于小天体探测器受到多种摄动力干扰,加上小天体本身的自转、形状不规则等特性,探测器获取的小天体影像通常会发生旋转、平移、缩放以及光照变化等,要实现影像间可靠的匹配,连接点的提取必须能克服尺度、光照变化等不利因素的影响。

SIFT是计算机视觉领域的经典算法[6]。该算法在空间尺度中寻找极值点,并提取出其位置、尺度和旋转不变量。SIFT算法对旋转、尺度缩放、亮度变化保持不变性,并对影像的仿射变换、视角变化、噪声等影响也有一定的稳定性。本文采用该算法,对小天体影像的特征进行提取和描述。为保证提取的特征在整张影像上均匀分布,采用影像金字塔逐级提取,并利用GPU加速算法,提高特征提取的运算速度。

图1为“黎明号”探测器2011年在2 420 km测绘轨道(Survey orbit)获取的VESTA小行星的连续4帧影像(分辨率为1 024×1 024)。从图中可以看出,小行星表面影像信息量较为贫乏,且由于小行星的自转和探测器的运动,造成地物形状和光照变化较大。图2为利用SIFT算法提取的特征点,图中可以看出,虽然尺度变化较大,但即使在纹理细节不丰富的平坦区域,还是能提取出大量的特征点(实验中,每帧影像提取特征点数量约3 000个),这为下一步特征的匹配和相机位置姿态的估计奠定了基础。

图1 VESTA小行星原始影像Fig.1 Raw sequence images of VESTA

图2 VESTA影像特征提取结果Fig.2 Extracted feature points from VESTA

1.2 随机KD树搜索策略匹配连接点

在数量众多的序列小天体影像上大量特征点高维空间中确定同名特征,效率较低。SIFT的特征描述符为128维,虽然可以通过PCA(Principal component analysis)算法对描述符做降维处理,但会损失特征描述符的一些信息,会导致信息原本就相对贫乏的小天体影像特征匹配的准确性下降。因此,为加快特征匹配的速度,本文引入随机KD树搜索[7](Randomized KD-tree algorithm)策略实现在全维特征空间的特征匹配。随机KD树收集多个树的子节点对各个类别的投票,然后选择获得最多投票的类别作为最终的判断结果,通过计算树的所有子节点的平均值来解决回归问题。随机树建立时的基本子系统是决策树。随机KD树算法便于并行计算,因此容易提升计算性能。对于特征向量间匹配的测度,采用目前基于尺度不变特征匹配最常用的匹配测度——欧氏距离测度欧氏距离最小等价于N维空间点Y与点X之间的距离最小。

1.3 RANSAC误匹配点剔除

采用第1.2节方法得到的匹配点只是粗匹配的结果,其中含有误匹配点,特别是针对纹理信息较少重复特征较多的小天体影像,误匹配率非常高,这对后续的相机位置姿态的求解是非常不利的。本文采用随机抽样一致性算法(Random sample consensus,RANSAC)剔除粗差[8]。RANSAC算法的基本思想是随机选择两个点,这两个点确定了一条直线,在这条直线一定距离范围内的点称为这条直线的支撑。随机选择重复数次,具有最大支撑特征集的直线被确认为是样本点集的拟合。在拟合的误差距离范围内的点称为内点,反之则为外点。根据算法描述可知,随机选取的包含外点的初始点集确定的直线不会获得很多特征点支撑,如图3所示,点A作为外点被剔除。拟合直线只是一种简单的变换模型,对于小天体影像而言,拟合仿射或者透视变换模型更合适,其基本原理与直线的拟合是一样的,本文在使用RANSAC算法剔除误匹配点时使用透视变换作为其拟合的模型。

图3 RANSAC剔除粗差点示意图Fig.3 Outliers removal using RANSAC

对第1.1节中提取的VESTA影像SIFT特征进行随机KD树策略匹配,并采用拟合透视变换模型的RANSAC对误匹配点进行剔除。表1为三幅影像粗匹配点以及正确匹配点的数量统计结果。图4为剔除误匹配后的结果。

表1 粗匹配和粗差剔除结果Table 1 Initial image matching and outliers removal

从实验结果可以看出,由于小天体表面形态的特点,存在相当高的误匹配点,必须进行剔除。经过RANSAC算法剔除误匹配点后,特征点的匹配可靠性得到了极大的提高。这些可靠匹配的同名像点,将用于下一步估计相机位置和姿态参数。

图4 特征点匹配结果Fig.4 Feature point matching results

2 稀疏光束法平差求解相机位置姿态

小天体序列影像位置和姿态的解算,是基于图像重建小天体三维形状的关键。文献[2- 3]的小天体三维重建方法,需要利用初始的位置和姿态数据。本文在无影像初始位置和姿态数据的条件下,利用序列影像内在的几何约束关系,采用稀疏光束法平差[9]方法求解影像相对位置和姿态参数。

2.1 SBA算法原理

稀疏光束法平差是计算机视觉领域在光束法平差的基础上发展的一种利用Levenberg—Marquardt算法求解相机参数和三维点坐标的方法。

可以用7个参数来描述一个相机成像模型,分别是3个旋转参数、3个平移参数和1个相机焦距,如图5所示。

从图中可以看出,每一对匹配连接特征点都对应一个三维点,稀疏光束法平差的误差方程为所有三维点在对应影像上的像点误差的平方和,即

其中,I代表所有影像,x(i)为影像i所对应的三维点列表。设点Xj在影像i上的重投影像点为uij,所对应的特征点为mij,则像点误差为

采用Levenberg-Marquardt算法求解的迭代方程为

式中,Φ是相机参数和特征点坐标的向量;r是残差向量;σ为3倍的中误差;J=∂r/∂Φ是一个M×N的矩阵;M为影像个数;N是待求解未知数的个数;Cp为协方差矩阵。

图5 影像特征点与三维点关系示意图Fig.5 Relationship between image point and 3D point

由于矩阵J大部分元素是0(只有当是对应影像或特征点的参数时导数才不为0),所以可以直接求解JTJ中的元素。由于矩阵JTJ是一个稀疏阵,稀疏矩阵降低了相机参数与特征点坐标之间的相关性。利用稀疏光束法平差先求解相机参数,然后再求解特征点相对坐标,同时得出摄影时相机的相对平移位置和姿态数据。

2.2 稀疏光束法平差策略

为保证平差精度,本文针对小天体影像的特点,采用如下策略:

1)平差初始影像对选择具有较多连接点的两张影像,以该像对的左片像空间坐标系为整个小天体的局部坐标系,通过相对定向和前方交会方法计算该像对相对的位置、姿态以及连接点在上述局部坐标系中的位置初值,并利用稀疏光束法平差方法对这些参数进行最优化求解。

2)添加一幅新影像,利用已求得的三维点坐标以及相应影像之间的连接点计算影像的相机参数,并利用前方交会计算对应特征点的相对坐标,再进行光束法平差优化。

3)继续添加其他影像,直到所有影像都完成,如果一幅影像与其他任一幅影像的匹配特征点太少,则该影像不处理。

3 多视密集匹配小天体三维形状重建

在求解序列影像相对位置、姿态参数的基础上,小天体三维形状重建的关键步骤是序列影像的密集匹配。但是由于小天体影像纹理信息相对贫乏,大量的陨石坑内光照角度不一致,采用传统的匹配方法会导致大量的误匹配点。为达到稳健的密集匹配结果,必须增加相应的约束条件。本文将半全局匹配算法(Semi-global matching,SGM)[10]引入到小天体影像的密集匹配过程,综合采用核线几何约束等方法增加匹配的稳健性,并利用影像金字塔分层匹配提高匹配效率。在此基础上,利用多视最小二乘匹配实现同名像点的高精度匹配,最终重建小天体的三维形状。具体流程如图6所示。

图6 密集匹配重建小天体3D形状流程Fig.6 Process of 3D shape reconstruction

半全局匹配算法基于互相关信息的逐像素匹配思想,通过多个一维约束构成近似全局的二维平滑约束。设参考影像像素p灰度为Ibp,对应待匹配影像的同名点q灰度为Imq。函数q=ebm(p,d)表示匹配影像上对应于参考影像像素p的核线,核线参数是d。对于核线影像,匹配影像位于参考影像右侧, d为视差,则有如下公式

匹配时一个重要的因素,是待匹配区域的尺寸与形状。待匹配区域越大,匹配的稳健性越高。匹配代价函数计算基于互相关信息(MI),互相关信息对亮度变化不敏感。其定义是从两幅图像的熵H(即信息含量)出发

熵由相应图像灰度的概率分布P计算得出

本文利用上述密集匹配方法和策略,对图1中4帧影像进行密集匹配,图7为得到的深度图,图8为4帧影像地形提取结果。

图7 密集匹配后得到的深度图Fig.7 Disparity map after dense image matching

图8 四帧影像地形重建结果Fig.8 Reconstruction results using four images

从图7可以看出,在相机位置姿态估计的基础上,利用密集匹配方法精细地重建了重叠区域的深度图,从最终构建的三维地形(图8)来看,较小的陨石坑也得到了重建。

4 实验结果与分析

为验证上述方法的有效性,从“黎明号”探测器在测绘轨道(Survey orbit)获取的VESTA小行星的1 000多帧影像中,筛选了80张影像(部分影像如图9所示),进行VESTA小行星三维重建试验。

从原始的影像中不难看出,拍照的角度和光照变化都比较大,影像基本处于无序状态。原始影像经过SIFT特征提取匹配以及误匹配点剔除后得到连接点,经过稀疏光束法平差恢复连接点的三维位置如图10所示。图中包含了53 255个连接点。同时还恢复了相机的相对位置和姿态,相对位置的平差精度小于0.4个像素。

利用这些连接点,已经基本重建了VESTA三维形状,但是缺乏表面细节,许多地方还有较大的空洞,因此,还需要利用多帧重叠影像的密集匹配,得到稠密点云。图11为经过多视密集匹配后得到稠密的点云,并经过三角化得到VESTA三维形状重建的最终结果。

局部区域(如陨石坑)的三维重建结果如图12所示,该陨石坑为VESTA上最大的陨石坑Marcia,其表面细节地形得以精确重建。

小天体表面地形精确重建后,可以制作小天体表面的正射影像图,图13为VESTA小行星南极区域的正射影像图。

图9 VESTA序列影像Fig.9 Sequence images of VESTA

图10 连接点三维重建结果Fig.10 3D shape reconstruction using tie points

图11 VESTA三维形状重建的结果Fig.11 3D shape reconstruction results of VESTA

图12 Marcia陨石坑三维重建结果Fig.12 3D shape reconstruction result of rarcia crater

图13 VESTA小行星南极区域的正射影像图Fig.13 Orthophoto in south pole of VESTA

5 结论

小天体的三维重建是小天体探测的重要内容,也是小天体着陆探测的基础。本文针对小天体影像的特点,在无初始位置和姿态的条件下,研究了小天体相对形状的重建问题,提出了一种基于序列影像的小天体三维形状重建方法,并利用“黎明号”探测器的绕飞段实测影像,对灶神星(VESTA)的三维形状重建进行了实验。下一步,拟结合SPICE库提供的初始位置和姿态数据对VESTA小行星三维重建结果进行精化。

[1] Cheng Y,Miller J K.Autonomous landmark based spacecraft navigation system[C]∥13th AAS/AIAA Spaceflight Mechanics Meeting.Ponce,Puerto Rico:[s. n.],2003.

[2] Kersten K D,Matz F,Preusker F,et al.High resolution Vesta HAMO atlas derived from DAWN FC images[C]∥43rdLunar and Planetary Science Conference.The Woodlands,TX,USA:[s.n.],2012.

[3] Preusker F,Scholten F,Matz K D.Topography of Vesta from DAWN FC stereo images[C]∥43rd Lunar and Planetary Science Conference,The Woodlands.TX,USA: [s.n.],2012.

[4] 邵巍.基于图像信息的小天体参数估计及探测器自主导航研究[D].哈尔滨:哈尔滨工业大学,2009.[Shao W.Image based small body parameters estimation and spacecraft autonomous navigation[D].Harbin:Harbin Institute of Technology,2009.]

[5] 崔平远,邵巍,崔祜涛.绕飞过程中小天体三维模型重构及探测器运动估计研究[J].宇航学报,2010,31(5):1381- 1389. [Cui P Y,Shao W,Cui H T.3-D Small body model reconstruction and spacecraft motion estimation during flyaround[J].Journal of Astronautics,2010,31(5):1381 -1389.]

[6] Lowe D G.Distinctive image features from scale invariant key points[J].International Journal of Computer Vision, 2004,60(2):91- 110.

[7] Silpa A C,Hartley R.Optimised KD-trees for fast image descriptor matching[C]∥IEEE Conference on Computer Vision and Pattern Recognition.[s.l.]:IEEE,2008.

[8] Fischler M A,Bolles R C.Random sample consensus:a paradigm for model fitting with application to image analysis and automated cartography[J].Communication Association Machine,1981,24(6):381- 395.

[9] Lourakis M,Argyros A.SBA:A software package for generic sparse bundle adjustment[J].ACM Transactions on Mathematical Software,2009,36(1):1- 30.

[10] Hirschmüller H.Stereo processing by semi-global matching and mutual information[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,2008,30(2):328- 341.

通信地址:郑州市陇海中路66号测绘学院航测教研室(450052)

电话:(0371)81635393

E-mail:lan_cz@163.com

[责任编辑:高莎]

3D Shape Reconstruction for Small Celestial Body Based on Sequence Images

LAN Chaozhen1,2,GENG Xun1,XU Qing1,CUI Pingyuan2
(1.Institute of Surveying and Mapping,Information Engineering University,Zhengzhou 450052,China;2.Institute of Deep Space Exploration Technology,Beijing Institute of Technology,Beijing 100081,China)

A method for 3D shape reconstruction of small celestial body based on sequence images is proposed. Using the internal geometric constraints of sequence images,the relative position and attitude parameters of sequence images are estimated based on sparse bundle adjustment.In order to reconstruct the 3D shape of the small celestial body,such solutions as epipolar line geometric constraints,semi-global matching strategy and multi-view least squares matching are combined to perform dense image matching.In order to solve the mis—matching problem caused by the different scale and lack of texture details,the search strategy of random KD tree and RANSAC method with perspective model are used to detect outliers.Then,the 3D shape of VESTA is reconstructed using the method above by sequence images obtained from the Dawn mission.The results demonstrate our method's effectiveness.

small celestial body;3D reconstruction;deep space exploration;dense matching;sequence images

P234

:A

:2095-7777(2014)02-0140-06

蓝朝桢(1979—),男,讲师,博士,主要研究方向:数字摄影测量,三维可视化。

2013-11-01;

2013-12-30

国家重点基础研究发展计划(973计划)(2012CB720000);中国博士后基金资助(2013M540870)