近年来, 我国已经成为全球电动自行车生产和销售第一大国, 同时, 电动自行车也已经成为人们出行的主要交通工具之一. 据统计, 2018年我国电动自行车的保有量已经达到了2.5亿^[1], 随着GPS定位技术的不断成熟, 越来越多的电动自行车被安装了GPS设备, 由此获取到的海量轨迹数据为分析城市交通状况和用户出行规律提供了很好的基础^{[2, 3]}. 地图匹配算法的任务是将GPS记录的轨迹匹配到交通工具实际经过的道路上, 这是对轨迹数据进行深度分析和有效利用的必要步骤. 但是, 电动自行车上普遍使用的成本低廉的GPS传感器无法提供高精度的定位, 因此往往无法从原始的轨迹数据中直接得知实际经过的准确路线.

目前, 学者对于地图匹配的研究已经有20多年的历史, 专注于解决地图匹配问题的文献也有上百篇之多^[4]. 然而, 现有的各类算法大多是研究机动车轨迹数据的匹配问题. 对于电动自行车轨迹的地图匹配算法的相关研究还比较少. 相比于机动车的轨迹数据, 电动自行车轨迹数据具有以下3个特点.

特点1: 电动自行车轨迹上存在大量停留点.

特点2: 电动自行车轨迹采样频率高, 行驶速度较慢, 导致采样点距离较小, 轨迹点密度相应增大.

特点3: 电动自行车行驶路段更多, 由于路网精度的限制, 部分轨迹可能无法匹配到对应的路网上.

具体而言, 特点1会导致地图匹配准确率的下降. 停留点指轨迹上一组连续且相距较近的轨迹点, 其产生原因分为两类: (1)由于停止或运动速度极慢; (2)由于车辆反复行驶过某区域^[3]. 这两类原因均会导致文献[4]中提到的不必要迂回的出现, 从而降低匹配的准确率. 例如, 在图1中, 当车辆在轨迹点 $ {p_2} $ 与 $ {p_3} $ 处行驶速度极慢, 为两个停留点, 由于GPS的测量误差^[5], $ {p_2} $ 与 $ {p_3} $ 在道路A上的投影点前后顺序与道路A的方向相反, 导致地图匹配算法计算出 $ {p_2} \to {p_3} $ 的路线为道路A→道路B→道路C→道路D→道路A, 与事实情况不符.

文献[6]对GeoLife^{[7, 8]}轨迹数据集中不同出行方式的轨迹进行了分析, 认为自行车与步行出行方式的轨迹数据有更高的方向变化概率与停止概率. 在真实场景中, 机动车只有当等待红绿灯或处于交通拥堵路段时才会导致轨迹上出现停留点, 但是电动自行车用户拥有更大的行驶自由, 例如可以随意在路边停车、在某较小区域内反复行驶, 从而产生更多的停留点. 因此对电动自行车轨迹上的停留点进行预先识别与处理十分必要.

特点2会显著地降低地图匹配算法的效率. 电动自行车车速较慢, 因此轨迹点更加密集, 高峰时段下轨迹点密度会进一步提高. 文献[9]将采样时间间隔在1 s至10 s之间的采样方式定义为高频采样, 而时间间隔高于2 min的采样方式定义为低频采样. 然而, 这种定义方式只考虑了时间间隔, 忽略了采样点在空间上的距离. 当交通工具速度不同时, 即使采样频率相同, 也可能导致采样点空间距离相差很大. 图2(a)显示了一段郑州电动自行车轨迹, 相邻轨迹点平均间距约为60 m、采样时间间隔约为12 s, 图2(b)显示了一段CRAWDAD数据集中旧金山出租车的轨迹^[10], 相邻轨迹点平均间距约为600 m、采样时间间隔约为30 s. 尽管电动自行车的采样频率约为机动车的3倍, 但是轨迹点密度却达到了机动车的10倍. 在许多开源轨迹数据集中, 机动车轨迹上的轨迹点密度相对较低. 例如, 开源数据集T-Drive^{[11, 12]}采集了北京市的一万余辆出租车一周内的轨迹数据, 平均采样时间间隔177 s, 轨迹点间隔623 m, 文献[13]中使用了2010–2011年间采集的北京市6万余辆出租车轨迹数据, 平均采样时间间隔70 s, 轨迹点间隔560 m.

鉴于大量的轨迹数据会带来极大的存储和计算开销, 许多研究人员关注于如何压缩轨迹以提升地图匹配速度^[14]. 基于DP算法^[15]的轨迹简化方法应用最为广泛, 其基于下述假设: 轨迹点的位置记录是交通工具在该时刻位置的真实记录. 由于GPS测量存在误差, 当交通工具位于高层建筑旁、桥梁旁时, 误差会急剧增大^[5]. 基于DP的轨迹简化算法会计算出与原始轨迹形状差异最小的简化轨迹, 因此倾向于保留GPS误差较大的轨迹点^[4], 进而可能导致匹配准确率下降. 例如, 图3中高层建筑旁的轨迹点GPS误差很大, 若使用简化后的轨迹进行地图匹配, 将会得到错误匹配结果. 由于卡尔曼滤波算法能够用于最优状态估计, 因此使用卡尔曼滤波算法与轨迹简化算法相结合, 能够在简化过程中不断修正轨迹点. 同时, 现有基于DP算法^[15]的轨迹简化算法^[16–20]依赖于人为设置的距离阈值, 只有简化过程结束后才能计算出轨迹的简化比例.

对于电动自行车轨迹数据的特点3, 我们使用无效轨迹片段来表示轨迹上不存在于给定路网的部分. 当路网精度降低时, 电动自行车的无效轨迹片段会相应增加. 在实际的路网中, 电动自行车能够在一些机动车不能驶入的路段行驶. 例如, 在图3中, 星形轨迹点不存在于图中的路网上, 这些轨迹点组成了一个无效轨迹片段. 一些全局匹配算法使用折线间距离度量方法如弗雷歇距离来计算路网上的各条可能路线与待匹配轨迹的相似度^[21]. 当轨迹上存在无效轨迹片段时, 这些算法仍会将相似度最大的路线作为匹配结果. 基于状态转移模型的地图匹配算法^{[14, 22–25]}则会因无效轨迹片段的存在导致匹配过程中断^[4], 无法对轨迹的有效部分计算出匹配路线. 上述两类算法都会因无效轨迹片段的存在而使匹配准确率降低.

为解决现有地图匹配算法因电动自行车轨迹数据的上述特点导致的匹配准确率下降、匹配效率不高的问题, 本文提出了一种可自适应路网精度的电动自行车轨迹地图匹配方法KFTS-AMM. 该方法融合了基于卡尔曼滤波算法的轨迹简化算法(KFTS), 以及分段的隐马尔可夫模型(hidden Markov model, HMM)的地图匹配算法(AMM). 本文的贡献主要集中在以下3个方面.

(1) 提出了基于时空聚类的停留点处理算法, 能够准确地识别与合并停留点, 减少轨迹上的迂回, 从而提高地图匹配算法的准确率.

(2) 提出了基于卡尔曼滤波算法的轨迹简化算法KFTS, 改进了DPhull算法, 使其不依赖人为设置的距离阈值, 仅使用轨迹简化比例作为简化过程参数. 由于轨迹简化算法倾向于保留误差较大的轨迹点^[4], 因此在轨迹简化过程中加入滤波操作对轨迹点加以修正, 能够提升后续地图匹配过程的准确率.

(3) 使用分段的HMM的自适应地图匹配算法AMM, 能够在匹配的过程中不断检测出无效的轨迹片段, 同时将有效的轨迹片段匹配到给定精度的路网图上, 在不同精度的路网上均能取得良好的准确率.

本文第1节回顾了地图匹配算法、轨迹简化算法与轨迹停留点识别的相关研究工作. 第2节介绍相关的定义以及问题描述. 第3节详细介绍了本文提出的算法. 第4节展示了在真实数据集上的对比实验结果. 第5节总结了现有工作以及未来工作的重点.

1 相关工作

在本节中, 我们回顾了地图匹配、轨迹简化以及轨迹停留点识别的相关研究工作.

1.1 地图匹配

文献[26]提出根据轨迹匹配过程中的采样点范围将地图匹配算法分为局部/增量算法^[27–29]与全局算法^{[21, 30]}. 局部/增量算法使用轨迹的局部特征进行匹配, 计算速度快, 适用于高频率采样的轨迹, 但当路网密集度较高时, 准确率普遍较低. 全局算法使用整条轨迹进行匹配, 为轨迹寻找一条相似度最高的路线, 计算复杂度较高.

Quddus等人根据地图匹配算法所使用的技术将地图匹配算法分为以下4类^[26]: 基于几何关系的匹配算法^[31]、基于拓扑关系的匹配算法^[32]、基于概率统计的匹配算法^[33]和其他匹配算法(例如: 使用了诸如扩展卡尔曼滤波器^[34]、模糊逻辑^[35]、证据理论^[36]和贝叶斯推理^[37]等技术). 但是, 随着近年来又有许多新的方法被提出, 这一分类标准已经不再适用. Chao等人在文献[4]中提出根据核心匹配模型将地图匹配算法分为4个类别: 相似度模型^{[21, 27]}、状态转移模型^{[9, 14, 23, 24, 38]}、候选进化模型^[39]和评分模型^[40]. 同时, Chao等人还列举了3种给地图匹配带来挑战的数据质量问题, 包括不必要的迂回、异常点和路网中道路的密度, 这3种问题均会降低地图匹配算法的准确率. 在基于状态转移模型的地图匹配算法中, 基于隐马尔可夫模型的地图匹配算法^{[14, 22–25, 38, 41]}是最流行的, 同时也被证实具有较高的准确率^[42].

一些研究人员致力于提高地图匹配效率. 地图匹配效率的提升可以通过以下4种方式^[14]: 使用空间索引技术^{[43, 44]}以加快对某个轨迹点的近邻点与近邻边的查找速度、避免路网图中最短路径的重复计算^{[14, 45]}、使用分布式与并行计算技术^{[22, 46–48]}、压缩轨迹以减少参与计算的轨迹点.

上述地图匹配算法设计的实验所使用的数据集多为机动车轨迹数据集, 由于电动自行车轨迹上存在大量停留点、采样频率高、轨迹点密集且存在大量无效轨迹片段的特点, 已有算法不适宜直接应用在电动自行车轨迹数据上.

1.2 轨迹简化

轨迹简化算法是目前最主流的轨迹数据压缩处理方法^[49], 将轨迹数据视为一条连续的折线, 通过删除部分对折线形状影响较小的轨迹点, 实现简化后折线与原始折线的近似拟合, 从而完成轨迹数据的压缩. 根据轨迹简化算法是否适用于实时计算, 可以将现有轨迹简化算法划分为在线轨迹简化算法与离线轨迹简化算法^[2].

