基于影像數(shù)據(jù)的馮·卡門(mén)撞擊坑形貌分析

鄭晨，姚鴻泰

（1. 河南大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院應(yīng)用數(shù)學(xué)研究所，開(kāi)封 475000；2. 中國(guó)科學(xué)院 國(guó)家天文臺(tái)，北京 100012）

0 引 言

從古時(shí)的夜觀天象到當(dāng)今的“嫦娥”飛天，人類一直未停止過(guò)對(duì)浩瀚宇宙的觀測(cè)步伐。而距離人類最近的天體——月球，也成為了人類探測(cè)的重要對(duì)象。從20世紀(jì)50年代美國(guó)和前蘇聯(lián)發(fā)射的探測(cè)器飛掠月球，到實(shí)現(xiàn)人類首次登月的“阿波羅號(hào)”，再到我國(guó)開(kāi)展的“嫦娥系列”探月工程，人類對(duì)月球的認(rèn)知在不斷深入和細(xì)化。

2013年12月15日，我國(guó)“嫦娥3號(hào)”著陸器在虹灣地區(qū)順利著陸，成為了世界上第3個(gè)成功實(shí)施月球軟著陸的國(guó)家[1-2]，也實(shí)現(xiàn)了我國(guó)探月工程“繞、落、回”中“落”的目標(biāo)。但是，目前成功軟著陸于月球的人類探測(cè)器主要位于月球正面，還沒(méi)有在月球背面成功著陸的先例。2014年以來(lái)，中國(guó)探月與航天工程中心經(jīng)多方論證，確定了“嫦娥4號(hào)”著陸月球背面這一任務(wù)目標(biāo)和相應(yīng)的技術(shù)方案[3-4]，有望實(shí)現(xiàn)人類首次月球背面著陸及勘察。著陸區(qū)域則初步選定為月球背面南極–艾肯（South Pole-Aitken，SPA）盆地內(nèi)的馮·卡門(mén)（Von Kármán）撞擊坑，該區(qū)域?qū)儆诟叩氐匦蝃2]，為“嫦娥4號(hào)”著陸器的軟著陸帶來(lái)了困難。

著陸區(qū)地形地貌的分析是影響著陸設(shè)計(jì)的一個(gè)重要方面。為了輔助“嫦娥4號(hào)”探月工程的開(kāi)展，本文根據(jù)已有馮·卡門(mén)撞擊坑的影像觀測(cè)數(shù)據(jù)，利用基于隨機(jī)場(chǎng)的聚類分析方法，從聚類的角度對(duì)該區(qū)域的地形地貌進(jìn)行了分析與表示，以期為“嫦娥4號(hào)”著陸區(qū)的選擇提供有意義的參考，為探測(cè)任務(wù)的順利開(kāi)展提供更多的技術(shù)支持。

1 數(shù)據(jù)與方法

1.1 實(shí)驗(yàn)數(shù)據(jù)

主要采用中心點(diǎn)坐標(biāo)為（44.8°S，175.9°E），40～48°S，172°E～180°～178°W范圍內(nèi)的LOLA（Lunar Orbiter Laser Altimeter）高程DEM數(shù)據(jù)[5]和LROC（Lunar Reconnaissance Orbiter Camera）的高分辨率影像數(shù)據(jù)[6]（圖1），對(duì)馮·卡門(mén)撞擊坑的內(nèi)部形貌進(jìn)行分析。

馮·卡門(mén)撞擊坑屬于SPA內(nèi)的典型地貌類型，其內(nèi)部的月球物質(zhì)成分具有代表性[2,7-8]，對(duì)月球火山及月殼活動(dòng)的研究具有重要價(jià)值。同時(shí)，它位于月球背面，是開(kāi)展低頻射電天文觀測(cè)的理想地點(diǎn)。而且，從LOLA高程數(shù)據(jù)和LROC影像數(shù)據(jù)的直方圖（見(jiàn)圖2）可以發(fā)現(xiàn)，馮·卡門(mén)撞擊坑的高程數(shù)據(jù)（圖2（a））主

