基于多平面分割和矩阵变换的航摄边坡点云滤波算法

An Aerial Photography Slope Point Cloud Filtering Algorithm Based on Multi-Plane Segmentation and Matrix Transformations

摘要

点云滤波是实现地面点与非地面点分离，获取最真实地面点云的重要处理手段，是进行公路边坡点云灾害识别的基础.为了克服传统点云滤波算法在边坡场景下处理速度慢、结果精度低和误差大等问题，提出了基于多平面分割和矩阵变换的航摄边坡点云滤波算法.该方法首先利用基于曲率的区域生长算法对边坡进行多平面分割，得到多个边坡子点云；其次拟合得到边坡子点云平面模型，再利用旋转矩阵将边坡子点云进行空间转换至水平面；然后通过模拟布料下沉并设置距离阈值分离出非地面点；最后再次利用旋转矩阵的逆矩阵进行空间位置还原，进而得到滤波后的边坡点云.利用精细贴近仿地航线设计方法获取高精度点云模型以在不同的边坡场景下进行算法测试，并与其他传统滤波算法结果进行对比，结果表明：在所有的试验中，本文算法的总误差分别为7.11%、4.15%、1.45%、4.41%，在所有测试算法中最小；Kappa系数分别为0.77、0.90、0.96、0.90，在所有测试算法中最大.本文提出的算法在面对地形和植被情况复杂的边坡情景下，有着较高的准确性以及较强的适用性，为公路边坡点云滤波提供了新的解决方案.

Abstract

Point cloud filtering is a crucial processing technique for separating ground points from non-ground points and obtaining the most accurate ground point cloud. It serves as the foundation for landslide identification in highway slope point clouds. To address issues such as slow processing speed， low result accuracy， and high error rates encountered by traditional point cloud filtering algorithms in slope scenarios， an airborne slope point cloud filtering algorithm based on multi-plane segmentation and matrix transformation is proposed. This method initially employs a region growing algorithm based on curvature for multi-plane segmentation of the slopes， resulting in multiple sub-point clouds of the slopes. Subsequently， it fits plane models for these sub-point clouds and uses rotation matrices to spatially transform them onto a horizontal plane. Non-ground points are separated by simulating fabric settling with a distance threshold. Finally， the inverse of the rotation matrix is applied for spatial position restoration， yielding the filtered slope point cloud. High-precision point cloud models are obtained using a fine approximation flight route design method for algorithm testing in various slope scenarios. Results are compared with those of other traditional filtering algorithms， demonstrating that the algorithm in this paper outperforms the others with total errors of 7.11%， 4.15%， 1.45%， and 4.41% in all experiments， respectively. Furthermore， the Kappa coefficient values are 0.77， 0.90， 0.96， and 0.90， all of which are the highest among all tested algorithms. The proposed algorithm exhibits high accuracy and applicability， particularly in complex slope scenarios characterized by varying terrains and vegetation cover. It offers a new solution for point cloud filtering in highway slope applications.

Graphical abstract

关键词

公路边坡 / 无人机 / 航摄点云 / 点云滤波 / 多平面分割 / 矩阵变换

公路交通是我国交通运输体系的重要组成，对国家和社会的发展产生深远影响. 在《国家公路网规划（2013年—2030年）》的指导下，我国公路建设加快推进，截至2022年，全国公路里程535.48万公里，其中，高速公路里程17.73万公里^［1］. 在公路建设取得瞩目成就的同时，公路边坡的安全问题也愈发突出. 近年来，无人机倾斜摄影技术^［2］因具有可达性高、安全性强、成本小等优点而被广泛应用于安全监测. 不少学者利用无人机倾斜摄影技术与三维重建算法获取边坡点云数据，从三维空间来解译边坡位移变化以实现边坡灾害监测^［3］. 利用点云数据实现边坡位移监测时，关键的一步数据处理是点云滤波，即滤除边坡点云中的非地面点，保留最真实的地面. 尽管国内外学者早已着手于点云滤波的研究，但多数聚焦于简单的地形地貌或针对特定情境，这些滤波方法应用于边坡点云时效果并不理想.