离线轨迹简化算法没有输入轨迹规模的限制, 能够得到全局最优的简化结果. Bellman算法^[50]是最早用于解决轨迹简化问题的算法, 利用了动态规划技术, 但是计算复杂度较高, 达到了 $ {\rm{O}}({n^3}) $ . DP算法^[15]是另一种经典的轨迹简化算法, 在最坏情况下的时间复杂度为 $ {\rm{O}}({n^2}) $ . 许多研究人员为了提升DP算法的效率提出了改进方法, 如DPhull算法^[16]引入了凸包技术, 使算法在最坏情况下的时间复杂度下降到了 $ {\rm{O}}(n\log n) $ . 现有基于DP算法的轨迹简化算法大多依赖于人为设置的距离阈值, 无法直接确定最终的简化比例. 文献[51]提出了不使用距离阈值作为参数的简化方法, 使用递归的方式, 动态调整距离阈值以控制简化比例, 但是易造成轨迹上不同片段简化比例的不平衡. 文献[52]提出的简化算法首先依据速度对轨迹分段, 再结合DP算法对轨迹简化, 以保留轨迹时空特征.

在线场景中, 需要实时对轨迹数据进行简化, 因此必须根据计算设备的存储与计算性能限制轨迹的输入规模. Sliding window (SW)算法^[17]和opening window (OW)算法^[18]都使用了能够动态调整大小的滑动窗口作为轨迹点缓冲区, 利用两个轨迹点作为起止点标记一个滑动窗口, 当窗口内某轨迹点与窗口首尾轨迹点连线的欧式距离超过设定阈值时, 将窗口的起始轨迹点向前移动. 这两种算法不同之处在于, SW算法将起始轨迹点移动到终止轨迹点的前一个轨迹点处, 而OW算法则将其移动到距离首尾连线距离最大的轨迹点处. STTrace算法^[19]与SQUISH算法^[20]设置了固定大小的缓冲区, 同样根据窗口内各轨迹点距首尾轨迹点连线的距离大小来决定是否将窗口向前滑动.

1.3 停留点识别

停留点是指在某一段连续的时间内, 位置变化小于一定距离的一组轨迹点. 当车辆运动速度过于缓慢(例如用户推行电动自行车)、暂时停止行驶, 或在较小区域内反复行驶时, 就会产生停留点.

停留点识别方法有3类: 基于差异的方法^{[8, 53]}、基于历史数据与额外位置信息的方法^[8]和基于聚类的方法^{[54, 55]}. 其中, 基于差异的方法需要通过设置速度、距离或者时间阈值来识别停留点, 这种方法在拥有足够先验知识的情况下才能够有较高的准确率; 基于历史数据与额外位置信息的方法通过用户历史出行数据挖掘用户的出行模式, 进而结合关键位置信息推测停留点; 基于聚类的方法使用聚类算法寻找轨迹点聚类, 但是这种方法同样依赖输入的参数, 如距离和时间阈值. 当距离阈值与时间阈值设置合理时, 基于聚类的方法在没有历史数据和位置信息的情况下仍能达到较高的准确率. 文献[8]设置了时间与距离阈值, 线性扫描轨迹以识别停留点, 然后以求均值的方式合并停留点. 然而这种识别方法只考虑了相邻位置的时空关系, 当轨迹在一段区域内反复迂回, 则无法有效识别出停留点.

DBSCAN算法^[56]是一种基于密度的聚类算法, 能够识别任意形状的聚类, 且不需要预先确定聚类个数. 但是由于时空数据的时间属性和空间属性尺度不同, 因此该算法不能很好地解决时空数据的聚类问题. 为此, Birant等人提出了ST-DBSCAN算法^[57], 该算法使用了空间邻域半径和时间邻域半径作为聚类内轨迹点的距离阈值, 能够很好地解决时空数据的聚类问题. 本文StayPointsProcess算法即使用ST-DBSCAN算法完成停留点的识别.

1.4 卡尔曼滤波

卡尔曼滤波算法^[58]是一种利用线性系统状态方程, 迭代地对系统每一时刻状态做出最优估计的算法. 许多研究人员对其进行了创新, 使卡尔曼滤波算法也可应用于非线性系统, 例如使用泰勒级数展开将非线性系统线性化的扩展卡尔曼滤波算法(extended Kalman filtering, EKF)^{[59, 60]}和结合了无迹变换(unscented transform, UT)的无迹卡尔曼滤波算法(unscented Kalman filtering, UKF)^[61]. 卡尔曼滤波算法在许多领域都有巨大的应用价值, 例如车辆状态预测^[62]、动态目标跟踪^[63]、传感器信息融合^[64]等. 卡尔曼滤波算法能够融合多种传感器数据, 减小仅使用GPS传感器对车辆定位时的误差. 文献[34]使用车内里程表、陀螺仪与GPS传感器数据, 通过卡尔曼滤波算法构建车体的运动模型, 然后结合地图匹配算法对车辆精确定位. 文献[65]提出了将卡尔曼滤波算法与地图匹配算法融合入动态的贝叶斯网络中, 同样使用了多种传感器数据对车辆运动模型建模以达到对车辆的精确定位.

2 基本概念和问题描述

在本节中, 我们给出了本文中出现的一些关键变量和符号的说明, 如表1所示. 同时, 我们给出了本文相关概念的定义以及所要解决的问题的描述.

定义1. 轨迹点. 轨迹点p为GPS传感器的一条记录, 代表了电动自行车的时空状态, 可以由一个三元组表示, $ p: < lon, lat, t > $ 其中, $ lon $ 和 $ lat $ 分别表示电动自行车所处位置的经度和纬度, t表示这条记录的时间戳. 由于传感器存在测量误差, 因此轨迹点p的位置测量值并不完全等于真实位置.

定义2. 轨迹. 轨迹 $ Tr:{p_1} \to {p_2} \to \ldots \to {p_n} $ 是由一组按照时间顺序排列的轨迹点组成的序列. 图4中 $ {p_1} \to {p_2} \to \ldots \to {p_{14}} $ 展示了一段轨迹的示例.

定义3. 路网图. 路网图 $ G(V, E) $ 是一个有向图, 其中, (1) V是路网图中的顶点集合, 顶点 $ v \in V $ , 表示道路的端点; (2) E是路网图中有向边的集合, 有向边 $ e \in E $ , 在路网中表示道路段, 可表示为一个二元组 $ e: < e.start, e.end > , \forall e.start, e.end \in V $ . 在图4中, 顶点 ${v_1}, {v_2}, \ldots , {v_6}$ 为路网图上的道路端点, 顶点间的有向边为路网图中的道路段. 顶点与顶点间的有向边共同组成了路网图.

定义4. 路线. 路线 $ R:{e_1} \to {e_2} \to \ldots \to {e_n} $ 是由一组路网图中的相连道路段组成的序列, $ \forall {e_i} \in E, {e_{i + 1}}.start = {e_i}.end, i \in [1, n - 1] $ . 在图4中, $ {e_1} \to {e_2} \to {e_3} $ 为图中轨迹 $ {p_1} \to {p_2} \to \ldots \to {p_{14}} $ 实际经过的路线.

定义5. 停留点. 轨迹 $ Tr:{p_1} \to {p_2} \to \ldots \to {p_n} $ 上, 存在一条子轨迹 $T{r_{\rm sub}}:{p_a} \to {p_{a + 1}} \to \ldots \to {p_b}, 1 \leqslant a \leqslant b \leqslant n$ , 且 $ \forall a \leqslant i \leqslant j \leqslant b, dist({p_i}, {p_j}) \leqslant threshold $ , 其中 $ dist({p_i}, {p_j}) $ 为轨迹点 $ {p_i} $ 和 $ {p_j} $ 的距离, $ threshold $ 为停留点距离阈值, 则称 $ T{r_{\rm sub}} $ 上的轨迹点为停留点. 在图4中, 虚线框中的轨迹点 ${p_9} \to {p_{10}} \to \ldots \to {p_{14}}$ 为一组停留点的示例.

