基于马尔科夫随机场的三维激光雷达路面实时分割

朱 株，刘济林

（浙江大学 信息与电子工程学系，浙江 杭州310027）

精确感知三维环境是地面自主智能车辆的一项重要工作.为保证安全行驶，自主车辆通常装有激光雷达和相机组成的传感器系统，用于探测障碍区域和地面区域.

早期自主机器人平台一般采用单线激光雷达作为距离传感器，具有数据量小及响应速度快的优点，但是存在探测范围有限的缺点.随着VELODYNE 64线激光雷达的问世，这种三维激光雷达开始替代单线激光雷达，广泛用于自主车平台［1-3］.

现有的三维激光雷达地面分割算法主要有如下2类.

1）基于栅格地图的地面分割算法.这类算法将数据栅格化，投影到一个栅格地图中，分析每个栅格中的三维点云特性来确定该栅格的障碍物属性.Kammel等［4］提出了一种基于栅格中点云最大高程差的地面检测方法，如果栅格最大高程差小于某个阈值，栅格将被标记为地面区域.Himmelsbach等［5］将三维数据以极坐标的方式栅格化，并且用非参数的方法来拟合地平面，分离地面和障碍物.Guo等［6］将栅格图以无向图的方式处理，用马尔科夫随机场将栅格分为4类，从中选取地面.虽然基于栅格的检测方法较为稳定，但是该方法的检测精度取决于栅格大小（往往是20 cm×20 cm），相比于原始三维点云数据精度（0.2 cm），栅格算法精度较小；受限于计算速度，栅格地图范围不大，典型的250×150栅格图相当于50 m×30 m的检测范围，相比于原始数据200 m×200 m的检测范围，大量数据丢失；由于三维数据分布的非均匀性，大量栅格内并没有数据，造成了存储和处理的浪费.

2）基于图的地面检测.三维激光雷达数据结构的特点：任意三维点与其在同一扫描线上相邻角度的点以及前后扫描线相同角度的点组成4邻域系统，可以将三维点云以无向图的方式组织起来.Montemerlo等［7］使用三维数据环之间的距离来判断单个三维点是否属于地面区域（三维激光雷达地面反射点在二维平面上呈现出环的样式，障碍物点的环间距要远小于地的面点间距）.Moosmann等［8］使用局部凹凸性分割不同物体.Douillard等［9］采用Mesh方式构造图，并且通过计算梯度来确定地面点.

相比于栅格类的分割算法，基于图的算法精度较高，并且能处理全部雷达数据，但是现有算法均使用单个三维点作为图的节点，特征单一，容易受到噪声干扰，在起伏道路环境中鲁棒性较差.

本文提出一种新的基于图的算法.将扫描线χ-y平面上的投影曲线进行分段，构建以线段为节点的马尔科夫随机场，通过对线段特征分析建立能量函数，并采用图分割算法求全局最优解.算法不仅能在100 ms内完成，而且在不同场景下分割效果都优于现有基于图的分割算法.

1 算法概述

观察雷达数据在χ-y平面上的投影可以发现，地面数据和障碍物数据在线型上的明显区别.地面区域如图1（a）环状图案所示，单条扫描线障碍区域放大后如图1（b）、（c）所示：在同一扫描线上呈断开状或是折角状.据此，本文提出的地面分割算法采用线段特征.

图1 扫描线在x-y平面上投影的样式Fig.1 Projection of scan lines onχ-y plane

图2 地面分割算法流程图Fig.2 Flow chart of ground segmentation algorithm

图2为地面分割算法的流程框图，该流程分为3个步骤：扫描线投影分段；线段端点精确定位；建立无向图马尔科夫随机场并且求解.

首先输入激光雷达原始数据，并且进行预处理，包括异常距离三维点滤除以及采用GPS进行点云位姿修正.

考虑到在χ-y平面的投影中每条扫描线处于同一直线上的相邻点属性相同，算法首先对每条扫描线进行快速分割.由于存在噪声，分割后相邻线段端点（角点）可能出现定位错误，通过角点检测，准确选取端点，增加算法的精确度.采用线段作为图节点，不仅大大减少了节点数量（从每帧大约130 000个点减少到20 000个点左右）加快图分割的处理速度，而且能获取多个基于线段的特征，增加了算法的鲁棒性.最后通过图分割算法，获取全局优化后的分割结果，提升算法的性能.

2 扫描线投影分割

