在高速公路测量中，车载激光扫描作为一种先进的测量手段^[1]，在不与测量物接触、不影响交通的情况下，能够采集大面积高精度、高密度的路面、交通设施、植被等公路及附属设施的表面信息。通过采集的点云数据可以获取高精度的路面地形，为道路改扩建、路面沉降和破损分析等应用提供基础数据，而这些深入应用的前提是如何快速精确地从激光点云数据中提取出路面点云^[2]。

目前，已有不少国内外学者对激光点云数据的地面滤波方法进行了研究。在机载方面，Axelsson^[3]提出了TIN三角网加密方法处理机载点云数据。算法首先选择局部区域内的最低点作为种子点来构建初始稀疏三角网，然后通过迭代对TIN三角网进行加密，从而实现滤波；Vosselman^[4]提出了一种基于坡度变化的滤波算法，其核心思想为：当局部范围内的任意点与待判点的坡度大于给定的坡度阈值时，待判点即为非地面点，算法的关键在于选择合适的坡度阈值；Zhang等^[5]描述了一种渐进式形态学滤波算法来过滤地物点，即在运用数学形态学方法时，依据地形坡度和窗口尺寸计算相应的高差阈值。在车载方面，Sherif ^[6]提出基于点云空间特征向量进行聚类分析，通过计算法线方向和特征值进行人行道、路面、路缘石等分类；Anttoni^[7]利用强度与高程图像自动分类道路标示和路缘石点云，并通过Delaunay三角网构建道路表面；Yang Bisheng^[8]提出一种通过逐条扫描线的处理从车载激光点云中提取路面和路缘石点云的方法，主要考虑道路路坎的形状，在扫描线上设置滤波窗口的高差。

目前，车载道路点云数据处理中，研究对象一般为城市规则道路，通常包含明显的路缘石，缺少针对高速公路(无路缘石)路面提取的研究，由于城区道路与高速公路的差异，相关算法(如基于扫描线的方法)并不能适应高速公路的点云特征；同时一些算法过多地依赖先验条件，增加处理难度。虽然在机载点云数据处理中地面点云滤波研究较多，但由于两者点云的扫描距离、目标、点密度及精度不同，所以不能直接沿用机载的方法。本文针对高速公路原始点云数据，提出了一种路面初始种子点的提取方法，然后基于坡度滤波思想，通过局部格网邻域坡度对比分析提取全部路面点。

1 高速公路路面点滤波算法

车载激光扫描系统可以获取高速公路路面及两侧地物点云，路面滤波的目的在于去除路面的车辆、交通设施、植被及两侧地物等，只保留所有路面点云。滤波算法流程图如图 1所示。

图 1 滤波算法流程图 Fig. 1 Flow Chart of Filtering Algorithm

图选项

高速公路路面坡度较小，点云分布平滑。对路面进行三维格网化处理，传统方法通常是将每个格网中高程最低点作为该格网的初始地面种子点。为提高初始地面种子点的提取概率，本文提出了近似平面约束法，算法原理如下：假定单个格网内部存在路面点，那么在格网内应存在一个平滑的地表平面，则至少需要三个非共线地面点。近似平面约束模型如图 2所示，首先将单个格网分割成大小相等的四个子格网，获取每个子格网的最低点，然后从中随机抽取三个点。三点空间关系满足如下条件：①三点构成的三角形的边长大于距离阈值(1/2格网步长)；②任意两点坡度小于区域路面最大坡度阈值；③三点构成的三角形尽可能接近正三角形。当格网内不存在符合条件的三个点时，则去除格网；若符合条件，则选择其中最低点作为疑似路面点种子点。

图 2 平面约束种子点(虚拟格网俯视图) Fig. 2 Seed Point by Plane Constraint

图选项

车载移动测量系统沿着道路采集数据，由于道路两侧地物多样，当扫描到表面近似水平且面积较大的平整地物时，如低矮、整齐的植被表面，格网最低点会被错误判读为路面种子点。针对此问题，本文首先利用有序最小二乘估计窗口坡度阈值，然后运用多尺度窗口迭代方法对格网种子点进行多次迭代过滤，进一步去除种子点中的噪声点。迭代过程如下。

