1. 引言
越来越多的残疾人依靠电轮椅出行。然而电力不足，无法返程会带来很多不便。很多导航系统开始为驾驶电轮椅的残疾人进行路线规划，找出最节能的路线，比如eNav导航系统(eNav navigation system) [1] 。有研究指出路线的耗能与地面的高度呈指数关系 [1] ，一个高精度三维道路地图的输入是路线成功规划的关键。开源地图(Open Street Map, OSM)由于其开源的特性在导航系统中应用广泛，然而道路的高度信息在OSM中非常稀疏，且精度不高(方差高于5米)，有越来越多的研究者关注高精度三维地图的生成。
本文以机载激光雷达点云为基础(均方差仅为20 cm)拓展OSM，生成可供eNav使用的高精度三维道路地图。为了解决机载激光雷达高空作业使得部分地面区域被遮盖而无点云采样点的问题，提出了基于投影的霍夫变换的方法(Projection and Hough Transform, PHT)。
2. 三维道路地图的生成
为整个城市生成三维道路地图的文献较少，大部分论文集中在城市路面环境的分割和分类 [2] [3] 。由于整个城市路面情况复杂，且点云和OSM数据量大，因此，先对局部复杂的环境分割出道路，再将该方法应用于所有区域，是获得整个城市路面高度的便捷方法。但这也对算法鲁棒性提出了更高的要求。
基于激光雷达点云的道路的检测与分割的方法主要可分为以下6类。1) 基于体素的方法；2) 形态学方法；3) 基于道路特征的方法；4) 分割法；5) 基于机器学习和深度学习的方法；6) 拟合平面的方法。
基于体素的方法 [4] [5] 比较简单，往往都假设道路的高度是最低的，这种假设不适用复杂的路面情况，比如斜坡，高架桥等。形态学方法 [6] 一般用来提取道路的边界。它只能粗略估计道路路肩，需要其它后续方法的辅助。道路和其他环境可以由路肩区分开，因此，一些研究 [7] [8] 通过识别路肩点来分割路面。以上方法对于参数的大小比较敏感，假设过于简单，不能应用于整个城市。
基于分割 [9] 的方法需要先验知识的辅助，对于处在边缘的点并不鲁棒。Dzafic [1] 用OPTICS算法分割出道路点集合；然而OPTICS算法对于点云中的缺失和噪声问题不够鲁棒。近几年，基于深度学习的方法在很多公开数据集竞赛上取得了骄人的成果。Caltagirone Luca [10] 应用全卷积网络检测路面并在KITTI自动驾驶数据集上取得第一名的成绩，但是深度学习方法要求数据集场景相似且训练数据量大，测试数据集小。
常见的拟合平面的方法有随机采样一致性算法(Random Sample Consensus Algorithm, RANSAC) [11] [12] [13] ，霍夫变换(Hough Transform, HT) [14] [15] 。RANSAC只能拟合出一个模型，多条道路存在会干扰结果。而HT可以一次拟合多个模型，更适合复杂的城市环境。Xiangyun Hu [15] 应用了WHT提取路面中心线的弧度单元，结果表明对每个点分配不同的权重来投票可以提高路面分割的结果，但是准确计算权重也增加了算法复杂度和对点云预处理的要求。
拟合平面的方法是基于局部路面可近似于平面的假设，不依赖于路肩，护栏等特征。它可以方便的拓展到城市多变的场景中，同时可以根据拟合出的平面直接计算出被遮挡路面的高度，PHT就是基于拟合平面的方法。
3. 算法描述
PHT是逐个处理OSM中的道路点，提取其邻域点云。这样运行内存比载入全部点云更加友好。该方法是基于局部路面可用平面拟合的假设。算法流程如图1所示。可分为3个步骤：1) 基于正交投影将三维路面投影成二维直线，然后利用HT拟合该直线，找出道路候选点集合。相比拟合平面，用HT拟合直线提高了算法的效率；投影也自动过滤了干扰面；HT能一次拟合多条直线，可以很好的适应多条道路存在的情况。2) 利用RANSAC和LSM根据得出的候选点集合拟合出道路平面。3) 为了解决室内通道无点云覆盖，提出了利用相关联的室外通道的高度值以及投影距离加权的方法估计室内通道的高度值。4) 针对坡度异常的道路，利用带权重的霍夫变换(Weighted Hough Transform, WHT)进行再计算。
针对每一个OSM中的道路点，根据XY平面内的距离提取它在激光雷达点云中的邻域点。如果邻域太小，由于OSM本身位置的偏差以及遮挡等问题，目标路面可能不存在于邻域点云中。另一方面，较大的邻域半径会带来更多噪声和干扰面，计算效率也会下降。因此，根据提取出的平面的梯度来实时调整半径。具体为梯度超过25˚ (当地政府提供的先验知识：路面坡度不会超过25˚)，邻域半径将会增加。
获取到的点云不可避免会存在干扰点和异常值，采用基于平均距离的统计分布方法滤除离群点；由于局部路面的法向量变化平缓，将邻域法向量变化程度大的作为干扰点去除。
预处理完点云后，PHT主要由室外通道处理，室内通道处理和坡度异常再计算三个阶段组成。
3.1. 室外通道的处理技术
3.1.1. 正交变换
点云预处理后，下一步是找到合适的投影面使得三维路面可以被投影成直线。如图2所示，首先投影向量应该尽可能垂直于目标路面的法向量，其次要寻找的投影面必须垂直于XY平面。
Figure 2. 3D road surface projected into a straight line
图2. 三维路面投影成直线示意图
利用PCA提取点的法向量时，由于协方差距离矩阵是对称的，因此PCA得到的三个特征向量是相互正交且独立的。最小的特征值对应的特征向量就是点的法向量。将离OSM中目标点p在XY平面上最近的k个点的法向量的平均值近似为p所在路面上的法向量。所以可以利用其余两个特征向量作为一组正交基计算出投影向量n，这样可以保证投影向量尽可能垂直于平面法向量。
(1)
为上文提到的k个点分别利用PCA分解得到的k个特征向量的和，分别对应特征值最大和特征值第二大。a，b为系数。
另一方面，现实中，因为路面是平坦的，因此PHT要寻找的投影面必须垂直于XY平面，即投影向量是平行于XY平面的。
(2)
为向量
的z维度分量。
因为投影向量n的长度不影响投影，为了简化计算，假设a为1，则计算结果为
(3)
图3(a)是OSM中某点的谷歌地图全景。按照计算出的投影向量(图3(b)中绿色向量所示)投影后结果如图3(c)所示。可以发现许多干扰面比如屋顶立面经过投影后是无序的，并不能形成直线。而两个道路平面成功投影成两条直线(图3(c)蓝色直线所示)。图3(d)显示的是k个点的法向量，某些法向量存在轻微偏移目标路面法向量的情况，但是向量的平均操作，提高了算法的鲁棒性。
3.1.2. 利用霍夫变换获取道路候选点
HT的思想是让每个点为预定义的它可能属于的结构投票。投票在所有点上累积，投票最多的模型即为最终拟合出的模型。
通过梯度的约束 [16] 来提高HT的效率，每个点只对梯度附近的参数投票同时忽略梯度异常的投票值。图4是某实际道路点的霍夫空间表示。结果表明，滤除干扰点并利用梯度投票的霍夫空间投票更集中，为后续寻找投票极大值奠定基础。
(a) OSM中某点的谷歌地图全景，纬度为50.7985768，经度为6.0618652；(b) 激光雷达点云，绿色向量为计算出的投影向量，蓝色向量为k个点法向量；(c) 点云投影结果，蓝色直线为两个道路平面投影后的点拟合出的直线；(d) k个点法向量细节图。
Figure 3. Dimensionality reduction of point cloud
图3.点云降维的结果
(a) 没有滤除干扰点的霍夫空间；(b) 滤除干扰点霍夫空间；(c) 滤除干扰点并利用梯度投票的霍夫空间。
Figure 4. Hough Space
图4. 点(纬度为50.7363403，经度为6.1011415)霍夫空间表示
使用高斯滤波器平滑霍夫空间可以去除一些干扰的极大值点，同时也使目标直线的参数的投票值更突出。图5显示了高斯滤波器对于HT的影响，图5(a)和图5(b)的对比说明了滤波有效去除了霍夫空间内有干扰的极大值。
(a) 高斯滤波后的霍夫空间投票数，蓝色点表示极值点；(b) 高斯滤波前的霍夫空间投票数；(c) 高斯滤波后投影，蓝色直线表示极值点对应的直线；(d) 高斯滤波前投影。
Figure 5. The effect of the Gaussian filter
图5. 高斯滤波器的作用
采用文献 [17] 提出的扫描平面的方法去寻找局部极大值。其中邻域太小会使算法更易检测到多余的直线，邻域太大，可能会遗漏目标直线的检测。PHT倾向于选择小的邻域，因为多余的直线并不会影响最终直线的选择，如图5中的(c)和(d)，它们最终拟合的目标直线都是绿色的直线。
3.1.3. RANSAC改善最小二乘法拟合道路平面
在霍夫空间寻找到极大值后，可根据极大值追寻投票给该参数的点的集合。再使用RANSAC将这些点分为内点和外点。内点是满足估计参数的点，外点是一些干扰点。最后，在内点集合中，使用最小二乘法拟合出道路平面。当霍夫空间存在多个极大值，即存在多个平面时，选择XY平面内与离OSM中目标点p点距离最近的极值点。
图6(b)中，紫色的点是道路点集合。利用该先验知识，计算HT拟合的直线和LSM结合RANSAC拟合的平面的投影的误差，分别为44厘米和14厘米。由此可见，RANSAC可被用于去除外点的影响，提高拟合平面的准确性。
拟合出平面后，即可根据p点的x和y值求出对应高度z值。
(a) 道路点(纬度为50.774883，经度为6.0822265)激光雷达点云投影，红色点为全部邻域点；(b) LSM拟合平面投影，绿色直线为LSM拟合的平面的投影，蓝色线为HT拟合的直线。
Figure 6. RANSAC and LSM improve the accuracy
图6. RANSAC与LSM提高了拟合的精度
3.2. 室内通道的处理技术
室内通道也是OSM道路的重要组成部分，然而室内区域并没有点云覆盖，因此这部分道路的高度需要依据相连的室外通道来计算。
Figure 7. The method to calculate the height of indoor ways
图7. 室内通道高度计算示意图
搜索OSM中每一条室外通道，如果只有一个点被标识是室内的，则说明这条室外通道与一条室内通道相连。每一条室内通道都有起始点。起点和终点分别与两条室外通道相连。如图7所示，将室内道路点投影到起始点的连线上，按投影比例和起始点高度差hse即可算出室内每一点的高度。
3.3. 基于权重霍夫变换的异常坡度的处理技术
在处理完OSM中所有道路点后，每条道路相邻的两点都可以根据自身的高度和距离算出坡度。当坡度大于25˚时，该道路一定有某些点高度估计不准确。坡度异常的原因主要有以下3点：1) 道路被遮挡，由于OSM中的目标道路点在XY平面一定落入遮挡物范围内，致使高度估计值实际是遮挡物的高度，因此会造成坡度异常。2) OSM经纬度存在误差。比如图9(a)，图9(b)在OSM都是桥上的两点，结果一点的经纬度在谷歌地图上显示在桥下。这样这两点就存在错误的高度差。3) 阶梯一般也会造成大的坡度，符合预期；PHT使用阶梯旁边的斜坡来代替阶梯的高度。WHT主要针对前两点改善算法。图8是WHT的流程图。首先生成所有异常边集合。当一条道路两条异常边有交点时，将交点加入异常点集合。
对于隧道，立交桥等OSM中有标识的路面，因为HT已经拟合出所有平面，则试着将其它平面的高度值赋给该异常点直到坡度小于25˚。如图9(b)，图9(c)所示，被桥遮挡的路面高度值从桥的高度变成自身高度。
但是，对于没有标识的情况，比如图9(e)所示，由于茂密的树林的遮挡，点云集合中只有极少部分点来自于目标路面，HT依然会为干扰点拟合出一个模型，所以图9(e)中会出现目标点的高度处在树冠中间。但同时异常的坡度也预示着当前提取的模型的不合理，因此，将这些支持当前模型的点的权重减少，再次拟合模型，重复以上步骤直到坡度小于25˚。当有权重小于0时但坡度依然异常，此时增加邻域的半径重新计算。图9(e)是由于树冠的遮挡，高度值估计错误的情况，经过半径的自动增加如图9(g)所示，目标道路点的增多以及图9(f)非道路点权重的降低，使得最终得到了合理的道路估计值，如图9(g)中蓝色点所示。
4. 实验结果与分析
本文与德国亚琛工业大学计算机科学与技术嵌入式实验室合作，利用德国科隆市政府提供的亚琛市的机载激光雷达点云，为亚琛市生成城市级的三维道路地图。
由于获取路面的真实高度需要昂贵的高精度GPS设备，并且由于实验覆盖了整个城市，因此手动测量路面高度进行算法的评估是不现实的。基于Dzafic利用OPTICS生成亚琛市道路高度地图的方法 [1] ，对比了PHT与OPTICS算法的结果，结果表明，与OPTICS算法相比，8.7%的节点两个算法评估的高度有差异且差异大于1米。图10是这些高度差的直方图分布，高度差主要集中在1米到3米之间，占据了86.2%。随机采样500个高度差在1米到3米之间的点，100个高度差在3米到7米的点，以及所有大于7米的点，总计是62个。这样分类的依据是根据现实生活中物体的高度，一般小于3米是一些灌木，集
(a) OSM桥面上一点的谷歌地图显示(纬度为50.7984121，经度为6.1456789)；(b) OSM桥面上一点在谷歌地图上显示在桥下(纬度为50.7981679，经度为6.1460355)；(c) OSM中被桥遮挡的路面点q (纬度为50.7910241，经度为6.1339308)；(d) q的邻域点云，绿色平面为最终拟合的路面；(e) 由于邻域半径过小，点云中目标路面的采样点过少，被树冠遮挡(纬度为50.787027，经度为6.0819526)；(f) (e) 的投影，其中红色点为权重减少的点，紫色点为权重不变的点；(g) (e) 中的点，邻域自动增加，绿色平面是WHT拟合出的路面。
Figure 9. The abnormal incline and the effect of WHT
图9. 坡度异常以及WHT的作用
装箱，公交车等，而3米到7米是一些大树，大于7米一般是建筑物等，这样可以方便分析原因。
随机采样点后，按照场景的类型进行分类，然后根据实地勘察或者谷歌地图和其三维点云人工判断哪个算法更加合理。
4.1. 高度差小于3米的结果分析
在随机采样的500点里，87.3%的点本论文提出的算法评估的更准确。这些点里，OPTICS算法更倾向于选择路面上的一些干扰物作为路面的高度，比如灌木丛或者小汽车的高度，如图11(c)所示。这些情况占据了70%，还有15%是因为存在斜坡的情况也会导致OPTICS估计错误。
PHT估计错误的比例是12.7%，主要是因为这些点在OSM中并不在路面经纬度范围内，如图11(e)所示；而公共汽车的高度不够高，使得相应的路面坡度为23.6˚；因此，WHT不能被调用。
Figure 10. Histogram of height difference which is greater than 1 meter
图10. 高度差大于1米的直方图
(a) 正确率对比；(b) OPTICS预估错误的原因；(c) OPTICS受干扰点影响，紫色点为OPTICS估计值，蓝色点PHT的估计值；(d) 图(c)中点谷歌地图显示(纬度为50.7667512，经度为6.1071091)；(e) 目标点落在了公共汽车顶面上，导致该算法估计错误，紫色点为OPTICS估计值，蓝色点为PHT的估计值；(f) 图(e)中点的谷歌地图显示(纬度为50.7994947，经度为6.0729785)。
Figure 11. Analysis of height difference which is smaller than 3 meters
图11. 高度差小于3米结果分析
4.2. 高度差大于3米小于7米的结果分析
高度差在3米和7米之间的点有486个点，随机采样了100个点。预测成功了85%的点，这些点中，
(a) 正确率对比；(b) OPTICS预估错误的原因；(c) 图(d)中点谷歌地图显示(纬度为50.7905282，经度为6.0721709)；(d) OPTICS受干扰点影响，紫色点为OPTICS估计值，蓝色点为PHT的估计值；(e) 图(f)中点的谷歌地图显示(纬度为50.776214，经度为6.0439491)；(f) 目标点落在了地面上而非目标路面上，导致PHT估计错误，紫色点为OPTICS估计值，蓝色点为PHT的估计值。
Figure 12. Analysis of height difference which is between 3 meters and 7 meters
图12. 高度差大于3米小于7米结果分析
(a) 正确率对比；(b) OPTICS预估错误的原因；(c) 图(d)中点谷歌地图显示(纬度为50.7872251，经度为6.0816284)；(d) OPTICS受干扰点影响，紫色点为OPTICS估计值，蓝色点为PHT的估计值；(e) PHT可以正确估计被建筑物遮挡的路面高度，紫色点为OPTICS估计值，蓝色点为论文提出的算法的估计值；(f) 图(e)中点的谷歌地图显示(纬度为50.7676544，经度为6.0920514)。
Figure 13. Analysis of height difference which is greater than 7 meters
图13. 高度差大于7米结果分析
路面被遮挡的情况约占90%，大部分是被树和建筑物遮挡。即使如图12(d)所示，路面采样完整，并没有被树冠影响，OPTICS依然会错误估计，被干扰点的高度影响。PHT对遮挡情况和干扰点更加鲁棒。
在其余15%的情况下，由于OSM的经纬度偏差，使得目标点在XY平面上与其它路面更加相近。比如如图12(f)所示，目标点正好落在上下平面的交界处，与下面的地面更加近，造成了PHT估计值偏低。
4.3. 高度差大于7米的结果分析
高度差大于7米的点一共有62个，大约87%的情况下，PHT预测的更准确。而OPTICS表现不佳的原因主要是这些区域大部分是在植被茂密的区域或者室内区域。比如如图13(d)所示，OPTICS估计的高度值是树的高度。相反，如图13(d)、图13(e)所示，PHT可以准确估计被树冠或者建筑物如火车站顶棚遮挡的路面的高度。
有8个点PHT估计不准确，主要原因是一条路上只存在一条边坡度异常，而到底是边的起点还是终点估计错误，该算法不能处理这种情况。不过，这样的情况仅仅只有0.3%的点。
5. 结论
本论文提出了生成城市级三维道路地图的通用性算法PHT解决了无点云覆盖问题。该方法的创新型在于：1) 将三维点云投影到二维，使路面近似投影成一条直线。降维提高了霍夫变换的效率，同时也过滤了一些干扰面的影响，便于后续拟合局部道路平面；2) 使用WHT解决了OSM经纬度有误差的问题；3) 自动调节点云邻域大小，适应点云密度的变化。
结果表明，相比OPTICS算法，准确率更高。该方法可操作性强，对于点云被遮盖的情况以及干扰点具有鲁棒性。
致谢
对东南大学生物医学学院和德国亚琛工业大学计算机科学与技术学院嵌入式实验室表示感谢，对给予转载和引用权的资料、图片、文献、研究思想和设想的所有者，表示感谢。