鄭晨，姚鴻泰

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

0 引 言

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

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

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

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

1.1 實驗數(shù)據(jù)

主要采用中心點坐標為（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），對馮·卡門撞擊坑的內(nèi)部形貌進行分析。

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

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

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

1.2 實驗方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

⑤根據(jù)公式（8）和最大后驗準則，更新像素處的標記

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

2 實驗結(jié)果

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

2.1 LROC實驗結(jié)果

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

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

表2 K = 3時影像各類別的定量指標Table 2 Quantitative indexes of images with K = 3

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

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

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

表3 K = 7時影像各類別的定量指標Table 3 Quantitative indexes of images with K = 7

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

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

2.2 LOLA實驗結(jié)果

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

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

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

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

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

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

3 結(jié) 論

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

致謝

本文實驗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/下載，在此表示感謝！

參考文獻

[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é)英，等. 月球背面地形對軟著陸探測的影響分析[J].深空探測學(xué)報，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號”月球背面軟著陸任務(wù)設(shè)計[J].深空探測學(xué)報，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.