1) 窗口坡度阈值估计。在近似平面约束处理之后，从目标格网一定尺度窗口的格网中提取每个格网的最低点，然后通过最小二乘拟合局部空间平面，以平面的坡度作为局部近似坡度阈值^[9]。由于初始种子点可能存在粗差点，因此采用一种有序最小二乘法拟合平面，流程如下:①将窗口格网中的种子点按照高程进行排序，取小于中值的点拟合平面，设局部区域的空间平面方程如式(1)所示，利用最小二乘法计算a、b、c 的初始值，然后计算平面坡度s。

②利用式(2)，按排序顺序计算大于中值的种子点到拟合平面的距离d，将距离d与阈值σ进行比较，小于阈值，则加入点集，再次拟合平面，获取平面参数估计值；大于阈值，则去除该点及以后的所有点。

式(2)和式(3)中，d_i为点到平面的距离；s_i为坡度；l_i为邻域窗口的对角线长。③利用所有保留下来的点重新计算参数 a、b、c，即可获得最佳平面参数，以最终平面的坡度作为局部近似坡度阈值。

2) 格网种子点判定。对于未知格网A，由路面点云的连续性可以认为路面上8邻域窗口的格网数不小于3个，否则为孤立地物，需去除；若不是孤立点，由格网A的最低点和8邻域格网中的最低点构成以A为顶点的扇形三角网，如图 3所示，约束三角形边长大于1/2格网步长，且为锐角三角形，计算每一个三角面相对于水平面的坡度^{[3, 4, 10]}，同时比较三角面三个顶点的高程关系，然后取坡度排序后的中值与窗口坡度估计阈值进行比较，若大于阈值，且A点高程较高，则认为A格网内的所有点为非路面点，需去除，被去除的格网不再参与下一步滤波。

图 3 扇形三角坡度估计 Fig. 3 Slope Estimation of Fanning Triangle

图选项

3) 多尺度窗口迭代分析。如果A格网及8邻域窗口格网的最低点均为非路面点，且高程相近，如图 4所示，A格网及周围灰色区域高程均为3 m，通过8邻域大小的窗口进行一次分析只能去除A格网及周围部分高程为3 m的种子点，无法去除所有非路面种子点，需要进行多尺度迭代滤波。其流程如下：首先利用相同尺度窗口对上一次邻域分析的结果进行迭代处理，直到前后两次保留的格网数之差在阈值范围内；然后扩大迭代分析的窗口，按照上述步骤再进行相同处理，通过多尺度窗口迭代逐层将道路非路面格网种子点和空间孤立的格网点去除，尺度最大窗口小于区域内最大的非路面平整区域，由于车载系统沿道路采集，窗口大小由道路两侧低矮、平滑的植被的表面大小决定，一般区域宽度小于10 m。

图 4 非地面平整区域示意图 Fig. 4 Diagram of Non-Ground Level Area

图选项

获取路面格网种子点之后，根据路面点云的连续性，本文将局部坡度滤波的范围控制在2 m内，根据格网步长转换为实际窗口邻域大小。以8邻域窗口为例，即由格网内的待定点与8邻域格网中的路面种子点参照上文格网种子点判定方法构成扇形三角网，计算三角面的坡度，取待定点为最高点的三角面中的最大坡度值与8邻域的坡度估计阈值进行比较，在阈值范围内的为路面点。

2 实验分析

本文利用 C++实现了上述算法，采用山东科技大学和青岛秀山移动测量公司研发的Vsurs I型车载式三维空间移动测量系统采集的点云数据。选择两处比较典型的高速公路路段，实验区域点云类型包括道路、交通附属设施、植被和汽车，两段高速公路长度约为150 m，路面点位精度在5 cm左右，路面点云平均距离小于0.1 m，点云平均密度大于100点/m²，其中直道点云总数为836 501个，平均点距为0.08 m，路面的坡度小于4.5%，边坡坡度大于65%，点云最大高差为8.01 m；弯道点云总数为687 823个，平均点距为0.07 m，路面的坡度小于5%，边坡坡度大于65%，点云最大高差为7.38 m。

在建立三维格网时，为了保留路面细节，格网不宜设置较大，但格网较小会显著降低处理效率。本文建立了三维格网步长与坡度、局部路面高差的关系为：

