2. 61618部队, 北京 102102
2. Unit 61618, Beijing 102102, China
在小行星探测过程中,形貌测绘与表征技术作为认识、探测小行星目标的重要内容,是保障小行星探测器自主导航、安全着陆的前提,逐渐成为深空探测领域的研究热点[1-2]。小行星探测器所携带的科学仪器主要有面阵相机、成像光谱仪、激光测距仪等,通过这些仪器可获取大量的影像数据和激光测距数据。激光测距仪能够实时测量探测器与目标小行星之间的距离,结合探测器的位姿信息,可以重建出小行星上目标点的三维坐标,但该方法因获取点数十分有限而不被采用。成像传感器具有范围大、功耗低、分辨率高等特点,被广泛应用于深空探测任务中,利用其影像数据重建小行星的三维地形成为首选方式。但是,摄影条件不严格,影像重叠率较为随意,甚至处于无序状态,加之影像中光照条件变化大,小行星表面色调相对单一,这给小行星三维形状的重建造成很大困难。而利用立体相机同时获取小行星目标区域的光学影像,可以增加影像间的重叠率,有助于提高目标区域三维地形测量的实时性。
本文以立体相机为数据获取设备,设计了一种基于系列影像的小行星三维模型自动重建方法,即利用处于无序状态的系列影像来构建小行星整星的三维形貌模型;提出了一种立体视觉在线测绘小行星局部形貌的方法,即在探测器抵近过程中实时测量探测区域的三维地形。在已知目标区域地形信息的情况下,提出了一种基于等值线分析的撞击坑特征提取方法,可以从获取的三维地形中识别出撞击坑特征,进而提高了探测器对小行星表面形貌的感知能力,提高了小行星探测器的自主导航与安全着陆能力。
1 多视匹配小行星全球模型重建近年来,国内外学者对序列影像的小行星形状重建技术开展了卓有成效的研究。例如,Helfenstein利用Galileo探测器所获取的光学影像重建了951 Gaspra小行星的三维模型[3];Kersten、Preusker等对“黎明号”(DAWN) 探测器所获取的近万幅Vesta小行星表面立体影像进行摄影测量处理,获得了该小行星的全球DTM模型,水平定位精度达到8 m,高程定位精度达到6 m[4-5]。崔平远、邵巍等提出了一种绕飞过程中小行星三维模型重构方法,采用PCA-SIFT(principle component analysis-sacle invariant feature transform) 算法有效提高了三维模型重建的速度[6]。从文献来看,传统小行星全球模型是通过多视匹配、区域网平差等摄影测量处理而重建出来的,但国内外对如何从处于无序状态的序列小行星影像构建三维模型的流程、策略的整体研究较少。
在无初始位置和姿态的条件下,需要首先确定每幅小行星表面影像的方位,其一般做法是:1) 用SIFT(scale invariant feature transform) 或者SURF (speeded up robust features) 匹配算法获取重叠影像中的同名像点,并以此作为连接点;2) 结合稀疏光束法平差算法解算出各个影像的相对位置与姿态;3) 再综合利用核线几何约束、多视最小二乘匹配、半全局匹配等策略进行密集匹配,同时引入影像金字塔分层策略提高匹配效率,最终重建出小行星的三维地形模型,其流程如图 1所示。
为了验证上述方法的有效性,本节从DAWN探测器所获取的Vesta小行星立体影像中随机筛选了80张影像(部分影像如图 2所示) ,来重建Vesta小行星的全球地形模型。拍摄时轨道高度约800 km,从系列影像中可以看出,拍摄角度、光照变化都比较大,影像基本处于无序状态。
利用SIFT特征提取与匹配算法以及RANSAC误匹配点剔除方法,共获取了53 255个连接点,同时恢复了各幅影像的相对位置和姿态,相对位置的平差精度达到了0.4个像素。经过密集匹配与平差处理,重建出一组由130多万个三维点构成的Vesta小行星全球模型,如图 3所示。
实验结果表面,该方法适合批量影像的多视匹配与联合平差处理,能够较好地重建小行星全球模型,但是其重建效率尚不能满足探测器在抵近过程中实时观测目标小行星的要求。
2 抵近过程中立体视觉在线测量在小行星探测器抵近着陆的过程中,常采用光学相机和激光雷达的导航方法来引导,如近地小行星交会探测器(NEAR) 结合光学影像与激光测距信息来估计探测器的运动轨迹;“隼鸟号”探测器(Hayabusa) 借助光学CCD导航相机和激光高度计对Itokawa小行星进行了全球测绘,并在光学系统制导下成功着陆小行星。但是,仅用二维影像来判别三维形貌特征是不稳定的,即便结合稀疏的激光点云仍难实现探测器对小行星上探测区形貌的感知能力。因此,Johnson、Mathies 等提出了一种利用机器视觉和图像处理技术进行软着陆导航运动参数估计的方法[7];崔平远、崔跃军等提出了一种基于双目相机的小行星软着陆导航方法[8],为了满足实时处理的要求,仅对立体影像中的稀疏特征点进行了定位与导航分析。实际上,随着计算机技术与立 体视觉技术的快速发展,完全能够实现立体影像实时重 建目标三维信息,在此基础上本文提出了一种立体视觉 在线实时测量方法,即边抵近飞行边测量三维地形。
该方法是采用立体相机分别从两个视点来观察目标小行星,以获取其在不同视角下的影像,经过相机标定、核线纠正、影像匹配以及前方交会等步骤,可计算出影像中目标点的三维信息。其中,影像匹配采用了半全局匹配算法,对立体影像进行了逐像素匹配以获取密集的同名点。
通过序列立体影像对中的前后影像,可将重建出的独立模型连接成一个统一的形貌模型,进一步实现探测器在抵近过程中边获取立体影像边重建探测区三维地形,并按照顺序将各模型规划之统一坐标系中,其难点在于精确确定相邻两组立体影像对之间的几何关系。本文对相邻立体影像对进行稀疏匹配,提取出对应同名点的目标点坐标并以此作为模型连接点,通过模型间的空间变换关系将各个独立模型归化至统一坐标系中,其实现流程如图 4所示。
SURF算法具有稳定性好、速度快且正确率高的特点,因此本文采用该算法获取了前后影像的稀疏连接点。利用这些连接点可以计算出它们之间的空间变换关系,进而将各个独立模型归化至同一坐标系中,而这种关系的求解采用了空间相似变换公式
$\left[ {\begin{array}{*{20}{c}} {{X_{\rm{T}}}}\\ {{Y_{\rm{T}}}}\\ {{Z_{\rm{T}}}} \end{array}} \right] = \lambda \left[ {\begin{array}{*{20}{c}} {{a_1}}&{{a_2}}&{{a_3}}\\ {{b_1}}&{{b_2}}&{{b_3}}\\ {{c_1}}&{{c_2}}&{{c_3}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} X\\ Y\\ Z \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{X_0}}\\ {{Y_0}}\\ {{Z_0}} \end{array}} \right]$ | (1) |
其中:
由式(1) 可知,式中含有7个未知参数,而一对连接点的方程个数为3个,因此解算空间相似变换参数至少需要3个不在一条直线上的连接点。实际中,为了保证精度与可靠性,常需要4个或4个以上的连接点来计算空间相似变换参数。由式(1) 的线性化方程就可计算出前后立体模型间的变换参数,而由这些参数亦可将后一个模型归化至前一个模型中。因此,深空探测器在抵近小行星的过程中可将实时获取的独立模型逐个归化至一个统一的坐标系中,形成整个探测区域完整的三维地形模型。
3 基于等值线分析的撞击坑特征快速识别小行星表面的地形形貌复杂多样,国内有学者参照IAU官方网站给出了小行星形貌特征及其描述信息,并归纳了其表面普遍存在的5种典型形貌特征,指出了撞击坑是小行星表面最为常见且最为明显的形貌特征,其自动识别方法逐渐成为深空探测领域的研究热点。例如,Kim等将坡度指数应用到撞击坑识别中,使得中大型撞击坑的识别率达到75%左右[9-10]; 同时,Bue等也指出仅仅用深度图或坡度指数提取撞击坑还存在应用范围的限制[11]。Salamuniccar等提出了采用融合CCD影像和DEM数据进行撞击坑的自动提取及识别的方法[12-13]。但是,上述撞击坑特征提取方法都需要逐点对DEM数据进行分析,耗时较长,且受岩石、皱褶等地形影响较大[13-14]。而本节所提方法以线特征代替点特征,大大提高了形貌特征提取的速度,避免了形貌中杂点的干扰。
首先,小行星表面撞击坑特征的空间结构具有以下共同特点[15]:
1) 撞击坑特征的外层等高线常表现为近似圆形或椭圆形;
2) 撞击坑特征的外层等高线内侧可能存在多条嵌套的相似等高线,而外侧等高线形状不定、距离较远;
3) 撞击坑特征的中央剖面线可用下凹的抛物线形态拟合。
依据撞击坑特征的空间结构特点,本文利用等值线的形状与分布规律来快速提取小行星表面的撞击坑特征,其具体流程如图 5所示。
其中,小行星表面局部地形受整个星体曲面模型影响较大,需要拟合局部地形形貌的基本面,即能够描述整个局部区域的光滑曲面。将局部地形形貌信息投影至基本面上,重新生成特征更加明显的形貌信息,以减少局部地形坡度变化对撞击坑特征提取正确率的影响。
针对拟合生成的DEM数据,本文应用移动曲面拟合法来逐点内插栅格点坐标,以获取等间距的规范数据集。再用矢量法来内插局部地形的等值线,即按逐条等高线的走向边搜索边内插点。而撞击坑特征的识别过程如下:
步骤1:保留符合撞击坑特征的等值线。先依据撞击坑内等值线的特点保留所有闭合曲线,再应用最小二乘法拟合这些闭合曲线并获取其圆度。当曲线圆度满足阈值要求时,认为该等值线可能为撞击坑上的等值线,同时记录其中心坐标与半径信息。
步骤2:分析等值线之间的关系,聚合嵌套的同坑等值线。结合等值线的半径及其距离关系,将同一撞击坑的等值线聚合在一起,而且记录它们之间的相互关系。
步骤3:判断等值线区域是否为撞击坑。为了避免突出地形或石块的影响,需要对该区域再进行凹凸性评价判断,拟采用抛物线方程来拟合该区域的切面 信息。
步骤4:应用数学形态学方法提取撞击坑边缘。以判别撞击坑的外层等值线为依据,采用数学形态学方法通过膨胀运算检测撞击坑的实际边缘。
步骤5:精确计算撞击坑特征的位置信息。
4 仿真实验情况为了模拟深空探测器抵近测量小行星表面地形、识别其形貌特征,本文通过软硬件相结合的方式设计了一套小行星表面抵近探测仿真原型系统,硬件包括运动控制模块、探测数据获取模块、数据实时处理模块以及小行星局部地形模型等,其结构示意图如图 6所示。
其中,小行星地形模型是参考了Vesta小行星的局部形貌用细沙堆积而成的,模型包括山地、平原以及撞击坑特征等,整个沙盘的大小为4 m×4 m,如图 7(a)所示。探测数据获取模块主要由光学立体观测相机构成,高度范围为2~3 m,基线长度小于10 cm,如图 7(b)所示。结合张正友标定法,获取了立体相机的内参数及其间的相对位姿关系,标定结果如表 1和表 2所示。
利用光学立体观测相机可连续获取小行星地形模型的立体影像,这里列出其中一条航带共18对立体影像,如图 8所示。通过对立体影像进行核线纠正、密集匹配、空间前方交会处理,可由各独立立体影像生成所摄区域的相对三维模型;依据相邻的立体像对,先获取其同名像点及其对应的地形坐标,再利用空间相似变换将各个独立模型归化至统一坐标系中,生成整片区域的地形模型,如图 9所示。
其中,纠正每帧立体影像耗时约为18 ms;SGM算法中深度范围设置为0~64个像素,密集匹配耗时为120~150 ms;而由同名像点交会对应的目标点坐标耗时仅为30~35 ms;相邻前后影像SURF匹配平均耗时为211 ms。实验中共获取立体影像为51对,包含3条航带且每个航带17对立体影像,从立体影像对获取、独立模型生成至模型坐标的归一化共花费时间约61 s。
为了评价地形形貌重建的精度,实验采用激光扫描仪获取了该模拟地形的三维数据,并随机抽取了10个相对离散的检查点。结果显示,对应点最大高程误差为31.91 mm,平均误差为3.38 mm。误差原因可能是立体相机在运动过程中有抖动现象,使得匹配误差有所放大,加之各立体模型在连接过程中又存在误差传递现象,使得地形重建的误差越来越大。
4.2 由等值线分析法识别模型中撞击坑特征以获取的小行星局部地形为数据,经过预处理后该区域的栅格大小为600×200,应用等值线分析法从地形模型中自动识别其表面撞击坑特征,识别结果如图 10所示。
由上图可知,该地形模型共存在6个撞击坑特征,应用等值线分析法从模型中共提取了5个,没有出现提取错误的撞击坑特征,只遗漏了1个撞击坑特征,整个撞击坑特征提取所耗费的时间约为2.16 s。因此,本次撞击坑特征提取实验的正确率为83.3%,耗时远小于探测器着陆时的悬停时间。遗漏原因是该撞击坑特征正好位于山顶位置,在自动提取过程中受山体等值线的影响大而未被识别。
5 总 结论文针对小行星影像的特点,提出了一种基于序列影像的小行星三维重建方法。通过对序列影像进行特征匹配、稀疏光束法平差、核线几何约束、半全局匹配、多视最小二乘匹配等处理,并用Vesta小行星影像构建了其三维形貌信息。为了满足深空探测器对目标小行星表面地形的实时测量与表征需求,本文结合光学立体观测相机,提出了一种适应探测器抵近观测的立体视觉在线测量方法,有效克服了传统摄影测量处理方法的滞后性;并提出了一种基于等值线分析的撞击坑特征提取方法,与其他撞击坑特征提取方法相比,具有一定的稳定性和实用性。
本文所提的立体视觉测量方法还存在一定的精度误差,尤其是在独立模型连接过程中地形测量误差会越来越大,进而影响小行星探测器的导航与定位。后续研究中,拟在实验系统中加入激光测距仪,可直接测量目标点与探测器之间的距离,并将获取的测距信息与影像信息联合平差处理,以求得到更高精度的小行星地形形貌。
[1] | 崔祜涛, 崔平远. 软着陆小行星的自主导航与制导[J].宇航学报, 2002, 23(5): 1-4. |
[2] | 徐青, 耿迅, 蓝朝桢, 等. 火星地形测绘研究综述[J].深空探测学报, 2014, 1(1): 28-35. |
[3] | Helfenstein P, Veverka J, Thomas P C, et al. Galileo photometry of asteroid 951 Gaspra[J]. ICARUS, 1994, 107(1): 37-60.doi: 10.1006/icar.1994.1005 |
[4] | Kersten K D,Matz F,Preusker F,et al. High resolution Vesta HAMO atlas derived from DAWN FC images[C]//43rd Lunar and Planetary Science Conference. The Woodlands,TX,USA:[s.n.],2012. |
[5] | Preusker F,Scholten F,Matz K D,et al. Global shape of (4) Vesta from Dawn FC stereo images[C]//Vesta in the light of Dawn:First exploration of a protoplanet in the asteroid belt.[S.l.]:[s.n.],2014.http://cn.bing.com/academic/profile?id=9346a5ec3b60a3755b1183f4ccadfa79&encoded=0&v=paper_preview&mkt=zh-cn |
[6] | 崔平远, 邵巍, 崔祜涛. 绕飞过程中小天体三维模型重构及探测器运动估计研究[J].宇航学报, 2012, 31(5): 1381-1389. |
[7] | Johnson E A, Mathies H L. Precise image-based motion estimation for autonomous small body exploration[J].Robotics and Automation in Space, 1999: 627-637. |
[8] | 崔跃军. 基于双目相机的小行星软着陆导航方法及其仿真平台研究[D]. 哈尔滨:哈尔滨工业大学,2009. Cui Y J. Algorithms based on binocular vision for softlanding asteroid[D]. Harbin:Harbin Institute of Technolog,2009. |
[9] | Kim J R, Muller J P, Gasselt S V, et al. Automated crater detection:a new tool for Mars cartography and chronology[J]. Photogrammetric Engineering and Remote Sensing, 2005, 71(10): 1205-1217.doi: 10.14358/PERS.71.10.1205 |
[10] | Michael G G. Coordinate registration by automated crater recognition[J].Planetary and Space Science, 2003, 51(9): 563-568. |
[11] | Bue B D, Stepinski T F. Machine detection of martian impact craters from digital topography data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45(1): 265-274.doi: 10.1109/TGRS.2006.885402 |
[12] | Salamuniccar G, Loncaric S, Grumpe A,et al. Hybrid method for detection of lunar craters based on topography reconstruction from optical images[C]//Conference:Image and Signal Processing and Analysis.[S.l.]:ISPA,2011. |
[13] | 罗中飞, 康志忠, 刘心怡. 融合嫦娥一号CCD影像与DEM数据的月球撞击坑自动提取和识别[J].测绘学报, 2014, 43(9): 178-185. |
[14] | 王栋, 邢帅, 徐青, 等. 一种基于三维形貌的深空星体表面撞击坑自动提取方法[J].测绘科学技术学报, 2015, 32(6): 619-625. |
[15] | 王栋, 邢帅, 徐青, 等. 一种类球型小行星表面撞击坑自动提取方法[J].宇航学报, 2016, 37(4): 470-479. |
[16] | 蓝朝桢, 耿迅, 徐青, 等. 基于序列影像的小天体三维形状重建方法研究[J].深空探测学报, 2014, 1(2): 140-145. |