基于Apollo图像的月表撞击坑自动提取

李 超，王心源，骆 磊，吉 玮

(1.中国科学院对地观测与数字地球科学中心数字地球重点实验室，北京 100094;2.中国科学院研究生院，北京 100049)

基于Apollo图像的月表撞击坑自动提取

李 超1，2，王心源1，骆 磊1，2，吉 玮1

(1.中国科学院对地观测与数字地球科学中心数字地球重点实验室，北京 100094;2.中国科学院研究生院，北京 100049)

撞击坑是月表的主要形貌单元。研究月表撞击坑的形状与分布特征对描述月表撞击过程与演变规律有很大帮助。根据月表撞击坑的真实形状，提出了一套椭圆形撞击坑自动提取方法，其主要包括3个步骤:图像预处理、边缘检测与约束参数的Hough多椭圆分层检测。与人工提取结果比较，月表撞击坑自动提取的总体准确率在76%以上(其中比较大的撞击坑的检测率可以达到83%)。统计检测结果表明，当撞击坑的长半轴D＜160 m时，其累积大小频率分布曲线撞击坑个数N与D－3.8有明显的线性关系。

月表撞击坑;自动提取;Hough多椭圆变换;边缘

0 引言

撞击坑(craters)对于研究类地行星(卫星)的形貌、相对地质年龄、岩石构造以及行星起源有着重要意义。Shoemaker等［1］认为，在行星表面上的撞击坑数量可以反映不同年代的撞击现象。这个原理是基于撞击坑在空间与时间随机形成的假设，常用以评价星体表面形成的绝对或相对年龄。这种理论是通过统计撞击坑的累积大小频率分布(sizefrequency distribution，SFD)来体现的。以前的学者［2－4］大部分是通过人工方法提取撞击坑进行分析的，但由于这种方法具有耗时长、缺乏客观性、对操作者专业知识要求高等弱点，使自动提取撞击坑的要求逐渐强烈起来。近年来，许多学者从不同方面提出撞击坑自动检测方法［5－9］，形成了各自的特征识别系统，包括模式匹配、纹理分析、基因算法与面向对象等多种提取方法。这些方法都是将撞击坑抽象为圆形。但是，由于特定的撞击角度以及传感器拍摄视角等原因，撞击坑的形状更符合椭圆形。本文利用Apollo 15图像数据，根据月表撞击坑的真实形状，选取受到3个大型撞击坑(哥白尼(Copernicus)、阿里斯塔克(Aristarchus)和开普勒(Kepler)撞击坑)强烈影响、内部有很多次生撞击坑(secondary craters)的区域作为研究区，旨在提出一种新的椭圆形检测方法，为提取与分析月表次生撞击坑提供方法支持。

1 研究区与数据源

