由于城乡道路建设飞速发展，道路信息快速获取与更新成为亟待解决的问题。传统道路数据采集方式包括野外测量、高分遥感影像等，其中高分影像以其高分辨率特征成为主要技术手段。然而，高分影像价格昂贵，难以开源获取，影像中的噪声干扰使得提取算法复杂。尽管高分影像提高了时间分辨率，但仍不能实时获取，难以满足人们对道路数据低成本、高实时性的要求^[1]。随着移动传感技术的发展和众源地理信息的出现，产生了海量的时空轨迹大数据。其中，车辆轨迹蕴含了丰富的道路信息，直接反映了道路网络的几何特征和运动状态^[2]，其实时、低廉、众源特性相比遥感影像等传统数据源更具优势。如何从海量、高噪音的轨迹数据中快速提取道路信息，成为交通地理信息领域感兴趣的问题。

目前，学者们围绕众源轨迹提取道路数据集中在道路中心线提取^[1-3]。依据几何模型和数学统计方法的差异，可分为4类：一是轨迹聚类方法，包括轨迹点、轨迹线聚类，如文献[4-5]运用K-Means聚类算法并结合高斯模型从车辆轨迹中提取道路中心线；文献[6]通过路口转向判断模型对轨迹分类，然后聚类轨迹提取路网；文献[7]根据轨迹线的空间邻近与方向关系聚类，提取路网数据。聚类方法适合高精度、低噪音轨迹数据。二是轨迹增量提取方法，如文献[8]基于力学思想，根据轨迹运动方向增量融合提取道路网；文献[9]基于轨迹点与候选路网之间的空间、语义关系来实现路网增量化生成；文献[10]提出了符合人类空间认知规律的轨迹线增量融合生成路网方法。三是基于图论的方法，文献[11-12]分别从GPS轨迹中提取道路交叉点、道路边，并基于图论的方法构建道路地图。四是栅格化方法，如文献[13-14]将轨迹数据栅格化，利用图像形态学方法提取道路面与骨架线；文献[15-16]将栅格化方法与聚类、机器学习、统计等方法相结合提取道路数据。以上模型方法适于道路中心线提取，难以识别道路边界。当前，道路边界提取仍以高分遥感影像^[17]、点云数据^[18]为主，但高分影像与众源轨迹相比成本高、实时性差，且提取算法复杂、专业要求较高，难以满足生产需求。从众源轨迹线集中提取道路边界的研究相对较少，其主要有两类方法：一是将轨迹数据栅格化^[13-14]，利用矢量化算法、自动矢量化工具如ArcScan提取道路面。二是图论方法，如文献[19]运用约束Delaunay三角网提取道路面域粗轮廓。以上方法都受轨迹密度差异影响大、难以精确提取道路边界，尤其对轨迹密度差异显著、道路结构复杂、不同时间跨度等特殊情形的道路边界提取适应性不强。

已有研究多运用高精度(采样间隔1～5 s)轨迹数据提取道路中心线，对于提取道路边界研究较少。其原因在于：①轨迹采样稀疏与高噪音问题对道路边界识别增加难度；②轨迹密度悬殊问题，导致轨迹高密度区域道路边界被拓宽，轨迹低密度区域道路难以提取。如何基于轨迹点(线)几何特征差异建立高效的道路边界识别模型成为该问题研究的关键。本文引入约束Delaunay三角网及其对偶Voronoi图表征轨迹密度变化差异，针对不同路面车辆分布频率悬殊、复杂道路网结构、多时间跨度等困难情形下的轨迹线数据处理建立针对性方法。

1 众源轨迹线几何特征分析

