张 兵，武殿文

（1.内蒙古电力勘测设计院有限责任公司，内蒙古 呼和浩特 010020；2.内蒙古能源规划设计研究院，内蒙古 呼和浩特 012100）

0 引言

风电场地形图测绘如今普遍采用无人机低空航测技术，以便获取测区高精度数字正射影像图(digital orthophoto map，DOM)和数字地表模型(digital surface model，DSM)。通过功能强大的无人机数据后处理软件一般能够得到较为理想的测区DOM，即使在非地面点复杂地区也不会有太多影响其精度的地方。但DSM对于非地面点的反应却十分敏感。对于植被十分茂密的风电场测区，往往由于树木高度的不确定性很难得到精确地形数据，进而影响风电场地形图测绘的准确度。所以必须剔除非地面点对数字高程模型(digital elevation model，DEM)的影响。剔除非地面点的解决方法，传统主要采用手工采集，然后减去树木、房屋等地物点的高度，得到地面点的高程点绘制等高线，或在每个区域内部采取野外实测树高，或者在航测内业采用立体测图的办法获取树高，对每个植被覆盖区域采取减去固定高程差值的办法来得到整个风电场区域的DEM。上述传统方法采取人工干预的方式来剔除测区内植被或者房屋等非地面点，人工工作量很大、效率很低。另一种途径是借鉴激光雷达(light detection and ranging，LiDAR)点云的处理方法对多视影像匹配点云处理，将已获取的点云数据中地面点与非地面点分离，获取剔除非地面点后的DEM，用于风电场矢量地形图绘制。这种办法人工参与少、作业效率高，在实际应用中反复验证点云滤波后DEM精度情况显得尤为重要。

1 LiDAR点云滤波技术的发展

点云数据滤波方法最早来自激光点云的数据处理，一直是激光点云应用最前沿和最难解决的棘手问题，机载LiDAR点云数据的滤波算法仍相对落后，许多学者对点云滤波算法进行了研究，并提出了多种滤波算法，可概括为2种：①利用点云数据的几何特性(三维坐标)进行的滤波； ②利用点云数据的回波特性(回波强度、多回波信息) 进行的滤波[1-2]。LiDAR点云滤波技术通过近些年的发展，出现很多基于点云数据的几何特性(数据本身)的滤波方法，有基于坡度、基于曲面拟合、基于数学形态学和基于布料模拟等。基于数据本身的滤波方法，这一过程是利用从数据内部提取的有用信息建立判别规则，实现过滤非地面点，获取只包含地面点的点云数据的目标，这种滤波方法的对象不只是滤除植被，还包括了人工地物、移动物体等表面噪声。基于数据本身的滤波目的是滤除点云数据或DSM的表面噪声，在滤除高低不一的建筑物或孤立植被方面具有较好效果，但在大范围植被覆盖地区或成片森林，立体影像匹配生成的点云数据通常不包括地面点，因此不能满足这类算法所要求的局部邻域内最低点就是地面点的初始条件，无法有效滤除植被。但是针对多视影像匹配点云的滤波只能采用基于数据本身的滤波方法。目前基于数据本身的LiDAR点云滤波技术主要使用的方法包括基于坡度的滤波算法、基于不规则三角网滤波方法、基于移动曲面拟合的滤波算法、基于数学形态学滤波算法以及基于布料模拟的滤波算法[3-5]。本文将分别论述以上几种滤波方法，并重点研究基于布料模拟的滤波算法。将基于布料模拟滤波算法研究成果应用于风电场地形图等高线绘制。

2 多种基于数据本身LiDAR滤波技术的算法原理与分析

基于坡度的滤波算法由VOSSELMAN首先提出，其思想是依据坡度的差异区别出地面点和非地面点，但对整个点云区域采用固定坡度阈值的方式不够合理[6]。基于坡度的滤波算法原理认定地形表面为平缓光滑曲面，局部区域地形急剧变化很大可能由非地面点引起，相邻点间高差值和距离值的滤波函数被定义为地面点归类函数，如式(1)所示。

式(1)中点pi与pj的差值不大于Δhmax(d)即可将pi判定为地面点，否则为非地面点。

图1中基于不规则三角网滤波算法是通过初始地面点p1、p2、p3与待判定点p构建不规则三角网，分别计算d、a1、a2、a3与设定的阈值比较，如果小于阈值将p归类为地面点并加入初始三角网，通过反复迭代计算直到没有新的地面点结束。

图1 不规则三角网略图

基于移动曲面拟合算法强调地形的平滑性，将点云数据二位排序获取地面点种子区域，计算种子点曲面方程并将临近点代入计算拟合高程值，若不低于设定阈值，判定为地面点并吸收到二次曲面方程中，通过反复迭代计算直到所有点云判定结束，在改变网格尺寸和设置不同的窗口领域大小时，为了有效地将地面点和非地面点分开，需要设置不同的高差阈值[7]，如式(2)所示。

