一种顾及地形特征的山区LiDAR高精度DEM提取算法

彭检贵 邢元军 宋亚斌 王宗跃

摘要：针对目前常规DEM插值算法无法很好地表达断裂区域地形细节的问题，提出了一种基于地形特征约束的高精度DEM插值算法。首先利用离散点构建无约束TIN，然后利用边缘检测算法提取地形断裂特征，并嵌入地形断裂线作为约束，构建约束TIN，在此基础上制作高精度DEM。最后，利用黄土高坡数据进行了DEM提取实验。结果证明，该算法较传统算法在地形细节的保留上更具优势。

关键词：地形特征；数字高程模型；激光雷达；约束三角网

DOIDOI：10.11907/rjdk.171700

中图分类号：TP312文献标识码：A文章编号：16727800（2017）010002704

0引言

数字高程模型（Digital Elevation Model，简称DEM）是用一组有序数值阵列形式表示地面高程的一种实体地面模型[1]，是数字地形模型（Digital Terrain Model，简称DTM）的一个分支[2]。DEM是描述包括高程在内的各种地貌因子的基础，在测绘、水文、气象、地貌、地质、土壤、工程建设、通讯、军事等国民经济和国防建设领域有着广泛应用[3]。DEM数据源的获取方式多种多样，如：传统地面测量、摄影测量、干涉合成孔径雷达（InSAR）和机载激光雷达（Light Detection and Ranging，简称LiDAR）等[4]。地面测量精度高，但耗时又费力，不适合大规模应用；被动方式的摄影测量在植被覆盖区域和不利气象条件下效果不理想[5]；机载激光雷达和InSAR技术由于具有较高精度，并能进行较低成本的大范围应用，使其成为了工程应用的首选[6]。

机载激光雷达集激光扫描系统、全球定位系统（GPS）和惯性导航系统（INS）三种技术于一身，无需大量地面控制点，即能快速准确地获取地表高密度、高精度的高程信息[7]，已逐渐成为制作高精度DEM产品的一种优质数据源。从机载LiDAR点云数据到DEM产品，需要一定处理流程，其中最重要的步骤是点云滤波和DEM插值。滤波是指从点云中分离地面点、非地面点（包括人工目标（如建筑物、桥梁）和自然地物表面（如：灌木、树木等）的过程。作为制作DEM的重要步骤，滤波一直是国内外学者的研究重点。目前，已形成了多种代表性算法，如数学形态学的滤波方法[810]、基于分割后拓扑重建的滤波方法[11，12]、基于坡度的滤波算法[1314]、渐近三角网加密算法（PTD）[1517]等。通过ISPRS第三工作组对几种常用滤波算法的实验与分析[18]，Axelsson认为渐近三角网加密算法在地形保留和误差抑制上做的最好。而专门针对DEM插值的研究并不多，目前常用的插值算法有移动曲面拟合法、径向基函数法、反距离加权法、克里金插值法、线性三角网插值算法等[19]。

以上几种算法是在地形连续的假设上利用数学公式拟合地形，在点云密度较高的地区能取得较好效果，但在地形起伏较大、点云分布稀疏的区域则存在较多的信息损失。在植被密布的山区，LiDAR穿透力有限，经过滤波处理后，植被下方的地面点云密度可能非常低，此时传统插值方法已不再适用，需要顾及地形的骨架信息，尽可能地保留地形细节。

针对上述问题，本文提出一种顾及地形骨架的山区LiDAR高精度DEM提取方法。首先，利用普适性较好的渐近加密三角网算法（PTD）滤波获取全部地面脚点，同时获取无约束TIN；然后利用边缘检测算法提取地形骨架线（包括山脊线、山谷线以及断裂线），将地形骨架线作为约束条件，构建约束TIN，并基于约束TIN制作高精度DEM；最后，利用黄土高坡的点云进行了DEM提取实验。实验结果证明，本文提出的插值算法在地形起伏较大的山区能够获得精度更高、更加貼近真实地形的DEM。

1顾及地形特征的山区LiDAR高精度DEM提取

