李 超，王心源，駱 磊，吉 瑋

(1.中國科學院對地觀測與數(shù)字地球科學中心數(shù)字地球重點實驗室，北京 100094;2.中國科學院研究生院，北京 100049)

基于Apollo圖像的月表撞擊坑自動提取

李 超1，2，王心源1，駱 磊1，2，吉 瑋1

(1.中國科學院對地觀測與數(shù)字地球科學中心數(shù)字地球重點實驗室，北京 100094;2.中國科學院研究生院，北京 100049)

撞擊坑是月表的主要形貌單元。研究月表撞擊坑的形狀與分布特征對描述月表撞擊過程與演變規(guī)律有很大幫助。根據(jù)月表撞擊坑的真實形狀，提出了一套橢圓形撞擊坑自動提取方法，其主要包括3個步驟:圖像預處理、邊緣檢測與約束參數(shù)的Hough多橢圓分層檢測。與人工提取結(jié)果比較，月表撞擊坑自動提取的總體準確率在76%以上(其中比較大的撞擊坑的檢測率可以達到83%)。統(tǒng)計檢測結(jié)果表明，當撞擊坑的長半軸D＜160 m時，其累積大小頻率分布曲線撞擊坑個數(shù)N與D－3.8有明顯的線性關(guān)系。

月表撞擊坑;自動提取;Hough多橢圓變換;邊緣

0 引言

撞擊坑(craters)對于研究類地行星(衛(wèi)星)的形貌、相對地質(zhì)年齡、巖石構(gòu)造以及行星起源有著重要意義。Shoemaker等［1］認為，在行星表面上的撞擊坑數(shù)量可以反映不同年代的撞擊現(xiàn)象。這個原理是基于撞擊坑在空間與時間隨機形成的假設，常用以評價星體表面形成的絕對或相對年齡。這種理論是通過統(tǒng)計撞擊坑的累積大小頻率分布(sizefrequency distribution，SFD)來體現(xiàn)的。以前的學者［2－4］大部分是通過人工方法提取撞擊坑進行分析的，但由于這種方法具有耗時長、缺乏客觀性、對操作者專業(yè)知識要求高等弱點，使自動提取撞擊坑的要求逐漸強烈起來。近年來，許多學者從不同方面提出撞擊坑自動檢測方法［5－9］，形成了各自的特征識別系統(tǒng)，包括模式匹配、紋理分析、基因算法與面向?qū)ο蟮榷喾N提取方法。這些方法都是將撞擊坑抽象為圓形。但是，由于特定的撞擊角度以及傳感器拍攝視角等原因，撞擊坑的形狀更符合橢圓形。本文利用Apollo 15圖像數(shù)據(jù)，根據(jù)月表撞擊坑的真實形狀，選取受到3個大型撞擊坑(哥白尼(Copernicus)、阿里斯塔克(Aristarchus)和開普勒(Kepler)撞擊坑)強烈影響、內(nèi)部有很多次生撞擊坑(secondary craters)的區(qū)域作為研究區(qū)，旨在提出一種新的橢圓形檢測方法，為提取與分析月表次生撞擊坑提供方法支持。

1 研究區(qū)與數(shù)據(jù)源

