基于數字高程模型的月海地貌信息重建

禹 健，郝秀芳，安永泉

(1. 山西大學 自動化系，山西 太原 030013； 2. 中北大學 信息與通信工程學院，山西 太原 030051)

美國1994年1月25日發射Clementine探測器[1-3]，標志著月球探測熱潮的二次興起，歐空局和日本等也加入了月球探測的行列[4，5]. 歐空局的智能1號(Smart-1)、 日本的月亮女神(Selene )、 中國的嫦娥1號至4號(Chang’e-1至4)號、 印度的月球初航1號至3號(Chandrayaan-1至3)等先后飛向月球，獲得了大量寶貴月球影像、 地形等數據. 其中高程類數據反映了月表地形的結構特征.

月海地貌重建考慮的主要特征為各年齡撞擊坑. 關于月表撞擊坑自動識別，日本學者Sawabe等2006年基于Apollo照片，得到了約80%的撞擊坑[6]； 中科院遙感所的岳宗玉[7](2008年)提出了面向對象分類方法的月表撞擊坑識別； 李超等(2012年)提出了一種新的橢圓形檢測方法[8]； 魯宇航等(2013年)提出用最大類間方差法對月面影像分類，并擬合邊緣[9]； 江泓昆等(2013年)提出一種基于特征空間的撞擊坑自動識別的自適應算法[10]； 王棟等(2016年)提出一種將等值線分析與球面窗口掃描相結合的識別算法[11].

國內學者的坑識別方法所得結果，根據坑個數和區域地質年代之間的經驗公式進行求證，基本合理[12]. 其中有兩個問題尚待解決，一是同一撞擊坑由于光照角度不同可能被分識成幾個撞擊坑而多次統計； 二是大型不規則撞擊坑的識別仍需要在ArcGIS中手動修訂，即需把屬于某大坑的目標體合并.

在本文的工作中，基礎數據為Chang’e-2激光高度計所獲取的、 中國科學院國家天文臺后期處理的LEVEL0A數字高程數據(DEM)，精度為500 m×500 m. 提出的基點彌散法，實現不規則撞擊坑的尋找，能夠識別坑唇細節，給出坑的位置，坑深，長軸直徑，短軸直徑，給出了不同地形類型的統計規律，結合坡度坡向約束，重建海地貌仿真地形. 與NASA公布的結論對比，運用經驗公式進行驗證，該技術可用于探索更多月表區域.

1 DEM數據預處理

“嫦娥二號”衛星攜帶有8套24件科學探測儀器，包括： CCD立體相機、 激光高度計、 伽馬X射線譜儀、 微波探測儀、 太陽高能粒子探測器和低能離子探測器，在月球探測測控系統和地面應用系統支持下，這些有效載荷返回了月表地貌、 物質成分、 內部結構和天文環境的科學數據. 本文采用的激光高度計數字高程數據，是進行源包排序、 優化拼接、 去重復、 去源包包頭后的有效載荷科學數據塊，包含一定分辨率的月面地形信息. 首先將數據塊進行投影變換、 插值細化和拼接的預處理工作，然后得到坡度坡向信息，進而進行隕石坑的識別.

1.1 投影校正

投影變換是對非水平基準地貌遙測數據必不可少的預處理工作.

源數據為墨卡托(Mercator)，即正軸等角圓柱投影. 將球面展開成平面，建立二維平面和三維球面的對應關系. 所有數據塊變換成正軸投影面的正射投影. 即投影平面切于月球面上一點，視點在無限遠，投影光線為互相平行的直線，與投影平面相垂直. 所有緯線圈無長度變形，投影中心一定范圍內經線圈可確保不失真.

采用同一坐標系(球面坐標系)下正解變換法克服投影變形，如圖1 所示.

圖1 參考基準球坐標變換Fig.1 The reference spherical coordinates transform

已知球面上兩點(lon1,lat1),(lon2，lat2). 此兩點所在大圓的方位角為(alon,alat)，使用的坐標均為地理坐標系坐標.

tan(lat)=-cos(alon-lon)/tan(alat),

(1)

式中：AN=tan(b),ON=1/cos(b), tan(c)=AM/OA,AM=AN/cos(∠NAM).

方位角參數的方程為

(2)

已知點(lon1,lat1)，在方位角為(alon,alat)的方位投影坐標系中的極坐標參數

ctg(aplon)=ctg(alon-lon1)/sin(alat).

(3)

需要分段來處理，分4段.

tg(aplon)=tg(alon-lon1)*sin(alot),

ctg(aplon)=ctg(alon-lon1)/sin(alat),

(4)

