多视角多尺度陨石坑特征检测识别方法

高锡珍，汤 亮，崔平远

(1. 北京控制工程研究所,北京 100094； 2. 空间智能控制技术重点实验室,北京 100094；3. 北京理工大学深空探测技术研究所,北京 100081)

0 引 言

着陆探测及采样返回是未来深空探测的主要发展方向。为了获得更有价值的科学数据，未来希望探测器具备在复杂地形区域中精确避障及位姿确定的能力。然而复杂地形区分布着岩石、沟壑及陨石坑等地貌，给探测器安全探测带来很大威胁。

陨石坑是太阳系天体上常见的自然形貌，随时间推移在不同光照条件和相机位姿下具有一致的几何轮廓(圆或椭圆)和明暗区域特性，因此陨石坑是可用于导航的观测特征。此外，陨石坑边缘地势陡峭，是障碍地形之一，且极具科研价值，需要对其进行规避及跟踪。Bandeira等针对陨石坑特有明暗结构强化边缘，分析边缘形状生成相应的模板，实现陨石坑模板匹配提取。文献[4-6]利用陨石坑明暗特性生成多尺度模板进行陨石坑明暗区域提取与配对，为减少模板数量，根据光亮-阴影区域所占比不变性生成固定模板检测陨石坑。但该算法只能获取陨石坑的中心点，不能确定陨石坑半径。Maass等针对各种不同形状的图像区域，利用最大稳定极值区域法(MSER)提出了尺度、位置以及光照不变的陨石坑检测方法。上述方法提取陨石坑后可利用基于区域的匹配算法实现不同图像间陨石坑匹配。Cheng等提出了基于边缘的陨石坑检测算法，利用Canny算子提取图像中的边缘曲线，筛选并配对符合陨石坑形态特征的曲线，然后拟合配对曲线为椭圆来表征陨石坑。在陨石坑匹配方面，通过构建陨石坑椭圆曲线间的仿射不变信息，提出基于仿射不变量的陨石坑匹配算法。Yu等提出了一种结合明暗区域信息与图像边缘的陨石坑检测与匹配算法，基于椭圆曲线中仿射不变量采用WTA(Winner takes all)策略进行匹配，最后采用三维陨石坑模型消除误匹配。针对陨石坑阴影不明显的行星着陆导航图像，陈建清等提出了基于图像灰度值特征的陨石坑自主检测方法，解决了一般陨石坑检测方法对图像太阳高度角的依赖问题。针对3D特征检测与匹配问题，结合数字高程图，骆磊等通过计算特征坡度和曲率完成了月球陨石坑的检测。Salamuniccar等利用数字地形图取代光学图像，通过进行梯度及形态学分析，利用模糊边缘检测器实现了行星表面陨石坑特征的提取。

除上述直接利用陨石坑形态特性、模板和相对几何构型进行陨石坑识别外，文献[13-14]研究了基于特征描述符进行陨石坑检测与匹配的方法，首先利用高斯金字塔和ELSD算法提取陨石坑圆弧，在支撑区域内构造圆弧带描述符，结合最近邻距离比(NNDR)和欧式距离约束实现曲线匹配，通过仿真验证了该算法具有缩放、旋转、光照不变性。但是该方法存在陨石坑边缘误检测及重复检测的问题。针对陨石坑检测存在不充分和不准确的问题，随着机器学习的发展，Kang等和Wang团队分别提出了基于监督学习、半监督学习以及主动学习的天体表面陨石坑检测算法。在保证样本充足情况下，机器学习算法对不规则退化陨石坑都能达到良好识别效果，但目前来说该算法仅适用于在轨观测和离线分析，并不适用于着陆等实时性要求较高的过程。

地外探测着陆与大范围巡视过程中，探测器位置及视角动态变化导致目标图像特征尺度、形状等外观模型变化，且序列图像间存在信息冗余，分别检测降低计算效率，容易造成目标漏检测。因此，考虑到障碍规避和精确导航的需要，有必要快速有效地检测出陨石坑。本文提出了多视角多尺度陨石坑特征检测识别方法，基于明暗区域信息和形态学处理实现陨石坑图像特征边缘粗提取，利用快速傅里叶变换计算模板与图像特征的匹配度，从而实现图像特征所在区域粗定位；在此基础上，通过引入惯性测量信息计算图像间重叠面积，更新模板形状，并预测搜索区域，解决序列图像特征视角尺度变换及信息冗余带来的特征无法检测和识别效率低的困难，最终在单一陨石坑所在区域内实现其精细检测。

