月表卫星影像的陨石坑几何特征研究

黄 磊，王梦雪，李 淼

（成都理工大学 地球探测与信息技术教育部重点实验室，四川 成都610059）

0 引 言

陨石坑是月表最明显的形貌特征，它对人们研究月球演化、太阳系陨石撞击历史和空间环境等都有重要意义，在陨石坑的研究中，陨石坑影像几何特征的信息提取是一种重要手段［1］。随着科技的进步，很多国家发射了大量的月球探测器，如何利用月球探测器的高分辨率图像研究月表陨石坑的几何特征和月表年龄关系，已成为一个研究的难点与热点问题。许多学者从不同角度进行了分析研究，如使用陨石坑大小频度法估测月表年龄，建立表暴露年龄与单位面积中陨石坑数量关系等［2，3］。但是，这些研究只是定性描述陨石坑几何特征，没有对陨石坑几何特征进行系统的分类以及定量描述陨石的某一几何特征。

成都理工大学大学空间信息研究室在此研究基础上，利用 “嫦娥一号”卫星CCD 立体相机获取的月表遥感影像，研究陨石坑的几何特征，提出了一种 “六位数字十级分级法”，用0－9个数字量化描述陨石坑的几何特征，为描述陨石坑的影像和几何特征建立模型和标准。六位数字分别是直径、边缘清晰度、边缘圆度、坑唇宽度、坑底结构复杂度和高程。本文是研究边缘圆度模型的建立以及边缘圆度的评价算法，依据圆度误差的定义和评定方法，分别用最小区域圆法、最小二乘圆法、最小外接圆法和最大内接圆法对 “嫦娥一号”CCD 卫星影像的陨石坑边缘进行圆度误差判定，然后利用改进的BP神经网络，组合4种评价算法得到的圆度误差的结果作为圆度等级的评价指数。最后，把圆度评价模型加载到自主开发的数字月球平台上做结果展示，圆度等级设为0－9级，对陨石坑的边缘形变做10级量化。

1 基本思路

1.1 算法综述

圆度误差是指被测实际轮廓相对理想圆的变动量，是一项重要的形位精度指标，是评定许多其它形位误差的基础，关于圆度误差的评定方法研究很多，本文所用的圆度考量因子是圆度误差。根据国际标准化组织新颁发的标准规定了有5种圆度误差的评定方法，分别是最小二乘圆法（the least square circle method，LSCM）、最小区域圆法（the minimum zone circle method，MZCM）、最小外接圆法（the minimum circumscribed circle method，MCCM）、最大内接圆法（the maximum inscribed circle method，MICM）和切比雪夫拟合圆法［4］。本文提出一种基于圆度误差的陨石坑边缘圆度评定算法，该算法对图像增强和降噪平滑后的图像做边缘提取，然后用圆度误差计算方法进行判定，得出边缘圆度量化评价结果。对圆度误差计算方法的选择，本文选用的评定方法分别是最小二乘法、最小区域法、最小外接圆法和最大内接圆法。最后，用标准数据对神经网络进行训练，将4种评价结果作为网络的输出，得到量化的规则度评价等级，算法流程如图1所示。

图1 边缘圆度评价模型算法流程

1.2 边缘圆度的定义

边缘圆度的判定是对月球陨石坑的CCD 卫星图像的边缘进行圆度误差的判定，圆度误差作为陨石坑形变的定性判定，误差越大，边缘圆度就越差，误差越小，边缘圆度就越好。以 “六位数字十级分级法”的角度出发，从定量来看，分为0 到9 级，圆度误差逐渐减小，边缘圆度逐渐提高。

1.3 图像处理方法

1.3.1 图像预处理“嫦娥一号”CCD 图像的预处理对边缘圆度的判定有影响，本文图像的预处理主要采用图像增强和双边滤波。图像增强是根据局部灰度值的统计数据对选取原始图像进行增强处理，增大图像的对比度，有选择地图像中陨石坑的几何特征，同时抑制图像中的其它信息。

为消除图像噪声对月表陨石坑的影响，采用非线性、非迭代的双边滤波方法处理图像。双边滤波在处理相邻像素时，既要考虑到距离上的邻近关系，还要考虑灰度上的相似度，通过对这两者的非线性组合，可以自适应地进行滤波。图像在滤除噪声的同时，能很好地保留图像的边缘信息。预处理结果如图2所示。

