CN101950436A - 利用激光高度计数据制作数字高程模型的方法 - Google Patents

Abstract

本发明公开了一种利用激光高度计数据制作数字高程模型的方法，该方法包括步骤：利用激光高度计获取高程测量数据；对激光高度计的高程测量数据进行滤波；经过滤波，对剔除奇异点之后的激光测高数据进行插值处理，制作全月球规则格网的数字高程模型DEM。利用本发明的方法制作的月球DEM模型，由于使用了全月数据，并且通过使用滤波和插值过程，使得制作的DEM模型在数据精度和分辨率方面都有显著提高。

Description

利用激光高度计数据制作数字高程模型的方法

技术领域

本发明涉及数字地形图制作方法，特别是涉及一种利用激光高度计数据制作数字的方法。

背景技术

早期的航空航天激光高度计系统包括70年代装载在Apollo飞行器上的激光高度计，近些年发射的激光高度计还有SLA(Shuttle Laser Altimeter)和MOLA(Mars Orbital LaserAltimeter)，它们可以获取米级精度的地球和火星的地形数据。在月球地形测绘中，高精度的月表高程测量是通过激光测距技术获取的。目前只有美国Clementine、Apollo 15、Apollo 16、Apollo17、日本SELENE和我国的“嫦娥一号(CE-1)”卫星对月球表面进行过激光测距；只有SELENE和CE-1的激光探测数据覆盖了全月面。此外Archinal等人利用Clementine激光高度计和UVVIS图像数据制作了全月DEM模型，即ULCN2005DEM模型。ULCN2005全球DEM模型由ULCN2005月球控制网提供，日本SELENE的全球DEM模型数据目前还没有发布，当前部分结果发表在H.Araki等人的论文。

中国于2007年10月24日发射的环绕月球的“嫦娥一号(CE-1)”卫星上搭载的激光高度计(LAM)和CCD立体相机对月球进行了一年的探测，获取了月面地形数据。其中，激光高度计的科学探测任务是获取沿星下点轨迹方向卫星到月面的距离测量数据，用于解算该方向上月表高程测量数据，利用高程散点数据高效获取月球数字高程模型(DEM)，用于月表地形制图和研究。

但是，由于仪器系统噪声、月表地形起伏、卫星轨道与姿态测量等多个环节都存在误差，LAM原始测高数据中不可避免地含有误差。从原始数据到DEM处理过程中数据的误差会被传播和放大，导致测高数据中存在明显的高程奇异点，甚至出现严重的失真和扭曲现象。

现有技术在制月球DEM过程中，一方面使用的数据仅是月球正面数据，另一方面未考虑数据中存在的各种误差，使得制作的月球DEM不精确不全面。

发明内容

(一)要解决的技术问题

为了克服现有技术的上述缺陷，本发明提供一种利用激光高度计产品数据制作月球DEM的方法，其使用了全月数据，并且在制作DEM过程中考虑了数据误差情况。

本发明的利用激光高度计产品数据制作月球DEM的方法包括步骤：利用激光高度计获取高程测量数据；对激光高度计的高程测量数据进行滤波；对滤波后的数据进行检查；经过滤波，对剔除奇异点之后的激光测高数据进行插值处理，制作全月球规则格网的数字高程模型DEM。

优选地，高程测量数据是全月球数据。

优选地，对存在偶然误差的高程测量数据进行滤波。

优选地，滤波处理为先进行单轨滤波，再进行区域滤波。

优选地，在单轨滤波过程中，逐点统计各高程点与其前后几个数据点的高程差和标准差：

$\stackrel{\‾}{h}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{h}_i$

$\mathrm{\σ}=\sqrt{\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(h_i-\stackrel{\‾}{h})}^2}{N-1}}$

式中，

为高程平均值，σ为标准差，N为点集中高程点个数，h_i为高程点的月面高程测量值，如果被统计的高程点满足

则该高程点被当作奇异点剔除，其中M为高程阈值设定参数。

优选地，月海区域取M＝2，高地区域取M＝4。

优选地，在单轨滤波之后，将高程散点数据合并，再进行区域滤波，区域滤波包括：1)选取高程点集，采用最小二乘拟合的方法，确定这些点集决定的一个最小二乘平面h_i＝a₀+a₁x_i+a₂y_i，即初始拟合面，计算每个点的拟合高程值

并计算初始拟合面的高程标准差σ，根据高程标准差σ设定一阈值Mσ，计算每个点高程值h_i与拟合高程值

