基于PTD和改进曲面拟合的高山区水电工程机载激光雷达点云滤波方法

朱依民，田林亚，毕继鑫，林 松

(1.河海大学地球科学与工程学院，江苏 南京 211100； 2.浙江华东测绘与工程安全技术有限公司，浙江 杭州 310014)

我国水能资源空间分布不均衡，具有明显的空间聚集现象。水利水电工程多位于西南地区的四川、云南、贵州、西藏等地，这些地方具有山高坡陡、河谷狭窄、植被茂密、气候复杂多变、测区地类地物复杂多样等特点，致使传统对地观测技术无法及时满足水电工程的勘测设计需求[1]。随着无人机技术[2-3]的发展，基于无人机倾斜摄影测量技术的地形图测绘得到了很好的应用[4]，但该技术仅适用于地表较为裸漏的测区，而对于植被茂密的地区，植被下方无法通过影像匹配得到点云数据。机载激光雷达(light detection and ranging，LiDAR)测量技术作为一种主动式对地观测技术，集激光测距技术、差分GPS技术、惯性测量技术于一体，具有受外界环境影响小、自动化程度高、采集数据时间短、测量精度高等优点。与无人机倾斜影像立体匹配技术相比，机载LiDAR测量具有较强的穿透力，能快速采集包含大量地物三维坐标信息的点云数据，特别是其独特的扫描式测量方法使其在获取高山区水电工程的丰富地表信息方面有着得天独厚的优势，为快速、精准、自动生成水电工程地表模型提供了可能[5]。

机载LiDAR点云数据滤波是重建水电工程地表模型的重要环节，但高山区水电工程区域植被茂密、山区高度差大、河流宽度窄，致使基于数学形态学的滤波算法[6]、基于坡度变化的滤波算法[7]、基于渐进三角网加密(progressive triangulated irregular network densification, PTD)的滤波算法[8-9]、基于曲面拟合的滤波算法[10]以及基于分割的滤波算法[11]在高山区机载LiDAR点云数据滤波中的应用效果不佳，后期仍需大量人工干预三角网建网成果。李乐林等[12]采用金字塔策略按分辨率从高到低的方式对点云数据进行格网分层，通过局部统计分析的方法确定高差阈值，进而利用薄板样条插值法从上至下对不同层格网控制点的高程进行插值，最终得到最下层格网的地面点云集合，但该滤波方法的精度依赖于顶层格网种子点的选取以及通过统计分析得到的高差阈值的合理性。杨铭等[13]通过改变八叉树模型节点的尺寸，将复杂地形的点云分解成大量的山坡地形点云，然后对每个分解部分利用坡度进行地面点滤波和地面点集合的合并，该方法中距离和坡度阈值的设置以及八叉树的分割标准将直接影响最终的滤波效果。孙蒙等[14]利用统计学中的偏度与峰度思想实现点云滤波，首先对原始点云进行基于偏度与峰度变化曲线的初次滤波，然后对初次滤波获取的地面点进行曲面拟合，最后对曲面拟合高差进行统计和判断，以此完成点云二次滤波获得地面点，该方法需要持续将高程最大点作为非地面点进行剔除直至偏度小于等于零为止，因此该方法仅适用于平坦区域的点云滤波，并不适用于起伏较大的山地点云滤波。郑缉涛等[15]提出一种结合可变半径滚动圆和B样条拟合的点云滤波方法，该方法利用一个可变半径圆环从面向地心侧的方向滚过由点云构成的扫描线，保留扫描线上被圆环滚过的点并作为地面点，对这些地面点进行均匀差值后利用B样条拟合方法进行曲面拟合，最后对曲面拟合高差进行阈值判断完成点云二次滤波获得地面点集合，但该方法在滤波过程中进行的多次插值和采样会丢失原始点云数据信息。针对上述问题，本文提出一种基于PTD和改进曲面拟合的机载LiDAR点云滤波算法，算法思想是先对点云数据进行预处理，再使用PTD算法进行初次滤波获取部分地面点集合，然后基于格网划分对每一个格网进行曲面拟合实现二次滤波，最后通过点云的精细化处理得到精确地面点的集合，研究成果为机载LiDAR测量技术在高山区水电工程中的应用提供技术支撑。

1 基于PTD和改进曲面拟合的高山区点云滤波方法

1.1 点云数据的预处理

1.1.1点云去噪

