激光雷达(Light Detection and Ranging，LiDAR)技术是目前应用最广泛的高分辨率测量技术之一^[1,2]，根据搭载平台的不同，可分为星载、机载、车载和地面LiDAR^[3]。其中，机载LiDAR系统及其数据处理技术发展较为成熟，主要用于获取大范围高精度数字地面模型和数字表面模型等，而合适的点云滤波算法是获取高精度数字模型的前提^[4]。根据滤波策略不同，激光点云滤波算法主要分为基于表面调节滤波、基于形态学滤波、基于不规则三角网滤波、基于分割与分类滤波、基于统计分析滤波与基于多尺度滤波^[4]。基于表面调节滤波主要是根据不同的误差特性，通过循环迭代操作逐步精化数字地面模型（Digital Terrain Model，DTM）得到高精度DTM，从而分离出地面点^[5,6]，但该方法容易忽略地形细节，将接近地面的地物点误分为地面点^[6]。基于形态学滤波则是通过设定描述地形差异的特定斜率运算符实现点云滤波操作^[7,8]，在陡峭地区具有较强的鲁棒性，并且能够去除较小的非地面地物，保留地形细节,但该方法参数设置依赖于地形特征，在地形复杂区域的适应性需进一步研究^[6]。基于不规则三角网滤波的可靠性和准确性已经被大量的研究所证实，但在处理不连续地形时滤波效果较差^[9,10,11]。基于分割与分类滤波能够充分利用额外的几何与纹理等特征进行点云滤波^[12,13]，可以有效地去除不同尺寸和形状的地物，非常适用于城市地区，但难以有效地应用于植被覆盖密度较高的地区。基于统计分析滤波借助统计学理论对点云数据进行处理，避免了人为调节参数的不确定性，具有较强的鲁棒性^[14,15]，在地物单一的平坦地区，滤波效果较好，但其稳定性低于其他有参方法。基于多尺度比较的点云滤波则借助多种滤波策略对非地面点进行识别，对于一般平坦区域效果较好，但在地形急剧变化的区域效果较差^[16,17]。

以上几类方法通常对特定地形具有良好的滤波效果，但获取的机载LiDAR点云数据通常来源于地物多样、地形复杂的场景，如何提高复杂场景下点云数据的滤波效率，减少参数人为设置提高滤波的自动化程度，提出更为有效且稳健的滤波方法，一直以来都是机载LiDAR滤波研究者关注的热点与难点^[4,18]。综合现有的几种滤波策略，可以在一定程度上提高传统单一滤波方法在复杂场景下的滤波精度，如何寻找有效的算法组合也成为了当前研究的热点^[4]。考虑到布料模拟滤波方法具有参数设置少、可去除微小地物及在地形平坦城区滤波效果好等特点，本文在顾及地形多样性特征的前提下，综合高程归一化处理流程，对布料模拟滤波算法进行改进，有效解决了该方法在应用中由于内力作用而产生的地形适应问题，提高了算法的精度和鲁棒性，具有更为广泛的地形适用性。

布料模拟滤波算法（Cloth Simulation Filtering，CSF）是一种基于表面调节的滤波算法，它巧妙地将物理模型中的“布料模拟”算法应用于LiDAR点云滤波处理^[19]，即通过点云反转使地形倒置，然后使用“弹簧—质点”布料模型^[20]（图1）模拟点云滤波过程。“布料”质点基于重力作用进行下落（图2(b)、2（e）），质点之间的内力作用使质点向上位移（图2（d）、2（g））。在质点的重力作用下，当下落到该质点的最邻近激光点的最低高度以下时，标记为不可移动质点，并定位到该最低高度（图2（c）、2（f））。质点间的内力作用主要通过邻近激光点的上拉作用进行位置修正。通过这种上下物理过程的迭代模拟，确定每一个地形“布料”质点的最终位置，并得到地面点。具体原理如图2所示。