由于地面不平整以及激光雷达本身具有的检测噪声，扫描线在χ-y平面上投影含有噪声.采用基于最大模糊线段的方法对扫描线投影进行快速分割.

模糊线段（blurred segment，BS）定义为介于直线l1：aχ+by=μ以及直线l2：aχ+by=μ+ω之间的有限离散二维点集，如图3所示.ν为BS的宽度：

图3 模糊线段Fig.3 Blurred segment

一个二维点集模糊线段的参数可以通过计算该点集凸包络的支持线来获取［10］，计算复杂度与点的数目成线性关系.分割算法流程如下：

1）将单条扫描线的映射点顺序输入，如果当前点和前一个点的距离大于阈值T p-d，则将前一个点作为当前线段的终止端点，当前点作为新一条线段的起始端点.相邻2个点的属性都标记为分离点.

2）若当前点和前一个点的距离小于阈值Tp-d，则进行模糊线段参数计算.如果模糊线段的宽度大于阈值Tv，则将前一个点作为当前线段的终止端点，当前点作为新一条线段的起始端点.相邻2个点的属性都标记为连接点.

三维激光雷达在不同距离下有不同的分辨率，Tp-d和Tv均为距离相关的自适应参数：

式中：D为当前点在χoy平面的投影到原点的距离，μ1、μ2为比例因子.由于 VELODYNE HDL64E雷达在不同距离，扫描线在地面上的点间隔不同，对此阈值使用比例因子与距离的乘积.考虑到在10 Hz运行下，激光雷达角分辨率最大为0.18°，2个连续扫描点间距理论上为D×0.18π／180.作为阈值将该值适当放大，取3倍值，计算中μ1=0.009 42.式（3）用于确定分割线段的最大宽度.该阈值不仅受到距离的影响，也受到不同地形环境的影响.本文使用经验值，平坦地形下μ2=0.003，在颠簸环境下μ2=0.005.

分割结果见图4，相邻的线段用2种不同的灰度表示.

3 角点检测

图4 模糊线段分割结果Fig.4 Segmentation by blurred segment

模糊线段有一定的宽度，当连接点很多时不能准确定位在局部曲率最大的角点上.可以通过角点检测算法重新定位连接点.

实空间R上的算子A被定义为以下形式：

二维曲线Γ定义为参数化方程a（s），核函数φ定义为高斯函数，其中u为均值，t为方差：

将a（u）在u=0附近进行泰勒展开：

可得：

积分算子在χ处的值在曲线法方向nχ上正比于曲线曲率κχ，选择局部曲率最大的点作为角点.角点定位前后分割结果如图5（a）、（b）所示.

图5 角点定位前后不同的分割结果Fig.5 Segment results before and after corner detection

4 建立无向图马尔科夫随机场

马尔科夫随机场框架建立了全局标记优化模型，常被用于解决场景分类问题［11］.定义G= （V，E）为无向图，节点v i∈V；边 （vi，v j）∈E对应于相邻的节点对.每条边（vi，vj）∈E被赋以非负权重W（v i，v j），用于描述相邻节点vi及v j间的相似度.

标记推测算法能为每个节点找到概率最大的标记，例如图分割［12］.定义L为一个有限的标记集合，f为标记赋值函数，f（v i）为每个节点赋以标记.假设标记在区域内变化缓慢，而在区域边缘变化剧烈.标记赋值的准确度由能量方程来衡量：

式（8）中等号右侧第一项是数据项，是对节点本身的约束；第二项是平滑项，是对能量方程解的平滑性约束［13］.

本文算法利用激光雷达原始数据结构，构建无向图.把每条线段当作一个节点，所有包含与当前节点线段相邻三维点的线段，与该节点线段组成一个相邻集合.定义与节点相邻的同一扫描线线段为侧邻线段（SN），定义与节点相邻的前一扫描线线段为上邻线段（UN），定义与节点相邻的后一扫描线线段为下邻线段（DN）.根据定义，每条线段只有左右2条侧邻线段，可以有多条上下邻线段.

根据线段内点数依据阈值Tp-n将线段分为2类：长线段（Sl）和短线段（Ss）.由于地面较为平坦，长线段大多属于高概率地面区域（Agh），只有在以下情况长线段被定义为高概率的障碍区域（Aoh）：

1）长线段端点标记为分离点，且该点距雷达原点的距离与侧邻线段分离点距原点距离的差大于阈值Td.

