数字高程模型的构建方法、装置、设备及存储介质与流程

1.本申请属于地理信息处理技术领域，尤其涉及一种数字高程模型的构建方法、装置、设备及计算机可读存储介质。


背景技术：

2.数字高程模型(digital elevation model，dem)，是用一组有序数值阵列形式表示地面高程的一种实体地面模型，通过有限的地形高程数据实现对地面地形的数字化模拟。随着数字中国的进程加快发展，数字高程模型已经被应用于各种领域。
3.然而，经本申请的发明人研究发现，现有方法构建的数字高程模型普遍存在与实际地形偏差较大的问题。


技术实现要素：

4.本申请实施例提供一种数字高程模型的构建方法、装置、设备及计算机可读存储介质，能够减小构建的数字高程模型与实际地形之间的偏差，提高构建的数字高程模型的准确性。
5.第一方面，本申请实施例提供一种数字高程模型的构建方法，方法包括：
6.获取目标区域的点云数据集，点云数据集包括目标区域中多个位置点的空间坐标，空间坐标包括位置点的平面坐标和高程；
7.对点云数据集进行闭运算，得到去除低异常值的第一点云数据集，低异常值为高程小于预设范围的最小边界值、且高程与最小边界值的差值大于预设第一阈值的第一位置点的空间坐标；
8.基于第一点云数据集，构建目标区域的数字高程模型。
9.第二方面，本申请实施例提供了一种数字高程模型的构建装置，装置包括：
10.获取模块，用于获取目标区域的点云数据集，点云数据集包括目标区域中多个位置点的空间坐标，空间坐标包括位置点的平面坐标和高程；
11.滤波模块，用于对点云数据集进行闭运算，得到去除低异常值的第一点云数据集，低异常值为高程小于预设范围的最小边界值、且高程与最小边界值的差值大于预设第一阈值的第一位置点的空间坐标；
12.构建模块，用于基于第一点云数据集，构建目标区域的数字高程模型。
13.第三方面，本申请实施例提供了一种电子设备，电子设备包括：处理器、存储器及存储在存储器上并可在处理器上运行的计算机程序，计算机程序被处理器执行时实现如第一方面提供的数字高程模型的构建方法的步骤。
14.第四方面，本申请实施例提供了一种计算机可读存储介质，计算机可读存储介质上存储计算机程序，计算机程序被处理器执行时实现如第一方面提供的数字高程模型的构建方法的步骤。
15.经本申请的发明人研究发现，导致构建的数字高程模型与实际地形偏差较大的主
要原因是：现有技术在利用点云数据集构建dem的过程中，低异常值(例如因故障导致的高程极低值)会被作为正常点云数据参与到数字高程模型的构建中，造成数字高程模型与实际地形偏差较大。有鉴于此，本申请实施例的数字高程模型的构建方法、装置、设备及计算机可读存储介质，获取目标区域的点云数据，点云数据包括目标区域中多个位置点的空间坐标，空间坐标包括位置点的平面坐标和高程；对点云数据进行闭运算，得到去除低异常值的第一点云数据，低异常值为高程小于预设第一阈值的第一位置点的空间坐标；基于第一点云数据，构建目标区域的数字高程模型。由于本申请实施例通过闭运算去除了点云数据集中的低异常值，故低异常值不会参与到数字高程模型的构建中，从而使得构建的数字高程模型与实际地形之间的偏差较小，进而提高构建的数字高程模型的准确性。
附图说明
16.为了更清楚地说明本申请实施例的技术方案，下面将对本申请实施例中所需要使用的附图作简单的介绍，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
17.图1为本申请一实施例提供的数字高程模型的构建方法的流程示意图；
18.图2为图1所示的数字高程模型的构建方法中s101的流程示意图；
19.图3示意性示出了激光雷达坐标系；
20.图4为图1所示的数字高程模型的构建方法中s102的流程示意图；
21.图5为本申请又一实施例提供的数字高程模型的构建方法的流程示意图；
22.图6为本申请又一实施例提供的数字高程模型的构建方法的流程示意图；
23.图7为本申请又一实施例提供的数字高程模型的构建方法的流程示意图；
24.图8是本申请一实施例提供的数字高程模型的构建装置的结构示意图；
25.图9是本申请一实施例提供的电子设备的结构示意图。
具体实施方式
26.下面将详细描述本申请的各个方面的特征和示例性实施例，为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及具体实施例，对本申请进行进一步详细描述。应理解，此处所描述的具体实施例仅意在解释本申请，而不是限定本申请。对于本领域技术人员来说，本申请可以在不需要这些具体细节中的一些细节的情况下实施。下面对实施例的描述仅仅是为了通过示出本申请的示例来提供对本申请更好的理解。
27.需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括
……”
限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。
28.在对本申请实施例进行说明之前，为了便于对本申请实施例理解，下面首先对于本申请实施例涉及到的技术名词进行解释。
29.数字高程模型(digital elevation model，dem)：dem是通过有限的地形高程数据实现对地面地形的数字化模拟(即地形表面形态的数字化表达)，它是用一组有序数值阵列形式表示地面高程的一种实体地面模型，是数字地形模型(digital terrain model，dtm)的一个分支，其它各种地形特征值均可由此派生。一般认为，dtm是描述包括高程在内的各种地貌因子，如坡度、坡向、坡度变化率等因子在内的线性和非线性组合的空间分布，而dem是零阶单纯的单项数字地貌模型，其他如坡度、坡向及坡度变化率等地貌特性可在dem的基础上派生。
30.激光雷达点云：在逆向工程中通过测量仪器得到的产品外观表面的点数据集合也称之为点云，通常使用三维坐标测量机所得到的点数量比较少，点与点的间距也比较大，叫稀疏点云；而使用三维激光扫描仪或照相式扫描仪得到的点云，点数量比较大并且比较密集，叫密集点云。当一束激光照射到物体表面时，所反射的激光会携带方位、距离等信息。若将激光束按照某种轨迹进行扫描，便会边扫描边记录到反射的激光点信息，由于扫描极为精细，则能够得到大量的激光点，因而就可形成激光雷达点云。
31.数学形态学(mathematical morphology)：数学形态学是一门建立在格论和拓扑学基础之上的图像分析学科，是数学形态学图像处理的基本理论。基本的运算包括：腐蚀和膨胀、开运算和闭运算、骨架抽取、极限腐蚀、击中击不中变换、形态学梯度、top
‑
hat变换、颗粒分析、流域变换等。
32.开运算：在数学形态学中，开运算被定义为先腐蚀后膨胀。
33.闭运算：在数学形态学中，闭运算被定义为先膨胀后腐蚀。
34.腐蚀：腐蚀(通常用表示)，是数学形态学中两个基本的算子之一(另一个是膨胀)。它最初是为二值图像定义的，之后被扩展到灰度图像以及随后的完全格。
35.膨胀：膨胀(通常用
⊕
表示)，是数学形态学中基本的算子之一。它最初是为二值图像定义的，之后被扩展到灰度图像以及随后的完全格。
36.德洛内(delaunay)三角网：点集的三角剖分(triangulation)对数值分析(比如有限元分析)以及图形学来说，是极为重要的一项预处理技术。尤其是delaunay三角剖分，由于其独特性，关于点集的很多种几何图都和delaunay三角剖分相关，如voronoi图，emst树，gabriel图等。delaunay三角剖分有最大化最小角，“最接近于规则化的“的三角网和唯一性(任意四点不能共圆)两个特点。
37.克里金法(kriging)：是依据协方差函数对随机过程/随机场进行空间建模和预测(插值)的回归算法。在特定的随机过程，例如固有平稳过程中，克里金法能够给出最优线性无偏估计(best linear unbiased prediction，blup)，因此在地统计学中也被称为空间最优无偏估计器(spatial blup)。
38.在阐述本申请实施例所提供的技术方案之前，为了便于对本申请实施例理解，本申请首先对现有技术中存在的问题进行具体说明：
39.数字高程模型dem的构建过程大致分为：1)获取点云数据集；2)对点云数据集进行滤波，滤除点云数据集中的非地面点；3)基于滤除非地面点的点云数据集生成dem。
40.现有技术在利用点云数据集构建dem的过程中，没有滤除低异常值，低异常值会被作为正常点云数据参与到数字高程模型的构建中，这会使得地面点滤波过程产生严重偏差，最终会造成数字高程模型与实际地形偏差较大。其中，低异常值可以理解为在高程上明
显低于点云数据集中正常点云数据的异常空间坐标。低异常值主要是因多次反射的脉冲或激光测距仪的故障而产生，低异常值会影响对邻近地面点的误判，若不滤除会造成数字高程模型与实际地形偏差较大。
41.为了解决现有技术问题，本申请实施例提供了一种数字高程模型的构建方法、装置、设备及计算机可读存储介质。
42.本申请实施例的技术构思在于：通过闭运算去除点云数据集中的低异常值，避免低异常值参与到数字高程模型的构建中，从而使得构建的数字高程模型与实际地形之间的偏差较小，进而提高构建的数字高程模型的准确性。
43.下面首先对本申请实施例所提供的数字高程模型的构建方法进行介绍。
44.图1为本申请实施例提供的数字高程模型的构建方法的流程示意图。如图1所示，该方法可以包括以下步骤：
45.s101、获取目标区域的点云数据集，点云数据集包括目标区域中多个位置点的空间坐标，空间坐标包括位置点的平面坐标和高程；
46.s102、对点云数据集进行闭运算，得到去除低异常值的第一点云数据集，低异常值为高程小于预设范围的最小边界值、且高程与最小边界值的差值大于预设第一阈值的第一位置点的空间坐标；
47.s103、基于第一点云数据集，构建目标区域的数字高程模型。
48.本申请实施例的数字高程模型的构建方法，首先获取目标区域的点云数据，点云数据包括目标区域中多个位置点的空间坐标，空间坐标包括位置点的平面坐标和高程；然后对点云数据进行闭运算，得到去除低异常值的第一点云数据，低异常值为高程小于预设第一阈值的第一位置点的空间坐标；最后基于第一点云数据，构建目标区域的数字高程模型。由于本申请实施例通过闭运算去除了点云数据集中的低异常值，故低异常值不会参与到数字高程模型的构建中，从而使得构建的数字高程模型与实际地形之间的偏差较小，进而提高构建的数字高程模型的准确性。
49.下面介绍上述各个步骤的具体实现方式。
50.首先介绍s101、获取目标区域的点云数据集。
51.在本申请实施例中，点云数据集可以包括目标区域中多个位置点的空间坐标，空间坐标可以包括位置点的平面坐标和高程，即三维坐标。目标区域可以是任意区域，本申请不限于此。
52.具体地，目标区域的点云数据集可以根据激光雷达的原始数据、激光雷达所在载体的轨迹pos文件以及激光雷达相对于载体的安装信息而获得。其中，激光雷达可以包括机载激光雷达(airborne light detection and ranging，lidar)。lidar是一种主动式航空遥感对地观测系统，它集成激光测距仪、全球定位系统gps和惯性测量单元imu于一身。lidar在三维空间信息的实时获取方面产生了重大突破，为获取高时空分辨率的地球空间信息提供了一种全新的技术手段。机载lidar通过发射和接收激光脉冲，能直接、快速获得地表密集的高精度三维点坐标(即机载lidar点云数据)。获取的点云数据精度高、密度大，能很好地表现地表细节。
53.在本申请实施例中，目标区域的点云数据集例如可以根据传统lidar系统采集的空间数据获得，也可以是根据全波形lidar系统采集的空间数据获得。与传统lidar系统相
比，全波形lidar系统具有两点优势：一方面接收的回波信号包含更加全面的地面几何信息。这意味在同一个光斑内，全波形lidar系统采集的数据比传统lidar系统采集的数据能更详细地描述目标地物的粗糙度、水平分布和垂直结构等信息(尤其在林区)；另一方面，全波形lidar系统采集的数据(下文简称全波形数据)经过波形分解以后可以得到更丰富的高密度、高精度三维点云数据。全波形数据分解的基本目的是实现从一维信号到三维点数据的转化，即根据一维波形信号获取三维地形信息。采用高斯分解方法可以提取额外波形信息(脉冲宽度和振幅等)，获得的波形参数不仅使波形数据处理精度更高，而且能反映地物特性，在一定程度上弥补了传统离散lidar数据缺少光谱信息的不足。另外，全波形相较于离散lidar数据具有分辨率方面的优势。通常，离散lidar数据由于缺少光谱信息，在地物精分类应用上存在先天不足，无法识别低矮植被，滤波方法难以将非地面点完全剔除。而全波形lidar数据在测量过程中能够详细记录地物的所有回波信息，在后续的波形分解中除了能获取高精度的几何数据还能提取其他的波形参数，比如地物的散射特性，使得点云分类过程大大降低对影像数据的依赖。这个优点还可以避免诸如多源数据融合的问题，甚至无需影像信息就能达到地物的高精度分类。
54.如图2所示，s101具体可以包括以下步骤s1011和s1012。
55.s1011、获取激光雷达的遥感器的投影中心的空间坐标、投影中心的姿态参数、投影中心与多个位置点之间的距离、激光雷达与成像扫描周期内的中间像元间的夹角。
56.如图3所示，g点表示激光雷达的遥感器的投影中心，p点表示目标区域中的任意一个位置点，(x
g
,y
g
,z
g
)表示激光雷达的遥感器的投影中心的空间坐标，(x
p
,y
p
,z
p
)表示激光雷达的遥感器的投影中心的空间坐标，s表示投影中心与位置点之间的距离，(ω,κ)表示投影中心的姿态参数，θ表示激光雷达与成像扫描周期内的中间像元间的夹角。
57.s1012、根据投影中心的空间坐标、投影中心的姿态参数、投影中心与多个位置点之间的距离、激光雷达与成像扫描周期内的中间像元间的夹角，确定多个位置点的空间坐标。
58.计算任意一个位置点的空间坐标的表述式为：
[0059][0060][0061][0062]
其中，b'＝cosωsinαcosκ+sinκsinω。
[0063]
通过s101可以获得目标区域中多个位置点的空间坐标，得到目标区域的点云数据集。
[0064]
以上为s101的具体实现方式，下面介绍s102的具体实现方式。
[0065]
s102、对点云数据集进行闭运算，得到去除低异常值的第一点云数据集。
[0066]
在本申请实施例中，低异常值为高程小于预设范围的最小边界值、且高程与最小边界值的差值大于预设第一阈值的第一位置点的空间坐标。如前所述，低异常值可以理解
为在高程上明显低于点云数据集中正常点云数据的异常空间坐标。预设范围可以看作是点云数据集中正常点云数据所处的范围，如[
‑
10米，20米]，低异常值为与预设范围的最小边界值(如
‑
10米)之间的差值大于预设第一阈值的空间坐标。其中，预设范围和预设第一阈值的大小具体根据实际情况灵活调整，本申请不限于此。
[0067]
闭运算包括膨胀和腐蚀两个过程。本申请实施例的思路是：先对点云数据集进行膨胀处理，当对低异常值进行膨胀运算时，低异常值的高程会通过膨胀以达到其邻域窗口定义的最高点；然后对点云数据集进行腐蚀处理，低异常点则恢复为原高程，将原高程与每个位置点的空间坐标的闭运算结果高程进行比较，将差值大于阈值的低异常值移除。
[0068]
具体地，在s102中，将膨胀和腐蚀这两种基本形态操作结合起来对低异常值进行去除。点云数据集f中待定空间坐标(x，y)的膨胀操作定义为：
[0069]
[δ
b
(f)](x,y)＝max{f(x+i,y+ji,j∈d
b
；(x+i),(y+j)∈d
f
)}
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0070]
其中，f表示点云数据集，b表示窗口，b在数学形态学中也称为结构元素；f和b所在区域分别为d
f
和d
b
。
[0071]
腐蚀是指从相邻点的高度获得最小值，表达式如下：
[0072]
[ε
b
(f)](x,y)＝min{f(x+i,y+ji,j∈d
b
；(x+i),(y+j)∈d
f
)}
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0073]
将膨胀和腐蚀相结合，将结构元素b在点云数据集f上的形态学闭运算定义为：
[0074][0075][0076]
其中，t1，t2和t3均为阈值，p
i
和p
j
均为位置点，z
i
和z
j
分别是位置点p
i
和位置点p
j
的高程值。由于极低高程和散射分布的低异常值特性，可以用上述表达式(7)中的准则来检测低异常值的集合l。如果位置点p
i
的闭运算结果大于阈值t1，则位置点p
i
相对于周围点具有局部极低的海拔。由于低异常点几乎没有近邻点，所以窗口内的相邻点之和被用来判断p
i
是否是散射的，其中，窗口是由实际数据中低离群值的分布预定义的。低异常值通常远低于周围点云，上述表达式(7)中的阈值可以通过试验和误差来确定。
[0077]
由此，通过对点云数据集进行闭运算，则可以将低异常值滤除。
[0078]
在形态学滤波中，滤波窗口(也称为数学形态学中的结构元素)，对滤波结果有着十分重要的作用。滤波窗口大小的选择对滤波性能非常敏感。一方面，只有在所采用的结构元素足够大的情况下才能过滤掉对象；另一方面，大型结构元素很可能会切断地形高度的突变，如悬崖、山脊和山峰。因此，为了处理不同尺寸的对象以及最大程度地滤除掉点云数据集中的低异常值，在s102中，可以迭代滤除低异常值，每次迭代中更新结构元素和结构元素对应的高差阈值，并利用更新的结构元素和高差阈值对点云数据集进行闭运算。
[0079]
具体地，如图4所示，s102可以包括如下步骤：
[0080]
s1021、利用结构元素对点云数据集进行膨胀处理，以改变每个位置点的高程，得到每个位置点改变后的高程。
[0081]
s1022、利用结构元素对膨胀处理后的点云数据集进行腐蚀处理，去除第一差值大于预设第二阈值的第一位置点的空间坐标，得到第二点云数据集，第一差值为位置点改变后的高程与改变前的高程之间的差值。
[0082]
其中，s1021和s1022的具体过程请参见上文描述，为了描述简洁，在此不再赘述。
[0083]
s1023、扩大并更新结构元素和预设第二阈值。结构元素与预设第二阈值一一对应，预设第二阈值又称高差阈值，随结构元素的改变而改变。在一些实施例中，例如可以采取指数递增的方式更新结构元素，如2，4，8，
……
，2
n
；也可以采取线性递增的方式更新结构元素，如1，2，3，
……
，n。两种方式相比，指数递增的方式所得的结构元素个数更少，计算速度更快；线性递增的方式所得窗口更多，所以后续的计算时间更长。
[0084]
s1024、判断是否满足预设第一条件，在满足预设第一条件的情况下，执行s1025；在不满足预设第一条件的情况下，返回s1021。其中，预设第一条件可以包括：更新结构元素的次数大于或等于预设次数，或者，结构元素的尺寸大于预设尺寸。简单来讲，预设第一条件可以包括已迭代的次数达到预设次数或者结构元素的尺寸已达到预设尺寸。
[0085]
s1025、将第二点云数据集确定为去除低异常值的第一点云数据集。
[0086]
通过迭代滤除低异常值，能够处理不同尺寸的对象以及最大程度上滤除掉点云数据集中的低异常值，使得最终构建的数字高程模型与实际地形之间的偏差较小，进而提高构建的数字高程模型的准确性。
[0087]
考虑到点云数据集的数据量较大，如果全部串行处理会造成计算紧张而造成的计算速度过慢，影响用户体验。因此，为了提高整体的运算效率，在对点云数据集进行闭运算之前，可以先对点云数据集进行分块处理，即，将点云数据集分成多个子点云数据集，然后分别对每个子点云数据集进行闭运算。与传统的分块方式不同的是，本申请实施例不是按照矩形或者某个规则的图形(即预设面积，如5米
×
5米)进行分块，而是按照点云的梯度变化进行分块，本申请采取的按照点云的梯度变化进行分块的方式不仅能够降低的整体的点云处理速度，也可以为后面的点云数据集处理(如点云数据分类和滤波)提供帮助。
[0088]
具体地，如图5所示，在s102之前，本申请实施例提供的数字高程模型的构建方法还可以包括以下步骤：
[0089]
s1011’、对于多个位置点中的每个位置点，根据位置点的空间坐标和至少一个初始位置点的空间坐标，计算位置点与至少一个初始位置点之间的坡度。在本申请实施例中，初始位置点为点云数据集中的任一空间坐标对应的位置点。
[0090]
具体地，首先可以从点云数据集的多个位置点中任选一个位置点作为初始位置点。然后，计算点云数据集中除该初始位置点之外的其他位置点与该初始位置点之间的坡度。
[0091]
s1012’、对于每个初始位置点，从点云数据集中筛选出坡度小于或等于预设第三阈值的位置点的空间坐标，得到每个初始位置点对应的子点云数据集。
[0092]
具体地，对于每个位置点，若该位置点与至少一个初始位置点之间的坡度小于或等于预设第三阈值，说明该位置点与初始位置点属于同一种类型，属于同一区块，则将该位置点的空间坐标归类至该初始位置点的空间坐标所在的子点云数据集中。而若满足预设数量的位置点与初始位置点之间的坡度均大于预设第三阈值，说明周围的位置点与初始位置点不属于同一分类，这时该次的分块结束，则更新初始位置点和子点云数据集。为了避免某个子点云数据集中的数据过多，即为了避免单块的点云数据过多，当子点云数据集中的空间坐标的数量大于预设数量时，也会更新初始位置点和子点云数据集。
[0093]
在更新初始位置点和子点云数据集之后，初始位置点和子点云数据集均已改变。此时，计算多个位置点中未归类的位置点与更新后的初始位置点的坡度，当坡度小于或等
于预设第三阈值时，将未归类的位置点的空间坐标归类至更新后的初始位置点的空间坐标所在的子点云数据集中，依次类推重复上述过程，直至将多个位置点的空间坐标均归类到对应的子点云数据集中，即直至将点云数据集的所有位置点的空间坐标分类至对应的子点云数据集中结束。
[0094]
相应地，在s102中，在对点云数据集进行闭运算时，可以对获得的多个子点云数据集进行闭运算，得到去除低异常值的第一点云数据集，从而提高了计算速度。
[0095]
在一些实施例中，为了使得最终得到的dem的准确性更高，在s102之前，还可以对点云数据集进行降噪，以去除点云数据集中的噪声点，这里的噪声点主要是高异常值。在本申请实施例中，高异常值为高程大于预设范围的最大边界值、且高程与最大边界值的差值大于预设第五阈值的空间坐标。高异常值可以理解为在高程上明显高于点云数据集中正常点云数据的异常空间坐标。预设范围可以看作是点云数据集中正常点云数据所处的范围，如[
‑
10米，20米]，高异常值为与预设范围的最大边界值(如
‑
10米)之间的差值大于预设第五阈值的空间坐标。其中，预设范围和预设第五阈值的大小具体根据实际情况灵活调整，本申请不限于此。
[0096]
具体地，如图6所示，在s102之前，本申请实施例提供的数字高程模型的构建方法还可以包括以下步骤：
[0097]
s1021’、根据每个位置点的空间坐标和多个目标位置点的空间坐标，计算每个位置点与多个目标位置点之间的距离的平均值。
[0098]
在本申请实施例中，每个目标位置点的空间坐标为点云数据中的任一空间坐标。比如，对于每个位置点而言，目标位置点为该位置点临近的位置点。
[0099]
这里，主要使用的方式是统计学滤波的方式。即，对于每个位置点而言，计算该位置点与多个目标位置点中的每个目标位置点之间的距离，得到多个距离；然后计算这多个距离的平均值。
[0100]
s1022’、去除点云数据集中平均值大于预设第四阈值的第二位置点的空间坐标，得到第三点云数据集。
[0101]
具体地，对于每个位置点而言，若该位置点与多个目标位置点之间的距离的平均值大于预设第四阈值，则说明该位置点的空间坐标为噪声点，则将该位置点的空间坐标去除。在去除点云数据集中平均值大于预设第四阈值的第二位置点的空间坐标之后，得到去除高异常值的第三点云数据集。
[0102]
相应地，在s102中，可以对第三点云数据集进行闭运算，得到第一点云数据集，从而使得最终得到的dem的准确性更高。
[0103]
为了进一步滤除高异常值，在s102之后，还可以对第一点云数据集进行开运算，进一步去除高异常值。
[0104]
具体地，在s102之后，本申请实施例提供的数字高程模型的构建方法还可以包括以下步骤：
[0105]
s103’、对第一点云数据集进行开运算，得到去除高异常值的第四点云数据集。
[0106]
开运算包括腐蚀和膨胀两个过程。首先通过腐蚀操作从第一点云数据集中除去比树木等结构元素的尺寸小的非地面点，然后通过膨胀操作恢复被腐蚀的建筑物边缘。开运算既可以移除比结构元素尺寸小的非地面点，又可以保留大于结构单元窗口尺寸的物体。
开运算可以过滤高异常值(植被、房屋等)。
[0107]
开运算的表达式如下：
[0108][0109]
腐蚀处理的表达式如下：
[0110]
(fθb)(i,j)＝z(i,j)＝min
z(s,t)∈w
(z(s,t))
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)
[0111]
膨胀处理的表达式如下：
[0112][0113]
其中，w表示结构元素的窗口，z(i,j)表示点云数据集中第i行第j列的位置点的高程。
[0114]
相应地，在s103中，可以基于进一步滤除高异常值的第四点云数据集，构建目标区域的数字高程模型，从而进一步提高构建的dem的准确性。
[0115]
需要说明的是，开运算也可以在s102之前执行，即开运算可以在闭运算之前执行，也可以是在闭运算之后执行，本申请不限于此。
[0116]
下面介绍s103的具体实现方式，s103、基于第一点云数据集，构建目标区域的数字高程模型。
[0117]
具体地，如图7所示，s103具体包括以下步骤：
[0118]
s1031、基于目标三角化算法对第一点云数据集进行三角化处理，得到包含多个三角网的目标三角网网络；
[0119]
s1032、获取目标区域的多个栅格，并确定每个栅格对应的三角网；
[0120]
s1033、利用目标内插算法和每个栅格对应的三角网的顶点的空间坐标，计算每个栅格内的目标点的高程，得到包含多个目标点的高程的数字高程模型。
[0121]
下面依次介绍s1031、s1032和s1033。
[0122]
构建dem模型的表达方式有很多种,而不规则三角网是其中最重要的方法之一，它是通过对离散的点数据构建许多个三角形来拟合真实的地面，它的特点是能够很好的表现地形起伏特点，特别是对山谷和山脊等特征线地区能够很好的表现，构建不规则三角网的算法和组织方式也有很多，在一些实施例中，本申请采用delaunay三角网对点云进行三角化。
[0123]
具体地，在s1031中，首先将第一点云数据集分割成多个子集，其中，每个子集包括预设数量(如每子集3或4个点)的位置点的空间坐标。然后对每个所述子集进行三角化，并利用局部优化算法lop对三角化的子集进行优化，得到每个子集对应的三角网。最后，将每个子集对应的三角网进行合并，得到目标三角网网络。
[0124]
具体过程如下：
[0125]
第一步、将第一点云数据集以横坐标为主、纵坐标为辅按升序排序。
[0126]
第二步、若第一点云数据集中点数大于阀值，则继续将第一点云数据集化为点个数近似相等的两个子集，并对每个子集做如下工作：
[0127]
①
获取每子集的凸壳；
[0128]
②
以凸壳为数据边界进行三角化，并用局部优化算法lop优化成delaunay三角网；
[0129]
③
找出连接左右子集两个凸壳的底线和顶线；
[0130]
④
由底线到顶线合并两个三角网。
[0131]
第三步、若第一点云数据集中点数不大于阀值，则直接输出三角剖分结果。
[0132]
在一些实施例中，数据点集划分可以采用递归分割快速排序法；子集凸壳的生成可以采用格雷厄姆算法；子集三角化可以采用任意方法，如子集最小到3或4个点则可直接三角剖分；子网合并则需先找出左右子集凸壳的底线和顶线，然后逐步合并三角剖分得到最终delaunay三角网。
[0133]
在s1032中，生成目标区域的多个栅格。栅格是根据用户的需求成图的精度形成的，如果需求的成图dem的精度是10米的，那么对应到的栅格就按照10米进行生成。比如，100
×
100的目标区域内，可以生成100个10
×
10的栅格。在生成目标区域的栅格之后，可以将生成的栅格映射到之前获得的目标区域的目标三角网网络之中，从而确定每个栅格对应的三角网。
[0134]
在s1033中，目标内插算法可以为逐点内插算法，逐点内插方法主要包括移动拟合法、加权平均法、voronoi图法、反距离加权法(inverse distance weighted，idw)等。在一些实施例中，s1033可以选用反距离加权法。反距离加权法插值类似于克里金插值，都是通过对已知采样点赋权得到未知点的值，不同之处在于克里金方法在考虑采样点间的距离同时还考虑采样点的整体空间分布情况，反距离加权插值仅是由未知采样点和已知采样点间的距离远近来决定权重，通过计算邻近区域离散点群的平均值来估计单元格点的方法为反距离加权插值，是一种逐点内插方法，其基本思想是当预测一个点的值时，接近预测点的值对预测点的影响更大，相反影响越小，即周围点对预测给出的权重与预测点的距离成反比。反距离加权插值作为一种简单有效的数据插值方法，其运算速度相对较快。
[0135]
通过反距离加权法和每个栅格对应的三角网的顶点的空间坐标，可以计算出每个栅格内的目标点(如栅格的中心点)的高程，进而可以得到包含多个目标点的高程的数字高程模型。
[0136]
利用反距离加权法计算每个栅格内的目标点的高程的表达式如下：
[0137][0138][0139]
其中，z0表示目标点的高程，(x0,y0)分别表示目标点的横坐标和纵坐标，(x
i
,y
i
,z
i
)表示三角网的任意一个顶点的空间坐标，p为任意正实数。
[0140]
由此，通过上述步骤可以得到数字高程模型dem。在得到数字高程模型dem，可以输出dem，dem可以按照标准格式tiff进行输出。
[0141]
基于上述实施例提供的数字高程模型的构建方法，相应地，本申请还提供了数字高程模型的构建装置的具体实现方式。请参见以下实施例。
[0142]
首先参见图8，本申请实施例提供的数字高程模型的构建装置800可以包括以下模块：
[0143]
获取模块801，用于获取目标区域的点云数据集，点云数据集包括目标区域中多个位置点的空间坐标，空间坐标包括位置点的平面坐标和高程；
[0144]
滤波模块802，用于对点云数据集进行闭运算，得到去除低异常值的第一点云数据
集，低异常值为高程小于预设范围的最小边界值、且高程与最小边界值的差值大于预设第一阈值的第一位置点的空间坐标；
[0145]
构建模块803，用于基于第一点云数据集，构建目标区域的数字高程模型。
[0146]
本申请实施例的数字高程模型的构建装置，获取模块用于获取目标区域的点云数据，点云数据包括目标区域中多个位置点的空间坐标，空间坐标包括位置点的平面坐标和高程；滤波模块用于对点云数据进行闭运算，得到去除低异常值的第一点云数据，低异常值为高程小于预设第一阈值的第一位置点的空间坐标；构建模块用于基于第一点云数据，构建目标区域的数字高程模型。由于本申请实施例通过闭运算去除了点云数据集中的低异常值，故低异常值不会参与到数字高程模型的构建中，从而使得构建的数字高程模型与实际地形之间的偏差较小，进而提高构建的数字高程模型的准确性。
[0147]
在一些实施例中，滤波模块802具体用于：利用结构元素对点云数据集进行膨胀处理，以改变每个位置点的高程，得到每个位置点改变后的高程；利用结构元素对膨胀处理后的点云数据集进行腐蚀处理，去除第一差值大于预设第二阈值的第一位置点的空间坐标，得到第二点云数据集，第一差值为位置点改变后的高程与改变前的高程之间的差值；扩大并更新结构元素和预设第二阈值；判断是否满足预设第一条件，在满足预设第一条件的情况下，将第二点云数据集确定为去除低异常值的第一点云数据集，预设第一条件包括：更新结构元素的次数大于或等于预设次数，或者，结构元素的尺寸大于预设尺寸；在不满足预设第一条件的情况下，返回执行利用结构元素对点云数据集进行膨胀处理，直至满足预设第一条件。
[0148]
在一些实施例中，本申请实施例提供的数字高程模型的构建装置800还可以包括：分块模块，用于对于多个位置点中的每个位置点，根据位置点的空间坐标和至少一个初始位置点的空间坐标，计算位置点与至少一个初始位置点之间的坡度，初始位置点为点云数据集中的任一空间坐标对应的位置点；对于每个初始位置点，从点云数据集中筛选出坡度小于或等于预设第三阈值的位置点的空间坐标，得到每个初始位置点对应的子点云数据集。滤波模块802具体用于：对获得的多个子点云数据集进行闭运算，得到去除低异常值的第一点云数据集。
[0149]
在一些实施例中，分块模块具体用于：确定一个初始位置点；当坡度小于或等于预设第三阈值时，将位置点的空间坐标归类至初始位置点的空间坐标所在的子点云数据集中；当满足预设第二条件时，更新初始位置点和子点云数据集，预设第二条件包括：满足预设数量的位置点的坡度均大于预设第三阈值，或者，子点云数据集中的空间坐标的数量大于预设数量；计算多个位置点中未归类的位置点与更新后的初始位置点的坡度，当坡度小于或等于预设第三阈值时，将未归类的位置点的空间坐标归类至更新后的初始位置点的空间坐标所在的子点云数据集中，直至将多个位置点的空间坐标均归类到对应的子点云数据集中。
[0150]
在一些实施例中，本申请实施例提供的数字高程模型的构建装置800还可以包括：降噪模块，用于根据每个位置点的空间坐标和多个目标位置点的空间坐标，计算每个位置点与多个目标位置点之间的距离的平均值，每个目标位置点的空间坐标为点云数据中的任一空间坐标；去除点云数据集中平均值大于预设第四阈值的第二位置点的空间坐标，得到第三点云数据集。滤波模块802具体用于：对第三点云数据集进行闭运算，得到第一点云数
据集。
[0151]
在一些实施例中，本申请实施例提供的数字高程模型的构建装置800还可以包括：开运算模块，对第一点云数据集进行开运算，得到去除高异常值的第四点云数据集，高异常值为高程大于预设范围的最大边界值、且高程与最大边界值的差值大于预设第五阈值的空间坐标。构建模块803具体用于：基于第四点云数据集，构建目标区域的数字高程模型。
[0152]
在一些实施例中，构建模块803具体用于：基于目标三角化算法对第一点云数据集进行三角化处理，得到包含多个三角网的目标三角网网络；获取目标区域的多个栅格，并确定每个栅格对应的三角网；利用目标内插算法和每个栅格对应的三角网的顶点的空间坐标，计算每个栅格内的目标点的高程，得到包含多个目标点的高程的数字高程模型。
[0153]
在一些实施例中，构建模块803具体用于：将第一点云数据集分割成多个子集，每个子集包括预设数量的位置点的空间坐标；对每个子集进行三角化，并利用局部优化算法对三角化的子集进行优化，得到每个子集对应的三角网；将每个子集对应的三角网进行合并，得到目标三角网网络。
[0154]
在一些实施例中，获取模块801具体用于：获取激光雷达的遥感器的投影中心的空间坐标、投影中心的姿态参数、投影中心与多个位置点之间的距离、激光雷达与成像扫描周期内的中间像元间的夹角；根据投影中心的空间坐标、投影中心的姿态参数、投影中心与多个位置点之间的距离、激光雷达与成像扫描周期内的中间像元间的夹角，确定多个位置点的空间坐标。
[0155]
图8所示装置中的各个模块/单元具有实现图1中各个步骤的功能，并能达到其相应的技术效果，为简洁描述，在此不再赘述。
[0156]
基于上述实施例提供的数字高程模型的构建方法，相应地，本申请还提供了电子设备的具体实现方式。请参见以下实施例。
[0157]
图9示出了本申请实施例提供的电子设备的硬件结构示意图。
[0158]
电子设备可以包括处理器901以及存储有计算机程序指令的存储器902。
[0159]
具体地，上述处理器901可以包括中央处理器(central processing unit，cpu)，或者特定集成电路(application specific integrated circuit，asic)，或者可以被配置成实施本申请实施例的一个或多个集成电路。
[0160]
存储器902可以包括用于数据或指令的大容量存储器。举例来说而非限制，存储器902可包括硬盘驱动器(hard disk drive，hdd)、软盘驱动器、闪存、光盘、磁光盘、磁带或通用串行总线(universal serial bus，usb)驱动器或者两个或更多个以上这些的组合。在一个实例中，存储器902可以包括可移除或不可移除(或固定)的介质，或者存储器902是非易失性固态存储器。存储器902可在综合网关容灾设备的内部或外部。
[0161]
在一个实例中，存储器902可以是只读存储器(read only memory，rom)。在一个实例中，该rom可以是掩模编程的rom、可编程rom(prom)、可擦除prom(eprom)、电可擦除prom(eeprom)、电可改写rom(earom)或闪存或者两个或更多个以上这些的组合。
[0162]
存储器902可以包括只读存储器(rom)，随机存取存储器(ram)，磁盘存储介质设备，光存储介质设备，闪存设备，电气、光学或其他物理/有形的存储器存储设备。因此，通常，存储器包括一个或多个编码有包括计算机可执行指令的软件的有形(非暂态)计算机可读存储介质(例如，存储器设备)，并且当该软件被执行(例如，由一个或多个处理器)时，其
可操作来执行参考根据本申请的一方面的方法所描述的操作。
[0163]
处理器901通过读取并执行存储器902中存储的计算机程序指令，以实现图1所示实施例中的方法/步骤s101至s103，并达到图1所示实例执行其方法/步骤达到的相应技术效果，为简洁描述在此不再赘述。
[0164]
在一个示例中，电子设备还可包括通信接口903和总线910。其中，如图9所示，处理器901、存储器902、通信接口903通过总线910连接并完成相互间的通信。
[0165]
通信接口903，主要用于实现本申请实施例中各模块、装置、单元和/或设备之间的通信。
[0166]
总线910包括硬件、软件或两者，将电子设备的部件彼此耦接在一起。举例来说而非限制，总线可包括加速图形端口(accelerated graphics port，agp)或其他图形总线、增强工业标准架构(extended industry standard architecture，eisa)总线、前端总线(front side bus，fsb)、超传输(hyper transport，ht)互连、工业标准架构(industry standard architecture，isa)总线、无限带宽互连、低引脚数(lpc)总线、存储器总线、微信道架构(mca)总线、外围组件互连(pci)总线、pci
‑
express(pci
‑
x)总线、串行高级技术附件(sata)总线、视频电子标准协会局部(vlb)总线或其他合适的总线或者两个或更多个以上这些的组合。在合适的情况下，总线910可包括一个或多个总线。尽管本申请实施例描述和示出了特定的总线，但本申请考虑任何合适的总线或互连。
[0167]
另外，结合上述实施例中的数字高程模型的构建方法，本申请实施例可提供一种计算机可读存储介质来实现。该计算机可读存储介质上存储有计算机程序指令；该计算机程序指令被处理器执行时实现上述实施例中的任意一种数字高程模型的构建方法。计算机可读存储介质的示例包括非暂态计算机可读存储介质，如电子电路、半导体存储器设备、rom、随机存取存储器、闪存、可擦除rom(erom)、软盘、cd
‑
rom、光盘、硬盘。
[0168]
需要明确的是，本申请并不局限于上文所描述并在图中示出的特定配置和处理。为了简明起见，这里省略了对已知方法的详细描述。在上述实施例中，描述和示出了若干具体的步骤作为示例。但是，本申请的方法过程并不限于所描述和示出的具体步骤，本领域的技术人员可以在领会本申请的精神后，作出各种改变、修改和添加，或者改变步骤之间的顺序。
[0169]
以上的结构框图中所示的功能块可以实现为硬件、软件、固件或者它们的组合。当以硬件方式实现时，其可以例如是电子电路、专用集成电路(application specific integrated circuit，asic)、适当的固件、插件、功能卡等等。当以软件方式实现时，本申请的元素是被用于执行所需任务的程序或者代码段。程序或者代码段可以存储在机器可读介质中，或者通过载波中携带的数据信号在传输介质或者通信链路上传送。“机器可读介质”可以包括能够存储或传输信息的任何介质。机器可读介质的例子包括电子电路、半导体存储器设备、rom、闪存、可擦除rom(erom)、软盘、cd
‑
rom、光盘、硬盘、光纤介质、射频(radio frequency，rf)链路，等等。代码段可以经由诸如因特网、内联网等的计算机网络被下载。
[0170]
还需要说明的是，本申请中提及的示例性实施例，基于一系列的步骤或者装置描述一些方法或系统。但是，本申请不局限于上述步骤的顺序，也就是说，可以按照实施例中提及的顺序执行步骤，也可以不同于实施例中的顺序，或者若干步骤同时执行。
[0171]
上面参考根据本申请的实施例的方法、装置(系统)和计算机程序产品的流程图
和/或框图描述了本申请的各方面。应当理解，流程图和/或框图中的每个方框以及流程图和/或框图中各方框的组合可以由计算机程序指令实现。这些计算机程序指令可被提供给通用计算机、专用计算机、或其它可编程数据处理装置的处理器，以产生一种机器，使得经由计算机或其它可编程数据处理装置的处理器执行的这些指令使能对流程图和/或框图的一个或多个方框中指定的功能/动作的实现。这种处理器可以是但不限于是通用处理器、专用处理器、特殊应用处理器或者现场可编程逻辑电路。还可理解，框图和/或流程图中的每个方框以及框图和/或流程图中的方框的组合，也可以由执行指定的功能或动作的专用硬件来实现，或可由专用硬件和计算机指令的组合来实现。
[0172]
以上，仅为本申请的具体实施方式，所属领域的技术人员可以清楚地了解到，为了描述的方便和简洁，上述描述的系统、模块和单元的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。应理解，本申请的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本申请揭露的技术范围内，可轻易想到各种等效的修改或替换，这些修改或替换都应涵盖在本申请的保护范围之内。