机载LiDAR扫描系统在扫描时会受诸如多路径效应、飞鸟等影响而产生低位误差和高位误差等噪声点[16]，绝大多数滤波算法都是将点云的最低点作为地面种子点进行处理，为了避免噪声点对后期点云数据滤波产生影响，应先对原始点云数据进行去噪处理。常用的去噪方法有局部平面拟合法、频率域法、高程分布直方图法等，由于噪声点一般较少且在高程空间分布上比较孤立，本文结合高程分布直方图和KD树对点云数据进行去噪。

1.1.2点云回波次数判断

大多数商业LiDAR系统都可以从点云数据中获取4次回波信息[17]。机载LiDAR测量系统发射的激光脉冲在高山区水电工程传播路径上会遇到植被、建筑物、地面等不同地物，进而产生不同的反射信息，且激光脉冲具有一定的穿透力，穿过植被的树叶或枝干会产生多次回波，但遇到建筑物或地面时只有一次回波信息[18]。

通过分析不同测量目标点云的回波信息，可以发现单次回波的点云一般是地面点、建筑物点以及植被的冠层点和枝干点等；多次回波的首次回波点云通常为高大植被的冠层点，中间次回波点云主要为高大植被的枝叶点以及低矮植被点，末次回波点云一般是植被的中间枝叶点和地表反射的激光点。综上，可以利用回波次数对点云数据进行预处理，在具有单次回波、多次回波的末次回波的点云集合中获取地面点，减少非地面点的数量，可提高后期滤波操作的可靠性和算法的效率。

1.2 PTD算法

如图1所示，P1、P2、P3为地面种子点，P4为待判断点，ρ为P4到P1、P2、P3形成的平面的距离，β1、β2、β3分别为P1、P2、P3和P4的连线与P1、P2、P3形成的平面的夹角。为了方便数据索引以及后期计算，要先对经过预处理的点云数据进行格网划分，格网边长为区域内建筑物最大尺寸，这可保证每一格网都含有地面点；其次，将每个格网内的高程最低点作为地面种子点，并以种子点构建初始PTD模型；遍历区域内所有待判断点，对其进行距离和角度的阈值判断，将符合条件的纳入至初始PTD模型中，并更新PTD模型，对剩余点进行迭代判断和PTD模型更新，PTD算法的具体流程见文献[19]。在本文算法中，由于渐进三角网加密是为下一步的曲面拟合算法提供一定数量的地面种子点，并且考虑到PTD算法的耗时较长，所以在此进行两次迭代计算，并把得到的地面点集合作为下一步改进的曲面拟合算法的种子点。

图1 PTD算法的原理

1.3 改进的曲面拟合算法

曲面拟合的思想是将任意复杂空间曲面表示为一个二次曲面函数，函数值代表拟合得到的对应平面点的高程值，再将拟合得到的高程值与实际的高程测量值相减便可得到拟合高程差。对于地面点而言，拟合高程差较小，而非地面点的拟合高程差较大，因此可根据这一特性设置一个合理的拟合高程差阈值来分离地面点和非地面点。

令xn、yn、zn为第n个地面点云的三维坐标值，设地表二次曲面拟合模型的方程为

Z=a0+a1X+a2Y+a3X2+a4XY+a5Y2

(1)

其中X=(x1,x2,…,xn)TY=(y1,y2,…,yn)T

式中：a0、a1、a2、a3、a4、a5为方程的待求参数，其误差方程矩阵形式为

(2)

l=(z1,z2,…,zn)T

当n≥6时，可解算最佳未知参数进而拟合该空间曲面，由于是同精度的独立观测，所以将由不同测量值的权构成的权阵P设为单位阵，根据VTPV取最小值的原则，通过下式求解参数值：

(3)

将解得的参数值代入式(1)便可得到拟合曲面方程，将点云的平面坐标值带入拟合方程中可计算拟合高程值，再通过比较拟合高程差与阈值的大小实现点云的二次滤波。

图2 基于地面种子点的格网划分

图3 单个格网的曲面拟合

1.4 地面点的精细化处理

经过前面的点云去噪、回波判断、渐进三角网加密初次滤波以及改进曲面拟合算法的二次滤波，原始的机载LiDAR点云数据被分离成地面点集合和非地面点集合。此时地面点集合还包含少量的低矮植被点，究其原因是低矮植被点与地面点紧紧相连，进行距离和坡度阈值判断时很难完全将其去除。为了进一步提高滤波的精度，本文利用点的几何特征对这部分点云进行进一步剔除，对于经过上述步骤处理得到的地面点集，遍历计算每一点及其K邻域内点的高度方差σ2，该值在植被处较高，而在地面及平面建筑物屋顶处取值较小，因此可以设定一个合适的高度方差阈值进一步剔除地面点集合中夹杂的低矮植被点。