投影變換后進行經典的雙線性內插，采用坐標點對齊結合邊界異常點剔除進行拼接.

1.2 坡面參數計算

采用D8法[13]給出所處理區域每一點e的最大坡降方向，如圖2 所示.

圖2 最大坡降方向求取窗口Fig.2 The window for calculating maximum gradient direction

在3×3的格網局部窗口中，設中心格網為e，只有其周圍8個格網點對應的8種可能坡降方向,每個方向間隔45°. 判斷條件:

max{k×(zc-zi)}i=1,2,…,8.

(5)

圖3 最大坡向矩陣確定算法流程Fig.3 The algorithm process of maximum slope directionmatrix calculation

2 月球坑識別

月表地形醒目的結構特征是布滿環形構造(環形體). 按照成因機制，分為月海穹窿和撞擊坑. 月海穹窿是月海地區呈環形或橢圓形出現、 中間向上稍隆起的一種表面平滑而坡度較小的正地形[14]. 撞擊坑遍布全月，是月面最主要的地貌形態，有相對平坦的底部和撞擊事件發生后大量濺射物堆積形成的坑唇. 為識別需要，本文給出月球坑的廣義定義，將月海穹窿作為一種特殊的，坑深為負值，坑唇為零的坑包含在識別結果內.

2.1 定義

月球坑定義為獨立洼地. 長寬比率φv=d1(v)/d2(v)大于閾值(0.6～0.8). 其中d1(v)為第v個撞擊坑進行橢圓擬合時橢圓長軸直徑，d2(v) 為第v個撞擊坑進行橢圓擬合時橢圓短軸直徑. 對于年輕坑和成熟坑，在某一高程平面上可有凸輪廓(坑唇).

圖4 月球地形影像數據視圖Fig.4 Image data view of lunar terrain

圖4 中： ①為撞擊坑(基斯Kies)(直徑45 km，比例尺： 1∶3 000 000)，②為月球表面的穹窿構造，識別結果將體現為坑深為負值的撞擊坑.

當d1(v)或d2(v)的值過大，則為一片低洼的地區，不是坑. φv值過大，則為溝狀地形，月谷或月溪，不是坑.

分辨率為500 m×500 m的DEM數據，讀出后顯示的視覺效果如圖5 所示.

圖5 月球地形DEM數據視圖Fig.5 DEM data view of lunar terrain

2.2 基點彌散法月球坑自動識別

月球坑識別采用基點彌散法進行自動逐組尋找.

基礎工作為兩項： 確定所處理區域的坑唇邊緣點集合； 確定所處理區域的最大坡降矩陣，即每點的最大坡降方向.

先尋找整個數據文件中的最低點，標記，作為第一個獨立洼地的尋找依據； 然后運用基點彌散尋找所有滿足閾值條件的點，終止條件是包含高程值變化量兩次過0； 識別過程中計算記錄此坑的參數(中心點，長軸半徑，短軸半徑，坑深等).

后續工作是坑的“隱藏”. 第一組已識別坑點被重新賦值后，進行第二組坑的識別. 尋找最低點，重復運用基點彌散法，逐組將坑識別出來. 整體技術構圖如圖6 所示.

圖6 基點彌散法組成Fig.6 Structure of the seed dispersion method

2.2.1 坑唇邊緣點確定

坑唇邊緣點定義為兩個相互正交的方向上，一個方向凸起，而另一個方向沒有凹凸性變化的點.

(6)

坑唇線是坑唇點的集合. 利用DEM數據提取地面的平面曲率及地面的正負地形，取正地形上平面曲率的最大值為坑唇點. 求取原始DEM數據層的最大高程值H，計算(H-DEM)得到與原來地形相反的DEM數據層，即反地形DEM； 若以坡向變率(SOA)表征平面曲率，具體步驟如圖7 所示，坑唇線求取結果如圖8 所示.

圖7 坑唇點確定算法結構框圖Fig.7 The algorithm structure of searching crater lip points

圖8 坑唇線求取結果示意Fig.8 The result of crater lip lines

2.2.2 非規則邊界月球坑的檢測

掃描尋找最低點，確定坑最深點(i0,j0). 以(i0,j0)為中心建立窗口，掃描窗口內的所有柵格點，將滿足下列條件的柵格點標上坑標記: 1) 柵格點高程小于其8鄰域，稱為特征點； 2) 與特征點相鄰； 3) 高程值不低于相鄰的特征點.

