道路是空间信息系统中的重要基础数据元素，利用计算机技术从遥感影像中自动提取道路是一个长期的研究热点^[1]。早期，受限于卫星传感器的分辨率，道路在遥感影像中表现为线状地物。因此，道路提取的研究多集中于线状目标提取^[2]，如基于种子点及基于最小二乘模板匹配的航空影像道路提取算法^[3]等。随着遥感影像的空间分辨率逐步提高，道路提取的方法变得更为复杂，如形态学方法^[4]、图割方法^[5]等。道路的边缘、几何和拓扑关系、尺度特性等也被用于道路识别^{[6, 7, 8, 9]}。机载激光雷达(LiDAR)作为一种新型传感器，不仅能快速、精确获取地面三维坐标信息，同时也能记录激光反射强度^[10]，能从空间坐标和辐射两个方面表现地物的特性，为空间地物的识别与提取提供新的手段^[11]。

利用LiDAR数据进行道路网络提取的研究已经较多。如Alharthy及Clode等通过点云数据滤波移除非地面点，并利用强度信息进行简单的阈值分割得到道路区域^{[12, 13]}。Rieger等利用高分辨率遥感影像来检测平行线，并与LiDAR融合得到道路^[14]。除了道路之外，Sithole等人还利用LiDAR点云提取桥梁信息^[15]。实际上，由于机载LiDAR数据滤波算法的成熟^[10]，通过LiDAR高程信息能够有效地消除建筑物等干扰道路提取的信息。然而LiDAR点云的强度信息由于噪声较大，并没有引起足够的重视。

研究发现，对LiDAR数据进行地形滤波、强度分割等处理即可初步提取道路区域，然而这些区域中包含较多非狭长形、无规律的假道路区域，严重影响了提取效果。为此，本文提出了一种基于双十字剖线及狭长性的道路区域判别方法，通过狭长性计算区域权值，剔除噪声与假道路区域，提取出狭长状、强度较一致的道路区域。通过对道路区域进行细化和平滑，计算道路中心线，形成道路网。最后，通过实验验证了本文所提方法的有效性和可行性。算法流程如图 1所示。

图 1 算法流程图 Fig. 1 Flowchart of the Algorithm

图选项

道路区域的初步提取首先利用成熟的地形滤波算法，依据LiDAR数据的高程信息滤除建筑物、植被、电力线等非地面点，得到地面点集；然后根据选取道路区域的强度值范围，对地面点进行强度阈值分割，分离地面点中道路与非道路点，从而提取出初步的道路区域。具体步骤如下。

1) 地形滤波。采用Axelsson^[16]提出的基于自适应三角网模型对地形曲率进行约束的方法，分离点云中的地面点与非地面点。该方法给定地形种子点，构建稀疏三角网，对初始地形曲率进行近似，通过选取与所在三角形曲率变化最小的点作为候选点，在满足曲率约束的条件下，添加候选点到三角网中，逐步迭代加密，直至无满足曲率约束条件的候选点为止，最终获取连续的地面点集。

2) 强度内插。采用基于三角网的内插方法，将地面点根据位置与强度信息内插成栅格强度图像。其中，内插栅格点的强度根据其落于三角网三角形的顶点强度，结合距离加权而计算得到。

3) 道路区域的初分割。LiDAR数据中道路因其自身特有的材质属性，在强度反射上与其他地物有着显著的不同。通过确定代表道路区域的强度范围，利用阈值分割选取出道路区域。在本文实验中，该阈值通过人工选取。

经过以上步骤，可初步提取出道路区域，如图 2所示。这些区域中包含大量的非道路区域，其来源主要有两个方面：一方面，室外地面停车场和操场等在高程、材质上与道路非常接近，无法通过以上步骤将其排除；另一方面，实际数据中，常会因LiDAR强度信息的光谱分辨率低和噪声大等因素导致有些地物与道路区域的强度混淆，影响了提取的效果。

图 2 道路初步提取结果示意 Fig. 2 Results of Initial Road Extraction

图选项

从图 2可以看出，正确的道路区域在几何形状上表现为有规则的狭长状。干扰区域如停车场、操场等非道路区域则表现为无规则的块状。故而，本文设计一种双十字剖线贯穿道路区域，统计区域中的每个点在4个方向上的规律性及狭长性，并分别通过规律因子和形状因子对其进行描述；再将规律因子和形状因子加权组合为该点的道路指数，最后通过道路指数的高低判断该点是否为道路。

