一种基于等高线的小型撞击坑识别方法*

刘宇轩，李春来，刘建军

(1. 中国科学院国家天文台，北京 100101；2. 中国科学院大学，北京 100049，3. 中国科学院月球与深空探测重点实验室，北京 100101)

行星表面通常分布着大量的撞击坑构造，它们是行星表面最显著的地貌特征，也是研究行星内部物质的窗口[1]。由于月球上没有大气，撞击坑保存比较完好，因此研究月球形貌、相对地质年龄、岩石构造以及月球的起源有重要意义。撞击坑识别的研究，从早年使用望远镜观察描绘撞击坑的形貌，到利用遥感图像数据目视识别撞击坑，再到近年来通过计算机自动识别撞击坑。随着月球影像与地形数据的大量获取，撞击坑的自动识别提取成为一个新的研究热点。近年来，国内外众多学者提出大量的研究方法，并建立多种尺度的撞击坑数据库[2-6]。通过对现有方法的分析总结，撞击坑识别的方法大致可以归纳为4类：人工目视识别、基于形态拟合、基于机器学习以及基于地学信息分析的方法。人工目视识别法：利用影像和地形数据，建立识别撞击坑的解译标识，在计算机环境下，采用人工交互的方式对撞击坑逐个识别，这种方法是目前最准确的撞击坑识别方法，该方法的优点是准确性高，但费时费力，不适合于全月球和小尺度撞击坑识别。基于形态拟合法：由于撞击坑在图像上有类圆形的几何特点，在图像分割的基础上以霍夫变换和二次曲线为代表的形状拟合算法用于模拟提取完整的撞击坑边界，其中文[7]利用组合霍夫变换算法，文[8]采用模糊霍夫变换算法对撞击坑进行识别，文[9]提出一种基于弦中点霍夫变换算法，文[10]采用最小二乘法拟合边缘椭圆的方法实现对撞击坑的提取。基于机器学习的方法：文[11]采用BP神经网络加训练样本的算法，文[12]利用决策树算法，文[13]利用面向对象方法实现对撞击坑的识别。基于地学信息分析方法：文[14]利用坡度作为识别撞击坑的要素，文[15-17]利用500 m分辨率的嫦娥一号数字高程模型数据和形态拟合的方法，进行撞击坑的识别提取和研究，文[18]利用坡度数据、纹理信息和剖面曲率作为识别依据，采用水淹法对火星撞击坑进行分析。

本文使用目前世界上覆盖最全、精度最高的嫦娥二号数字地形数据，提出基于等高线地形分析的方法，针对小型撞击坑进行识别提取研究，填补撞击坑数据库小型撞击坑数量不足和准确率不高的缺陷。

1 数 据

2008年2月，国务院批准探月工程二期立项，二期工程主要目标是实现月面软着陆，开展月面就位探测与自动巡视勘察。嫦娥二号作为二期工程的技术先导星，首要的任务是对全月球和嫦娥三号的月面着陆区进行详查，精细测绘全月球的地形地貌。嫦娥二号卫星于北京时间2010年10月1日18时59分57秒在西昌卫星发射中心成功发射。10月24日16时49分，嫦娥二号卫星上的两线阵CCD立体相机首次开机工作，并成功获取了前视与后视两个视角的影像数据。截至2011年5月20日，嫦娥二号卫星CCD立体相机完成全部任务，共获取607轨图像数据。全月球影像图制作过程中选取的100 km轨道高度成像的384轨数据(其中南北纬70°范围344轨，极区40轨)，用于全月球数字高程模型数据产品生产。数据处理流程包括：数据准备、测区数据组织、平差点准备、全球平差、测区初级产品生产、标准产品生产、全球数据镶嵌等。结合嫦娥二号卫星轨道与姿态测量数据、5个月面绝对控制点进行全球平差，经过4年的数据处理，成功制作了7 m、20 m、50 m分辨率的数字高程模型数据。嫦娥二号数字高程模型数据在空间分辨率、数据一致性和完整性镶嵌精度等方面优于国际同类全月球数字产品。本文选择7 m分辨率的数字高程模型数据进行小型撞击坑的识别提取。

2 方 法

2.1 撞击坑模型分析