基于数学形态学的滤波算法对点云数据先进性腐蚀运算再进行膨胀运算，将区域内最小和最大高程值填充到结构元素内部，然后运用开闭运算反复迭代出点云数据内地面点和非地面点。形态学滤波主要依靠开运算实现点云类属判断[8-9]。主要运算式如式(3)～式(6)所示。

基于布料模拟滤波算法是根据牛顿第二定律原理在模拟过程中内部力和外部力同时作用且随时间变化，假设布料足够柔软黏住地面则布料最终形状会是DSM，考虑将点云倒置给布料赋予硬度因素则布料形状贴近数字地面模型，物理表达式如式(7)所示。

在实际应用中，统计很多实际案例发现基于坡度的滤波算法对多数地形都具有适用性，需要根据不同地形设定坡度阈值，不适合复杂地形，坡度阈值的选定也较为困难。基于不规则三角网的滤波算法适用于地形起伏不大的地区，但是构建三角网运算量很大，距离和角度阈值设定很复杂且不易剔除植被。基于移动曲面拟合的滤波算法原理简单算法速度快，但函数高程阈值很难判断，窗口尺寸设置不好把握。基于数学形态学的滤波算法简单易于实现，但在处理原始点云数据中容易丢失数据，固定的滤波窗口和阈值使滤波精度不高。而基于布料模拟滤波算法在实际应用中适用于大多数地貌类型，对于各种地物类型的滤波效果也更为优异，已被更多应用于DSM的点云滤波。在吸收很多人关于点云相关研究的基础上人们建立起来一套大型跨平台开源编程点云库(point cloud library，PCL)，PCL实现了大量点云相关的通用算法，其中就包括各类点云滤波算法。利用PCL模块化的现代编程模板库实现点云滤波算法实践显得十分高效。

PCL可结合多种程序来编译滤波算法，实现点云数据的自动滤波，提供了大量数据类型和数据结构，以及许多最先进的算法对3D数据进行处理，包括滤波、模型重估、表面重建等。PCL中的所有模块和算法都是基于标准的第三方库Boost、Eigen、FLANN、VTK、CUDA、OpenNI、Qhull来传送处理数据的，避免多次复制系统中已存在的数据的需要。PCL利用OpenMP、GPU、CUDA等先进高性能计算技术，通过并行化提高程序实时性。从算法的角度，PCL是指纳入了多种操作点云数据的三维处理算法。每一套算法都是通过基类进行划分，把贯穿整个流水线处理技术的所有常见功能整合在一起，保持整个算法实现过程中的紧凑和结构清晰，提高代码的重用性、简洁性。

3 CSF算法原理与程序实现

基于布料模拟滤波(cloth simulation filter，CSF)算法相较于其他滤波算法具有算法精度高、适用面广等优点，本文把基于布料模拟的滤波算法作为无人机多源数据匹配生成DSM的主要滤波方法。布料模拟是物理过程的模拟，符合牛顿第二定律原理。模拟布料落到倒置的点云表面上，分析布料的节点与相应的点之间的相互作用，以约定的时间内模拟布料的形状，而布料节点在垂直方向上的位移依照式(8)所示。

算法的核心是根据PCL使用编程语言完成CSF算法编译。其大致的实施过程是将原始点云数据倒置，设定布料模拟参数和初始位置。然后将点云数据与栅格质点投影到同一平面，结合式(6)和式(7)计算模拟布料所有质点的垂直位移，使用位移后的布料质点与初始点云之间进行距离计算，并与设定阈值比较，距离小于分类阈值则将此点归为地面点，距离大于分类阈值则将此点归为非地面点。如图2所示是基于CSF算法简易流程图。

图2 CSF简易流程图

具体在PCL中，根据提供的多种开发语言的点云处理技术算法源码。首先在PCL中定义函数库以及需要使用的规范文件，读取和保存LAS格式点云数据。在进行CSF过程中关键在于滤波参数的设置，需要设置质点为不可移动，设定布料网格分辨率和布料硬度，设定布料模拟质点与原始点云的距离阈值以及重复计算最大迭代次数。以下是CSF基于LAS格式点云python简易源码。

# coding: utf-8

import laspy

import CSF

import numpy as np

import open3d as o3d

inFile = laspy.file.File(r” in.las”,mode=’r’)# read a las file

points = inFile.points

xyz=np.vstack((inFile.x,inFile.y,inFile.z)).transpose()# extract x,y,z and put into a list

csf = CSF.CSF()

# prameter settings

csf.params.bSloopSmooth =False

csf.params.cloth_resolution =0.2

csf.params.rigidness =5

csf.params.time_step =0.65