要分布在–6～–1 km之間，且撞擊坑內(nèi)部的起伏變化較??；而LROC的影像數(shù)據(jù)值呈近似正態(tài)分布（圖2（b）），坑內(nèi)的地形地貌整體也較為平坦。因此，該區(qū)域是“嫦娥4號(hào)”初步選定的著陸區(qū)之一[4]。但是，從圖1（b）可以觀測(cè)到，撞擊坑內(nèi)仍存有大量的局部中小型環(huán)狀撞擊坑區(qū)域，在直方圖中也反映出異常的峰值點(diǎn)（圖2紅色虛線框標(biāo)注）。由于LOLA高程數(shù)據(jù)和LROC影像數(shù)據(jù)的灰度值分布相對(duì)集中，造成局部形貌的數(shù)值變化較小。這種局部低對(duì)比度現(xiàn)象將直接影響人工判讀的精度，對(duì)分析撞擊坑內(nèi)部的形貌造成潛在的影響。

圖1 馮·卡門(mén)撞擊坑的實(shí)驗(yàn)數(shù)據(jù)Fig. 1 Test data of the Von Kármán crater

1.2 實(shí)驗(yàn)方法

為解決觀測(cè)數(shù)據(jù)中低對(duì)比度問(wèn)題，提高觀測(cè)數(shù)據(jù)的可讀性，本文研究了基于馬爾科夫隨機(jī)場(chǎng)模型（Markov Random Field，MRF）[9]的聚類分析方法。根據(jù)已有研究，聚類分析是提高深空探測(cè)數(shù)據(jù)局部低對(duì)比度的一種有效方法[10]，而MRF模型是一類利用概率圖模型進(jìn)行聚類的方法。它通過(guò)在觀測(cè)數(shù)據(jù)和聚類標(biāo)記間建立概率模型，不僅可以凸顯特征相近的不同地物形貌間的差異，而且還可以刻畫(huà)相鄰像素的空間關(guān)系，提高模型的紋理描述能力。

圖2 馮·卡門(mén)撞擊坑LOLA和LORC數(shù)據(jù)的直方圖Fig. 2 Histogram of LOLA and LROC data for the Von Kármán crater

具體而言，假設(shè)大小為的觀測(cè)影像數(shù)據(jù)定義在柵格位置集合上，其中表示像素處的觀測(cè)數(shù)值。標(biāo)記隨機(jī)場(chǎng)也定義在位置集合上，其中每個(gè)表示像素處的聚類類別。此處，是一個(gè)隨機(jī)變量，且，其中是聚類的類別總數(shù)。若是標(biāo)記隨機(jī)場(chǎng)的一個(gè)實(shí)現(xiàn)，那么，聚類問(wèn)題在MRF模型中就被轉(zhuǎn)化為最大后驗(yàn)概率的求解問(wèn)題，即

上式中是實(shí)現(xiàn)的集合，即根據(jù)貝葉斯公式，式（1）可表示為

其中第2個(gè)等號(hào)成立的原因是由于是已知的觀測(cè)數(shù)據(jù)，不影響最終的聚類結(jié)果。所以，在MRF模型的聚類過(guò)程中，需要首先確定似然函數(shù)和標(biāo)記場(chǎng)聯(lián)合概率分布的具體函數(shù)形式。

對(duì)于聯(lián)合概率MRF模型假設(shè)標(biāo)記隨機(jī)場(chǎng)具有空間馬氏性，即其中表示與像素空間相鄰的像素點(diǎn)組成的集合。由Hammersley–Clifford定理可知[9]，此時(shí)服從Gibbs分布

上式中是歸一化常數(shù)，是實(shí)現(xiàn)的能量函數(shù)。其中，是定義在勢(shì)團(tuán)上的勢(shì)函數(shù)。本文模型中，為了在保持聚類精度的前提下降低計(jì)算量，我們只考慮具有二階勢(shì)團(tuán)的多層邏輯模型[11-12]，此時(shí)勢(shì)函數(shù)為

對(duì)于似然函數(shù)假設(shè)在給定標(biāo)記隨機(jī)場(chǎng)實(shí)現(xiàn)的條件下，各像素點(diǎn)的觀測(cè)特征是相互獨(dú)立的，即根據(jù)圖2的直方圖，尤其是LROC的影像數(shù)據(jù)分布，本文模型采用正態(tài)分布來(lái)描述各像素點(diǎn)的似然函數(shù)分布，即在的條件下其中n是觀測(cè)數(shù)據(jù)的波段數(shù)，參數(shù)分別是類別的正態(tài)均值和方差。

