受益于GPS(global position system)定位技术的发展及普及，大量轨迹数据在源源不断地产生.这些轨迹数据中，车辆轨迹是非常重要的一部分.一方面，车辆沿着道路行驶，但由于GPS定位存在误差，GPS的读数有时候会偏离车辆真实的行驶路线.另一方面，尽管车辆在空间中是连续运动的，但GPS通常每隔一段时间产生一个读数，也就是说，由GPS点组成的轨迹是车辆真实位置的采样.当车辆行驶速度比较快而GPS采样频率较低时，2个相邻GPS点会相距较远，两点之间的位置信息无法直接获得，这称为GPS轨迹的不确定性.地图匹配将GPS点映射到路网上，能够找到车辆真实的行驶路线，减少轨迹的不确定性.如图1所示，给定一条轨迹tr=
p1,p2,p3,p4,p5
(称为GPS点轨迹)和路网G，地图匹配将轨迹tr转化成路网G中一个首尾相连的路段序列e1,e2,e3,e4,e5(称为匹配轨迹).
Fig. 1 Illustrationof Map-Matching
图1 轨迹地图匹配示例
作为轨迹数据挖掘最基本的操作之一，地图匹配在许多空间数据智能场景中都非常有用.例如，在轨迹压缩[1-2]应用中，先将GPS点轨迹转化成由路段ID组成的匹配轨迹，然后利用整数编码或者最短路径的方法对轨迹进行压缩，减少存储空间.文献[3]通过比较匹配轨迹和原始GPS点距离评估城市区域的GPS环境友好性，即哪些地方的GPS误差高，哪些位置的GPS误差低，从而指导部署辅助定位设备.文献[4]通过地图匹配操作实现对移动机器人的定位.文献[5]通过查询经过某一点的匹配轨迹，找到从该点出发一段时间内能够抵达的区域范围.在轨迹异常检测[6-7]中，通过比较匹配轨迹与预设路径的差异判断车辆是否偏离了给定路线.文献[8-9]查询经过给定路径的所有匹配后的轨迹，可以分析每条道路的交通流量[10-11]，或者更好地为车辆导航[12].
轨迹地图匹配非常困难，特别是在并行道路或者高架桥等路网密集的区域.目前应用最广泛的方法是基于隐马尔可夫模型(hidden Markov model, HMM)的地图匹配方法[13]，因为其不仅考虑了路网的连接性，还考虑了GPS点之间的相互关系，因此具有较高的准确率.然而，基于HMM的地图匹配算法的运行效率很低，难以满足实时性很高的场景需求，在大规模轨迹数据的场景下，单台计算机的资源有限，这个问题显得尤为突出.例如，在实时可达区域分析文献[5]中，需要实时匹配全城范围内的车辆，才能捕捉每条道路的交通情况，最终响应用户实时可达区域查询的需求；在危化品监管项目中，需要实时监控危化品车辆是否偏离了既定的路线，因为一旦发现异常，需要实时地对其进行拦截，防止其驶入居民区域，造成安全隐患.即使对于离线匹配的场景，也需要尽可能地加快地图匹配的效率，减少资源的消耗.
有大量的工作用来加速基于HMM的地图匹配效率.例如，通过构建空间索引[14-16]加速查找每个GPS点的候选路段，但仅采用这种方法对整体提速并不明显，因为基于HMM的地图匹配的瓶颈在于大量的最短路径计算.为此，一部分研究工作简化了最短路径计算，包括减少待匹配的GPS点个数[2,17-18]、略过较低概率匹配点之间的最短路径计算[19-20]、用近似值替代最短路径[21]等.然而，简化计算的方法通常需要一些额外的信息，包括GPS点方向、路网等级等，其通用性受限，还会降低地图匹配的准确率.因此，另一部分工作[22-25]直接加速最短路径的计算，但它们仍有较大的提升空间.还有一部分工作[18,26-29]利用并行计算或分布式计算技术，充分利用更多机器资源来提高计算效率，但是它们没有考虑到路网和轨迹之间的大小差异，按照采样轨迹对数据进行分区，容易造成数据不均衡，或者轨迹切分过细破坏轨迹完整性从而导致匹配准确率降低等问题.
针对已有工作的不足，面向大规模轨迹数据的场景，本文提出了一个基于层次收缩[30]的分布式地图匹配处理框架(contraction hierarchy-based map-matching, CHMM).首先，在大多数实际项目中，涉及到的路网数据非常小，因此，我们将路网数据广播到整个集群，对轨迹进行随机分区，保证了数据平衡和轨迹完整性.其次，针对基于HMM的地图匹配算法效率较低问题，我们提出了使用基于层次收缩的多对多最短路径查询算法，在保证匹配准确率不变的情况下，大大加速了图匹配过程.层次收缩算法最初被用来加速计算两点之间的最短路径，其思想是，先对路网进行预处理，增加一些捷径，并对每个道路顶点进行编号；然后采用双向搜索算法，利用编号限制搜索方向，利用捷径跳过一些节点的遍历，从而显著提升最短路径的计算效率.本文对基于层次收缩的最短路径算法加以改造，实现多对多(即多个出发点，多个目标点)的最短路径查询，从而减少重复计算，提升地图匹配效率.注意，本文提出的优化方案与其他加速方案[2,19]是正交的，即可以与其他加速方案同时使用.本文的贡献主要包含4个方面：
1) 提出了一个分布式地图匹配框架.采用了一种新的分区策略，即将路网广播、轨迹分段后再随机分区，从而实现负载均衡，保持轨迹完整性.
2) 提出了基于层次收缩的多对多最短路径搜索算法，从而加速地图匹配算法的效率;
3) 采用了4个真实的数据集进行实验，实验结果表明，本文提出的分布式算法能够大大提升地图匹配的效率，同时具有很强的可扩展性；
4) 本文提出的CHMM算法已经落回到了产品JUST[31-33]中，并支持了多个真实的项目落地.我们开源了CHMM算法的核心代码，并提供了一个在线演示环境供读者体验(1)http://chmm.urban-computing.com/.
1 先验知识
本节先给出一些形式化的定义，然后介绍Apache Spark的必要知识.
1.1 问题定义
定义1. 路网(road network).路网是一个有向图G=(V,E)，其中，V={v1,v2,…,vn}是节点集合，表示交叉路口；E={e1,e2,…,em}是有向边集，表示路段集合.每条路段e∈E是一个三元组e=(vstart,vend,l)，vstart,vend∈V表示e的通行方向为从vstart到vend，路段长度为l，称vstart和vend分别是e的起始节点和终止节点.为方便起见，令e=(vstart,vend)，l=len(vstart,vend).对于双向路段，可以用2条有向边来表示.
定义2. 轨迹(trajectory).轨迹tr=p1,p2,…,pn是由同一个移动对象产生的按照时间先后顺序组织的GPS点序列，其中GPS点pi，1≤i≤n，是一个三元组pi=(lati,lngi,timei)，表示该移动对象在timei时刻上报的位置坐标为(lati,lngi).注意，由于GPS定位存在误差，移动对象上报的位置有时并不是其真实所处的位置.
定义3. 地图匹配(map-matching).给定轨迹tr=p1,p2,…,pn和路网G=V,E，地图匹配将轨迹tr映射到路网G上，得到一条匹配后的轨迹tr′=(e1,time1),(e2,time2),…,(em,timem)，其中，(ei,timei)表示移动对象进入路段ei的时间为timei.对于任意(ei,timei)和(ei+1,timei+1)，1≤i<m，均有timei<timei+1，且ei的终止节点是ei+1的起始节点，即ei.vend=ei+1.vstart.这保证了匹配后的轨迹tr′在路网上是连续的.为简便起见，对于匹配后的轨迹，我们忽略其时间信息，即tr′=e1,e2,…,em.在接下来的章节中，我们统称原始轨迹tr和匹配后的轨迹tr′为轨迹.
1.2 Apache Spark
Apache Spark[34](简称Spark)是一个基于内存的分布式处理框架，它能够处理大规模的数据，并保证了容错性.Spark提供了一种数据抽象，称为RDD(resilient distributed dataset)，每个RDD包含多个分区(partition)，不同的分区可能分布在集群的不同机器中，每台机器可以有多个executor，每个executor可以并行处理多个分区，每个分区即为一个并行处理单元.RDD可以通过一些并行的操作(例如map,filter,reduce等)进行构建和转换.我们应当尽可能保证每个RDD分区内的数据量差不多，从而实现负载均衡，加快整体的效率.RDD可以缓存在内存中，或者持久化在磁盘中，以便重复使用，减少数据计算，这在循环迭代的场景中非常有用.在Spark中，可以将变量广播(broadcast)到所有executor中，或者通过zipPartitions操作将2个RDD合并成一个RDD，这样能够实现2个RDD之间的信息交换.Shuffle操作能够将数据在不同的分区中重新组织和分配，但是shuffle操作代价非常昂贵，因为它会造成不同分区甚至不同机器之间的数据转移，因此我们应该尽可能地减少shuffle的次数.
尽管CHMM算法是基于Spark实现的，但是我们可以很方便地将它移植到其他的分布式计算框架中，例如Apache Hadoop[35]或者Apache Flink[36].
2 基本框架
图2给出了基于层次收缩的分布式地图匹配算法CHMM的基本框架，共包含3个部分：轨迹数据预处理、路网数据预处理、轨迹地图匹配.
1) 轨迹数据预处理.如图2左上部分所示，我们首先读取轨迹数据，然后对轨迹进行分段，即将一条长的轨迹分段成多条短的子轨迹，一方面能够更好地组织轨迹数据，另一方面也能够加快后续地图匹配的效率(详见第3.1节).
2) 路网数据预处理.如图2左下部分所示，我们首先读入路网数据，然后利用层次收缩算法[30]进行收缩，最终得到收缩后的路网，用来加快后续地图匹配过程中的最短路径计算(详见第3.2节).
3) 轨迹地图匹配.如图2右半部分所示，这个过程接收层次收缩的路网和分段后的轨迹作为输入，基于这2个数据集，首先对轨迹数据进行分区，使得每个分区中的数据量尽可能相等，从而达到负载均衡.然后在每个分区里面并行执行地图匹配算法，得到最终的匹配轨迹(详见第3.3节).我们还将CHMM算法落回到了JUST中[31-33,37]，并提供了一个用户界面[38]，用户可以通过该系统触发地图匹配过程，并查看匹配结果(详见第4节).
Fig. 2 Framework of CHMM
图2 CHMM算法的基本框架
3 分布式地图匹配算法CHMM
在本节中我们将详细介绍CHMM的各个步骤，最后对算法性能进行分析.
3.1 轨迹数据预处理
轨迹数据是按照移动对象组织的，一个移动对象产生的GPS点可以视作一条轨迹.然而，GPS点是源源不断地产生的，若仅仅按照移动对象组织GPS点，那么一条轨迹就是一个面向未来无限延展的流，难以对轨迹进行存储、管理和挖掘.特别地，针对基于HMM的地图匹配情况，上述组织轨迹的方式会带来3个问题.1)基于HMM的地图匹配的时间复杂度为
其中n为一条轨迹中GPS点数目，|V|是路网顶点数目，|E|是路网的边数，k是每个GPS点的候选路段条数(见3.3节).轨迹越长，涉及到的最短路径搜索空间越大，导致地图匹配过程非常耗时.2)若移动对象长期停在一个地方(例如停车场)，但是未关闭GPS定位，这段时间产生的GPS点对地图匹配结果正确性影响不大，反而会降低地图匹配效率，应该将之去除.3)若2个相邻GPS点产生的时间间隔过长(比如一个小时)，那么它们之间地图匹配结果的正确性难以保证.
为解决上述问题，本文在地图匹配之前，先对轨迹进行分段.轨迹分段是指将一条长的轨迹分解成多条短的子轨迹.有多种轨迹分段的方法，例如：根据出租车的一次载客过程[40]或者共享自行车一次借车过程[41]对轨迹进行分段，根据快递员一次送货过程[26,42]对轨迹进行分段等.上述的分段方法都需要具体的业务知识，并不适用于所有轨迹数据.因此，本文采用基于驻留点和时间跨度的方法对轨迹进行分段，只利用轨迹的时空信息，适用范围更广(任何其他的轨迹分段方法都可用于CHMM算法).具体分成2步：驻留点检测和轨迹切分.
1) 驻留点检测.驻留点是移动对象在一段时间t内产生在一个区域范围d中的GPS点.我们采用文献[42]的驻留点检测方法，依次将每个GPS点作为锚点，检查其后继节点是否满足时间和空间条件.如图3所示，假设p1是当前的锚点，其后继节点p2,p3,p4到p1的距离都小于d，p5到p1的距离大于d.接着我们检查p4与p1的时间间隔是否大于t.若大于t，那么p1,p2,p3,p4都是驻留点.接着，我们以p2作为锚点重复前面的过程，直到所有的GPS点都检查完毕.在本文中，我们令d=100 m，t=100 s.该取值不但能够排除十字路口等红绿灯的情况(通常红绿灯的等待时间在100 s以内)，也能够避免GPS的定位误差而带来的影响(GPS的定位误差通常不超过100 m)，还能够较好地捕捉长时间停留在某处的情形.
Fig. 3 Illustration of trajectory segmentation
图3 轨迹分段示例
2) 轨迹切分.获得所有的驻留点后，我们以驻留点作为其中一类分割点.另一类分割点考虑相邻2个GPS点的时间间隔，若间隔大于t′，那么这2个GPS点之间也会形成一个分割点，如图3的p6,p7所示.轨迹分段利用这些分割点，将长轨迹分割成多段短的子轨迹.图3中将一条长轨迹分成了3段子轨迹.注意，t′不能设置得过低，因为如果设置得过低，将会产生过多的子轨迹，每条子轨迹的起/止节点无法充分利用前/后GPS点的信息，影响匹配结果的准确性.在本文中，我们令t′=60 min，因为如果2个连续GPS点超过了这个时间，它们之间的距离通常会相距非常远，它们相互之间参考价值不大，其间的地图匹配结果的正确性也难以保证.
Fig. 4 Distributed trajectory processing based on Spark
图4 基于Spark的分布式轨迹预处理
分布式实现.利用Spark，可以很容易实现轨迹数据预处理.如图4所示，我们首先利用textFile操作从外部存储中读取轨迹数据，形成轨迹RDD.然后调用Spark的map操作识别出轨迹的驻留点，得到新的轨迹RDD，最后继续调用flatMap操作对轨迹切分(因为一条长轨迹将会切分成多条子轨迹)，得到分段后的子轨迹RDD.注意，整个轨迹预处理过程不涉及到数据的shuffle，利用Spark的并行处理能力，能够高效地预处理海量的轨迹数据.
3.2 路网数据预处理
为了加快基于HMM的地图匹配算法效率，一种最有效的方式是提升最短路径的查询效率.基于层次收缩的最短路径查询[30]权衡了预处理时间、存储空间、查询效率3方面的因素，能够以较少的预处理时间、较少的存储空间、高效地支持两点之间的最短路径.为此，在本节，我们对路网数据进行层次收缩，得到收缩后的路网；同时讨论当路网非常大时，我们应当如何处理.在3.3节，我们阐述如何利用层次收缩后的路网加速地图匹配.
3.2.1 路网层次收缩
假设给定路网G=V,E,其节点V={v1,v2,…,vn}按照某种“度量标准”排序.层次收缩就是依照这种“度量标准”，依次对节点进行收缩.具体过程为：依次将节点v及其相邻边从原来的路网中删除，同时增加一些边，使得剩下的图节点之间的最短路径长度与删除v之前各节点间的最短路径长度保持一致，直到所有的节点都删除完毕.这里，删除的次序称为该节点的层数(level)，添加的边称之为捷径(shortcut).文献[30]指出，从直观上看，应该尽可能地让收缩后的路网的边数更少.因此，他们发现边差(edge difference)是一个非常有效的排序度量标准.节点v的边差是指，删除节点v后增加的捷径数与删除的边数之差.注意，给定一个度量标准，找到最优解非常耗时，因此可以使用启发式方法找到次优解，但对结果的正确性没有影响.
接下来用一个实际例子来阐述层次收缩过程.给定如图5(a)所示的路网，图中每条边都是双向边，边上的数字表示路段长度，节点旁边的数字表示该节点的边差.例如，节点v2的边差是1-3=-2，因为假设删除v2，我们会删除v2的3条相邻边(v1,v2),(v2,v3),(v2,v5)，但是为保持所有节点对的最短路径长度不变，我们需要增加1条长度为6的捷径(v1,v3)，否则节点v1到v3的最短路径长度无法得以保留.我们随机选择最小边差的一个节点v2进行收缩，并对v2进行层数编号，得到如图5(b)所示的路网，虚线表示已经删除的边和节点，边(v1,v3)表示增加的捷径，节点内的数字表示层数编号.同时，我们更新剩余节点的边差.接着，我们收缩节点v1，得到如图5(c)所示的路网.以上过程不断迭代重复，直到所有的节点都收缩完毕，最后得到如图5(d)所示的收缩路网.
Fig. 5 Example of contraction hierarchy
图5 层次收缩的例子
还有其他的排序度量标准，在文献[30]中有更详细的讨论.算法1给出了路网层次收缩算法的伪代码.注意，整个过程是在原来的路网图上进行修改.我们用一个优先队列q来保证编号的顺序，并初始化当前层次数(行①).接着，依次检查q中每个节点，直到队列q为空.在每次检查时，对q中的最小度量节点v出队，并对v进行编号(行③～④).然后，检查删除v是否影响剩余图的最短路径(我们只需要检查v的两两邻居之间的路径即可).如果v是其2个邻居u,w的最短路径上的节点，共分成2种情况：1)若(u,w)早已存在，更新原有边(u,w)的长度(行⑩)；2)若(u,w)不存在，我们添加捷径(u,w)到图G中(行
).在处理完节点v后，我们需要更新v的所有邻居的排序度量(行
)，保证下一次迭代中，检查的是度量数最小的节点.
算法1. 路网层次收缩算法.
输入：原始路网G=V,E,某排序度量标准;
输出：收缩后的路网G.
① 初始化一个优先队列q,并将所有的v∈V入队，按照给定的度量标准排序;level=1;
② while q≠null do
③ v←q.pop();
④ v的层次编号设为level;level++;
⑤ for each(u,v)∈E并且u尚未编号do
⑥ for each(v,w)∈E并且w尚未编号do
⑦ if u→v→w是u到w的最短路径then
⑧ l←len(u,v)+len(v,w);
⑨ if(u,w)已经存在E中then
⑩ 更新(u,w)的长度为l;
else
将(u,w,l)加入到E;
end if
end if
end for
end for
更新v的所有邻居的排序度量;
end while
return G.
3.2.2 超大路网处理方案
一般而言，城市中的路网数据更新周期较长，例如按月或者按季度更新，因此，上述路网收缩过程可以周期性执行，并将收缩结果存储起来，再次使用时直接加载收缩后的路网即可.此外，在实际落地项目中，通常是以城市甚至区域作为交付的基本单元，而一个城市中的路网通常不会很大，例如，作为中国大都市的北京，其路网占用的磁盘大小约为100 MB，任何一台标准的服务器均可以进行处理.对于无法完全加载到单台机器内的情形，我们可以先对路网切割，例如按照行政区域进行切割或者采用图分割算法[43]，将一张大的路网切割成多个小的子路网，然后对每个子路网单独进行层次收缩，最后对落在每个子路网内的轨迹数据，依次采用本文提出的分布式地图匹配算法(见第3.3节).
Fig. 6 Handling trajectories and road networks at boundaries
图6 处理边界的轨迹和路网
然而，如果对路网进行了切割，我们需要慎重考虑子路网边界处的路段和轨迹数据.如图6所示，分割线(注意实际情况中，分割线未必是直线)将路网划分成了2个子路网G1和G2，轨迹tr=p1,p2,p3,p4,p5横跨了这2个子路网.由于GPS点的定位误差，靠近图分割线的GPS点的真实位置可能落在了所在子路网之外，例如，尽管p3落在了图G1中，但是其候选路段可能在图G2中.此外，GPS点p3的匹配路段也受p4的匹配结果的影响.为了提高匹配精度，我们将每个子路网向外扩张α米，得到扩张路网，然后再对扩张路网执行层次压缩操作.对于边界处的轨迹，我们移除落在扩张路网之外的GPS点，只保留原始路网内GPS点的匹配结果.文献[44]指出，当α=500 m时，对匹配精度几乎没有影响.
例如图6中，R1+R2表示子路网G1的区域，R3+R4表示子路网G2的区域，R1+R2+R3表示G1的扩张路网
的区域，R2+R3+R4表示G2的扩张路网
区域.当处理内的轨迹数据时，我们去除p5，只对子轨迹p1,p2,p3,p4进行地图匹配，但只保留p1,p2,p3的匹配结果.同样，当处理内的轨迹数据时，我们去除p1和p2，只对子轨迹p3,p4,p5进行地图匹配，但只保留p4和p5的匹配结果.最后将边界处的p3和p4的匹配点，通过最短路径拼接成轨迹tr的最终匹配结果.
Fig. 7 Distributed road network preprocessing based on Spark
图7 基于Spark的分布式路网预处理
3.2.3 分布式实现
图7展示了利用Spark实现路网预处理的过程.首先，通过textFile操作从外部文件中读取路网数据，形成路网RDD.注意这里的路网图是以路段的形式存储的，在不同的分区之间以点分割的形式进行分割.接着，我们通过partitionBy操作对路网重新分区.我们实现了一种按照城市范围的自定义分区方法，使得落在同一个城市扩张范围内的路段重新分配到同一个分区中.若路段与多个城市扩张范围相交，则会分配到多个分区中.在每个新的分区中，我们可以得到一个子路网，通过mapPartitions操作实现对子路网的层次收缩，采用mapPartitions的原因是单个分区内的所有路段都会参与整体的收缩运算.最后，同样通过调用mapPartitions操作，将每个分区内的层次收缩子路网存储为单个文件，同时记录了每个文件的路网空间范围，方便后续地图匹配过程使用.
3.3 轨迹地图匹配
本文采用基于HMM的地图匹配算法[13]，因为其具有较高的匹配精度，而且只使用轨迹的时空信息，无需其他的业务知识，适用性广，是目前使用最多的地图匹配算法之一.给定轨迹tr=p1,p2,…,pn和子路网G，基于HMM的地图匹配算法共包含3步：寻找候选点、构建转移图和确定匹配路径.
1) 寻找候选点.我们首先为每个子路网构建空间索引(例如R-trees[14]或者Quad-trees[15])，方便快速找到轨迹中每个GPS点pi周围β范围内(β<α)的候选路段.注意，这里不包括压缩路网中的捷径.令β=100 m，因为GPS点的定位误差通常不会超过100 m.然后，对于每个候选路段，找到pi的候选点.假设
是GPS点pi的第j条候选路段，那么pi在上的候选点是中距离pi最近的位置点，即：
(1)
其中，d(pi,q)表示两点pi,q之间的欧氏距离.pi在候选路段
的候选点，要么是pi到的垂直投影点，要么是的端点.如图8(a)所示，pi的候选路段为
对应的候选点为
Fig. 8 Steps of Map-Matching based on HMM
图8 基于HMM的地图匹配步骤
2) 构建转移图.在获得每个GPS点的候选点后，我们构建一张转移图，如图8(b)所示.转移图是一张有向加权层次图，按轨迹GPS点的时间先后顺序进行层次组织，相邻GPS点的候选点之间进行全连接.图中每个节点和边赋予一个概率值.针对GPS点pi的某个候选点
，其概率符合标准正态分布：
(2)
其中，σ是GPS定位误差的标准差，我们采用文献[39]的取值，令σ=20 m，d(
,pi)是两点,pi之间的欧氏距离.对于每条边
其概率为
(3)
其中，d(pi,pi+1)表示2个GPS点的欧氏距离，
表示2个候选点的最短路径距离.分子取较小值、分母取较大值是为了防止概率大于1.
3) 确定匹配路径.这一步使用Viterbi算法[45]找到转移图中概率最大的路径.形式化而言，就是确定每个GPS点的最终候选点
使得式(4)最大：
(4)
确定每个GPS点的候选点后，通过计算相邻2个候选点的最短路径，即可获得最终匹配结果.
上述基于HMM的地图匹配算法，最耗时的地方是第2步中计算相邻2个候选点的最短路径距离
当路网比较密集时，每个GPS点的候选点很多，情况更加糟糕.我们做了相关实验，发现即使使用A*算法[46]计算两点间的最短路径，其耗时也占整个地图匹配耗时的95%以上.基于此，本文采用层次收缩算法加快最短路径的计算，从而加速地图匹配的过程.
算法2. 基于层次收缩的最短路径搜索算法.
输入：层次收缩路网G=(V,E),起始点s，终止点t;
输出：从s到t的最短路径长度dmin.
① if s和t在同一条路段上then
② return s到t的路段距离;
③ end if
④ 初始化优先队列qs和qt，用来存储待检查的节点，其键为对应节点到s或t的最短路径距离;
⑤ 初始化哈希表hs和ht，其键为已检查的节点，值为对应节点到s或t的最短路径距离;
⑥ 将s的所有出口节点加入qs，将t的所有入口节点加入qt；dmin←∞;
⑦ while qs≠null或者qt≠null do
⑧ if qs≠null do /*正向走一步*/
⑨ (di,vi)←qs.pop();
⑩ if di>dmin do /*正向扩张结束*/
清空qs;
else if hs[vi]=null do
hs[vi]=di;
if ht[vi]≠null do
/*表示在中间相遇*/
dmin←min{dmin,di+ht[vi]};
end if
for each(vi,vk)∈E且vk.level>
vi.level do
将(di+len(vi,vk),vk)加入qs;
end for
end if
end if
if qt≠null do /*反向走一步*/
(dj,vj)←qt.pop();
if dj>dmin do /*反向扩张结束*/
清空qt;
else if ht[vj]=null do
ht[vj]=dj;
if hs[vj]≠null do
/*表示在中间相遇*/
dmin←min{dmin,dj+hs[vj]};
end if
for each(vk,vj)∈E且vk.level>
vj.level do
将(dj+len(vk,vj),vk)加入qt;
end for
end if
end if
end while
return dmin.
在第3.2节，我们已经对路网进行了层次压缩.给定压缩后的路网G，起始点s，终止点t，基于层次收缩的最短路径搜索过程[30]类似于双向Dijkstra搜索算法.如算法2所示，若s和t(可能不是路网节点)在同一个路段，则直接计算s到t的路段距离并返回(行①～③).否则，我们通过3个步骤获得最终结果：
1) 初始化(行④～⑥).我们借助于2个优先队列来分别控制正向和反向的搜索顺序，借助于2个哈希表来分别记录某个节点是否已经被访问过以及它们到s或t的最短距离，并将s的出口节点(即从s出发能够抵达的第1个道路节点)和t的入口节点(即能够抵达t的最后一个道路节点)分别加入到2个优先队列中.dmin表示当前解，初始化为正无穷大.
2) 循环迭代(行⑦～).不断从起止两端依次单步扩张，直到2个队列均为空为止.共包含2个部分：正向扩张(行⑧～)和反向扩张(行～).正向扩张中，首先从优先队列qs中弹出一个离s最近的节点vi，然后验证它离s的距离是否超过了当前解dmin.若已经超过了dmin，那么qs中所有节点离s的距离均大于dmin，因此我们清空qs以结束正向扩张.否则，验证vi是否在正向扩张中被访问过了，只有vi是首次被访问我们才从vi继续扩张，以防止重复计算.从vi继续扩张时，我们首先将vi加入hs，用以记录vi是否被访问过以及它到s的最短距离.然后检查其是否被反向扩张遍历到了(称为相遇).若相遇，则表明找到了一条经由vi的候选路径，更新dmin.接着，验证vi的所有后继邻居节点vk，只有vk的层级大于vi的层级，才将vk加入到优先队列qs.反向扩张的逻辑类似，特别注意它是沿着路段的反方向进行扩张，且边的起始节点vk的层级要大于终止节点vj的层级.
3) 结果返回(行).返回dmin作为最终结果.
例如，要找到图5(d)中v1到v6的最短路径长度，其过程如表1所示，每一步执行完毕后，各变量的取值在相应位置给出.例如，在第4轮扩张时，正向扩张从优先队列qs中取出v5，它到s的距离为6，小于dmin，因此不会手动清空qs.此外，它尚未被正向扩张访问过，因此将v5加入hs中.由于它已被反向扩张访问过，说明此时已经找到了一条经由v5的路径，更新dmin.然后检查v5的邻居，发现不存在层级比v5大并且未访问过的节点.此时，2个优先队列均为空，搜索完毕.最后得到v1到v6的最短距离为7.注意，整个过程中，我们都没有访问节点v2.
相对于双向Dijkstra算法，基于层次收缩的最短路径算法主要有2点优化：1)利用捷径，能够绕过一些中间节点，从而减少迭代次数；2)限制了搜索方向，即正向搜索沿着层次增加的方向，反向搜索沿着层次减少的方向.因此，基于层次收缩的最短路径搜索算法更快.注意，基于层次收缩的最短路径算法虽然与A*算法[46]有相似之处，即通过启发式方法控制了搜索方向，但两者的启发式方法不能直接结合在一起，因为对于层次收缩算法而言，每条边的层级数都不一样，也就无法先按照层级数扩张，再按照A*算法的启发函数来扩张；在路网收缩前，我们也无法事先得知起止点，也就无法利用A*的启发函数来对路网进行收缩.
Table 1 Example of Shortest Path Search Based on CH
表1 基于层次收缩的最短路径搜索示例
步骤qshsqthtdmin初始化(0,v1)空(0,v6)空∞第1轮迭代(4,v4),(6,v3)hs[v1]=0(1,v5)ht[v6]=0∞第2轮迭代(6,v3),(6,v5)hs[v1]=0hs[v4]=4ht[v6]=0ht[v5]=1∞第3轮迭代(6,v5)hs[v1]=0hs[v4]=4hs[v3]=6ht[v6]=0ht[v5]=1∞第4轮迭代hs[v1]=0hs[v4]=4hs[v3]=6hs[v5]=6ht[v6]=0ht[v5]=17
然而，针对地图匹配的情形，算法2仍有较大的优化空间.如图8(b)所示，假设GPS点pi有m个候选点，pi+1有n个候选点，在转移图中，pi到pi+1共有m×n条有向边，因此需要触发m×n次最短路径的计算.一种最直观的方式，是调用m×n次算法2，每次计算一组候选点对的最短路径.然而，这样会带来很多冗余的计算.给定如图9(a)所示的路网，节点从左往右层级逐渐升高.要计算A到D和F的最短路径，采用一对一的最短路径计算方法，其搜索过程如图9(b)所示，路径中的节点表示访问的先后顺序.可见，图9(b)中虚线部分的搜索过程是冗余的.
Fig. 9 Redundant computing for one-to-one shortest path
图9 单对单最短路径冗余计算
为了尽可能地减少冗余计算，本文提出了基于层次收缩的多对多最短路径搜索方法.其大体思想是，给定一个层次收缩路网G=(V,E)，出发节点集合S={s1,s2,…,sm}和目标节点集合T={t1,t2,…,tn}，我们依次从每个出发节点si和目标节点tj往外扩张一步.其中，对于每个出发节点，我们沿着有向边且层级增加的方向扩张；对于每个目标节点我们沿着有向边的反方向且层级增加的方向扩张.若某个出发节点(或目标节点)与所有目标节点(或出发节点)确定了最短路径距离，那么该出发节点(目标节点)便停止扩张.
算法3给出了基于层次收缩的多对多最短路径搜索算法的伪代码.其基本框架与算法2类似，不同之处在于，1)在初始化阶段，我们为每个起始点和终止点都分配一个优先队列和一个哈希表，当前解是一个m×n的数组Dmin；2)在循环迭代阶段，终止条件改成了2个优先队列集合是否均为空，即表示是否所有起止点之间都找到了最短路径.每次迭代时，我们先让所有起始点的扩张正向走1步，然后让所有终止点的扩张反向走1步.在正向扩张中，行⑨表示si节点已经找到了到所有目标节点的最短路径，因此优先队列qsi不应该再检查，直接从优先队列集合中删除.在行～，我们更新所有si和tj相遇时的最短路径.反向扩张中也作了类似的改动.3)最后返回的是所有起止点对的最短路径距离Dmin.
算法3. 基于层次收缩的多对多最短路径搜索算法.
输入：层次收缩路网G=V,E,起始点集合S={s1,s2,…,sm}，终止点集合T={t1,t2,…,tn};
输出：所有起止节点si,tj之间的最短路径长度
① 初始化2个优先队列集合QS={qs1,qs2,…,qsm}和QT={qt1,qt2,…,qtn}，qsi∈QS，qtj∈QT的排序键分别为对应节点到si和tj的最短路径距离;
② 初始化2个哈希表集合HS={hs1,hs2,…,hsm}和HT={ht1,ht2,…,htn}，hsi∈HS，hti∈HT的键为已检查的节点，值为该节点到si或tj的最短路径;
③ 初始化一个数值集合
每个元素初始化为∞.
表示si到tj当前状态下的最短路径距离;
④ 对于si∈S和tj∈T在同一路段的情形，直接计算si到tj的最短距离，更新
⑤ 将si∈S的所有出口节点加入qsi，将tj∈T的所有入口节点加入qtj；
⑥ while QS≠∅或者QT≠∅ do
⑦ for each qsi∈QS do
/*所有起始点正向走一步*/
⑧ (di,vi)←qsi.pop();
⑨ 
/*si正向扩张结束*/
⑩ 从QS中移除qsi;
else if hsi[vi]=null do
hsi[vi]=di;
for each htj∈HT且htj[vi]≠null
do
end for
for each(vi,vk)∈E且vk.level>
vi.level do
将(di+len(vi,vk),vk)加入qsi;
end if
end if
end for
for each qtj∈QT do
/*所有终止点反向走一步*/
(dj,vj)←qtj.pop();
/*tj反向扩张结束*/
从QT中移除qtj;
else if htj[vj]=null do
htj[vj]=dij;
for each hsi∈HS且hsi[vj]≠null
do
end for
for each(vk,vj)∈E且vk.level>
vj.level do
将(dj+len(vk,vj),vk)加入qtj;
end if
end if
end for
end while
return Dmin.
图10展示了利用Spark进行地图匹配的处理过程.在第3.1节，我们获得了分段后的轨迹RDD.在第3.2节，我们获得了多份收缩后的子路网，并存储在了外部磁盘中.在地图匹配过程中，我们首先计算分段轨迹的空间范围range1，然后在Spark的Driver进程中依次读取与range1相交的每个子路网文件，匹配完某个子路网内的轨迹后，销毁该路网变量，再读取下一个路网.由于每个子路网不会很大，因此可以在Driver进程的内存中存储起来.对于每个子路网，我们首先向集群中广播其空间范围.分段后的轨迹RDD调用filter方法，过滤掉路网空间范围外的轨迹数据.由于经历了轨迹分段和过滤，分区间的轨迹数目可能会相差比较大，因此，需调用repartition方法，一方面调整分区的数目，使得耗时的地图匹配能够充分利用集群资源，另一方面能缓解数据不平衡的问题.最后通过调用map方法结合广播的子路网对每条轨迹进行地图匹配.
Fig. 10 Distributed trajectory map matching based on Spark
图10 基于Spark的轨迹地图匹配
3.4 性能分析
轨迹预处理过程中，假设一条原始轨迹的长度为n，共有m条轨迹，驻留点检测的最差时间复杂度为
轨迹分段的时间复杂度为
因此轨迹预处理的时间复杂度为
然而，驻留点检测通常会更快，因为移动对象会不断移动，在寻找每个锚点的后继点时，不会扫描所有的后续节点.故最终的时间复杂度可以表示为
其中c1为常数.
假设原始路网中边和顶点的数目分别为|E|和|V|.在路网预处理过程中，优先队列底层的数据结构为最小堆，将节点入队实际上是构建最小堆的过程，采用自底向上方法构建最小堆的时间复杂度为
接着，需要对优先队列中的每个节点出队，单个节点出队后需要调整队列，每次调整的时间复杂度为
共|V|个节点，因此时间复杂度为
对每个节点编号后，需要判断当前节点v是否在其邻居节点的最短路径上.假设v的度为k，其邻居对的数目为k(k-1)，采用Dijkstra算法计算两点间的最短路径的时间复杂度为
因此计算v的邻居两两之间的最短路径时间复杂度为
然而实际中会远远小于此值，因为v的邻居节点u,w之间通常相距很近，且在搜索过程中，我们可以限定搜索的空间范围为len(u,v)+len(v,w)，每次搜索开销可以看成是一个常数c2.因此，最终路网预处理的时间复杂度为
地图匹配时，构建索引、查找候选点的开销可以忽略不计.最耗时的部分在于最短路径计算.假设共有m条轨迹，每条轨迹有n个GPS点，每个GPS点有k个候选点，则其时间复杂度为
与路网预处理中不同，这里最短路径搜索耗时取决于2个相邻GPS点的距离，因此不能忽略.
尽管从时间复杂度上看CHMM算法与相对于传统基于HMM的地图匹配算法相同，但是通过以下3个优化，大大提升了地图匹配效率：1)分布式计算充分利用了多台机器资源；2)通过层次收缩路网限制了最短路径的搜索方向，跳过了大量中间节点的访问；3)多对多最短路径搜索减少了大量冗余计算.
4 系统展示
本文提出的基于层次收缩的地图匹配算法已经落回了JUST系统[31-33,37]中，用户通过如下的SQL语句，就能够实现快速地图匹配功能.
SELECT
st_trajMapMatchToCompletePath(t1.traj,
t2.rn)
FROM trajTable AS t1,
(SELECT
st_makeRoadNetwork(collect_list(r))
AS rn
FROM rnTable) AS t2
其中，trajTable是分段后轨迹表的名字，rnTable是层次压缩后的路段表的名字，通过collect_list函数将路段收集到了Spark Driver程序中，st_makeRoadNetwork将路段组装成了路网rn.JUST会自动将rn广播到集群中，并触发分布式的地图匹配算法.
图11展现了JUST系统的用户界面，用户在SQL面板中输入上述地图匹配SQL语句，就能在地图面板中查看地图匹配的结果，如图中绿色线所示.本文提出的算法已经支持了危化品车辆行为分析[7]、地图路网修复[48]等多个项目的落地.
Fig. 11 UI of JUST for Map-Matching based on CH[38]
图11 JUST实现基于层次收缩的地图匹配界面[38]
5 实验及分析
本节先介绍实验相关设置，然后分析实验结果.
5.1 实验设置
5.1.1 数据集
我们采用4个真实的数据集来验证本文提出的地图匹配方法的性能.
1) 西安轨迹数据集.采用滴滴盖亚数据开放计划(2)https://outreach.didichuxing.com/提供的轨迹数据，原始数据共包含中国西安市从2018-10-01—2018-11-30这2个月的出租车轨迹数据，其空间范围为经纬度坐标(108.92,34.20)到(109.01,34.28)，包含了2 846 251 982条GPS点，未分段的轨迹有3 960 100条，磁盘占用为217 GB，相邻2个GPS点的采样时间间隔约为3 s.分段后平均每条轨迹包含42个GPS点.我们随机采样一些GPS点，数据分布如图12(a)所示，由图12可知，轨迹在空间上分布非常不均匀.
2) 西安路网数据.从OSM(open street map)(3)https://www.openstreetmap.org/网站中下载了西安轨迹空间范围向外扩了5 km的路网数据，这样确保边界处的轨迹匹配结果正确.截取的路网共包含了71 906条边，33 370个顶点，磁盘占用8.47 MB.
3) 北京轨迹数据集.T-Drive[49-50]数据集包含中国北京市从2008-02-02—2008-02-08一周内10 366辆出租车的轨迹数据，原始数据共有17 662 984个GPS点，平均采样时间间隔为177 s，磁盘占用为752 MB，去除明显噪声(坐标为(0,0)的GPS点)后，其空间范围为经纬度坐标(116.24,39.83)到(116.48,40.02)，分段后的轨迹数目为314 086条.图12(b)展示了T-Drive采样数据的分布.
4) 北京路网数据.从OSM中下载北京轨迹数据集对应空间向外扩张5 km范围内的路网数据，共有43 856条边，19 665个顶点，磁盘占用3.89 MB.
前2个数据集统称为西安数据集，后2个数据集统称为北京数据集.默认情况下，我们采用西安数据集，因为其轨迹数据量更大、采样时间间隔更小、数据质量更好.
Fig. 12 Data distribution
图12 数据分布
5.1.2 衡量标准
本文在保证与基于HMM的地图匹配算法准确率不变的情况下，主要聚焦在地图匹配效率的提升上.因此，我们主要比较不同方法的预处理时间、地图匹配的时间.
5.1.3 对比方法
我们比较以下方法的地图匹配效率.注意对于不同方法，使用的都是分段后的轨迹.所有对比方法的代码均随着我们的方法一起开源.
1) CHMM.本文提出的方法，即采用广播路网的方式、轨迹随机分区的方法、多对多的层次压缩最短路径算法来加速地图匹配过程.
2) FMM[22].FMM方法首先计算每条路段到一段范围内所有其他路段的最短路径，然后通过查表的方式来计算最短路径，从而加快地图匹配的效率.我们实现了其分布式版本.在预计算过程中，采用单源Dijkstra算法计算δ范围内的最短路径并存储在一个二维哈希表中.在地图匹配阶段，若在二维哈希表中存在2条路段的最短距离，则直接返回，若不存在，则使用Dijkstra算法计算最短距离.
3) PartitionSpace[44].文献[44]首先对轨迹数据进行采样，在Spark Driver中采用四叉树的方式对空间进行划分，使得每个子空间的轨迹样本数目相近.接着对全量轨迹数据集按照此空间进行划分，使得每个分区内的轨迹数目尽可能相同.它还考虑了边界路网的问题，对于每个分区，往外扩张500 m路网.文献[15]采用的是一对一Dijkstra计算最短路径，但为了公平起见，我们在每个分区内采用了CHMM算法，即这个方法与我们的方法不同之处在于数据分区方式的不同.
4) Dijkstra.采用传统的采用单向Dijkstra算法计算最短路径距离，例如文献[13,51].
5) Bi-Dijkstra[23].文献[23]采用双向单对单Dijkstra算法计算最短路径距离，用来加速地图匹配.
6) m-n-Dijkstra[24-25].文献[24-25]采用双向多对多Dijkstra算法计算最短路径距离，用来加速地图匹配.
7) A*.文献[39]采用A*算法计算最短路径距离.
8) CH.采用单对单CH算法计算最短路径距离.
注意，实验结果中未展示双向A*算法[46]，因为文献[52]指出，双向A*算法在双向搜索相遇时，得到的只是近似结果.为得到准确结果，仍会继续搜索很长的时间，导致算法极慢.实际上，我们也实现并开源了双向A*搜索算法，终止条件为：任何一个双向搜索开放队列为空，或者可能的最短路径大于当前已经找到的最优值.在我们的实验环境中，即使匹配一条只有14个GPS点、采样间隔为24 s的轨迹，双向A*算法花费的时间也超过了130 min.我们不比较多对多A*算法以及双向多对多A*算法，因为A*算法的启发式函数是针对单个起止点的，当涉及到多个起止点时，难以设计启发式函数，这将是一个有趣的研究点.
5.1.4 参数设置
表2显示了实验过程中的各参数变量取值，其中若无特殊说明，默认的参数用粗体表示.我们通过查看机器数目和轨迹集大小的变化验证我们方法的
Table 2 Parameter Settings
表2 参数设置
参数取值机器数目1台,2台,3台,4台,5台轨迹采样间隔3s,9s,24s,39s,60sGPS点个数5,10,15,20,25,30
注：默认参数用黑体表示.
高可扩展性；通过不同采样率和GPS点个数的轨迹验证我们方法的广泛适用性.对于所有方法，在寻找候选点时，我们先对路网构建R-tree索引.
5.1.5 实验环境
集群由5台物理机器组成，每台机器有24个物理CPU核、128 GB内存，操作系统为CentOS 7.4.我们部署了Hadoop 2.7.3和Spark 2.3.3，编程语言使用Java和Scala，JDK的版本为1.8.在实验过程中，我们分配了5核和5 GB内存给了Spark Driver进程，为每个机器启动了30个executor，每个executor分配了5核和16 GB内存.
5.2 实验结果
Fig. 13 Preprocessing performance of FMM (Xi’an)
图13 FMM预处理性能 (西安)
5.2.1 路网数据预处理性能对比
我们首先比较不同方法的路网数据预处理性能.注意我们没有比较轨迹的预处理性能，因为所有方法都是基于分段后的轨迹进行地图匹配的，他们的轨迹预处理过程相同.在这些方法中，只有CHMM,CH和FMM需要进行路网数据预处理，而CHMM和CH方法都是对路网进行同样的层次压缩，因此我们只比较CHMM与FMM的路网预处理性能.路网数据预处理都是在Spark Driver程序中进行的.图13展示了FMM算法对西安路网数据的每一条路段，采用Dijkstra算法预计算其到δ空间范围内所有路段间的最短路径所需的时间和结果占用空间.由图13可知，FMM算法的预处理时间和预处理结果大小都随着δ增大呈指数增长，当δ=3 km时，预处理结果大小达到了1.67 GB，当加载到内存时，占用空间会更大，不适合分布式场景下广播.而CHMM的路网预处理时间，需要的时间仅为13 s，压缩后序列化的路网大小仅为19 MB，因此适合广播.综上，CHMM方法仅需要少量的预处理开销，非常适合分布式的场景.
5.2.2 地图匹配性能随采样间隔变化
由于地图匹配非常耗时，在这一组实验中，我们选取了西安数据集预处理后前10 000条轨迹进行分布式匹配，用来比较不同方法的匹配时间随着采样时间间隔的变化.原始轨迹的采样时间间隔为3 s一个点，9 s表示每3个点保留一个GPS点，以此类推.
表3展示了实验结果，其中括号中的数字表示对应采样时间间隔下轨迹的平均GPS点个数.从表3中我们可以得出以下7个结论：
1) 随着采样时间间隔从3 s增加到60 s，几乎所有方法的地图匹配时间先减少，然后缓慢增加.因为采样时间间隔越大，一条轨迹中包含的GPS点越少，需要计算的最短路径点对的数目越少，因此地图匹配时间越少.但是随着采样间隔越大，2个相邻的GPS点距离也就越大，计算最短路径的时间也越长，导致采样时间间隔从24 s增加到60 s时，地图匹配时间有所增加.
2) 对于FMM算法，随着δ的增加，从缓存中能够查到的最短路径记录命中率也就越大，因此，δ=1 000 m时的地图匹配效率比δ=500 m时更快.但是在我们的实验环境中，当δ=2 000 m时，路网预计算结果在内存中的空间大于2 GB，Spark无法广播大于2 GB的变量，因此我们未比较δ≥2 000 m的情形.
3) 当采样间隔较低时，FMM算法比Dijkstra更快，因为我们可以在预计算结果中，直接查询到某些最短路径，从而避免耗时的Dijkstra计算.例如，当δ=1 000 m、采样时间间隔为3 s时，命中率约为60%，避免了大部分的计算开销.随着采样时间间隔增加，FMM的用时与Dijkstra的差距越小，因为命中率会变小.
4) Bi-Dijkstra算法、m-n-Dijkstra算法、A*算法、CH算法、CHMM算法都比原始的Dijkstra更快，而本文提出的CHMM最快，因为通过压缩路网，CHMM算法能够限定搜索的方向，通过捷径能够跳过一些中间节点，而且通过多对多搜索，减少了重复计算.
5) Bi-Dijkstra算法比A*算法快，甚至是Dijkstra算法性能的20倍左右，因为Bi-Dijkstra算法减少了许多无效的路网扩张；此外，在实现中，我们采用了文献[53]的方法，一旦双向扩张相遇后，只扩张队列内的路段，进一步限定了搜索范围.
6) m-n-Dijkstra与CHMM类似，共享了部分中间结果，减少了大量的重复计算过程，这非常适合基于HMM的地图匹配算法，因此A*算法不如m-n-Dijkstra算法.
7) PartitionSpace方法内部也采用了CHMM算法，但是采用四叉树根据轨迹分布对轨迹进行分区(在我们的实验中，最大分区数目为150，因为集群的总核数为150).由于轨迹的空间分区非常不均衡，通过采样对空间进行分区，仍会造成分区不均衡的问题，而本文提出的CHMM能够保证分区均衡，因此比PartitionSpace方法更快.
Table 3 Map-Matching Time VS Sampling Interval on Xi’an Datasets
表3 西安数据集地图匹配时间随着采样间隔变化情况 s
算法采样间隔(平均GPS点个数)3s(46)9s(23)24s(16)39s(14)60s(12)FMM(500m)28251291114812441448FMM(1000m)19989698548901115PartitionSpace10553394348Dijkstra34301556127412981446Bi-Dijkstra14075716778m-n-Dijkstra4432332536A*22381036884868918CH5936403140CHMM2927232324
注：黑体值为最优结果.
5.2.3 地图匹配性能随GPS点个数变化
在这组实验中，我们展示了采用北京数据集的实验结果.西安数据集的实验结果类似，因此结果未展出.我们选取了北京数据集中GPS点个数大于30的10 000条轨迹，每条轨迹的采样时间间隔大约是24 s.对于每条轨迹，分别截取前5,10,15,20,25,30个GPS点，形成新的轨迹，用来验证不同方法随着GPS点个数增加的耗时变化情况.
由表4可知，随着GPS个数的增加，所有方法所需要的耗时都增加，因为我们需要匹配更多的点，计算更多的最短路径.所有方法中，CHMM算法效率最高.m-n-Dijkstra算法表现也很不错，因为它共享了部分中间结果，减少了重复计算.在不需要多次执行地图匹配算法或者路网更新较频繁的场景，可以使用m-n-Dijkstra算法，因为它无需对路网数据预处理.有趣的是，FMM(500 m)比FMM(1 000 m)效率更快，这与常识以及表3的结果相反.经分析，本组实验中FMM算法的耗时都不高，因此其广播变量的时间不可以忽略.当δ=500 m时，预计算的最短路径结果的磁盘占用为14 MB；而当δ=1 000 m时，预计算的最短路径结果的磁盘占用为48 MB，所需的广播耗时更大.
Table 4 Map-Matching Time VS Number of GPS Points on Beijing Datasets
表4 北京数据集地图匹配时间随着GPS点个数变化情况 s
算法GPS点个数51015202530FMM(500m)3250647994110FMM(1000m)3857738899119PartitionSpace182239485760Dijkstra366189115142164Bi-Dijkstra192428343943m-n-Dijkstra171820212223A*3450647692104CH202225293135CHMM121417171921
注：黑体值为最优结果.
结合表3和表4来看，匹配北京数据集所需的时间比匹配西安数据集所需的时间小很多，因为北京数据集的路网更稀疏，每个GPS点的候选点更少，需要计算的最短路径更少.
各步骤的耗时占比.表5和表6分别统计了各个算法在地图匹配过程中，每个步骤的时间占比以及耗时情况.这里选用了西安数据集预处理后GPS点个数在30～50的1 731条轨迹，轨迹采样时间间隔为24 s.由于本组实验主要查看各部分耗时占比，我们采用单机版单线程环境运行了实验，各步骤的耗时是匹配这1 731条轨迹的时间之和.注意，PartitionSpace算法与CHMM算法只是分布式分区方法不同，在单机版单线程环境中，与CHMM算法一致，因此其结果未列出.所有的方法中，构建候选图的耗时最大，因为里面涉及到耗时的最短路径查询操作；确定匹配路径耗时最少(接近于0%)，这得益于Viterbi算法的使用.例如，对于传统的Dijkstra算法，其构建候选图的耗时高达99.19%，即使作了一些优化，最佳对比方法(即m-n-Dijkstra)构建候选图的耗时也需要占用92.65%.而CHMM是所有方法中构建候选图耗时唯一少于90%的方法，但也占用89.87%.将表3和表5结合起来看，构建候选图占比越小的方法整体运行得更快.因此，加快最短路径的计算对于基于HMM的地图匹配至关重要.另一个需要说明的现象是，表3与表6中各算法的性能之比是不同的，因为它们所使用的数据集不一样，而且执行环境不同，但他们之间的相对性能是一致的.
Table 5 Time Ratio of Different Steps in Map-Matching on Xi’an Datasets
表5 西安数据集地图匹配各步骤时间占比 %
地图匹配步骤算法FMM(500m)FMM(1000m)DijkstraBi-Dijkstram-n-DijkstraA*CHCHMM寻找候选点1.761.820.632.725.721.432.587.78构建候选图97.6397.7999.1996.6592.6598.2096.5289.87确定匹配路径0.010.010.000.010.010.000.010.03
Table 6 Time Cost of Different Steps in Map-Matching on Xi’an Datasets
表6 西安数据集地图匹配各步骤时间 s
地图匹配步骤算法FMM(500m)FMM(1000m)DijkstraBi-Dijkstram-n-DijkstraA*CHCHMM寻找候选点7.2885.9076.276.4735.786.376.3496.145构建候选图403.23317.83982.943230.02993.576437237.4371确定匹配路径0.0230.0170.0160.0150.0120.0050.0220.022总耗时41332599123810144524679
5.2.4 扩展性比较
为了验证CHMM方法的扩展性，我们做了2组实验.在第1组实验中，我们对西安数据集中100万条轨迹进行地图匹配，分析在不同集群节点数目情况下，所需要的匹配时间.我们从表3中选取了较快的4个算法CHMM,CH,m-n-Dijkstra,PartitionSpace进行对比，因为其他的算法效率太低，匹配100万条轨迹数据，即使使用5台节点的集群，运行一整晚均未出结果.由图14(a)可知，随着节点数目的增加，所有的方法所需的时间逐渐减少，然后趋于平稳.因为节点数目越多，计算资源越多，因此会更快，但当节点更多时，瓶颈将会由计算开销转化成通信开销，因此会趋于平稳.CHMM算法的运算效率比其他的方法更高，而且可以更快地稳定下来.
Fig. 14 Scalability of different methods (Xi’an)
图14 不同方法的扩展能力(西安)
图14(b)展示的是不同方法的执行时间随着数据集大小的变化情况.其中，100%的数据表示西安数据集中约为500万条轨迹.由图14可知，随着数据集的增大，所有方法所需的时间开销呈线性增长，但CHMM方法所需的时间少得多，而且增长率更低，CHMM优势更明显.整体而言，CHMM算法所需的时间是PartitionSpace算法的1/3.因此，CHMM的可扩展性更强，效率更高.
6 相关工作
对轨迹地图匹配问题的研究由来已久，其中基于HMM的轨迹地图匹配[13]应用最为广泛，因为其具有较高的准确率，能够处理噪声和低采样率轨迹的情形.但它也面临着效率较低的问题，因此，有大量加速轨迹地图匹配的研究工作.这些工作大体可分成4类：构建索引、简化计算、加速最短路径计算、并行计算或者分布式计算.
1) 构建索引.基于HMM的轨迹地图匹配的第一步是寻找每个GPS点周围一定范围内的路段，用于确定候选点.大多数地图匹配工作[13,39,44]都对路网构建空间索引，以加快对周围路段的查询.许多空间索引可供选择，例如：R-trees[14],Quad-trees[15],k-d trees[16]、网格索引或者它们的变种[54].尽管构建空间索引能够一定程度上加快地图匹配效率，然而，候选点的确定并不是基于HMM的轨迹地图匹配的最大瓶颈，仅仅通过空间索引对全局轨迹地图匹配的提速是非常有限的.
2) 简化计算.研究发现，基于HMM的轨迹地图匹配的瓶颈在于确定匹配路径[22]，其效率跟轨迹点的个数、路网的稠密度、最短路径的计算有关.为提升地图匹配的效率，直观方法是减少匹配GPS点的个数、减少候选路段的条数、简化最短路径的计算.①减少匹配GPS点的个数.当轨迹采样率非常高时，可以通过对轨迹进行降采样或者轨迹化简的方式，减少待匹配的GPS点从而加快地图匹配.文献[19]设计了一种考虑空间信息(例如路段长度和回转方向)的加权函数用来进行轨迹化简，并提出了3种轨迹化简方法来减少轨迹中的GPS点的个数.文献[20]首先考虑GPS点的移动方向对轨迹进行聚类，然后将同一个聚类内的GPS点匹配到同一条路段中，因为聚类的个数比GPS点的个数少很多，因此能够显著提升匹配效率，同时这种方法还能够缓解噪声点带来的影响.实时地图匹配[2,17-18]通常只考虑最近一段时间的GPS点数据，因此通常比离线地图匹配更快.②减少候选路段的条数.ST-Matching算法[39]只保留当前为止得分最多的l个候选点.SD-Matching算法[2]首先对每个候选点计算一个综合概率，它同时考虑了空间概率和角度概率的.然后只选择k个综合概率最大的点作为最终候选点，从而加快地图匹配效率.③简化最短路径的计算.SD-Matching算法[2]通过运动方向找到一个搜索区域，然后只考虑在搜索区域内的路径，大大简化了最短路径的计算.文献[17]针对手机基站数据，改造了HMM算法，提出了一个快速实时的维特比解码器，具有线性的空间复杂度和时间复杂度，其大体思路是消除了乘以0的路径组合.针对高频轨迹，文献[21]直接用GPS点的欧氏距离替代路网距离.简化计算的方法通常需要额外的信息，例如GPS点的方向，其通用性受限.此外，这类方法在一定程度上会降低匹配精度.
3) 加速最短路径计算.与简化计算不同，这类方法能保持匹配结果与原生HMM类的地图匹配算法一致.针对Dijkstra算法计算2个候选点的最短路径非常耗时的问题，文献[39]采用A*算法[46]或ALT算法[55]加速最短路径的计算.文献[23]采用双向Dijkstra算法来加速，而文献[24-25]采用多对多Dijkstra算法进行加速.前面这些提速方法需要在原始的路网上进行搜索，其加速效果仍然有较大提升的空间.文献[22]预计算并缓存了一定空间范围内任何一对路段之间的最短距离，当进行匹配时，通过查表的方式求得2个候选点之间的最短距离.这种方法非常消耗内存，在路网稠密的地方，常常容易内存溢出.
4) 并行计算或者分布式.另外一种方式通过使用更多的CPU来加速地图匹配过程.文献[51]采用多线程的方式并行化构建路网索引，同时并行化地图匹配过程，大大提升了匹配效率.文献[27]使用Hadoop[35]实现大规模GPS轨迹的分布式地图匹配.针对Hadoop执行效率较低的问题，涌现了很多基于Spark的分布式地图匹配工作[18,21,28-29,56].例如，文献[21]根据轨迹点采用了Quad-tree进行分区，同时，对于分区边界的路网，进行扩张和拷贝，从而提升边界点的匹配准确率.它还提出了一种分批处理的方法，防止数据过大时内存开销太大.但是根据轨迹进行分区的方法，在空间范围不大并且轨迹数目较多时，可能会面临着路网拷贝过多、单个分区内轨迹数目较多的问题，当划分空间非常小，轨迹也会被切分成很小，极端情况下，一个分区内的轨迹点数目只有1～2个，严重影响了地图匹配准确率.因此它的扩展性受限，效率不高，准确率可能会受影响.
本文提出的CHMM方法利用了空间索引加速候选点的查找，基于层次收缩的路网采用多对多的方式加速最短路径计算，同时提出了一个分布式处理框架，在保持高匹配率的同时，具有更高的效率和更强的可扩展性.
7 总结与展望
本文提出了一个基于层次收缩的分布式地图匹配算法框架，能够实现海量轨迹数据的快速匹配.具体来说，提出了一个简单但非常有效的数据分区策略，能够保证分布式环境下任务的负载均衡；提出了利用层次收缩的方法加速地图匹配效率，并在此基础上提出了基于层次收缩的多对多最短路径查询算法.实验结果表明我们提出的方法具有较高的效率和可扩展性.例如，在我们的实验环境中，对于500万条轨迹进行地图匹配，CHMM算法的耗时仅为PartitionSpace算法的1/3.所提的方法已经支持了多个真实空间数据智能项目的落地.我们还提供了一个在线演示系统，开源了核心代码，期待能够推动相关技术在更多的空间数据智能场景中落地.未来工作中，我们期待增加更多过滤条件在层次路网搜索中，同时考虑将本文提出的方法与其他优化方法例如文献[57]相结合，更进一步地提升地图匹配效率.
作者贡献声明：李瑞远负责想法的提出、文章的撰写以及rebuttal等事宜；朱浩文、王如斌负责执行实验、将算法落回产品；陈超负责论文算法讨论和整体质量把控；郑宇给我们提供了工作、实验环境，提出问题及大体方向，确定论文框架.
[1]Zuo Yimeng, Lin Xuelian, Ma Shuai, et al. Road network aware online trajectory compression[J]. Journal of Software, 2018, 29(3): 734-755 (in Chinese)(左一萌, 林学练, 马帅, 等. 路网感知的在线轨迹压缩方法[J]. 软件学报, 2018, 29(3): 734-755)
[2]Chen Chao, Ding Yan, Xie Xuefeng, et al. TrajCompressor: An online map-matching-based trajectory compression framework leveraging vehicle heading direction and change[J]. IEEE Transactions on Intelligent Transportation Systems, 2019, 21(5): 2012-2028
[3]Ma Liantao, Wang Yasha, Peng Guangju, et al. Evaluation of GPS-environment friendliness of roads based on bus trajectory data[J]. Journal of Computer Research and Development, 2016, 53(12): 2694-2707 (in Chinese)(马连韬, 王亚沙, 彭广举, 等. 基于公交车轨迹数据的道路GPS环境友好性评估[J]. 计算机研究与发展, 2016, 53(12): 2694-2707)
[4]Fang Zheng, Tong Guofeng, Xu Xinhe. A robust and effcient algorithm for mobile robot localization[J]. Acta Automatic Sinica, 2007, 33(1): 48-53 (in Chinese)(方正, 佟国峰, 徐心和. 一种鲁棒高效的移动机器人定位方法[J]. 自动化学报, 2007, 33(1): 48-53)
[5]Li Ruiyuan, Bao Jie, He Huajun, et al. Discovering real-time reachable area using trajectory connections[C] //Proc of the Int Conf on Database Systems for Advanced Applications. Berlin: Springer, 2020: 36-53
[6]Mao Jiangyun, Wu Hao, Sun Weiwei. Vehicle trajectory anomaly detection in road network via markov decision process[J]. Chinese Journal of Computers, 2018, 41(8): 1928-1942 (in Chinese)(毛江云, 吴昊, 孙未未. 路网空间下基于马尔可夫决策过程的异常车辆轨迹检测算法[J]. 计算机学报, 2018, 41(8): 1928-1942)
[7]Zhu Zheng, Ren Huimin, Ruan Sijie, et al. ICFinder: A ubiquitous approach to detecting illegal hazardous chemical facilities with truck trajectories[C] //Proc of the 29th ACM SIGSPATIAL Int Conf on Advances in Geographic Information Systems. New York: ACM, 2021: 37-40
[8]Li Ruiyuan, Ruan Sijie, Bao Jie, et al. Efficient path query processing over massive trajectories on the cloud[J]. IEEE Transactions on Big Data, 2018, 6(1): 66-79
[9]Li Ruiyuan, Ruan Sijie, Bao Jie, et al. Querying massive trajectories by path on the cloud[C] //Proc of the 25th ACM SIGSPATIAL Int Conf on Advances in Geographic Information Systems. New York: ACM, 2017: 1-4
[10]Zhang Junbo, Zheng Yu, Qi Dekang, et al. Predicting citywide crowd flows using deep spatio-temporal residual networks[J]. Artificial Intelligence, 2018, 259: 147-166
[11]Zhang Junbo, Zheng Yu, Qi Dekang, et al. DNN-based prediction model for spatio-temporal data[C] //Proc of the 24th ACM SIGSPATIAL Int Conf on Advances in Geographic Information Systems. New York: ACM, 2016: 1-4
[12]Yuan Haitao, Li Guoliang. A survey of traffic prediction: From spatio-temporal data to intelligent transportation[J]. Data Science and Engineering, 2021, 6(1): 63-85
[13]Newson P, Krumm J. Hidden Markov map matching through noise and sparseness[C] //Proc of the 17th ACM SIGSPATIAL Int Conf on Advances in Geographic Information Systems. New York: ACM, 2009: 336-343
[14]Guttman A. R-trees: A dynamic index structure for spatial searching[C] //Proc of the 1984 ACM SIGMOD Int Conf on Management of Data. New York: ACM, 1984: 47-57
[15]Finkel R A, Bentley J L. Quad trees a data structure for retrieval on composite keys[J]. Acta Informatica, 1974, 4(1): 1-9
[16]Bentley J L. Multidimensional binary search trees used for associative searching[J]. Communications of the ACM, 1975, 18(9): 509-517
[17]Algizawy E, Ogawa T, El-Mahdy A. Real-time large-scale map matching using mobile phone data[J]. ACM Transactions on Knowledge Discovery From Data (TKDD), 2017, 11(4): 1-38
[18]Wang Hongyu, Li Jin, Hou Zhenshan, et al. Research on parallelized real-time map matching algorithm for massive GPS data[J]. Cluster Computing, 2017, 20(2): 1123-1134
[19]Li Hengfeng, Kulik L, Ramamohanarao K. Spatio-temporal trajectory simplification for inferring travel paths[C] //Proc of the 22nd ACM SIGSPATIAL Int Conf on Advances in Geographic Information Systems. New York: ACM, 2014: 63-72
[20]Li Hengfeng, Kulik L, Ramamohanarao K. Robust inferences of travel paths from GPS trajectories[J]. Int Journal of Geographical Information Science, 2015, 29(12): 2194-2222
[21]Dogramadzi M, Khan A. Accelerated map matching for GPS trajectories[J/OL]. IEEE Transactions on Intelligent Transportation Systems, 2021 [2021-09-01]. https://ieeexplore.ieee.org/abstract/document/9312432
[22]Yang Can, Gidofalvi G. Fast map matching, an algorithm integrating hidden Markov model with precomputation[J]. International Journal of Geographical Information Science, 2018, 32(3): 547-570
[23]Koller H, Widhalm P, Dragaschnig M, et al. Fast hidden Markov model map-matching for sparse and noisy trajectories[C] //Proc of the 18th Int Conf on Intelligent Transportation Systems. Piscataway, NJ: IEEE, 2015: 2557-2561
[24]Wei Hong, Wang Yin, Forman G, et al. Fast Viterbi map matching with tunable weight functions[C] //Proc of the 20th Int Conf on Advances in Geographic Information Systems. New York: ACM, 2012: 613-616
[25]Jagadeesh G R, Srikanthan T. Fast computation of clustered many-to-many shortest paths and its application to map matching[J]. ACM Transactions on Spatial Algorithms and Systems (TSAS), 2019, 5(3): 1-20
[26]Ruan Sijie, Fu Xi, Long Cheng, et al. Filling delivery time automatically based on couriers’ trajectories[J/OL]. IEEE Transactions on Knowledge and Data Engineering, 2021 [2021-09-01]. https://ieeexplore.ieee.org/abstract/document/9497691
[27]Huang Jian, Qiao Shaoqing, Yu Haitao, et al. Parallel map matching on massive vehicle GPS data using mapreduce[C] //Proc of the 10th Int Conf on High Performance Computing and Communications & 2013 IEEE Int Conf on Embedded and Ubiquitous Computing. Piscataway, NJ: IEEE, 2013: 1498-1503
[28]Ruan Sijie, Li Ruiyuan, Bao Jie, et al. Cloudtp: A cloud-based flexible trajectory preprocessing framework[C] //Proc of the 34th Int Conf on Data Engineering (ICDE). Piscataway, NJ: IEEE, 2018: 1601-1604
[29]Francia M, Gallinucci E, Vitali F. Map-matching on big data: A distributed and efficient algorithm with a hidden Markov model[C] //Proc of the 42nd Int Convention on Information and Communication Technology, Electronics and Microelectronics (MIPRO). Piscataway, NJ: IEEE, 2019: 1238-1243
[30]Geisberger R, Sanders P, Schultes D, et al. Contraction hierarchies: Faster and simpler hierarchical routing in road networks[C] //Proc of Int Workshop on Experimental and Efficient Algorithms. Berlin: Springer, 2008: 319-333
[31]Li Ruiyuan, He Huajun, Wang Rubin, et al. JUST: JD urban spatio-temporal data engine[C] //Proc of the 36th Int Conf on Data Engineering (ICDE). Piscataway, NJ: IEEE, 2020: 1558-1569
[32]He Huajun, Li Ruiyuan, Bao Jie, et al. JUST-Traj: A distributed and holistic trajectory data management system[C] //Proc of the 29th ACM SIGSPATIAL Int Conf on Advances in Geographic Information Systems. New York: ACM, 2021: 403-406
[33]Li Ruiyuan, He Huajun, Wang Rubin, et al. Trajmesa: A distributed nosql storage engine for big trajectory data[C] //Proc of the 36th Int Conf on Data Engineering (ICDE). Piscataway, NJ: IEEE, 2020: 2002-2005
[34]Zaharia M, Chowdhury M, Das T, et al. Resilient distributed datasets: A fault-tolerant abstraction for in-memory cluster computing[C] //Proc of the 9th USENIX Symp on Networked Systems Design and Implementation (NSDI 2012). San Jose, CA: ACM, 2012: 15-28
[35]Dean J, Ghemawat S. MapReduce: Simplified data processing on large clusters[J]. Communications of the ACM, 2008, 51(1): 107-113
[36]Carbone P, Katsifodimos A, Ewen S, et al. Apache flink: Stream and batch processing in a single engine[J/OL]. Bulletin of the IEEE Computer Society Technical Committee on Data Engineering, 2015, 36(4). [2021-09-01]. https://reurl.cc/Gor3kD
[37]Li Ruiyuan, He Huajun, Wang Rubin, et al. TrajMesa: A distributed NoSQL-based trajectory data management system[J/OL]. IEEE Transactions on Knowledge and Data Engineering, 2021 [2021-09-01]. https://ieeexplore.ieee.org/abstract/document/9430714
[38]JUST. JD urban spatio-temporal data engine[OL]. [2021-07-30]. https://just.urban-computing.cn/ (in Chinese)(JUST. 京东城市数据时空引擎[OL].[2021-07-30] https://just.urban-computing.cn)
[39]Lou Yin, Zhang Chengyang, Zheng Yu, et al. Map-matching for low-sampling-rate GPS trajectories[C] //Proc of the 17th ACM SIGSPATIAL Int Conf on Advances in Geographic Information Systems. New York: ACM, 2009: 352-361
[40]Yuan Jing, Zheng Yu, Zhang Liuhang, et al. Where to find my next passenger[C] //Proc of the 13th Int Conf on Ubiquitous Computing. New York: ACM, 2011: 109-118
[41]He Tianfu, Bao Jie, Ruan Sijie, et al. Interactive bike lane planning using sharing bikes’ trajectories[J]. IEEE Transactions on Knowledge and Data Engineering, 2019, 32(8): 1529-1542
[42]Ruan Sijie, Xiong Zi, Long Cheng, et al. Doing in one go: Delivery time inference based on couriers’ trajectories[C] //Proc of the 26th ACM SIGKDD Int Conf on Knowledge Discovery & Data Mining. Virtual Event. New York: ACM, 2020: 2813-2821
[43]Karypis G, Kumar V. Analysis of multilevel graph partitioning[C] //Proc of the 1995 ACM/IEEE Conf on Supercomputing (Supercomputing’95). Piscataway, NJ: IEEE, 1995: 29-29
[44]Peixoto D A, Nguyen H Q V, Zheng Bolong, et al. A framework for parallel map-matching at scale using Spark[J]. Distributed and Parallel Databases, 2019, 37(4): 697-720
[45]Viterbi A. Error bounds for convolutional codes and an asymptotically optimum decoding algorithm[J]. IEEE Transactions on Information Theory, 1967, 13(2): 260-269
[46]Hart P E, Nilsson N J, Raphael B. A formal basis for the heuristic determination of minimum cost paths[J]. IEEE Transactions on Systems Science and Cybernetics, 1968, 4(2): 100-107
[47]Cormen T H, Leiserson C E, Rivest R L, et al. Introduction to Algorithms[M]. London, England: MIT Press, 2009
[48]Ruan Sijie, Long Cheng, Bao Jie, et al. Learning to generate maps from trajectories[C] //Proc of the AAAI Conf on Artificial Intelligence. Palo Alto, CA: AAAI Press, 2020, 34(1): 890-897
[49]Yuan Jing, Zheng Yu, Xie Xin, et al. Driving with knowledge from the physical world[C] //Proc of the 17th ACM SIGKDD Int Conf on Knowledge Discovery and Data Mining. New York: ACM, 2011: 316-324
[50]Yuan Jing, Zheng Yu, Zhang Chengyang, et al. T-drive: Driving directions based on taxi trajectories[C] //Proc of the 18th SIGSPATIAL Int Conf on Advances in Geographic Information Systems. New York: ACM, 2010: 99-108
[51]Song Renchu, Lu Wei, Sun Weiwei, et al. Quick map matching using multi-core cpus[C] //Proc of the 20th Int Conf on Advances in Geographic Information Systems. New York: ACM, 2012: 605-608
[52]Kaindl H, Kainz G. Bidirectional heuristic search reconsidered[J]. Journal of Artificial Intelligence Research, 1997, 7: 283-317
[53]Jin Xiaoqiang. Bi-directional Dijkstra algorithm and the acceleration of intermediate list[J]. Computer Simulation, 2004, 21(9): 78-81 (in Chinese)(靳晓强. 双向Dijkstra算法及中间链表加速方法[J]. 计算机仿真, 2004, 21(9): 78-81)
[54]Yang Keyu, Ding Xin, Zhang Yuanliang, et al. Distributed similarity queries in metric spaces[J]. Data Science and Engineering, 2019, 4(2): 93-108
[55]Goldberg A V, Harrelson C. Computing the shortest path: A search meets graph theory[C] //Proc of SODA. Philadelphia, PA: SIAM, 2005: 156-165
[56]Almeida A M R, Lima M I V, Macedo J A F, et al. DMM: A distributed map-matching algorithm using the MapReduce paradigm[C] //Proc of the 19th Int Conf on Intelligent Transportation Systems (ITSC). Piscataway, NJ: IEEE, 2016: 1706-1711
[57]Hu Gang, Shao Jie, Liu Fenglin, et al. If-matching: Towards accurate map-matching with information fusion[J]. IEEE Transactions on Knowledge and Data Engineering, 2016, 29(1): 114-127