近年来，国内外众多学者从不同理论角度对点云滤波进行了针对性的研究. 常用的点云滤波算法主要分为六类：基于坡度的滤波算法、基于形态学的滤波算法、基于曲面拟合的滤波算法、基于不规则三角网的滤波算法、基于分割的滤波算法以及基于机器学习的滤波算法. 基于坡度的滤波算法主要是利用地面点与地面点之间的坡度较小，而地面点与非地面点之间的坡度较大这一特性将非地面点分离出来. 由于该算法仅依赖单一的坡度阈值进行判断，因此在复杂地形地区精确度较低，尤其是坡度较大的斜坡地形中，地面点容易被误判为非地面点而被剔除^［4］. Zhang等^［5］提出了一种经典的形态学的滤波算法，该算法利用地面点与非地面点进行形态学开运算之后的高程变化不同，通过设置高程变化的阈值来分离非地面点. 试验证明，该算法在陡峭地区具有较强的鲁棒性，但其参数设置依赖于地形特征，在地形复杂区域的适应性需进一步研究. 对于地形起伏较大的地形特别是边坡地区，极易出现地面点误判为非地面点的情况^［6］. 基于曲面拟合的滤波算法是通过插值方法拟合出地形曲面，并根据一定准则逐步剔除不满足条件的非地面点. Kraus等^［7］利用线性函数并设置权重，逐步滤除非地面点，在森林区域取得了良好效果，但在其他地形条件下适应性较差^［8］. 基于不规则三角网的滤波算法通过选取地面种子点建立三角网格，迭代判断并将满足条件的点加入三角网格. 尽管其可靠性和准确性较高，但存在内存占用大、计算时间长的问题，且对种子点的选择极为敏感^［9］. 基于分割的滤波算法首先对点云进行分割处理，然后采用特定方法进行滤波. Lin等^［10］使用区域增长算法对点云进行分割，并通过建立不规则三角网分离地面点，有效降低了分类误差. 然而，该方法对分类准确性要求较高，在植被密度较高的地区应用较为困难^［11］. 基于机器学习的点云滤波算法将点云视为二分类问题，利用机器学习算法进行训练和分类. Hu等^［12］采用卷积神经网络进行点云滤波，取得了低误差和高准确率，但该方法需要大量标记的点云样本，涉及大量人力、物力和时间投入^［13］. 此外，Zhang等^［14］提出的布料模拟滤波算法也是常用滤波算法之一. 该算法通过模拟布料来拟合地面地形，根据点与模拟布料的距离来分离出非地面点. 该算法能够高效地分离地面点与非地面点，但在复杂场景以及密集点云中的滤波效果并不理想^［15］.

上述算法多针对简单的大平面场景，植被密度较小、点云密度较小的场景开发，而在面对地形起伏大、植被密度高的复杂地形，特别是公路边坡时，滤波效果往往不尽如人意. 鉴于传统滤波算法在边坡场景中无局限性，研究针对性的边坡点云滤波算法显得尤为重要.

高精度滤波实现的前提是生成高质量的边坡点云模型. 目前，学者们正积极致力于优化无人机倾斜摄影的航线规划，以期提升三维重建模型的质量. 例如，余加勇等^［3］探索了采用五向飞行航线获取边坡影像用于三维重建，但该方法面临作业量大、照片冗余度高以及建模周期长等问题. 传统的航线规划在处理边坡地形时，同样存在作业成本高、飞行风险大等问题^［16］，且获取的图像可能因地面分辨率差异显著而导致信息缺失和模型精度下降. 因此，研究专门针对边坡地形的无人机倾斜摄影航线规划方法显得尤为重要.

本文首先采用精细的仿地航线设计，实现对边坡地形的无人机精准航线规划，从而获取高质量的边坡影像数据. 随后，利用三维重建算法，生成高质量的边坡点云数据，以供后续处理使用.其次，利用多平面分割和矩阵变换的航摄边坡点云滤波算法对点云进行分割.再次，利用旋转矩阵进行变换，分离出地面点后还原边坡点云. 最后，将本文提出的滤波算法应用于长沙市马桥河路边坡试验中，并与传统滤波算法的结果进行对比.

1 边坡精细贴近仿地航线设计方法

根据目标边坡附近的地理环境，在确保无人机飞行安全及高效率的前提下，设计飞行高度、飞行速度、地面分辨率、航向重叠率、旁向重叠率、作业区域等参数，提前规划一组适用于飞行作业的航线. 以下将着重阐述飞行高度和重叠率参数的计算理论. 飞行高度，亦称为航高，是指无人机进行航摄时的飞行高度. 地面分辨率（Ground Sample Distance， GSD）是航摄的影像中每一个像素所对应的地面实际尺寸，是确定航高的关键参数. 通常，在其他条件保持不变时，相对航高越低，则地面分辨率越高. 在确定了相机参数和地面分辨率要求之后，根据摄影成像原理即可得到飞行高度H：

H=f\times \frac{R_{\mathrm{G}\mathrm{S}\mathrm{D}}}{a}（1）

式中：R_GSD为地面分辨率，m/pixel；a为像元尺寸；H为飞行高度，m；f为相机焦距，mm.

航线的重叠率主要分为航向重叠率和旁向重叠率. 保证航摄影像的重叠率是实现高精度的三维模型重建的关键. 航向重叠率是指无人机在沿一条航线飞行时，第一次拍摄的照片与第二次拍摄照片重叠的概率，如图1所示. 旁向重叠率是指无人机在第一条航线和第二条航线所拍摄的照片重叠的概率，如图2所示.