图2 原始图像和处理后效果图像

1.3.2 边缘提取

边缘提取的精度对边缘圆度的判定有影响，本文分析多种边缘检测的方法后，采用Canny算法提取月球陨石坑的边缘。由于Canny算子的计算量大，实时性较差，检测边缘的精度没有达到单像素级，因而采用改进型Canny边缘检测算法提取边缘及细化［5］。值得注意的是，因为太阳高度角、CCD 镜头拍摄角度等原因，陨石坑内经常会出现高亮区和阴影区，二者同时出现且它们之间的连线与太阳光入射矢量基本保持平行，如图2（a）所示。特别是在小型陨石坑内，由于半径小，当太阳光高度角接近90°时，会出现大面积的高亮区和阴影区，几乎占到坑底的一般面积，并覆盖大部分坑唇。提取的陨石坑边缘中既包括了真实的边缘信息，又包括了陨石坑内部高亮区和阴影区的伪边缘信息以及其它地形特征的边缘信息，如图3（a）所示。通过检测陨石坑图像的梯度方向，发现在陨石坑内部的伪边缘与陨石坑的真实边缘的梯度方向正好相反，真实陨石坑的边缘梯度方向与光源方向的夹角小于90°，但是伪边缘梯度方向与光源方向的夹角大于90°，即真实边缘满足式（1）

图3 伪边缘去除效果

2 边缘圆度评价算法

边缘圆度的评价实质上是根据圆度误差的定义构造函数模型，然后进行函数的优化求解［6，7］。最小区域圆法（MZCM）是以最小区域圆（MZC）作为评定的基准圆的方法，按此方法求得的圆度误差值ΔZZ

式中： （Ri）max、 （Ri）min——陨石坑边缘上取得的样点相对最小区域圆（MZC）的最大、最小偏离值；Ri——陨石坑边缘上取得的样点相对最小区域圆（MZC）的偏离值。同样将最小二乘圆（LSC）、最小外接圆（MCC）、最大内接圆（MIC）作为评定的基准圆求得的圆度误差分别为ΔZs、ΔZc、ΔZi。

2.1 圆度误差的数学模型

根据圆的标准方程，最小区域圆、最小二乘圆、最小外接圆和最大内接圆的一般方程见式（3）

（1）最小二乘圆的标准方程见式（3），在陨石坑的边缘上取得n （n ＞3） 个样点，各个点的坐标记为xi、yi（i＝1，2，…，n） ，则各个点到最小二乘圆的圆心的距离为

根据最小二乘圆法构造的目标函数见式（5），该函数的实质是求解圆心的坐标xc、yc和半径Rc的最优化问题，求解结果、、即为最小二乘圆的圆心坐标值和半径值

（2）最小区域圆的标准方程见式（3），陨石坑的边缘上各个样点到最小区域圆的圆心的距离为

根据最小区域圆法构造的目标函数见式（7），该函数的实质是求解圆心的坐标xc、yc的最优化问题，求解结果、即为最小区域圆的圆心坐标值

式中： ［Ri］max、 ［Ri］min——陨石坑边缘上取得的样点与圆心距离的最大值和最小值。

（3）最小外接圆的标准方程见式（3），陨石坑的边缘上各个样点到最小外接圆的圆心的距离为Ri见式（6）。根据最小外接圆法构造的目标函数见式（8），该函数的实质是求解圆心的坐标xc、yc的最优化问题，求解结果、即为最小外接圆的圆心坐标值

（4）最大内接圆的标准方程见式（3），陨石坑的边缘上各个样点到最大内接圆的圆心的距离为Ri见式（6）。根据最大内接圆法构造的目标函数见式（9），该函数的实质是求解圆心的坐标xc、yc的最优化问题，求解结果、即为最大内接圆的圆心坐标值

2.2 圆度误差的求解

圆度误差的求解过程为：

（1）在陨石坑的边缘上取得的n （n ＞3） 个样点，各个点的坐标记为xi、yi（i＝1，2，…，n） ，选取其中3个点，分别记为P1（x1，y1）、P1（x2，y21）、P1（x3，y3），将它们带入方程组 （10），可以得到圆心坐标初值x0、y0和半径初值R0