csf.params.class_threshold =0.5

csf.params.interations =200

#moredetailsaboutparameter:

csf.setPointCloud(xyz)

ground = CSF.VecInt()#

non_ground = CSF.VecInt()#

csf.do_filtering(ground,non_ground)#

outFile = laspy.file.File(r”ground.las”,

mode=’w’,header=inFile.header)

outFile.points = points[ground]# extract ground points,and save it to a las file.

outFile.close()#

4 基于DSM点云滤波方法的风电场地形图应用

无人机低空摄影制图技术用于提供风电场地形图需要的DOM和DSM。多视影像匹配技术的发展带来更为精细化的2D数据，针对无人机2D数据的风电场矢量地形图绘制已经十分成熟。但是在内蒙古区域内依然存在不少地物植被众多区域，因此在处理风电场矢量地形图时不可避免对非地面点进行过滤和DEM重建。将LiDAR点云滤波技术引入到多视影像匹配DSM滤波，能够提高后续制图成果的效率和精度。如图3所示是加入点云滤波步骤后的风电场地形图制作流程图。

图3 风电场地形图制作流程

数据滤波前后针对滤波区域进行了精度检校。图4为无人机低空摄影测量成果制作的DSM，图5为采用布料滤波算法进行点云滤波后重新构建生成的DEM。图6为采用原始DSM制作的等高线图形，图7为采用DEM制作的风电场等高线图形。通过加入DSM滤波技术，能够显著提高风电场地形图等高线制作的效率和精度。通过以上案例发现，在对原始DSM进行自动滤波后，大部分非地面点都被有效剔除，滤波后重构DEM能够贴合真实地形(图5)，但是在实际检测中发现部分正常的地形点数据经过滤波后被带偏移，部分地区高程点数值改变超过限差以上。针对以上问题在对现有数据进行滤波后再与原有DEM数据进行一次融合算法。该算法原理是首先对2次数据进行匹配，计算机自动锁定区域内所有滤波区域然后计算出每个滤波区域的闭合界线。设定在闭合线里边的数据采用滤波后数据，而在闭合线外部的数据替换成原有数据。在获得更新数据之后，又进行了DEM精度检校。采取以上办法可以很好地解决部分正常数据被滤波数据带偏的问题。最终再次重构DEM并生成等高线数据(图7)，该等高线数据能够真实体现风电场地形。

图4 数字地表模型DSM

图6 原始DSM生产等高线

图7 重构高程模型生产等高线

针对点云滤波精度评定方法有目视定性评估法和统计学定量评价法：①目视定性评估法：是观察滤波后地面点数据是否符合风电场地形图制图要求；②统计学定量评价法：点云滤波结果时定量分析第I类误差、第II类误差和总误差。设定a为地面点误分为非地面点的数量，b为非地面点误分为地面点的数量，c为正确划分地面点的数量，d为正确划分非地面点的数量。误差评价公式如式(9)～式(11)所示。

点云滤波算法一般要求在控制Ⅱ类误差的基础上尽量降低Ⅰ类误差，即在保证分离出的地面点最大程度上是地面点的同时，尽可能分离出较多的地面点[10]。以上3类误差比率系数值越小表明点云滤波效果越好。在多个工程实际应用中统计表明，能够满足风电场测图需求的点云滤波3类误差统计结果应满足Ⅰ类误差≤0.1，Ⅱ类误差≤0.05，总误差≤0.05。

5 结论

LiDAR点云滤波技术发展到今天，很多前人的经验在经过发展之后能更好地对LiDAR点云数据进行滤波。将LiDAR点云滤波技术用于多视影像匹配获取的DSM点云滤波，由于无人机航测数据的局限性只能在一定程度上解决植被等非地面点的干扰来获取DEM，并且通过实际检校发现很难达到更高精度。仅在风电场地形图无人机航测实际应用中通过多次的滤波成果与工测数据对比，发现大多数能够满足风电场地形图测绘精度要求。因此对于植被丰富地区的地形图测绘，为了达到厘米级精度要求，想通过DSM滤波技术方法得到高精度DEM不现实：一方面需要广泛运用无人机机载LiDAR技术采集测区激光点云，通过机载LiDAR多次回波技术直接穿透植被获取真地表地形数据，才是一次性解决高精度真地形模型的办法；另一方面由于单一一种LiDAR点云滤波技术都存在自身算法和阈值设定的局限性，因此很难确定到底哪种方法能够完美适用于各种地形风电场地形图的DSM滤波。对于植被丰富地区的地形图测绘，无法满足更高精度要求，需要更加深入地将多种LiDAR点云滤波技术合并使用，同时积极采纳图像学、人工智能学最新研究成果，提供适用于所有地形场景风电场DSM数据的点云滤波，以期达到更高精度和滤波阈值设定更智能化。