其中，slope为坡度；Δh为局部路面高程差值(局部范围小于单车道宽度)；Δx和Δy分别是x、y方向的步长。

实验路段路面最大坡度为5%，局部(2 m范围)路面点最大高差为0.1 m，根据式(4)计算格网的对角线距离为m，确定方形格网步长为1 m×1 m；根据垂直距离对点云进行分层，依据格网索引只对最下面一层进行处理，提高效率，本文根据点云最大高差，同时保证路面点均位于最下一层格网中，设置垂直步长为2 m。

表 1为路面提取算法运行结果的统计数据，可以看出，格网内近似平面约束法和多尺度窗口 分析法能有效提高地面初始种子点的概率，同时大大减少了所需处理的格网数，提高了运算效率。以直道为例，均匀格网处理后，有7 223个格网需要进一步处理，近似平面约束使格网数 降低为5 613，窗口迭代分析使格网数进一步降低为5 074，两步运算共去除了直道2 149个非路面格网种子点。

弯道路面运用本文的算法去除非路面种子点的前后对比如图 5所示。可以看出，在原始数据中，道路两侧存在大量非地面点，经过两步处理被成功去除；两段道路的路面滤波三维视图结果(路面为真彩点云)如图 6和图 7所示。可以看出，两侧与中间植被、交通设施等非路面点被准确地滤除，路面点云完整，行车道边缘的少部分平坦的路肩地面被提取出来，不影响后期应用；图 7中弯道路面中间有护栏遮挡，由于采集系统在左侧行驶，右侧路面点云比较稀疏，在滤波时进行格网点密度约束，经过滤波只保留车载系统行驶一侧的高密度路面点云。

图 5 路面初始种子点提取过程 Fig. 5 Extraction Process of Road Surface Seed Points

图选项

图 6 直道路面滤波结果 Fig. 6 Filtering Results of Straight Road

图选项

图 7 弯道路面滤波结果 Fig. 7 Filtering Results of Curve Road

图选项

借鉴机载LiDAR数据滤波评价的误差指标^[10]，将误差error(I)定义为路面点误分为非路面点的误差，则：

其中，m表示路面点分类正确的点数；n表示路面点分类成非路面点的点数。

误差 error(II)定义为非路面点误分为路面点的误差，则：

其中，k表示非路面点分类成路面点的点数；l表示非路面点分类正确的点数。

为综合评估算法提取路面的可行性，本文对路面提取的准确率和完整率两项指标进行评估，路面准确率(pavement accuracy rate,PAR)定义为正确提取路面点面积与自动提取全部路面点面积的比值；路面完整率(pavement integrity rate,PIR)定义为正确提取路面面积与实际道路路面面积的比值。

其中，η表示提取正确路面点的面积；μ表示提取非路面点的面积；β表示实际路面面积。

按上述误差计算方法，采用人工提取的方法从直道数据中提取出误差统计值m、n、k、l和评估统计值η、μ、β，其中实际道路路面面积β为4 704.63 m²。表 2列出了本文算法与其他几种滤波算法(扫描线法^[2]、数学形态^[5]、斜率法^[9])的对比，其中第II类误差要小于以上三种方法，第I类误差稍微偏大，但是为了保持较高的数据精度，实验中主要需要控制第 II类误差；本文算法提取路面的准确率为95.74%，优于其他方法，完整率为98.11%，缺失区域主要为车辆遮挡，可见，本文提取的路面点云误差较小，准确率高。

车载激光扫描系统能快速获取高速公路路面点云数据，本文通过对高速公路点云数据的空间分布特征进行了分析，充分考虑路面点云平滑、连续分布的特点，运用近似平面约束法和窗口迭代分析有效提取了初始路面种子点；在估计局部坡度时，提出了有序最小二乘拟合法。算法结合了数据统计特性，通过中值提取增强了滤波的稳健性。实验结果和误差分析证明，本文提出的路面点滤波方法的适应性强，稳定性高，能有效提取高速公路路面高密度点云，通过对直道结果进行分析，得出I类与II类误差分别为2.8%、1.08%；准确率和完整率分别达到95.74%、98.11%。另外，本文算法需要一些人工经验作为阈值参数(如迭代次数控制阈值)，自适应性有待进一步研究。下一步将利用真彩信息进行少量路肩去除。