(1)多尺度曲率滤波算法(MCC)。MCC定义一个包含所有点x,y,z坐标的向量Z(s),x和y为平面坐标,z为点云高程。利用薄板样条插值法Z(s)得到连续栅格表面。在栅格表面传递一个3×3的格网定义一个新的向量X(s),其z为3×3栅格的平均高程; 比较Z(s)中各点的高程与X(s)高程和曲率容差t之和的大小,当Z(s)<X(s)+t时,为地面点云。然后Z(s)重新定义为未分类的点进行不断迭代,直到所有的点都低于阈值。MCC结果主要受点云间距和曲率容差阈值影响^[21]。

(2)扩展窗口高程阈值算法(ETEW)。ETEW根据相邻地面点的高程差和地面点与非地面点高程差的差异进行点云分类,通过对比预设半径内点的高程差与预先设置的阈值,剔除非地面点^[22]。首先将点云数据细分成方形单元,保留单元内最低高程的点; 其次计算单元中高程最低点和其他点的高程差,高程差小于阈值时视为地面点,大于阈值则为非地面点。随着单元格和阈值的增大,重复此过程,直到不再剔除非地面点。ETEW分类的效果和单元格大小、坡度因子、迭代次数有密切关系^[22]。

(3)适应三角网滤波算法(ATIN)。ATIN以点到三角形不规则网络(TIN)表面的距离为主要约束条件,进行滤波计算。首先将整体区域划分成小的单元格网,每个单元格只保留高程最低的点(种子点)构建初始稀疏的TIN; 然后计算每个点到最近三角形表面的距离以及点到三角形顶点的直线与三角形表面的夹角。当距离和角度都小于预先设定的阈值时,则将该点归类为地面点; 后续迭代中使用这些基本点来构造新的TIN,计算其余点到表面的距离和角度,并与阈值进行比较,直到未分类的点都不可添加到地面点集时,迭代终止^[18]。该算法分类结果和单元格的选取、初始三角网的大小关系最密切。

(4)渐进形态学滤波算法(PM)。数学形态学滤波算法基于集合理论从图像中提取特征,采用腐蚀和膨胀两个基础算子组成“开”运算对LiDAR点云数据分类。Zhang等^[16]在此基础上开发了PM,其核心原理是逐渐增大窗口尺寸及对应的高差阈值,在保留地面数据的同时,去除不同尺寸的非地面目标的测量值。首先将矩形格网覆盖在点云数据集上,每个单元格网包含一个测量点,其为单元格内的最低高程点。单元格中测量点的高程构成一个初始近似曲面; 然后进行“开”运算,获得次生曲面,将第i次的曲面与第(i-1)次的曲面之间高程差值和阈值进行比较,判断单元格网中的点是否为非地面点^[23]。单元格和起始窗口大小对PM结果影响较大。

(5)基于坡度的滤波算法(SBF)。SBF通过坡度差来分离地面点和非地面点。算法原理是将一个矩形网格覆盖在点云数据集上,每个激光点分配到一个单元网格,如果同一单元格中有多个点下落,则选择高程最低的点作为数组元素; LiDAR点在给定半径范围内,该点与任意其他点之间斜率的最大值小于预先设定的阈值时,被划分为地面点^[15]。单元格大小、搜索半径、最大坡度对该算法影响较大。

(1)参考数据获取。通过对比参考点云数据与不同算法的分类结果,确定滤波精度。参考数据可通过纯人工滤波、GNSS/全站仪测量、自动+手动滤波等方法获取^[22]。纯人工滤波是作业人员根据对研究区的先验知识,通过手动分类实现地面点与非地面点的区分,分类精度受作业人员知识、经验等影响较大。GNSS/全站仪测量所得地面点数据密度较低,且野外工作量大,不易实现。“自动+手动”方法可兼顾点云精度和工作量,已被广泛用于获取参考点云数据^[18]。因此,本研究选用“自动+手动”方式,在Terrasolid 2016软件中处理得到参考点云。首先,通过Terrasolid 2016内置滤波算法剔除杂点和噪声点,初步分离地面点和非地面点。然后,结合无人机LiDAR获取的董庄沟高分辨率照片,开展目视判读,修正植被覆盖度、地形陡峭等敏感区域的分类结果,获取高精度的点云分类数据。

(2)精度评价。点云分类精度评价采用国际摄影测量与遥感学会推荐的方法。该方法计算3类误差,分别为Ⅰ类误差、Ⅱ类误差和总误差^[17]。其中Ⅰ类误差表示地面点云错分为非地面点云的比例,Ⅱ类误差表示非地面点云错分为地面点的比例,总误差代表总体点云错分比例。

Ⅰ类误差=b/(a+b)(1)

Ⅱ类误差=c/(c+d)(2)

总误差=(b+c)/(a+b+c+d)(3)

式中:a为正确分类的地面点数量; b为地面点误分为非地面点的数量; c为非地面点误分为地面点的数量; d为正确分类非地面点的数量。