在確定了似然函數(shù)和標(biāo)記場(chǎng)的聯(lián)合概率后，通過(guò)不動(dòng)點(diǎn)的思路來(lái)逐像素迭代更新（2）式的后驗(yàn)概率，實(shí)現(xiàn)最終的聚類分析，即

具體的算法請(qǐng)見(jiàn)表1。其中，關(guān)于式（7）中的正態(tài)分布參數(shù)，可以根據(jù)EM算法[13]估計(jì)如下

式中，表示類別為的像素點(diǎn)的數(shù)目。

表1 MRF算法流程Table 1 Algorithm of the MRF model

③根據(jù)公式（9），以及標(biāo)記實(shí)現(xiàn)更新(7)式的參數(shù)及似然函數(shù)

④根據(jù)公式（5）和標(biāo)記實(shí)現(xiàn)更新（6）式的標(biāo)記聯(lián)合概率

⑤根據(jù)公式（8）和最大后驗(yàn)準(zhǔn)則，更新像素處的標(biāo)記

⑥更新標(biāo)記實(shí)現(xiàn)則輸出作為聚類結(jié)果；否則，令t=t+ 1，返回步驟③。

2 實(shí)驗(yàn)結(jié)果

在上述的MRF模型，定義鄰域集合為8鄰域，分別對(duì)LROC影像數(shù)據(jù)和LOLA高程數(shù)據(jù)進(jìn)行了聚類分析，具體結(jié)果如下。

2.1 LROC實(shí)驗(yàn)結(jié)果

圖3 LROC影像數(shù)據(jù)K = 3時(shí)的聚類結(jié)果和對(duì)應(yīng)直方圖Fig. 3 Clustering results of LROC data with K = 3 and its corresponding histogram.

對(duì)于LROC影像數(shù)據(jù)，首先在K= 3的類別數(shù)目下進(jìn)行了聚類分析，結(jié)果如圖3（a）所示。此時(shí)， MRF模型的結(jié)果主要分為了撞擊坑內(nèi)較平坦的區(qū)域、向光區(qū)域和陰影區(qū)域這3種類別，它們對(duì)應(yīng)的影像灰度值分別位于圖2（b）中的峰值和兩邊的尾部區(qū)域，如圖3（b）的直方圖和表2的定量數(shù)值所示。由于馮·卡門(mén)撞擊坑內(nèi)的中小型撞擊坑具有明顯的向光和陰影區(qū)域，因此根據(jù)該結(jié)果，可以更直觀快速地判讀出這些局部撞擊坑。然而，不同撞擊坑間的差異卻難以得到更精細(xì)的表示。

表2 K = 3時(shí)影像各類別的定量指標(biāo)Table 2 Quantitative indexes of images with K = 3

為了進(jìn)一步凸顯不同撞擊坑間的差異，考察了K=7時(shí)的聚類結(jié)果，如圖4所示。此時(shí)聚類結(jié)果中每種類別對(duì)應(yīng)的地形地貌被進(jìn)一步細(xì)化，其中坑內(nèi)的平坦區(qū)域被細(xì)分為兩類（橙色和粉色區(qū)域）、向光區(qū)域被細(xì)分為3種不同程度的類別（青色、黃色和綠色區(qū)域）、陰影區(qū)域則被細(xì)分為藍(lán)色和紫色兩種類別。同時(shí)，相比于K= 3，此時(shí)每種類別內(nèi)的標(biāo)準(zhǔn)差會(huì)相應(yīng)降低，各類別的灰度均值則分別位于圖2（b）中的一些局部異常峰值點(diǎn)附近，如表3所示。隨著聚類的細(xì)化，圖4的聚類結(jié)果可以描述出不同中小型撞擊坑間的差異。

圖4 LROC影像數(shù)據(jù)K = 7時(shí)的聚類結(jié)果Fig. 4 Clustering results of LROC data with K = 7

