1 引言

激光雷达(LiDAR)技术是一种主动遥感技术，与被动遥感技术相比，有着明显的优势，它可以快速、准确、高效地获取地球表面的三维结构信息^[1]。近年来，LiDAR技术已广泛应用于数字地面模型提取、三维建筑模型重建、道路检测、电力线提取、森林参数估计等地球空间信息学科的众多领域^[2-5]。

根据搭载平台的不同，LiDAR系统包括车载LiDAR、机载LiDAR、地面LiDAR以及背包LiDAR。相较于其他平台的LiDAR系统，背包LiDAR具有易操作、低成本、高精度、高效率、数据查看简洁方便的优点^[6-7]。背包LiDAR设计更加轻巧便捷，操作过程相对简单，且受外界的影响较小，无需考虑飞行航线布控，在无GPS等高成本系统条件下也实现厘米级精度，通过手机、平板电脑等移动端实时显示采集到的点云数据和运动轨迹^[8-9]。在森林环境中，背包LiDAR的优势最为明显。森林中崎岖不平的地形限制了传统平台的移动性，而背包LiDAR的轻巧灵活性使得复杂地区的森林数据采集有了新的突破，目前，背包LiDAR已广泛应用于数字地形模型提取、森林资源调查以及山路检测等领域^[10-12]。

点云滤波是背包LiDAR技术在以上领域应用的关键步骤。利用背包LiDAR系统获取的地表点云数据不仅包括地面点，还包含建筑物、桥梁、电线和树木等地物点^[13]。因此，需要对其进行滤波处理，以获得表征地表形态的地面点云。现有的滤波方法主要分为4类：1)基于插值的滤波算法；2)基于坡度的滤波算法；3)基于聚类分割的滤波算法；4)基于形态学的滤波算法。

基于插值的滤波算法是通过各种插值算法对选取的地面种子点进行插值计算，拟合出近似地形表面，再将点云数据的高程与其估计值进行比较，得到地面点云数据。常用的方法有薄板样条(TPS)插值法、最小二乘插值法、最近邻点插值法和不规则三角网插值法。Evans和Hudak^[14]采用多尺度曲率分类 (MCC) 方法将点云数据分类，并将超过曲率阈值的点判为地物点，否则判为地面点。Mongus等^[3]利用TPS插值法得到地形表面，并对点云进行逐级滤波。Chen等^[15]优化了TPS插值方法中种子点选择和滤波判断条件，提高了整体精度。但是，基于插值的滤波算法通常较耗时，对于陡峭和高度变化较大的地形区域的点云数据处理具有一定难度。Zhang等^[16]采用模拟布料滤波(CSF)方法来构建近似地面曲面，该方法先将地形翻转，再将一块柔软的虚拟布料根据重力作用落到地面上方，布料最终下落的形状为数字地面模型。

基于坡度的点云滤波方法最早由Vosselman^[17]提出，利用点云数据间的坡度区分地面点和非地面点，坡度阈值的设置是影响滤波效果的关键因素。Meng等^[18]提出一种多方向的地面滤波(MGF)方法，结合基于邻域和方向扫描方法的优点，检测不同方向上细微的高程变化，但需要设置坡度阈值和高差阈值。Sithole^[19]提出一种改进的方法，根据实际地形的坡度变化设置自适应坡度阈值。Susaki^[20]通过构建多层次的数字地面模型(DTM)计算变化的坡度阈值。虽然现有的一些基于坡度的算法可以根据实际地形动态调整阈值，但是基于坡度的算法不适用于植被较低的崎岖地形。

基于分割的点云滤波方法的原理是依据点云的特征进行聚类，然后对聚类分割的结果按照一定的规则进行点云滤波。Yang等^[21]根据分割结果的边界点与DEM的高差对分割结果进行滤波。Chen等^[22]采用基于区域增长的分割方法获取地面种子点，提出一种改进的基于分割的多分辨率分级分类滤波算法。综上所述，基于分割的滤波算法的识别效果在很大程度上取决于分割的质量。LAStools开源软件所使用的基于不规则三角网的DTM建立方法最早由Axelsson^[23]提出，该算法首先分割出多个点云区域，再从中依次选取初始地面点，构建不规则三角网(TIN)作为参考面，最后依据距离和夹角的阈值划分出地面点和地物点。

