道路是关系国民经济的重要基础设施, 道路信息的及时、准确获取与更新对交通管理、自动车辆导航及应急事务处理均具有重要意义.多回波激光雷达(Light detection and ranging, LiDAR)可以快速获取地表高精度的3D点云数据及回波、强度等辅助信息, 该技术的出现使得道路3D信息的快速提取成为可能.

已有经典的面向LiDAR数据的道路提取方法中主要采用的数据结构有2D格网^[1-10]、不规则三角网^[11-12]和点云^[13-26].其中, 前两者的同一平面坐标只能对应一个高程值, 用其表达3D LiDAR数据存在的信息损失会影响提取结果的完整性且其提取结果为2D形式, 难以直接构建3D道路网、反映道路的空间结构; 点云是LiDAR数据的原始表达, 可完全保留原始的3D信息, 但其并未明晰表达内部各LiDAR点的邻接信息, 由此造成了算法设计的困难.为了克服上述限制, 本文提出了一种基于灰度体元模型(Grayscale voxel model, GVM)的3D道路提取算法.

GVM是将离散、非规则分布的空间点集表示为以体元为基本单元的规则空间数据结构, 并依据体元内LiDAR点的平均强度的离散化程度对体元赋值.其优势在于: 3D数据结构更有利于综合利用多次回波信息、易于移植2D图像处理算法或在此基础上开发新算法; 各体元间隐含有几何拓扑关系, 易于操作和处理; 同时融合了LiDAR数据的高程和强度信息, 可用于高程或强度信息不同的目标分类.本文提出的算法即很好地利用了上述优势, 基于道路"空间结构上具有连通性、反射强度上表现出一致性"的特性、利用3D连通区域构建的理论构建道路提取模型.算法设计简单且道路提取结果—3D道路体可直接用作灰度3D道路模型, 相较于传统的道路提取算法具有明显优势.

1.   GVM构建

包括将场景空间依据体元分辨率划分成3D体元格网、依据体元内LiDAR点的强度值为体元赋值等步骤, 如图 1所示.另外, LiDAR数据中通常包含的异常数据为非真实目标信息、呈现噪声特性, 其存在会影响GVM构建的准确性及效率, 因此, 必须首先予以剔除.

图 1  GVM构建流程图

Fig. 1  Flowchart of GVM construction

下载: 全尺寸图片 幻灯片

1.1   异常数据剔除

异常点通常占少数且位于局部高程最高或最低的位置.因此, 若将场景空间以LiDAR点的平均点间距为格网单元划分为3D格网, 异常点将位于包含LiDAR点的最高或最低的格网中.因此, 可首先定位上述格网并将其作为候选异常数据集.然后, 对各候选格网, 比较其和给定邻域内(邻域必须足够大, 如15 m$\times $15 m)包含LiDAR点的格网的平均高程的高差, 若其绝对值大于给定阈值$T_{ed}$, 则对应的格网内的LiDAR点判作异常点, 予以剔除; 否则予以保留. $T_{ed}$可取各个候选异常格网和其周围格网的高差绝对值的3倍标准差.剔除异常数据集记做$Q$, $Q$ = {$q_{i}(x_{i}$, $y_{i}$, $z_{i})$, $i$ = 1, 2, $\cdots$, $n$}, 其中, $i$是LiDAR点的索引, $n$是LiDAR点数, $q_{i}$代表第$i$个LiDAR点, ($x_{i}$, $y_{i}$, $z_{i})$为第$i$个LiDAR点的坐标.

1.2   场景空间范围确定

由包含剔除异常数据集的轴向平行包围盒(Axis-aligned bounding box, AABB)确定. AABB = {($x$, $y$, $z)$|$ x_{\rm min}$ $\le x$ $\le x_{\rm max}$, $y_{\rm min}$ $\le y$ $\le y_{\rm max}$, $z_{\rm min}$ $\le z$ $\le z_{\rm max}$}, 其中, $x_{\rm max}$ ($y_{\rm max}$, $z_{\rm max})$ = max {$x_{i}$ ($y_{i}$, $z_{i})$, $i$ = 1, 2, $\cdots$, $n$}, $x_{\rm min}$ ($y_{\rm min}$, $z_{\rm min})$ = min {$x_{i}$ ($y_{i}$, $z_{i})$, $i$ = 1, 2, $\cdots$, $n$}分别代表$Q$中$x$、$y$和$z$坐标的最大值和最小值.