1 陨石坑检测识别流程

本文提出的陨石坑提取算法将结合惯性信息重点解决序列图像间视角尺度变换及信息冗余带来的特征无法检测和识别效率低的困难，算法的总体流程如图1所示。首先引入惯性测量信息计算序列图像间重叠面积，在此基础上，基于形态学处理实现陨石坑图像特征边缘粗检测，利用快速傅里叶变换计算模板与图像特征的匹配度，从而实现图像特征所在区域粗定位，在此基础上，考虑到序列图像间视角尺度变换，结合惯性信息更新模板形状，并预测搜索区域，最终在单一陨石坑所在区域内实现其精细检测。

图1 陨石坑特征检测识别方案流程Fig.1 Scheme flow of crater feature detection and recognition

2 多视角多尺度陨石坑特征粗识别

2.1 序列图像重叠区域计算

陨石坑作为障碍地形具有重大科研价值，需要对其跟踪并有效规避，陨石坑需要在序列图像的重叠区域内进行选择，避免对整幅图像进行检测，减少了搜索区域。因此，寻找相交的多边形并计算重叠区域是实现候选区域确定的关键。

(1)

(2)

(3)

(4)

由式(3)计算得到探测器运动后的相机视场与星表平面交点如式(5)所示。

(5)

为了计算重叠区域，需要考虑四边形顶点以及四边形边界之间交点的位置，算法步骤如下：

(1)对于第一个四边形的每个顶点，判断它是否包含在第二个四边形里面。如果是，存储该点坐标；

(2)对于第二个四边形的每个顶点，判断它是否包含在第一个四边形里面。如果是，储存该点坐标；

(3)对于其中一个四边形的每条边，判断其与另一个四边形的每条边是否相交。如果是，计算并储存交点坐标；

(4)计算每个点对应的图像点位置，连接每个像点形成凸多边形；

(5)基于相机透视成像原理确定两幅图像中重合区域，在此范围内进行陨石坑候选区域粗检测。

基于上述步骤能够有效地实现序列图像间重合区域确定，从而为陨石坑候选区域粗检测提供前提条件。

2.2 候选区域粗检测

陨石坑形貌近似为凸出地表的环形，在良好的光照条件下，陨石坑边缘存在明显的明暗变化。如果这些明暗区域被认为是可能的陨石坑的标志，那么所有与边缘明暗变化等同的区域都应与图像背景分离并增强。因此，考虑将明暗区域信息引入到边缘检测中，以得到陨石坑边缘所在候选区域。

对于某一像素点(,)，基于区域的边缘检测方法选择以其为中心，大小为×的图像窗口，首先计算窗口内灰度均值()，然后寻找窗口内像素点最大灰度值max()和最小灰度值min()，并计算窗口内灰度均值()分别与max()和min()的差值，最后将该值分配给原图像像素点。具体计算过程如式(6)所示。

(,)=max[()-min(), max()-()]

(6)

为了最小化噪声影响，取=3。将带有×元素的窗口以原始图像的每个像素为中心重复式(6)所示的操作过程，从而创建一个新图像。

根据图像的灰度值范围，计算阈值，如式(7)所示。

=[max()-min()]+min()

(7)

式中：为常数。

利用阈值，对图像进行如下二值化处理，得到陨石坑边缘候选区域的二值图像。

(8)

2.3 候选区域定位

通过区域检测得到了一系列陨石坑边缘候选区域，但光照和星表地形复杂多样导致候选区域不可能全为陨石坑边缘所在区域。在边缘提取过程中，直接对候选区域进行处理存在较大偏差，如果候选区域为虚假区域，会浪费大量计算资源，甚至使结果不可理解。因此，本节利用模板图像和候选区域进行相似性比较，筛选候选区域，从而确定陨石坑边缘所在真实区域。