根据基本形态，撞击坑可分为简单撞击坑(碗型小型撞击坑)、复杂撞击坑(距中央峰的较大型撞击坑)、多环撞击坑(撞击盆地)3类[19]。如图1，一个理想状态下的复杂撞击坑模型主要由坑唇边缘(rh)、坑壁(w)、坑底(f)、撞击培岩(m)和中央峰或中央环(p)5种形态要素组成。撞击溅射堆积物由坑唇边缘(rh)、丘状沉积物(rr)和辐射线沉积物(rc)构成。撞击溶岩(m)覆盖在其他形态要素之上。简单碗型撞击坑(c)一般分布在丘状沉积物(rr)之上。撞击坑识别是一种特征识别，对影像或地形数据进行分析处理，如图2，识别撞击坑的坑唇边缘(rh)，并提取撞击坑的直径(D)、中心位置和高程等信息。

图1 理想的的复杂撞击坑模型(Don wilhelms, 1987)
Fig.1 Ideal complex impact crater model

等高线是地形图上高程相等的相邻各点连成的闭合曲线，是最为常用的表达地形形态特征和起伏情况的形式，如图3。等高线从数字高程模型中提取出来，科学地反映月表的高程 、坡度、坡向、高地走势等基本月表地形地貌。嫦娥二号数字高程模型数据的空间分辨率很高，月表高程在二维地理空间上得到了连续表达，利用该数字高程模型数据较为方便、精细地提取月球撞击坑地形等高线，表达月球撞击坑形态，并对其进行分析和处理。

从等高线的形态分析，如图4撞击坑具有3个典型的空间特点：(1)撞击坑坑唇边缘等高线通常为近似圆形或椭圆形；(2)撞击坑坑唇边缘等高线内侧通常存在多条嵌套的相似等高线，多为同心圆或椭圆，且间距较小，外侧等高线形状不定，且间距较远，与坑唇边缘有明显区别；(3)撞击坑坑底高程小，坑边缘高程大，其剖面线可以用下凹的抛物线形态拟合。因此，通过撞击坑地形等高线的特征以及与其他等高线之间的相互关系，可以分析并识别月表的撞击坑。因此，只要选取合适的等高距，撞击坑坑唇边缘的闭合近圆形等高线就可以逼近真实边缘，可以准确识别并提取撞击坑的地理位置、大小和深度等信息。

图2 简单撞击坑和复杂撞击坑剖面图
Fig.2 Impact crater structure

图3 等高线示意图
Fig.3 Contour diagram

图4 撞击坑等高线示意图
Fig.4 Crater contour diagram

2.2 指标体系和参数的确定

撞击坑等高线往往表现为相互嵌套的类圆环状。因此，初步判别撞击坑特征的思路是根据等高线的圆度、矩形度、体态比等指标参数保留满足条件的等高线，再分析等高线之间的关系，聚合嵌套的等高线并标示外层等高线；然后通过剖面线的凹凸剔除穹隆等正向地形区域，保留负向地形的撞击坑区域；最后通过霍夫变换拟合撞击坑的边缘，得到撞击坑的中心位置和大小。因此，指标体系和参数包括圆度、矩形度、体态比、剖面线凹凸度。

(1)圆形度是图像处理中重要的概念，用于特征提取与描述，其值反映被测量边界的复杂程度，公式为

e=4πA0/P2,

(1)

其中，A0为面积；P为周长。这个特征对圆形面域取最大值1，坑唇边缘越接近1的等高线越符合条件。

(2)矩形度是常用的目标物体矩形拟合因子，其公式为

e1=A0/AR，

(2)

其中，A0为面积；AR为最小外接矩形的面积；e1表示某物体对其最小外接矩形的充满程度。对于矩形物体，e1取最大值1；对于圆形物体，e1取0.79；对于其他细长物体，e1取值变小。矩形拟合因子的值限定在0～1之间。

(3)体态比定义为

e2=W/L，

(3)

其中，W为最小外接矩形的宽；L为最小外接矩形的长。圆形面域的体态比取值为1，这个参数可以把细长不规则的等高线剔除。

(4)剖面线凹凸性：因为月表还存在其他隆起的地形，如穹隆也会产生近圆形的闭合等高线。潜在撞击坑还存在，对此情况的处理方法是判断潜在撞击坑所在的数字高程模型剖面线是否为下凹形态，如果剖面线为下凹形态，则是撞击坑；若剖面线为上凸形态，则不是撞击坑。根据内等高线与等高线的数值大小关系确定是否为真实撞击坑。

2.3 撞击坑边界的拟合

对满足条件的等高线，利用霍夫变化的方法拟合其边界，最终得到中心点和半径大小。

