劉宇軒，李春來(lái)，劉建軍

(1. 中國(guó)科學(xué)院國(guó)家天文臺(tái)，北京 100101；2. 中國(guó)科學(xué)院大學(xué)，北京 100049，3. 中國(guó)科學(xué)院月球與深空探測(cè)重點(diǎn)實(shí)驗(yàn)室，北京 100101)

行星表面通常分布著大量的撞擊坑構(gòu)造，它們是行星表面最顯著的地貌特征，也是研究行星內(nèi)部物質(zhì)的窗口[1]。由于月球上沒(méi)有大氣，撞擊坑保存比較完好，因此研究月球形貌、相對(duì)地質(zhì)年齡、巖石構(gòu)造以及月球的起源有重要意義。撞擊坑識(shí)別的研究，從早年使用望遠(yuǎn)鏡觀察描繪撞擊坑的形貌，到利用遙感圖像數(shù)據(jù)目視識(shí)別撞擊坑，再到近年來(lái)通過(guò)計(jì)算機(jī)自動(dòng)識(shí)別撞擊坑。隨著月球影像與地形數(shù)據(jù)的大量獲取，撞擊坑的自動(dòng)識(shí)別提取成為一個(gè)新的研究熱點(diǎn)。近年來(lái)，國(guó)內(nèi)外眾多學(xué)者提出大量的研究方法，并建立多種尺度的撞擊坑數(shù)據(jù)庫(kù)[2-6]。通過(guò)對(duì)現(xiàn)有方法的分析總結(jié)，撞擊坑識(shí)別的方法大致可以歸納為4類：人工目視識(shí)別、基于形態(tài)擬合、基于機(jī)器學(xué)習(xí)以及基于地學(xué)信息分析的方法。人工目視識(shí)別法：利用影像和地形數(shù)據(jù)，建立識(shí)別撞擊坑的解譯標(biāo)識(shí)，在計(jì)算機(jī)環(huán)境下，采用人工交互的方式對(duì)撞擊坑逐個(gè)識(shí)別，這種方法是目前最準(zhǔn)確的撞擊坑識(shí)別方法，該方法的優(yōu)點(diǎn)是準(zhǔn)確性高，但費(fèi)時(shí)費(fèi)力，不適合于全月球和小尺度撞擊坑識(shí)別?；谛螒B(tài)擬合法：由于撞擊坑在圖像上有類圓形的幾何特點(diǎn)，在圖像分割的基礎(chǔ)上以霍夫變換和二次曲線為代表的形狀擬合算法用于模擬提取完整的撞擊坑邊界，其中文[7]利用組合霍夫變換算法，文[8]采用模糊霍夫變換算法對(duì)撞擊坑進(jìn)行識(shí)別，文[9]提出一種基于弦中點(diǎn)霍夫變換算法，文[10]采用最小二乘法擬合邊緣橢圓的方法實(shí)現(xiàn)對(duì)撞擊坑的提取?；跈C(jī)器學(xué)習(xí)的方法：文[11]采用BP神經(jīng)網(wǎng)絡(luò)加訓(xùn)練樣本的算法，文[12]利用決策樹(shù)算法，文[13]利用面向?qū)ο蠓椒▽?shí)現(xiàn)對(duì)撞擊坑的識(shí)別。基于地學(xué)信息分析方法：文[14]利用坡度作為識(shí)別撞擊坑的要素，文[15-17]利用500 m分辨率的嫦娥一號(hào)數(shù)字高程模型數(shù)據(jù)和形態(tài)擬合的方法，進(jìn)行撞擊坑的識(shí)別提取和研究，文[18]利用坡度數(shù)據(jù)、紋理信息和剖面曲率作為識(shí)別依據(jù)，采用水淹法對(duì)火星撞擊坑進(jìn)行分析。