的差值，并根据差值与该阈值的关系来滤除奇异点；2)在步骤1)之后，利用最小二乘原理确定剔除后的点集决定的最小二乘二次曲面

计算拟合曲面的高程标准差σ，根据该高程标准差σ设定一阈值Mσ，如果某个点高程值与拟合高程值之差超过该阈值，则该点被滤除，其中a₀，a₁，a₂为最小二乘平面系数值，b₀，b₁，b₂，b₃，b₄，b₅为最小二乘曲面系数，(x_i，y_i)高程点的月面位置。

优选地，采用克里金法进行插值处理。

优选地，制作DEM时规则格网的大小为3km。

优选地，对全月制作DEM模型，根据地形特点将全月表面划分为南极4个区域，利用克里金法将高程数据分别插值成空间分辨率为3km的规则格网DEM模型，4个区域的数据经镶嵌后得到全月球空间分辨率为3km的DEM模型。

利用本发明的方法制作的月球DEM模型，由于使用了全月数据，并且通过使用滤波和插值过程，使得制作的DEM模型在数据精度和分辨率方面都有显著提高。

附图说明

图1为本发明利用激光高度计数据制作数字高程模型方法的一实施例的流程图；

图2为本发明滤波处理过程的流程图；

图3为滤波前后的地形晕渲图；

图4为各种插值方法地形晕渲图；

图5为插值方法地形剖面线对比图；

图6为不同测距数据覆盖率情况下月球3km格网DEM模型高程测量结果图；

图7a为利用CE-1卫星激光高度计探测数据制作的全月面DEM模型；

图7b为利用CE-1卫星激光高度计探测数据制作的月球两极地区DEM模型；

图8a为使用CE-1激光测高DEM数据，利用本发明的方法制作的风暴洋-雨海-澄海地区三维景观图；

图8b为使用ULCN2005激光测高DEM数据，利用本发明的方法制作的风暴洋-雨海-澄海地区三维景观图；

图9a为使用CE-1激光测高DEM数据，利用了本发明的方法制作的月表中低纬度区域地形晕渲图；

图9b为使用ULCN2005激光测高DEM数据，利用了本发明的方法制作的月表中低纬度区域地形晕渲图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

参照图1所示的流程图，本发明的利用激光高度计产品数据制作月球DEM的方法包括：

(1)获取高程测量数据。

搭载在探月卫星上的激光高度计是实现获取月球表面三维影像的一个重要载荷，并且该类激光高度计一方面可测量卫星到月球表面星下点的距离，另一方面可提供月球表面数字高程模型DEM径向高度数据。

激光高度计获取数据的基本原理是从卫星上发射一束大功率的窄激光脉冲到月球表面，并接收月球表面后向散射的激光信号，通过测量激光往返延迟时间来计算卫星到月表的距离，也就是激光高度计获得的高程数据。

本申请使用的高程数据优选地使用CE-1探月卫星携带的激光高度计所测量的高程数据，也就是测高数据，该数据已经公布于http://moon.bao.ac.cn。

(2)对激光高度计的原始测高数据进行滤波。

由于仪器系统噪声、月表地形起伏、卫星轨道与姿态测量等多个环节都存在误差，激光高度计LAM原始测高数据中不可避免地含有误差。从原始数据到DEM处理过程中数据的误差会被传播和放大，导致测高数据中存在明显的高程奇异点，甚至出现严重的失真和扭曲现象。因此，为了得到更为精确的DEM，去除测高数据中的高程奇异点是有必要的。

原始测高数据误差按照性质可分为偶然误差(也就是随机噪声)、系统误差和粗差(也就是错误值)。系统误差与激光高度计LAM硬件设备有关，对测高结果的影响具有规律性，在数据预处理阶段通过系统校正已经得到了纠正。粗差是由于仪器工作状态不正常而导致的错误的输出结果，粗差的存在导致DEM地形重建发生严重失真和扭曲。粗差数据在激光高度计LAM数据预处理过程中已作为无效数据被剔除掉。偶然误差是激光高度计在数据采集过程中由月表地形起伏、卫星平台不稳定等随机因素引起的，对观测结果的影响没有任何规律。

在本发明中，通过对高程数据进行滤波来处理偶然误差。