道路边界作为道路面内外的分界线，其两边的轨迹点密度差异显著，这种密度差异是识别道路边界的重要特征。要分析轨迹密度差异、识别道路边界需要解决以下3个难点问题。①高噪音特征。由于GPS信号漂移、采样间隔稀疏、道路弯曲复杂度^[6]等导致轨迹数据高噪音、缺失等问题^[19]，使得无法直接从原始轨迹中提取道路数据。②不同道路结构下的密度差异。在城市主干道、交通枢纽区域轨迹点密集，而支干路轨迹点稀疏,这种密度差异导致提取的主干道边界被拓宽，支路无法完整提取。特别对于复杂道路结构，如环形道路交叉口、平行道路结构，道路间间隔越小，道路边界越模糊，道路与非道路区域难以区分，有效处理复杂道路结构下的轨迹密度差异性问题是本文研究的目标之一。③不同时间尺度下的密度差异。轨迹作为时间序列数据，不同时间长度下的数据量具有质的变化。时间尺度过小，轨迹数据量小，则提取的道路面域小于真实道路面；时间尺度过大，轨迹数据量大，提取的道路面大于真实道路面。顾及行车频率小的低等级道路，需要选取时间跨度长的轨迹，然而对于高等级道路会造成轨迹点过于密集以至于邻近道路混为一团，难以区分其间的边界，选取多长时间段的轨迹线作为道路提取的数据对象是一个难点问题。

针对众源轨迹线提取道路信息的难点问题，本文重点关注轨迹密度变化差异，引入约束Delaunay三角网及其对偶Voronoi图模型，计算三角形边长度与Voronoi面域大小等几何特征表达轨迹点密度变化差异，建立众源轨迹点(线)边界识别模型和道路边界提取方法。

2 Delaunay三角网支持下的轨迹数据提取道路边界

从轨迹线中提取道路边界的核心是识别轨迹点集分布的边界^[20]，本质上是一个空间邻近分析问题。Delaunay三角网是空间邻近分析的有力工具，广泛用于空间数据聚类^[21]、空间邻近关系探测^[22-23]。故本文运用约束Delaunay三角网(constrained Delaunay triangulation，CDT)计算轨迹密度变化率、边长距离，建立道路边界识别模型，提取道路边界。

2.1 道路边界识别指标计算

(1) 密度变化率指数。道路边界作为道路面内外的分界线，两边的轨迹点密度差异显著，密度变化率成为道路边界识别的重要指标。文献[23]认为每个点目标的存在都要获取其生存空间，相邻点目标对分布空间的竞争结果相当于对空间做等分式划分(不考虑点性质差异)，并用点的Voronoi图面积的倒数表示该点的密度。在由轨迹线生成的Voronoi图中(图 1(a))，道路内部Voronoi图面积小，非道路区域Voronoi图面积大，道路边界两侧的轨迹点密度差异大。由于Delaunay三角网与Voronoi图互为对偶，故将轨迹密度变化差异融入Delaunay三角形边中。三角形边将空间分为左右区域，如果三角形边的左右两侧轨迹点密度差异大，则为道路边界区域。因此，三角形边相对左右两点的密度之比即是该边的密度变化率，该比值称为密度变化率指数(density change rate index，DCRI)，计算公式为