轨迹简化问题描述: 对于给定的原始轨迹 $ Tr:{p_1} \to {p_2} \to \ldots \to {p_n} $ , 使用轨迹简化算法将计算出一条新的轨迹 $T{r{'}}:p'_1 \to p'_2 \to \ldots \to p'_m$ , $ m < n $ , 是原始轨迹Tr的子序列, 即轨迹Tr'上的轨迹点也存在于原始轨迹Tr上, 称Tr'上的轨迹点为特征点. 轨迹简化比例 $ ratio = m/n, ratio \in (0, 1] $ 是指简化后的轨迹Tr'与原始轨迹Tr的轨迹点数量之比, $ ratio $ 越大, 轨迹简化过程保留的轨迹点越多.

地图匹配问题描述: 给定一条轨迹 $ Tr:{p_1} \to {p_2} \to \ldots \to {p_n} $ 和路网图 $ G(V, E) $ , 地图匹配算法将计算出轨迹Tr实际经过的路线 $ R:{e_1} \to {e_2} \to \ldots \to {e_m} $ . 在图4中, 路线 $ {e_1} \to {e_2} \to \ldots \to {e_m} $ 是轨迹 $ {p_1} \to {p_2} \to \ldots \to {p_{14}} $ 真实经过的路线.

3 基于卡尔曼滤波的轨迹简化与自适应地图匹配算法(KFTS-AMM)

在本节中, 我们详细地描述了本文提出的可自适应路网精度的电动自行车轨迹地图匹配方法KFTS-AMM. 其整体的架构图如图5所示, 主要由4个部分组成, 分别是轨迹提取、停留点识别与合并、轨迹简化与地图匹配, (1)轨迹提取: 删除GPS日志中存在属性缺失、冗余和存在明显错误的记录. 对GPS日志按照电动自行车设备ID与记录时间进行排序, 通过设置时间阈值分割提取出轨迹数据; (2)停留点处理: 使用时空聚类算法识别停留点, 结合旋转卡壳算法(rotating calipers algorithm)^[66]对停留点进行合并; (3)轨迹简化: 使用分段的卡尔曼滤波算法对轨迹进行修正, 然后使用基于DPhull算法^[16]改进的DPhull-ratio算法进行轨迹简化; (4)地图匹配: 使用分段的隐马尔可夫模型(HMM)地图匹配算法进行地图匹配.

本文提出的算法涵盖了电动自行车轨迹数据的清洗、轨迹提取、停留点处理、轨迹压缩简化以及地图匹配过程. 在实际应用场景中, 能够兼顾地图匹配的效率和准确率. 同时, 简化后的轨迹可以有效降低数据的存储空间.

3.1 轨迹提取 3.1.1 轨迹数据清洗

本文使用的是郑州市采集的真实电动自行车行驶数据. 数据传输过程中存在部分数据丢失和部分数据传输错误的情况, 因此需要对初始采集数据进行清洗. 原始日志字段包括电动自行车ID、数据采集时间、数据产生时间、经度与纬度. 其中, 电动自行车设备ID和数据产生时间可以唯一确定一条数据. 轨迹数据清洗的过程包括以下3个步骤: (1)删除存在属性缺失的记录; (2)对电动自行车设备ID和时间戳属性均相同的重复记录, 只保留其中一条; (3)由于已知郑州市的经纬度范围为东经112.42°至东经114.14°, 北纬34.16°至北纬34.58°, 因此将经度或纬度超出这一范围的异常记录删除.

3.1.2 轨迹数据的分割

交通部门使用安装在电动自行车上的边缘设备采集轨迹数据. 首先按照电动自行车设备ID对日志数据进行分割, 然后再按照数据产生时间进行排序, 通过这种方式可获取每辆电动自行车按照时间顺序排列的GPS日志序列. 为了从排序后的GPS日志中提取电动自行车的轨迹数据, 还需要按照时间间隔对其进行分割. 当电动自行车长时间静止后, GPS传感器会停止采集数据. 因此, 我们使用了固定的时间阈值来分割不同的轨迹, 我们将时间阈值设置为3 min. 若排序后的GPS日志中相邻记录的数据采集时间超过3 min时, 可以认为其属于不同轨迹.

3.2 停留点识别与合并

本节详细分析了电动自行车轨迹的停留点分布特点, 并据此特点提出了电动自行车轨迹上对停留点的处理算法StayPointsProcess.

停留点是指在某一段连续的时间内, 位置变化小于一定距离的一组轨迹点. 当车辆运动速度过于缓慢(例如用户推行电动自行车)、暂时停止行驶或者在某一较小区域内反复行驶时, 就会产生停留点. 在实际场景中, 用户会在骑行电动自行车的途中停车或下车推行并保持一定时间, 此时会随之产生停留点. 图6显示了一段真实轨迹与其GPS测量记录. 当电动自行车速度接近于零时, 出现了一组停留点, 由于GPS测量值在真实位置附近呈混合高斯分布, 因此停留点的GPS记录分布在图中椭圆形虚线框范围内. 图7展示了StayPointsProcess算法的过程, 图中左侧椭圆形虚线框内为一组停留点, 右侧多边形框为停留点的凸包, 凸包直径上的点为停留点合并得到的新轨迹点, 该点位于凸包的直径上. 在图7中, 使用均值计算得到的轨迹合并点, 受停留点密度影响, 未分布在这一区域的中心位置, 若以此点作为合并点, 则可能导致合并后的轨迹不平滑.

StayPointsProcess算法首先使用ST-DBSCAN算法^[57]完成停留点的识别, 然后通过旋转卡壳算法^[66]求出停留点的凸包并计算凸包直径, 最后使用电动自行车的平均速度在直径上取出若干点代替停留点, 从而形成新的轨迹. 若轨迹中停留点的个数为 $ p $ , 则构造凸包的时间复杂度为 ${\rm{O}}(p\log p)$ , 设凸包上点个数为 $ q $ , $ q \leqslant p $ 使用旋转卡壳算法计算此凸包直径的时间复杂度为 ${\rm{O}}(q)$ , 则停留点的合并过程时间复杂度为 ${\rm{O}}(p\log p)$ . 若轨迹上的轨迹点个数为 $n,\; n \geqslant p$ , 则ST-DBSCAN的平均时间复杂度为 ${\rm{O}}(n\log n)$ . 综上, StayPointsProcess算法的整体时间复杂度为 ${\rm{O}}(n\log n)$ .

算法1 StayPointsProcess描述了停留点的处理过程, 其输入包括原始轨迹Tr、空间邻域半径 $ ep{s_1} $ 、时间邻域半径 $ ep{s_2} $ 、停留点最小邻居数minPts以及距离计算指标 $ metric $ , 默认的距离计算指标为曼哈顿距离(Manhattan distance). 算法输出为新轨迹 $T{r{'}}$ . 算法 1第11行 $Diamete{r_{\rm stay}}$ 为使用旋转卡壳算法 $ RC $ 计算出的停留点凸包直径, 第12行 $merge(Diamete{r_{\rm stay}}, {v_{\rm avg}})$ 为停留点的合并过程, 其中 ${v_{\rm avg}}$ 为轨迹的平均速度.

算法1 StayPointsProcess的具体步骤如下.

• 第1步(第1–4行): 初始化当前停留点集合、输出轨迹为空. 构造聚类模型并拟合数据, 标记每个轨迹点是否为噪声点. 计算整个轨迹的全局平均速度.

• 第2步(第6行): 判断当前点是否为噪声点, 即该点是否不属于任何聚类, 若该点为噪声点, 则执行步骤3; 若不为噪声点, 则执行步骤4.

• 第3步(第7–14行): 获取当前聚类的标签, 向后遍历轨迹点, 将属于同一聚类的轨迹点加入当前停留点集合. 然后使用旋转卡壳算法计算停留点凸包直径, 在直径上取若干个点作为合并后的轨迹点, 加入输出轨迹. 清空停留点集合.

• 第4步(第16行): 将非噪声轨迹点加入输出轨迹.

算法1. 停留点处理算法StayPointsProcess.

输入: 原始轨迹 $Tr:{p_1} \to {p_2} \to \ldots \to {p_n}$ , 空间邻域半径 $ep{s_1}$ , 时间邻域半径 $ep{s_2}$ , 停留点最小邻居数 $minPts$ , 距离计算方法 $metric$ ;

输出: 停留点识别合并后的轨迹Tr'.

1. $ {s}_{\rm stay}\leftarrow \varnothing , T{r}'\leftarrow \varnothing ; $ // 停留点集合 ${s_{\rm stay}}$ 初始化为空集合

2. $model \leftarrow ST\text{-}DBSCAN\left( {ep{s_1}, ep{s_2}, minPts, metric \leftarrow \rm{`Manhattan\text{'}}} \right);$ // 构建聚类模型

3. $labels \leftarrow model.fit\left( {Tr} \right);$ // 使用聚类模型拟合数据

4. ${v_{\rm avg}} \leftarrow \sum\limits_{i = 1}^{n - 1} {dis({p_i},{p_{i + 1}})} /({p_n}.t - {p_1}.t)$ ; // 计算轨迹的全局平均速度

5. for $i = 1$ to $n$ do

6. 　 if $labels\left[ i \right] \ne - 1$ do // 判断该点是否属于某个簇类

7. 　　 $labe{l_{\rm cur}} \leftarrow labels\left[ i \right];$

8. 　　 while $labels\left[ i \right] = labe{l_{\rm cur}}$ do

9. 　　　 ${s_{\rm stay}}.add\left( {Tr\left[ i \right]} \right);$

10.　　　 $i \leftarrow i + 1;$

11.　　 end while

12.　　 $Diamete{r_{\rm stay}} \leftarrow RC\left( {{s_{\rm stay}}} \right);$ // 使用旋转卡壳算法计算停留点凸包的直径

13.　　 $Tr'.add\left( {merge\left( {Diamete{r_{\rm stay}}, {v_{\rm avg}}} \right)} \right);$ // 在凸包直径上取若干轨迹点, 加入新轨迹

14.　　 ${s_{\rm stay}}.clear()$ ; // 清空当前停留点集合

15.　 else

16.　　 $Tr'.add\left( {Tr\left[ i \right]} \right);$

17.　 end if

18. end for

3.3 轨迹简化

本节详细介绍了本文所提出的基于分段的卡尔曼滤波的轨迹简化算法KFTS. 该算法首先使用分段的卡尔曼滤波算法修正轨迹曲线, 然后使用基于DPhull算法改进的DPhull-ratio算法进行轨迹简化, 得到简化后的轨迹.

3.3.1 基于轨迹简化比例的轨迹简化算法DPhull-ratio

现有的离线与在线轨迹简化算法大多使用距离阈值作为轨迹简化过程的参数^[49], 只有当简化过程终止时, 才能计算出轨迹简化比例. 但是, 距离阈值的设置依赖于先验知识和对轨迹数据的初步统计结果, 若设置不够合理, 则无法得到期望的简化比例. 文献[51]使用递归的方式, 动态调整距离阈值以控制简化比例, 易造成轨迹上不同片段简化比例的不平衡. 本文在DPhull算法的基础上, 提出了基于轨迹简化比例的轨迹简化算法DPhull-ratio, 能够仅使用轨迹简化比例作为参数, 得到简化结果. 利用优先队列存储轨迹上的子轨迹段, 以达到不同片段简化比例的平衡.

DP算法的计算过程如图8所示, 设待简化的轨迹为 $Tr:{p_1} \to {p_2} \to \ldots \to {p_6}$ , 距离阈值为 ${D_{\rm threshold}}$ . 首先, 将轨迹起点 $ {p_1} $ 与终点 $ {p_6} $ 设置为特征点, 并做连线, 由于 $ {p_4} $ 与首尾连线垂直距离最远, 且超过了距离阈值 ${D_{\rm threshold}}$ , 因此 $ {p_4} $ 被作为轨迹的特征点. 然后, 在 $ {p_4} $ 处分割原轨迹 $ Tr $ , 得到子轨迹段 $ {p_1} \to {p_2} \to {p_3} \to {p_4} $ 与 $ {p_4} \to {p_5} \to {p_6} $ . 对这两段子轨迹分别递归执行上述过程, 最终确定特征点为 $ {p_1} $ 、 $ {p_2} $ 、 $ {p_4} $ 与 $ {p_6} $ , 简化后轨迹为 $T{r^{'}}:{p_1} \to {p_2} \to {p_4} \to {p_6}$ . 文献[51]在DP算法的基础上, 实时计算当前简化比例, 若达到输入简化比例, 则简化过程结束. 这种方式易造成轨迹不同片段简化比例的不平衡. DPhull算法引入了凸包的结构, 能够快速查找距离首尾连线距离最大的轨迹点, 使算法的最坏时间复杂度下降到了 ${\rm{O}}(n\log n)$ .

当需要获得某一设定简化比例的简化结果, 使用DP算法(包括DPhull)需要对距离阈值参数多次尝试, 导致计算资源的浪费. 为了解决此问题, 本文提出了DPhull-ratio算法, 一次运行即可得到期望的简化比例. DPhull-ratio算法基于贪心的思想, 维护了一个优先队列, 队列中每个节点记录了一段子轨迹和距离子轨迹首尾连线垂直距离最大的轨迹点, 称其为分割点, 节点的优先级即为分割点距离首尾连线的距离. 这样设置优先级的原因在于, 轨迹段上距离首尾连线垂直距离更大的轨迹点更能代表轨迹段的轮廓. DPhull-ratio算法每次取出一个队头节点, 在分割点处对子轨迹二分, 并计算分割后两个子轨迹的分割点、优先级用于构建新的队列节点, 最后将新节点加入优先队列.

DPhull算法在递归过程中, 每次选出距子轨迹段首尾连线距离最大, 且超过距离阈值的轨迹点作为特征点, 这一选择依据与DPhull-ratio算法的优先级设置方式相似. 因此, 对于同一条轨迹, 当DPhull-ratio算法与DPhull算法得到同一简化比例的结果时, 所选取的特征点集合应具有较大相似度, 与原始未简化轨迹的误差也应该较为接近.

在图9中, 待简化轨迹为 $Tr:{p_1} \to {p_2} \to \ldots \to {p_{10}}$ , 设轨迹简化比例为50%. 优先队列中初始只存在一个节点, 节点代表的轨迹段为 ${p_1} \to {p_{10}}$ , 分割点为 ${p_6}$ , 表示为 $ < {p_1} \to {p_{10}}, {p_6} > $ . 然后, 将其作为队列头节点取出, 令 ${p_6}$ 为特征点, 并在此处分割. 分割后得到两个新节点, 分别为 $ < {p_1} \to {p_6}, {p_4} > $ 与 $ < {p_6} \to {p_{10}}, {p_8} > $ , 将其加入队列. 迭代上述过程, 当选取的特征点个数达到总轨迹点个数的50%时, 简化过程结束. 算法2描述了上述计算过程.

DPhull-ratio算法增加了维护优先队列的消耗, 设优先队列长度为 $ n $ , 则优先队列的插入与删除操作时间复杂度为 ${\rm{O}}(\log n)$ . 若轨迹简化比例为 $ ratio $ , 则用于维护优先队列的时间复杂度为 ${\rm{O}}(ratio \times n\log (ratio\times n))$ . 每次取出节点后, 需要使用凸包计算分裂后节点的分割点, 时间复杂度为 ${\rm{O}}(\log n)$ , 因为共需取出 $ratio\times n$ 次队头结点, 总时间消耗为 ${\rm{O}}(ratio\times n\log n)$ . 综上, DPhull-ratio算法的时间复杂度为 ${\rm{O}}(ratio\times n\log (ratio\times n) + ratio\times n\log n)$ , 即 ${\rm{O}}(ratio\times n\log n)$ .

算法2. 基于压缩率的轨迹简化算法DPhull-ratio.

输入: 轨迹 $Tr:{p_1} \to {p_2} \to \ldots \to {p_n}$ , 轨迹简化比例 $ratio$ ;

输出: 简化后的轨迹 $ T{r{'}} $ .

1. $Reserve \leftarrow \{ 0\} ;$ // $Reserve$ 数组表示每个轨迹点是否为特征点, 数组长度为n

2. $ Reserve[0] \leftarrow 1, \;Reserve[n - 1] \leftarrow 1$ ; // 初始化首尾轨迹点为特征点

3. $cnt \leftarrow 2;$ // $cnt$ 为特征点个数

4. $Init\left( {Heap} \right)$ ; // $Heap$ 为优先队列, 对其进行初始化

5. while $cnt/n < ratio$ do

6. 　 $T{r_{\rm sub}}, {p_{\rm split}} \leftarrow Heap.top();$ // 选取队头轨迹段, 分割点

7.　 $Reserve[{p_{\rm split}}.index] \leftarrow 1, cnt \leftarrow cnt + 1$ ; // 分割点设置为特征点, 特征点数加1

8. 　 $Heap.pop();$ // 将队头元素弹出

9. 　 $T{r_{\rm sub\_1}}, T{r_{\rm sub\_2}} \leftarrow split\left( {T{r_{\rm sub}}, {p_{\rm split}}} \right)$ ; // 在分割点处二分轨迹段

10. 　 ${p_{\rm split\_1}} \leftarrow {{hull}}\left( {T{r_{\rm sub\_1}}} \right);$ // 使用凸包计算子轨迹段 $T{r_{{\rm{sub}}\_1}}$ 的分割点

11. 　 ${p_{\rm split\_2}} \leftarrow {{hull}}\left( {T{r_{\rm sub\_2}}} \right);$ // 计算子轨迹段 $T{r_{{\rm{sub}}\_2}}$ 的分割点

12. 　 $< T{r_{\rm sub\_1}}, {p_{\rm split\_1}} > .priority \leftarrow Dis\left( {{p_{\rm split\_1}}, \overline {T{r_{\rm sub\_1}}} } \right);$ // 节点的优先级为分割点与轨迹段首尾连线的垂直距离

13. 　 $< T{r_{\rm sub\_2}}, {p_{\rm split\_2}} > .priority \leftarrow Dis\left( {{p_{\rm split\_2}}, \overline {T{r_{\rm sub\_2}}} } \right);$ // 计算第2个节点的优先级

14.　 $Heap.push\left( { < T{r_{\rm sub\_1}}, {p_{\rm split\_1}} > } \right);$ // 将新节点加入优先队列

15. 　 $Heap.push\left( { < T{r_{\rm sub\_2}}, {p_{\rm split\_2}} > } \right);$

16. end while

DPhull-ratio算法的具体步骤如下.

• 第1步(第1–4行): 初始化Reserve数组, 将首尾节点设置为特征点, 特征点个数设置为2, 初始化优先队列.

• 第2步(第5行): 判断特征点个数与总轨迹点个数比例是否达到轨迹简化比例 $ ratio $ , 若达到, 则算法结束, 否则执行步骤3.

• 第3步(第6–9行): 取出队头节点, 将分割点设置为特征点, 并在分割点处分割轨迹.

• 第4步(第10–15行): 利用凸包分别计算步骤3得到的子轨迹段的分割点, 构造两个新节点, 并计算两个新节点的优先级, 然后加入优先队列.

3.3.2 卡尔曼滤波轨迹修正算法

卡尔曼滤波算法可以为每一时刻的状态进行最优估计, 其基于两个方程, 分别是系统状态的预测方程:

$ {X_k} = A{X_{k - 1}} + {q_k} $

(1)

和系统状态的观测方程:

$ {Z_k} = H{X_k} + {r_k} $

(2)

在公式(1)中, ${X_k}$ 为系统的状态向量, A为系统的状态转移矩阵(表示对系统状态转移的理想估计), ${q_k}$ 为系统的预测过程噪声, 服从均值为0, 协方差矩阵为Q的正态分布. 在公式(2)中, ${Z_k}$ 为系统的观测向量, H为系统的观测矩阵, ${r_k}$ 为观测噪声, 服从均值为0, 协方差矩阵为R的正态分布. 本文所使用的电动自行车轨迹只包含GPS数据, 因此卡尔曼滤波轨迹修正算法仅使用轨迹点的经纬度信息用于构建预测与观测方程.

卡尔曼滤波算法迭代地进行预测过程与更新过程. 预测过程使用2个公式:

$ \overline {{x_k}} = A{x_{k - 1}} $

(3)

$ \overline {{P_k}} = A{P_{k - 1}}{A^{\rm{T}}}+ Q $

(4)

公式(3)中, 先验状态向量 $\overline {{x_k}} $ 为上一时刻的后验状态向量 ${x_{k - 1}}$ 与状态转移矩阵A的乘积, 公式(4)计算先验卡尔曼误差估计协方差矩阵 $\overline {{P_k}} $ , ${P_{k - 1}}$ 为上一时刻的后验卡尔曼误差估计协方差估计矩阵. 更新过程使用3个公式, 分别为:

$ {K_k} = \overline {{P_{k - 1}}} {H^{\rm{T}}}{\left( {H\overline {{P_{k - 1}}} {H^{\rm{T}}} + R} \right)^{ - 1}} $

(5)

$ {x_k} = \overline {{x_k}} + {K_k}\left( {{Z_k} - H\overline {{x_k}} } \right) $

(6)

$ {P_k} = \left( {I - {K_k}H} \right)\overline {{P_k}} $

(7)

公式(5)计算得到的卡尔曼增益矩阵 ${K_k}$ 决定了最优估计更接近预测值还是更接近观测值, 公式(6)计算得到的后验状态向量 ${x_k}$ , 是经过卡尔曼滤波修正后的状态, 公式(7)计算后验卡尔曼误差估计协方差矩阵 ${P_k}$ .

算法3 KF描述了对输入轨迹段进行卡尔曼滤波轨迹修正的过程, 算法的输入包括轨迹段 $T{r_{\rm sub}}$ 、坐标预测误差方差 $va{r_{\rm predict}}$ 、坐标观测噪声方差 $va{r_{\rm gps}}$ , 输出为修正过的轨迹 $Tr'_{\rm sub}$ . 算法2中第9行、第10行、第11行和第13行分别对应公式(3)、公式(4)、公式(5), 算法第13行对应公式(6), 算法第14行对应公式(7).

算法3. 卡尔曼滤波轨迹修正算法KF.

输入: 轨迹 $T{r_{\rm sub}}:{p_1} \to {p_2} \to \ldots \to {p_n}$ , 预测误差方差 $va{r_{\rm predict}}$ , GPS观测误差方差 $va{r_{\rm gps}}$ ;

输出: 修正后轨迹 $Tr'_{\rm sub}$ .

1. $Tr'_{\rm sub} \leftarrow \varnothing ;$

2. $Init\left( {P, A, H, Q, R} \right);$ // 初始化 $ P、A、H、Q与R $

3. for $i = 1$ to $n$ do

4. 　 if $i = 1$ do

5. 　　 $X \leftarrow {\left( {\begin{array}{*{20}{c}} {T{r_{\rm sub}}\left[ 1 \right].x}&{T{r_{\rm sub}}\left[ 1 \right].y} & 0 &0 \end{array}} \right)^{\rm{T}}};$ // 初始状态向量

6. 　 else

7. 　　 $\vartriangle t \leftarrow T{r_{\rm sub}}\left[ i \right].t - T{r_{\rm sub}}\left[ {i - 1} \right].t;$

8. 　　 $A\left[ 0 \right]\left[ 2 \right] \leftarrow \vartriangle t, A\left[ 1 \right]\left[ 3 \right] \leftarrow \vartriangle t;$ // 更新状态转移矩阵

9. 　　 ${X_{\rm prior}} \leftarrow AX;$

10.　　 ${P_{\rm prior}} \leftarrow AP{A^{\rm{T}}} + Q;$ // 计算先验卡尔曼误差估计协方差矩阵

11.　　 $K \leftarrow {P_{\rm prior}}{H^{\rm{T}}}{\left( {H{P_{{\rm{prior}}}}{H^{\rm{T}}} + R} \right)^{ - 1}};$ // 计算卡尔曼增益

12.　　 $Z \leftarrow {\left( {\begin{array}{*{20}{c}} {T{r_{\rm sub}}\left[ i \right].x}&{T{r_{\rm sub}}\left[ i \right].y} \end{array}} \right)^{\rm{T}}};$ // 当前时刻轨迹点的GPS测量值

13.　　 $X \leftarrow {X_{\rm prior}} + K\left( {Z - H{X_{\rm prior}}} \right);$ // 计算后验状态向量

14.　　 $P \leftarrow \left( {I - KH} \right){P_{\rm prior}};$ //计算后验卡尔曼误差估计协方差矩阵

15.　　 $Tr'_{\rm sub}.add {\left( {X\left[ 0 \right], X\left[ 1 \right]} \right)} ;$ // 当前轨迹段中加入修正后的点

16.　 end if

17. end for

状态向量X为(x, y, v_x, v_y)^T, 其中x为电动自行车在该时刻的经度在GCS_WGS_1984 (world geodetic system 1984, GPS使用的坐标系统)投影坐标系下的横坐标, y为电动自行车在该时刻的纬度在GCS_WGS_1984投影坐标系下的纵坐标, x和y的单位为m, v_x与v_y分别为电动自行车在该时刻x方向与y方向上的速度分量, v_x和v_y单位为m/s. 由于状态向量有4个维度, 因此后验估计协方差矩阵P的维度为4×4, 初始值设置为零矩阵.

KF算法使用匀速直线运动作为电动自行车的运动模型, 因此状态转移矩阵A设置为4×4的单位矩阵. 观测值有2个, 分别是电动自行车的经度方向位置x和纬度方向位置y, 因此观测矩阵H设置为2×4, 初始值设置为 $\left( {\begin{array}{*{20}{c}} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \end{array}} \right)$ . 若x与y坐标预测误差为 $va{r_{\rm predict}}$ , 则系统预测过程噪声的协方差矩阵Q为 $\left( {\begin{array}{*{20}{c}} {va{r_{\rm predict}}} & 0 & 0 & 0 \\ 0 & {va{r_{\rm predict}}} & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{array}} \right)$ . 若 $ x $ 与 $ y $ 坐标测量误差为 $va{r_{\rm gps}}$ , 则系统观测过程噪声的协方差矩阵R为 $\left( {\begin{array}{*{20}{c}} {va{r_{\rm gps}}} & 0 \\ 0 & {va{r_{\rm gps}}} \end{array}} \right)$ .

3.3.3 电动自行车转弯判断

KFTS算法假设电动自行车匀速直线运动, 但是电动自行车在行驶过程中可能会在路口转弯, 这会影响卡尔曼滤波的状态预测方程的准确度. 因此KFTS算法将卡尔曼滤波算法分段进行, 即将转弯前后的轨迹段分别进行卡尔曼滤波轨迹修正.

存在轨迹段 ${p_a} \to {p_{a + 1}} \to \ldots \to {p_b}$ , 其平均速度定义为: $ \overrightarrow {{v_{a \to b}}} = \dfrac{{dist\left( {{p_a}, {p_b}} \right)}}{{{p_b}.t - {p_a}.t}} $ , 速度的方向由 ${p_a}$ 指向 ${p_b}$ , $dist\left( {{p_a}, {p_b}} \right)$ 为 ${p_a}$ 与 ${p_b}$ 的距离. 对于轨迹 $Tr:{p_1} \to {p_2} \to \ldots \to {p_n}$ , 设窗口大小为 $ w $ (表示窗口内轨迹点的数量). 若某窗口内轨迹段为 ${p_1} \to {p_2} \to \ldots \to {p_w}$ , 下一窗口内轨迹段为 ${p_{w + 1}} \to {p_{w + 2}} \to \ldots \to {p_{2\times w}}$ , 当两窗口平均速度的夹角大于45°, 且窗口内轨迹速度方向变化小于45°时, 可认为电动自行车在两窗口间完成了转弯. 判断是否转弯的条件定义如下: $angle\left( {\overrightarrow {{v_{1 \to w}}} , \overrightarrow {{v_{w + 1 \to 2\times w}}} } \right) > 45^\circ $ , 且 $angle\left( {\overrightarrow {{v_{i \to i + 1}}} , \overrightarrow {{v_{i + 1 \to i + 2}}} } \right) < 45^\circ $ , $i \in \left[ {1, w - 2} \right] \cup \left[ {w + 1, 2\times w - 2} \right]$ , 其中, $angle\left( {\overrightarrow {{v_1}} , \overrightarrow {{v_2}} } \right)$ 表示速度 $\overrightarrow {{v_1}} $ 与速度 $\overrightarrow {{v_2}} $ 间的夹角. 下文算法KFTS描述中的 $Turn\left( {{p_i}, {p_{i + 1}}} \right)$ 表示判断电动自行车是否在两窗口间, 即轨迹点 ${p_i}$ 与轨迹点 ${p_{i + 1}}$ 间完成了转弯.

3.3.4 基于分段卡尔曼滤波算法的轨迹简化算法

在图10中, $ v $ 和 $ {v{'}} $ 分别代表窗口1和窗口2内轨迹的平均速度. 由于 $ v $ 与 $ {v{'}} $ 的夹角接近90°, 且窗口1与窗口2内的轨迹方向变化很小, 因此可认为电动自行车在两个窗口间完成了转弯. 对转弯前后的轨迹分别进行卡尔曼滤波可以减小转弯对卡尔曼滤波的预测方程的影响.

算法4 KFTS描述了基于分段卡尔曼滤波的轨迹简化算法. KFTS算法的具体步骤如下.

• 第1步(第1行): 初始化输出轨迹 $T{r_{\rm new}}$ 为空, 当前轨迹片段 $T{r_{\rm sub}}$ 为空.

• 第2步(第3行): 将当前轨迹点加入轨迹片段 $T{r_{\rm sub}}$ .

• 第3步(第4行): 判断当前轨迹片段 $T{r_{\rm sub}}$ 轨迹点数是否超过窗口大小, 若超过, 则执行步骤4.

• 第4步(第5–10行): 判断轨迹是否完成转弯, 若完成, 则使用算法2 KF将当前轨迹片段 $T{r_{\rm sub}}$ 进行卡尔曼滤波轨迹修正. 修正后的轨迹片段加入输出轨迹 $T{r_{\rm new}}$ , 清空 $T{r_{\rm sub}}$ ; 若未完成转弯, 则将当前轨迹点加入 $T{r_{\rm sub}}$ .

• 第5步(第13行): 使用DPhull-ratio算法对修正后的轨迹进行简化.

KFTS算法首先使用卡尔曼滤波修正算法KF对轨迹进行修正, KF算法能够在线性时间内对一段轨迹滤波, 时间复杂度为 ${\rm{O}}(n)$ . 然后使用DPhull-ratio算法对轨迹简化, 由于DPhull-ratio的时间复杂度为 ${\rm{O}}(ratio\times n\log n)$ , 因此KFTS算法整体的时间复杂度为 ${\rm{O}}(ratio\times n\log n)$ .

算法4. 轨迹简化算法KFTS.

输入: 轨迹 $Tr:{p_1} \to {p_2} \to \ldots \to {p_n}$ , 窗口大小w, 预测误差方差 $va{r_{\rm predict}}$ , GPS观测误差方差 $va{r_{\rm gps}}$ , 轨迹简化比例 $ratio$ ;

输出: 简化后轨迹 $T{r_{\rm new}}$ .

1. $T{r_{\rm new}} \leftarrow \varnothing , T{r_{\rm sub}} \leftarrow \varnothing ;$

2. for $i = 1$ to $n$ do

3. 　 $T{r_{\rm sub}}.add\left( {Tr\left[ i \right]} \right);$

4. 　 if $T{r_{\rm sub}}.{si\textit{z}e}() \geqslant {{w}}$ do

5.　　 if $Turn\left( {Tr\left[ i \right], Tr\left[ {i + 1} \right]} \right)$ do // 判断轨迹是否在此处完成转弯

6. 　　　 $T{r_{\rm new}}.add\left( {KF\left( {T{r_{\rm sub}}, va{r_{\rm predict}}, va{r_{\rm gps}}} \right)} \right);$ // 将当前轨迹段进行卡尔曼滤波

7. 　　　 $T{r_{\rm sub}}.clear();$ // 清空当前轨迹段

8.　　 else

9.　　　 $T{r_{\rm sub}}.add\left( {Tr\left[ i \right]} \right);$

10.　　 end if

11. 　 end if

12. end for

13. $T{r_{\rm new}} \leftarrow$ $\mathrm{DPhull}\text{-}\mathrm{ratio}\left( {T{r_{\rm new}}, ratio} \right);$ // 使用DPhull-ratio算法进行轨迹简化

3.4 地图匹配

本节介绍了地图匹配部分的内容, 首先介绍了候选点的选择和候选图的构建方法, 其次介绍了如何解决无效轨迹片段导致的候选图中断问题, 最后给出了分段的隐马尔可夫模型(HMM)的地图匹配算法AMM, 其能够自适应路网精度, 解决因路网精度不足导致的匹配准确率下降的问题. 地图匹配算法的应用场景包括离线场景和在线场景^[4], 但是大部分电动自行车由于成本问题, 未安装计算设备, 因此本节中提出的地图匹配算法AMM只针对离线场景.

3.4.1 候选点选择和候选图构建

我们首先给出候选图的定义, 并且对候选点与候选边进行详细说明.

定义6. 候选图. 候选图 $ {G_c}(C, {E_c}) $ 是一个有向图, C为候选点集合, $ {E_c} $ 为候选点间相连边的集合. 候选图为一个分层的有向图, 有向边只存在于相邻层的节点之间, 且方向由上一层节点指向下一层节点. 图11为候选图的示例.

定义6中的候选点是对轨迹点实际位置的推测, 候选点所在的路段称为候选边, 通常取轨迹点在候选边上的投影点为候选点, 且轨迹点与候选点距离应小于固定的距离阈值, 该阈值与GPS误差范围有关. 候选点选择过程基于K近邻算法(K-nearest neighbor, KNN)实现, 可以表示为: $CS\left( p \right) = \mathrm{KNN}{\left( p \right)_{k, r}}$ , 其中 $ CS(p) $ 为轨迹点p的候选点集合, $ k $ 为近邻数, r为轨迹点和候选点的距离阈值, 对于轨迹 $Tr:{p_1} \to {p_2} \to \ldots \to {p_n}$ , 每一个轨迹点都存在若干候选点, 同一个轨迹点的候选点集合构成了候选图中的候选点层, 因此每一个轨迹点均在候选图中有一个候选点层与之对应. 候选点层 ${G_c}.layer\left[ i \right]$ : $\left( {c_i^1, c_i^2, \ldots , c_i^k} \right)$ , 表示轨迹点 $ {p_i} $ 的候选点集合. 候选点 $ c $ 有3个属性, 分别为所在候选边 $ c.e $ 、累积概率 $ c.accum $ 与前驱点 $ c.prev $ . 有向边 ${e_{a \to b}} \in {E_c}$ , 起点 $ a $ 与终点 $ b $ 分别为相邻候选层的候选点, $a \in {G_c}.layer\left[ i \right]$ , $b \in {G_c}.layer\left[ {i + 1} \right]$ . 在图11中, ${p_1} \to {p_2} \to {p_3}$ 为一段轨迹示例, 轨迹点与候选边的距离阈值为r, 道路段 $ {e_1} $ 与道路段 ${e_2}$ 为轨迹点 ${p_1}$ 的两条候选边, 轨迹点 ${p_1}$ 在 ${e_1}$ 与 ${e_2}$ 上的投影点 $c_1^1$ 与 $c_1^2$ 为其候选点.

文献[9]认为在候选图中, 连接候选点的边 ${e_{a \to b}}$ 的权值包括空间分析值与时间分析值, 由于本文使用的路网数据没有道路速度的相关信息, 因此本文只使用空间分析值代表 ${e_{a \to b}}$ 的权值. 空间分析值的定义与文献[9]中一致, 表示为: ${F_s}\left( {c_{i - 1}^t \to c_i^s} \right) = N\left( {c_i^s} \right)\times p\left( {c_{i - 1}^t \to c_i^s} \right)$ , 空间分析值为候选点 $c_{i - 1}^t$ 移动到候选点 $c_i^s$ 的概率值, 等于当前候选点的观察概率 $N\left( {c_i^j} \right)$ 与候选点 $c_{i - 1}^t$ 到候选点 $c_i^s$ 的转移概率 $p\left( {c_{i - 1}^t \to c_i^s} \right)$ 的乘积. 观察概率定义为: $N\left( {c_i^j} \right) = \dfrac{1}{{\sqrt {2\text{π} } \sigma }}{{\rm{e}}^{ - {{\left( {x_i^j} \right)}^2}/2{\sigma ^2}}}$ , 其中, 候选点 $c_i^j$ 的观察概率服从方差为 $\sigma $ 的正态分布, $x_i^j$ 为轨迹点 ${p_i}$ 和候选点 $c_i^j$ 的距离. 转移概率定义为: $p\left( {c_{i - 1}^t \to c_i^s} \right) = \dfrac{{{{\min}}\left( {{d_{i - 1 \to i}}, {w_{\left( {i - 1, t} \right) \to \left( {i, s} \right)}}} \right)}}{{\max ({d_{i - 1 \to i}}, {w_{\left( {i - 1, t} \right) \to \left( {i, s} \right)}})}}$ , 其中, ${d_{i - 1 \to i}}$ 表示轨迹点 ${p_{i - 1}}$ 与轨迹点 ${p_i}$ 的距离, ${w_{\left( {i - 1, t} \right) \to \left( {i, s} \right)}}$ 表示候选点 $c_{i - 1}^t$ 到候选点 $c_i^s$ 的最短路径距离. 由于车速限制, 最短路径距离不能过大, 否则连接两候选点的有向边不存在. 文献[14]论证了文献[9]中的转移概率定义可能会由于密集的轨迹数据与较大的邻域半径而出现概率大于1的问题, 因此本文使用文献[14]中的转移概率定义方式. 候选点的累积概率为上一层候选点的累积概率与其到该候选点的空间分析值之和的最大值, 计算方法为: $c_i^j.accum = \mathop {\max }\limits_{t \in \left[ {1, k} \right]} \left( {c_{i - 1}^t.accum + {F_s}\left( {c_{i - 1}^t \to c_i^j} \right)} \right)$ . 候选点的前驱点的计算方法为: $c_i^j.prev = \mathop {{{\arg\max}}}\limits_{{{c}} \in {G_c}.layer\left[ {i - 1} \right]} \left( {c.accum + {F_s}\left( {{{c}} \to c_i^j} \right)} \right)$ . 在图12中, 当前候选点 $c_2^1$ 的前驱点为上一层候选点中累积概率与其到 $c_2^1$ 的空间分析值之和的最大的候选点.

3.4.2 无效轨迹片段导致的候选图中断

无效轨迹片段是指, 在交通工具行驶的路线中, 某一段不存在于给定路网中的行驶轨迹. 一些使用全局匹配的地图匹配算法, 例如使用弗雷歇距离度量轨迹与路线相似度的算法^[21], 当存在无效轨迹片段的轨迹时, 仍会计算出一条路线, 但是这明显与实际情况不符. 使用基于HMM的地图匹配算法^{[14, 23–25, 37, 38]}时, 虽然会避免将无效轨迹片段匹配到路网上, 但是会引起候选图中断, 从而导致匹配过程中断, 最终部分有效轨迹片段也无法成功匹配. 在图13中, 轨迹片段 ${p_5} \to {p_8}$ 为无效轨迹片段, 该片段上的每个轨迹点的候选点都为空, 当使用基于HMM的地图匹配算法时, 在回溯寻找匹配路线时会导致算法中断.

3.4.3 分段的隐马尔可夫模型的地图匹配算法

无效轨迹片段可导致候选图中断, 并且会随着路网精度的降低而不断增多, 降低匹配准确率. 本文提出了分段的HMM的匹配算法AMM, 使用了分段匹配机制, 通过候选图区分有效与无效轨迹片段, 对连续候选图匹配路线.

算法5给出了自适应路网精度的地图匹配算法AMM的详细描述, 输入包括轨迹Tr, 路网图 $ G(V, E) $ , 轨迹点邻域半径r, 轨迹点近邻数k. 若输入轨迹 $ Tr $ 存在无效轨迹片段, 则输出为轨迹Tr实际经过的路线集合 ${S_{\rm path}}$ ; 若输入轨迹Tr不存在无效轨迹片段, 则 ${S_{\rm path}}$ 只包含一条路线. AMM算法的具体步骤如下.

• 第1步(第2行): 使用KNN算法计算当前轨迹点 $ Tr[i] $ 的候选点集合 $CS(Tr[i])$ .

• 第2步(第3行): 判断当前轨迹点的候选点集合是否为空, 若不为空, 执行第3步; 若为空, 执行第6步.

• 第3步(第4行): 判断当前候选图是否为空, 若不为空, 执行第4步; 若为空, 执行第5步.

• 第4步(第5–14行): 计算候选图上一层节点与候选图当前层节点间有向边的权值、当前层候选点的累积概率和前驱节点.

• 第5步(第16–19行): 当前层候选点为候选图的第一层候选点, 无前驱节点, 计算候选点的累积概率.

• 第6步(第22–24行): 为当前候选图匹配路线, 将所得路线加入路线集合 ${S_{\rm path}}$ , 清空当前候选图. 第23行中PathInference为路线推断算法, 详见算法6.

算法5. 地图匹配算法AMM.

输入: 轨迹 $Tr:{p_1} \to {p_2} \to \ldots \to {p_n}$ , 路网图 $G\left( {V, E} \right)$ , 轨迹点邻域半径r, 轨迹点近邻数k;

输出: 路线集合 ${S_{\rm path}}$ .

1. for $i = 1$ to $n$ do

2. 　 $CS\left( {Tr\left[ i \right]} \right) \leftarrow KNN{\left( {Tr\left[ i \right]} \right)_{k, r}};$ // KNN算法寻找轨迹点的候选点

3. 　 if $CS\left( {Tr\left[ i \right]} \right) \ne \varnothing$ do

4. 　　 if ${G_c} \ne \varnothing$ do

5. 　　　 ${C_{\rm last}} \leftarrow {G_c}.layer\left[ { - 1} \right];$ // 取出候选图最后一个候选点层

6. 　　　 foreach $c \in CS\left( {Tr\left[ i \right]} \right)$ do

7. 　　　　 ${G_c}.C.add\left( c \right);$ // 候选点加入候选图

8. 　　　　 foreach $v \in {C_{\rm last}}$ do

9. 　　　　　 ${e_{v \to c}} \leftarrow {F_s}\left( {v \to c} \right);$ // 计算空间分析值、有向边权值

10.　　　　　 ${G_c}.{E_c}.add\left( {{e_{v \to c}}} \right);$

11.　　　　 end foreach

12.　　　　 $c.accum \leftarrow \mathop {\max }\limits_{v \in {C_{\rm last}}} \left( {v.accum + {F_s}\left( {v \to c} \right)} \right);$ // 计算累积概率

13.　　　　 $c.prev \leftarrow \mathop {\arg \max }\limits_{v \in {C_{\rm last}}} \left( {v.accum + {F_s}\left( {v \to c} \right)} \right);$ // 寻找候选点的前驱点

14.　　　 end foreach

15.　　 else

16.　　　 foreach $c \in CS\left( {Tr\left[ i \right]} \right)$ do

17.　　　 ${G_c}.C.add\left( c \right);$ // 加入最后一个候选点层

18.　　　 $c.accum \leftarrow N\left( c \right);$ // 第1层候选点的累积概率等于观察概率

19.　　　 end foreach

20.　　 end if

21.　 else

22.　　 $R \leftarrow \mathit{PathInference}\left( {{G_c}} \right)$ ; // 路线推断

23.　　 ${S_{\rm path}}.add\left( R \right);$ // 将路线加入输出路线集合

24.　　 ${G_c}.clear();$ // 清空当前候选图

25.　 end if

26. end for

算法5第6步中使用的PathInference算法为路线推断算法, PathInference算法的输入为候选图 $ {G_c} $ , 输出为该段轨迹对应的路线R.

算法6. 路线推断算法PathInference.

输入: 候选图 ${G_c}$ ;

输出: 路线R.

1. ${c_{\rm cur}} \leftarrow \mathop {\arg \max }\limits_{c \in {G_c}.layer\left[ { - 1} \right]} \left( {c.accmu} \right);$ // 选出候选图 ${G_c}$ 最后一层 ${G_c}.layer\left[ { - 1} \right]$ 中累积概率最大的候选点 ${c_{\rm cur}}$

2. for $i = {G_c}.layers.si{\textit{z}}e()$ to 1 do // 从候选图的最后一层回溯到第1层

3. 　 $R'.add\left( {{c_{\rm cur}}.e} \right);$ // 将候选点所在路段e加入逆序路线 $R'$

4. 　 ${c_{\rm cur}} \leftarrow {c_{\rm cur}}.prev;$ // 当前节点置为当前候选点的前驱

5. end for

6. $R \leftarrow R'.reverse();$ // 将逆序路线R'逆置并输出路线R

图14为路线推断示例, 图中箭头由候选点指向其前驱节点. 首先从轨迹点 $ {p_4} $ 的候选点中选出累积概率最大的候选点 $c_4^3$ , 然后不断回溯寻找前驱点, 最后可得到路线 $c_1^4.e \to c_2^2.e \to c_3^2.e \to c_4^3.e$ , 即为当前候选图对应的实际路线.

3.4.4 地图匹配算法AMM复杂度分析

设输入轨迹的轨迹点个数为 $ n $ , 每个轨迹点候选点个数为 $ k $ . AMM算法的时间复杂度由以下几个部分组成: (1)为每个候选点查找候选点以构建候选图, 时间复杂度为 ${\rm{O}}(kn)$ ; (2)计算候选图中相邻层候选点间的转移概率, 本文使用了文献[14]中预先计算路网图顶点最短路径并查表的方式避免重复计算, 候选图层数等于轨迹点个数 $ n $ , 每层有 $ k $ 个候选点, 因此时间复杂度为 ${\rm{O}}({k^2}n)$ ; (3)路线推断算法PathInference迭代寻找候选点的前驱点, 时间复杂度为 ${\rm{O}}(n)$ . 综上, AMM算法的整体时间复杂度为 ${\rm{O}}(nk + {k^2}n + n)$ , 即 ${\rm{O}}({k^2}n)$ .

4 实　验

本节介绍了实验所使用的路网数据和轨迹数据, 给出了实验中地图匹配准确率与地图匹配速度的定义. 同时, 本节设计了8组试验, 其中, 实验1用于测试基于简化比例的轨迹简化算法DPhull-ratio的误差; 实验2用于验证轨迹简化算法KFTS能否有效修正轨迹; 实验3用于测试不同超参数对地图匹配算法AMM的匹配准确率的影响; 实验4使用多个现有地图匹配算法与本文提出的地图匹配算法AMM进行匹配准确率比较; 实验5用于测试本文提出的轨迹简化算法KFTS对地图匹配准确率的影响; 实验6用于验证经过KFTS算法简化后的轨迹数据能否提高地图匹配速度; 实验7用于测试停留点处理算法StayPointsProcess对地图匹配准确率的影响; 实验8对3个真实电动自行车轨迹案例进行了分析. 实验运行在Windows 10操作系统计算机上, 计算机处理器为Intel Core i5-4590 CPU, 主频为3.30 GHz, 内存(RAM)为14 GB.

4.1 数据介绍 4.1.1 路网数据

实验使用了郑州市的开放街道地图(OpenStreetMap, OSM. 来源: https://www.openstreetmap.org)的开源路网数据集, 其统计信息如表2所示, 包含了OSM采集到的所有级别的道路, 包含了机动车道路、非机动车道路以及人行小路. OSM为各类道路设置了不同的标签, 表3展示了各类标签与国内道路类型的大致对应情况以及道路数量占比. 路网图拓扑结构如图15所示.

4.1.2 轨迹数据

原始GPS日志共有2687214282条记录, 删除存在缺失值、冗余或经纬度有明显错误的记录后, 剩余2515378612条记录, 经提取得到10651448条轨迹, 涉及1349641辆电动自行车, 覆盖的时间范围为2019年9月16日至2019年9月30日, 传感器设置采样时间间隔为10 s.

提取的轨迹的长度分布如图16(a)所示, 约75%的轨迹长度分布在1–4 km, 这说明使用大部分用户倾向于在中距离路程中使用电动自行车. 轨迹上相邻点的平均距离频率分布如图16(b)所示, 近半数的轨迹点间距分布在50–100 m, 由此可以看出约半数用户骑行电动自行车的速度在18–36 km/h之间. 本文选取提出轨迹中的1000条轨迹作为轨迹子集, 对其进行人工标记, 完成后续实验.

4.2 地图匹配准确率以及速度评价标准

本文实验使用的匹配准确率定义如公式(8)所示:

$ \mathop \sum \limits_i^{\left| {GT} \right|} compare\left( {GT\left[ i \right], MR\left[ i \right]} \right)/\left| {GT} \right| $

(8)

其中, $GT$ 表示被标记过的轨迹, $\left| {GT} \right|$ 表示被标记过的轨迹个数, MR表示轨迹匹配路线, $compare\left( {GT\left[ i \right], MR\left[ i \right]} \right)$ 用于比较轨迹正确结果 $GT\left[ i \right]$ 与匹配结果 $MR\left[ i \right]$ , 若完全一致, 则结果为1. 公式(8)衡量了成功匹配的轨迹个数占轨迹总数的比例. 匹配速度定义为匹配总时间除于轨迹总个数: ${S_{\rm match}} = {T_{\rm match}}/\left| {GT} \right|$ .

4.3 实验1: 基于简化比例的轨迹简化算法DPhull-ratio误差测试

实验1用于验证本文提出的DPhull-ratio算法能否在得到期望简化比例的前提下, 将简化后轨迹与原始轨迹的误差控制在一定范围内. 我们使用了最大垂直欧式距离 $PE{D_{\rm max}}$ (max perpendicular Euclidean distance)以衡量轨迹简化的误差. 我们将本文提出的DPhull-ratio算法与DPhull算法进行了对比, 比较了在轨迹集合简化比例相同时两种算法的 $PE{D_{\rm max}}$ , 同时我们使用了杰卡德相似系数(Jaccard similarity coefficient)衡量了两种算法在简化过程中选取的特征点集合的相似度.

设原始轨迹为Tr: ${p_1} \to {p_2} \to \ldots \to {p_n}$ , 简化后的轨迹为 ${{Tr':p}}_{\text{1}}{'} \to {{p}}_{\text{2}}{'} \to \ldots \to {{p}}_{{m}}{'}$ , 数组Reserve: $\{ 0, 1|i \in \left[ {1, n} \right]\} $ 表示原始轨迹 $ Tr $ 上的轨迹点是否存在于简化后轨迹 $T{r{'}}$ 上, Reserve数组元素取值为0或1, Reserve $\left[ i \right] = 1$ 表示轨迹点 ${p_i}$ 在轨迹简化过程中被保留, 否则被舍弃. $PE{D_{\max}}$ 为被简化的轨迹点距离被简化后的轨迹的最大距离, 表示了轨迹简化的误差上界. 最大垂直欧式距离 $PE{D_{\max}}$ 的计算方法如公式(9)所示:

$ PE{D_{\max}} = \mathop {\max }\limits_{i \in \left[ {1, n} \right]} dist\left( {{p_i}, \overline {{p_{i\_prev}} \to {p_{i\_next}}} } \right)\times\left( {1 - {{Reserve}}\left[ i \right]} \right) $

(9)

其中, ${p_{i\_prev}}$ 表示轨迹点 ${p_1}$ 至轨迹点 ${p_{i - 1}}$ 中距离轨迹点 ${p_i}$ 最近的被保留的轨迹点, ${p_{i\_next}}$ 表示轨迹点 ${p_{i + 1}}$ 至轨迹点 ${p_n}$ 中距离轨迹点 ${p_i}$ 最近的被保留的轨迹点, $\overline {{p_{i\_prev}} \to {p_{i\_next}}} $ 为这两个轨迹点相连得到的线段, $dist( {{p_i},} {\overline {{p_{i\_prev}} \to {p_{i\_next}}} } )$ 表示轨迹点 ${p_i}$ 距离该线段的垂直距离. 设两种算法计算得到的轨迹保留数组分别为Reserve₁和Reserve₂, Jaccard相似系数的计算方法为: $\left| {{{Reserve}_1} \cap {{Reserve}_2}} \right|/\left| {{{Reserve}_1} \cup {{Reserve}_2}} \right|$ .

我们对DPhull-ratio算法与DPhull算法在轨迹集合简化比例相同时的 $PE{D_{\max}}$ 指标进行了测试. 具体做法是: 首先, 使用某一距离阈值作为DPhull算法的参数, 此时轨迹集合中每条轨迹的简化比例是不同的; 然后, 对于轨迹集合中的每条轨迹, 使用DPhull-ratio算法以相同的简化比例对其进行简化, 最终轨迹集合的平均简化比例也相同. 轨迹集合的平均垂直欧式距离 $PE{D_{\max}}$ 随轨迹集合平均简化比例的变化如图17所示. 从结果可以看出: (1)两种算法得到的 $PE{D_{\max}}$ 指标都会随着平均简化比例的上升而下降, 这是因为简化比例越大, 简化过程保留的轨迹点越多, 因此简化后轨迹会更接近原轨迹的形状; (2)简化比例较大时, DPhull-ratio算法与DPhull算法得到的 $PE{D_{\max}}$ 指标十分接近, 而当简化比例较小, 为10%与20%时, DPhull-ratio算法的效果更好. 这说明了DPhull-ratio算法通过优先队列来选择特征点的机制是有效的, 同时使用队列节点的分割点与节点包含轨迹段的垂直距离作为节点优先级能够帮助选出更能描绘轨迹轮廓的轨迹点.

DPhull-ratio算法与DPhull算法选取的特征点集合的Jaccard相似系数随简化比例的变化如图18所示, 从结果可以看出: (1)两种算法选取的特征点集合的Jaccard相似系数较高, 简化比例为10%时相似系数也有近80%; (2) Jaccard相似系数随简化比例增大而不断提高, 这与图17中两种算法得到的 $PE{D_{\max}}$ 的相似度变化趋势一致; (3)简化比例较低时, Jaccard相似系数较低, 联系图17可以说明, 在简化比例较低时, 根据合理设置的优先级选取特征点的机制更有效, 因此较低简化比例下Jaccard相似系数也较低.

4.4 实验2: 轨迹简化算法KFTS效果评价

实验2用于验证本文提出的轨迹简化算法KFTS是否可以有效地修正轨迹点和平滑轨迹, 分别测试了KFTS算法与DPhull-ratio算法在轨迹简化比例 $ ratio $ 为0.2、0.5与0.8的情况下, 简化后轨迹的平均误差 $Erro{r_{\rm mean}}$ 与最大误差 $Erro{r_{\max}}$ . 轨迹简化比例 $ ratio $ 表示简化后轨迹点个数与原轨迹点个数的比例. $ ratio $ 越大, 表示保留的轨迹点越多. 平均误差 $Erro{r_{\rm mean}}$ 定义为: $Erro{r_{\rm mean}} = \displaystyle\mathop \sum \limits_{i = 1}^n dist\left( {pr{j_i}, {p_i}} \right)/n$ , 其中 $pr{j_i}$ 为轨迹点 ${p_i}$ 在真实路线上的投影点, $dist\left( {pr{j_i}, {p_i}} \right)$ 表示投影点 $pr{j_i}$ 与轨迹点 ${p_i}$ 的欧式距离. 最大误差 $Erro{r_{\max}}$ 定义为: $Erro{r_{\max }} = \mathop {\max }\limits_{i \in \left[ {1, n} \right]} dist\left( {pr{j_i}, {p_i}} \right)$ .

实验结果如表4所示, 从中可以看出: (1)在使用不同 $ ratio $ 的情况下, KFTS算法均能得到更低的最大误差 $Erro{r_{\max}}$ 与平均误差 $Erro{r_{\rm mean}}$ ; (2)随着 $ ratio $ 不断提高, KFTS算法相比于DPhull-ratio算法的评估指标的提升比例不断下降, 这是因为基于DP算法^[15]的轨迹简化算法倾向于保留GPS误差较大的轨迹点. 当 $ ratio $ 较小时, 大部分GPS误差较小的轨迹点因简化而被舍弃, 所以KFTS算法相较于DPhull-ratio在评估指标上提升较为明显; 当 $ ratio $ 较大时, GPS误差较大的轨迹点所占比重减小, 因此KFTS算法相比DPhull-ratio算法的评估指标的提升比例降低; (3)当ratio相同时, 最大误差相较于平均误差的提升比例更高, 这是因为KFTS算法可以较好地修正GPS误差大的轨迹点. 同时, 评估指标的提升比例均未超过10%, 原因可能是KFST算法中的分段卡尔曼滤波模型不够精细, 仅使用轨迹点的时间与位置属性不能得到更精确的估计值.

总之, 实验2结果证明, 本文提出的轨迹简化算法KFTS能够在简化轨迹的同时, 较好地修正轨迹, 有效解决了基于DP算法^[15]的轨迹简化算法得到的简化后的轨迹误差较大的问题.

4.5 实验3: 地图匹配算法AMM超参数选择

本文提出的地图匹配算法使用的超参数有2个, 分别为轨迹点邻域半径r和候选点个数k. 本文设计了实验3用于测试不同超参数对于地图匹配算法AMM匹配准确率的影响. 实验3结果如图19所示.

从实验结果中可以观察到, 随着邻域半径r的扩大与候选点个数k的增加, 匹配准确率不断提高. 当邻域半径r选择50 m、候选点个数k选择4时, 匹配准确率基本达到最高. 继续扩大邻域半径r与增加候选点个数k所带来匹配准确率的提高很小, 反而会降低算法的匹配效率. 文献[5]指出, 民用GPS传感器的测量误差在3 m之内, 且在城市路网中, 单方向中可供电动自行车行驶的平行路段最多不会超过2条, 因此邻域半径r选择50 m、候选点个数选择4时, 不会遗漏候选点. 为了达到最佳实验效果, 后续进行的实验所使用的邻域半径r选择50 m、候选点个数选择4.

总之, 实验3结果证明, AMM算法的两个超参数轨迹点邻域半径r与候选点个数k, 分别在选择4与50 m时, 可以兼顾准确率与效率, 达到最优效果.

4.6 实验4: 地图匹配算法匹配准确率测试

OSM对不同级别的道路设置了相应的标签. 在表5中, 我们给出了几种常见道路的标签, 道路的标签代表了道路的级别. 表5中的道路的级别从motorway到residential依次降低. 一般情况下道路的精细度较低时, 包含的道路级别也相对较高, 例如只保留城市路网中的主干道路以组成路网的骨架结构. 本文通过选择不同道路标签来模拟不同精度的路网, 其中: (1)高精度路网包含motorway、trunk、primary、secondary、tertiary、unclassified、service与residential等8种道路标签的道路; (2)中精度路网包含motorway、trunk、primary、secondary、tertiary、unclassified与service等7种道路标签的道路, 路网的密集程度相对于高精度路网来说有所降低; (3)低精度路网包含motorway、trunk、primary、secondary、tertiary与unclassified这6种道路标签的道路, 路网的密集程度最低.

图20展示了郑州市东经113.7756°至东经113.8254°, 北纬34.7741°至北纬34.8134°范围内的不同精度的路网图, 道路密度由高精度路网图到低精度路网图不断降低.

实验4选取5种已有地图匹配算法作为对比算法, 分别是FMM算法^[14]、ST-Matching算法^[9]、Segmented-MM^[29]、IVMM算法^[24]以及AIVMM算法^[25]. 其中FMM算法、ST-Matching算法、IVMM算法与AIVMM算法均为基于隐马尔可夫模型的地图匹配算法, Segmented-MM算法属于局部/增量算法. IVMM算法与AIVMM算法在HMM模型的基础上加入了互动投票机制, 利用了候选点间的相关性. AIVMM算法在IVMM的基础上加入了两个特定情况的约束, 以减少匹配路线的迂回. FMM算法为了解决匹配轨迹存在反向移动的问题, 设计了惩罚机制, 使用惩罚因子( $ pf $ )来表示惩罚的程度, $ pf $ 参数取值范围为0–100, 取值越高表示对反向移动的修正更强. 本文选取了 $ pf $ 分别等于0、50与100这3种情况下的FMM算法进行对比, 分别表示为FMM-1、FMM-2与FMM-3. 输入轨迹数据均经过了停留点处理.

实验4结果如图21所示, 从中可以看出: (1)随着路网精度的下降, 本文提出的AMM算法与其他方法的匹配准确率都在下降, 但是其他方法的匹配准确率下降更为明显, 说明了本文提出的AMM算法对不同精度路网的适应度更强. 当路网精度下降时, 电动自行车轨迹数据中无法匹配到给定路网的无效轨迹增多, 造成基于HMM模型的算法匹配过程中断. Segmented-MM算法同样无法对无效轨迹进行匹配. 因此对比方法的匹配准确率下降较快; (2) AMM算法匹配准确率也会随路网精度降低而下降, 原因在于当路网精度不断降低、无效轨迹不断增多时, 候选图中连续层数过少, 导致隐马尔可夫模型不够精确, 因此准确率降低; (3)由于输入轨迹已进行了停留点处理, 因此Segmented-MM针对高频率采样下, 路口处轨迹匹配进行的优化对准确率的提升不明显; (4) FMM-1、FMM-2与FMM-3在相同精度的路网下, 匹配准确率相差不大, 原因在于本文所使用的电动自行车轨迹数据采样频率较高, 相邻轨迹点距离较近, 而文献[14]提到的反向移动问题多出现在相邻轨迹点距离较远的情况下.

总之, 实验4结果表明, 本文提出的地图匹配算法AMM相较于其他对比方法在不同精度路网上的匹配准确率更高, 同时随着路网精度的下降, 匹配准确率的下降幅度更小, 表现出更强的适应度, 有效解决了因路网精度不足而产生的无效轨迹片段使地图匹配算法准确率降低的问题.

4.7 实验5: 轨迹简化算法KFTS对地图匹配算法准确率的影响

本文提出的轨迹简化算法KFTS算法包含了卡尔曼滤波轨迹修正与轨迹简化两个步骤. 实验5用于验证KFTS算法中的两个关键步骤对地图匹配算法准确率的影响. 首先, 为了验证卡尔曼滤波轨迹修正步骤对于地图匹配准确率的影响, 我们在实验中将KFTS-AMM算法与其变体DPhull-ratio-AMM进行了对比, KFTS-AMM算法在轨迹简化前加入了卡尔曼滤波轨迹修正步骤, 而DPhull-ratio-AMM则没有轨迹修正步骤. 其次, 为了验证轨迹简化步骤对于地图匹配准确率的影响, 我们评估了KFTS-AMM与DPhull-ratio-AMM在不同轨迹简化比例下的地图匹配准确率, 其中轨迹简化比例 $ ratio $ 设置为0.2、0.5、0.8与1. 轨迹简化比例 $ ratio $ 表示简化后轨迹点个数与原轨迹点个数的比例, $ ratio $ 越大, 则保留的轨迹点越多. 当 $ ratio $ 为1时, KFTS-AMM与DPhull-ratio-AMM均不对轨迹进行简化. 路网数据使用高精度路网图.

实验5结果如后文图22所示, 从中可以看出: (1)在轨迹简化比例 $ ratio $ 相同时, KFTS-AMM比DPhull-ratio-AMM的匹配准确率较高, 说明经过卡尔曼滤波修正后的轨迹可以提高地图匹配的准确率; (2)随着 $ ratio $ 不断下降, DPhull-ratio-AMM的匹配准确率不断下降, 这是因为随着 $ ratio $ 的下降, GPS误差大的轨迹点所占比例不断增大, 因此对匹配准确率的影响越来越明显; (3)当 $ ratio $ 由1下降到0.8、0.5时, KFTS-AMM的匹配准确率变化很小, 可以忽略不计, 说明对采样频率较高的轨迹进行适度简化, 不会降低地图匹配准确率; (4)当 $ ratio $ 下降到0.2时, KFTS-AMM的匹配准确率下降较大, 原因是当 $ ratio $ 下降到0.2时, GPS误差较大的轨迹点所占比例较大, 极大地影响了匹配准确率.

综上, 实验5结果证明本文提出的轨迹简化算法KFTS中的卡尔曼滤波轨迹修正步骤能够在一定程度上提高地图匹配的准确率, 同时在对轨迹进行适度的简化, 能够在不降低准确率的前提下提升地图匹配效率.

4.8 实验6: KFTS算法对地图匹配效率的影响

实验6用于测试经轨迹简化算法KFTS处理后的轨迹能否给地图匹配算法带来效率上的提升. 实验6使用KFTS算法在不同轨迹简化比例 $ ratio $ 下得到轨迹数据作为地图匹配算法AMM的输入, 分别测试匹配速度. $ ratio $ 分别选取0.2、0.4、0.6、0.8与1, $ ratio $ 越大, 则保留的轨迹点越多. 实验6结果如图23所示, 平均每条轨迹的匹配时间随着 $ ratio $ 不断下降而缩短, 说明了经过轨迹简化算法KFTS处理后的轨迹可以加快地图匹配的速度.

4.9 实验7: 停留点处理算法对匹配准确率的影响

实验7用于评估本文提出的停留点处理算法StayPointsProcess对地图匹配算法准确率的影响. 实验7分别测试了在不同精度路网上, AMM算法与未使用StayPointsProcess算法的AMM算法的地图匹配准确率. 实验7结果如后文图24所示, 在高、中、低3种精度的路网中, 相比于未使用停留点处理的AMM算法, 包含了停留处理过程的AMM算法均能取得更高的地图匹配准确率. 实验7结果表明, 本文提出的停留点处理算法StayPointsProcess可以有效识别与合并停留点, 从而解决了电动自行车轨迹数据中存在大量停留点所导致的地图匹配准确率下降的问题^[4], 能够有效提升电动自称车轨迹地图匹配算法的性能.

4.10 实验8: 案例分析 4.10.1 案例1: 存在停留点轨迹的匹配

图25中白色网格为路网图, 红色点为电动自行车轨迹. 图25(a)为一段包含停留点的电动自行车轨迹, 椭圆虚线框内为一组停留点. 当未对停留点处理时, 使用AMM算法的匹配结果如图25(c)所示. 由于AMM算法设置了候选点间最短路径的距离上限, 虽然避免了不必要的迂回^[4], 但是匹配过程则因停留点的候选点间不存在最短路径而中断. 使用StayPointsProcess算法对停留点识别与合并后, AMM算法匹配结果如图25(d)所示, 与正确结果一致.

4.10.2 案例2: 存在GPS误差较大轨迹点的轨迹的匹配

图26中白色网格为路网图, 红色点为电动自行车轨迹. 图26(a)中圆形虚线圆框内的轨迹点GPS误差较大. 图26(c)与图26(d)中的虚折线分别为使用DPhull-ratio算法和KFTS算法简化后的轨迹, 红色折线分别为匹配后的路线. 图26(c)中, 由于误差较大的轨迹点得到了保留, 导致匹配路线错误. KFTS-AMM对GPS误差较大的轨迹点进行了修正, 因此得到了正确匹配结果.

4.10.3 案例3: 存在无效轨迹片段的轨迹的匹配

图27中白色网格为路网图, 红色点为电动自行车轨迹. 图27(a)中椭圆虚线框内为无效轨迹片段. 图27(c)中为FMM算法^[14]的匹配结果, 由于无效轨迹片段存在, 导致匹配过程中断, 因此只有部分轨迹得到了正确的匹配结果. KFTS-AMM算法的匹配结果如图27(d)所示, 与正确结果一致.

5 总结与展望

本文提出了一种可自适应路网精度的电动自行车轨迹地图匹配方法KFTS-AMM. 该方法融合了基于卡尔曼滤波算法的轨迹简化算法(KFTS), 以及分段的隐马尔可夫模型的地图匹配算法(AMM). 本文提出的算法解决了现有地图匹配算法因电动自行车轨迹数据存在停留点、采样频率高、轨迹密度大且受限于路网精度较低而存在无效轨迹片段等特点导致的匹配准确率下降和匹配效率不高的问题. 基于真实的郑州市电动自行车轨迹数据的实验结果证明, 本文提出的算法与现有方法相比, 在不同精度路网图上的适应度更强, 匹配准确率更高, 且简化后的轨迹数据能够显著提升匹配速度. 未来的研究工作包括以下几个部分: (1) 解决地图匹配在隧道、高架桥与高楼附近等多路径效应明显的区域内地图匹配准确率较低的问题; (2) 尝试利用更多传感器数据进行数据融合, 使卡尔曼滤波模型更精确, 进一步提高轨迹修正效果; (3) 本文所提方法不适用于海量轨迹数据处理场景, 在未来的工作中尝试结合大数据与并行化技术, 提升算法的计算效率; (4) 地图匹配过程中获取的无效轨迹片段能够用于更新路网数据, 完整的轨迹对于深入挖掘用户出行规律十分重要, 我们将会在未来的工作中对无效轨迹片段深入地分析与挖掘.