高程数据滤波的效果与地形的起伏程度直接相关。月球表面大体上可划分为月海和月陆(也称月球高地)两大主要的地理单元，其中月海地区撞击坑密度较低，地势整体上比较平坦，地形变化较缓慢；月陆地区撞击坑密度较高，地形变化幅度较大，尤其在陆海交界和大型撞击坑坑缘，地形起伏最大。为了有效地滤除高程数据中的奇异点数据，优选地将全月表面划分为南极(S60～S90°)、北极(N60～N90°)、中纬度月海(S70～N70°，W95～E95°)和中纬度月陆(S70～N70°，W180～W85°，E85～E180°)4个区域，采用不同的高程阈值分别进行了滤波处理。

滤波过程包括单轨滤波和区域滤波。

单轨滤波是采用地学统计法，逐点统计各高程点与其前后几个数据点的高程差和标准差：

$\stackrel{\‾}{h}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{h}_i---\left(1\right)$

$\mathrm{\σ}=\sqrt{\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(h_i-\stackrel{\‾}{h})}^2}{N-1}}---\left(2\right)$

式中，

为高程平均值，σ为标准差，N为点集中的高程点个数，N优选为5，h_i为高程点的月面高程值。如果足印点高程满足下式条件，该点将作为奇异点被剔除。

$|h_i-\stackrel{\‾}{h}|>\mathrm{M\σ}---\left(3\right)$

式中，M为高程阈值设定参数，在本发明中M是根据实验结果确定的滤波效率，对于月海区域，优选取M＝2，对于高地区域，优选取M＝4。

对测高数据进行单轨滤波后，将高程散点数据合并，再进行区域滤波。区域滤波采用移动曲面拟合法，分两步完成。

首先，以DEM格网点为中心，其一定半径范围内的高程点集，采用最小二乘拟合的方法，确定这些点集决定的一个最小二乘平面，即初始拟合面，计算每个点的拟合高程值，并计算初始拟合面的高程标准差，判断每个点高程值与拟合高程值的差值是否满足阈值的要求，以剔除高程奇异点；然后，利用最小二乘原理确定剔除后的点集决定的最小二乘二次曲面，计算拟合曲面的高程标准差。如果某个点高程值与拟合高程值之差超过给定的阈值，则该点被滤除。在上面两个步骤中，高程阈值的设定与单轨滤波原则一致，即月海区域以2倍高程标准差作为阈值，高地区域以4倍作为阈值。

计算公式如下：

h_i＝a₀+a₁x_i+a₂y_i(4)

$h_i=b_0+b_1{x}_i+b_2{y}_i+b_3{x}_i^2+b_4{x}_i{y}_i+b_5{y}_i^2---\left(5\right)$

$\stackrel{\‾}{h}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{h}_i---\left(6\right)$

$\mathrm{\σ}=\sqrt{\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(h_i-\stackrel{\‾}{h_i})}^2}{N-1}}---\left(7\right)$

式中，h_i为点集中高程点的月面高程值，a₀，a₁，a₂为最小二乘平面系数值，b₀，b₁，b₂，b₃，b₄，b₅为最小二乘曲面系数，(x_i，y_i)是高程点的月面位置，为根据式(6)得到的拟合高程值，也就是高程平均值，σ为高程标准差，N为点集中的高程点个数，优选为5。。如果足印点高程满足下式条件，该点将作为奇异点被剔除。

$|h_i-\stackrel{\‾}{h_i}|>\mathrm{M\σ}---\left(8\right)$

高程阈值的设定与单轨滤波原则一致，在本发明中是根据实验结果确定的滤波效率，月海区域取2倍高程标准差作为阈值，高地区域采用4倍高程标准差。

(3)对滤波后数据进行检查。

由于数据本身的中断和丢失等多种原因，输入数据会产生断点，断点与伪高程附近均有正确高程被滤除。于是在滤波过程中需要增加对数据的检查，通过与保留的高程数据插值对比，对滤除掉的高程数据进行确认，将其中误滤除的数据重新保留。对剔除伪高程点的前后效果对比如图3所示，左中右分别是直接用未剔除伪高程点的数据、Apollo地形数据和剔除伪高程点的数据制作的渲染地形，图中箭头标识出了原数据中的伪高程点。尽管存在少量系统偏差，但滤波后地形数据域Apollo地形间的整体一致性体现出伪高程点得到了正确滤除。这样，通过对步骤(2)中被滤除的数据进一步检查，也就是进一步判断其是否是奇异点，如果是，则将其删除，如果不是，则保留。

上述数据滤波和检查过程可参照图2所示的流程图。

(4)通过以上滤波后，对剔除奇异点之后的激光测高数据进行插值处理，制作全球规则格网的DEM模型。