基于形态学的滤波算法最早由Lindenberger^[24]提出，其基本思想是对图像进行腐蚀与膨胀，通过处理前、后图像的特征形状尺寸变化来识别地物点。Zhang等^[25]提出一种滤波窗口逐渐增大时阈值不断发生变化的渐进形态滤波方法。Hui等^[26]将形态学运算和克里格插值方法结合，计算各层次局部区域的地形起伏度，并对滤波准则进行优化。但是，基于形态的滤波算法在复杂场景下很难设置合适的窗口尺寸，对于地形起伏较大的地区很容易误分类，故基于形态学的滤波算法在地形起伏较大地区的鲁棒性仍需提高。

综上所述，虽然已有的滤波方法能够在平坦区域获得较好的滤波效果，但点云滤波研究仍然存在以下难点与挑战：1)对于背包LiDAR点云滤波研究较少，现有方法无法在复杂山区获得良好的滤波效果；2)滤波方法往往需要复杂的参数调节才能获得较高的滤波精度，滤波方法自动化程度较低；3)滤波方法无法有效地保护地形细节，识别精度较低。针对以上问题，本文提出一种基于背包LiDAR点云多尺度渐进DTM (MPD)构建方法。该方法首先通过构建多尺度渐进滤波窗口，实现地面种子点的渐进提取。针对各尺度所获得的地面种子点，采用最小二乘方法构建出拟合DTM，再采用形态学开运算探测地物点，最后根据自适应变化的阈值筛选出地面点，渐进构建出准确的DTM。实验结果表明，MPD方法能够有效地实现背包LiDAR技术在复杂地形区域的高精度DTM提取。

2 基本原理

基于背包LiDAR点云多尺度渐进DTM构建方法的流程如图1所示。首先，根据研究区域内最大地物尺寸设置最大窗口尺寸，选取该窗口尺寸下的最低点作为地面种子点，并用最小二乘法(LS)构建地面曲面模型D_{DTM_LS}；然后，对由当前尺度下地面点插值计算方法生成的数字地面模型D_DTM与D_{DTM_LS}的差值进行形态学开运算，利用随窗口自适应变化的高差阈值T_H获取地物点；最后，利用去除地物点的地面点重新构建D_DTM，并设置随地形起伏自适应变化的高差阈值S_th，再对点云数据的实际高程值与D_DTM中对应的拟合高程值进行差值判断，更新地面点。滤波窗口由大到小变化，直至窗口尺寸达到最小时滤波结束，输出地面点。具体包括以下三个步骤：1)渐进DTM构建；2)多尺度形态学运算地物点探测；3)自适应地面点提取。

图 1. 多尺度渐进DTM构建方法的流程

Fig. 1. Flow of multi-scale progressive DTM construction method

下载图片 查看所有图片

2.3 渐进DTM构建

利用地面种子点建立数字地面模型DTM，其中种子点选取的准确性对于地面模型的建立有重要影响。将滤波窗口内最低点作为地面种子点，图2所示的圆点为筛选出的地面种子点。当窗口尺寸渐进变化时，利用所提的MPD算法获取地面种子点，并利用种子点拟合出一个粗糙的曲面模型。曲面拟合公式为 φ(x，y)＝f(x，y)＋JLS，(1)式中：φ(x，y)为构建的拟合曲面；f(x，y)为真实曲面；J_LS为误差矩阵。

图 2. 多尺度渐进拟合曲面。(a)初始拟合曲面；(b)渐进拟合曲面；(c)最终拟合曲面

Fig. 2. Multi-scale progressive fitting surfaces. (a) Initial fitting surface; (b) progressive fitting surface; (c) final fitting surface

下载图片 查看所有图片

采用拟合方法中经典的最小二乘法来拟合地面曲面φ(x，y)，将各窗口尺寸下获取的地面种子点代入曲面拟合公式，可得到各点(x_i，y_i)的近似高程值Z(x_i，y_i)，i＝1，…，n，n为地面种子点数量。当误差∑||φ(x，y)-f(x，y)||²为最小值时，拟合的地面曲面模型与实际地面模型最为接近。通过(2)式和(3)式可求出构建曲面的参数A。 minφ∑i＝1nJi2＝∑i＝1n||φ(xi，yi)－f(xi，yi)||2，(2)φ(x，y)＝XTA，(3)式中：J为拟合曲面与真实曲面间的误差矩阵，当误差平方和最小时，可获得最佳拟合曲面；X为点云坐标；A为待求曲面参数。