如图1所示，研究区位于月海内的欧拉(Euler)撞击坑与丢番图(Diophantus)撞击坑之间。研究所用数据为 Apollo 15 图像数据(http://wms.lroc.asu.edu/apollo/browse?camera=M＆mission=15)。该图像编号为M－1845，采用Mapping/Metric相机近垂直角度拍摄，拍摄高度为102.99 km，太阳高度角为7°;图像空间分辨率为6.909 m，中心坐标为(N 25.73°，W 37.97°)，跨度为160 km ×160 km。

图1为采用GCS_Moon_2000空间参考系正射纠正后的Apollo图像。

图1 正射纠正后的Apollo图像Fig.1 Apollo image after orthographic projection correction

2 撞击坑检测方法

研究区内撞击坑特征十分复杂，为了便于计算机处理，在研究区内选择6个小范围的、具有便于识别撞击坑特征的感兴趣区(area of interest，AOI)(图1中红色矩形标识的区域)。在获取图像中的AOI之后，通过图2所示的3个步骤自动提取撞击坑。

图2 撞击坑自动检测流程图Fig.2 Flow chart of automatic detection of craters

2.1 数据预处理

根据图像局部灰度值的统计数据对选取的AOI进行增强处理，增大图像的对比度，有选择地突出图像中感兴趣的撞击坑轮廓特征，同时抑制(掩盖)图像中其他信息，使图像与计算机或人的视觉响应特性相匹配。

为消除图像噪声对月表撞击坑检测的影响，采用非线性、非迭代的双边滤波(bilateral filtering)方法［10］对图像进行处理。这种双边滤波在处理相邻像素时，不仅考虑到距离上的邻近关系，同时还考虑到灰度上的相似性，通过对这二者的非线性组合，可自适应地进行滤波。这样图像在滤除噪声的同时，也能很好地保留图像的边缘信息。Apollo原始图像与预处理结果如图3所示。

图3 原始图像与预处理结果Fig.3 Original image and preprocessing result

2.2 边缘提取

分析Apollo图像特点及综合比较各种边缘检测方法后，最终采用Canny算法［11］提取月表撞击坑边缘。鉴于传统Canny算子的计算量比较大、实时性较差，并且检测边缘的精度仍没有达到单像素级，同时受到多种因素的影响造成部分伪边缘，为此本文采用改进型Canny边缘检测算法［12］提取月表撞击坑边缘并进行细化，结果如图4(a)所示。由于获取的Apollo图像数据的太阳高度角为7°，会使撞击坑高出的坑壁出现阴影。与Apollo原始图像(图3(a))比较可以看出，提取的月表撞击坑边缘(图4(a))中既包含了撞击坑的真实边缘信息，又包括了因光照条件和地形因素造成的撞击坑内部光亮区与阴影区的伪边缘信息以及其他地形特征的边缘信息。通过检测撞击坑图像的梯度方向［13］，发现在撞击坑内部的伪边缘与撞击坑的真实边缘的梯度方向正好相反，真实撞击坑的边缘梯度方向与光源方向之间的夹角小于90°，而伪边缘梯度方向与光源方向之间的夹角大于90°，即真实边缘满足

式中:g为边缘梯度矢量;n为光源方向矢量;|g|为梯度矢量的模长;|n|为光源方向矢量的模长。

通过上述方法选择边缘梯度方向与光源方向夹角小于90°的边缘，就可以有效剔除伪边缘(图4(b))。最后，按照边缘的方向进行一定邻域内的连接，同时考虑噪声以及下文论述的Hough变换的计算量，通过设置一定的阈值将较小的边缘从中剔除。

图4 边缘检测Fig.4 Edge detection

2.3 椭圆拟合

现有的椭圆拟合算法大致分为5类:①基于经典的Hough变换及其改进算法［14－15］;②有效统计算法(statistically efficient method)［16］; ③遗传算法(genetic algorithm)［17］; ④ 基于椭圆 几何 特性 算法［18－19］;⑤基于广义正交概念的 K随机采样(K－RANSAC)及其改进算法［20］。椭圆的自由参数包括椭圆的中心坐标(x，y)，长、短半轴的长度(D，B)、长轴的偏转角度(θ)。在实际椭圆检测应用中，对于多个椭圆的检测往往存在计算耗时长、内存需求大等问题，因此如何降低椭圆的空间维度以便加快检测速度已成为椭圆检测的主要问题。

Hough变换是一种在图像中定位形状的技术［21］。它是基于“证据收集”方法，通过模板匹配过程的重新描述来实现的(其中，证据是在累加器数组中的投票)。Hough变换的基本思想是利用点—线的对偶性，对图像进行坐标变换，使之在另一个坐标系(累加器空间)中的特定位置上出现峰值，检测曲线的问题就简化为寻找峰值位置的问题。这种变换的优点是受噪声与曲线间断的影响较小，但最大的缺点是计算量非常大，且需要比较大的内存。为此，改进的Hough变换被先后提出，比如概率Hough变换(PHT)、随机Hough变换(RHT)等。但是，这类方法主要针对单个或者简单的椭圆，直接应用到月表撞击坑的提取，效果还不理想。本文结合最小二乘与Hough变换，通过选取AOI，采取基于局部主成分分析(PCA)感兴趣参数约束Hough多椭圆分层检测的思路［22］，方法步骤如下:①设定长度阈值T1与T2(T1＞T2);②在边缘图像中标记并提取长度大于T1的连续曲线段;③对各线段做PCA方向分析，确定是否属于有效曲线段;④对所有感兴趣曲线段按照标记顺序依次利用最小二乘椭圆拟合方法得到感兴趣椭圆粗略参数;⑤根据拟合结果进而模糊约束Hough变换参数搜索范围，得到精确椭圆参数，并输出这5个自由参数;⑥利用检测结果更新图像空间，删除已经检测到的椭圆;并按照一定的步长迭代进行，直到连续曲线段长度小于T2，所有椭圆检测完毕。图5为按照上述步骤得到的6个AOI最终的椭圆检测结果。

图5 AOIs区域的椭圆(＞7个像元)拟合结果Fig.5 Results of the ellipse fitting larger than seven pixels in AOIs

3 结果评价与分析

3.1 检测精度评价

由于研究区所采用Apollo图像的空间分辨率为6.909 m，并且在这样的空间分辨率下没有已知对应的撞击坑数据库，因此为了评价自动提取月表撞击坑的准确性，需在研究区内人工提取撞击坑(如图5中的绿色椭圆)。准确性评价采用1993年Shufelt［23］提出的评价方法，即

式中:Pd为检测比率;Fb为分歧比率;Pq为质量比率;Pt为该区域内人工提取与自动提取都检测出的撞击坑数量(正确提取的撞击坑数量);Pf表示因拟合算法的缺陷造成的人工未提取出，而自动提取结果中却存在的撞击坑数量(误提取的撞击坑数量);Nf表示人工提取出，但因图像质量等造成拟合算法未能提取的撞击坑数量(遗漏提取的撞击坑数量)。具体评价结果见表1。

表1 撞击坑自动检测结果Tab.1 Results of automatic detection of craters

从表1可看出:由于采用了分层的椭圆提取方法，所以大、小撞击坑都可提取出来;研究区内大撞击坑(长半轴D＞25个像元)数量上比较少，而小撞击坑(8≤D≤25)数量上比较多;大撞击坑由低太阳高度角造成的黑白边缘特征更加明显与清晰，检测比率略高于小撞击坑，总体检测质量略好;对于全部的撞击坑，检测质量可以达到76%以上。

3.2 撞击坑大小累积分布特征

撞击坑的数量特征往往通过累计大小频率分布(SFD)来体现［24］。SFD描述了撞击坑的长半轴(D)与大于该直径的撞击坑的个数(N)的关系，满足N=kD－b的幂指数关系(其中－b是对2个坐标轴求对数之后的斜率)。本文利用所提取的椭圆(Pt)的长半轴，统计得出研究区内撞击坑的SFD曲线(图6)。

图6 撞击坑的累积大小频率分布Fig.6 SFD of the craters

在提取的撞击坑中，长半轴D=160 m(图6中的红色标记点)是一个分界点。当D＜160 m时，N与D取对数后，有明显的线性关系，即该曲线满足N=kD－3.8( －3.8 为取对数后的斜率)。这与 Shoemaker等学者［25－28］通过统计得出的小撞击坑(D ＜175 m)的幂指数在［3.5，4.0］范围内的结论相一致，说明这种关系具有普遍性，同时也从侧面说明了本文提取方法的准确性与可行性。

4 结论

1)本文利用Apollo 15图像数据，根据月表撞击坑的真实形状，提出一种撞击坑提取的新算法。通过对比度增强、双边滤波消除噪声等图像预处理;利用改进型Canny边缘检测和伪边缘剔除等进行边缘提取;利用最小二乘椭圆拟合与基于约束参数的Hough多椭圆变换等多个层次的拟合方法，在不受环形山的大小、歪曲等因素的影响下，稳定、高精度将撞击坑自动提取出来。