沿某一方向作剖线，该剖线将灰度图像上这一方向的一系列点相连，形成一个点串，如图 3(a)所示；在该点串上，若存在由连续的道路区域点构成的子串，则保留穿过该子串的线段，如图 3(b)所示；若存在由连续的背景区域点构成的子串，则删除穿过该子串的线段，同时，考虑到存在一定的噪声，若存在剖线端点距离非常近的两根剖线，则需要将它们相连，如图 3(c)所示；由此形成了若干根剖线段，如图 3(d)所示。对灰度图像中该方向的所有系列的点做剖线，形成这个方向上的剖线组，同时在该方向的垂直方向上做另一剖线组，形成十字剖线组。同理，变换一个角度(45°)后，再做一个十字剖线组，与之前的十字剖线组构成双十字剖线组。图 4和图 5分别示意了规则、细长的道路区域，以及不规则的非道路区域的双十字剖线及与剖线L垂直相交的剖线组L₀,L₁,L₂,…,L_n－1。

图 3 剖线段的形成过程 Fig. 3 Formation Process of a Cut Profile

图选项

图 4 道路区域的双十字剖分线与剖线组 Fig. 4 Cross Profiles and a Group of Profiles of the Road Region

图选项

图 5 非道路区域的双十字剖分线与剖线组 Fig. 5 Cross Profiles and a Group of Profiles of the Non-road Region

图选项

权值计算包括两个方面，剖线L的权值Q_L的大小体现了该剖线所在区域在剖线方向上的狭长性与规律性；点P的权值Q_P的大小体现了该点位于道路区域中的概率大小。通过计算这两方面的权值，可以全面地判别道路与非道路区域。

双十字剖线组权值的具体计算方法如下。

1) 剖线权值Q_L的计算。

剖线组里的剖线段L，在其垂直方向的剖线组里，存在着若干根与之相交的剖线组S={L₀,L₁,L₂,…,L_n-1}，如图 4所示。每根剖线的权值Q_L=αQ_A^T。其中α为权值系数矩阵，表征规律性与狭长性的比重。Q_A为权值因子矩阵，由权值因子Q_i和Q_j构成。

权值因子Q_i表征剖线组S内剖线长度的变化大小(规律因子)，体现了区域的规律性。其中，d_Li为剖线组S中L_i的长度。因为道路一般宽度不变，所以若剖线L在道路区域中，剖线组S中剖线长度变化较小(图 4)，则Q_i值较大，说明具有较强的规律性。若剖线L在噪声中，则Q_i值较小，说明规律性不明显。Q_i表达为：

权值因子Q_j表征剖线L的长度与剖线组S长度均值的比值(形状因子)，体现了区域的狭长性。通常道路区域表现为狭长状，在垂直相交的两个方向上具有差异很大的剖线长度。若Q_j值较大，表明此区域在此方向上较狭长，符合道路特征。若Q_j值较小，表明此区域非狭长区域，有可能为停车场等非道路区域。Q_j表达如式(4)：

2) 点的权值Q_P的计算。

对于某个区域而言，过其中一点P存在由双十字的4条剖线组成的剖线组M=L₁,…,L_n|(1≤n≤4)。剖线组M中的每条剖线的权值代表了各个方向上该点所在区域的狭长性与规律性。当某条剖线与该道路走向的方向一致时，可认为该剖线具有最强的规律性和狭长性，它的权值为过点P的所有剖线中最大的。该权值真实地体现了该区域与道路区域的相似程度。故点P的权值Q_P为过点P所有剖线里剖线权值最大的那个值。Q_P表达为：

根据十字剖分线对区域形状的统计，区域中的每个点均含有权值。权值的高低表示该点所在区域在4个方向中最接近道路方向上的狭长性和规律性。道路区域的狭长性和规律性决定了道路点的权值较非道路区域高。为此，依据得到的统计权值，进行阈值分割，分离道路与非道路区域。

设道路区域分割阈值为ρ_R，则对于区域A里的每一个点，若点P的权值Q_P大于ρ_R，则为道路点，否则为非道路点。通过比较后，得到道路点点集R=P_i∈A:Q_i>ρ_R.本实验中，道路区域分割阈值ρ_R为人工选取。

道路中心线的提取实现比较简单：首先通过图像的细化，得到中心线条；然后将这些点连起来形成连续的曲线；再进行平滑和简化，得到道路中心线。

