一种基于变化检测的LiDAR点云DEM更新方法研究

王杰栋 俞涵

（浙江省第二测绘院，浙江 杭州 310012）

当前，测绘生产部门利用机载LiDAR 数据生产DEM 主要采用自动滤波加后期人工编辑的方式[1]。基于单一LiDAR 点云数据进行滤波的结果往往存在大量的错分和漏分点云现象，后期需要大量的人工编辑来提高DEM精度，大大增加了DEM的生产周期和生产成本。

考虑到许多数据管理部门已经存有测区DEM 的实际情况，在获得最新的点云数据后，可以利用这些数据直接检测地形的变化区域（变化区域一般远小于整个测区范围），实现DEM 的快速更新[2]。因此，本文提出一种基于变化检测的LiDAR 点云DEM 更新方法。

1 算法原理

1.1 预处理

预处理主要包括LiDAR 点云数据与已有DEM 坐标系统一和LiDAR 点云滤波。本文以已有DEM 所定义的坐标系为参考坐标系，将LiDAR 点云数据统一到参考坐标系中。利用商用软件Terrasolid 对点云数据进行初始滤波处理，获得初始地面点及DEM。

1.2 候选变化区域获取

针对预处理获得的DEM 数据及已有的DEM 数据，利用公式（1）获取高程差值模型（dDEM）。

公式（1）中，t1、t2表示前后两个不同时相。通过对dDEM模型进行判断，可以定位测区内的变化区域。初始的dDEM除了表示空间形态变化的大区域范围外，还存在少量细小的图斑。这些细小的图斑是因为点云自动滤波获取的数据存在错分及漏分现象，导致生成的DEM 在部分区域存在偏差。为了有效去除干扰区域，确定真正发生变化的区域，本文采用高程差阈值处理和形态学滤波算法对dDEM模型进行进一步处理。

1.3 伪变化区域去除

1.3.1 高差阈值处理

通过设定两期DEM 相对高差阈值，可以消除大部分伪变化区域。

公式（2）中，P 表示像素值，（c，r）代表像素点所在位置，T 为高程差阈值。

根据机载LiDAR 数据后处理技术规定，比例尺为1∶2000 的DEM 在平地和丘陵地的高程中误差分别为0.4m 和0.5m。本文采用0.4m 作为高程差阈值。

1.3.2 多尺度形态学开重建滤波处理

由于滤波结果存在少量突起或凹陷错分点、DEM插值误差等问题，进行高差阈值去噪后的dDEM（0.4）中依旧有一部分伪变化区域。针对这些伪变化区域，本文通过将dDEM（0.4）转化为二值图像，采用多尺度形态学开重建滤波算法进行去除。

形态学开重建方法是一种基于测地膨胀的形态学重建方法[3]。测地膨胀不需要直接选择特定的结构元素，而是基于标记影像和掩膜影像进行迭代计算。假定f 为标记影像，p 为掩膜影像（Df＝Dp），B 表示结构元素，则形态学膨胀可定义为：

基于上述形态学膨胀重建的定义，形态学开重建可定义为：

公式（4）中O 表示形态学开运算。形态学开重建是在开运算的基础上增加了重建过程。

由于伪变化区域的形态大小各不相同，如果采用单一的结构元素进行滤波，往往不能有效去除不同尺度的伪变化区域噪声。本文采用有限个不同大小的结构元素参与影像滤波运算，利用不同尺度结构元素能去除对应尺度噪声的特点，构建多尺度形态学滤波器。多尺度形态学滤波定义如下：

公式（5）中，nB 表示尺度为n 对应的结构元素。

1.4 区域滤波

利用已有的DEM数据辅助变化区域点云进行滤波，主要采用优化初始地面种子点的选取方式。其基本流程如下：

（1）对变化区域内点云数据建立规则格网索引，格网尺度根据经验设置为5m。

（2）设置高程差阈值 Hthreshold，判断每个点与已有DEM的高程差值，将满足阈值的点云记为待定地面点。

（3）遍历格网内的待定地面点，选择高程值最小的待定地面点作为初始地面点。

（4）从已有DEM 中提取断裂线区域，对落入断裂线区域内的格网进行分裂处理，获得小格网（小格网尺度为格网的一半），将小格网内高程值最小的待定地面点加入到初始地面点集中。

（5）利用初始地面点数据，采用TIN 加密滤波算法获得变化区域的地面点云及DEM。

在进行区域滤波之后，还需对变化区域进行少量人工交互编辑修正，从而获得精确的DEM。

2 试验与分析

2.1 试验数据说明

试验数据取自吉林省长春市。前一期DEM 数据，如图1（a）所示，由采集自2009 年5 月的LiDAR 点云数据生成，DEM 分辨率是1m；后一期点云数据，如图1（b）所示，采集自2010 年10 月，数据范围约为500×600 m2，平均点密度约为6.20 个/m2。2009 年，该地区很多地块处在开发初级阶段，2010 年，随着大量房地产项目的启动，地表变化剧烈，地面出现了大量土堆和较深的建筑物地基。图1（c）是利用商业软件Terrasolid 进行自动滤波处理后生成的概略DEM。

2.2 试验结果与分析

利用两期DEM 数据相减，可以获取该时间段的地形地貌变化情况，如图2（a）所示；图2（b）是对图2（a）进行0.4m 高程差阈值处理后的结果。

图1 试验数据

图2 dDEM

针对部分高程差值和真正变化区域类似的伪变化区域，采用多次度形态学开重建进行去除，如图3 所示，其中图3（a）是图3（b）的二值化结果，图3（b）是形态学运算后的结果。

图3 形态学运算结果

选取图3（b）中红色区域作为变化区域，采用区域滤波算法获得的变化区域DEM，如图4（a），并将早期DEM 中的变化区域替换为新生成的区域DEM，从而实现DEM 更新，如图4（b）。在地形变化区域采用基于DEM 辅助的滤波方法能有效提高滤波的准确度，减少错分和漏分点云现象，降低后期的人工编辑量，且地形变化区域范围远小于整个测区范围，导致大部分测区可以使用原有的DEM 数据，提高了DEM 的生产效率，降低了生产成本。

图4 DEM

3 结束语

本文提出了一种基于变化检测的LiDAR 点云DEM更新方法，该方法充分利用已有的DEM，通过比较分析两期DEM 数据提取地形变化区域，再对这些区域进行滤波和人工编辑生成区域DEM，实现DEM的快速更新。试验表明，该方法在两次量测期间测区地形没有发生重大变化的情况下，利用点云数据和已有DEM 生成新的DEM，能有效提高生产及更新DEM 的速度， 降低生产成本。