CSF滤波是一种自上而下的表面迭代精化过程，在平地地区具有较好的滤波结果^[1,4]，可将低矮植被、桥梁和建构筑物等地物去除，并且能够保持较强的地形细部特征。

在CSF算法中，内力作用分为相邻两可移动质点间的内力作用和相邻可移动质点与不可移动质点间的内力作用。相邻两个可移动质点间的内力作用需将两质点调整为平行状态（图2）。而相邻可移动质点与不可移动质点间的内力作用所产生的质点位移仅依赖于设定的硬度值，硬度参数越大，内力作用对物体所产生的位移越大。例如，可移动质点与不可移动质点的高程差为d，当硬度参数为1时，上移1/2 d；硬度参数为2时，上移3/4 d（图3）。

CSF算法在单一研究区域中，将硬度参数设置为唯一值，对于不同的地形条件，质点间的相同高程差所产生的内力作用相同，没有考虑到地形的多样性。因而，当迭代次数设置过大或硬度参数设置过小时，低矮建筑物易被误分为地面点；迭代次数设置过小或硬度参数设置过大时，山顶信息被削除的概率增大，算法的地形自适应差。尽管增加滤波算法的迭代次数可减少山顶信息被去除现象，但同时会影响算法在平地地区的滤波效果，出现山顶信息被削除及包含较多地物点等错误（图4），同时增加大量的算法时间^[21]。

张吴明等^[19]为了解决这一问题，利用相邻质点间高程变化小的特性对该算法进行了后处理，可在一定程度上改进陡坡处的滤波精度，但该算法过度依赖于邻近不可移动质点的位置，导致在山顶区域、大量植被覆盖的山地区域应用效果仍不理想。

由于布料模拟滤波的过程是水平“弹簧-质点”模型简单的物理过程模拟，通过反复的循环迭代，最终确定质点位置，每次布料迭代的结果都是上次迭代结果的精化。若布料的初始形态能够大致地反映地形形态，在满足保留山顶信息及滤波精度的同时，迭代时间将会大大缩减。布料初始形态的确定一定程度上等价于对LiDAR点云进行高程归一化预处理，因而，本文针对布料模拟算法的特性以及山地地区的地形特点，研究了顾及地形特征的布料模拟点云滤波算法。该方法使用能够反映地形起伏变化的粗DTM对LiDAR点云数据进行高程归一化预处理，得到归一化高程的LiDAR点云。去除地形特征的LiDAR点云基本处于平面内，利用CSF算法在平地地区滤波效果好的优势，对该高程归一化LiDAR点云进行地形“布料”的物理过程模拟，确定布料形态，算法如图5。最终通过确定的阈值进行地面点的提取。

本文采用局部高程插值的方法进行点云粗差的剔除^[22]。通过搜索某个点（P）的邻近点采用距离加权平均插值法内插出该点的高程值Z，然后计算P点内插高程Z与测量高程的差值，如果该差值绝对值大于阈值ε就剔除P点。

复杂场景地区地形起伏大、地物复杂多样，为了使粗DTM能够较为真实地反映地形起伏特征，本文采用多级格网法。基于坡度动态确定格网尺寸，以格网内最低点作为地面种子点，进行曲面拟合得到粗DTM。一般来说，地形坡度较大的区域，基本没有较大建筑物分布，格网尺寸设置相对较小；而对于地形坡度较小的城市等区域，建筑物分布较为密集，格网尺寸需设置相对较大。根据建筑物基于地形坡度的分布特点，通过初始DTM的坡度数据，动态确定粗DTM格网尺寸的方法，既可以保留基本的地形起伏特征，又不会引入误差影响后期处理精度。多级格网种子点获取的具体流程见图6。

通过能够反映地形基本起伏特征的粗DTM对LiDAR点云进行高程归一化处理，可以避免CSF滤波算法的地形自适应问题，提高算法在大范围复杂场景区域的适应性和滤波精度，同时可以减少迭代次数，降低算法的时间消耗，提高算法效率。