式中，D(p_left)、D(p_right)分别表示三角形边E_i相对的左、右两点的密度，DCRI的值域为(0,1］。根据轨迹点的密度(图 1(a))和DCRI值分布特征(图 1(c))，在道路内部和非道路区域轨迹密度差异小，三角形边DCRI值大；道路边界区域轨迹密度差异大，DCRI值小；DCRI值越小，轨迹密度差异越大，是道路边界的可能性越大。设定DCRI阈值(用DCRI_Value表示)，将三角形边分为两类(图 1(c))，如果三角形边DCRI≤DCRI_Value，则该边位于道路边界区域，可作为道路边界，简称突变边，用De表示；反之称为普通边，用Pe表示。

(2) 边长距离。Delaunay三角网中边的长度表达了点目标间在空间分布上的远近关系，边越长轨迹稀疏，密度低；反之轨迹密集，密度高(图 1(b))。由于车辆轨迹沿道路网聚集分布，在由轨迹线构建的Delaunay三角网中(图 1(b))，道路内部区域三角形边长度小，非道路区域边长度大。设定边长阈值，将三角形边分为两级(图 1(b))，即可区分道路与非道路区域，识别道路边界。根据Delaunay三角网中边长的统计特征^[22]，计算公式如下

式中，Len_Mean(DT)表示三角网DT平均边长；Len_Variation(DT)表示三角网边长变异；α表示调节系数。α值越大(默认为1)，整体约束越宽松，反之则越严格。根据边长度分布特征(图 1(b))，如果边长度Len≥Len_Value，可作为道路边界(图 1(e))，简称长边，用Le表示；反之，称为短边，用Se表示。

上述两个指标从不同侧面反映轨迹数据的密度变化差异和道路边界特征，但由于道路网的复杂性与轨迹数据分布差异性，使得单一指标具有不同的控制领域与适应条件。如图 1(d)所示，DCRI指标确定的突变边界在空间分布上具有混合性、不均匀性、不连续性，使得突变边识别的道路边界不完整、道路边界拓扑断开。尤其对于轨迹点稀疏的支干道路边界识别不完整(图 1(d)中A)。DCRI指标适合于轨迹数据密集、密度变化差较大的情形，如主干道(图 1(d)中B)、道路交叉口(环形道路、转盘)等复杂道路结构。边长指标适于轨迹稀疏的支干路(图 1(e)中A)、道路结构简单下的轨迹数据，这能解决DCRI指标提取道路边界不完整的问题。但对于轨迹密集、道路结构复杂的情形，用单一边长指标难以顾及全局与局部的密度差异，提取的边界会出现道路面被拓宽(图 1(e)中D)、环形闸道与空白区域无法区分的问题，而DCRI指标可解决该问题。故本文将两个识别指标融合集成构建道路边界识别模型，用于道路边界提取。

2.2 基于Delaunay三角网的道路边界提取

2.2.1 轨迹数据预处理

运用文献[18]中方法对原始轨迹进行常规(包括时间错误、越界等)清洗并生成轨迹线。由于轨迹稀疏采样，直接对原始轨迹线构建约束Delaunay三角网，则破坏三角网邻近特性。故本文根据轨迹转向角与轨迹段长度自适应加密轨迹线解决该问题。按照文献[24]中方法计算每个轨迹点的转角θ和该轨迹点相邻轨迹段的长度L。通过试验表明，如果θ为(30°,90°)且L>300 m，则该点相邻轨迹段加密多、加密步长较小，反之加密少、加密步长较大。

2.2.2 道路边界提取方法

根据2.1中所述，本文集成边界识别指标构建道路边界识别模型，如图 2(a)所示，对于每一条三角形边，如果是长边(Le)或者突变边(De)，该边为道路边界，称为“障碍边”。根据“障碍边”的数目将Delaunay三角形分为4类：只有1条障碍边的为Ⅰ类三角形；有两条障碍边的为Ⅱ类三角形；有3条障碍边的为Ⅲ类三角形；没有障碍边的为Ⅳ类三角形(图 2(b))。道路边界提取算法思想为：在保持拓扑连通性的前提下，将任意Ⅳ类三角形作为“种子点”，运用“种子点”区域增长算法寻找多边形范围。由与种子点相连的任意一个三角形出发，向3方向扩展，对于一个三角形来说，搜索路径为一边进入、两边输出，因此采用二叉树广度优先遍历的方法遍历三角网，一旦输出边为障碍边,则停止该障碍边方向上的搜索。Ⅱ类三角形可看作叶子节点，只有一边进入，没有输出；Ⅰ类三角形是拥有一个孩子节点的非叶子节点；Ⅳ类三角形是拥有两个孩子节点的非叶子节点；Ⅲ类三角形为障碍区，不搜索。

根据边与三角形的邻接关系持续扩展直到扩展边均到达障碍边,所有三角形范围构成了封闭区域,完成道路多边形生成(图 2(c))。由于真实环境中路网结构复杂、轨迹密度差异显著，在算法实际应用时还需要根据三角形类型、上下文关系来决策确定道路边界(障碍边)。

(1) 孤立障碍边与障碍区。由于道路内部、不同等级道路连接处的轨迹密度差异，在搜索扩展中出现零散分布的孤立障碍边、障碍区(图 3(a)中C)。障碍区由少量Ⅱ类三角形、Ⅲ类三角形组成，形成了道路多边形区域内的障碍区“空洞”。这些孤立障碍边、障碍区的邻接区域都已扩展搜索，处于活跃状态。对于该情形，把障碍边当作非障碍边继续扩展。对于障碍区还需要设定面积阈值与环形道路中的非道路区域区分开。

(2) 扩展“瓶颈”与边界平滑。在搜索过程中出现单一障碍边(图 3(b)中A)、单一Ⅲ类三角形(图 3(b)中B)阻断扩展路径，无法完成搜索。对于单一障碍边，如果障碍边左右邻接的三角形都处于扩展活跃状态，需跨越该障碍边继续扩展。对于Ⅲ类三角形，如果该Ⅲ类三角形的3个邻接三角形中有两个以上三角形处于扩展活跃状态，将与活跃状态三角形邻接的边作为非障碍边继续扩展。另外，城市支路上出现障碍边与非障碍边间隔分布(图 3(c))情形。针对该情形，判断每个Ⅱ类三角形障碍边的邻接三角形类型。如果是Ⅱ类、Ⅰ类、Ⅳ类，则跨越障碍边继续扩展；如果是Ⅲ类三角形，则以障碍边为边界终止扩展。扩展完成后的道路边界凹凸不平(图 3(b))，需对道路边界平滑化简，最后完成道路多边形提取。

2.3 指标参数取值及探讨

在道路边界提取过程中，选择合适的Len_Value(调节系数α值)与DCRI_Value阈值显得尤为关键。根据2.1节中对参数值分布规律的分析，参数取值即是找出那些数量少的长边Le和突变边De。本文根据参数值分布规律，利用ArcGIS自然间断点分级法进行取值试验。分析不同参数取值的试验结果精度，发现DCRI最佳取值范围为(0.2,0.6)，边长调节系数α最佳取值范围为(0.5,2)。在各自指标的最佳取值范围内，取不同的组合值进行试验，分析组合取值对道路提取结果精度的影响，即参数取值的敏感度分析(图(4))，其敏感度变化符合正态分布曲线，如图 4(c)所示。在图 4(c)中，当两个参数组合取值时，其最佳取值区间有微小的变化，DCRI变为[0.3,0.5]，α为[0.5,2]，在该范围内的结果精度达可到60%左右。当指标取值位于(DCRI=0.4,α=1)附近时，结果精度到达80%以上。在最佳取值范围内，通过调节、组合参数取值使其适应不同情形下的道路边界提取。

3 试验与分析 3.1 数据与试验环境

车辆轨迹数据为北京市2012年11月3日(星期日)一天的出租车轨迹。出租车轨迹数据包括车辆ID、GPS时间、GPS经纬度等信息。轨迹点采样间隔为10～120 s不等，共有轨迹点32 882 436个，生成轨迹线为144 205条，预处理后为45 768条，本试验选取北京市望京社区一定范围内作为试验区域(图 5)。本试验在P4/8GB/2GB/Win7环境下，基于ArcGIS平台、C#编程语言进行道路面、线数据提取的算法实现。

3.2 试验结果分析

对预处理加密的轨迹线构建约束Delaunay三角网，计算DCRI与边长阈值，根据道路边界识别模型，运用“种子”点区域增长算法提取道路边界，生成道路面域数据(图 5)。道路面域二次构建Delaunay三角网，提取道路中心线(图 5)。将提取的道路数据与影像地图、OSM电子地图叠加进行定性评价。如图 5所示，道路面、道路线数据的整体结果基本覆盖试验区的道路，并且准确度较高。如图 5(b)、(e)，本方法对于复杂道路交叉口(样区A)区域能较好区分环形闸道中的空白区域。如图 5(c)、(f)，本方法能够处理不同道路等级(样区B)下的轨迹密度差异，区分空间上邻近道路的道路边界。但在城市内部的部分低等级道路面无法提取或提取精度低(图 5(c))，原因在于轨迹线太少。

以北京市标准道路矢量数据为参考，将本文结果与轨迹数据栅格化后利用ArcGIS矢量化的结果对比分析。采用多边形叠置方法评价道路面提取的有效性，并借鉴文献[3]中方法计算两种常用评价指标，即准确率(precision)和召回率(recall)来评价叠置结果。准确率反映试验结果的准确度，召回率反映试验结果的完整度。对于道路线数据，采用文献[25]中提出的缓冲区检测方法与文献[4]中方法得到的试验果对比评价，即分别建立2 m、5 m、8 m半径的缓冲区进行道路长度的准确率、召回率统计比较。结果统计评价如表 1。

表 1中P表示准确率,R表示召回率。从表中可看出，本文方法提取的面域数据在准确率、完整度上比栅格化方法提高10%以上。其主要原因在于栅格化方法对轨迹密度差异适应差，低密度区域道路无法提取(图 6(a)中A、B处)；高密度区域的结果大于真实道路面、邻近道路边界区分度不高(图 6(a)中C、D处)。由于GPS定位误差，不论哪种方法的提取结果都存在误差，难以与真实边界完全重合(图 6(a)、(b))。基于道路面生成的道路线数据在2 m缓冲区的精确度比栅格化方法有约40%的提高，完整度有约35%的提高。在拓扑正确性方面，栅格化方法出现较多孤立的道路面、线拓扑断开、道路边界z字形问题(图 6(a))，而本文的“种子点”扩展方法解决了该问题。

本文针对不同等级道路(不同行车频率)、复杂道路交叉口、不同时间跨度等3种特殊情形下的轨迹数据进行试验分析：

(1) 不同道路等级。城市区域不同等级道路具有不同交通流量、不同行车频率，导致轨迹密度悬殊。如图 7(a)所示，本文方法相比栅格化方法对于空间上邻近的双干道(图 7(a)中A)、低等级道路(图 7(a)中C)识别更完整。栅格化方法对密度低值区域无法提取或提取不完整(图 7(a)中C)，导致边界不平滑或拓扑断开，对高密度区域如双干道边界、环形匝道区域(图 7(a)中A、B)识别不如本文方法。

(2) 复杂道路结构。道路交叉口区域轨迹密集、数据量大，环形匝道、主干道等不同等级道路在空间上邻近，导致栅格化方法难以完整提取环形匝道内空白区域边界(图 7(b)中A、B处)。本文方法在边界提取中顾及了轨迹密度差异性，对于环形闸道与空白区域的区分较好(图 7(b)中边界)、边界提取更完整。

(3) 不同时间长度。分别选取1小时、1天、1周不同时间长度的轨迹线进行试验分析(图 7(c))。当处理1周轨迹数据量时(图 7(c)中D)，栅格化方法已经无法区分出环形道路中的空白区域，将空间上邻近的多条道路识别为一条道路。本方法对于1周数据量不进行轨迹线加密、选取约束性严格的参数值(DCRI=0.58，α=0.4)，提取的边界精度可达到73.1%，道路中的环形空白区也能较好识别(图 7(c)中D)。对于一条道路，轨迹量太少或太大(时间跨度过长)，即使通过调节参数也难以达到理想结果。

综上，本方法相比栅格化方法更适于从不同行车频率、不同时间跨度、复杂道路结构下的轨迹线集中提取道路边界。

4 结 论

本文针对众源轨迹线集提取道路边界问题，引入约束Delaunay三角网及其对偶Voronoi图几何模型，通过DCRI和边长距离表征轨迹密度变化差异，将两种不同几何维数的控制条件集成构建决策模型，运用“种子点”区域扩展方法实现道路边界的精确提取。通过北京市出租车GPS轨迹数据试验分析，证明了该方法的有效性。本研究贡献包括：运用Delaunay三角网及Voronoi图几何模型建立了基于轨迹点(线)几何特征差异的道路边界识别模型，运用决策模型、“种子点”区域扩展算法实现道路边界提取。该方法更适于轨迹密度差异悬殊情形下的道路边界提取，解决了从众源轨迹中提取道路边界问题。

同时，仍还有诸多内容需深入完善：①本文仅提取了道路几何数据，还需要提取道路属性数据如路名、车道、交通规则等，如何将多源异构的时空轨迹数据融合提取道路基础设施，感知城市交通动态将是下一步的研究重点。②本文方法对于复杂道路交叉口、立交桥等具有立体三维层次的道路面、线识别提取还需要作进一步研究。