高精度的地形获取需要完整的地面点云,以便有效监测微地形变化。Ⅰ类误差导致的地面点云缺失难以在后续数据处理中消除,是地形信息准确表达的主要难点,而Ⅱ类误差(非地面点云误分为地面点云)则相对容易去除^[17]。因此,本文滤波中,在尽可能保证地形信息完整的同时,降低各类误差,各参数初始取值为允许范围的下限,以固定步长逐步增加参数取值,直至参数取值上限,对比不同参数取值下的滤波误差,获取参数的最优取值(表2)。

表2 各滤波算法主要参数的最优取值

基于各样本数据的地面点云参考数据,生成0.5 m分辨率的DEM,计算不同样本区域的平均坡度(表1),分析坡度对不同滤波算法精度的影响。结果表明,各算法的Ⅰ类误差与平均坡度呈显著正相关(图4)(p<0.05),而Ⅱ类误差和总误差与坡度的相关性不显著(p>0.05)。MCC和PM的Ⅰ类误差随坡度上升速率最小,ETEW的Ⅰ类误差随坡度上升速率最大,SBF和ATIN的Ⅰ类误差随坡度上升速率居中。坡度小于45°的9个样本区域中,6个样本区域的MCC的Ⅰ类误差最小,坡度最高的3个样本区域(48.54°,49.72°,54.95°)PM的Ⅰ类误差最低。

Ⅰ类误差较Ⅱ类误差和总误差对坡度变化更敏感。原因在于各算法倾向于将陡峭区域整体误分为非地面点,造成Ⅰ类误差的迅速上升,而陡峭区域点云中包含的部分真实非地面点被正确分为非地面点云,缓慢降低了Ⅱ类误差。各算法Ⅰ类误差对坡度变化敏感性差别的原因可能在于,ETEW仅依据高差阈值判断归类点云,不能有效考虑陡峭区域的影响,SBF和ATIN通过选择合适的坡度(或夹角)阈值一定程度上降低了坡度对滤波结果的影响,而MCC利用曲率容差有效考虑了坡度的影响,因此Ⅰ类误差随坡度变化较小^[21]。PM通过能够有效考虑了陡峭区域的地形变化(图1),因此对坡度变化不敏感。

基于参考数据中植被点云与原始点云点数量的比值,计算各样本区域植被覆盖度^[19],分析植被覆盖对滤波结果的影响。结果表明,Ⅰ类误差与植被覆盖度呈显著负相关(p<0.05),而Ⅱ类误差和总误差与植被覆盖度相关性不显著(p>0.05)(图5)。MCC和PM的Ⅰ类误差随植被覆盖度上升下降速率最低,ETEW的Ⅰ类误差随植被覆盖度增加降低速率最高,SBF和ATIN居中。植被覆盖度较低85%的区域,MCC和PM的Ⅰ类误差较ATIN小,而在植被覆盖度大于85%时,ATIN的Ⅰ类误差小于PM,但植被覆盖度达到90%时,ATIN的Ⅰ类误差与MCC接近。

Sithole等^[17]等研究认为Ⅰ类误差随植被覆盖度上升而增大,而本文得出与之相反的结论。本研究区域地形复杂,地形不仅直接影响滤波结果,而且显著影响植被分布,植被覆盖度高的区域比较平缓,Ⅰ类误差较小。而植被覆盖度低的区域,地形陡峭,坡度因素占主导,导致Ⅰ类误差比较高。同时,地形对植被覆盖的影响解释了各算法的Ⅰ类误差对植被覆盖度变化的敏感性差异,也解释了Ⅰ类误差较Ⅱ类误差对植被覆盖度的变化更为敏感。

对各样本区域原始的点云数据进行抽稀,获得不同抽稀密度(5%,10%,25%,50%,100%)的样本数据,利用滤波算法对不同点云密度数据进行滤波,分析点云密度对不同滤波算法精度的影响。结果表明,随着点云密度增加,MCC,ATIN,PM的Ⅰ类误差无明显变化,Ⅱ类误差和总误差缓慢下降,在5%和100%点云密度下的Ⅰ类误差略高于10%～50%点云密度; SBF和ETEW的Ⅰ类误差增加,Ⅱ类误差和总误差下降(图6)。

图4 坡度与滤波算法误差的关系

图5 植被覆盖度与滤波算法误差关系与坡度和植被覆盖度的关系

图6 滤波算法误差随点云密度的变化

通常情况下,低点云密度下裸地与植被不易分离,导致Ⅱ类误差增大^[22-23],因此随着点云密度增加各算法的Ⅱ类误差都呈下降趋势。SBF和ETEW的Ⅰ类误差与点云密度正相关,原因是随着点云密度的升高,固定邻域内的点数增加,大于高程或坡度阈值而被归于非地面点的数量增加,这部分点中包含真实地面点云,它们的剔除提升了Ⅰ类误差。其次MCC,ATIN,PM误差随点云密度变化较为缓慢,原因可能是此3种滤波算法均需利用点云构建连续表面,连续表面的质量影响算法的精度^[15],本文获取的原始的点云密度较高(平均密度26～56个/m²),在其基础上抽稀对连续表面质量影响不明显。