（2）将求到的圆心坐标初值x0、y0和半径初值R0代入构造函数式（5），同时将圆心坐标初值x0、y0分别代入构造函数式 （7）、式 （8）、式 （9），然后进行迭代求解，分别得到满足迭代求解精度的函数解、、和、，即最小二乘圆的圆心坐标值、半径值和最小区域圆、最小外接圆、最大内接圆的圆心坐标值。

（4）分别将求得的最小二乘圆法、最小区域圆法、最小外接圆法、最大内接圆法相对应的Ri（i＝1，2，…，n） 的值代入式（2），即可分别得到圆度误差ΔZs、ΔZZ、ΔZc、ΔZi。

3 基于BP网络的边缘圆度等级判定

上述的4种圆度评定方法由于使用不同的理想圆，评定结果会有所差别，考虑到月表的复杂性和评价结果的局部误差，只对4种结果进行加权求和做归一化处理得到的圆度的等级缺乏准确性。如何建立一种标准的圆度等级划分模型对检测样本进行判定，本文借助于BP神经网络良好的复杂非线性问题处理和纠错能力，设计一种适用于月表情况BP神经网络判定算法［8，9］。

3.1 建立学习样本

“嫦娥一号”CCD 相机拍摄的遥感影像在月球的中低纬度地区的成像质量更高，本文的实验对象选择这个区域。首先选取10个标准陨石坑为样本，样本的标准为陨石坑周围地形起伏不大，坑唇较窄，坑底较为平缓，然后对这些样本进行人为的模糊处理，划分成圆度好、一般、差3 个圆度等级。通过4种评价方法分别对样本进行评价，再将不同圆度等级的样本评价结果通过线性内插为10 个评价值，对应相应的0－9级。然后将得到的100个评价值用以建立学习样本S＝s i，（ ）j4×100，神经网络的输入数据应进行归一化，再把神经网络的期望输出设置为圆度Y ＝［0，0.1…，0.9，…，0，0.1，…，0.9］1×100。传统的最大最小值法易使数据归一化为传递函数Sigmoid的极大值0和极小值1，导致连接权值应足够大才能使神经网络的输出值与之相匹配，进而使训练次数增多和训练时间变长。为了避免这个现象，可以把区间数据规范到［0.15，0.85］区间，本文使用改进型归一化公式，见式（11），最小二乘圆法评价作为逆指标，归一化公式见式（12）

3.2 神经网络训练

BP神经网络分为三层，分别为输入层、中间层（或隐含层）、输出层［10］。其中，输入层X1、X2、X3、X4分别是4种评价算法的值，中间层（或隐含层）的各个神经元通过Wij（i＝1，2，…，n；j＝1，2，…，p） 和输入层建立权重连接，通过Vjk（j＝1，2，…，p；k＝1，2，…，q）和输出层建立连接，输出层Yi（i＝0，1，…，9）是最终的评级结果，如图4所示。

从图4可知，神经网络中各层单元都与上一层或下一层建立连接，但是各层之间并无连接。神经网络学习的基本步骤为：

（1）正向传播，将输入层输入的学习样本通过权值连接送入中间层，经中间层的传递函数激活，再通过连接权值送入到输出层，经输出层的传递函数激活，的得到输出结果；

（2）反向调整，BP 神经网络是一种前馈型神经网络，定义了一种误差函数E，若E 不满足预定精度，则需要不断调整连接权值和阀值减少误差；

（3）循环调整，调整连接权值和阀值，直到误差函数E满足要求且需要设置迭代上限避免学习陷入死循环。

图4 神经网络

如图5所示，其中，图5（a）是神经网络的误差变化曲线，误差曲线表明训练了972步后，输出误差小于误差上限0.001，图5（b）是神经误差的误差检验结果图，是为了检验神经网络的准确性，再次把学习样本输入网络进行检验，得到的散点图，散点图表明圆度的评价等级检验的绝对误差都小于0.1，因而神经网络的准确性得到了保证。

图5 神经网络的训练结果

3.3 评价结果展示

本文的实验选取了30个陨石坑，为了使圆度判定的结果与主观的判定相互印证，从中又选取了10幅能够主观上判断圆度等级的影像在本文中展示，相关图像和判定数据见图6和表1。表1中的神经网络测试结果对应圆度判定的等级，根据学习样本和测试集的神经网络测试结果，神经网络测试结果的值与边缘圆度判定等级划分的对应表见表2。

图6 测试图片集 （按降序圆度减小）

