随着社会经济水平的发展，交通道路的更新频率也呈现加速之势，而建立现势性的道路网络是建设智能交通、数字城市、智慧城市的基本前提。传统的测绘手段成本高昂，且更新周期长，如何实现道路网快速有效的更新，已成为智能交通、智慧城市建设的关键问题之一。针对该问题，研究者从不同角度提出了两类方法，一是从高分辨率遥感影像中提取各类专题（包括道路网）信息^[1-2]；二是基于道路相关的其他类别数据（如公交轨迹），通过几何计算来建立道路网^[3-4]。通过遥感图像提取道路的研究可归纳为几何形态、层次特征及基于知识的提取方法。其中，基于道路几何形态特征的提取是比较流行的一种方法，即利用道路本身的几何、光谱、拓扑等特征来进行提取^[5]。此外，由于道路本身存在层次等级，可以借助影像分级从不同的粒度（对象层次）实现分层次特征的道路信息提取^[2]。针对道路信息复杂的情况，Chen等^[6]和Yin等^[7]提出了基于知识的分类方法，即利用现有的道路和其他辅助数据进行自动或半自动的道路提取。

尽管遥感影像为道路网更新提供了可行的途径，但受影像质量、影像分类算法耗时和精度以及部分高分辨率影像的价格因素影响，实践中该方法应用范围较窄，主要应用于一些重点区域。而基于免费地理数据，特别是众源地理数据，通过几何计算从中提取道路网的方法则受到更多关注^{[3, 8]}，如从出租车全球定位系统(Global Positioning System, GPS）轨迹数据中提取道路网，已是近年研究的热点之一^{[4, 8]}。

GPS轨迹数据记录了道路信息，从GPS轨迹中提取交通道路，按是否考虑轨迹点的上下文联系进行分类，可分为独立轨迹计算和轨迹上下文计算。独立轨迹计算方法将所有的轨迹点作为独立的个体对待，典型方法如轨迹点聚类方法^[9]、轨迹点密度方法^[10-12]等。这类方法基于轨迹点的距离分析，计算速度快，但因为没有考虑轨迹点关系，且要求轨迹点必须高密度分布，否则效果较差。而轨迹上下文计算充分考虑轨迹点间的时序和空间联系，通过对上下文关系的统计建模来构建交通网络，典型方法如轨迹合成算法、地图匹配（map-matching）方法等^[13-14]，其特点是充分考虑了轨迹的前后关系，对轨迹点的密度要求不高，但具有较高的计算复杂性，对数据精度要求高，否则很容易形成道路偏移现象，可能产生重复路段^[15]。文献[16]提出利用隐马尔可夫模型（hidden Markov model, HMM）提取失配轨迹的方法。该方法利用了状态转移概率模型，充分考虑了相邻轨迹点间的连接特征，利用历史道路网的几何网络功能，较精确地区分新旧道路上的轨迹点，能够实现增量道路更新，在轨迹数据量不是很大的情况下精度较好；但在提取失配轨迹点时，需要多轮循环，计算量很大，在大规模道路网更新中效率低下。

基于上述对历史道路网更新方法的分析，本文提出了一种新的道路网自动更新算法，即结合距离分析和GPS轨迹的上下文关系建模优势，将最新的GPS轨迹数据提取的道路合并到历史道路网中，完成道路更新。

1 道路网相关定义

道路网（road network, RN）是交通道路组成的具有网络拓扑结构的有向图。道路网的增量更新目标可描述为：某历史时刻点t₀的道路网表示为RN_t0，利用t₀~t_n间（t_n > t₀）采集的GPS轨迹数据，将RN_t0更新到t_n时刻，构建t_n时刻的现势道路网RN_tn。

定义1   轨迹线（trajectory）。用符号T表示，定义为由存在时间偏序关系的n个相邻轨迹点 < p₁, p₂…p_n > 构成的坐标序列，其中轨迹点p_i=(s_i, t_i)，包括空间位置属性s_i和时间属性t_i(1≤i≤n)两个基本内容，且t₁ < … < t_i-1 < t_i < t_i+1 < … < t_n。T中任意相邻的m个（m≥2）轨迹点形成T的一个轨迹段（trajectory section, TS）。

道路网眼（road network mesh, RNM）是由道路网中若干条路段围成的最小面状封闭区域。任意一条轨迹线T都可被相邻的若干个RNM切分成一个或多个首尾相连的轨迹段，轨迹段上的轨迹点都隶属于唯一的网眼（落在网眼内）；对任意一个RNM，落在该网眼内的所有轨迹段是轨迹分析的基本对象。图 1(a)显示了由若干个封闭RNM组成的一个区域，其中某网眼RNM_p由4条路段（R₁、R₂、R₃、R₄）形成。

定义2   失配轨迹点。任意一个轨迹点p_i可能隶属于历史道路网RN_t0，也可能属于新的道路。判断其能否匹配到RN_t0的过程，称为失配轨迹点判别（简称失配点判别）。若p_i不能够匹配到RN_t0，则称其不隶属历史道路网，即失配轨迹点。

定义3   候选点（candidate point, CP）。对轨迹T上隶属于RNM_p的一轨迹点p_i（即p_i在RNM_p内），定义一个距离阈值r，若p_i到组成RNM_p的任意路段的几何距离不大于r，则定义该道路上存在p_i的候选点，且为与p_i几何距离最近时对应的位置点。如图 1(b)所示，在距离阈值r下，组成RNM_p的路段中仅有R₁、R₄与p_i的几何距离小于r，且两道路上对应的cp_i¹、cp_i²点距离p_i最近，则cp_i¹、cp_i²点分别称为p_i在道路R₁、R₄上的候选点。

轨迹点p_i的候选点cp_i是该轨迹点在对应的历史道路上最可能出现的真实位置，若产生多个候选点，则与道路网数据精度和GPS定位精度有关。如果道路网数据精度足够高，采集的p_i的空间坐标位置足够精确，且p_i来源于现有道路，则该轨迹点能精确匹配到道路网RN_t0的某条道路上。由于数据精度原因，实际p_i与道路匹配点cp_i之间存在一定的位置偏差，即使轨迹点不落在道路上，也不能确定其为失配点，因此必须通过后续失配点判别计算来确定。为了减少计算量的同时保证计算结果的可靠性，常在预设的距离阈值r下寻找最多β个最近距离的候选点，称r为最远候选点距，β为最大候选点数。

2 失配轨迹点提取

失配轨迹点提取是从轨迹数据中构建新道路的关键信息，也是利用轨迹数据增量更新道路网算法的核心。本文提出的失配轨迹点的判别方法是基于几何距离分析和HMM建模方法的优势，利用轨迹点p_i所在的道路网眼距离及轨迹线关联上下文特征进行综合分析，判定其是否为失配点。若无法通过距离判别，则通过HMM建模计算。

2.1 轨迹点距离分析

轨迹点与道路网的几何距离是判断轨迹失配的重要依据。假设GPS采样误差符合正态分布，则某轨迹点p_j是否隶属于历史道路网，与该点到道路网的几何最短距离d_j密切相关。文献[17]对道路网眼内侧进行双缓冲，即预设两个距离阈值——最远候选点距r和最邻近距γ（r > γ），d_j > r，则认为轨迹点是失配的; 若d_j < γ，则轨迹点一定为非失配。该方法在进行双缓冲判别时，仅基于单个轨迹点计算，未利用轨迹点的上下文信息（如前轨迹点和后轨迹点），在个别轨迹点误差超过设置的阈值情况下，将对整个轨迹段造成错判。

本文在几何距离分析时，考虑了轨迹点上下文的距离关系，采用基于轨迹点上下文的综合距离判别方法，对于道路网眼内包含p_j的轨迹段TS_i={p_j-m, p_j-m+1…p_j…p_j+m-1, p_j+m}（m为邻接数，一般可取1或2，且j > m），TS_i集合内共计有大于等于2m+1个轨迹点。若有大于m个轨迹点满足几何最短距离d > r，则认为轨迹点p_j是失配的；同样，若有大于m个轨迹点满足d < γ，则p_j是非失配的。对未能通过几何距离判别的轨迹点，进一步采用HMM进行判别。

2.2 HMM建模

HMM为轨迹点失配性判别提供了一套概率统计方法。对于一个动态系统，可用转移概率表示从一种状态至另一种状态转换过程的可能性，若转移概率依赖于其相邻的前一系统状态，则称该过程为马尔可夫过程，可用马尔可夫模型来模拟。一般情况下，每个状态对于观察者来说是直接可见的，若该过程同时隐含未知状态时，则可用HMM的观测序列推算隐含状态序列，其产生观测结果的概率最大^[16]。HMM具有表现概率矩阵、状态转移矩阵和初始状态概率向量3要素^[18]。

在轨迹点的失配判别中，由于轨迹点p的观测坐标与其真实值存在误差，可假设该误差服从正态分布。同时假设p隶属于历史道路网上的某条道路，该道路上对应的候选位置为cp（隐藏状态），则包含了p的轨迹段对应的轨迹点系列形成观测状态序列，每个轨迹点对应的候选点cp形成的序列即为隐藏状态序列，通过判断隐藏状态序列（即cp候选点序列）能否产生得到的观测序列（产生概率）即可判断p的失配性。由轨迹点采样误差服从正态分布的假设，p点对应cp的表现概率（发射概率）可表示为^[16]：

${\rm{Pr}}({\rm{cp}} \to p) = \frac{1}{{\sqrt {2{\rm{ \mathit{ π} }}} {\sigma _z}}}{{\rm{e}}^{ - 0.5{{(\frac{{{{\left\| {{\rm{cp}} - p} \right\|}_{{\rm{great}}\;{\rm{circle}}}}}}{{{\sigma _z}}})}^2}}}$

(1)

式中，${{{\left\| {{\rm{cp}} - p} \right\|}_{{\rm{great}}\;{\rm{circle}}}}}$为cp点和p点的球面距离，一般可以用平面几何距离代替；${{\sigma _z}}$为GPS测量标准差，常选取15~25 m^{[16, 19]}。

状态转移矩阵由状态转换概率元素组成，假设p_i和p_i+1为道路网眼RNM内的轨迹段上相邻的两个轨迹点，对应隐藏状态候选点为cp_i和cp_i+1，通过正态分布计算隐藏状态转换概率Pr_b为：

$\begin{array}{l} {\rm{Pr}}\_{\rm{b}}({p_{i + 1}} \leftarrow {\rm{cp}}_{i + 1}^{s2}|{p_i} \leftarrow {\rm{cp}}_i^{s1}) = \\ \;\;\;\;\;\;\left\{ \begin{array}{l} (0,\frac{{||{\rm{cp}}_{i + 1}^{s2} - {\rm{cp}}_i^{s1}|{|_{{\rm{net}}}}}}{v} > t\\ {P_{\left( {{\rm{cp}}_i^{s1},{\rm{cp}}_{i + 1}^{s2}} \right)}},\frac{{||{\rm{cp}}_{i + 1}^{s2} - {\rm{cp}}_i^{s1}|{|_{{\rm{net}}}}}}{v} \le t \end{array} \right. \end{array}$

(2)

式中，${P_{\left( {{\rm{cp}}_i^{s1},{\rm{cp}}_{i + 1}^{s2}} \right)}} = \frac{1}{{\sqrt {2{\rm{ \mathit{ π} }}} }}{{\rm{e}}^{ - 0.5{{\left( {\Delta d} \right)}^2}}}$，且$\Delta d = {\left\| {{p_i} - {p_j}} \right\|_{{\rm{great}}\;{\rm{circle}}}} - ||{\rm{cp}}_i^{s1} - {\rm{cp}}_{i + 1}^{s2}|{|_{{\rm{net}}}}$；Pr_b为条件概率，表示在${p_i} \leftarrow {\rm{cp}}_i^{s1}$的前提下${p_{i + 1}} \leftarrow {\rm{cp}}_{i + 1}^{s2}$的概率，结合观测序列p_i和p_i+1，用于判断历史道路网上连通两个候选点的可能性（概率）；t为p_i到p_i+1的时间间隔；上标s1、s2分别表示p_i、p_i+1对应的第s1、s2个候选点；${||{\rm{cp}}_{i + 1}^{s2} - {\rm{cp}}_i^{s1}|{|_{{\rm{net}}}}}$为${{\rm{cp}}_{i + 1}^{s2}}$到${{\rm{cp}}_i^{s1}}$的路网距离；v为实际行车的平均速度。当候选点的路网距离的理论最小所需时间比t大时（即$\frac{{||{\rm{cp}}_{i + 1}^{s2} - {\rm{cp}}_i^{s1}|{|_{{\rm{net}}}}}}{v} > t$），表明不可能通过s1、s2候选点实现从p_i到p_i+1的连通，因此状态转换概率Pr_b为0。在存在可能的候选连接点的情况下，Δd用于计算观测点球面距离和候选点路网距离的差异，差异越小，则转换概率越大，由候选点序列实现观测点序列的可能性越高。该概率计算可由Viterbi算法^[18]实现。

应用HMM分析GPS轨迹点的失配性，其原理如图 2所示。其中，轨迹点p_i即为HMM中的观测（可见）状态，cp_i为对应的HMM的隐藏状态，为每一个观测轨迹点p_i依照表现概率所对应的若干个基于最短距离的最可能候选点cp_i¹, cp_i²…cp_iⁿ（n≤β）。计算每一个轨迹点p_i的候选点到下一个轨迹点p_i+1的候选点cp_i+1¹, cp_i+1²…cp_i+1ⁿ（n≤β）的转换概率（图 2(a)），文献[19]详细列举了4种轨迹线断裂的情形，图 2(b)显示了其中的两例，对于任意的p（p=1, 2…β），cp_i^p都无法经过cp_i+1^p通达至cp_i+2^p，即认为从p_i至p_i+1或p_i至p_i+2为失配段。

3 基于道路网眼的道路网更新算法

道路网更新算法通过对GPS轨迹数据的质量控制，考虑轨迹点上下文的距离分析及HMM建模，提取出失配轨迹点和轨迹段，将失配轨迹进行骨架线重构，实现道路重建和更新。

3.1 GPS轨迹点质量控制

对GPS轨迹数据进行质量控制，去除冗余轨迹点数据，形成符合逻辑的轨迹线，在历史道路网及GPS轨迹数据基础上，为生成道路网眼提供数据。

轨迹点质量控制以轨迹线为基本处理单元，对轨迹线内的轨迹点进行去冗余和轨迹分析，从而减少冗余数据量，并为失配轨迹点提取提供高质量的轨迹线数据。对于任意的轨迹线T，设p_i-1、p_i和p_i+1为T上的3个相邻轨迹点，对于给定的预设距离dis_a，若p_i-1和p_i+1的直线距离dis(p_i-1, p_i+1)≤dis_a，则p_i为冗余轨迹（即在p_i-1至p_i间或从p_i至p_i+1间可能存在停车状态，或途中GPS采样丢失），故将p_i删除。此外，对于T上任意相邻的两个轨迹点p_i和p_i+1，设dis_b为一个预设距离，若p_i和p_i+1间直线距离dis(p_i, p_i+1) > dis_b，则认为在p_i和p_i+1间存在行车GPS轨迹点丢失，不能将p_i和p_i+1作为一个轨迹线上的连续轨迹点，而需要将轨迹线T进行分割，生成两个独立的轨迹线并进行后续的失配轨迹分析。

3.2 道路网眼生成

基于GPS轨迹数据和历史道路网形成最大凸包面^[20]，形成的凸包面需要确保其空间上覆盖了历史道路网和所有GPS轨迹点数据。该凸包的边界和历史道路网（线）数据建立由道路网划分的封闭拓扑多边形提取最小面域，即道路网眼。图 3显示了建立道路网眼的过程，3条轨迹线（点）和4条道路线共同建立了一个最大区域面（最大凸面），将最大凸面的边界线与道路网合并进行面拓扑重建，形成最小封闭面，即道路网眼。采用最大凸面保证了生成的道路网眼覆盖整个研究区，从而所有的GPS轨迹点都有机会参与失配点判别。

3.3 失配轨迹提取

对于道路网眼RNM_j（j=1, 2…m）以及落在RNM_j中的k个轨迹段TS_ji（j=1, 2…m, i=1, 2…k），采用几何最短距离分析和HMM分析，提取出所有道路网眼内的失配轨迹段。在道路网眼RNM_j内，若p为轨迹段TS_ji内的任一轨迹点，且p为失配点，则无需对该轨迹段里的其余轨迹点进行失配检验，可直接判定TS_ji为完全失配段，从而极大地提高了失配轨迹的提取效率。设TS_j.={TS_j.}为RNM_j内的所有失配轨迹段集合，提取所有的道路网眼失配轨迹段，获取研究区失配轨迹段的总集合TSS={TS_j., j=1, 2…m}。

3.4 轨迹段骨架线重构

道路网眼内的轨迹段密度决定了新增道路的可能性，新增加的道路轨迹段覆盖密度一般较高。因此，可以通过轨迹段密度方法来实现对轨迹段骨架线的重构。本文采用轨迹段缓冲和密度判别的方法，具体方法是：先对网眼内发现的所有失配轨迹段进行双侧γ缓冲，将缓冲区进行叠加，计算轨迹段的积累密度，提取出轨迹段密度“山脊线”，作为新增道路。对“山脊线”的提取等同于分水线的提取，后者可以看作水流的起源点，常用的地理信息系统工具（如ArcGIS）提供了地表汇流（flow accumulation, FA）模块，通过FA模块可方便地计算出汇流量的空间分布，无汇流的区域即为分水线。

采用道路网眼作为算法分析的基本单元时，由于采样本身的稀疏性，当存在重构的骨架线任意一端无法连接到围成道路网眼的路段时，采用骨架线延长（长度γ）的方法，将其与网眼路段的交点作为融入历史道路网的结合点。若无法产生交点，则认为该新增道路为未接入历史道路网的悬挂路段。

4 实验与分析 4.1 数据来源及参数设置

本实验采用ECML PKDD(European Conference on Machine Learning & Principles and Practice of Knowledge Discovery in Databases)提供的数据集^[21]，具体为葡萄牙波尔图市（8.66°W~8.61°W, 41.18°N~41.22°N）的出租车轨迹。该数据集包括了2013-07-01至2014-06-30期间442辆出租车采集的78万条轨迹线、3 689万个轨迹点，相邻轨迹点的采样时间间隔为t=15 s，按照最大行驶速度120 km/h计算，相邻采样点间的最大距离应该不超过500 m，若大于500 m，则说明存在GPS轨迹点采集丢失，应对所在轨迹线进行切分处理。此外，作为志愿者地理信息的典型代表，OpenStreetMap（OSM, http://www.openstreetmap.org）数据可用来分析城市道路信息^[22]。为了与GPS轨迹数据配合，从OSM中提取了2013-01-01的历史道路数据。作为比对，也利用了2016-12发布的谷歌地图进行可视化结果显示。OSM道路网和谷歌地图比对后, 道路显示结果为基本匹配；与实际历史道路网相比，OSM存在部分道路缺失，但不会对实验结果产生影响，因此采用OSM数据是可行的。

按照§3.1的GPS轨迹点质量控制方法对车辆轨迹数据进行质量控制。定义dis_a=30 m，dis_b=500 m，对数据中的所有轨迹线进行分析。若相邻p_i-1、p_i、p_i+1点间的p_i-1和p_i+1的直线距离小于等于dis_a，则删除p_i。如果相邻的两个轨迹点p_i、p_i+1的间距大于500 m，则认为存在轨迹缺失问题。道路网眼的生成利用了OSM提取的城市道路及ECML PKDD轨迹数据，采用§3.2所述方法实现。几何距离分析中，设计失配轨迹点判别的邻接数m=2，根据文献[16, 21]，设置GPS测量误差标准差=20 m，最远候选点距r=100 m，最近距离γ=15 m，最大候选点数β=4。由于OSM数据集中没有提供道路行车速度数据，实验中的行车平均速度v设置为120 km/h，该取值保证了HMM中对失配点检测不会造成遗失（但增加了部分计算量）。

4.2 实验结果与分析

通过历史道路网和GPS轨迹点数据空间叠加，形成研究区的凸包，由历史道路网将研究区凸包分割成连续道路网眼，并覆盖所有轨迹点，形成了失配点提取的最小操作单元，一个典型的道路网眼如图 4(a)所示。对该网眼内所有轨迹段的轨迹点进行基于几何距离的失配点判别。首先，通过道路网眼γ和r内侧缓冲，分布标记网眼内的轨迹点距离特征；其次，对于任意路段，采用基于轨迹点上下文的综合判别方法。对于轨迹点p_i及其前、后相邻的2m个邻接轨迹点，判别p_i是否为失配点或非失配点，图 4(a)显示的是道路网眼内通过距离判别获取的失配轨迹点（三角形）的实例；对于无法通过基于轨迹点上下文距离判别的轨迹点，应用本文所提的HMM建模完成失配轨迹点提取（图 4(a)方形点所示）。提取所有失配轨迹点所对应的轨迹段，结果如图 4(b)所示；对失配段进行缓冲，形成空间密度分布图，对轨迹段进行骨架线重构（图 4(c)）；利用Douglas-Peucker算法^[23]将生成的骨架线进行平滑处理，形成连续但存在首尾悬挂端点的道路段（图 4(d)）；对悬挂道路段进行端线延伸和匹配，达到失配轨迹点到增量道路的转换以及与历史道路网的融合（图 4(e)）；在影像地图上进行验证，发现更新的道路网与实际情况基本一致（图 4(f)）。

本文提出的基于道路网眼的失配轨迹提取方法，综合利用了道路网眼内侧缓冲、轨迹点上下文的几何距离综合判别和基于HMM模型的方法来判别失配轨迹点和失配轨迹段。实验发现，在设置了最邻近距γ=15 m和最远候选点距r=100 m的情况下，通过内侧双缓冲（γ、r）及轨迹点上下文统计分析，计算得到的失配点占总量的绝大部分（90%），而通过HMM得到的失配轨迹点平均仅占10%。由此可见，大部分失配点都可以通过缓冲分析和对轨迹点上下文统计分析得到，相比单纯使用HMM提取失配点的方法，可极大减少计算量。

对于无法通过最短距离提取的剩余小部分轨迹点（约占10%），则利用HMM方法来实现失配点提取的准确性，本文方法的主要计算时间源于此。但与其他类似方法相比（如文献[16]），本文提出的失配轨迹获取方法是基于独立网眼单元进行计算的，即网眼单元可分别进行失配轨迹提取，更易于利用并行和分布式运算来实现，对于大规模的道路网数据能大幅度提高计算效率，实现实时计算。

5 结语

开源数据及GPS轨迹数据为增量更新道路网提供了经济可行的途径。采用增量更新道路网的优势是不必处理全部采集记录数据，而仅需和新增道路密切相关的轨迹数据生成最新道路网。本文首先提出了利用道路网眼控制下的GPS轨迹点更新历史道路网的方法，结合几何距离判别和HMM模型的优势快速提取失配轨迹点；在几何距离判别中利用了轨迹点与现有道路的几何距离以及轨迹点本身的上下文信息；然后在HMM判别中利用HMM的隐藏态转移概率、隐藏态到观测结果的表现概率，判别由隐藏状态（历史道路网的候选点）序列产生观测状态（GPS轨迹点）序列的可能概率，提取出失配轨迹点；最后基于失配轨迹的密度分布，通过道路骨架线提取方法，完成了新道路获取，并以此更新历史道路网。

基于道路网眼进行单元控制，充分发挥了距离分析和HMM提取失配轨迹的优势。HMM在提取失配轨迹时考虑了道路网的连通关系，对轨迹点的数量依赖性低，但计算复杂度高；而轨迹点上下文距离分析法计算简单，且在GPS点集密集的情况下具有高精度的特点。两者结合，相比现有方法（如文献[16]），可大量减少轨迹候选点的数量，降低计算量，保证提取的失配轨迹的完整性。

采用道路网眼的方法提取的失配轨迹更加连贯和准确，在道路网眼内通过轨迹段延伸，使提取的新道路能够很好地融合到历史道路网中。

以道路网眼为基本处理单元，使得算法可以支持并行计算，从而可以设计在线道路更新模式，为道路导航等应用提供实时的道路网更新结果。实验证明，本文方法在具有一定规模的GPS轨迹点区域具有良好的应用效果。

城市道路更新是数字城市和智慧城市发展和应用的基本要求，其更新的频率和准确性对包括自动车辆导航、城市公共服务应用领域都具有重要的应用价值。本文仅从GPS轨迹点和历史道路网入手，提出基于道路网眼的增量道路更新方法，未来将集成包括车辆轨迹点和其他社会感知网络（如人口流动）、高分辨率遥感影像等数据在内的多源空间数据，同时研究更智能化的数据融合分析算法和道路信息提取模型，以提高城市道路网更新的可靠性和实时性。