霍夫变换算法是图像处理中识别几何图形以及其特征的一种常用方法。霍夫变换的思想是：对二值图像I上任意一个边缘点I(x,y)，如果它在以(a,b)为圆心，r为半径的圆上，则满足：

(x-a)2+(y-b)2=r2.

(4)

从(4)式可以看出，每个边缘点I(x,y)映射成空间 (a,b,r)上一个二次曲面。将(4)式写成参数 形式：

a=x-rcosθ,b=y-rsinθ，

(5)

其中，θ∈[0, 2π)。设待检测圆的半径r∈[R1,R2]，则按照霍夫变换，对任意一个边缘点I(x,y)，θ，r分别以Δθ，Δr的步长遍历[0, 2π)和[R1,R2]，得到参数空间(a,b,r)的子空间，并对相应的累加器单元AC(a,b,r)加1。对图像上所有的点遍历之后，得到累加器数组AC，其中任意一点AC(a,b,r)表示在圆(a,b,r)上的边缘点的数目。相同数目的边缘点在不同半径的圆上被认为是圆的可能性不同，半径越小，可能性越大。因此定义圆度的概念：

p=n/r，

(6)

其中，n为圆的边缘点的数目 ；r为圆的半径；p为圆度。将累加器数组AC中的元素按(6)式计算圆度，得到三维圆度数组C。圆形目标对应于圆度数组的局部峰值，即对应的C(a,b,r)为局部最大值。

3 结 果

3.1 识别流程

在嫦娥二号7 m全月数字高程模型数据中随机选择G227区、G228区进行实验，基于ENVI和IDL，按照章节2中的原理和方法，编写程序进行识别提取，如图5。首先，设置合理的等高距在嫦娥二号数字高程模型数据中提取等高线；然后，通过圆形度、矩形度、体态比、剖面凹凸性筛选出初选撞击坑区域；最后，通过霍夫变换提取撞击坑边界信息。将提取的信息导入ArcGIS，确定撞击坑的地理位置和半径大小，并通过高程差计算得到撞击坑的深度信息，生成具有位置、大小、深度的小型撞击坑数据库。

图5 IDL程序计算演示
Fig.5 IDL program diagram

3.2 识别结果

基于数字高程模型数据等高线自动识别结果如图6，具体统计结果见表1。

本文的方法在G227区、G228区分别识别撞击坑728个和989个。为了评价本文方法的精度，使用Shufelt[20]提出的算法评价因子：检测百分比为(Detection percentage)D=100TP/(TP+FN)，分支系数为B=FP/TP，质量百分比为Q=100TP/(TP+FP+FN)。TP表示正确识别的撞击坑数量；FN表示遗漏识别的撞击坑数量；FP表示错误识别的撞击坑数量。目前最好的评价方法是采用人工识别的结果作为对比参照。

图6 撞击坑识别结果示意图
Fig.6 The crater detection result diagram

表1撞击坑识别结果统计表
Table1Craterdetectionresultstatistics

TP/个FN/个FP/个D/%B/%Q/%G227697533192.94.489.2G228918817191.97.781

从表1可知，本文方法在两个实验区的评价因子D分别为92.9%、91.9%，均值为92.4%，B分别为4.4%、7.7%，均值为6.1%，Q分别为89.2%、81%，均值为85.1%。从评价因子可以看出，本文的方法对于小型撞击坑有非常好的检测结果，正确识别率高，且错分率低，存在小部分遗漏识别的情况。

3.3 讨 论

本文提出一种基于高精度数字高程模型数据的派生数据等高线，通过计算机自动识别提取小型撞击坑的方法。相对于传统的逐点分析坡度指数，提取图像边缘转换成对地形线的分析，可以有效避免地形皱褶和碎石群的干扰，识别流程更为简单高效，识别速度较快。设置不同的等高距，分批次分区间识别不同尺度的撞击坑，其正确识别率高，遗漏识别率低，说明该方法是可行有效的。不过需要注意的是，对于部分退化严重的撞击坑，坑缘破坏严重，坑壁平缓，等高线不能明显表达和准确勾画出撞击坑的形貌轮廓，无法被准确识别，导致漏分情况的发生。但通过对比同位置的影像数据，通过人工识别分析是可以分辨的。因此，今后的研究方向和重点在于如何识别退化严重的撞击坑，这也是所有撞击坑识别方法中的薄弱环节和有待解决的问题。在实际应用中，可用该方法的识别结果补充现有月球撞击坑数据库中的小尺度撞击坑，丰富现有的月球撞击坑数据库。