2)将撞击坑的自动提取结果与人工提取结果比较，总体检测精度在76%以上，直径比较大的撞击坑的检测精度则更高，可达83%。通过统计研究区撞击坑的累积大小频率分布(SFD)结果发现，在撞击坑长半轴D＜160 m时，幂指数为－3.8。

3)由于受到获取图像时的光照条件以及月表形貌复杂特征的影响，存在部分伪边缘检测现象。同时也应该注意到，在撞击坑检测精度评价中没有统计自动提取的椭圆长短轴与人工提取的椭圆长短轴的差别，评价方法还有待改进。

4)由于本文研究区处于受3个大型撞击坑(哥白尼(Copernicus)、阿里斯塔克(Aristarchus)、开普勒(Kepler)撞击坑)强烈影响的区域，区内尚有很多次生撞击坑。下一步的工作是将该方法应用于嫦娥二号数据，并尝试结合嫦娥一号数据，找到区分次生撞击坑以及原生撞击坑的因子，并分析次生撞击坑的父坑位置等。

［1］Shoemaker E M，Hackman R J，Eggleton R E.Interplanetary Correlation of Geologic Time［J］.Adv Astronaut Sci，1962，8:70 －89.

［2］Nagumo K，Nakamura A M.Reconsideration of Crater Size－frequency Distribution on the Moon:Effect of Projectile Population and Secondary Craters［J］.Adv Space Res，2001，28(8):1181 －1186.