考虑到陨石坑尺度稀疏特性，会产生大量冗余模板，进而导致大量无意义计算。因此，首先需要分析图像中陨石坑边缘候选区域大小。本节利用连通域分析陨石坑边缘候选区域大小，并在此基础上最终确定陨石坑边缘区域。为了保证候选区域检测充分性，采用四连通域分析候选区域大小。四连通域如图2所示。

图2 四连通域Fig.2 Four connected domain

通过对候选区域的四连通域进行分析，得到候选区域上相距最远的两点和，分别计算其水平方向和垂直方向上的距离，取最大值定义为模板图像中圆环外径。模板图像圆环外径确定如式(9)所示。

=max(|max-min|,|max-min|)

(9)

式中：max,max,min和min分别表示两点和在图像坐标系下的横纵坐标值。

创建二值化模型作为搜索模板，由于陨石坑自然形状不规则，搜索模板设定为圆环，其内部半径是外部半径的一半，周围是一个黑色正方形的背景。同时，陨石坑特征尺度大小不一，需要生成半径大小不同的模板图像。不同大小模板图像如图3所示。

图3 模板形状Fig.3 Templates

在此基础上，利用快速傅里叶变换(Fast Fourier transform, FFT)进行相关性分析，可确定陨石坑边缘区域，避免产生大量冗余模板，提高了计算效率。

2.4 模板更新

在实际中，大多数的陨石坑在星表以圆形形态分布。然而，在序列图像中，由于星载相机并不一直垂直指向行星表面，陨石坑图像并不呈现为圆形。探测器运动过程中位置和指向的持续变化进一步导致了图像畸变。图4显示了视角尺度变换对模板的影响，原始的特征模板为圆环，修改后的特征模板标记为椭圆环。半径为的圆环，经过视角尺度变换(+1,+1)形成椭圆环，变换后的椭圆可以由参数长轴、短轴和长轴倾角表示。因此在序列图像中应用相同的基于圆环的特征检测模型可能会造成不可估量的精度损失，需要特征模板更新。

图4 视角尺度变换后的模板变化Fig.4 Template change after view and scale transformations

序列图像间像点的几何关系可以描述为：

(10)

式中：表示相机焦距；和表示第幅图像中像点坐标；+1和+1表示第+1幅图像中像点坐标；,和表示探测器在位置处时星表平面点在相机坐标系中的三维坐标；+1表示探测器在位置+1处时星表平面点在相机坐标系中的轴坐标。

初始模板图像表示如式(11)所示。

(,,1)(,,1)=0

(11)

式中：表示初始模板图像的二次型矩阵。

结合式(10)和式(11)计算得到视角尺度变换后模板图像表达式为：

(,,1)+1(,,1)=0

(12)

式中：+1表示变换后模板图像的二次型矩阵。

将式(12)改写成一般解析式为：

+++++=0

(13)

式中：参数,,,,和分别表示模板图像一般方程的系数。

为了更加直观地描述变换后模板图像，此时图像为椭圆,且满足≠,将式(13)改写成中心点(,)、长短轴,和长轴倾角参数表示：

(14)

=

(15)

(16)

(17)

同时探测器运动导致图像中目标运动，而运动后的目标往往只在图像中一个非常小的区域出现，所以在搜索匹配时，搜索位置可通过估计值(,)确定，比较当前及其邻域相关系数大小，找出合适区域。在行星着陆任务中，位姿估计信息可以通过位姿确定系统(如惯性测量信息)提供，据此可求得下降图像中当前帧相对前一帧图像的变换参数，从而更新特征模板，实现陨石坑候选区域确定。需要指出的是，本节利用惯性敏感器件提供初始状态信息，后续状态信息在陨石坑特征匹配基础上可通过视觉导航系统修正和提供。

3 陨石坑边缘精细提取

利用边缘检测算法直接对整幅图像处理存在一定局限性，因为设定较低的阈值可以保证边缘的完整性，但是也必然引入其它陨石坑边缘和大量虚假边缘，增加曲线拟合的困难和计算量，同时也给陨石坑的匹配带来困难，而较高的阈值会漏掉那些灰度变化不剧烈的边缘，从而造成边缘检测不完全、不充分。因此，本节利用上节提取的陨石坑候选区域，把每一个候选区域看作一个稀疏型陨石坑分布，然后精细提取陨石坑的边缘。