相比于LROC影像數(shù)據(jù)的直接判讀，結(jié)合圖3和圖4的聚類結(jié)果，馮·卡門(mén)撞擊坑內(nèi)的地形地貌，尤其是撞擊坑等局部低對(duì)比度區(qū)域，可以得到更合理的表示，如圖5所示。具體而言，對(duì)于該區(qū)域內(nèi)較大型的撞擊坑，它的向光與陰影區(qū)域是漸變的，因此在K= 3的聚類結(jié)果中會(huì)從左至右地表現(xiàn)出向光區(qū)域與陰影區(qū)域的交替（一般都為兩次，分別對(duì)應(yīng)坑的邊緣和內(nèi)部光照變化）；而在K= 7的聚類結(jié)果中，大型撞擊坑會(huì)含有多種類別的聚類結(jié)果，如在圖5（a）中示例的撞擊坑就由兩種陰影區(qū)域類別和3種向光區(qū)域類別共同構(gòu)成。對(duì)于坑內(nèi)的中型撞擊坑，由于坑的邊緣不是很明顯，在K= 3的聚類結(jié)果中通常只表現(xiàn)出坑內(nèi)部的一次光照變化情況，如圖5（c）所示，而K= 7的聚類結(jié)果中出現(xiàn)的類別數(shù)目也會(huì)少于大型撞擊坑。對(duì)于小型和微型撞擊坑，由于其直徑很小，導(dǎo)致在K= 3的聚類結(jié)果中只表現(xiàn)出小塊的陰影類別，其在K= 7的聚類結(jié)果中也一般只含有兩類表示陰影區(qū)域的類別，如圖5（b）所示。除了上述3種不同大小的環(huán)形撞擊坑，馮·卡門(mén)撞擊坑內(nèi)還具有一些長(zhǎng)條型的凹陷區(qū)域，本文稱之為長(zhǎng)型撞擊坑。這種類型的撞擊坑在光照變化方面與中型撞擊坑類似，一般只含有一次的陰影與光照區(qū)域交替，所以類別特征的分布是相似的，其差別主要體現(xiàn)在坑的形狀上。根據(jù)上述特征，馮·卡門(mén)撞擊坑內(nèi)典型形貌的統(tǒng)計(jì)結(jié)果總結(jié)如表4所示（注：因小型撞擊坑數(shù)目過(guò)多，因此采用程序自動(dòng)地統(tǒng)計(jì)了聚類結(jié)果，其數(shù)目顯示為3 994，因程序統(tǒng)計(jì)可能存在一定誤差，所以下表記為約4 000個(gè)）。

表3 K = 7時(shí)影像各類別的定量指標(biāo)Table 3 Quantitative indexes of images with K = 7

圖5 不同類型撞擊坑在K = 3和K = 7時(shí)的特征Fig. 5 Characteristics of different craters with K = 3 and K = 7

表4 馮·卡門(mén)撞擊坑內(nèi)各典型形貌的統(tǒng)計(jì)結(jié)果Table 4 Statistical results of various craters in the Von Kármán crater

2.2 LOLA實(shí)驗(yàn)結(jié)果

對(duì)于LOLA高程數(shù)據(jù)，由于馮·卡門(mén)撞擊坑內(nèi)數(shù)據(jù)差異較小，導(dǎo)致整個(gè)區(qū)域的可視化特征非常接近，如圖6（b）所示。這種低對(duì)比度的影像數(shù)據(jù)不能有效地凸顯撞擊坑內(nèi)不同地物的差異，而一些在LROC數(shù)據(jù)中可以觀測(cè)到的典型形貌也難以在高程數(shù)據(jù)中進(jìn)行判讀。為了描述不同地物形貌在LOLA高程數(shù)據(jù)中的差異，我們?nèi)圆捎肕RF模型對(duì)其進(jìn)行聚類分析，具體結(jié)果請(qǐng)見(jiàn)圖6。