首先將一個基點入集，并在標志矩陣中做已處理標示； 若集不空，對頭格網點出集，在標志矩陣中做坑標示； 否則判斷該格網點的8鄰域格網點，如果最大坡降方向指向該格網，并且尚未處理過，則入集，并在標志矩陣中做已處理標示;循環增加標志矩陣，終止條件為： 1) 窗口包含足夠多高程值為0的點(經驗閾值)，因為古老撞擊坑邊緣被破壞，坑唇線檢測不出； 2) 窗口包含足夠多坑唇點(經驗閾值)，記錄窗口大小.

2.2.3 已識認坑的隱藏

完成坑的定位，參數記錄，邊界確定工作后，將此坑隱藏，以尋找下一個坑.

1) 獨立坑的隱藏

對于檢測出的所有坑特征點，將其高程值賦以鄰域格網的最小高程值，便可將所有單點坑隱藏. 將區域內高程值低于坑唇點閾值高程值的所有點的高程用坑唇點的高程代替.

2) 復合坑的隱藏

當出現復合坑，坑套坑時，利用復合坑間的指向關系以及在坑檢測過程中所得到的坑矢量特征，將柵格操作與矢量操作結合起來進行處理[11].

首先用鄰接表1 來表示坑間的鄰接關系，其間指向關系可通過查表確定，有坑合并時更新.

根據檢測過程中得到的坑間指向關系，復合坑中構成環狀指向的坑被測出，如圖9 中I, II, III坑； 將環中的坑合并為一個新坑(如圖10 中的洼地V),生成的新坑是獨立坑(沒有指向其他坑)，則可按照處理獨立坑的方法將其填平； 尋找、 合并至復合坑處理完畢.

表1 坑更新前后的鄰接關系表

圖9 復合坑及其指向關系Fig.9 Composite craters and their mutual pointing relation

圖10 一次合并后的復合坑及其指向關系Fig.10 Composite craters and their mutual pointing relationafter once merged

2.3 坑參數計算

坑參數包括坑中心點經緯度坐標(x,y)，因為帶中央峰的撞擊坑最低點并不是中心點，所以坑中心點根據坑長短直徑中點確定. 表征坑位置，坑深depth，坑長軸直徑a，坑短軸直徑b，坑唇傾斜率θ(指坑識認出后坑唇擬合橢圓與當地正北方向的傾斜角，如圖11 所示).

圖11 擬合橢圓的傾斜角(月球坑俯視圖)Fig.11 The tilt Angle of fitting ellipse(Vertical viewof lunar craters)

3 坑識別結果

對任意區域撞擊坑進行自動識別，并統計撞擊坑個數. 盡管所得到的坑很不規則，但對于個數統計不造成影響.

3.1 識別結果

月球坑分批識別，如果彈出的坑滿足要求，可通過調節坑深閾值中斷，如果窮盡識別，最終可以識別出所有可認的小坑、 碎坑. 視覺上可能會連成片，即每一片不規則形狀都是坑的情況. 但數據庫中是逐一記錄的，區分清晰，結果為新的圖層，其上任意坑內單擊，啟動數據庫鏈接，得到該坑的特征參數計算結果如圖12～圖15 所示. 坑的特征參數包括坑的名稱、 深度、 長軸直徑、 短軸直徑、 坑中心點坐標、 坑唇傾斜角.

圖12 坑識別結果示意Fig.12 Result of Craters identification

圖13 坑組Ⅰ與坑組Ⅲ的識別Fig.13 Identification of pit group Ⅰ and Pit group Ⅲ

圖14 坑組Ⅱ識別和組內的某單坑識別Fig.14 Identification of pit group Ⅱ and one single crater in it

若如坑組Ⅰ中坑套坑的情況，參數計算時會出現如圖15 的多個中心點橫坐標centerx、 中心點縱坐標centery、 坑深depth、 長軸直徑a、 短軸直徑b、 長軸直徑傾斜角angle.

圖15 坑組Ⅰ的參數識別結果Fig.15 Result parameters of pit group Ⅰ

此時，坑套坑的情況能夠識別出來，但參數區分很小的坑也將被區分成兩個，使坑的個數變多. 通過調節坑唇點被包含閾值合并組坑.

表2～表4 給出部分識認出的坑指標，位置特征相近度，直徑特征相近在±3 km，推定以知名坑為識別結果驗證基準.

表2 識別結果1

表2 中識別的坑參數以Abbe坑(1970年命名)為驗證基準坑(緯度57.3S，經度175.2E，最長直徑66 km).

表3 識別結果2

表3 中識別的坑參數以Grimaldi坑(1935年命名)為驗證基準坑(緯度5.5S，經度68.3W，最長直徑172 km).

表4 識別結果3

表4 中識別的坑參數以Ventris坑(1970年命名)為驗證基準坑(緯度4.9S，經度158.0E，最長直徑95 km).