2）长线段端点标记为连接点，且它与侧邻线段夹角与直角的差小于阈值Ta.

3）长线段与上邻线段的梯度大于阈值Tg.梯度定义为线段包含三维点平均高度差与线段在χ-y平面上的平均距离的比值.

短线段大多属于高概率障碍区域，只有在以下情况长线段被定义为高概率的地面区域：

1）短线段端点标记为连接点，与侧邻定义为高概率地面区域的长线段夹角与直角的差值大于阈值Ta.

2）短线段与上邻定义为高概率地面区域的长线段的梯度小于阈值Tg.

本文采用线段点数阈值的经验值Tp-n=6.Td、Tg均为障碍物阈值，障碍物会造成同一扫描线上线段的分离以及相邻扫描线上线段间产生梯度.这2个阈值的取值与地形或者障碍物区分敏感度有关.实验中，城市道路的Td=40，Tg=tan 12°（0.212）；乡村颠簸道路的Td=60，Tg=tan 25°（0.466）.Ta主要用于区分城市道路中的马路牙子和地面，实验中Ta=30.

根据以上准则，可得：

η是符合判断准则下，标记为地面的先验概率值，实验中η=0.8.

为保证局部区域标记的一致性，定义平滑项：

φ（Δh）是根据相邻线段在三维空间中的平均高度差Δh来度量相似性的函数，k为经验系数（实验中k=0.01）：

最后使用图分割的方法（Graph-Cut）优化能量方程，得出标记结果.

5 实验过程及其结果

实验平台如图6所示.三维激光雷达安装在车顶部，除少部分盲区外，能较为全面地获取车体周围的三维环境信息.

图6 三维激光雷达与自主车Fig.6 3D LIDAR and autonomous landing vehicle

分别在城市平坦路面场景和乡村起伏道路场景采集数据，将本文算法结果与局部凹凸性算法［8］、Mesh梯度算法［9］以及人工标记的结果进行对比.图7中第1行是城市平坦路面场景；第2行是乡村起伏道路场景；图7（a）、（f）是对2种场景的相机图片；图7（b）、（g）是三维激光雷达数据人工标记的结果，在本文中该图为灰度图，实际中该图的其中地面区域标记为绿色，障碍区域标记为红色，黑色区域为未知区域；图7（c）、（h）为局部凹凸性算法的结果；图7（d）、（i）为 Mesh梯度算法结果；图7（e）、（j）为本文算法的结果.

在城市平坦路面场景下，Mesh梯度算法（图7（d））结果和本文算法（图7（e））的结果都较好，接近人工标记的结果（图7（a））.局部凹凸性算法（图7（c））结果较差，以局部平面的法向量方向作为判断准则，用单个三维点的邻接4点拟合平面.因此，在构建图时，删除了长度大于某个阈值的边，以防止远距离点用于平面拟合，造成障碍物边缘判断错误.由于在远距离，激光雷达相邻点本身相距很远，导致删除长边后的图出现大量未知区域，影响了检测效果.

在乡村起伏道路场景下，地面起伏变化大，使得基于单个三维点梯度特征的Mesh梯度算法（图7（i））检测效果下降，出现了大量地面区域漏检点；路边大量低矮植物造成了地面区域误检点.局部凹凸性算法（图7（h））不仅在远距离存在未知区域；在近距离，由于地面不平整，使得局部平面的法向量变化很大，导致地面区域漏检，路边大量低矮植物又造成局部平面的法向量变化缓慢，产生地面区域误检点.本文算法（图7（j））效果优于Mesh梯度算法和局部凹凸性算法：由于本算法使用线段的梯度特征，减少噪声对单点梯度的影响，线段夹角特征和分离点距离特征的判断增加了地面区域判断的可靠性，最后使用马尔科夫随机场进行全局优化使得结果具有更好的鲁棒性.

图7 2种场景的分割结果Fig.7 Segmentation results of two kinds of scenes

式中：NTP为地面点被正确标记的数目，NFP为障碍点被标记为地面点的数目，NFN为地面点被标记为障碍点的数目，NTN为障碍点正确标记的数目.

城市场景下的算法准确率如图8所示，横坐标为连续帧的编号，纵坐标为准确率.Mesh梯度算法和本文算法都有较高的准确率.

乡村起伏道路场景下的算法准确率如图9所示，本文算法相比于Mesh梯度算法和局部凹凸性算法，不仅准确率高于这2种算法，而且稳定性也要好于这2种算法.在起伏道路环境中，地面梯度变化剧烈，基于单点局部特征的算法稳定性差不适用于这类场景.本文算法通过多特征全局优化，获得了较高的稳定性，更适合乡村起伏道路场景.

