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

禹 健，郝秀芳，安永泉

(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)不同年龄的月球坑用双抛物线旋转函数模拟，仿真重建了符合坑分布统计规律的月海地貌，可用于着陆安全概率分析和区域选择.