姚杰,侯至群,朱大明,熊建华,付志涛,杜思雨,袁付松
(1.昆明理工大学国土资源工程学院,云南 昆明 650093; 2.昆明市测绘研究院,云南 昆明 650051; 3.昆明市子午环测绘咨询服务有限责任公司,云南 昆明 650000)
机载LiDAR技术已广泛应用于地面数字高程模型(DEM)获取[1]、道路提取[2]、电力线提取[3]、森林参数估计[4]、三维城市模型建立[5]等领域,滤波是这些应用中一项必要的处理环节,是一个过滤地物点,保留地面点的过程,目的是为了突出重要信息。近些年来,许多专家学者对点云滤波算法进行了针对性的研究,主要有基于邻域点三维关系的滤波算法,如渐进加密三角网滤波(PTD)[6]、坡度滤波[7]、形态学滤波[8];基于平面拟合滤波算法,如最小二乘拟合滤波[9]、RANSC滤波[10];基于物理过程模拟和滤波,如布料滤波等[11]。各种滤波算法经过不断改进,在各自应用领域都实现了较好的滤波效果。但是在面对大面积综合复杂地形区域时滤波效果仍不理想。如在地形起伏较大区域,PTD滤波、CSF滤波、形态学滤波难以保留陡坡和山顶等地形特征信息,坡度滤波能够较好地保留地形特征信息但又难以过滤植被和房屋分布密集区域。针对这一问题,本文提出一种组合滤波方法,首先利用改进渐进加密三角网滤波算法对LiDAR点云数据进行整体处理,再通过邻域高差滤波对缺失的陡坡和山顶进行补全,最后通过平面拟合滤波去除近地点噪声,取得了较好的滤波效果,得到了完整且光滑的地面点云数据;并利用同名点较差的方法对地面点的整体精度进行了评定。
组合滤波由三种滤波算法组合而成,本文重点介绍三种滤波算法的原理和处理结果,预处理、裁剪等基础步骤不再详述,算法流程如图1所示。
图1 组合滤波算法流程图
渐进加密三角网算法(PTD)[12]最初是Axelsson提出的,算法通过预先生成一个稀疏的三角格网,再通过逐级迭代加密的方式进行分类,直至所有点云分类完毕。Zhao等人对PTD算法进行了改进[13],提出了IPTD滤波算法(Improved Progressive TIN Densification,IPTD)。该算法包括三个步骤:①利用形态学方法获取潜在的地面种子点;②获得准确的地面种子点;③建立基于TIN的模型,对TIN进行迭代致密化。量取最大建筑物尺寸作为初始格网大小,对点云数据进行格网化处理,对于不含建筑物的点云数据,以默认值作为初始格网大小。取格网内的高程最低点作为起始种子点。遍历所有待分类点,查询各点水平面投影所落入的三角形,计算点到三角形底面的距离及点到三角形三个顶点与三角形所在平面形成夹角的最大值,将其分别与迭代距离与迭代角度进行比较,如果小于对应阈值,则将此点判定为地面点,并加入三角网中;如果大于对应阈值,则去除。重复此过程,直至所有地面点分类完毕。
但是,按照格网最低点先向下逐级致密化再向上逐级致密化的算法原理,如图2所示,山顶点和坡顶点需事先已知。如果初始网格中没有包含这些点,那么陡坡点就不会被保留下来。尽管可以通过调整迭代角度和停止迭代三角网边长来逼近坡顶,但即使调整参数也不能完全还原陡坡点和山顶点,且过度的调整参数会使得其他类别的误差增多。
图2 IPTD向下致密化原理
邻域高差滤波是利用种子点与邻域点高差阈值来获取小于阈值目标点的原理进行滤波的算法,通过控制高差阈值来控制点云的坡度大小。从坡度计算公式(式1)可知,坡度是通过高差与水平距离的比值计算所得,那么通过控制一定半径范围内高差阈值,就可实现保留陡坡点云、过滤非地面点云的滤波效果,算法原理如图3所示。根据邻域高差滤波原理可知,在邻域搜索半径大于最小点间距的情况下,搜索半径越小,则对非地面点的过滤效果越好。但是如果邻域搜索半径小于最小点间距,那么该算法会失去滤波功能,所以滤波结果会存在搜索半径距离的近地点噪声,并且大面积的建筑物屋顶点云因坡度满足阈值也不会被过滤。综上所述,邻域高差滤波作用的区域有限,无法大面积应用,但是其保留陡坡和山顶处点云的特点却能够很好地补充IPTD滤波算法的不足。
图3 邻域高差滤波原理
坡度α=arctan(高差/水平距离)
(1)
平面拟合滤波算法的主要功能是对近地点噪声进行过滤,并使点云表面光滑。其基本步骤包括:①以种子点为参考原点,在原点周围一定搜索半径范围内搜索邻域点,将邻域点利用最小二乘方法进行平面拟合;②计算点云到拟合平面的投影距离;③根据实际需求设置距离阈值,该阈值代表将点云向拟合平面调整的最大距离。平面拟合滤波可以让点云表面更加光滑、与地面更加趋近一致;尤其是对表面粗糙的地面点处理后,点云排列会更加规则。因此,该方法能够有效过滤IPTD滤波和邻域高差滤波获取结果中的近地点噪声。
实验区位于昆明市长水机场附近,面积约为 3.75 km2,如图4、图5所示,点云数据由大棕熊100固定翼飞机搭载RIEGLVQ-1560i机载激光雷达设备于2020年采集所得。总点数4.9亿个,最大点密度 2 616 per/m2,最小点密度 1 per/m2,平均点密度为 131 per/m2。本文使用LiDAR360、ArcGIS等软件对点云数据进行相关处理和分析。
图4 研究区全景
图5 局部细节展示
(1)IPTD滤波
首先采用IPTD算法对原始LiDAR数据进行初次滤波,分离地面点与非地面点。经过反复实验,确定参数设置为:最大建筑物尺寸 400 m,最大值迭代角度10°,地形坡度角为88°。当边长小于 5 m时减小迭代角度,小于 1 m时停止构建三角形。经过滤波后,得到地面点云如图6所示。
图6 IPRD滤波后得到的地面点云
从结果看出,IPTD点云滤波后陡坡和山顶处点云容易被当作非地面点错误过滤掉,如图7所示。为了正确描述缺失处的地形特征,在小比例尺地形测绘时,目前通常做法是利用插值[14]方法进行近似补全,但是,在大比例尺地形测绘时,插值方法会导致特征地形处出现失真。采用本文提出的方法,可以有效减少这类错误的出现。
图7 部分陡坡等地形变化特征点
(2)补全缺失点云并进行邻域高差滤波
通过对初次滤波数据的检查,矢量裁剪了15块缺失区域(如图8所示),一类为自然区域,如自然土坡、山体陡坎边缘、山顶等区域,如图9中的2、15部分;另一类为人工区域,如道路护坡、挡土墙等区域,如图9中的7、9部分。这些区域在构建DEM时非常重要,因此有必要对这些区域进行补全并进行第二次滤波。
图8 15块缺失区域
图9 部分缺失区域补全示意图
根据每块缺失地形的最大坡度计算高差阈值,经邻域高差滤波后不仅过滤了大部分非地面点数据,还完整的保留山坡的全部信息。但是由于大量的邻域搜索计算,整体滤波计算时间较长,且近地表面搜索半径范围内存在一定的噪声,如图10所示。上文2.2中提到,这是由于邻域高差滤波算法对小于搜索半径的点不敏感所致。针对这一问题,提出利用平面拟合滤波的方法对近地点噪声进行过滤。
图10 邻域高差滤波结果表面噪声
(3)平面拟合滤波
在前两轮滤波之后,地面上大部分的地物点都已经被过滤去除。针对近地点噪声问题,提出利用平面拟合滤波算法,首先通过一定范围内的邻域点拟合平面,并计算邻域内所有点到拟合平面的距离;再根据实际噪声高度,设置相应的距离阈值,修改噪声点高程;最终得到了光滑的地面点数据,如图11所示。
图11 平面拟合滤波平滑表面
经过上述三轮滤波处理后,得到了平滑的地面点数据。但是,补全后的地面数据是分散的,故需将15处缺失区域地面点和初次滤波得到的地面点数据进行融合,得到完整的地面点数据。因与初次滤波得到地面数据运用的算法不同,融合后会出现点云密度不均匀现象,如图12所示,再利用重采样的方法即可获得均匀且完整的地面点云数据。
图12 补全区域与初始地面点密度不一致
因LiDAR激光穿透力有限,建筑物和茂密植被区域仍然会出现空洞,按照《激光雷达数据处理技术规范》CH/T 8023-2011中 1∶500地形图的DEM格网间距为 0.5m×0.5 m的要求设置,利用反距离权重插值[15]的方法内插出空洞区高程值,格式为TIFF文件。
点云滤波评价方法有定量分析、定性分析等评定方法,也有通过构建DEM后选点与实测点进行交叉验证、同名点较差等评定方法。因本实验数据为实测数据,故采用构建DEM后采取同名点较差的方法更加精确可靠。
为了对整体测区精度进行准确评定,将测区划分为平地、丘陵、山地等区域,在不同区域内均匀选取一定数量的高程点。通过计算选取的DEM高程点与实测检核点的高程残差△h的中误差来评定DEM模型的精度,从而达到评定组合滤波结果精度的目的。实验区检核点由GNSS系统采集,经过平面、高程基准转换统一后,在DEM相同位置选取同名点共189个,如图13所示,按照中误差计算公式(式2)计算。
图13 同名点分布
式中,mh为公共点高程中误差,单位为(m);△hi为检核点高程较差,单位为米(m);n为参与精度评定的点数为参与精度评定的点数。经过计算,不同区域同名点较差的中误差如表1所示,其中平地、丘陵地和山地的精度都远高于标准限差值,所得DEM数据准确。
表1 同名点较差中误差
从误差结果分布图14可以看出:①总体误差值均优于规范限差并呈正态分布;②平地和山地中误差较小,平地精度又优于山地,平地主要为农田和空地,山地主要为植被;丘陵地中误差较大,丘陵地中主要地物为厂房和居民房。丘陵地区存在较大的误差值,主要是本文提出的组合滤波效果在房屋分布密集区域并没有独特的优势和针对性改进,并且构建居民区和厂区DEM时运用了内插算法,且内插面积较大,也是导致误差值较大的原因之一。因此,本文提出的组合滤波算法在较多建筑物区域处理效果方面还有待进一步提高。
图14 误差分布
点云滤波是LiDAR点云数据处理过程中最重要的一环,从地物影响因素的角度出发,不仅要考虑如地形、建筑物、植被等因素影响,还要考虑多种因素的联合影响,在单一滤波算法中加强对一种因素影响的控制同时会减弱对另一种影响因素的控制,难以适应大面积地形复杂的区域。本文提出一种组合滤波的方案,首先对初始LiDAR点云数据利用IPTD滤波算法去除了大部分的地物点,通过邻域高差滤波的方法对IPTD滤波在陡坡、山顶等区域的缺失点进行了补全,最后通过平面拟合滤波的方法对数据进行了平滑,消除了近地点噪声;经过融合、重采样、内插等处理之后,制作了实验区 1∶500地面DEM,通过同名点较差的方法进行精度评定,精度均满足规范要求,证明组合滤波是一种有效可靠的滤波算法,并取得了如下结论:
(1)相比以往单一滤波方法,组合滤波方法的适应性更强,可以做到各滤波算法之间取长补短,有效降低了单一滤波算法的不足对结果的影响。
(2)组合滤波方法适用于大面积综合场景下的LiDAR数据,对自然和人工陡坡、山顶等地形的滤波效果好;但是也存在处理步骤多、智能化程度低、对建筑物密集地区的过滤能力不足等问题。