为了减少噪声响应，提高边缘检测正确率以及限制单个边缘点对于亮度变化的定位，对候选区域图像采用Canny边缘检测方法。首先利用高斯滤波器进行图像去噪声，利用sobel算子包含横向及纵向两组3×3的矩阵，与候选区域图像作平面卷积运算，得到纵向及横向的灰度差分值，近似为像素梯度，最后非极大值抑制多响应边缘，使边缘的宽度尽可能为1个像素点。在完成了陨石坑真实边缘的提取后，利用提取出的边缘采取投票策略拟合为椭圆，确定椭圆参数。

边缘精细提取具体步骤如下：

(1)利用Canny边缘检测方法提取某一候选区域图像中的陨石坑边缘，以候选区域中心为陨石坑中心；

(2)随机抽取5个像素并拟合出对应的曲线,计算拟合出的椭圆中心，判断椭圆中心点与候选区域中心的距离，如果距离在一定阈值范围内，判定曲线为备选陨石坑边缘曲线；

(3)设定阈值，计算该条曲线上所有剩余的点与该椭圆的距离，如果符合，判定该点近似的在该条曲线上，并记录；

(4)判断点的个数，若多于设定的最少数目，认定该拟合椭圆是正确的,并用全部在椭圆曲线上的点重新计算椭圆方程，计算新获得的椭圆与全部点距离的平方和，以此作为拟合误差,最后在正确的拟合中挑选拟合误差最小的一组参数作为曲线参数的估计结果;

(5)对每一个候选区域，重复完成上面的步骤(1)～(4)，拟合得到所有候选区域的陨石坑边缘参数，从而得到在原图像中的陨石坑边缘参数。

4 仿真校验

本节采用美国“火星全球勘测者”(Mars Global Surveyor)在绕火星飞行时利用火星轨道相机(MOC)所拍摄的火星表面图像进行实验分析。该相机获取的图像分辨率较高，覆盖范围广，包含的形貌特征多样，能够反应星表陨石坑图像的特点。

基于2.2节所述方法对原始图像(图5(a))进行处理，常数值取0.25，陨石坑边缘候选区域检测结果如图5(c)所示。

图5 陨石坑边缘候选区域Fig.5 Crater edge candidate area

对于图5(a)所示的火星表面局部图像，利用边缘检测算法提取其中的陨石坑边缘。这里分别采用基于一阶微分算子、基于二阶微分算子、基于Canny算子以及Ostu-Canny算子的边缘检测方法处理图像并对处理结果进行对比分析。图像检测结果如图6所示。

图6 不同算子边缘检测结果Fig.6 Feature detection results of different operators

可以看出上述方法对于边界较为明显的陨石坑都获得了良好的检测结果，一阶微分算子对灰度渐变区域处理效果较好，二阶微分算子灰度突变定位准确，Canny边缘检测不易受噪声干扰，边缘定位准确，Ostu-Canny边缘检测阈值自适应，但都存在边缘检测不完全、伪边缘过多的缺点，无法精确确定陨石坑边缘轮廓信息。而利用区域信息保留了所有陨石坑所在区域，同时有助于后续进行模板匹配。

图7给出了MOC在不同视角下拍摄的火星表面图像。基于所研究的陨石坑识别方法，首先利用惯性测量信息计算序列图像间重合区域，如图8所示。

图7 火星表面序列图像Fig.7 Mars surface sequence images

图8 星表图像重合区域Fig.8 Overlapping area between neighbored scenes

利用基于区域的边缘检测方法检测重合图像中的陨石坑特征。为了消除细小特征并平滑特征边界，在陨石坑特征区域检测的基础上对特征进行形态学闭操作，陨石坑边缘候选区域检测结果如图9所示，其中图9(a)表示利用惯性测量信息计算得到的图像间重合区域，即待处理图像，图9(b)和图9(c)分别表示按照式(6)和式(8)对图9(a)处理的结果，图9(d)表示对图9(c)进行形态学处理的结果。