將直徑大于60 km的524個知名坑作為基準驗證坑，坑中心點經度偏離最大15.63%，最小0.77%，坑中心點維度偏離最大13.16%，坑直徑長邊偏離最大4.39%.

基點彌散方法識別月隕坑的誤差與誤判來源于3個方面： ① 定義導致的誤差. 將月球坑定義為獨立洼地，根據擬和橢圓長短軸半徑與坑深的比值分辨坑的年齡和形態. ② 數學擬和算法導致的誤差. 坑參數(坑所在的基礎地形坡度、 坑數學中心、 坑深、 坑大小參數、 坑唇等細節特征)需要采用數學擬和為橢圓的假設，所建立模型與地勢復雜的坑存在誤差. ③ 機器實施算法導致的誤差. 月球真實數字高程模型分辨率在百米級，用于著陸安全性分析的月球數字高程模型分辨率達到厘米級，動態范圍大，PC機工作所涉及的采樣頻率范圍大.

3.2 統計分布規律

虹灣地區屬于平緩月海區域，NASA給出此類地形撞擊坑的直徑與大于該直徑的坑的數目之間的統計關系[13，14]

lgN=-2×lgD-1,D≤40,

lgN=-3×lgD+0.602,D>40,

(7)

式中：N為單位平方米的撞擊坑數目；D為撞擊坑臨界直徑.

統計圖如圖16 所示.

(a) NASA給出的坑統計分布規律

(b) 算法得到的坑統計分布規律圖16 坑分布規律Fig.16 Distribution of lunar craters

橫坐標為撞擊坑直徑，單位為km，縱坐標為相應直徑的坑的個數(N). 圖16(a)為NASA給出的坑統計分布規律. 特征點動態供給法對虹灣地區坑直徑與個數的統計分布規律如圖16(b) 所示.

4 月海地貌仿真重建

4.1 撞擊坑的模型

采用雙拋物線擬合并繞z軸旋轉來構造月球坑模型. 坑唇間的坑底部分采用式(8)計算

(8)

坑唇部分的模型采用式(9)計算

(9)

式中：x為撞擊坑模型上各點x值；y為撞擊坑模型上各點y值；D為坑直徑；d為坑唇寬度；h為坑深；H為坑唇高；a為拋物線系數，如圖17 所示.

圖17 撞擊坑的模型圖Fig.17 Model of craters

當撞擊坑的直徑小時，其截面形狀描述為

z(rc)=

當直徑大時，其截面形狀描述為

(11)

式中：u,p和?為經驗性常數參數.

4.2 地貌仿真生成

選定顯示區域范圍，生成區域趨勢面，根據上文得到的撞擊坑分布規律，石塊(坑深為負值)分布規律，設定石塊和撞擊坑的最小直徑，最大直徑.

根據撞擊坑統計分布規律計算在該區域中撞擊坑的個數，在選定區域內隨機生成這些坑的中心點坐標，調用典型撞擊坑的數學模型，分別算出每個撞擊坑高程值. 用戶輸入參數界面如圖18 所示.

最后，生成該區域的地勢高程，再調用生成撞擊坑函數和生成石塊函數，每個點的高程值進行疊加. 算出該區域的三維高程數據. 三維視圖如圖19 所示.

圖18 生成地形輸入項圖Fig.18 The input interface of terrain generation

圖19 月海地形仿真三維顯示圖Fig.19 Three dimensional display of lunar mare simulation terrain

圖19 中所有坑均為雙拋物線擬合坑，直徑分布遵循月海地區撞擊坑分布規律，著陸器4個足墊的平均坡度均為，是最優的著陸情況.

5 結 論

撞擊坑的直徑分布范圍很寬，小的只有幾十厘米或更小. 直徑大于10 km的撞擊坑的總面積約占整個月球表面積的7%-10%. 基于基點彌散法識別月球撞擊坑，進而進行地貌仿真重建. 得到結論：

1) 單個月球撞擊坑的自動識別，可給出坑中心點經緯度，坑深，坑擬合橢圓的長軸直徑，短軸直徑，坑唇線與當地水平面的傾斜角.

2) 全局月球坑統計分布規律. 可給出任意指定區域內，撞擊坑直徑和大于此直徑的坑的個數之間的統計關系. NASA給出的月球坑統計分布規律在月海區描述與分析結果基本一致，高地區月球坑數量比本文分析結果少10%左右.

3)不同年齡的月球坑用雙拋物線旋轉函數模擬，仿真重建了符合坑分布統計規律的月海地貌，可用于著陸安全概率分析和區域選擇.