根据约束条件的有无，三角剖分可分成常规三角剖分（Delaunay Triangulation，简称DTIN）和约束三角剖分（Constrained delaunay Triangulation，简称CDTIN）。利用离散点云构建三角网时，不仅对三角形的形状有要求，而且离散数据本身也会影响三角网的局部合理性。通过对一些特殊地物或地形的点（如断裂、山脊、山谷、堤坝等）进行组合，对TIN进行强制性约束，从而使构建的三角网更符合实际地形。因此，基于约束三角网内插的DEM也更加贴近真实地形。本文在无约束TIN中嵌入地形结构线，构建约束TIN，最后基于约束TIN制作高精度DEM。

1.1基于PTD的地面点提取与非约束三角网构建

Axelsson提出的PTD算法以三角网为基础，通过由粗到细的过程，利用一定的光滑条件，逐步获取精密的地面点。由于该算法具有较好的普适性，已在商业化的LiDAR数据处理软件（如芬兰的Terra Solid、中国的LiDAR_Suite）中实现，算法的基本步骤如下：①对LiDAR点云构建格网索引，所需的参数一般需要人为设定；②对于格网的每个分块，搜索其最低点作为初始地面点，并构建稀疏的地形TIN；③根据待判断点与所在三角形的角度、距离、镜像点原则等条件作判断（见图1），确定余下LiDAR脚点是否能够加入三角网，若满足条件，则将其加入三角网中；④重复步骤③，直到所有点都已被判定为地面点或非地面点。

图1PTD算法加密过程

PTD作为经典的滤波算法之一，已得到了广泛研究，其具体细节在文献[15][17]中有着非常翔实的描述，本文不再赘述。本文选择PTD算法滤波获取地面脚点，首先考虑其极好的普适性，在地形起伏大的山区也能有效保留地形骨架并较好地抑制噪声；另一方面，PTD算法在加密滤波的同时也生成三角网，极大简化了后续插值工作。

1.2嵌入地形特征构建约束TIN

地形特征线，主要包括山脊线、山谷线以及断裂线。作为地形的骨架，地形特征线不仅决定了地形地貌的几何形态和基本走势，也有特定的物理意义。本文利用国产激光雷达数据处理软件LiDAR_Suite基于离散的地面点云提取断裂线、山脊线和山谷线，将上述3种地形特征线作为约束条件嵌入TIN模型，构建约束Delaunay三角网。endprint

约束边嵌入的算法很多，如分割合并法[20]、Shell三角化法[21]以及对角线交换算法[22]等。综合考虑实验难度以及插入效率，本文选择对角线交换算法构建约束TIN。对角线交换法的基本思想是，从起始点出发，判断每一条对角线的可交换性，若可以则交换，反之则继续迭代，直到约束边完整嵌入三角网为止。为了避免在影响域比较复杂时，判断和交换也随之变得复杂从而导致迭代失败的情况，本文采取m+2多边形对角线交换法，通过约束边将影响域分成2个区域，然后删除这两个多边形区域内的三角形，重新在该区域生成三角形，具体流程如下：

（1）假设有颜色区域为断裂区域（见图2），BL（浅灰色部分）表示高程较高的平坦区域，BR（深灰色部分）表示断裂线下方区域，L1L2为断裂线段（约束线段）。将约束线段的各个顶点看成离散点，将约束线段顶点加入到获得的CDT中。利用LOP算法对三角网进行调整，使之满足LiLj（i，j

（2）根据约束边确定影响域L。假设正在处理的约束线段为L1L2，则与L1L2相交的三角形构成的区域称为约束线段L1L2的影响域MT={T＼-1，T＼-2，…，T＼-n}，其中MT中三角形的外边组成了该影响域的影响多边形B={L1，a，b，c，……，L2}（图2中有颜色的区域）。该影响多边形有3个特点：①B是简单多边形，L1L2為B的一条对角线，将B分成BL和BR两部分，并且BL和BR也必须是简单多边形（见图2）；②BL和BR也可以进行Delaunay三角剖分；③ 对于任意Lk，假如Lk是距离L1L2最近的点（Lk≠L1，Lk≠L2），可确定LiLk∈B，LkLj∈B。删除L1L2上下两侧的多边形BL和BR内的三角形，重新进行Delaunay 三角剖分，从而将约束线段L1L2嵌入到三角网中。以BL为例，实现步骤如下：①建立一个空的堆栈stack，首先将L1L2的边放入stack中；②在stack中弹出一边，设为LcLiLj，遍历BL顶点，寻找离LiLj最近的点Lc（Lc≠Li，Lc≠Lj），形成三角形LcLiLj；③假如LcLi（或LcLj）是BL的一条边，则不入堆栈；反之则将LcLi（或LcLj）放入堆栈stack中；④若堆栈不为空，则返回步骤②，否则针对BL的操作结束；⑤ 使用LOP算法对BL和BR的剖分三角形进行局部优化处理，使之满足Delaunay三角形的两个基本性质。