從圖6可以發(fā)現(xiàn)，當(dāng)K值設(shè)定較小，如圖6（c）中K= 3時(shí)，除了坑中部小型凸峰的峰值區(qū)域外，整個(gè)撞擊坑區(qū)域被劃分成一種類別，此時(shí)聚類結(jié)果主要表現(xiàn)出數(shù)據(jù)的宏觀特征，即坑內(nèi)區(qū)域整體較平坦。當(dāng)進(jìn)一步增大K值時(shí)，如圖6（d）中K= 5時(shí)，坑內(nèi)的典型地物，如中部的凸峰就可以得到較完整的顯示。當(dāng)K= 10時(shí)，聚類結(jié)果中可以顯示出該區(qū)域內(nèi)整體的高程變化和一些較深的中小型撞擊坑。而當(dāng)K= 22時(shí)，一些從圖6（b）中難以目視判讀的小型撞擊坑，也可以在結(jié)果中得到表現(xiàn)；同時(shí)，整個(gè)馮·卡門(mén)撞擊坑內(nèi)的高程變化也可以細(xì)分為左上部分和右下部分兩種類別，這和圖4中橙色和粉色的平坦區(qū)域劃分是相對(duì)應(yīng)的。最后，當(dāng)K= 28時(shí)，坑內(nèi)右下部分中一塊略有凹陷的區(qū)域可以得到表現(xiàn)。

縱觀上述結(jié)果，通過(guò)MRF模型的聚類分析，LOLA高程數(shù)據(jù)中的低對(duì)比度區(qū)域可以得到更合理的表示，一些典型地物形貌也能在聚類結(jié)果中得以彰顯。

因此，相比于目視判讀，基于聚類的表示方法可以有效地提高低對(duì)比度區(qū)域的可視化水平，能輔助相應(yīng)區(qū)域的地形分析。同時(shí)，本文使用的馬爾科夫隨機(jī)場(chǎng)模型因其有效的空間描述能力，可以進(jìn)一步優(yōu)化聚類的結(jié)果。為了突顯本文模型的特點(diǎn)，在圖7中對(duì)比了本文模型和經(jīng)典Kmeans算法[10]在LOLA高程數(shù)據(jù)的聚類結(jié)果，可以看出，在K= 5的條件下，經(jīng)典Kmeans算法因?yàn)閮H考慮了各像素的高程數(shù)據(jù)，沒(méi)有考慮不同像素間的空間關(guān)系，所以它的聚類結(jié)果（圖7（a））沒(méi)有完整地識(shí)別出中央的凸峰。本文采用的馬爾科夫隨機(jī)場(chǎng)模型通過(guò)似然函數(shù)和聯(lián)合概率分別考慮了影像數(shù)據(jù)中的像素高程數(shù)據(jù)和各像素間的空間關(guān)系，更充分地使用了影像信息，所以聚類結(jié)果圖7（b）中較完整地識(shí)別出了中央的凸峰。

圖6 LOLA高程數(shù)據(jù)聚類結(jié)果Fig. 6 Clustering result of LOLA data

圖7 本文方法和Kmeans方法的聚類結(jié)果對(duì)比Fig. 7 Comparison of clustering results between our method and the Kmeans

3 結(jié) 論

本文基于LOLA高程數(shù)據(jù)和LROC高空間分辨率影像數(shù)據(jù)，通過(guò)馬爾科夫隨機(jī)場(chǎng)模型對(duì)月球背面SPA內(nèi)的馮·卡門(mén)撞擊坑進(jìn)行了地形地貌的聚類表示和分析。本文的貢獻(xiàn)主要體現(xiàn)在以下兩點(diǎn)：①模型的似然函數(shù)刻畫(huà)了影像數(shù)據(jù)的概率分布，而標(biāo)記隨機(jī)場(chǎng)則描述了地形地貌間的空間關(guān)系，結(jié)合這兩點(diǎn)的本文模型通過(guò)聚類可以有效地突顯出影像數(shù)據(jù)中低對(duì)比度區(qū)域內(nèi)地形地貌的差異；②通過(guò)分析聚類結(jié)果可知，馮·卡門(mén)撞擊坑的地形從整體來(lái)看，除了在中上部有一個(gè)凸峰外，其余區(qū)域是相對(duì)平坦的，高低起伏也較緩慢；從局部來(lái)看，該區(qū)域內(nèi)仍存在大量的環(huán)狀撞擊坑。而這些撞擊坑又可以根據(jù)聚類結(jié)果分為大型、中型、小型和長(zhǎng)型撞擊坑，其中大型撞擊坑不僅深度較大，而且坑的邊緣還有明顯的地勢(shì)起伏，但數(shù)量相對(duì)較少；中型和長(zhǎng)型撞擊坑的邊緣相對(duì)平緩，但是撞擊坑自身仍具有一定的深度，且分布于馮·卡門(mén)撞擊坑內(nèi)部各處；而小型撞擊坑的高程變化較小，但數(shù)量卻很多，廣布于擬著陸區(qū)域。在工程可控的前提下，著陸區(qū)域可考慮選擇凸峰西側(cè)的下方。因?yàn)榇颂幉粌H相對(duì)平坦，而且還接近各類型撞擊坑，便于展開(kāi)后續(xù)的科學(xué)探測(cè)工作。