最后，将本文设计的边缘圆度的评价模型加载到数字月球平台，评价模型可以对月表的陨石坑进行边缘圆度判定，结果展示如图7所示。

表2 圆度等级划分的数值区间

图7 边缘圆度评价示例

4 结束语

本文基于圆度误差的理论和求解方法，提出了一种月表陨石坑边缘圆度评价算法，通过BP神经网络综合了最小二乘圆法、最小区域圆法、最小外接圆法和最大内接圆法4种算法的评价结果，进而得到边缘圆度等级。通过边缘圆度算法建立评价模型，并且对实验样本进行了评价，实验结果和人工对比结果基本一致。本文的边缘圆度评价结果最终可以和其它5个判定因子一起，建立了全月性的月表的陨石坑几何影像特征库，为研究月表地质和探月工程的后续工作提供了一个研究平台。由于运用的 “嫦娥一号”全月球影像数据分辨率为120m，不能精确的对中小型和小型的陨石坑的进行边缘圆度等级判定，未来的研究工作是利用 “嫦娥二号”7m 分辨率影像数据，根据陨石坑大小的级别，分别建立陨石坑的边缘圆度等级判定模型。

［1］LI Kun，LIU Jianjun，MOU Lingli，et al.Crater retention age calculation of lunar geological unit using small－scale craters［J］.Journal of Jilin University（Earth Science Edition），2012，42 （2）：452－459 （in Chinese）.［李坤，刘建军，牟伶俐，等.利用小型撞击坑测算月球地质单元撞击年龄 ［J］.吉林大学学报地球科学版，2012，42 （2）：452－459.］

［2］Morota T，Haruyama J，miyamoto H，et al.Formation age of the lunar crater Giordano Bruno［J］.Meteorities ＆Planetary Science，2009，44 （8）：1115－1120.

［3］Melosh H J.Planetary surface processes ［M］.Cambridge University Press，2011.

［4］FANG Qinlin.Research on roundness evaluation algorithm and programming ［D］.Wuhan：Huazhong University of Science and Technology，2007 （in Chinese）.［方沁林.圆度误差评定的算法研究与软件设计 ［D］.武汉：华中科技大学，2007.］

［5］FENG Junhua，CUI Hutao，CUI Pingyuan，et al.Autonomous crater detection and matching on planetary surface ［J］.Chinese Journal of Aeronautics，2010，31 （9）：1858－1863 （in Chinese）.［冯军华，崔祜涛，崔平远，等.行星表面陨石坑检测与匹配方法［J］.航空学报，2010，31 （9）：1858－1863.］

［6］ZHANG Yumei.Research on roundness on－machine measurement based on multiscale analysis［D］.Changchun：Jilin University，2009 （in Chinese）.［张玉梅.基于多尺度分析的圆度误差在线检测研究 ［D］.长春：吉林大学，2009.］

［7］ZHANG Chunyang，LEI Xianqing，LI Jishun，et al.Method for roundness error evaluation based on geometry optimization［J］.Journal of Mechanical Engineering，2010，46 （12）：8－12（in Chinese）.［张春阳，雷贤卿，李济顺，等.基于几何优化的圆度误差评定算法 ［J］.机械工程学报，2010，46 （12）：8－12.］

［8］WEN Chao.Research on image－based cultural relic retrieval and ontology－based annotation ［D］.Xi’an：Northwest University，2012 （in Chinese）.［温超.基于图像的文物检索与本体标注方法研究 ［D］.西安：西北大学，2012.］

［9］WAN Dingsheng，HU Yuting，REN Xiang.Applied research of BP network with feedback input［J］.Computer Engineering and Design，2010，31 （2）：398－400 （in Chinese）.［万定生，胡玉婷，任翔.带反馈输入BP 神经网络的应用研究 ［J］.计算机工程与设计，2010，31 （2）：398－400.］

［10］WANG Qianqian，HUANG Zhiwen，LIU Kai，et al.Classification of plastics with laser－induced breakdown spectroscopy based on principal component analysis and artificial neural network model［J］.Spectroscopy and Spectral Analysis，2012，32 （12）：3179－3182 （in Chinese）.［王茜倩，黄志文，刘凯，等.基于主成分和人工神经网络的激光诱导击穿光谱塑料分类识别方法研究 ［J］.光谱学与光谱分析，2012，32 （12）：3179－3182.］