2 实例验证

2.1 数据源

为了验证提出的滤波算法在高山区水电工程机载LiDAR点云数据滤波中的可行性，选取DJI-M600搭载的HS-600机载LiDAR测量系统(图4)的实测数据进行试验，测区面积约323 072 m2，共扫描点云474 339个，点云密度1.47个/m2，测区地形起伏大、植被茂密，且存在部分建筑物(图5中黄色方框区域)和水域(图5中红色方框区域)等地物类型，完全满足试验的需要，原始点云数据如图5所示。

图4 DJI-M600搭载HS-600机载LiDAR测量系统

图5 高山区水电工程原始点云数据

2.2 试验结果及分析

首先对原始点云数据进行去噪和回波次数判断处理，进而基于PTD算法对经过去噪和回波次数判断的点云进行初次滤波。在PTD算法中需要对点云进行虚拟格网划分，且格网的大小要大于等于测区内最大建筑物的边长才能保证将建筑物全部滤除，由于本试验实际测得测区内建筑物的最大边长为24.5 m，所以试验中将格网边长设置为25 m。根据测区的实际情况，在试验中通过设置不同的阈值来观察PTD的效果，试验发现，当设置的三角网加密距离阈值ρ0=0.5 m，角度判断阈值β0=10°时，采用PTD算法既能提取足够数量的地面点集合，又能得到一个分布比较均匀的地面种子点集合，这有利于为接下来曲面拟合的二次滤波提供有效的地面种子点。利用初始种子点生成的三角网如图6所示，经过两次迭代加密得到的地面点如图7所示。

图6 利用初始种子点生成的三角网

图7 基于两次三角网加密处理的地面点

图8 基于改进的曲面拟合算法滤波后的地面点

经过上述步骤，原始的机载LiDAR点云数据被分为地面点集合和非地面点集合，但此时地面点集合仍包含少量的低矮植被点，设置0.5 m的高度方差阈值剔除低矮植被点，实现高山区水电工程地面点的精细化处理，精细化处理结果如图9所示，得到的数字高程模型如图10所示。

图9 基于精细化处理的地面点云

图10 基于精细化地面点云的数字高程模型

目前，定量评价滤波算法误差的指标主要有两个：①国际摄影测量与遥感学会(ISPRS)提出的滤波评价准则，即点云误分率；②滤波前后点云的均值与标准差的差值。本文采用常用的点云误分率对提出的算法滤波效果进行定量分析，第Ⅰ类误差E1指将地面点分类为非地面点的误差，第Ⅱ类误差E2指将非地面点分类为地面点的误差。E1和E2反映了算法的适应性，总误差Eall反映了算法的可行性，E1、E2、Eall的计算公式为

E1=b/(a+b)

(4)

E2=c/(c+d)

(5)

Eall=(b+c)/(a+b+c+d)

(6)

式中：a为滤波算法中正确分类的地面点数量；d为正确分类的非地面点数量；b为地面点误分类成非地面点的数量；c为非地面点误分类成地面点的数量。

为了进行横向对比分析，采用传统的形态学滤波算法、区域生长滤波算法、PTD算法和曲面拟合滤波算法分别对同一山地点云数据进行滤波并进行精度评价，计算结果如表1所示。

表1 各滤波算法的精度计算结果 %

由表1可知，与PTD算法、曲面拟合滤波算法、区域生长滤波算法和形态学滤波算法相比，本文提出的基于结合PTD和改进曲面拟合的滤波算法的E1、E2和Eall的最大减小幅度分别达6.25%、1.82%和2.93%，表明本文提出的滤波算法更适用于高山区水电工程机载LiDAR点云数据的滤波处理。

3 结 语

本文提出了一种综合点云回波特性、PTD算法、改进的曲面拟合算法和地面点精细化处理的高山区深植被点云滤波方法，采用DJI-M600搭载HS-600的机载LiDAR测量系统的实测数据进行试验，利用ISPRS提出的点云误分率对滤波的精度进行评定。结果表明，与传统的PTD算法、曲面拟合滤波算法、区域生长滤波算法和形态学滤波算法相比，基于PTD和改进的曲面拟合的滤波算法的第Ⅰ类误差、第Ⅱ类误差和总误差的最大减小幅度分别为6.25%、1.82%和2.93%，验证了基于PTD和改进曲面拟合的滤波算法在高山区水电工程机载LiDAR点云滤波中的适用性。