LiDAR测量技术现在广泛用于三维数字城市建设^[1–3]、道路安全设施检测^[4,5]、公路勘察设计等方面^[6,7], 车载LiDAR测量具有数据采集的高精度、高效率和高密度等特性^[8,9], 但是车载LiDAR移动速度较慢, 影响道路使用, 而且当路况不好的时候, 车辆无法驶入, 对于特殊状况想要对道路进行检测则力不从心; 机载LiDAR测量则是有快速度、高效率, 受天气、路况影响小, 数据生产周期短等特性, 该测量方式可分为高空扫描和低空扫描, 相比于高空扫描, 低空扫描点云数据间隙小, 精度更高^[10], 更适合用于对道路的养护. 因此本文使用低空飞行的无人机对特定的道路进行扫描.

道路数据的提取方法, 文献[9]中提出了使用曲线拟合道路数据的算法, 但是只将路面的坡度走势进行拟合, 并不能提高路面的提取精度, 文献[11]提出以点云的区域分布特性为基础, 通过对数据进行主成分分析(PCA), 构建点云平面基元检测的新模型, 虽然能将平面进行分割, 但是由于道路细节近似拱形, 所以效果并不理想, 文献[12]使用种子区域生长算法提取道路, 但是需要手动确定种子点, 而且提取精度较粗, 不能将路面的坑槽有效提取.

综上所述, 通过对上述算法的汲取和改进, 主要创新部分在于使用了聚类的算法来确定道路数据的种子点, 改进了半自动确定道路的方法, 并且用动态拟合的方式提取道路.

本文使用基于无人机的道路点云分割提取算法, 主要内涵在于无人机低空扫描提高数据精度和动态拟合分割提高提取精度, 其整体流程如图1所示, 使用无人机低空扫描路面获取精度较高的点云数据, 计算点云数据的法向量, 丰富数据信息, 为后续的聚类算法提供多维样本数据, 之后使用聚类算法根据点云数据的高程信息和法向量信息确定种子数据, 再对断面进行拟合, 利用数据动态拟合分割提取路面数据和非路面数据, 最后对非路面数据使用区域生长算法进行分割.

1 丰富点云数据信息

本文采用无人机装载激光扫描设备的方法, 无人机低空沿道路方向飞行, 扫描指定路段, 为方便后续说明, 记y轴方向为道路方向, x轴方向为道路横截面方向, z轴为道路高程信息, 道路三维点云数据作为集合P.

后续聚类算法需要使用点云的法向量信息, 就要确定其法向量, 本文使用主成分分析法(Principal Component Analysis, PCA)^[13–17], 用来计算点云数据的法向量. 对于道路的点云数据集合P, 由于它是无序的, 所以在对集合P进行处理之前, 需要为其建立拓扑关系, 本文使用K-D树(K-Dimensional树的简称)^[18–20]数据结构建立点云之间的拓扑关系, 之后再用PCA计算法向量. PCA计算方法如下:

设集合P中的点p_i的k邻域点集为N_k(p_i), 构造关于点p_i的3×3协方差矩阵如式(1)所示:

$\begin{array}{*{20}{c}} {C = {{\left[ {\begin{array}{*{20}{c}} {{p_1} - \overline p }\\ \vdots \\ {{p_k} - \overline p } \end{array}} \right]}^{\rm{T}}}\left[ {\begin{array}{*{20}{c}} {{p_1} - \overline p }\\ \vdots \\ {{p_k} - \overline p } \end{array}} \right],}&{{\text{其中}}\;\overline p = \dfrac{1}{k}\displaystyle\sum\limits_{{p_j} \in {N_k}\left( {{p_i}} \right)} {{p_j}} } \end{array}$

(1)

矩阵C表示了点集N_k(p_i)的点的分布情况, 即邻域点与质心点p的偏离程度, p_k为k邻域点集中的点, p为数据点p_i及其k邻域所构成集合的质心.

矩阵C为对称半正定矩阵, 有3个非负实数特征值, 设为λ₀, λ₁, λ₂, 且λ₀≤λ₁≤λ₂, 其对应的特征向量分别为v₀, v₁, v₂. 最小特征值λ₀对应的特征向量v₀即可近似为点p_i的法向量n_i, 为方便表示, 记单位法向量n_i= (x, y, z).