精细贴近仿地航线设计主要包含两部分：交叉环绕航线和贴近仿地航线. 交叉环绕航线所采集的影像数据主要用于还原边坡的真实三维空间位置，以实现边坡三维模型的重建. 在进行交叉环绕航线规划时，无人机实际上是对目标边坡区域进行环绕倾斜拍摄，无人机的镜头俯仰角一般设置为-60°~-45°，然后结合无人机的镜头参数和地面分辨率要求确定航线的飞行高度. 确定飞行高度后，根据预设的重叠率、飞行作业的区域范围以及镜头俯仰角，确定每一个航点的具体位置.

贴近仿地航线采集的影像数据旨在提升边坡模型的精度，以实现高精度的边坡三维模型构建. 相较于传统航线，贴近仿地航线解决了在起伏地形区域使用传统航线所面临的分辨率降低和重叠度不足等问题. 如图3所示，传统航线规划保持恒定飞行高度，但随着地形的变化，图像数据的重叠率逐渐下降，最终导致建模失败. 相反，如图4所示，贴近仿地航线根据地形变化调整航点的绝对高度，确保每个航点与坡面的相对高度不变. 在边坡贴近仿地航线规划中，我们使用“井”字飞行航线规划路径，无人机向下正射获取影像. 飞行高度是根据地面分辨率、相机镜头参数和边坡高程计算确定的. 随后，根据设置的重叠率、飞行区域和镜头俯仰角，确定每个航点的位置.

2 多平面分割和矩阵变换的点云滤波

本文所提出的基于多平面分割和矩阵变换的点云滤波方法步骤如图5所示，主要包括边坡多平面分割、基于旋转矩阵的空间转换、边坡地表点云识别3个步骤.

2.1 边坡点云多平面分割

利用无人机倾斜摄影所获取的边坡图像数据生成边坡点云数据^［3］，首先基于区域增长理论^［17］根 据不同的坡度对边坡点云数据进行分割处理，得到不同坡度的边坡子点云，为点云的空间转换做好准备.

边坡多平面分割的具体步骤如下.

1）计算法线和曲率，选择种子点.依据确定的搜索范围， 针对点云中的每个点计算其曲率和法线，依据点的曲率进行升序排序，选择曲率最小的点作为种子点放入种子点序列，以确保区域从曲率最小的点开始增长，提高分割效率.

2）法线夹角与曲率判断. 设置法线夹角阈值和曲率阈值，搜索当前种子点的邻域点，计算邻域点的法线与当前种子点的法线之间的夹角，若小于法线夹角阈值的邻域点，则放入到当前平面并标记为已分类，再检查当前平面点的曲率；若小于曲率阈值，则放入种子点序列中，并删除当前种子点.

3）平面区域生长. 提取出种子点序列中的下一个种子点，继续重复步骤2中的法线夹角和曲率判断，直至种子点序列清空，则完成当前平面区域分割.

4）多平面分割. 剔除边坡点云中标记为已分类的点，将剩下的点云重复上述步骤1至步骤3，直至剩余的点云数目小于聚类个数的阈值，至此完成边坡点云的多平面分割. 利用曲率对边坡点云坡面进行分割之后如图6所示，得到边坡子点云后方便进行边坡点云的空间转换.

2.2 基于旋转矩阵的点云空间转换

对每一个边坡子点云进行空间转换，转换主要分为两个步骤.

1） 边坡子点云平面提取.首先随机选取3个点，根据平面方程A_{\mathrm{0}}x+B_{\mathrm{0}}y+C_{\mathrm{0}}z+D_{\mathrm{0}}=\mathrm{0}，计算出模型参数. 用剩余的点计算到该平面方程的距离，将距离与设定的阈值做比较，如果距离小于阈值，则该点为内点，否则为外点. 统计该参数模型下内点的个数. 继续执行模型参数计算和内点识别，若当前模型的内点数量大于已经保存的最大内点数量，则变更新的模型参数，保证模型参数始终是内点数量最多的模型参数. 重复上述步骤，不断迭代，直到达到迭代阈值，找到内点个数最多的模型参数，最后利用内点再次对模型参数进行估计，从而得到最终的边坡子点云 P_n 平面模型，如图7所示.

2）根据提取的平面的参数计算旋转矩阵，并对边坡子点云进行旋转的空间转换. 边坡子点云平面的法向量\mathit{n}=(A,B,C)，在XOZ平面的投影为\mathit{a}=(A,\mathrm{0},C)，在YOZ平面的投影为\mathit{b}=(\mathrm{0},B,C). 平面XOZ法向量\mathit{p}=(\mathrm{0,1},\mathrm{0})，平面YOZ法向量\mathit{q}=(\mathrm{1,0},\mathrm{0}). 计算得到向量\mathit{a}与平面YOZ的夹角，如式（2）所示.

\theta =\frac{\mathrm{3}\mathrm{\pi }}{\mathrm{2}}+\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\mathit{a}\times \mathit{q}}{|\mathit{a}|\times |\mathit{q}|}\right)（2）

则边坡子点云 P_n 绕Y轴逆时针旋转的旋转矩阵为：