用于DEM模型的插值算法主要有克里金(kriging)法、反距离权重、最小曲率、不规则三角网等，在本发明中，通过对各种插值算法的内符合精度、高程渲染图、地形剖面等进行了实验和比较，如图4，5所示，本发明优选地使用克里金插值法对全月球激光测高数据进行处理，克里金法的内符合精度可达到200m以内。为了分析“嫦娥一号”激光测高数据的成图能力，确定全月球DEM的空间分辨率，经过实验发现，不同格网大小情况下高程数据的测量误差和DEM的精度不同，结果表明，在网格尺寸大于3km的情况下，网格中存在真实测距数据的概率大于50％，DEM高程误差优于58.6m(1σ)，小于高度计的高程测量实际误差60m(1σ)，因此，本发明优选使用3km作为规则格网大小进行全球DEM制图，其能够满足高程测量的误差。

数据覆盖率对DEM空间分辨率具有一定影响，这种影响可通过实验进行检验，例如，选择南半球高纬度地区、赤道附近和北半球高纬度地区3个实验区，分别对30％、50％和100％等数据覆盖率下高程测量值的相关性和数值误差进行分析，对3个实验区的分析结果相似，图6给出了赤道附近一个实验区的试验结果。

如图6所示，纵坐标为以某覆盖率制作的DEM模型上所有格网点的高程值，横坐标为以100％覆盖率制作的DEM模型上所有格网点的高程值。本图分析区域位于月球赤道附近西经178.5°～西经157.5°，南纬10.5°～北纬10.5°范围，面积约367802km²，原始测高数据共49467个。本图仅显示了30％、50％和100％覆盖的实验结果。对于50％数据覆盖率的情况，采样数为25567个，相关系数为ρ＝0.997，相对于100％覆盖率的高程测量误差范围为-976m～654m，算术平均值为-10m，标准差为σ＝59m。对于30％数据覆盖的情况，采样数为16749个，相关系数为ρ＝0.984，相对于100％覆盖率的高程测量误差范围为-1737m～1603m，算术平均值为-50m，标准差为σ＝251m。

实验分析表明，30％数据覆盖率引起的高程误差达251m(1σ)，远远超过了激光高度计的高程测量误差(60m(1σ))，而50％数据覆盖率引起的DEM高程误差为58.6m(1σ)，小于高度计的高程测量误差，因此认为当激光足印点对月面区域数据覆盖率达到50％以上时，制作的DEM模型的精度能够满足高度计月表高程的测量精度。通过实验分析，在3km格网下全月球的数据覆盖率达到51.6％。因此CE-1的激光高度计数据能够制作空间分辨率为3km的全月DEM模型。

在制作全月DEM地形图时，我们仍然采用滤波时的策略，根据地形特点将全月表面划分为南极(S60～S90°)、北极(N60～N90°)、中纬度月海(S70～N70°，W95～E95°)和中纬度月陆(S70～N70°，W180～W85°，E85～E180°)4个区域，利用克里金法将高程数据分别插值成空间分辨率为3km的规则格网DEM模型，4个区域的数据经镶嵌后得到全月球空间分辨率为3km的DEM模型，结果如图7a和图7b所示。

图7a是根据本发明的方法，利用CE-1卫星激光高度计探测数据制作的全月面DEM模型。原始探测数据获取于2007年11月28日至2008年12月04日，共约912万个探测数据，覆盖西经180°到东经180°，南、北纬90°之间的范围；投影方式采用摩尔维德等面积伪圆柱投影，中央经线采用西经270°，左半部为月球背面，右半部为月球正面；月球坐标系采用均地极轴坐标系，高程参考面采用半径为1737.4km的月球正球体表面，参考原点为月球质量中心，空间分辨率为3km，全球最大高差为19.807km，最高点位于(-158.656°E，5.441°N，+10.629km)Korolev和Dirichlet-Jackson盆地中间的隆起构造上，最低点位于(-172.413°E，70.368°S，-9.178km)南极艾肯特盆地(South Pole-Aitken Basin)内。

图7b是根据本发明的方法，利用CE-1卫星激光高度计探测数据制作的月球两极地区DEM模型。左图为月球南极60°S～90°S地区的月表DEM模型，右图为月球北极60°N～90°N地区的月表DEM模型，投影方式采用极球面等角投影，以极点为中心，70°纬线为等变形线；月球坐标系、高程参考面和参考原点与图6a相同，空间分辨率为3km。