如图2（a）所示，假如直接基于无约束三角网进行线性插值（具体细节见1.3节），由于三角形穿越断裂线，BL部分插值得到的高程将明显低于真实高程，而BR部分的高程则可能大于真实值。而在嵌入地形特征线构建约束三角网之后，地形结构线将不再穿越三角网（见图2（c）），线性插值获取的高程值更加贴近真实地形。

1.3基于约束TIN提取顾及地形特征的DEM

建立Delaunay三角网（TIN）之后，即可基于TIN计算区域内任意一点的高程，在点密度均匀的情况下，可以通过线性内插的方式获取各网点高程，即令三角形三点确定的倾斜平面作为该微小面元的地表。算法原理如图3所示。

对于待插值的格网点P（x，y），首先根据P点的平面坐标确定P隶属的分块格网号i，然后依次计算待插值点P与格网内其它点距离的平方：

D2i=（x-xi）2+（y-yi）2（2）

求取距离最小的点，设为Q1，然后依次取出Q1为顶点的三角形，判断点P是否位于该三角形内。若是，停止判断，反之则继续判断。若Q1为顶点的三角形都不包含点P，则继续判断距离次近的点Q2，如此迭代直到找到点P的外包三角形。

第三步：若P（x，y）的外包三角形为ΔQ1Q2Q3，三顶点的坐标分别是（x1，y1，z1）、（x2，y2，z2）和（x3，y3，z3），然后通过线性插值法确定P点的高程：

Z=Z1-（x-x1）（y21z31-y31z21）+（y-y1）（z21x31-z31x21）x21y31-x31y21（3）

2实验与结果分析

为了考察本算法的实际运行效果，本文基于国产激光雷达平台LiDAR_Suite进行二次开发，在VC++的环境下，实现了本文算法，并进行了实验与结果分析。实验数据选自黄土高坡地区，该区域地形复杂、地表形态丰富。数据采集于2011年9月，点云总数是1 562 831，点云平均间距约0.7m，垂直误差约0.15m。

分析图6可知，基于无约束TIN线性插值获取的DEM在地形平缓的区域，能够较好地保留地表细节，但在断裂区域附近，存在较明显的“平滑”现象。与此同时，如图7所示，基于约束TIN插值获取的DEM，其地表断裂处的地形细节都保留得更为完整。从视觉上分析，本文插值算法比传统不顾及地形骨架的方法更适用于山区DEM的插值。

为了更好地评价本文插值算法的效果，将基于GPS RTK量测获得的断裂地形附近的三维坐标作为检查点，定量评价基于移动曲面拟合、反距离加权法、线性三角网（无约束）、线性三角网（有约束）获取的DEM精度，其结果如表1所示。

表1中，A＼-1～A＼-6为地形较平缓区域所量测的检查点，A＼-7～A＼-11为断裂地形附近量测的检查点。分析表1可知，几种传统插值算法获取的DEM精度差异并不大，中误差在0.5～0.6m之间，约为点云高程误差的3～4倍，满足国家机载LiDAR数据后处理规范山区1：2 000 DEM的成图要求。另一方面，地形突变区域（A＼-7～A＼-11）的误差明显大于平坦地区（A＼-1～A＼-6）的误差，可见传统插值算法在地形细节保真方面有较大的提高空间。相对而言，基于约束TIN插值获取的DEM，在地形平缓区域的误差和传统插值结果接近，在地形突变结构附近区域的误差却明显小于传统插值结果，可见顾及地形骨架信息之后，DEM的提取精度有了明显提升。

3结语