1.3   体元分辨率的确定

体元分辨率是GVM构建过程中最关键的参数, 其值过高将导致数据的过度表达及目标的过度分割; 该值过低将导致数据的欠表达及目标的欠分割.为此, 需选择合适的体元分辨率以保证后续道路的提取质量.考虑到LiDAR系统理想采样时各LiDAR点呈间隔均匀的格网分布, 此时$x$、$y$方向的体元分辨率($\Delta x$, $\Delta y)$可依据LiDAR点的平均点间距确定^[27], 即:

其中, $A_{xy}$是LiDAR点水平投影所得2D点集的沿轴向最小外接矩形的面积.同理, $z$方向的体元分辨率$\Delta z$可将上述理论扩展至$XZ$、$YZ$平面, 依式(2)确定,

其中, $A_{xz}$ ($A_{yz})$是LiDAR点的$XZ$ ($YZ$)平面投影所得2D点集的沿轴向最小外接矩形的面积.式(2)中取最小值是因为其代表构建的体元模型和LiDAR数据间存在更少的精度损失.

1.4   场景空间划分为3D体元格网

给定体元分辨率即可将AABB划分为$r $行、$c $列、$l $层的3D体元格网, 可用一个3D阵列表示.设$V$是3D阵列中的体元集合, $V $= {$v_{j }(r_{j}$, $c_{j}$, $l_{j})$, $j$ = 1, 2, $\cdots$, $m$}, 其中, $j$是体元索引, $m$是体元数, $v_{j}$是3D阵列中第$j$个体元, ($r_{j}$, $c_{j}$, $l_{j})$是第$j$个体元的坐标.

1.5   点云映射及体元赋值

依式(3), 将$Q$中的各LiDAR点映射到$V$中, 并对含有LiDAR点的体元赋值LiDAR点强度均值、不含有LiDAR点的体元赋值0.若上述体元值的灰度等级非256, 则离散化到$\{0, {\cdots}, 255\}$得到体元灰度值.其中, 0值体元对应背景, 称作背景体元, 其余体元对应目标, 称作目标体元.构建所得即为GVM, 其本质为3D灰度图像, 后续的道路提取应用将基于该图像进行.

另外, $Q$中各LiDAR点的强度可能存在异常值, 其亦不能反映真实的目标信息, 必须在体元赋值前予以剔除.剔除方案为:统计$Q$中各个LiDAR点的强度值, 并以直方图的形式可视化显示统计结果; 进而目视判定真实目标的最高和最低强度值, 分别记做$T_{h}$和$T_{l}$, 将强度值高于$T_{h}$和低于$T_{l}$的LiDAR点判定为强度异常数据.

2.   基于体元的道路体提取

由于单条道路为连续、光滑的曲面, 道路点云的激光反射强度相近且与周围邻近地物具有明显的区别^[9], 而GVM是对各目标空间位置及其激光反射强度的3D离散化表示, 所以属于单条道路的体元在GVM中不但是3D连通的而且这些体元的灰度值接近(备注:若机动车等目标遮挡导致道路部分区域缺失LiDAR点, 则3D连通区域中间可能出现空洞, 但不影响道路目标的3D连通区域的构建).因此, 本文提出基于3D连通区域标记算法提取道路目标.首先, 依据道路的灰度特性选取种子体元; 进而, 将与种子体元3D连通且灰度差小于设定阈值的体元标记为道路体元; 最后, 利用数学形态学优化道路提取结果获得道路体元集.

2.1   种子体元选取

以顶视图的方式可视化显示GVM, 手工选取道路种子(一条道路选择一个种子即可.但实际应用中, 为了避免错误选择一个小方框内的体元作为种子), 种子体元集记作$V_{s }$ = {$v_{k}(r_{k}$, $c_{k}$, $l_{k})$, $k $ = 1, 2, $\cdots$, $t$}, 其中, $k$是种子的索引, $t$是种子数.

2.2   标记道路种子及其3D连通区域

3D连通区域构建方案为:对任一种子体元$V_{s}$, 依次从$V_{s}$的6邻域(或18、26及其他邻域尺度, 如图 2所示)内体元出发, 遍历所有邻接体元, 直至GVM中和$V_{s}$有路径相通且灰度差小于阈值的体元都被标记.算法详细步骤如下.

图 2  邻域尺度