本文使用目前世界上覆蓋最全、精度最高的嫦娥二號(hào)數(shù)字地形數(shù)據(jù)，提出基于等高線地形分析的方法，針對(duì)小型撞擊坑進(jìn)行識(shí)別提取研究，填補(bǔ)撞擊坑數(shù)據(jù)庫(kù)小型撞擊坑數(shù)量不足和準(zhǔn)確率不高的缺陷。

1 數(shù) 據(jù)

2008年2月，國(guó)務(wù)院批準(zhǔn)探月工程二期立項(xiàng)，二期工程主要目標(biāo)是實(shí)現(xiàn)月面軟著陸，開(kāi)展月面就位探測(cè)與自動(dòng)巡視勘察。嫦娥二號(hào)作為二期工程的技術(shù)先導(dǎo)星，首要的任務(wù)是對(duì)全月球和嫦娥三號(hào)的月面著陸區(qū)進(jìn)行詳查，精細(xì)測(cè)繪全月球的地形地貌。嫦娥二號(hào)衛(wèi)星于北京時(shí)間2010年10月1日18時(shí)59分57秒在西昌衛(wèi)星發(fā)射中心成功發(fā)射。10月24日16時(shí)49分，嫦娥二號(hào)衛(wèi)星上的兩線陣CCD立體相機(jī)首次開(kāi)機(jī)工作，并成功獲取了前視與后視兩個(gè)視角的影像數(shù)據(jù)。截至2011年5月20日，嫦娥二號(hào)衛(wèi)星CCD立體相機(jī)完成全部任務(wù)，共獲取607軌圖像數(shù)據(jù)。全月球影像圖制作過(guò)程中選取的100 km軌道高度成像的384軌數(shù)據(jù)(其中南北緯70°范圍344軌，極區(qū)40軌)，用于全月球數(shù)字高程模型數(shù)據(jù)產(chǎn)品生產(chǎn)。數(shù)據(jù)處理流程包括：數(shù)據(jù)準(zhǔn)備、測(cè)區(qū)數(shù)據(jù)組織、平差點(diǎn)準(zhǔn)備、全球平差、測(cè)區(qū)初級(jí)產(chǎn)品生產(chǎn)、標(biāo)準(zhǔn)產(chǎn)品生產(chǎn)、全球數(shù)據(jù)鑲嵌等。結(jié)合嫦娥二號(hào)衛(wèi)星軌道與姿態(tài)測(cè)量數(shù)據(jù)、5個(gè)月面絕對(duì)控制點(diǎn)進(jìn)行全球平差，經(jīng)過(guò)4年的數(shù)據(jù)處理，成功制作了7 m、20 m、50 m分辨率的數(shù)字高程模型數(shù)據(jù)。嫦娥二號(hào)數(shù)字高程模型數(shù)據(jù)在空間分辨率、數(shù)據(jù)一致性和完整性鑲嵌精度等方面優(yōu)于國(guó)際同類全月球數(shù)字產(chǎn)品。本文選擇7 m分辨率的數(shù)字高程模型數(shù)據(jù)進(jìn)行小型撞擊坑的識(shí)別提取。

2 方 法

2.1 撞擊坑模型分析

根據(jù)基本形態(tài)，撞擊坑可分為簡(jiǎn)單撞擊坑(碗型小型撞擊坑)、復(fù)雜撞擊坑(距中央峰的較大型撞擊坑)、多環(huán)撞擊坑(撞擊盆地)3類[19]。如圖1，一個(gè)理想狀態(tài)下的復(fù)雜撞擊坑模型主要由坑唇邊緣(rh)、坑壁(w)、坑底(f)、撞擊培巖(m)和中央峰或中央環(huán)(p)5種形態(tài)要素組成。撞擊濺射堆積物由坑唇邊緣(rh)、丘狀沉積物(rr)和輻射線沉積物(rc)構(gòu)成。撞擊溶巖(m)覆蓋在其他形態(tài)要素之上。簡(jiǎn)單碗型撞擊坑(c)一般分布在丘狀沉積物(rr)之上。撞擊坑識(shí)別是一種特征識(shí)別，對(duì)影像或地形數(shù)據(jù)進(jìn)行分析處理，如圖2，識(shí)別撞擊坑的坑唇邊緣(rh)，并提取撞擊坑的直徑(D)、中心位置和高程等信息。