［3］Florensky C P，Basilevsky A T，Grebennik N N.The Relationship Between Lunar Crater Morphology and Crater Size［J］.Earth，Moon，and Planets，1976，16(1):59 －70.

［4］Moore H J，Boyce J M，Hahn D A.Small Impact Craters in the Lunar Regolith:Their Morphologies，Relative Ages，and Rates of Formation［J］.The Moon and the Planets，1980，23(2):231 －252.

［5］Vinogradova T，Burl M，Mjolsness E.Training of a Crater Detection Algorithm for Mars Crater Imagery［C］//Proceeding of the IEEE Aerospace Conference，2002，7:3201 －3211.

［6］Honda R，Iijima Y，Konishi O.Mining of Topographic Feature from Heterogeneous Imagery and Its Application to Lunar Craters［C］//Proceeding of the Progress of Discovery Science London:Springer－ Verlag，2002:395 －407.

［7］Cheng Y，Johnson A E，Matthies L H，et al.Optical Landmark Detection for Spacecraft Navigation［C］//Proceedings of the 13th AAS/AIAA Space Flight Mechanics Meeting.Ponce，Puerto Rico，2003:203－235.

［8］Kim J R，Muller J P.Impact Crater Detection on Optical Image and DEM.ISPRS WG IV/9:Extraterrestrial Mapping Workshop“Advances in Planetary Mapping 2003”［EB/OL］.2003.http://astrogeology.usgs.gov/Projects/ISPRS/MEETINGS/Houston2003/index_houston.html.

［9］Wan C，Cheng W M，Zhou Z P，et al.Automatic Extraction of Lunar Impact Craters from Chang’E － 1 Satellite Photographs［J］.Science China Physics，Mechanics and Astronomy，2012，55(1):162－169.

［10］Tomasi C，Manduchi R.Bilateral Filtering for Gray and Color Ima-ges［C］//Proceedings of the 6th IEEE International Conference on Computer Vision.Bombay:IEEE，1998:839 －846.

［11］Canny J.A Computational Approach to Edge Detection［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1986，8(6):679－698.

［12］余洪山，王耀南.一种改进型canny边缘检测算法［J］.计算机工程与应用，2004，20(1):27 －29.Yu H S，Wang Y N.An Improved Canny Edge Detection Algorithm［J］.Computer Engineering and Applications，2004，20(1):27 －29(in Chinese with English Abstract).

［13］冯军华，崔祜涛，崔平远，等.行星表面陨石坑检测与匹配方法［J］.航空学报，2010，31(9):1858 －1863.Feng J H，Cui H T，Cui P Y，et al.Autonomous Crater Detection and Matching on Planetary Surface［J］.Acta Aeronautica Et Astronautica Sinica，2010，31(9):1858 －1863(in Chinese with English Abstract).

［14］Wei Y.Circle Detection Using Improved Dynamic Generalized Hough Transform［C］//Proceedings of IEEE International Geoscience and Remote Sensing Symposium.Seattle，WA:IEEE，1988，2:1190－1192.

［15］Zhang S C，Liu Z Q.A New Algorithm for Real－time Ellipse Detection［C］//Proceedings of IEEE International Conference on Machine Learning and Cybernetics.XI’an，China:IEEE，2003，1:602 －607.

［16］Ji Q，Haralick R M.A Statistically Efficient Method for Ellipse Detection［C］//Proceedings of 1999 International Conference on Image Processing.Kobe，Japan:IEEE，1999:730 －734.

［17］Kawaguchi T，Nagata R I.Ellipse Detection Using a Genetic Algorithm［C］//Proceedings of the 14th International Conference on Pattern Recognition Brisbane，Qld:IEEE，1998，1:141 －145.