Fig. 2  Adjacency size

下载: 全尺寸图片 幻灯片

1) 初始化, 设置存储种子体元的初始栈, 并标记种子体元为道路体元;

2) 从初始栈栈顶弹出一个元素, 获取其空间邻域内未标记的目标体元与该元素的灰度差, 若差值绝对值小于给定阈值$T_{i}$, 则标记为道路体元并存入栈中;

3) 如果栈中为空, 终止程序; 否则, 返回2).

在上述标记过程中应用不同的邻域尺度、$T_{i}$会得到不同的道路提取结果.邻域尺度及$T_{i}$过小会导致道路体提取结果不完整; 邻域尺度及$T_{i}$过大则影响效率及准确性.最佳邻域尺度及$T_{i}$将在实验中确定.

算法基于"相邻且灰度值接近的目标体元可视为同类"的思想, 是传统2D连通区域标记算法的3D扩展应用.算法综合利用了高程及灰度信息, 且在利用高程信息时, 将高程接近条件转换为了空间邻域的搜索, 很好地利用了体元结构内部隐含的邻域关系, 设计简单、高效.

2.3   3D道路体结果优化

实际生活中, 道路与停车场等目标邻接且高程及灰度均非常接近, 因而道路提取结果中往往存在一些非道路数据.此处利用数学形态学去除非道路体元.因为道路具有一定的方向性和宽度、呈线状, 因此, 设计线型结构元素检测并提取道路.具体方案为:

1) 将第2.2节所得道路体元投影至$XY$平面, 并将其转换为二值图像$I_{2}$, 其中道路体元对应的像素赋值1, 非道路体元对应的像素赋值0;

2) 采用2 $\times $ 2的方形结构元素对$I_{2}$进行腐蚀得到影像图$I_{3}$, 以断开细小的连接;

3) 针对非道路区域类型, 在$I_{3}$中采用线型结构元素(大小如3或5, 角度如70°)对图像进行开运算, 得到影像图$I_{4}$;

4) 统计$I_{4}$中连通区域的面积, 通过设定面积阈值将连通域小于设定阈值的区域去掉, 得到去除局部非道路区域的影像图$I_{5}$;

5) 将$I_{5}$中的道路点投影到GVM中得到3D道路体.

3.   实验数据