{\mathit{R}}_{ny}=\left[\begin{array}{cccc}\mathrm{c}\mathrm{o}\mathrm{s}\theta & \mathrm{0}& \mathrm{s}\mathrm{i}\mathrm{n}\theta & \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ -\mathrm{s}\mathrm{i}\mathrm{n}\theta & \mathrm{0}& \mathrm{c}\mathrm{o}\mathrm{s}\theta & \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]（3）

向量\mathit{b}与平面XOZ的夹角φ为：

\phi =\frac{\mathrm{\pi }}{\mathrm{2}}-\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\mathit{b}\times \mathit{p}}{|\mathit{b}|\times |\mathit{p}|}\right)（4）

则边坡子点云 P_n 绕X轴逆时针旋转的旋转矩阵为：

{\mathit{R}}_{nx}=\left[\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{c}\mathrm{o}\mathrm{s}\phi & -\mathrm{s}\mathrm{i}\mathrm{n}\phi & \mathrm{0}\\ \mathrm{0}& \mathrm{s}\mathrm{i}\mathrm{n}\phi & \mathrm{c}\mathrm{o}\mathrm{s}\phi & \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]（5）

通过计算得到旋转后与平面XOY平行的边坡子点云如式（6）所示.

{\mathit{P}}_n^{\mathrm{\text{'}}}={\mathit{R}}_{ny}{\mathit{R}}_{nx}{\mathit{P}}_n（6）

2.3 边坡地表点云识别

经过多平面分割和空间转换后，边坡点云场景得到了简化，布料模拟算法在处理大平面场景时具有处理效果好、速度快的特点，因此利用模拟布料算法识别简化后的边坡地表点云. 如图8所示，利用坐标转换，使边坡子点云{\mathit{P}}_n^{\mathrm{\text{'}}}沿着水平面XOY进行上下翻转；基于布料模拟理论^［14］，设置初始布料网格的大小，通常将布料的位置设置在最高点的上方；将所 有的边坡子点云{\mathit{P}}_n^{\mathrm{\text{'}}}的点与布料网格的节点投影至同一个水平面，保留其投影之前的高程信息，构建Kd-tree以找到每一个布料网格点在边坡子点云{\mathit{P}}_n^{\mathrm{\text{'}}}中的最邻近点；计算所有的可移动的布料网格点受重力作用产生的位移，并与当前布料网格点对应边坡子点云{\mathit{P}}_n^{\mathrm{\text{'}}}中最邻近点的高程信息进行比较，如果当前布料网格点的高程低于或者等于最邻近点的高程，则把当前布料网格点的高程设置为最邻近点的高程，并将其设为不可移动点；计算每个可移动的布料网格点受到其邻近节点内力作用所产生的位移；重复上述过程，直到最大高程变化足够小或者达到所预设的迭代次数；计算布料网格点与边坡子点云{\mathit{P}}_n^{\mathrm{\text{'}}}之间的高程差，将高程差小于设定高度阈值的点云判断为地面点，反之则判断为非地面点. 通过布料模拟分离出边坡子点云{\mathit{P}}_n^{\mathrm{\text{'}}}的地面点{\mathit{G}}_n^{\mathrm{\text{'}}}.

利用旋转矩阵还原边坡子点云. 将分离出的边坡子点云{\mathit{P}}_n^{\mathrm{\text{'}}}的地面点{\mathit{G}}_n^{\mathrm{\text{'}}}通过旋转矩阵还原至初始位置，得到原始边坡地面子点云公式为：

{\mathit{G}}_n={\mathit{R}}_{nx}^{-\mathrm{1}}{\mathit{R}}_{ny}^{-\mathrm{1}}{\mathit{G}}_n^{\mathrm{\text{'}}}（7）

最后，将所有获得的边坡地面子点云G_n^{}进行拼接，得到最后滤波的结果，即边坡地面点云合集，如式（8）所示.

G=\{{\mathit{G}}_{\mathrm{1}}^{},{\mathit{G}}_{\mathrm{2}}^{},{\mathit{G}}_{\mathrm{3}}^{},\cdots ,{\mathit{G}}_{m-\mathrm{1}}^{},{\mathit{G}}_m^{}\}（8）

3 边坡试验

3.1 试验区域与仪器设备

为了验证基于多平面分割和矩阵变换的边坡点云滤波算法的有效性，开展了边坡测试工作，所获取的数据为后文的研究奠定了基础. 如图9所示，本次试验数据的采集对象为长沙市望城区马桥河路附近的两处边坡（北纬28.271 5°，东经112.845 6°）. 两处边坡分别命名为马桥河路1号边坡、马桥河路3号边坡. 马桥河路1号边坡长约232 m，高约32 m，坡度约为43.58°，分类为陡坡和超高边坡；马桥河路3号边坡长约69 m，高约13 m，坡度约为27.21°，分类为中等坡和高边坡. 为了验证不同植被环境下滤波算法的有效性，分别采集了两期数据，同一边坡两期数据的植被情况差异较大，第一期采集日期为2023年4月1日，第二期采集日期为2023年7月8日.

本次边坡数据采集设备为大疆M210-RTK无人机和大疆精灵4 RTK， M210-RTK无人机可实时接收全球定位系统和北斗导航卫星系统信号，带有实时动态测量（Real Time Kinematic，RTK）模块，可实现厘米级定位精度. 搭载云台相机为禅思X5s，采用M4/3传感器，像素高达2 080万，本试验使用的镜头为 15 mm.精灵4 RTK使用网络RTK系统，具备厘米级导航定位系统和高性能成像系统，相机型号为FC6310R，像素高达2 000万.

3.2 高清边坡影像获取

为了高效获取高质量的边坡影像，此次研究利用精细贴近仿地航线规划方法获取影像. 将交叉环绕航线的各个航点位置和执行动作生成KML文件导入飞行器遥控器上，仔细检查边坡周围环境和航点位置，确认安全后通过遥控器将文件传输至无人机系统，无人机开始执行航摄任务，采集边坡高清影像. 将获取的边坡影像导入三维重建软件生成两个边坡的数字表面模型文件（Digital Surface Model， DSM），再将DSM文件导入航线规划软件，根据边坡的地形确定贴近仿地航线的各个航点位置以及飞行路线，出于安全考虑本次试验保持无人机始终在边坡正上方20 m处执行贴近摄影，预期GSD为0.6 cm/pixel. 生成的精细贴近仿地航线如图10所示，具体的飞行参数如表1、表2所示.

3.3 高精度三维点云模型重建

将所有获取的边坡影像导入三维建模软件中进行三维重建. 首先，对所有的边坡影像数据进行空三解算，使用边坡区域的KML文件来确定生产的三维模型范围，减少冗余数据的产生，提高建模效率，再利用运动恢复结构算法以及密集匹配算法计算获取边坡的高精度三维点云模型.

4 试验数据处理与分析

4.1 试验数据预处理

通过倾斜摄影和实景三维重建获取的边坡点云冗余数据较多，存在噪声点，直接用于边坡点云滤波试验会出现计算量大、耗时长以及处理结果不准确等问题. 因此，需要对获取的原始边坡三维点云数据进行降采样、噪点剔除等数据预处理操作后，才能将其用于边坡点云滤波试验.

常用的点云降采样方法有空间降采样、体素降采样、均匀降采样、随机降采样等. 本文根据边坡点云滤波试验要求选择空间降采样方法，在保留边坡形态特征的同时降低点云密度. 对试验使用的马桥河路1号边坡两期点云数据采用空间5 cm间隔降采样处理，即在降采样后的点云数据中，任意两点的距离都大于5 cm，其中马桥河路1号边坡一期点云数由71 774 682减少至4 665 195，马桥河路1号边坡二期点云数由42 388 111减少至2 730 459. 对试验使用的马桥河路3号边坡两期点云数据采用空间5 cm间隔降采样处理，其中马桥河路3号边坡一期点云数由28 001 360减少至1 683 701，马桥河路3号边坡二期点云数由20 969 774减少至1 121 968. 对比原始边坡点云，降采样后的点云数据在保留了边坡形态特征的同时减少了点云数据量.

通过倾斜摄影测量建模获取的三维点云数据，在生成过程中可能会因传感器的测量误差、环境的干扰等因素而包含一些无规则、无意义的噪点，这些噪点会直接影响边坡点云数据的质量和后续处理效果. 本文采用统计滤波算法对两期马桥河路边坡进行去噪处理，即对边坡点云中的每个点，计算其周围一定半径内的点的平均值和标准差，若该点到平均值的距离大于一定倍数的标准差，则判定该点为噪点，将其剔除. 处理结果如图11所示，统计滤波算法能有效地去除边坡点云中的噪点.

4.2 边坡点云滤波

在完成边坡点云数据预处理的基础上，将处理后的边坡点云数据运用本文提出的多平面分割和矩阵变换的点云滤波方法进行处理. 本文提出的方法需要设置的参数有最小聚类点数S_min、最大聚类点数S_max、曲率阈值c、网格分辨率G、迭代次数n、高度阈值H.

由于马桥河路1号边坡为多阶梯式边坡，边坡范围大，区域植被较少且多为矮小植物，同时二期植被较一期更多更高，为了把地面点与非地面点精确分开，经过多次试验，马桥河路1号边坡一期数据处理参数设置为：S_min=500 000、S_max=3 000 000、c=1、G=0.1 m、n=500、H=0.2 m；马桥河路1号边坡二期数据处理参数设置为：S_min=500 000、S_max=1 000 000、c=1、 G=0.2 m、n=500、H=0.3 m. 由于马桥河路3号边坡范围较小，植被多为低矮植物但植被数量较多，没有呈现明显的阶梯状，同时二期植被较一期更加密集高大.马桥河路3号边坡一期数据处理参数设置为：S_min=200 000、S_max=1 000 000、c=1、G=0.1 m、n=1 000、H=0.1 m；马桥河路3号边坡二期数据处理参数设置为：S_min=200 000、S_max=1 000 000、c=1、G=0.2 m、n=1 000、H=0.3 m. 马桥河路3号边坡二期数据处理结果如 图12所示，可以看出本文算法能有效地将研究区域内地面点与非地面点完全分离.

除了使用本文提出的基于多平面分割和矩阵变换的点云滤波方法进行处理，还引入了几种常用的点云滤波算法进行滤波处理以进行对比分析. 这些算法具体包括：布料模拟滤波算法（Cloth Simulation Filter， CSF）、渐进形态学滤波算法（Progressive Morphological Filter， PMF）^［5］、基于坡度滤波算法^［18］、移动曲面滤波算法（Hierarchical Moving Curved Fitting Algorithm， HMCFA）^［19］.

4.3 结果定量分析

交叉表格法^［20］是国际摄影测量与遥感学会（ISPRS）于2003年提出的一种常用的评价点云滤波精度的方法. 为了定量评价本文提出的点云滤波算法的精度和可靠性，采用人工识别的方式提取点云地面点作为参考数据，并应用交叉表法来评价点云滤波精度，从而实现对滤波效果好坏的量化评估. 如表3所示，评价体系将误差分为Ⅰ类和Ⅱ类误差. Ⅰ类误差指的是将地面点误识别为非地面点的误差，Ⅱ类误差指的是将非地面点误识别为地面点的误差. 计算两类误差的公式如式（9）、式（10）所示，其中T_g为地面点正确分类的数量；F_g为地面点误判为非地面点的数量；F_u为非地面点误判为地面点的数量；T_u为非地面点正确识别的数量；R_G和R_U分别表示真实地面点和非地面点数量；C_G和C_U分别表示滤波算法识别的地面点和非地面点的数量；N为点云总数. 总误差的计算如式（11）所示.

除了利用交叉表法评价分类精度，还可根据Kappa系数（C_Kappa）进行定量分析，Kappa系数计算见式（12）. Kappa系数一般位于0~1范围内，当Kappa系数位于0~0.20时，认为分类精度极低；位于0.21~0.40时，认为分类精度较差；位于0.41~0.60时，认为分类精度中等；位于0.61~0.80时，认为分类精度较高；位于0.81~1.00时，认为分类精度极高^［21］.

\mathrm{Ⅰ}\mathrm{类误}\mathrm{差}{E}_{\mathrm{Ⅰ}}=\frac{F_{\mathrm{g}}}{R_{\mathrm{G}}}（9）

\mathrm{Ⅱ}\mathrm{类误}\mathrm{差}{E}_{\mathrm{Ⅱ}}=\frac{F_{\mathrm{U}}}{R_{\mathrm{U}}}（10）

\mathrm{总误}\mathrm{差}{E}_{\mathrm{总}}=\frac{F_{\mathrm{u}}+F_{\mathrm{g}}}{N}（11）

C_{\mathrm{K}\mathrm{a}\mathrm{p}\mathrm{p}\mathrm{a}}=\frac{N\times \left(T_{\mathrm{g}}+T_{\mathrm{u}}\right)-\left(C_{\mathrm{G}}\times R_{\mathrm{G}}+C_{\mathrm{U}}\times R_{\mathrm{U}}\right)}{N^{\mathrm{2}}-\left(C_{\mathrm{G}}\times R_{\mathrm{G}}+C_{\mathrm{U}}\times R_{\mathrm{U}}\right)}（12）

根据交叉表精度评价模型和Kappa系数，计算得到不同算法试验结果的误差统计和Kappa系数值如图13所示.

如图13（a）所示，马桥河路1号边坡Ⅰ期数据处理结果中，本文算法的Ⅰ类误差为7.56%，Ⅱ类误差为4.86%，总误差为7.11%，尽管其Ⅰ类误差略高于其他几种滤波算法，但Ⅱ类误差远低于其他几种算法，且总误差在所有算法中最小. 此外，Kappa系数为0.77，是所有处理算法中最大的.

如图13（b）所示，马桥河路3号边坡Ⅰ期数据处理结果中，本文算法的Ⅰ类误差为3.50%，Ⅱ类误差为5.92%，总误差为4.15%，均小于其他几种对比算法，Kappa系数更是达到了0.90，大于其他所有算法.

如图13（c）所示，马桥河路1号边坡Ⅱ期数据处理结果中，本文算法的Ⅰ类误差为0.65%，Ⅱ类误差为3.40%，总误差为1.45%，3类误差均远低于其他算法处理结果，Kappa系数达到0.96，远大于其他算法，说明本文算法分类效果极佳.

如图13（d）所示，马桥河路3号边坡Ⅱ期数据处理结果中，本文算法的Ⅰ类误差为1.45%，Ⅱ类误差为9.38%，总误差为4.41%，3类误差在所有算法中均最小，Kappa系数为0.90，在所有处理算法中最大.

本文提出的滤波算法在所有试验中表现出色，相较于其他算法，其总误差最小且Kappa系数最大. 特别是在坡度更大的马桥河路1号边坡试验中，本文算法相较于其他算法在滤波精度上有显著提升. 同时，本文算法受植被环境复杂度的变化影响较小. 因此，本文基于平面分割和矩阵变换的点云滤波算法，面对不同植被环境及不同等级的边坡具有较强的适用性和稳定性，能够精确地将边坡区域内的非地面点剔除.

4.4 结果定性分析

为了进一步说明本文算法的有效性，对试验结果进行定性分析. 选取边坡试验结果中典型地形断面进行多种滤波算法的效果对比，处理结果如图14所示，其中边坡蓝色条纹处为选取断面位置.

图14（a）为马桥河路1号边坡Ⅰ期的试验结果，HMCFA、基于坡度滤波算法以及PMF滤波算法在此断面几乎无法滤除植被，同时HMCFA在边坡顶部还存在错误识别地面点为非地面点的情况；CSF算法能够有效滤除边坡底部的植被，而难以做到有效滤除边坡变坡处的植被；本文算法采用多平面分割和空间转换机制，能够有效克服边坡高差变换对于滤波算法的影响，能够最大限度地滤除植被，从而获取边坡的真实地形.

图14（b）为马桥河路1号边坡Ⅱ期的试验结果，CSF、PMF以及基于坡度滤波算法均无法识别边坡底部以及两个变坡处的植被，滤波效果欠佳；虽然HMCFA算法能识别出边坡底部和顶部变坡处的植被，但滤波效果不理想；本文算法虽然在边坡中部变坡处存在部分植被与地面连接处误判现象，但总体来说能够准确识别出此断面上的所有植被，滤波效果相对较好.

图14（c）为马桥河路3号边坡Ⅰ期的试验结果，HMCFA、基于坡度滤波算法、CSF以及PMF滤波算法均能识别边坡底部植被，但这些算法在识别植被与地面连接处时均存在不同程度的误判；在边坡中部，4种对比的滤波算法均有不同程度的将地面点误判为植被的现象；在边坡变坡处，4种对比的算法均未能识别出植被，同时在边坡顶部HMCFA同样出现将地面点误判为植被的现象；本文算法在底部植被和地面连接处也存在小部分误判现象，但底部植被识别效果优于其他所有对比算法.总体而言，本文算法在植被识别方面表现最佳，能够获取最接近真实地面的信息.

图14（d）为马桥河路3号边坡Ⅱ期的试验结果，CSF、PMF、基于坡度的滤波算法以及HMCFA均能识别到此断面的三处植被，但都存在不同程度的滤除不彻底现象，特别是在植被与地面连接处，其中HMCFA算法在变坡顶部还存在将地面点误判为植被的现象；本文算法能有效地识别三处植被并彻底滤除，处理结果与真实地面高度一致.

综上，与其他传统算法相比，本文提出的算法在边坡变坡处、非地面点与地面点连接处以及植被密集处具有更低的误识别率. 在不同情形下，该算法都能够有效滤除植被，还原最真实的边坡地形，是一种在边坡地形中具有较强适用性和较高准确性的点云滤波算法.

5 结 论

本文提出了基于多平面分割和矩阵变换的航摄边坡点云滤波算法，通过试验验证了算法在边坡复杂情形下的精确度和有效性，得到如下主要结论.

1）构建了包括交叉环绕航线、贴近仿地航线的边坡工程无人机精细航线规划方法，能够采集到高清边坡影像，进而实现了高精度边坡点云数据的获取.

2）提出了包括边坡多平面分割、旋转矩阵空间转换、边坡地表点云识别3个步骤的边坡点云滤波方法，实现了边坡点云高效滤波降噪，显著降低了传统算法在地面点与非地面点连接处的高误判率.

3）建立了统计空间降采样和统计滤波降噪算法，有效降低了点云的数据量以及噪声点对算法的干扰，为提升后续点云滤波处理的效率和精度提供了有力支持.

4）本文算法能够准确地提取地面点，在马桥河路两期的两个边坡试验中，分类结果的总误差分别为7.11%、4.15%、1.45%、4.41%，Kappa系数分别为0.77、0.90、0.96、0.90. 与其他传统滤波算法相比，总误差最小，Kappa系数最高，提取地面点的精度最高.

参考文献

[1] 交通运输部. 2022年交通运输行业发展统计公报［N/OL］. 中国交通报，（2023-06-16）［2023-10-08］.

[2] 陈仁朋，吴熠文，余加勇，等 ．基于无人机影像序列的城市精细化三维模型精度评估［J］．湖南大学学报（自然科学版），2019，46（11）：172-180．

[3] CHEN R P， WU Y W， YU J Y，et al ．Method accuracy evaluations of building urban detailed 3D model based on the unmanned aerial vehicle image sequences and its accuracy evaluatios［J］．Journal of Hunan University （Natural Sciences），2019，46（11）：172-180．（in Chinese）

[4] 余加勇，薛现凯，陈昌富，等 ．基于无人机倾斜摄影的公路边坡三维重建与灾害识别方法［J］．中国公路学报，2022，35（4）：77-86．

[5] YU J Y， XUE X K， CHEN C F，et al. Three-dimensional reconstruction and disaster identification of highway slope using unmanned aerial vehicle-based oblique photography technique［J］．China Journal of Highway and Transport， 2022， 35（4）： 77-86．（in Chinese）

[6] CHEN C F， GUO J J， WU H M，et al ．Performance comparison of filtering algorithms for high-density airborne LiDAR point clouds over complex LandScapes［J］．Remote Sensing，2021，13（14）：2663．

[7] ZHANG K Q， CHEN S C， WHITMAN D，et al ．A progressive morphological filter for removing nonground measurements from airborne LIDAR data［J］．IEEE Transactions on Geoscience and Remote Sensing， 2003， 41（4）： 872-882．

[8] TAN Y M， WANG S， XU B，et al ．An improved progressive morphological filter for UAV-based photogrammetric point clouds in river bank monitoring［J］．ISPRS Journal of Photogrammetry and Remote Sensing，2018，146：421-429．

[9] KRAUS K， PFEIFER N ．Determination of terrain models in wooded areas with airborne laser scanner data［J］．ISPRS Journal of Photogrammetry and Remote Sensing，1998，53（4）：193-203．

[10] LI F， ZHU H H， LUO Z W，et al ．An adaptive surface interpolation filter using cloth simulation and relief amplitude for airborne laser scanning data［J］．Remote Sensing，2021，13（15）：2938．

[11] NAHR S T， SAADATSERESHT M ．Utility-pole detection based on interwoven column generation from terrestrial mobile laser scanner data［J］. The Photogremmetric Record， 2021，36（176）：402-424．

[12] LIN X G， ZHANG J X ．Segmentation-based filtering of airborne LiDAR point clouds by progressive densification of terrain segments［J］．Remote Sensing，2014，6（2）：1294-1326．

[13] 张刚，刘文彬，张男 ．基于区域特征分割的密集匹配点云渐进形态学滤波［J］．地球信息科学学报，2019，21（4）：615-622．

[14] ZHANG G， LIU W B， ZHANG N ．Progressive morphological filtering method of dense matching point cloud based on region feature segmentation［J］．Journal of Geo-Information Science，2019，21（4）：615-622．（in Chinese）

[15] HU X Y， YUAN Y. Deep-learning-based classification for DTM extraction from ALS point cloud［J］. Remote Sensing， 2016，8（9）： 730．

[16] LI B W， LU H， WANG H，et al ．Terrain-net： a highly-efficient，parameter-free，and easy-to-use deep neural network for ground filtering of UAV LiDAR data in forested environments［J］．Remote Sensing， 2022， 14（22）： 5798．

[17] ZHANG W M， QI J B， WAN P，et al ．An easy-to-use airborne LiDAR data filtering method based on cloth simulation［J］．Remote Sensing，2016，8（6）：501．

[18] ŠTRONER M， URBAN R， LIDMILA M，et al ．Vegetation filtering of a steep rugged terrain：the performance of standard algorithms and a newly proposed workflow on an example of a railway ledge［J］．Remote Sensing，2021，13（15）：3050．

[19] 王涛， 阳利永， 曾举， 等 ．不同航线规划模式下单体建筑物三维建模对比研究［J］．测绘通报， 2019（）： 298-302．

[20] WANG T， YANG L Y， ZENG J，et al ．Comparative study on 3D modeling of single buildings in different route planning modes［J］．Bulletin of Surveying and Mapping，2019（Sup.1）： 298-302．（in Chinese）

[21] 李仁忠， 刘阳阳， 杨曼， 等 ．基于改进的区域生长三维点云分割［J］．激光与光电子学进展， 2018， 55（5）： 325-331．

[22] LI R Z， LIU Y Y， YANG M，et al ．Three-dimensional point cloud segmentation algorithm based on improved region growing［J］．Laser & Optoelectronics Progress，2018，55（5）：325-331．（in Chinese）

[23] VOSSELMAN G ．Slope based filtering of laser altimetry data［J］．International Archives of Photogrammetry & Remote Sensing，2000， 33： 935-942．

[24] WEI S U， SUN Z P， ZHAO D L，et al. Hierarchical moving curved fitting filtering method based on LIDAR data［J］．International Journal of Remote Sensing，2009，13（5）： 827-839．

[25] SITHOLE G， VOSSELMAN G. Experimental comparison of filter algorithms for bare-Earth extraction from airborne laser scanning point clouds［J］. ISPRS Journal of Photogrammetry and Remote Sensing， 2004， 59（1/2）： 85-101．

[26] TANG W， HU J， ZHANG H，et al. Kappa coeffcient： a popular measure of rater agreement［J］. Shanghai Archives of Psychiatry，2015， 27（1）： 62-67．

基金资助

湖南省水利科技项目（XSKJ2021000-46）; 江苏省水利科技项目（2021074）