高程归一化的过程主要包括查询、差值计算两个步骤。通过激光脚点的x、y值，查询该点所在的DTM栅格的行列数，并获得该位置处的高程值HDTM，进而与对应像元的高程值h_DTM作差计算，最终确定该激光脚点的归一化高程值Hn。

通过以上步骤得到的高程归一化的LiDAR点云，一定程度上去除了地形起伏特性。将所有点云的高程转换为对应粗DTM的高差，通过重力及内力作用进行上下调节位置，进行地形布料的物理过程模拟，有效地利用CSF算法在平地地区的优势。

通过循环迭代重力、内力作用，直到最大的质点位移变化足够小或迭代次数达到给定值时停止迭代，得到“布料”质点的最终位置。使用点到面的点云间距离计算法，计算高程归一化处理后的点到布料点云的距离，距离小于阈值的归一化点云被分为地面点，其他点被分为非地面点。

将高程归一化的地面点及非地面点与粗DTM叠加运算，即可完成滤波。

本文采用某大学附近3 000 m×1 500 m的复杂场景区域作为实验区域（图7），区域内有平地、河流、山堆及山包等多种地形，且被大量植被覆盖，建筑物及其他地物形态各异、大小不一、分布不均。

本试验数据结合遥感影像，使用渐进加密三角网算法及人机交互方法对研究区的机载点云数据进行分类处理，并通过实地调查验证，保证数据99.99%的分类精度。LiDAR点云数据参数如表1所示。

为了验证算法的精度，按照以下滤波过程进行滤波处理，得到该区域的滤波结果。

数据预处理：地表点云粗差剔除的关键是计算时间和精度，由于时间是点云数量的二次函数^[21]，所以对这种大范围复杂区域，根据研究区域的点云密度，采用5 m×5 m的正方形格网进行内插计算，ε取0.05 m。数据预处理后的点云不仅去除了粗差较大的点，同时去除了局部范围内较难人为识别的粗差数据。

粗DTM的建立：主要采用多级格网最低点作为种子点进行区域内插得到。格网密度越高，采集的最低点数据密度越高，越容易保留地物点作为格网内最低点。为了提高不同地形条件下归一化运算精度，根据研究区数据特征，选取60 m×60 m方格网内最低点作为初始种子点，运用最邻近法进行1 m×1 m格网插值，并对插值后的DTM生成坡度图像。根据初始DTM的坡度图进行坡度分级，适当选取格网最低点，采用多级种子点方式可以提高粗DTM精度，目的是为了保留主要的地形变化特征，而CSF算法在细部地形变化区域滤波效果较好，所以不宜过分地追求粗DTM的精度。根据获得的格网最低种子点，采用线性插值法内插得到分辨率为1 m的粗DTM。

LiDAR点云进行高程归一化处理：内插后的粗DTM用于计算点云的高程归一化，可有效减少高程归一化处理后的点云锯齿化问题，且在1 m格网范围内的地形数据起伏小，对布料模拟过程影响可以忽略。高程归一化处理结果如图8所示。

布料模拟滤波：对高程归一化处理后的点云，采用布料模拟算法进行滤波处理，得到地面点及非地面序列。

地面点及非地面点提取：将高程归一化处理后的地面点与非地面点数据，与粗DTM叠加计算，得到最终分类结果（图9）。可以看出，本文算法有效地去除了平地、山地地区建筑物，地面点都得到了有效的保留。

为了说明算法的精度和可靠性，对滤波结果进行定量及定性评价。在定量分析中，I类误差为将地面点分为非地面点的误差，II类误差为将非地面点分为地面点的误差。表2中e和f表示参考数据中地面点和非地面点的个数；g和h分别是滤波后数据中地面点和非地面点的个数；n为扫描数据点总个数；b/n×100%为第I类误差；c/n×100%为第II类误差。