采用国际摄影测量与遥感协会(International Society for Photogrammetry and Remote Sensing, ISPRS)第三工作组(http://www.itc.nl/isprswgIII-3/filtertest/)提供的包含不同路网复杂度的城区LiDAR数据(CSite2和CSite3, 如图 3所示)作为实验数据检验算法的有效性和可行性.

图 3  实验数据

Fig. 3  Experiment data

下载: 全尺寸图片 幻灯片

数据由Optech ALTM系统获取. CSite2和CSite3的大小分别为0.5 km$\times 0.39$ km和0.42 km$\times\, 0.63$ km、平均点密度分布为1.54个/m$^{2 }$和1.59个/m$^{2}$, 点数分别为408 921和320 477, 回波次数为2次. CSite2中道路高程起伏较大, 道路与停车场、水泥路面相连, 道路中分布有立交桥. CSite3为路况较简单的单层路网, 道路两侧主要为建筑物和水泥路面.利用商业软件, 采用手工分类的方式从实验数据中分离出道路点集, 将其作为标准数据定量评价提出算法精度.

4.   结果与讨论

本文提出的算法的提取结果是道路体元, 而标准数据中则是LiDAR点.为了定量评价提出的算法的精度, 实验首先统计提取的道路体元内包含的LiDAR点, 然后与标准数据对比, 采用式(4)中完整率($R_{com})$、准确率($R_{cor})$、总误差($R_{e})$以及总体质量($R_{q})$等指标定量评价提出的算法的精度^[28].

其中, $TP$是正确提取的道路LiDAR点数, $ FP$是误识为道路的LiDAR点数, $FN$是未提取出的道路LiDAR点数.

4.1   参数的敏感性分析

在相同条件下, 不同邻域尺度及$T_{i}$的道路提取算法的提取结果的总误差见表 1.

表 1  各不同邻域尺度及$T_{i}$下的道路提取算法总误差

Table 1  Total errors of the proposed algorithm with different adjacent sizes and $T_{i}$

CSite2/CSite3提取结果总误差(%)
$T_{i}$	6邻域	18邻域	26邻域	56邻域	64邻域
1	46.30/40.64	32.26/31.87	25.89/24.15	20.47/17.94	17.98/12.44
2	39.57/32.32	28.34/25.02	19.07/18.87	12.81/8.29	18.73/15.65
3	34.25/30.06	26.16/19.21	20.83/10.33	15.58/13.84	23.16/20.45
4	37.29/33.12	30.21/20.04	24.56/12.87	26.74/16.51	29.86/25.31

下载: 导出CSV 

| 显示表格

当$T_{i}$不变时, 对比不同邻域尺度下的道路提取结果的总误差可知:邻域尺度的增大并不意味精度的必然提高.以6邻域为例, 种子的信息只能向其6个基本方向扩展, 而18、26等其他邻域则可同时向更多的方向进行扩展, 随着邻域尺度的增加, 一些非道路体元也有可能被错判为道路体元, 从而在提取结果中引入误差, 这可以解释为何部分64邻域对比56邻域精度反而有所降低.

当邻域尺度不变时, 对比不同$T_{i}$下的道路提取结果的总误差可知: $T_{i}$的增大并不意味精度的必然提高. $T_{i}$过大则可能将非道路LiDAR点错判为道路, 从而在提取结果中引入误差.

对比不同邻域尺度、$T_{i}$下的道路提取结果的总误差可知: 56邻域为最优邻域尺度、2为最优灰度差阈值.

4.2   道路提取结果及其在道路网构建中的应用

最优邻域尺度及$T_{i}$下的3D道路提取结果如图 4所示.

图 4  道路提取结果

Fig. 4  Road extraction results

下载: 全尺寸图片 幻灯片

由图 4可知:提出的算法的道路提取结果为3D道路体元, 其是一种新形式的3D道路表示形式—3D道路体; 图 4(a)中3处圆圈圈出的区域中有3座桥梁, 桥梁与道路的灰度值差别较大, 由此导致道路提取结果在该3处发生了断裂; 所有的道路均被成功提取, 但是提取结果中仍然存在一些与道路相连十分紧密的非道路区域未被剔除.

提出的算法所得3D道路体可进一步处理以构建3D道路网.如图 5(a)和图 5(b)所示, 对3D道路体结果进行空洞填补、提取道路骨架线^[29]、去除毛刺短分支、多余通路和环形通路后所得的道路网是一种新形式的基于体元的3D道路网.由图 5(c)和图 5(d)中道路网和GVM叠加顶视图可知, 道路网构建完整.由图 5(a)中圆圈所示, 道路网发生了明显的断裂.断裂的原因在于立交桥的桥下由于被遮挡而出现了LiDAR点的缺失.但是3D道路网的优势在于能明确这里是两条道路, 即立交桥结构.即使存在断裂, 也可以通过某些路径算法预测道路的连通性.倘若是基于2D图像构建道路网则无法判定立交桥结构的存在.这也体现了基于体元的道路提取算法的优势, 其可以以3D的形式直观地呈现3D立体结构.

图 5  道路网

Fig. 5  Road network

下载: 全尺寸图片 幻灯片

4.3   定量精度评价

提出的算法的定量精度评价如表 2所示.

表 2  提出的算法的精度

Table 2  The accuracy of the proposed algorithm

数据	$R_{com}$(%)	$R_{cor}$(%)	$R_{q}$(%)
CSite2	84.83	80.76	70.57
CSite3	88.71	81.50	73.84

下载: 导出CSV 

| 显示表格

表 2表明:算法提取的道路的总体质量可达70%以上, 平均完整率、正确率分别为86.77%、81.13%.影响完整率的关键因素在于邻域扩展过程中部分道路边缘点的灰度值变化较大而未被提取; 影响正确率的关键因素在于仅利用高程及灰度信息无法完全抑制道路向周围高程及灰度信息均接近的非道路目标(如停车场)的扩展.

4.4   算法对比实验

为了体现本文提出的算法的优势, 首先展开其与商用软件TerraSolid的对比实验.其中, 各实验数据的TerraSolid道路提取指标如表 3所示.

表 3  Terrasolid道路提取精度

Table 3  The road extraction accuracy of Terrasolid

数据	$R_{com}$ (%)	$R_{cor}$ (%)	$R_{q}$ (%)
CSite2	88.2	64.1	59.1
CSite3	93.8	56.9	54.8

下载: 导出CSV 

| 显示表格

表 3表明: TerraSolid提取的道路的平均完整率、正确率及总体质量分别为91%、60.5%及57.0%.对比表 2和3可知:对比本文提出的算法的道路提取精度, TerraSolid的完整率虽然具有绝对优势, 但其准确率较低, 由此导致其道路提取总体质量较差.原因在于: TerraSolid的道路提取是首先滤波得到地面点, 然后依据设定的道路的强度范围获得道路点.这种将满足某些强度特性的地面点判作道路点的道路提取方法虽然可以获得较高的完整率, 但不可避免地会将一些非道路点误判作道路点, 从而影响道路提取结果的正确率, 如图 6中错误分类的点所示(以CSite2为例).而本文提出的算法则综合利用了道路的几何和强度特性, 因此其道路提取的正确率及总体质量均较高.

图 6  CSite2的Terrasolid提取结果误差图

Fig. 6  Error distribution of Terrasolid for CSite2

下载: 全尺寸图片 幻灯片

另外, 鉴于传统的道路提取算法多基于2D图像设计(以下简称2D道路提取算法), 特开展其和本文提出的算法的性能对比实验.

实验所用2D道路提取算法的步骤如下: LiDAR点云数据读取; 异常数据剔除; 剔除异常数据集投影至$XY$平面得到2D点集; 求取2D点集的沿轴向最小外接矩形以确定其范围; 根据平均点间距确定格网分辨率; 依据格网分辨率对最小外接矩形进行划分得到2D格网; 将2D点集映射到2D格网, 并将包含LiDAR点的格网赋值LiDAR点的平均强度及平均高程, 否则赋值0.若格网强度值的强度等级非256级, 则将格网值离散化到$\{0, {\cdots}, 255\}$; 可视化显示2D灰度图像, 手工选取道路种子点(和本文提出的算法的道路种子选取相同区域); 标记与种子点2D连通、高程差小于设定阈值且灰度差小于设定阈值的像素作为道路像素, 完成基于2D灰度图像的道路提取.如图 7所示, 以CSite2为例, 在相当条件下(3D灰度图像邻域扩展时采用最优的56邻域、灰度差阈值2; 2D图像邻域扩展时采用20邻域、灰度差阈值为2, 高程差阈值为1米(1米与"3D灰度图像中的56邻域对应的高程差"相当)), 以相同平面位置的种子点为基础, 向周围做邻域扩展所得基于2D和3D图像的道路提取结果顶视图.

图 7  2D和3D道路提取算法结果对比

Fig. 7  Comparison of 2D and 3D road extraction algorithms

下载: 全尺寸图片 幻灯片

对比图 7(a)和7(b)可知: 2D和3D道路提取算法均提取出了所有主干道路, 但2D道路提取结果包含更多的噪声, 如图 7(a)中实线框所示, 而3D道路提取结果则噪声较少, 如图 7(b)中实线框所示; 对比图 7(a)中虚线框区域, 该区域对应为一座桥, 基于2D图像的道路提取结果则误判作道路, 而基于3D图像的道路提取结果则未误判; 对比图 7(b)中虚线框区域, 该区域包含立交桥, 立交桥的3D结构细节可直观呈现, 如图 4(a)所示; 但是立交桥的3D结构信息则无法在基于2D图像的道路提取结果中直观呈现.

5.   结论

本文针对机载LiDAR数据灰度体元模型构建及在此基础上的3D道路提取展开研究, 提出了一种基于体元的道路提取算法.该算法以3D连通区域构建理论为基础, 将LiDAR数据中的目标提取从点云聚类等传统方式转换成基于体元的空间邻域关系的搜索标记方式, 综合利用了GVM内部隐含的几何、拓扑关系及灰度信息, 算法设计简单.定量评价算法精度实验的结果表明:无论是复杂、简单的道路网, 算法提取的道路完整率、正确率均可达80 %以上, 总体质量可达70 %以上.算法的不足之处在于:需要人工交互选择道路种子, 自动化程度有待进一步提高; 算法仅利用了机载LiDAR数据的几何和强度信息, 由此可能导致的停车场等与道路邻接、高程且反射强度均接近的地物的错误提取.后续研究可尝试基于多值体元模型(以多光谱LiDAR数据为数据源, 体元赋值其内的LiDAR点的色彩信息)的3D道路提取研究, 以提高复杂场景下道路目标提取的准确性.