在各窗口尺寸下采用最小二乘法构建曲面模型，通过渐进式(由大到小)变化窗口尺寸，逐步拟合出准确的曲面模型。所提出的渐进DTM构建方法的滤波精度是由粗到精变化的。从图2可以看出，随着滤波窗口数量的逐渐减少，所提取的地面种子数量越多，所构建的拟合DTM曲面精度就越高。

2.4 多尺度形态学运算地物点探测

通过上述最小二乘法可得到各窗口尺寸下拟合的地面曲面模型D_{DTM_LS}以及各窗口下的非地物点插值生成的数字地面模型D_DTM。由于地形情况复杂，为降低地形坡度变化的影响，对D_DTM与D_{DTM_LS}的差值δ_DTM进行形态学开运算，将形态学开运算前、后高程变化较大的点视为地物点。其中，形态学开运算为先进行腐蚀运算，再对其结果进行膨胀运算^[27]。依据点云数据形态学开运算前、后的高差变化大小筛选出地物点。形态学开运算可表示为 (f∘g)＝(f⊖g)⊕g，(4)式中：º为开运算符号；⊖为腐蚀运算符号；⊕为膨胀运算符号。

在本实验中，最大滤波窗口需要大于实验区域内的最大建筑物尺寸，随后逐步减小滤波窗口尺寸，对背包点云数据进行滤波。这是因为若窗口尺寸过小，则不能剔除大型建筑物等地物点；若窗口尺寸过大，则会造成部分地面点被剔除，地形过平滑。为解决形态学开运算时窗口尺寸的设置问题，采用多尺度的滤波窗口来剔除不同大小的树木及建筑物等地物点。在大窗口尺寸下，可将突变的大型地物剔除，之后随窗口尺寸的减小，逐步剔除树木等小尺寸的地物点。

将各窗口尺寸下拟合的D_DTM与D_{DTM_LS}作差生成δ_DTM，通过对δ_DTM进行形态学开运算可生成δ_{DTM_open}，将δ_DTM与形态学开运算结果δ_{DTM_open}的差值和高差阈值进行比较。若形态学开运算前后地形表面高差大于高差阈值，判定为地物点；否则，判为地面点。在滤除地物点后，通过对空格网进行邻近点插值来保留地形。考虑到固定的高差阈值势必会造成一定的误差，本实验设定的高差阈值随窗口尺寸变化而自适应变化，在保留突变地形的同时，剔除地物点。计算公式为 δDTM＝DDTM－DDTM_LS，(5)TH＝η·w＋t，(6)式中：D_DTM为当前尺度下的DTM；D_{DTM_LS}为最小二乘法拟合的DTM；T_H为设置的高差阈值；w为当前形态学滤波窗口尺寸；t为最小地物识别阈值，在本实验中设置为0.2 m；η为滤波窗口尺寸的系数。为实现滤波阈值的自适应变化，本文将滤波阈值设置为窗口(w)的正相关函数，这是因为滤波窗口越大，表明要滤除的地物也越大，因此可以用较大的滤波阈值进行地物滤除。实验表明，η设置为0.1时，可实现地物的有效滤除。例如，当滤波窗口为20个格网尺寸时，滤波阈值则大于2 m，此阈值可以实现大型地物的有效滤除。最小地物识别阈值表征滤波窗口趋于0时的最小地物探测能力。由于本实验在此环节旨在尽可能地滤除地物，因此，凡是高程变化大于0.2 m的点都被识别为地物，并被剔除。

2.5 自适应地面点提取

