出租车轨迹数据地图匹配

时空轨迹数据是带有时间和空间信息的序列数据，数据获取成本低，覆盖范围广，且拥有时态特性，既可以进行微观个体活动模式的研究，也可以进行宏观活动系统的城市空间结构的研究。今天将分享一套出租车轨迹数据，该数据包含北京市2008年一周内10000+辆出租车、10分钟采样频率、1千万条以上的轨迹。

去看原文

时空轨迹数据的研究点

时空轨迹数据可进行城市动态景观分析，城市公共空间组织与规划，城市人口分布，城市交通可达性、土地利用与交通需求的依赖关系、城市路网脆弱性等研究。同时，利用配有GPS的浮动车数据(Floating Car Data)不仅能够动态追踪移动车辆的连续轨迹，而且由于浮动车辆的运动受限于路网，也能够动态感知城市道路的路网交通状态、拥堵状况、十字路口等重要交通节点的车流量状况，分析出租车轨迹的时空规律等。

例如，轨迹数据可用于提取人流、信息流的方向与强度，对于区域空间结构研究提供了新的视角。轨迹数据也可提取人口分布空间集聚动态规律，有助于交通问题的发现与解决，促进交通分析的精确化、科学化。

数据样例

1,2008-02-02 15:36:08,116.51172,39.92123
1,2008-02-02 15:46:08,116.51135,39.93883
1,2008-02-02 15:46:08,116.51135,39.93883
1,2008-02-02 15:56:08,116.51627,39.91034
1,2008-02-02 16:06:08,116.47186,39.91248

主要包含车辆ID，时间戳，经纬度等信息。

数据展示

路网及轨迹点展示：

一天的某量车辆轨迹动态展示：

某辆出租车多天轨迹2D平面展示：

某辆出租车多天轨迹3维展示：

地图匹配

选择部分轨迹序列，由于GPS的误差，出租车轨迹点通常不在路网中；然而，由于出租车是受路网约束的，因此需要进行地图匹配，即将在路网外的GPS点移动到路网上（不是简单的移动到最近的路段上），并且考虑轨迹的连贯性。

原始轨迹数据是不在路网上的：

匹配后的点一定是在路网上的。GPS点将在考虑前后轨迹点的同时，找到一个路段上最优的的一个候选点进行匹配（从下图可以看出并非直接移动到最近的路段上）：

地图匹配进行轨迹还原：

地图匹配准备

进行时空轨迹数据分析时，地图匹配是一项数据预处理部分重要的工作。此处地图匹配的算法可参考：如何把车辆轨迹数据匹配到路网？此匹配方法适用低频、高频轨迹数据地图匹配，并且考虑匹配轨迹的连贯性，而不是简单的将点移动到最近的路段上。此处给出地图匹配的Python实现和C#实现。

在进行地图匹配之前，首选需要准备路网数据（可从OSM下载），并对轨迹数据进行划分成不同的轨迹段。

以下是路网数据的预处理：

1. 路段筛选。从Open Street Map上下载的道路路网不仅包含车辆行驶道路，还包括自行车道、人行道、步行道等非机动道路。因此需要将这部分非机动车道路从路网中移除。
2. 路段划分。对于低等级的道路，需要在线的交点处进行划分，将路段分割成更小的单元。即在该交点处能够向左、右行驶，即真实道路十字路口的转向拓扑关系。即下图红色框内需要进行线分割。
   
3. 投影变换。从OSM下载的路网数据通常是WGS84地理坐标系，需要将其从WGS84地理坐标系投影为投影（平面）坐标系 ，以方便距离计算。
4. 生成网络数据集。提取交叉路口。
5. 空间连接并筛选。生产的网络数据集中存在junctions文件，进行空间连接，即将路网线数据和交叉路口连接起来。对于空间连接数大于2的，则用junctions进行打断。然后去掉空间连接数等于1的，即保证一条路段仅有两个junctions，即最小单元的路段的两端为交叉路口。
6. 生成拓扑表。根据空间连接结果生成路段两端点的坐标，长度等信息。然后生成拓扑表。

对于出租车轨迹数据，需要进行以下预处理：

1. 按车辆ID对数据进行划分，即同一辆车同一天的数据放在一起，并按时间戳字段进行排序；
2. 同车辆ID同天的数据进行轨迹段划分。因为出租车存在数据丢失或者一些其他原因，一辆出租车的轨迹数据并不是按采样间隔连续有数据，会存在很长一段时间段没有数据，因此需要进行轨迹段划分。对排好序的数据相邻两个轨迹点进行时间戳做差求时间间隔，若时间间隔大于给定的阈值，则认为是新的一段轨迹，并设置新的轨迹段ID。