圖1 理想的的復(fù)雜撞擊坑模型(Don wilhelms, 1987)
Fig.1 Ideal complex impact crater model

等高線是地形圖上高程相等的相鄰各點(diǎn)連成的閉合曲線，是最為常用的表達(dá)地形形態(tài)特征和起伏情況的形式，如圖3。等高線從數(shù)字高程模型中提取出來(lái)，科學(xué)地反映月表的高程 、坡度、坡向、高地走勢(shì)等基本月表地形地貌。嫦娥二號(hào)數(shù)字高程模型數(shù)據(jù)的空間分辨率很高，月表高程在二維地理空間上得到了連續(xù)表達(dá)，利用該數(shù)字高程模型數(shù)據(jù)較為方便、精細(xì)地提取月球撞擊坑地形等高線，表達(dá)月球撞擊坑形態(tài)，并對(duì)其進(jìn)行分析和處理。

從等高線的形態(tài)分析，如圖4撞擊坑具有3個(gè)典型的空間特點(diǎn)：(1)撞擊坑坑唇邊緣等高線通常為近似圓形或橢圓形；(2)撞擊坑坑唇邊緣等高線內(nèi)側(cè)通常存在多條嵌套的相似等高線，多為同心圓或橢圓，且間距較小，外側(cè)等高線形狀不定，且間距較遠(yuǎn)，與坑唇邊緣有明顯區(qū)別；(3)撞擊坑坑底高程小，坑邊緣高程大，其剖面線可以用下凹的拋物線形態(tài)擬合。因此，通過(guò)撞擊坑地形等高線的特征以及與其他等高線之間的相互關(guān)系，可以分析并識(shí)別月表的撞擊坑。因此，只要選取合適的等高距，撞擊坑坑唇邊緣的閉合近圓形等高線就可以逼近真實(shí)邊緣，可以準(zhǔn)確識(shí)別并提取撞擊坑的地理位置、大小和深度等信息。

圖2 簡(jiǎn)單撞擊坑和復(fù)雜撞擊坑剖面圖
Fig.2 Impact crater structure

圖3 等高線示意圖
Fig.3 Contour diagram

圖4 撞擊坑等高線示意圖
Fig.4 Crater contour diagram

2.2 指標(biāo)體系和參數(shù)的確定

撞擊坑等高線往往表現(xiàn)為相互嵌套的類圓環(huán)狀。因此，初步判別撞擊坑特征的思路是根據(jù)等高線的圓度、矩形度、體態(tài)比等指標(biāo)參數(shù)保留滿足條件的等高線，再分析等高線之間的關(guān)系，聚合嵌套的等高線并標(biāo)示外層等高線；然后通過(guò)剖面線的凹凸剔除穹隆等正向地形區(qū)域，保留負(fù)向地形的撞擊坑區(qū)域；最后通過(guò)霍夫變換擬合撞擊坑的邊緣，得到撞擊坑的中心位置和大小。因此，指標(biāo)體系和參數(shù)包括圓度、矩形度、體態(tài)比、剖面線凹凸度。

(1)圓形度是圖像處理中重要的概念，用于特征提取與描述，其值反映被測(cè)量邊界的復(fù)雜程度，公式為

e=4πA0/P2,

(1)

其中，A0為面積；P為周長(zhǎng)。這個(gè)特征對(duì)圓形面域取最大值1，坑唇邊緣越接近1的等高線越符合條件。

(2)矩形度是常用的目標(biāo)物體矩形擬合因子，其公式為

e1=A0/AR，

(2)

其中，A0為面積；AR為最小外接矩形的面積；e1表示某物體對(duì)其最小外接矩形的充滿程度。對(duì)于矩形物體，e1取最大值1；對(duì)于圓形物體，e1取0.79；對(duì)于其他細(xì)長(zhǎng)物體，e1取值變小。矩形擬合因子的值限定在0～1之間。