本文算法使用C++实现，在Intel酷睿i5双核2.8 GHz CPU、2 G内存的自主车平台上运行时间

图8 城市场景下算法准确率Fig.8 Accuracy of algorithms under urban scene

图9 乡村起伏道路场景下算法准确率Fig.9 Accuracy of algorithms under rural bumpy scene

在2种场景下运行100帧数据序列，通过正确率（R）对不同算法的结果与人工标记的真值进行量化比较.为每帧61 ms，其中原始激光雷达距离数据转换成三维点云25 ms，扫描线分割15 ms，马尔科夫随机场建立和求解21 ms，符合实时运行的要求.

6 结 语

本文提出了一种适用于多类场景的实时三维激光雷达地面分割算法.算法首先使用最大模糊线段方法对三维扫描先在平面上的投影曲线进行分割，然后采用角点检测算法重新精确定位线段端点，最后建立无向图马尔科夫随机场并求解标记地面区域.在城市平坦路面和乡村起伏路面下的实验结果表明：本文算法相比于其他算法具有较高的鲁棒性，能在起伏路面场景下稳定可靠地为自主车辆提供路面区域信息.

（

）:

［1］GLENNIE C.Calibration and kinematic analysis of the Velodyne HDL-64E S2 Lidar sensor［J］.Photogrammetric Engineering and Remote Sensing，2012，78（4）：339- 347.

［2］HASELICH M，BING R，PAULUS D.Calibration of multiple cameras to a 3D laser range finder［C］∥International Conference on Emerging Signal Processing Applications（ESPA）.Las Vegas：IEEE，2012：25- 28.

［3］杨飞，朱株，龚小谨，等.基于三维激光雷达的动态障碍实时检测与跟踪［J］.浙江大学学报：工学版，2012，46（9）：1565- 1571.YANG Fei，Zhu Zhu，GONG Xiao-jin，et al.Real-time dynamic obstacle detection and tracking using 3D Lidar［J］.Journal of Zhejiang University：Engineering Science，2012，46（9）：1565- 1571.

［4］KAMMEL S，PITZER B.Lidar-based lane marker detection and mapping［C］∥Intelligent Vehicles Symposium.Eindhoven：IEEE，2008：1137- 1142.

［5］HIMMELSBACH M，HUNDELSHAUSEN F，WUENSCHE H.Fast segmentation of 3D point clouds for ground vehicles［C］∥Intelligent Vehicles Symposium（IV）.San Diego：IEEE，2010：560- 565.

［6］GUO C，SATO W，HAN L，et al.Graph-based 2D road representation of 3D point clouds for intelligent vehicles［C］∥Intelligent Vehicles Symposium（IV）.Detroit：IEEE，2011：715- 721.

［7］MONTEMERLO M，BECHER J，BHAT S，et al.Junior：The stanford entry in the urban challenge［J］.Journal of Field Robotics，2008，25（9）：569- 597.

［8］MOOSMANN F，PINK O，STILLER C.Segmentation of 3D lidar data in non-flat urban environments using a local convexity criterion［C］∥Intelligent Vehicles Symposium.Xi’an：IEEE，2009：215- 220.

［9］DOUILLARD B，UNDERWOOD J，KUNTZ N，et al.On the segmentation of 3D LIDAR point clouds［C］∥International Conference on Robotics and Automation（ICRA）.Shanghai：IEEE，2011：2798- 2805.

［10］DEBLED RENNESSON I，FESCHET F，ROUYERDEGLI J.Optimal blurred segments decomposition in linear time［C］∥ Discrete Geometry for Computer Imagery.Berlin Heidelberg：Springer，2005：371- 382.

［11］BESAG J.Spatial interaction and the statistical analysis of lattice systems［J］.Journal of the Royal Statistical Society.Series B（Methodological），1974，36（2）：192- 236.

［12］BOYKOV Y，VEKSLER O，ZABIH R.Fast approximate energy minimization via graph cuts［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，2001，23（11）：1222- 1239.

［13］BOYKOV Y，KOLMOGOROV V.An experimental comparison of min-cut／max-flow algorithms for energy minimization in vision［J］.IEEE Transactions on Pattern Analysis and Machine Intelligenc e，2004，26（9）：1124- 1137.