为比较CE-1和ULCN2005的DEM质量，我们制作了月球正面同一地区的三维景观图8a，8b和中低纬度地区地形晕渲图9a，9b。

图8a，b是使用CE-1和ULCN2005激光测高DEM数据制作的风暴洋-雨海-澄海地区三维景观图比较。(a)CE-1数据景观图；(b)ULCN2005月面控制网数据景观图。图中A表示朱拉山脉；B表示阿尔卑斯山脉；C表示高加索山脉；D表示亚平宁山脉；E：虹湾。

图9使用CE-1和ULCN2005激光测高DEM数据制作的月表中低纬度区域地形晕渲图。图9a是CE-1数据晕渲图，图9b是ULCN2005月面控制网数据晕渲图。

图8a，b和图9a，b的数据覆盖范围、坐标系统、高程基准、投影方式完全一致，可以看出，CE-1的DEM模型能够清楚地辨别月表的地貌单元及其细节，而ULCN2005DEM模型的地形分辨能力明显较低，虽然能够识别大于100km的主要地貌单元，但不能分辨地形细节。

因此，CE-1获得的激光高度计DEM模型，与已有的DEM模型相比，数据的精度和分辨率得到了显著提高。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (11)

1.一种利用激光高度计数据制作数字高程模型的方法，该方法包括步骤：

步骤1，利用激光高度计获取高程测量数据；

步骤2，对激光高度计的高程测量数据进行滤波；

步骤3，经过滤波，对剔除奇异点之后的激光测高数据进行插值处理，制作全月球规则格网的数字高程模型DEM。

2.根据权利要求1所述的方法，其特征在于，步骤1进一步包括所述高程测量数据是全月数据。

3.根据权利要求2所述的方法，其特征在于，步骤2中对存在偶然误差的高程测量数据进行滤波。

4.根据权利要求3所述的方法，其特征在于，步骤2中进行滤波为先进行单轨滤波，再进行区域滤波。

5.根据权利要求4所述的方法，其特征在于：单轨滤波过程中逐点统计各高程点与其前后几个数据点的高程差和标准差：

$\stackrel{\‾}{h}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{h}_i$

$\mathrm{\σ}=\sqrt{\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(h_i-\stackrel{\‾}{h})}^2}{N-1}}$

式中，

为高程平均值，σ为标准差，N为点集中高程点个数，h_i为高程点的月面高程测量值，如果被统计的高程点满足

则该高程点被当作奇异点剔除，其中M为高程阈值设定参数。

6.根据权利要求5所述的方法，其特征在于在确定高程阈值时，月海区域取M＝2，高地区域取M＝4

7.根据权利要求6所述的方法，其特征在于：在单轨滤波之后，将高程散点数据合并，再进行区域滤波，区域滤波包括：

1)选取高程点集，采用最小二乘拟合的方法，确定这些点集决定的一个最小二乘平面h_i＝a₀+a₁x_i+a₂y_i，即初始拟合面，计算每个点的拟合高程值

并计算初始拟合面的高程标准差

根据高程标准差σ设定一阈值Mσ，计算每个点高程值hi与拟合高程值的差值，并根据差值与该阈值的关系来滤除奇异点；

2)在步骤1)之后，利用最小二乘原理确定剔除后的点集决定的最小二乘二次曲面

计算拟合曲面的高程标准差σ，根据该高程标准差σ设定一阈值Mσ，如果某个点高程值与拟合高程值之差超过该阈值，则该点被滤除，

其中a₀，a₁，a₂为最小二乘平面系数值，b₀，b₁，b₂，b₃，b₄，b₅为最小二乘曲面系数，(x_i，y_i)高程点的月面位置。

8.根据权利要求1-7任一项所述的方法，其特征在于，步骤3中采用克里金法进行插值处理。

9.根据权利要求8所述的方法，其特征在于：步骤3中在制作DEM时规则格网的大小为3km。

10.根据权利要求9所述的方法，其特征在于：步骤3进一步包括：对全月制作DEM模型，根据地形特点将全月表面划分为南极4个区域，利用克里金法将高程数据分别插值成空间分辨率为3km的规则格网DEM模型，4个区域的数据经镶嵌后得到全月球空间分辨率为3km的DEM模型。

11.根据权利要求8所述的方法，其中步骤2之后进一步包括检查步骤，用于确认被滤除掉的高程数据中是否存在被误滤除的高程点，如果存在，则将其重新保留。

Images