(3)體態(tài)比定義為

e2=W/L，

(3)

其中，W為最小外接矩形的寬；L為最小外接矩形的長(zhǎng)。圓形面域的體態(tài)比取值為1，這個(gè)參數(shù)可以把細(xì)長(zhǎng)不規(guī)則的等高線剔除。

(4)剖面線凹凸性：因?yàn)樵卤磉€存在其他隆起的地形，如穹隆也會(huì)產(chǎn)生近圓形的閉合等高線。潛在撞擊坑還存在，對(duì)此情況的處理方法是判斷潛在撞擊坑所在的數(shù)字高程模型剖面線是否為下凹形態(tài)，如果剖面線為下凹形態(tài)，則是撞擊坑；若剖面線為上凸形態(tài)，則不是撞擊坑。根據(jù)內(nèi)等高線與等高線的數(shù)值大小關(guān)系確定是否為真實(shí)撞擊坑。

2.3 撞擊坑邊界的擬合

對(duì)滿足條件的等高線，利用霍夫變化的方法擬合其邊界，最終得到中心點(diǎn)和半徑大小。

霍夫變換算法是圖像處理中識(shí)別幾何圖形以及其特征的一種常用方法?；舴蜃儞Q的思想是：對(duì)二值圖像I上任意一個(gè)邊緣點(diǎn)I(x,y)，如果它在以(a,b)為圓心，r為半徑的圓上，則滿足：

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

(4)

從(4)式可以看出，每個(gè)邊緣點(diǎn)I(x,y)映射成空間 (a,b,r)上一個(gè)二次曲面。將(4)式寫(xiě)成參數(shù) 形式：

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

(5)

其中，θ∈[0, 2π)。設(shè)待檢測(cè)圓的半徑r∈[R1,R2]，則按照霍夫變換，對(duì)任意一個(gè)邊緣點(diǎn)I(x,y)，θ，r分別以Δθ，Δr的步長(zhǎng)遍歷[0, 2π)和[R1,R2]，得到參數(shù)空間(a,b,r)的子空間，并對(duì)相應(yīng)的累加器單元AC(a,b,r)加1。對(duì)圖像上所有的點(diǎn)遍歷之后，得到累加器數(shù)組AC，其中任意一點(diǎn)AC(a,b,r)表示在圓(a,b,r)上的邊緣點(diǎn)的數(shù)目。相同數(shù)目的邊緣點(diǎn)在不同半徑的圓上被認(rèn)為是圓的可能性不同，半徑越小，可能性越大。因此定義圓度的概念：

p=n/r，

(6)

其中，n為圓的邊緣點(diǎn)的數(shù)目 ；r為圓的半徑；p為圓度。將累加器數(shù)組AC中的元素按(6)式計(jì)算圓度，得到三維圓度數(shù)組C。圓形目標(biāo)對(duì)應(yīng)于圓度數(shù)組的局部峰值，即對(duì)應(yīng)的C(a,b,r)為局部最大值。

3 結(jié) 果

3.1 識(shí)別流程

在嫦娥二號(hào)7 m全月數(shù)字高程模型數(shù)據(jù)中隨機(jī)選擇G227區(qū)、G228區(qū)進(jìn)行實(shí)驗(yàn)，基于ENVI和IDL，按照章節(jié)2中的原理和方法，編寫(xiě)程序進(jìn)行識(shí)別提取，如圖5。首先，設(shè)置合理的等高距在嫦娥二號(hào)數(shù)字高程模型數(shù)據(jù)中提取等高線；然后，通過(guò)圓形度、矩形度、體態(tài)比、剖面凹凸性篩選出初選撞擊坑區(qū)域；最后，通過(guò)霍夫變換提取撞擊坑邊界信息。將提取的信息導(dǎo)入ArcGIS，確定撞擊坑的地理位置和半徑大小，并通過(guò)高程差計(jì)算得到撞擊坑的深度信息，生成具有位置、大小、深度的小型撞擊坑數(shù)據(jù)庫(kù)。