处理好路网数据和轨迹数据，便可以进行地图匹配。

地图匹配Python实现

# -*- coding: utf-8 -*-
# @Author: 武辛
# @Email: geo_data_analysis@163.com
# @Note: 如有疑问，可加微信"wxid-3ccc"
def main():
  trajectory_data_rootdir = "轨迹数据根路径"
  network_path = "BJRoadNet/Beijing_Polyline.shp" # 路网shp文件路径
  bufferDis = 50 # 缓冲距离
  taix_id = 1 # 出租车辆id
  prex = "point_taxi_gps_20xx_by_id" # 保存文件的前缀

  # 加载轨迹数据以及基础数据
  datas = readTrajectory(trajectory_data_rootdir + "taxi_log_20xx_by_id/%s.txt" % (str(taix_id)))
  network = Network.loadNetwork("Beijing_NetworkLink", "NetworkLink")
  ArcPy_Util.bindFeatureLayer(network, network_path)

  isStop(datas, network_path, bufferDis) # 判断为静止的时段，就不考虑该轨迹点

  datas = datas[datas["Day"] == 2]
  datas.reset_index(drop = True)

  # 匹配部分代码
  candinate_points = match(network, datas, bufferDis)

  # 匹配结果整理，获取轨迹线、点信息
  print("Start to get the matched info...")
  route_polylines, route_points = getPolyline_Point_Info(candinate_points, network)

  # 匹配前结果存为shp文件
  points = [ArcPy_Util.lnglatToXY(coords["LNG"], coords["LAT"], ArcPy_Util.get_spatial_reference(network_path)) for index, coords in datas.iterrows() if not coords["Is_Stop"]]
  ArcPy_Util.save_points(points, ArcPy_Util.get_spatial_reference(network_path), "Matched_Taxi_Data/%s_%s.shp" % (prex, str(taix_id)))

  # 匹配后结果存为shp文件
  ArcPy_Util.save_points(route_points, ArcPy_Util.get_spatial_reference(network_path), "Matched_Taxi_Data/Mathched_%s_%s.shp" % (prex, str(taix_id)))
  ArcPy_Util.save_polylines(route_polylines, "Matched_Taxi_Data/%s_%s.shp" % ("polyline", str(taix_id)))

地图匹配C#实现

针对地图匹配的Python实现，采用了采样频率为5-10秒的轨迹数据进行了测试，路网规模的量级大概为（8094节点，18486条边），如下图，匹配6065个节点的时间大概需要3900+秒，可以看出这个匹配效率还是比较慢。

因此，本文采用C语言对该代码进行重写，并进行测试（同样采样频率数据，相同路网情况下）。如下图，可以看出，匹配100个点大概需要3秒，效率相对来说可以接受。

以下是某段轨迹的匹配结果：

T-Driver数据匹配结果输出

通过以上算法完成T-Driver数据的匹配，主要输出了GPS点划分不同轨迹匹配后的投影坐标以及相应的轨迹Shapefile文件，如下：


公众号后台回复“T-Driver”可获取匹配后的样例数据。

基于隐马尔科夫模型地图匹配实现

去看原文
以下是使用HMM匹配的结果：

V_ID,T_ID,LNG,LAT,TIMESTAMP,MATCHED_LNG,MATCHED_LAT
12,1,121.47647,31.281,2020-04-02 01:44:21,121.476453312,31.2801150372
12,1,121.47772,31.280,2020-04-02 01:44:31,121.477702482,31.2809243592
12,1,121.47836,31.282,2020-04-02 01:44:41,121.478281504,31.282192294
12,1,121.4785,31.283,2020-04-02 01:44:51,121.478443692,31.2835652779
12,1,121.4786,31.284,2020-04-02 01:45:05,121.478586156,31.285021706
12,1,121.47902,31.286,2020-04-02 01:45:11,121.47901228,31.2864233916
12,1,121.47945,31.288,2020-04-02 01:45:21,121.479431095,31.287804875
12,1,121.47994,31.289,2020-04-02 01:45:31,121.479948089,31.289306316
12,1,121.48073,31.290,2020-04-02 01:46:26,121.480746175,31.2906989314

更多内容，请关注地学分析与算法。