2 路面及资产分割提取 2.1 多滑动窗口的均值漂移聚类算法确定种子数据

利用聚类算法提取道路信息, 存在3个问题, 一是由于激光在扫描的时候存在噪声, 所以并不是所有的点的法向量都接近聚类中心, 结果只能提取道路上的部分数据; 二是如果对所有断面都采用聚类算法, 虽然可以提取道路, 但是非常耗时; 三是针对路面上高程信息不明显的数据, 同样需要提取, 而聚类算法并不能提取出来.

根据上述问题, 本文结合数据高程信息和法向量作为聚类算法的样本点, 多滑动窗口的均值漂移聚类算法^[21–23]确定拟合道路的种子数据, 之后利用动态拟合提取所有路面.

因为高程信息与法向量的信息比例不同, 所以需要将高程数据范围限定在0到1, 则使用高程数据在高程区间的相对值h做样本空间的高程信息, h值的计算公式如式(2)所示:

$ {{h_i} = \frac{{{{\rm{z}}_i} - \min (Z)}}{{\max (Z) - \min (Z)}}}, \;\;\;{ {{h_i} \in H}} $

(2)

其中, z_i为当前点的高程信息, Z为所有数据点的高程信息集合. 由式(1)得到数据点p_i的法向量n_i=(x, y, z), 将这两类数据投到特征空间, 则点p_i在特征空间的向量为q_i=(x, y, z, h).

得到标准化后的有n个样本点的特征空间Q=(q₁, q₂, …, q_n)后, 使用多滑动窗口的均值漂移算法确定中心点, 在特征空间Q中随机选取10个点作为初始中心点, 其滑动窗口半径为r, 半径r的取值越大则对路面数据的要求越宽泛, 越小则越严格, 本文依据数据点到距离其最近的10个点的距离来确定r值, 公式如式(3):

$r = \frac{1}{n}\sum\limits_{i = 1}^n {\max ({S_i})} $

(3)

其中, S_i为数据点q_i的最近邻10个样本点的欧氏距离的集合.

由式(3)得到高维区域的欧氏距离半径, 由高维区域内的所有样本点可以确定各个窗口的偏移均值, 其计算公式如式(4)所示:

$\begin{array}{*{20}{c}} {M({c_i}) = \dfrac{1}{{{m_i}}}\displaystyle\sum\limits_{{q_j} \in {S_r}} {({c_i} - {q_j})} },&{i = 1,2,\cdots,10} \end{array}$

(4)

其中, c_i为窗口中心, Sr为以c_i为中心, 半径为r的高维区域, q_j为区域内的点, m_i为区域内q_j的数量, M(c_i)为以c_i为窗口中心的偏移均值.

根据式对10个滑动窗口进行中心更新, 更新公式如式(5)所示:

$c_i^{t + 1} = M_i^{t} + c_i^t$

(5)

其中, $M_i^t $ 为t状态下求得的偏移均值, $c_i^t $ 为t状态下的中心, $c_i^{t+1} $ 为t+1状态下的中心.

当中心点不再变化时, 停止更新滑动窗口的中心, 则高维区域内点数最多的两个窗口内的点分别为人行道和机动车道的数据点. 通常, 人行道的高度要高于机动车路面, 因此两个密度中心可以用高程信息来区分, 高程信息低的是机动车路面, 高的是人行道路面.

2.2 对路面数据动态拟合

通过对道路数据点的观察并结合实际, 机动车路面一般是中间稍高于两侧, 这样利于排除雨水. 同样, 两边的人行道, 也存在一定弧度. 可以将路面粗略描述为光滑曲面, 针对这种形状, 采用四次多项式拟合路面横断面曲线, 使用最小二乘法计算多项式系数.

拟合函数定义如下:

$f(x) = {a_0} + {a_1}x + {a_2}{x^2} + {a_3}{x^3} + {a_4}{x^4}$

(6)

其中, a为多项式系数, x为路面横截横坐标. 由式(6)得到均方误差, 均方误差如式(7)所示:

$\begin{aligned} Q({a_0},&{a_1},{a_2},{a_3},{a_4}) \\ & = \sum\limits_{i = 1}^n \left({a_0} + {a_1}{x_i} + {a_2}x_i^2 + {a_3}x_i^3 + {a_4}x_i^4 - {y_i} \right)^2 \\ \end{aligned} $

(7)