为了验证该方法对参数设置的敏感特性，对原始CSF滤波算法及改进后算法进行对比分析。实验参数设置如表3。其中序号0为直接使用CSF滤波算法的滤波结果；1~10为进行高程归一化处理后使用CSF滤波算法的滤波结果。从表中序号0及序号7的结果中可以明显看出，在相同参数（即地形参数：relief；使用后处理：是；格网大小：2 m；迭代次数：1 000；阈值：0.5 m）下，滤波正确率由原来算法的88.90%提高到了95.19%；表中序号2及序号7两行数据除迭代次数不同外，其他参数均相同，两种情况下的滤波正确率相同，算法时间复杂度由原来的164 s降低到51.9 s。该结果表明：迭代次数大于一定值时，滤波正确率与迭代次数无相关性。通过对比表中序号1与2、3与4、5与6、7与8的结果可知，使用CSF滤波后的处理算法，比未经CSF滤波处理的算法高3%的正确率。

改进的滤波算法仅根据点云密度设置格网大小，选择布料模拟滤波后的处理算法和合适的阈值，无需设置更多的地形参数，减少了人为参数设置的影响。

除此之外，改进算法提高了滤波的正确率，降低了错分率、漏分率，同时算法迭代次数少，缩短了计算时间。这是由于高程归一化处理后，即便在没有完全高程归一化处理的山地，由于迭代过程中的质点之间内力的相互作用，使得该滤波算法很容易将山地与建筑物区别，降低了错分率。漏分率的降低主要在山地地区相当明显，在有限的迭代次数下，保留了更完整的山地信息。

为进一步说明改进方法的有效性，对实验结果进行定性分析。图10为直接使用布料模拟滤波的实验结果。该算法采用的后处理机制仅仅对山脚区域有些改善，山顶区域信息被削平的现象仍相当明显。对点云除去地形特征后再进行布料模拟滤波，有效地克服了CSF算法滤除山地信息的缺陷，进一步提高了算法的鲁棒性和地形适应性。通过对比真实DTM数据（图11）、改进算法DTM结果（图12）及原CSF算法DTM结果（图13），结果表明：粗DTM的引入可以极大程度地提高CSF算法的精度，增加对山地信息的保留，更好地保证滤波效果。因此，改进的滤波算法在保留平面地区滤波优势的同时，计算效率也得到了明显的提高。

然而，在坡度急剧变化区域，由于此区域最低种子点提取密度过小，使得基于最低种子点内插的粗DTM并不能反映该区域基本的地形起伏特征。为了提高该区域地面点的提取精度，在最低种子点选取过程中，可手动增加最低种子点的提取密度。

针对目前机载激光雷达点云滤波算法存在的不足，利用CSF算法在平坦地区的滤波优势，提出了基于高程归一化和布料模拟的LiDAR点云滤波算法。对算法进行了实验室验证，结果及分析表明：

（1）顾及地形特征的CSF滤波算法，在地形复杂且连续变化的大范围复杂场景区域能够获得很好的滤波效果，适应性较强。

（2）本文算法在剔除地表地物信息的同时，能够有效保留山顶信息。加入高程归一化预处理一定程度上等价于对布料位置的初始化。将LiDAR点云高程转换为粗DTM对应点的高差，将山地地形转为平地地形，充分利用了CSF算法在平坦地区滤波效果好的优势。

（3）该方法一定程度上依赖于粗差的去除及用于高程归一化处理的粗DTM。因此，建立的粗DTM应当能够表现出该区域主要的地形起伏特征，同时构建粗DTM有必要剔除机载点云地表下的噪点。在粗DTM提取过程中，可尝试采用其他滤波方法（如渐进加密三角网算法）提高粗DTM的提取精度。

（4）本文算法对于坡度变化较大且地形不连续的区域需要人工设置格网大小。