0 引言
陨石坑作为月球表面最典型且普遍的地貌单元和地质结构，已成为月球科学研究最多的地貌特征之一[1]。对月表撞击坑研究有助于加深对月表形貌的认识和理解，对探测器软着陆地点的修正具有重要意义[2]。由于陨石坑形态复杂多样，而人工检测非常耗时且难以处理大型数据集，因此需使用计算机来快速检测陨石坑[3]。目前应用较多的陨石坑检测方法主要有两种类型：基于图像特征的传统方法和基于机器学习的方法[4]。
基于图像特征的传统方法通常会对原始图像进行预处理以增强坑的边缘，并使用霍夫变换、椭圆拟合、遗传算法、Gist特征、分水岭变换、模式识别、径向一致性算法或上述方法的组合，将图像中的环形坑边缘识别为圆形或椭圆形特征。基于机器学习的方法主要使用可扩展模板模型、卷积神经网络、自适应增强和改进自适应增强方法、最小二乘支持向量机等[5]。Kang等人提出了一种使用定向梯度(HOG)特征直方图和支持向量机(SVM)分类器自动提取小规模撞击坑的方法[6]。PyCDA是一种基于神经网络的陨石坑检测器，由探测器、提取器和分类器组成。Silburt等人用U-Net网络对月球DEM图像数据集训练并输出预测目标,再用模板匹配算法提取出陨石坑的位置和直径信息 [7]。
然而，目前基于图像特征的方法虽然处理过程简单，但检测精度不高且不适用于多种尺寸陨石坑同时检测的情况；而基于机器学习的检测方法则包含繁杂的后期处理步骤，大大降低了推理速度，很难区分和训练。
为了更加快速和精确的实现对月表图像中陨石坑检测，本文提出一种用于检测多尺度月球陨石坑的新方法(C-Moon-Net)。该方法基于CenterNet网络[8]检测月球DEM图像数据集中的陨石坑，具有检测精度高和速度快的优势。首先，利用月球DEM图像制作数据集，并对数据集进行预处理，以适应深度学习网络。主干网络选择DLA(Deep Layer Aggregation)网络结构[9]，并使用TTFNet中提出的基本损失计算方法对数据集进行训练[10]。
1 数 据
本文使用与Silburt等人[7]相似的方法，选用的图像数据来自月球侦察轨道器(LRO)[11]和Kaguya合并数字高程模型的数字高程(DEM)图像[12]，如图1所示。该模型跨越纬度范围为±60°，经度范围是±180°，分辨率为512像素/度(59 米/像素)。该全球灰度图是简单圆柱投影，分辨率为184 320×61 440像素，位深度为16位/像素；本文将其下采样为92 160×30 720像素和8位/像素，以减少计算量。本文使用高程图像而不是光学图像，是因为在高程图像中陨石坑的外观不受入射阳光方向的影响。这样可以减少陨石坑之间的外观变化，从而更容易训练模型检测陨石坑。
图1 LRO和Kaguya合并数字高程模型的数字高程图像
本文数据集中每个DEM图像都是通过以下方式生成的：1)从图1全月数字高程图中裁剪出大量正方形图像，裁剪区域位置随机选择且分布均匀。2)将裁剪后的图像使用最近邻元法降采样为256×256像素，以降低运算量、增加感受野。3)使用Cartopy Python软件包将图像转换为正交投影[7]。原图像采用的简单圆柱投影越往两极变形越大，圆形的撞击坑在该投影中呈现出的是越来越扁的椭圆形，对后续检测的影响较大，而正射投影更接近陨石坑的真实形状，可以最小化图像失真。
对于每个输入图像，都会生成一个相应的地面真实目标，该目标蒙版也为256×256 像素。陨石坑在目标蒙版中编码为厚度为1像素的环，半径和中心取自坑在陨石坑目录中的物理位置和直径。用于构建真实目标的数据是通过合并两个学者建立的陨石坑目录获得的。对于5～20 km的陨石坑，使用Povilaitis等人收集的全球陨石坑数据集[13]；对于大于20 km的陨石坑，由使用Head等人发表的全月球陨石坑数据集[14]。将上述裁剪图像和真实目标等分成两个单独的数据集，用于训练和测试模型。这两个数据集是从月球的大小相等且互斥的部分采样的。训练集为跨越经度-180°至-60°的区域，测试集为经度60°至180°的区域。每个数据集包含3 000个DEM图像，每个DEM图像的坑口中位数为21。由于本文采用的人工生成的陨石坑数据集不完整，包含许多明显的遗漏标识。Fassett等人估计[15]，Head等人的数据集不完整性为12%。Povilaitis等人的数据集不完整性目前未知，但人工检验发现其不完整性更高。数据集不完备会导致错误的训练结果。为了解决真实标注不完备的问题，本文人工选取了数据集内真实标注最为完整、边缘更为准确的图像各3 000张，分别作为模型的训练集和测试集。
同时，为了适应CenterNet网络，将真实目标标注形状由环状转换为圆外切矩形，并记录矩形框的左上角和右下角位置坐标。如图2所示，图中红色圆形为陨石坑环形标注，黑色方框为圆外切矩形框。矩形框左上角和右下角坐标可由式(1)得到：
图2 真实目标标注形状转换图
(xmin≥0,ymin≥0)
(xmax≤256,ymax≤256)
(1)
其中，圆心坐标为(x0,y0)，直径为D；矩形框的左上角坐标为(xmin,ymin)，右下角坐标为(xmax,ymax)。
数据集内DEM图像、环形真实目标、转换后的数据标注框的示例如图3所示。该示例图的经度范围为(158.742 187 5, 162.410 156 25)，纬度范围为(14.882 812 5,18.550 781 25)。该示例图中陨石坑真实目标目录如表1所示，第1列表示目标的序号；第2～4列展示了目标在全月DEM图像中的中心经纬度坐标和直径(km)；第5～7列展示了目标在裁剪图像中的圆心坐标(x0,y0)和直径(像素)；第8～11列展示了目标由环形标注转换成矩形框标注后， 矩形框的左上角坐标(xmin,ymin)和右下角坐标(xmax,ymax)。
表1 示例图内陨石坑目标目录
序号直径/km纬度/(°)经度/(°)x0/pixy0/pix直径/pixxmin/pixymin/pixxmax/pixymax/pix08.6315.31161.67201.30225.7819.82191.40215.87211.21211.21149.2615.53160.29108.08211.10113.1451.51154.53164.65164.65235.3916.11161.75206.13170.6181.28165.49129.96246.77246.77338.1916.36159.6364.20153.1687.7120.35109.31108.05108.05413.0616.89161.86212.86115.7829.99197.86100.78227.85227.8558.0317.06159.9082.67104.2918.4473.4595.0791.8991.8967.9317.37161.44185.0282.5218.21175.9173.42194.12194.12714.4317.48162.13231.0574.9733.13214.4858.41247.62247.6285.3517.49158.798.7773.8612.292.6367.7114.9214.9299.7818.09159.2036.0432.1722.4624.8220.9547.2747.27106.5018.21161.69201.3823.8514.92193.9216.39208.84208.841115.4518.24160.93151.2321.8835.48133.494.14168.97168.97
(a) 输入DEM图像
(b) 目标蒙版图
(c) 矩形框标注图
图3 输入图像与真实标注示例图
2 C-Moon-Net算法
本文采用简化的CenterNet作为基础卷积神经网络(CNN)模型。CenterNet算法是Zhou等人在2019年4月提出的。它是在构建模型时将检测目标表示为一个点——目标边界框的中心点[8]。然后，其他属性(例如目标尺寸、3D范围、方向和姿势)直接从中心点的图像特征中回归得到。此时目标检测转换为标准的关键点估计问题。CenterNet简单地将输入图像馈送到生成热力图的完全卷积网络。此热力图中的峰值对应于目标中心。每个峰值处的图像特征可预测目标边界框的高度和宽度。该模型使用标准的密集监督学习进行训练。推理是单个前向预测网络，后期处理没有非最大值抑制。
首先将图片输入到主干网络中，生成一张热力图，热力图上的点与输出的预测目标框有关。输入图像和热力图如图4所示。然后将热力图上的点对应的图像特征输入预测网络中，直接预测目标框的尺寸大小，从而将检测问题简化为目标中心的预测。本文的C-Moon-Net使用DLA-34分割网络作为它的主干网络，该网络通过迭代深度聚合的特征融合方式从而实现多尺度目标的鲁棒性检测[9]。此外，使用TTFNet中建议的基本损失计算方法进行训练[10]，与CenterNet相比，TTFNet监督信号更多，模型的收敛速度更快。C-Moon-Net架构图如图5所示。假设输入图像大小为X×Y，网络首先对图像提取特征，并生成尺寸为
的特征图。然后，生成的特征图分别输入到两个卷积层，用于目标定位和边界框回归[16]。
(a) 输入图像
(b) 热力图
图4 输入图像与输出热力图对比
图5 本文算法架构图
第一个卷积层的卷积核大小为1×1，通道数为1，后面跟一个sigmoid层。令
表示这个sigmoid层的输出，其中
代表像素(i,j)属于一个目标中心点的概率。
第二个卷积层的卷积核的大小是1×1，通道数是4。假设第二个卷积层的输出为
如果像素(i,j)被认定为一个目标的中心点，
代表像素(i,j)到边界框4个边之间的距离。
在训练阶段，定位损失是用修改后的focal损失来计算的，回归损失用L1损失计算。
2.1 目标定位
给定第m个真实标签边界框，首先将其线性映射到特征图尺寸下。其中，下采样倍数R=4。然后采样二维高斯核生成
二维高斯核如式(2)所示：
(2)
式中，
是原始图像上第m个框的中心点，(x0,y0)代表特征图尺寸下第m个框的中心点，可以变换为(
。(ω,h)是特征图尺寸下第m个框的宽度和高度。最后，本文通过对Hm应用元素最大值来生成真值图 H。
高斯分布的峰值，也就是标注框中心点，认为是正目标，其他所有像素被视为负目标。使用修正的focal 损失[10]计算CenterNet中的定位损失。给定预测和真值H，定位损失Lloc如式(3)所示：
(3)
式中，N为训练图像中的目标数量。
2.2 尺寸回归
给定特征图尺寸下的第m个真实边界框，真实标签的回归定义为
给定在特征图尺寸下第m个标注框内的像素点(i,j)，回归目标Sij定义为从(i,j)到第m个标注框的四边的距离，归一化系数为4，表示为一个四维矢量(ωl,ht,ωr,hb)m。因此，原始图像在(i,j)处的预测框(x1,y1,x2,y2)可以表示为
x1=4i-16ωl,y1=4i-16htx2=4i+16ωr,y2=4i+16hb
(4)
将Hij>0的像素区域定义为Am，也称为子区域，将子区域中每个像素作为回归样本，回归损失定义为式(5)：
(5)
式中Nreg是Hij>0区域的像素数量，即回归样本的数量。Wij是样本权重，用于对不同尺寸下的边界框平衡损失，从而不会影响本文后续的讨论。
由于大尺寸图像的目标多样，大的目标可能产生成千上万个样本，而小目标可能只有几个样本。在归一化所有样本的损失后，小目标的损失可以忽略不计，这将影响小目标的检测性能。因此，样本权重Wij在平衡损失中发挥了重要作用。
假设(i,j)是子区域Am的第m个标记框，则Wij如式(6)所示：
Wij=
(6)
式中，Gm(i,j)是在(i,j)处产生的高斯概率，am是第m个框的面积。这样可以利用大目标中包含的标记信息，并保留小目标的标记信息。同时，它还能着重强调靠近目标中心的样本，减少模糊样本和低质量样本的影响。
最终的总损失可以表示为
(7)
本文中设定ωloc=1.0，ωreg=5.0。
在测试阶段，在定位分支输出的特征图中的峰值点对应的像素作为预测目标边框的中心点，而其他的像素的输出会被丢弃。与传统的锚点检测模型(如RetinaNet或Faster-RCNN)相比，本文采用的这种简化的无锚点模型有更快的测试速度(由于其推理速度更快)。
3 实验结果分析
3.1 网络训练
本实验使用的月表陨石坑数据集如第1节所述。训练集和验证集各由3 000张256×256像素的图像组成。在测试阶段，最大池化层kernel大小为3×3，用于提取定位分支中输出特征图的峰值点，并且所有峰值点的峰值只有大于Tconf才被认定为正目标。在本实验中Tconf设置为0.05。对应于回归分支输出结果的峰值点，通过公式(4)用于计算目标边界框的位置。
本实验基于PyTorch深度学习框架，在Anaconda上搭建PYTHON 3.6.9，PyTorch 1.3.0的虚拟环境。使用第2节所述的网络模型，所有训练和测试均在Windows 10系统、Tesla P100-PCIE-12GB显卡、CUDA10.1上进行。所有模型都经过随机梯度下降(SGD)算法训练。经过实验，在每次迭代中，设置batch size为32，模型精度和训练速度达到最优。所有模型都训练了150个epoch。实验采用余弦退火(cosine annealing)学习率调度策略，初始学习率设置为0.001。
3.2 模型评估标准
为了定量评估本文方法的有效性，使用标准PASCAL VOC评价指标[17]评估C-Moon-Net性能。
对于基于CNN的检测模型，使用特定的IoU(Intersection-over-Union )阈值，以低置信度筛选出检测结果。当阈值增加时，模型精度Pr会随之增加，而模型召回率Rr会随之减小。召回率Rr反映了所有真正为正例(positive targets)的样本中被分类判定为正例的比例，其定义如式(8)：
(8)
精度反映了被分类判定的正例中真正的正例样本的比例，其定义如式(9)：
(9)
其中，TP(True Positive)、FP(False Positive)，TN(True Negative)，FN(False Negative)的含义如表2所示。
表 2 FN、FP、TN、TP含义
指标真实标签分类结果FN正样本负样本FP负样本正样本TN负样本负样本TP正样本正样本
AP是目标检测算法的标准指标，它综合考虑了在不同置信度级别下模型的精度Pr和召回率Rr，AP定义如式(10)所示：
(10)
一个理想检测器的AP=1。本实验使用4种AP指标来评估结果，包括mAP，AP30，AP50和AP75。
3.3 实验结果
经过150轮的训练，损失基本达到平稳状态。表3是150轮训练后在测试集上的平均精度，IoU为0.5时表示整个网络的mAP，其值为67.53%。AP_small，AP_medium，AP_large分别表示像素面积小于322的小目标框、像素面积在322～962之间的中等目标框、像素面积大于962的大目标框的AP测量值。通过观察序号5～7和序号8～10发现，中等目标的检测效果最好，其次是大目标，最后是小目标。
表3 测试集AP指标
序号AP说明结果/%1Average Precision(AP)@[IoU=0.50:0.95 area=all]35.02672AP30@[IoU=0.30 area=all]74.62143AP50@[IoU=0.50 area=all]67.53104AP75@[IoU=0.75 area=all]29.58845AP30_small@[IoU=0.30 area=small]69.72766AP30_medium@[IoU=0.30 area= medium]80.51607AP30_large@[IoU=0.30 area= large]71.17668AP50_small@[IoU=0.50 area=small]61.82599AP50_medium@[IoU=0.50 area= medium]78.948810AP50_large@[IoU=0.50 area= large]69.7873
Silburt等人提出的DeepMoon方法也是基于神经网络对月球DEM图像数据集训练并测试。该方法的平均精度仅有56%，而本文的C-Moon-Net用AP50表示网络的平均精度为67.53%，比DeepMoon的精度上升11.53%。此外，DeepMoon无法可靠检测直径大于15个像素的陨石坑，仅限于像素半径小于15像素的坑[7]。而本文网络对多种尺度的陨石坑均可准确识别，并且对大目标的检测精度达到69.79%，远远高于DeepMoon。本文还与另一种基于神经网络的陨石坑检测器PyCDA进行了比较。该检测方法的精度为25%，只能检测出较少数量的大直径陨石坑。与之相比，本文检测方法精度有大幅提升，且对中小尺寸陨石坑检测结果更为突出，平均精度分别达到78.95%和61.83%。
众所周知，在图像中检测出的陨石坑坑数量越多、图像特征数量越多，对应于地形导航中位置估计的不确定性也就越低。图6是测试数据集的示例和真实标签图。图7展示了在同一DEM图像上的3个陨石坑检测方法的检测性能结果比较。这些图像显示了已检测到并成功与真实标签匹配的陨石坑。图7(a)是本文方法检测结果图，绿色框是真实标签框(ground truth)，红色框是检测结果框，蓝色框是FN假负例，百分比代表置信度。图7(b)是DeepMoon方法检测结果图，蓝色圆形即检测结果。图7(c)是PyCDA检测结果，红色圆形即检测结果。由结果可知，PyCDA对于DEM图像的检测效果最差，检测到的陨石坑数量最少，且只能检测出大直径陨石坑；DeepMoon检测的坑不完全，且只能检测中小尺寸的坑；而本文网络不仅能检测到各种尺寸的陨石坑，而且检测到的陨石坑数量最多。
(a) 输入图像
(b) 真实标签
图6 测试数据集
(a) 本文检测结果
(b) Deep-Moon检测结果
(c) PyCDA检测结果
图7 3种陨石坑检测器检测结果对比图
DeepMoon、PyCDA和本文方法C-Moon-Net均在Tesla P100-PCIE-12GB上的Python Keras中进行了测试。图8柱形图比较了3种陨石坑检测方法在测试集上对每张图片的平均检测时间。与其他基于神经网络的陨石坑探测方法相比，C-Moon-Net具有最短的检测时间，在0.01 s的级别。如表4汇总了3种检测方法的精度和预测时间比较。由实验结果可知，本文的陨石坑检测方法C-Moon-Net相较于DeepMoon检测速度提高了201倍，检测精度提升了11.5%，减少了网络的漏检率，且极大地提高了大目标的检测率。此外，与PyCDA检测方法比较，本文检测精度提升了42.5%，检测时间降低了95.8%，且对中小尺寸的陨石坑检测效果有显著提升。
图8 3种陨石坑检测方法对测试集每张图片的平均检测时间
表4 3种陨石坑检测方法检测时间和精度比较
检测方法预测时间/s训练时间/h平均精度/%C-Moon-Net0.0152.5867.5DeepMoon3.0276.4156PyCDA0.365.3625
各陨石坑探测方法处理月球图像和提取陨石坑所花费的时间和模型精度有较大差异，是由于其网络结构以及处理过程不尽相同。本文采用的CenterNet是利用关键点三元组来感知物体内部信息，使用DLA-34分割网络提取特征，该网络模型的深度较高，提取的特征非常多。其次，本文的检测方法采用的是端到端的网络模型，不需要提前设置锚框的超参数，没有非最大值抑制，更简洁更快速更精确，在速度和精度上达到了最好的权衡。而DeepMoon方法使用的是U-net网络。在网络结构方面，该网络包含4个卷积层和4个上采样层，它将每个卷积层得到的特征图融合到对应的上采样层[7]，网络模型深度不高，提取的特征数量远不如本文网络。在处理过程方面，该方法通过模板匹配来处理神经网络预测。模板匹配需要循环扫过所有可能的圆半径，从而延长了其运行时间。PyCDA方法由两个U-net神经网络组成，而不仅仅是一个神经网络，这也大大增加了它的运行时间。
4 结束语
现有的月球陨石坑检测方法需要繁杂的后期处理，难以满足实时陨石坑检测需求。针对该问题，本文提出一种用于检测月球多尺度陨石坑的新方法，C-Moon-Net，能够实现端到端实时检测。
1)与使用深度学习检测陨石坑的DeepMoon方法相比，本文所提方法的检测速度提升201倍，检测精度提升11.5%；本文方法检测速度与精度也优于其他基于神经网络的陨石坑探测方法。
2)本文所提方法首次采用简化的CenterNet网络作为训练模型，没有耗时且计算负担大的后期处理。该方法极大地提升了检测速度，以满足实时陨石坑检测的需求。
3)本文方法可以精确地检测各种尺寸的陨石坑并输出唯一的检测框，有效地克服了现有检测方法需要筛选过滤重复陨石坑目标的问题，提高了检测的准确率，具有强鲁棒性的特点。
在后期的研究工作中还需要降低位置估计误差并对边缘提取算法进行研究，以实现高精度的陨石坑检测。
[1] 李春来，刘建军，左维，等．中国月球探测进展(2011-2020年)[J]．空间科学学报，2021，41(1):68-75．
[2] 田亚骏，张明，林轻．月面大范围探测功能需求分析及其研究现状[J]．载人航天，2020，26(5):649-655．
[3] 程维明，刘樯漪，王娇，等．全月球形貌类型分类方法初探[J]．地球科学进展，2018，33(9):885-897．
[4] 王赤，张贤国，徐欣锋，等. 中国月球及深空空间环境探测[J]. 深空探测学报，2019，6(2)：105-118.
[5] WANG Y，WU B． Active Machine Learning Approach for Crater Detection from Planetary Imagery and Digital Elevation Models[J]. IEEE Trans on Geoscience and Remote Sensing, 2019,57(8):5777-5789．
[6] KANG Zhizhong, WANG Xingkun, HU Teng, et al. Coarse-to-Fine Extraction of Small-Scale Lunar Impact Craters from the CCD Images of the Chang’E Lunar Orbiters[J]. IEEE Trans on Geoscience and Remote Sensing, 2019,57(1):181-193.
[7] SILBURT A, ALI-DIB M, ZHU C C , et al. Lunar Crater Identification via Deep Learning[J]. ICARUS, 2019,317(1):27-38.
[8] ZHOU X Y, WANG D Q, KR
HENBÜHL P. Objects as Points [EB/OL]. 2019: arXiv: 1904.07850[cs.CV].https://arxiv.org/abs/1904.07850.
[9] YU F, WANG D Q, SHELHAMER E, et al. Deep Layer Aggregation [C]∥Proceedings of 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition, Salt Lake City, UT, USA:IEEE, 2018:2403-2412.
[10] LIU Z L, ZHENG T, XU G D, et al. Training-Time-Friendly Network for Real-Time Object Detection[J]. Proceedings of the AAAI Conference on Artificial Intelligence, 2020,34(7):11685-11692.
[11] CHIN G, BRYLOW S, FOOTE M, et al. Lunar Reconnaissance Orbiter Overview: The Instrument Suite and Mission [J]. Space Science Reviews, 2007, 129(4):391-419.
[12] KATO M, SASAKI S, TANAKA K, et al. The Japanese Lunar Mission SELENE: Science Goals and Present Status [J]. Advances in Space Research, 2008, 42(2):294-300.
[13] POVILAITIS R Z, ROBINSON M S, VAN DER BOGERT C H, et al. Crater Density Differences: Exploring Regional Resurfacing, Secondary Crater Populations, and Crater Saturation Equilibrium on the Moon [J]. Planetary and Space Science, 2018, 162:41-51.
[14] HEAD J W, FASSETT C I, KADISH S J, et al. Global Distribution of Large Lunar Craters: Implications for Resurfacing and Impactor Populations [J]. Science, 2010, 329(5998):1504-1507.
[15] FASSETT C I, HEAD J W, KADISH S J, et al. Lunar Impact Basins: Stratigraphy, Sequence and Ages from Superposed Impact Crater Populations Measured from Lunar Orbiter Laser Altimeter (LOLA) Data [J]. Journal of Geophysical Research: Planets, 2012, 117(E12):E00H06.
[16] YANG R, WANG R, DENG Y K, et al. Rethinking the Random Cropping Data Augmentation Method Used in the Training of CNN-Based SAR Image Ship Detector [J]. Remote Sensing, 2020, 13(1):34.
[17] EVERINGHAM M, WILLIAMS J,WINN J,et al. The Pascal Visual Object Classes Challenge 2007 (voc2007) Results[EB/OL].[2020-04-29].http://host.robots.ox.ac.uk/pascal/VOC/voc2007/.
