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

禹 健，郝秀芳，安永泉

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

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

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

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

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

1 DEM數(shù)據(jù)預(yù)處理

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

1.1 投影校正

投影變換是對非水平基準(zhǔn)地貌遙測數(shù)據(jù)必不可少的預(yù)處理工作.

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

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

圖1 參考基準(zhǔn)球坐標(biāo)變換Fig.1 The reference spherical coordinates transform

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

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

(1)

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

方位角參數(shù)的方程為

(2)

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

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

(3)

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

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

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

(4)

投影變換后進行經(jīng)典的雙線性內(nèi)插，采用坐標(biāo)點對齊結(jié)合邊界異常點剔除進行拼接.

1.2 坡面參數(shù)計算

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

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

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

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

(5)

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

2 月球坑識別

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

2.1 定義

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

圖4 月球地形影像數(shù)據(jù)視圖Fig.4 Image data view of lunar terrain

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

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

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

圖5 月球地形DEM數(shù)據(jù)視圖Fig.5 DEM data view of lunar terrain

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

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

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

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

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

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

2.2.1 坑唇邊緣點確定

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

(6)

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

圖7 坑唇點確定算法結(jié)構(gòu)框圖Fig.7 The algorithm structure of searching crater lip points

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

2.2.2 非規(guī)則邊界月球坑的檢測

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

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

2.2.3 已識認(rèn)坑的隱藏

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

1) 獨立坑的隱藏

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

2) 復(fù)合坑的隱藏

當(dāng)出現(xiàn)復(fù)合坑，坑套坑時，利用復(fù)合坑間的指向關(guān)系以及在坑檢測過程中所得到的坑矢量特征，將柵格操作與矢量操作結(jié)合起來進行處理[11].

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

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

表1 坑更新前后的鄰接關(guān)系表

圖9 復(fù)合坑及其指向關(guān)系Fig.9 Composite craters and their mutual pointing relation

圖10 一次合并后的復(fù)合坑及其指向關(guān)系Fig.10 Composite craters and their mutual pointing relationafter once merged

2.3 坑參數(shù)計算

坑參數(shù)包括坑中心點經(jīng)緯度坐標(biāo)(x,y)，因為帶中央峰的撞擊坑最低點并不是中心點，所以坑中心點根據(jù)坑長短直徑中點確定. 表征坑位置，坑深depth，坑長軸直徑a，坑短軸直徑b，坑唇傾斜率θ(指坑識認(rèn)出后坑唇擬合橢圓與當(dāng)?shù)卣狈较虻膬A斜角，如圖11 所示).

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

3 坑識別結(jié)果

對任意區(qū)域撞擊坑進行自動識別，并統(tǒng)計撞擊坑個數(shù). 盡管所得到的坑很不規(guī)則，但對于個數(shù)統(tǒng)計不造成影響.

3.1 識別結(jié)果

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

圖12 坑識別結(jié)果示意Fig.12 Result of Craters identification

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

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

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

圖15 坑組Ⅰ的參數(shù)識別結(jié)果Fig.15 Result parameters of pit group Ⅰ

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

表2～表4 給出部分識認(rèn)出的坑指標(biāo)，位置特征相近度，直徑特征相近在±3 km，推定以知名坑為識別結(jié)果驗證基準(zhǔn).

表2 識別結(jié)果1

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

表3 識別結(jié)果2

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

表4 識別結(jié)果3

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

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

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

3.2 統(tǒng)計分布規(guī)律

虹灣地區(qū)屬于平緩月海區(qū)域，NASA給出此類地形撞擊坑的直徑與大于該直徑的坑的數(shù)目之間的統(tǒng)計關(guān)系[13，14]

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

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

(7)

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

統(tǒng)計圖如圖16 所示.

(a) NASA給出的坑統(tǒng)計分布規(guī)律

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

橫坐標(biāo)為撞擊坑直徑，單位為km，縱坐標(biāo)為相應(yīng)直徑的坑的個數(shù)(N). 圖16(a)為NASA給出的坑統(tǒng)計分布規(guī)律. 特征點動態(tài)供給法對虹灣地區(qū)坑直徑與個數(shù)的統(tǒng)計分布規(guī)律如圖16(b) 所示.

4 月海地貌仿真重建

4.1 撞擊坑的模型

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

(8)

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

(9)

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

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

當(dāng)撞擊坑的直徑小時，其截面形狀描述為

z(rc)=

當(dāng)直徑大時，其截面形狀描述為

(11)

式中：u,p和?為經(jīng)驗性常數(shù)參數(shù).

4.2 地貌仿真生成

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

根據(jù)撞擊坑統(tǒng)計分布規(guī)律計算在該區(qū)域中撞擊坑的個數(shù)，在選定區(qū)域內(nèi)隨機生成這些坑的中心點坐標(biāo)，調(diào)用典型撞擊坑的數(shù)學(xué)模型，分別算出每個撞擊坑高程值. 用戶輸入?yún)?shù)界面如圖18 所示.

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

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

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

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

5 結(jié) 論

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

1) 單個月球撞擊坑的自動識別，可給出坑中心點經(jīng)緯度，坑深，坑擬合橢圓的長軸直徑，短軸直徑，坑唇線與當(dāng)?shù)厮矫娴膬A斜角.

2) 全局月球坑統(tǒng)計分布規(guī)律. 可給出任意指定區(qū)域內(nèi)，撞擊坑直徑和大于此直徑的坑的個數(shù)之間的統(tǒng)計關(guān)系. NASA給出的月球坑統(tǒng)計分布規(guī)律在月海區(qū)描述與分析結(jié)果基本一致，高地區(qū)月球坑數(shù)量比本文分析結(jié)果少10%左右.

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