本文针对传统DEM插值算法无法很好地表达地形起伏大区域地表细节的问题，提出了一种顾及地形骨架的高精度DEM插值算法。首先利用离散点构建无约束TIN，然后利用边缘检测算法提取地形骨架线，并嵌入地形断裂线作为约束，构建约束TIN，在此基础上利用约束TIN制作高精度DEM。实验结果证明，本文算法提取的山区DEM较传统算法精度更高，地形细节的保留更为逼真。

参考文献参考文献：

[1]汤国安，李发源，刘学军.数字高程模型教程[M].北京：科学出版社，2010.

[2]周启鸣，刘学军.数字地形分析[M].北京：科学出版社，2006.

[3]张瑞军，杨武年，刘汉湖，等.数字高程模型（DEM）的构建及其应用[J].工程勘察，2005（5）：6164.

[4]MAGUYA A S， JUNTTILA V， KAURANNE T. Adaptive algorithm for large scale dtm interpolation from lidar data for forestry applications in steep forested terrain[J]. Isprs Journal of Photogrammetry & Remote Sensing， 2013，85（11）：7483.

[5]高广，马洪超，张良，等.顾及地形断裂线的LiDAR点云滤波方法研究[J].武汉大学学报：信息科学版，2015，40（4）：474478.

[6]WILSON J P. Digital terrain modeling[J]. Geomorphology， 2012， 137（1）：107121.

[7]黃先锋，李卉，王潇，等.机载LiDAR数据滤波方法评述[J].测绘学报，2009，38（5）：466460.

[8]ZHANG K， CHEN SC， WHITMAN D，et al. A progressive morphological filter for removing nonground measurements from airborne LIDAR data[J]. Geoscience and Remote Sensing， IEEE Transactions on， 2003，41（4）：872882.

[9]CHEN Q， GONG P， BALDOCCHI D， et al. Filtering airborne laser scanning data with morphological methods[J]. Photogrammetric Engineering and Remote Sensing， 2007，73（2）：175185.

[10]PINGEL TJ， CLARKE KC， MCBRIDE WA. An improved simple morphological filter for the terrain classification of airborne LIDAR data[J]. ISPRS Journal of Photogrammetry and remote Sensing， 2013，77（1）：2130.

[11]SITHOLE G， VOOR GEODESIE NC. Segmentation and classification of airborne laser scanner data[M]. Nederlandse Commissie voor Geodesie， 2005.

[12]彭检贵，马洪超，高广，等.利用机载LiDAR点云数据提取城区道路[J].测绘通报，2012（9）：1619.

[13]SITHOLE G. Filtering of laser altimetry data using a slope adaptive filter[J]. International Archives of Photogrammetry Remote Sensing and Spatial Information Sciences， 2011，34（3/W4）：203210.

[14]SUSAKI J. Adaptive slope filtering of airborne LiDAR data in urban areas for digital terrain model （DTM） generation[J]. Remote Sensing， 2012，4（4）：18041819.

[15]AXELSSON P. DEM generation from laser scanner data using adaptive TIN models[J].International Archives of Photogrammetry and Remote Sensing， 2000，33（B4/1）：111118.

[16]左志权，张祖勋，张剑清.知识引导下的城区LiDAR点云高精度三角网渐进滤波方法[J].测绘学报，2012，42（2）：246251.

[17]隋立春，张熠斌，张硕，等.基于渐进三角网的机载LiDAR点云数据滤波[J].武汉大学学报：信息科学版，2011，36（10）：11591163.

[18]SITHOLE G， VOSSELMAN G. Experimental comparison of filter algorithms for bareearth extraction from airborne laser scanning point clouds[J]. ISPRS Journal of Photogrammetry and remote Sensing， 2004，59（12）：85101.

[19]张良，马洪超，邬建伟.联合机载LiDAR数据和潮汐数据自动提取潮位线[J].遥感学报，2012，16（2）：405416.

[20]LEE D T， B J SCHACHTER. Two algorithms for constructing a delaunay triangulation[J]. International Journal of Parallel Programming， 1980，9（3）：219242.

[21]PIEGL L A， RICHARD A M. Algorithm and data structure for triangulating multiply connected polygonal domains[J]. Computers & Graphics， 1993，17（93）：563574.

[22]刘学军，龚健雅.约束数据域的Delaunay三角剖分与修改算法[J].测绘学报，2001，30（1）：8288.

责任编辑（责任编辑：黄健）endprint