从仿真结果可以看出，在不同光照条件下,不同尺度的陨石坑边缘区域都得到了很好的检测，重合图像中10个较为清晰的陨石坑，其中有9个陨石坑区域获得了有效检测，包括轮廓不完整的陨石坑，如图9(d)所示，检测率为90.00%，验证了算法在运动过程中光照方向及探测器位置变化等复杂环境下的有效性。

图9 陨石坑边缘候选区域检测Fig.9 Crater edge candidate region detection

在序列图像中，由于星载相机可能不总是垂直指向行星表面，陨石坑图像并不呈现为圆形。探测器运动过程中指向变化进一步导致了陨石坑形状参数变化，如图10所示。

图10 视角尺度变化后图像Fig.10 Image with view and scale distortions

视角及尺度变化后图像中陨石坑边缘轮廓呈现出明显的椭圆形，利用惯性测量信息更新初始模板，通过FFT算法计算不同大小模板与候选区域之间的相关系数，搜索并确定相关系数最大位置。图11和图12分别给出了视角尺度变化后的陨石坑边缘确定区域相关系数和陨石坑边缘区域检测结果。仿真结果表明，利用模板匹配较为准确地确定了陨石坑边缘区域，说明检测方法不会因为运动过程中探测器姿态的改变而失效，确保了算法的有效性。

图11 FFT计算结果Fig.11 FFT calculation results

图12 边缘区域确定Fig.12 Edge area determination

确定所有陨石坑所在区域，对区域内单一陨石坑实现精细检测，精细检测过程参考文献[8]中Canny边缘检测算子提取图像中的边缘曲线。不同于Cheng等对全局图像进行边缘算子提取，此处仅针对模板确定的单一陨石坑所在区域进行Canny边缘检测，检测对比结果如图13所示，其中图13(a)表示单一陨石坑区域图像,图13(b)表示对全局图像进行边缘提取的结果，图13(c)表示直接对图13(a)进行边缘提取的结果。与图13(c)相比，图13(b)中伪边缘较多，光照复杂情况下难以筛选配对,影响后续拟合效果。因此,本文对单一陨石坑所在区域进行Canny边缘提取并进行拟合,检测结果如图14所示。

图13 两种陨石坑边缘检测算法结果Fig.13 Result of two crater detection algorithms

图14 陨石坑检测结果Fig.14 Crater detection results

图14表明轮廓较清晰完整的陨石坑，甚至其中稠密分布的陨石坑，都成功被检测出来，但由于设定模板大小以及图像噪声影响，原始图像中尺寸较小的陨石坑以及形状不规则的陨石坑没有被检测出来。在避障及导航应用中，通常选定那些形状较为规则，轮廓更为清晰明显的陨石坑进行匹配和跟踪，以获得更好的匹配和跟踪效果，从而保证探测器探测的安全性和可靠性。

为进一步验证本文所提算法的有效性，对不同类型陨石坑特征进行仿真，以行星图像数据为基础，结合动力学来模拟着陆序列图像，模拟图像及检测结果如图15所示。算法在不同类型陨石坑检测中同样取得了理想的检测效果,特征误检率低,这是由于本算法使用了区域与边缘联合信息进行检测识别，提高了陨石坑选取标准,因此剔除了更多的误检测,但由于检测模板大小设定原因,一定程度上降低了不明显陨石坑的检测率。

图15 陨石坑检测识别算法示例Fig.15 An example of the crater detection and recognition algorithm

5 结 论

针对星表形貌复杂类型多、特征尺度变化大问题，提出了多视角多尺度陨石坑特征检测识别方法。基于明暗区域信息和形态学处理实现了陨石坑图像特征边缘粗提取，避免了其它不规则纹理特征干扰，利用快速傅里叶变换计算模板与图像特征的匹配度，实现了图像特征所在区域粗定位，缩小了陨石坑检测时的搜索区域；在此基础上，通过引入惯性测量信息计算图像间重叠面积，更新模板形状，并预测搜索区域，解决了序列图像特征间视角尺度变换及信息冗余带来的特征无法检测和识别效率低的困难，最终在单一陨石坑所在区域内实现了其边缘提取。