致謝

本文實(shí)驗(yàn)LROC數(shù)據(jù)和LOLA數(shù)據(jù)分別從Java Missionplanning and Analysis for Remote Sensing（JMARS）https://jmars.mars.asu.edu/和網(wǎng)站http://imbrium.mit.edu/DATA/LOLA_GDR/CYLINDRICAL/FLOAT_IMG/下載，在此表示感謝！

參考文獻(xiàn)

[1]SUN Z Z，JIA Y，ZHANG H. Technological advancements and promotion roles of Chang’e-3 lunar probe mission[J]. Science China（Technological Sciences），2013，43（11）：2702-2708.

[2]李飛，張熇，吳學(xué)英，等. 月球背面地形對(duì)軟著陸探測(cè)的影響分析[J].深空探測(cè)學(xué)報(bào)，2017，4（2）：143-149.LI F，ZHANG H，WU X Y，et al. Influence analysis of terrain of the farside of the Moon on soft-landing[J]. Journal of Deep Space Exploration，2017，4（2）：143-149.

[3]YE P J，SUN Z Z，ZHANG H，et al.An overview of the mission and technical characteristics of Change’4 Lunar Probe [J]. Science China（Technological Sciences），2017，60（5）：658-667.

[4]吳偉仁，王瓊，唐玉華，等. “嫦娥4號(hào)”月球背面軟著陸任務(wù)設(shè)計(jì)[J].深空探測(cè)學(xué)報(bào)，2017，4（2）：111-117.WU W R，WANG Q，TANG Y H，et al. Design of Chang’e-4 lunar farside soft-landing mission [J]. Journal of Deep Space Exploration，2017，4（2）：111-117.

[5]SMITH D E，ZUBER M T，JACKSON G B，et al. The Lunar orbiter laser altimeter investigation on the lunar reconnaissance orbiter mission[J]. Space Science Reviews，2010，150 （1-4）：209-241.

[6]ROBINSON M S，BRYLOW S M，TSCHIMMEL M，et al. Lunar reconnaissance orbiter Camera（LROC）instrument overview [J]. Space Science Reviews，2010，150 （1-4）：81-124.

[7]PETERSON C A，HAWKE B R，BLEWETT D，et al. Geochemical units on the Moon：the role of South Pole-Aitken basin[C]//Lunar and PlanetaryScience XXXIII. USA：[s.n.]，2002.

[8]Peterson C A，Hawke B R，Lucey P G，et al. Anorthosite on the lunar farside and its relationship to South Pole-Aitken Basin[C]//Lunar and Planetary Science XXXI. USA：[s.n.]，2000.

[9]LI S Z. Markov random field modeling in computer vision. 3ed ed.[M]. New York：Springer-Verlag，2009.

[10]ZHENG C，PING J S，WANG M Y. Hierarchical classification for the topography analysis of Asteroid（4179）Toutatis from the Chang’E-2 images [J]. ICARUS，2016，278：119-127.

[11]DERINH，ELLIOTT H. Modeling and segmentation of noisy and textured images using gibbs random fields [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence，1987，9（1）：39-55.

[12]DERINH，COLE W S. Segmentation of textured images using Gibbs random fields [J]. Computer Vision，Graphics and Image Processing，1986，35：72-98.

[13]DEMPSTER A P，LAIRD N M，RUBIN D B. Maximum likelihood from incomplete data via the EM algorithm [J]. Journal of the Royal Statistical Society（series b），1977，39（1）：1-38.

[14]JAIN A K，MURTY M N，F(xiàn)LYNN PJ. Data clustering：A review [J].ACM Computing Surveys，1999，31（3）：264-323.