为求式(7)的极小值, 得到方程如式(8)所示:

${A^{\rm{T}}}A\left( {\begin{array}{*{20}{c}} {{a_0}}\\ {{a_1}}\\ \vdots \\ {{a_4}} \end{array}} \right) = {A^{\rm{T}}}\left( {\begin{array}{*{20}{c}} {{y_1}}\\ {{y_2}}\\ \vdots \\ {{y_n}} \end{array}} \right),{\text{其中}}{{A = }}\left( {\begin{array}{*{20}{c}} 1&{{x_1}}& \cdots &{x_1^4}\\ 1&{{x_2}}& \cdots &{x_2^4}\\ \vdots & \vdots & \ddots & \vdots \\ 1&{{x_n}}& \cdots &{x_n^4} \end{array}} \right)$

(8)

解方程得到a₀, a₁, a₂, a₃, a₄带入式(6)可以得到当前断面的拟合函数f (x), 用当前断面的拟合函数f (x)判断下一断面的点是否符合式(9).

$\left| {f(x) - y} \right| < a$

(9)

符合式(9)的则为道路数据点, 不符合的则为非道路数据点.

$a = 2d + \frac{v}{n} \times 11\% $

(10)

$d = \frac{1}{n}\sum\limits_{i = 1}^n {\max ({D_i})} $

(11)

其中, a的计算公式如式(10)所示, 因为需要考虑道路的坡度问题, 根据《城市道路工程设计规范》内容, 坡度最大不会超过11%, 所以对提取范围修正. d由式(11)计算求得, D_i为数据点p_i的最近邻10个点距离的集合, v为无人机的移动速度, n为激光扫描仪每秒扫描的次数, 11%为国家规定公路最大坡度.

所谓动态拟合, 是将所有符合式的道路数据点进行如式拟合, 得到新的拟合函数f (x), 对下一断面进行重复操作, 达到提取所有道路数据点的目的.

在利用拟合曲线提取机动车道和人行道路面数据时, 为避免使用拟合曲线时将周围环境信息误提, 需要对拟合曲线的区间进行约束.

通常情况下可以利用上一层断面的数据确定机动车路面和人行道路面的范围, 但是, 在有遮挡物的情况下计算的范围会受到遮挡物影响, 而且不止会使当前路面的提取受到影响, 还会波及后面的计算, 所以, 准确的计算出路面的范围是非常有必要的. 针对这种情况, 提出基于宽度的路面范围修正算法.

根据现实考虑, 因为左右两边人行道同时边沿被遮挡的几率很小, 所以可以利用左右两边人行道的宽度比修正左右两边, 利用靠墙点作为基准对宽度进行修正, 设计修正算法如下:

算法1. 路面区间修正算法

1) 记前一断面左边人行道数据集合Q_{left_old}, 右边人行道数据集合Q_{right_old};

2) 根据集合Q_{left_old}, Q_{right_old}计算对应的宽度W_left, W_right, 左边x区间的最小值x_leftmin, 右边x区间的最大最x_rightmax;

3) 判断左右宽度的比值 $\scriptstyle\omega $ 确定是否需要修正, 如果 $\scriptstyle\omega >\omega' $ , 则认为需要对右边进行修正, 如果 $\scriptstyle\omega < -\omega' $ , 则认为需要对左边进行修正, 利用式计算x, 之后使用人行道的拟合函数f(x)作为修正的y.

其中, $\omega' $ 的值为两边人行道宽度的敏感程度, 越小越敏感, 越大越不敏感, 这里取 $\omega' $ 为0.3, $ \omega$ 的确定公式如式(12)所示, x的计算公式如式(13)所示:

$\omega = \frac{{\left( {{W_{\rm{left}}} - {W_{\rm{right}}}} \right)}}{{max\left( {{W_{\rm{left}}},{W_{\rm{right}}}} \right)}}$

(12)