如圖1所示，研究區(qū)位于月海內(nèi)的歐拉(Euler)撞擊坑與丟番圖(Diophantus)撞擊坑之間。研究所用數(shù)據(jù)為 Apollo 15 圖像數(shù)據(jù)(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 撞擊坑檢測方法

研究區(qū)內(nèi)撞擊坑特征十分復雜，為了便于計算機處理，在研究區(qū)內(nèi)選擇6個小范圍的、具有便于識別撞擊坑特征的感興趣區(qū)(area of interest，AOI)(圖1中紅色矩形標識的區(qū)域)。在獲取圖像中的AOI之后，通過圖2所示的3個步驟自動提取撞擊坑。

圖2 撞擊坑自動檢測流程圖Fig.2 Flow chart of automatic detection of craters

2.1 數(shù)據(jù)預處理

根據(jù)圖像局部灰度值的統(tǒng)計數(shù)據(jù)對選取的AOI進行增強處理，增大圖像的對比度，有選擇地突出圖像中感興趣的撞擊坑輪廓特征，同時抑制(掩蓋)圖像中其他信息，使圖像與計算機或人的視覺響應特性相匹配。

為消除圖像噪聲對月表撞擊坑檢測的影響，采用非線性、非迭代的雙邊濾波(bilateral filtering)方法［10］對圖像進行處理。這種雙邊濾波在處理相鄰像素時，不僅考慮到距離上的鄰近關(guān)系，同時還考慮到灰度上的相似性，通過對這二者的非線性組合，可自適應地進行濾波。這樣圖像在濾除噪聲的同時，也能很好地保留圖像的邊緣信息。Apollo原始圖像與預處理結(jié)果如圖3所示。

圖3 原始圖像與預處理結(jié)果Fig.3 Original image and preprocessing result

2.2 邊緣提取

分析Apollo圖像特點及綜合比較各種邊緣檢測方法后，最終采用Canny算法［11］提取月表撞擊坑邊緣。鑒于傳統(tǒng)Canny算子的計算量比較大、實時性較差，并且檢測邊緣的精度仍沒有達到單像素級，同時受到多種因素的影響造成部分偽邊緣，為此本文采用改進型Canny邊緣檢測算法［12］提取月表撞擊坑邊緣并進行細化，結(jié)果如圖4(a)所示。由于獲取的Apollo圖像數(shù)據(jù)的太陽高度角為7°，會使撞擊坑高出的坑壁出現(xiàn)陰影。與Apollo原始圖像(圖3(a))比較可以看出，提取的月表撞擊坑邊緣(圖4(a))中既包含了撞擊坑的真實邊緣信息，又包括了因光照條件和地形因素造成的撞擊坑內(nèi)部光亮區(qū)與陰影區(qū)的偽邊緣信息以及其他地形特征的邊緣信息。通過檢測撞擊坑圖像的梯度方向［13］，發(fā)現(xiàn)在撞擊坑內(nèi)部的偽邊緣與撞擊坑的真實邊緣的梯度方向正好相反，真實撞擊坑的邊緣梯度方向與光源方向之間的夾角小于90°，而偽邊緣梯度方向與光源方向之間的夾角大于90°，即真實邊緣滿足

式中:g為邊緣梯度矢量;n為光源方向矢量;|g|為梯度矢量的模長;|n|為光源方向矢量的模長。

通過上述方法選擇邊緣梯度方向與光源方向夾角小于90°的邊緣，就可以有效剔除偽邊緣(圖4(b))。最后，按照邊緣的方向進行一定鄰域內(nèi)的連接，同時考慮噪聲以及下文論述的Hough變換的計算量，通過設置一定的閾值將較小的邊緣從中剔除。

圖4 邊緣檢測Fig.4 Edge detection

2.3 橢圓擬合

現(xiàn)有的橢圓擬合算法大致分為5類:①基于經(jīng)典的Hough變換及其改進算法［14－15］;②有效統(tǒng)計算法(statistically efficient method)［16］; ③遺傳算法(genetic algorithm)［17］; ④ 基于橢圓 幾何 特性 算法［18－19］;⑤基于廣義正交概念的 K隨機采樣(K－RANSAC)及其改進算法［20］。橢圓的自由參數(shù)包括橢圓的中心坐標(x，y)，長、短半軸的長度(D，B)、長軸的偏轉(zhuǎn)角度(θ)。在實際橢圓檢測應用中，對于多個橢圓的檢測往往存在計算耗時長、內(nèi)存需求大等問題，因此如何降低橢圓的空間維度以便加快檢測速度已成為橢圓檢測的主要問題。

Hough變換是一種在圖像中定位形狀的技術(shù)［21］。它是基于“證據(jù)收集”方法，通過模板匹配過程的重新描述來實現(xiàn)的(其中，證據(jù)是在累加器數(shù)組中的投票)。Hough變換的基本思想是利用點—線的對偶性，對圖像進行坐標變換，使之在另一個坐標系(累加器空間)中的特定位置上出現(xiàn)峰值，檢測曲線的問題就簡化為尋找峰值位置的問題。這種變換的優(yōu)點是受噪聲與曲線間斷的影響較小，但最大的缺點是計算量非常大，且需要比較大的內(nèi)存。為此，改進的Hough變換被先后提出，比如概率Hough變換(PHT)、隨機Hough變換(RHT)等。但是，這類方法主要針對單個或者簡單的橢圓，直接應用到月表撞擊坑的提取，效果還不理想。本文結(jié)合最小二乘與Hough變換，通過選取AOI，采取基于局部主成分分析(PCA)感興趣參數(shù)約束Hough多橢圓分層檢測的思路［22］，方法步驟如下:①設定長度閾值T1與T2(T1＞T2);②在邊緣圖像中標記并提取長度大于T1的連續(xù)曲線段;③對各線段做PCA方向分析，確定是否屬于有效曲線段;④對所有感興趣曲線段按照標記順序依次利用最小二乘橢圓擬合方法得到感興趣橢圓粗略參數(shù);⑤根據(jù)擬合結(jié)果進而模糊約束Hough變換參數(shù)搜索范圍，得到精確橢圓參數(shù)，并輸出這5個自由參數(shù);⑥利用檢測結(jié)果更新圖像空間，刪除已經(jīng)檢測到的橢圓;并按照一定的步長迭代進行，直到連續(xù)曲線段長度小于T2，所有橢圓檢測完畢。圖5為按照上述步驟得到的6個AOI最終的橢圓檢測結(jié)果。

圖5 AOIs區(qū)域的橢圓(＞7個像元)擬合結(jié)果Fig.5 Results of the ellipse fitting larger than seven pixels in AOIs

3 結(jié)果評價與分析

3.1 檢測精度評價

由于研究區(qū)所采用Apollo圖像的空間分辨率為6.909 m，并且在這樣的空間分辨率下沒有已知對應的撞擊坑數(shù)據(jù)庫，因此為了評價自動提取月表撞擊坑的準確性，需在研究區(qū)內(nèi)人工提取撞擊坑(如圖5中的綠色橢圓)。準確性評價采用1993年Shufelt［23］提出的評價方法，即

式中:Pd為檢測比率;Fb為分歧比率;Pq為質(zhì)量比率;Pt為該區(qū)域內(nèi)人工提取與自動提取都檢測出的撞擊坑數(shù)量(正確提取的撞擊坑數(shù)量);Pf表示因擬合算法的缺陷造成的人工未提取出，而自動提取結(jié)果中卻存在的撞擊坑數(shù)量(誤提取的撞擊坑數(shù)量);Nf表示人工提取出，但因圖像質(zhì)量等造成擬合算法未能提取的撞擊坑數(shù)量(遺漏提取的撞擊坑數(shù)量)。具體評價結(jié)果見表1。

表1 撞擊坑自動檢測結(jié)果Tab.1 Results of automatic detection of craters

從表1可看出:由于采用了分層的橢圓提取方法，所以大、小撞擊坑都可提取出來;研究區(qū)內(nèi)大撞擊坑(長半軸D＞25個像元)數(shù)量上比較少，而小撞擊坑(8≤D≤25)數(shù)量上比較多;大撞擊坑由低太陽高度角造成的黑白邊緣特征更加明顯與清晰，檢測比率略高于小撞擊坑，總體檢測質(zhì)量略好;對于全部的撞擊坑，檢測質(zhì)量可以達到76%以上。

3.2 撞擊坑大小累積分布特征

撞擊坑的數(shù)量特征往往通過累計大小頻率分布(SFD)來體現(xiàn)［24］。SFD描述了撞擊坑的長半軸(D)與大于該直徑的撞擊坑的個數(shù)(N)的關(guān)系，滿足N=kD－b的冪指數(shù)關(guān)系(其中－b是對2個坐標軸求對數(shù)之后的斜率)。本文利用所提取的橢圓(Pt)的長半軸，統(tǒng)計得出研究區(qū)內(nèi)撞擊坑的SFD曲線(圖6)。

圖6 撞擊坑的累積大小頻率分布Fig.6 SFD of the craters

在提取的撞擊坑中，長半軸D=160 m(圖6中的紅色標記點)是一個分界點。當D＜160 m時，N與D取對數(shù)后，有明顯的線性關(guān)系，即該曲線滿足N=kD－3.8( －3.8 為取對數(shù)后的斜率)。這與 Shoemaker等學者［25－28］通過統(tǒng)計得出的小撞擊坑(D ＜175 m)的冪指數(shù)在［3.5，4.0］范圍內(nèi)的結(jié)論相一致，說明這種關(guān)系具有普遍性，同時也從側(cè)面說明了本文提取方法的準確性與可行性。

4 結(jié)論

1)本文利用Apollo 15圖像數(shù)據(jù)，根據(jù)月表撞擊坑的真實形狀，提出一種撞擊坑提取的新算法。通過對比度增強、雙邊濾波消除噪聲等圖像預處理;利用改進型Canny邊緣檢測和偽邊緣剔除等進行邊緣提取;利用最小二乘橢圓擬合與基于約束參數(shù)的Hough多橢圓變換等多個層次的擬合方法，在不受環(huán)形山的大小、歪曲等因素的影響下，穩(wěn)定、高精度將撞擊坑自動提取出來。