［18］Xie Y H，Ji Q.A New Efficient Ellipse Detection Method［C］//Proceedings of the 16th International Conference on Pattern Recognition Quebec，Canada:IEEE，2002:957 －960.

［19］Elmowafy O M，Fairhurst M C.Improving Ellipse Detection Using a Fast Graphical Method［J］.Electronics Letters，1999，35(2):135－137.

［20］杨忠根，马 彦.使用广义正交概念的K－RANSAC椭圆提取［J］.自动化学报，2002，28(4):520 －526.Yang Z G，Ma Y.A New Method for Ellipse Detection Using K －RANSAC Based on Generalized Orthogonality Principle［J］.Acta Automatica Sinica，2002，28(4):520 － 526(in Chinese with English Abstract).

［21］Hough P V C.Method and Means for Recognizing Complex Patterns:US，3969654［P］.1962.

［22］牛晓霞，胡正平，杨 苏.局部PCA参数约束的Hough多椭圆分层检测算法［J］.计算机应用，2009，29(5):1365 －1368.Niu X X，Hu Z P，Yang S.Hierarchical Ellipse Detection Algorithm Based on Local PCA Hough Transform with Parameter Restraint［J］.Journal of Computer Applications，2009，29(5):1365 － 1368(in Chinese with English Abstract).

［23］Shufelt J，McKewon D M.Fusion of Monocular Cues to Detect Man － made Structures in Serial Imagery［J］.Computer Vision and Image Understanding，1993，57(3):307 －330.

［24］Crater Analysis Techniques Working Group.Standard Techniques for Presentation and Analysis of Crater Size － frequency Data［J］.Icarus，1979，37(2):467 －474.

［25］Shoemaker E M.Preliminary Analysis of the Fine Structure of the Lunar Surface in Mare Cognitum［M］//The Nature of the Lunar Surface.Baltimore:Johns Hopkins University Press，1965:23 －77.

［26］Ivanov B A.Notes About Secondary Crater SFD［C］//Workshop on Surface Ages and Histories:Issues in Planetary Chronology.LPI Contribution No.1320.Houston:Lunar and Planetary Institute，2006:28－29.

［27］Hartmann W K，Gaskell R W.Planetary Cratering 2:Studies of Saturation Equilibrium［J］.Meteorit Planet Sci，1997，32:109 －121.

［28］Hartmann W K.Martian Cratering 8:Isochron Refinement and the Chronology of Mars［J］.Icarus，2005，174(2):294 －320.

Automatic Detection of Lunar Elliptical Craters from Apollo Image

LI Chao1，2，WANG Xin － yuan1，LUO Lei1，2，JI Wei1
(1.Center for Earth Observation and Digital Earth，Chinese Academy of Sciences，Beijing 100094，China;2.Graduate University of Chinese Academy of Sciences，Beijing 100049，China)

Impact crater is the major geomorphologic unit of the lunar surface.Its shape and distribution characteristics can provide a large amount of information about the impact history and the lunar evolution process.In this paper，based on the real shape of craters，the authors proposed an automatic extraction method of elliptical craters.The method consists of three steps:image preprocessing，edge detection and ellipse fitting.The last step uses principal component analysis(PCA)of local interesting parameters and Hough transform to derive multiple ellipses.The total accuracy based on a comparison with manually－derived result is over 76%，with the result of larger craters even in excess of 83%.Statistics of the size－frequency distributions of the craters show that，when the length of the long half axis(D)is smaller than 160 meters，the crater number N and D－3.8has a significant linear relationship.

lunar crater;automatic detection;Hough multi－ellipse transform;edge

TP 751.1

A

1001－070X(2012)04－0071－05

2012－03－01;

2012－03－25

国家高技术研究发展计划(863计划)项目(编号:2010AA122202－02)和国家自然科学基金项目(编号:60972141)共同资助。

10.6046/gtzyyg.2012.04.12

李 超(1986－)，男，硕士研究生，主要从事遥感图像处理、InSAR技术研究。E－mail:dearth0531@126.com。

王心源(1964－)，男，博士，研究员，主要从事遥感考古、自然－文化遗产信息技术应用研究。E－mail:xywang@ceode.ac.cn。

(责任编辑: 邢 宇)