本文选取Optech 公司生产的激光扫描仪对尼亚加拉地区的居民区扫描生成的两块数据A、B作为实验数据。这两块数据拥有两次回波，平均点距约1.2 m。从该居民区的影像上看，这两块区域都拥有较多的房屋及植被，都存在宽阔的停车场。房屋、植被及在反射强度上类似于道路的停车场等都是提取道路过程中的干扰因素。

数据A经过滤波后，在经过基于三角网的强度内插后，得到了基于强度的灰度图像(图 6(a)、图 7(a))。由于现在的计算机运算速度非常快，且计算过程中先对点建立了空间索引，两组实验的三角网内插及剖分均在2 min之内实现。根据道路的反射特性，选择合适的阈值进行分割，得到道路区域(图 6(b)、图 7(b))。由图 6(b)、7(b)可知，阈值分割只能粗略地将道路区域分割出来，许多由植被或其他地物引起的，在反射强度上类似于道路的噪声及其他非道路区域如停车场等都无法很好地剔除。十字剖分线法能够通过区域形状的狭长性及规则性，提取出真正的道路区域。图 6(c)、图 7(c)是根据十字剖分线法得到的点的权值图，可以发现，稳定、狭长状的道路区域点拥有较高的权值，而不稳定的噪声区域点和非狭长状的停车场的点则权值较低。根据点的权值高低，选取权值较高的点作为道路区域，权值较低的点则剔除(图 6(d)、图 7(d))。之后，大部分噪声和停车场已被剔除。去掉一些噪声后得到真正的道路区域(图 7(e)、图 7(e))。使用基于细化的道路中心线提取法提取出道路区域的中心线(图 6(f)、图 7(f))。对道路中心线进行矢量化放入原始数据中，如图 8所示。可以看出，第一组实验中，除了右上角的一条支路未被提取出以外，其他的道路均被正确提取。该支路未被提取出的原因在于，相对于该数据中的其他道路，该道路较短小，长宽比不大，即不具有狭长性。第二组实验中，除了一些非常小的支路，几乎所有的道路均被提取出来。这充分说明了本文提出的方法在LiDAR点云数据中进行道路提取的正确性。

图 6 基于十字剖分线的道路提取过程(数据A) Fig. 6 Road Extraction Based on Cross Profiles (Data A)

图选项

图 7 基于十字剖分线的道路提取过程(数据B) Fig. 7 Road Extraction Based on Cross Profiles (Data B)

图选项

图 8 最终生成的道路矢量线 Fig. 8 Road Vector Diagram

图选项

通过对原图像与提取结果的对比可以看出，本文算法较好地提取出了实验数据中的道路，并得到了道路矢量中心线，对于部分强度与道路区域相似的干扰区域也能很好地剔除。以上实验结果表明：利用道路和假道路区域在形状上的差异，可以有效地剔除干扰，提高提取结果的准确性和完整性。虽然能提取出绝大部分道路，但LiDAR强度信息噪声大，某些连着干扰区域的道路将一并删除，破坏提取结果的完整性。通过分析断裂处两端道路线方向走向及高程变化等关系，可以得出断裂处存在道路的可能性，进而修复断裂线。实验过程中也发现，一些道路会和干扰区域相连，形成一块大区域。故而在利用十字剖线统计区域形状并删除干扰区域时，会将此区域中的道路一并删除，造成提取结果中部分道路线的断裂。此外，由于道路中心线、接头处的处理来自于基于强度图像的整体提取，在细节上仍然需要几何上的后处理，如平滑、邻近阶段融合等。

本文通过基于区域形状分析的双十字剖线法得出的区域权值来区分道路与非道路区域。权值的高低体现了该区域的狭长性与规律性。通过比较基于强度信息粗提取的道路区域的权值，剔除其中高程及强度形状与道路区域相似的假道路区域。最后，通过对道路区域进行细化和平滑，提取道路区域的中心像素，并对其矢量化，得到道路中心线。本文重点解决了基于强度信息的道路提取过程中的假道路区域的检测问题，为进一步精确地提取三维道路网奠定了基础。但本方法并不适合提取具有大量遮挡的道路，一方面数据缺失需要修补，另一方面强度信号也被严重干扰。下一步工作集中于：(1)面对可能存在的道路区域反射强度不一致导致预处理过程中道路区域分割不完全的问题，需要寻找更好的方法将所有的道路检测出来；(2)在进行道路与非道路区域判别时，采用自适应生成的判别阈值；(3)考虑结合其他多源数据，可以提取更多的路面信息^[17]。