$x = \left\{ \begin{array}{ll} {{x_{\rm{leftmin}}} + {W_{\rm{right}}}},&{w < - \omega } \\ {{x_{\rm{rightmax}}} - {W_{\rm{left}}}},&{w > \omega } \\ \end{array} \right.$

(13)

经过区间约束, 和路面拟合, 就可以利用上一断面的拟合函数和区间约束来提取当前路面.

取某一断面数据为例, 说明四次多项式可以很好的描述路面横断面的曲线. 做出如表1所示的以确定系数R-square作为对比标准的拟合效果对比表: .

通过表可以看出, 无论是路面还是左右两侧的人行道, 当拟合次数达到四次时, 都可以有很好的拟合效果, 而当拟合次数达到五次时, 虽然略有提高, 但是变化不大, 在满足拟合精确的前提下, 为保证算法运算速度, 防止出现过拟合现象, 确定使用四次拟合.

使用四次多项式拟合得到的拟合曲线, 为使图像明显, 对点云数据的横坐标比例进行了缩放, 如图2(a)、图2(b)所示, 可以看到曲线在点云区间拟合效果很好.

2.3 改进区域生长算法对非路面数据分割

通过上述动态数据拟合的算法, 可以将路面数据提取出来, 达到将路面数据与非路面数据的分割, 对路面与非路面数据分割后, 再使用基于改进的区域生长三维点云分割算法^[24]将非路面的各部分分别提取. 其过程为先计算各个点的平均曲率Kn, 计算公式如式(14)所示:

$K({p_i}) = \frac{1}{{4{A_{\min }}}} \times \sum\limits_{j \in N(i)} {\left( {\cot {\alpha _{ij}} + \cot {\beta _{ij}}} \right)\left( {{p_i} - {p_j}} \right) \times n} $

(14)

式中, n为法向量, A_min为p周围一个无限小的区域, α_ij、β_ij分别为连接p_i和p_j边的对角. 选取点云数据中曲率最小的点设置为种子点.

因为路面上的行人, 车辆以及资产之间距离较大, 所以仅用距离作为阈值就可以将各部分分割, 选用式(11)的d作为阈值, 判断邻近点与种子面之间的生长半径是否小于阈值d, 将小于阈值的点添加到当前区域, 当再无点可以添加时, 则开始分割下一区域.

3 实验结果分析

对西安市神州六路部分路段作为无人机数据获取的实验区, 利用无人直升机Scout B1-100搭载RIEGL VUX-1HA激光扫描仪系统对实验区进行扫描, 获取道路点云数据, 激光扫描仪的基本参数如表2所示.

将本算法与区域生长提取路面算法对比, 验证算法的效果. 该路段效果如图3所示, 红色部分为机动车道, 绿色部分为人行道, 蓝色部分为路上资产和路面坑包, 黑色部分为路边环境.

从图3中可以看出, 两个算法很好的将道路上的路灯, 树坑, 墙, 蓝色部分两边为遮挡墙, 蓝色竖线为路灯, 绿色人行道上的几个蓝色小块为树坑. 通过图3(a)路面上有些蓝色点, 人行道上有一片土包可以看出, 本文算法可以将路面上的一些细节提取出来, 而且对于曲率变化较小的土包也可以提取出来, 而图3(b)中的区域生长算法, 对路面上的细节却不敏感, 而且对于曲率变化较小的土包也无法提取. 用识别路面上的树坑, 土包, 凹槽数量作为对比说明两种算法的效果, 列出表3.

为验证该算法的适用性, 对西安市神州六路多处路段进行算法验证, 结果如图4所示, 同样, 红色部分为机动车道, 绿色部分为人行道, 蓝色部分为路上资产和路面坑包, 黑色部分为路边环境, 对于多处路段, 都可以很好的将路面进行提取, 并且受周围环境影响较小. 当路面上有车辆时, 如图4(b)、图4(d)无论车辆是在提取路面的初始断面, 还是在后续断面, 都不会影响到路面的提取.

综上所述, 本文算法可以很好的将路面数据进行识别, 将路面上的资产进行提取, 相比于传统的使用区域生长算法识别路面, 本文算法对路面细节更加敏感. 所以在路面养护和路上资产管理方面有很好的效果.

4 结语

本文在处理道路信息时, 比较依赖于扫描的起始位置, 虽然对遮挡物的抗干扰能力强, 但是, 对于遮挡物太多, 道路信息过少的情况下, 并不能很好的将道路数据的中心通过聚类算法聚类出来; 另外, 在对数据进行分割的时候, 对于杂糅在一起的数据也无法进行分割, 后续可以利用聚类算法, 或者神经网络对数据进行更好的分割; 最后对于和路面的高程数据差过小的数据也无法提取出来, 这种淹没在噪声中的有用数据提取难度很大, 需要后续仔细研究.