圖5 IDL程序計(jì)算演示
Fig.5 IDL program diagram

3.2 識(shí)別結(jié)果

基于數(shù)字高程模型數(shù)據(jù)等高線自動(dòng)識(shí)別結(jié)果如圖6，具體統(tǒng)計(jì)結(jié)果見(jiàn)表1。

本文的方法在G227區(qū)、G228區(qū)分別識(shí)別撞擊坑728個(gè)和989個(gè)。為了評(píng)價(jià)本文方法的精度，使用Shufelt[20]提出的算法評(píng)價(jià)因子：檢測(cè)百分比為(Detection percentage)D=100TP/(TP+FN)，分支系數(shù)為B=FP/TP，質(zhì)量百分比為Q=100TP/(TP+FP+FN)。TP表示正確識(shí)別的撞擊坑數(shù)量；FN表示遺漏識(shí)別的撞擊坑數(shù)量；FP表示錯(cuò)誤識(shí)別的撞擊坑數(shù)量。目前最好的評(píng)價(jià)方法是采用人工識(shí)別的結(jié)果作為對(duì)比參照。

圖6 撞擊坑識(shí)別結(jié)果示意圖
Fig.6 The crater detection result diagram

表1撞擊坑識(shí)別結(jié)果統(tǒng)計(jì)表
Table1Craterdetectionresultstatistics

TP/個(gè)FN/個(gè)FP/個(gè)D/%B/%Q/%G227697533192.94.489.2G228918817191.97.781

從表1可知，本文方法在兩個(gè)實(shí)驗(yàn)區(qū)的評(píng)價(jià)因子D分別為92.9%、91.9%，均值為92.4%，B分別為4.4%、7.7%，均值為6.1%，Q分別為89.2%、81%，均值為85.1%。從評(píng)價(jià)因子可以看出，本文的方法對(duì)于小型撞擊坑有非常好的檢測(cè)結(jié)果，正確識(shí)別率高，且錯(cuò)分率低，存在小部分遺漏識(shí)別的情況。

3.3 討 論

本文提出一種基于高精度數(shù)字高程模型數(shù)據(jù)的派生數(shù)據(jù)等高線，通過(guò)計(jì)算機(jī)自動(dòng)識(shí)別提取小型撞擊坑的方法。相對(duì)于傳統(tǒng)的逐點(diǎn)分析坡度指數(shù)，提取圖像邊緣轉(zhuǎn)換成對(duì)地形線的分析，可以有效避免地形皺褶和碎石群的干擾，識(shí)別流程更為簡(jiǎn)單高效，識(shí)別速度較快。設(shè)置不同的等高距，分批次分區(qū)間識(shí)別不同尺度的撞擊坑，其正確識(shí)別率高，遺漏識(shí)別率低，說(shuō)明該方法是可行有效的。不過(guò)需要注意的是，對(duì)于部分退化嚴(yán)重的撞擊坑，坑緣破壞嚴(yán)重，坑壁平緩，等高線不能明顯表達(dá)和準(zhǔn)確勾畫(huà)出撞擊坑的形貌輪廓，無(wú)法被準(zhǔn)確識(shí)別，導(dǎo)致漏分情況的發(fā)生。但通過(guò)對(duì)比同位置的影像數(shù)據(jù)，通過(guò)人工識(shí)別分析是可以分辨的。因此，今后的研究方向和重點(diǎn)在于如何識(shí)別退化嚴(yán)重的撞擊坑，這也是所有撞擊坑識(shí)別方法中的薄弱環(huán)節(jié)和有待解決的問(wèn)題。在實(shí)際應(yīng)用中，可用該方法的識(shí)別結(jié)果補(bǔ)充現(xiàn)有月球撞擊坑數(shù)據(jù)庫(kù)中的小尺度撞擊坑，豐富現(xiàn)有的月球撞擊坑數(shù)據(jù)庫(kù)。