2)將撞擊坑的自動提取結(jié)果與人工提取結(jié)果比較，總體檢測精度在76%以上，直徑比較大的撞擊坑的檢測精度則更高，可達83%。通過統(tǒng)計研究區(qū)撞擊坑的累積大小頻率分布(SFD)結(jié)果發(fā)現(xiàn)，在撞擊坑長半軸D＜160 m時，冪指數(shù)為－3.8。

3)由于受到獲取圖像時的光照條件以及月表形貌復雜特征的影響，存在部分偽邊緣檢測現(xiàn)象。同時也應該注意到，在撞擊坑檢測精度評價中沒有統(tǒng)計自動提取的橢圓長短軸與人工提取的橢圓長短軸的差別，評價方法還有待改進。

4)由于本文研究區(qū)處于受3個大型撞擊坑(哥白尼(Copernicus)、阿里斯塔克(Aristarchus)、開普勒(Kepler)撞擊坑)強烈影響的區(qū)域，區(qū)內(nèi)尚有很多次生撞擊坑。下一步的工作是將該方法應用于嫦娥二號數(shù)據(jù)，并嘗試結(jié)合嫦娥一號數(shù)據(jù)，找到區(qū)分次生撞擊坑以及原生撞擊坑的因子，并分析次生撞擊坑的父坑位置等。

［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，F(xiàn)airhurst M C.Improving Ellipse Detection Using a Fast Graphical Method［J］.Electronics Letters，1999，35(2):135－137.

［20］楊忠根，馬 彥.使用廣義正交概念的K－RANSAC橢圓提?。跩］.自動化學報，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參數(shù)約束的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

國家高技術(shù)研究發(fā)展計劃(863計劃)項目(編號:2010AA122202－02)和國家自然科學基金項目(編號:60972141)共同資助。

10.6046/gtzyyg.2012.04.12

李 超(1986－)，男，碩士研究生，主要從事遙感圖像處理、InSAR技術(shù)研究。E－mail:dearth0531@126.com。

王心源(1964－)，男，博士，研究員，主要從事遙感考古、自然－文化遺產(chǎn)信息技術(shù)應用研究。E－mail:xywang@ceode.ac.cn。

(責任編輯: 邢 宇)