在点云滤波过程中，固定的阈值不利于保护地形细节。为减少地形坡度变化对滤波算法的影响，实现滤波阈值的自适应变化，首先依据各窗口尺寸下生成的数字地面模型D_DTM，插值计算各点在不同地形坡度下的局部高程值，通过计算相邻格网间x方向和y方向的高程变化率，并依此设置各点随地形起伏变化的高差阈值。然后，计算每个点云数据的实际高程值Z和拟合高程值Z_DTM的差值。最后，依据判别公式，将差值结果与各点的高差阈值进行比较，提取地面点云数据。阈值S_th计算公式为 {Sx(x，y)＝[DDTM(x＋1，y)－DDTM(x－1，y)]/3Sy(x，y)＝[DDTM(x，y＋1)－DDTM(x，y－1)]/3,(7)Sth＝Sx(x，y)2＋Sy(x，y)2＋ε，(8)式中：S_x(x，y)、S_y(x，y)分别为格网(x，y)处的水平方向和垂直方向的高程变化率；S_th为随地形变化的高差阈值；ε为常数，表示地形坡度变化为0时的最小地面点云识别能力。在本文中，ε设置为0.6 m，即将高程变化小于0.6 m的点视为地面点。实验结果表明，此阈值能够实现地面点与地物点的分离。

若点云数据的实际高程值与拟合高程值大于高差阈值S_th，则将该点分类为地物点，标签为1；反之，将该点分类为地面点，标签为0。随着窗口尺寸的渐进变化，各窗口尺寸下生成的D_DTM和各点的高差阈值发生相应变化，对点云数据的类别标签进行更新。当窗口尺寸达到最小时，可得到较为准确的地面点标签。标签公式为 Llabel＝{1，Z(i)－ZDTM(i)>Sth(i)0，Z(i)－ZDTM(i)≤Sth(i)，(9)式中：Z(i)为第i个点的实际高程值；Z_DTM(i)为第i个点的拟合高程值；S_th(i)为第i个点对应的高差阈值；L_label为点云滤波的分类标签。

3 实验分析

3.1 实验数据

为验证所提MPD滤波算法在背包点云数据提取中的可行性，选取安徽省黄山市4组地形较为复杂的背包点云数据进行测试。图3所示为本实验采集的背包点云数据所处的地形场景，包括斜坡植被、突出地形、平坦公路等复杂因素。通过人工分类，将样本划分为地面点和非地面点，以确保类别标签的准确性。图4(a)～(d)分别为4组数据按RGB信息显示的点云数据。本实验的运行平台为Windows 10 64位操作系统，处理器i5-1035G1 CPU，1.19 GHz，内存为16 GB，算法采用MATLAB R2018a软件实现。

图 3. 研究区域场景。(a)场景1；(b)场景2

Fig. 3. Research area scenarios. (a) Scene 1; (b) scene 2

下载图片 查看所有图片

图 4. 原始点云数据。(a)样本1示意图；(b)样本2示意图；(c)样本3示意图；(d)样本4示意图

Fig. 4. Initial point cloud data. (a) Schematic of sample 1; (b) schematic of sample 2; (c) schematic of sample 3; (d) schematic of sample 4

下载图片 查看所有图片

3.2 评价指标

将点云分类精度评价指标分为I类误差(E_Type1)、Ⅱ类误差(E_Type2)、总误差(E_Total)和Kappa系数(I_Kappa)，其中：I类误差表示地面点被错分为非地面点的比例，即拒真误差；Ⅱ类误差表示非地面点被错分为地面点的比例，即纳伪误差。I类误差和Ⅱ类误差能够反映出算法的适用性。总误差为总体误差比例，能够反映算法的可行性；Kappa系数是统计学中度量一致性的指标，表征算法的鲁棒性。4个精度评价指标的计算公式为 EType1＝ba＋b，(10)EType2＝cc＋d，(11)ETotal＝b＋ca＋b＋c＋d，(12){IKappa＝p0－pc1－pcp0＝(a＋d)/epc＝(a＋b)(a＋c)＋(c＋d)(b＋d)/e2e＝a＋b＋c＋d,(13)式中：a为正确分类为地面点的数量；b为错误分类的地面点的数量；c为错误分类的地物点的数量；d为正确分类为地物点的数量。

3.3 实验结果

MPD滤波误差结果如表1所示。在4组样本中，MPD滤波算法的I类误差和II类误差的平均值分别为13.23%、6.11%，总误差约为7.88%，Kappa系数的平均值为79.48%。

3.4 对比分析

为进一步验证MPD滤波算法的有效性，选取基于渐进TIN(LAStools)、多尺度曲率分类(MCC)和布料模拟滤波(CSF)的算法分别对相同的背包点云数据进行滤波，并进行精度评价。

表2为MPD滤波算法与LAStools、MCC、CSF 3种算法的滤波结果。在4组复杂地形的背包点云数据中，MPD滤波算法的I类误差和II类误差的平均值分别为13.23%和6.11%。LAStools、MCC、CSF 3种滤波算法在4组数据中I类误差平均值分别为63.17%、51.31%和59.37%，均远大于MPD滤波算法的I类误差平均值。这是因为MPD滤波算法可根据地形坡度变化自动调整阈值，受地形坡度影响较小。MPD滤波算法的II类误差平均值大于其余3种滤波算法的II类误差平均值，这是因为LAStools、MCC、CSF 3种算法的I类误差过大，将大部分地面点错误分为地物点，从而导致纳伪误差较小，因此这3种滤波算法的II类误差较小。

图5(a)为4种滤波算法的总误差分布直方图，可以看到，相较于LAStools、MCC、CSF 3种滤波算法，MPD滤波算法的总误差最低。从Kappa系数分布直方图[图5(b)]可以看出，上述3种滤波算法的Kappa系数远低于MPD滤波算法的Kappa系数。LAStools、MCC、CSF 3种滤波算法在4组样本中的总误差平均值分别为16.03%、13.67%、15.72%，而MPD滤波算法的总误差平均值为7.88%，明显低于上述3种滤波算法的总误差平均值，说明本文滤波算法对于背包点云数据的滤波具有较高的可行性。LAStools、MCC、CSF算法的Kappa系数平均值分别为45.30%、55.90%、47.75%，而MPD滤波算法的Kappa系数平均值为79.48%，明显高于上述3种滤波算法的Kappa系数平均值，说明本文滤波算法对背包点云数据的滤波具有较高的鲁棒性。

图 5. 4种滤波算法的滤波误差精度及Kappa系数分布。(a)总误差对比；(b)Kappa系数分布

Fig. 5. Filtering error accuracy and Kappa coefficient distribution of four filtering algorithms. (a) Comparison of total error; (b) Kappa coefficient distribution

下载图片 查看所有图片

图6所示为4组样本的真实DTM与4种滤波算法拟合出的DTM。可以看到，基于背包LiDAR点云多尺度渐进DTM构建方法相较于LAStools、MCC、CSF滤波算法能较好地选出地面点云数据，并能在一定程度上保留地形细节。如图6所示，实线处为滤波效果较好的区域，虚线处为滤波效果较差的区域。可以看到，MPD滤波算法不仅能将公路上的行人和车辆等非地面点剔除干净，还能较好地滤除植被覆盖等地形细节。而在图6所示的虚线框处，LAStools、MCC、CSF 3种滤波算法虽然能较好地剔除非地面点，但对于森林覆盖的地面点分类误差较大，公路边缘大量的山体点云数据被视为地物点，表明上述3种传统的滤波算法在复杂地表处的细节特征保留较少，滤波效果较为粗糙。因此，将上述4种滤波算法拟合出的DTM与真实DTM进行比较，可以发现在复杂山路地形中MPD滤波算法拟合出的DTM与真实DTM最为相近，相较于其他3种滤波算法，MPD滤波算法的滤波结果更为细腻，保留的复杂地形信息更多。总体而言，本文算法相较于LAStools、MCC、CSF这3种滤波算法能更好地用于背包点云数据的滤波。

图 6. 不同滤波算法构建的DTM

Fig. 6. DTM constructed by different filtering methods

下载图片 查看所有图片

4 结论

针对传统的滤波算法不适用于复杂地形中背包点云数据处理的问题，提出一种基于背包LiDAR点云多尺度渐进DTM构建方法。该方法结合了形态学运算和曲面插值方法的优点，可根据窗口尺寸和地形动态调整阈值，实现自适应滤波，可适用于大多数复杂地形下背包点云的滤波，具有良好的鲁棒性和适用性。本文方法在多尺度窗口下，采用最小二乘法对地面种子点进行曲面拟合，渐进构建DTM，结合形态学开运算逐步获取准确地面点云数据；同时，实现自适应参数阈值的设置，高差阈值T_H和S_th分别随窗口及地形变化而自适应调整，滤波算法受地形变化影响较小，因此可在地形坡度剧烈变化地区获得较好的滤波效果。此外，利用本文方法对4组背包云数据进行滤波实验，针对这4组数据滤波的总误差平均值为7.88%，Kappa系数的平均值为79.48%，与LAStools、MCC、CSF 3种滤波算法相比，本文方法的平均总误差最小，平均Kappa系数最大。实验结果表明，本文方法对于背包点云数据可取得较好的滤波效果，能有效保护地形细节，在实验样本中均获得优于传统滤波算法的效果。