简介一套基于Python的武汉市出租车轨迹数据挖掘与分析项目源码主要面向计算机、交通、地理信息等相关专业的本专科生及研究人员可作为毕业设计、课程设计或科研预研的完整参考。项目实现了从原始轨迹数据读取、排序、插值平滑、路网匹配到轨迹可视化、上下车点热点分析、轨迹聚类、OD空间交互网络构建、目的地预测以及异常轨迹检测的完整技术链路并附有武汉行政区划矢量数据可支撑空间分析与结果展示。代码均经过运行测试按编号顺序组织配合README说明文档可快速复刻流程也便于在此基础上改造或扩展。压缩包共36个文件包含10个Python脚本、ArcGIS所需的shp/dbf/prj等地理数据文件及说明文档整体大小约11.58MB解压后即可查阅运行。目前已有213人学习下载适合希望系统掌握出租车轨迹挖掘方法并完成课设或毕设的开发者。1. 出租车轨迹数据挖掘这套流程到底在解决什么问题把某个“武汉市出租车轨迹”的压缩包解压后最先看到的通常是一堆 CSV 和 notebook。真正决定这套代码有没有用的不是聚类算法调得多漂亮而是从原始 GPS 报文到干净轨迹这条链路是否可靠。跑过这类数据的人都知道10 分钟原始数据里可能混入 2% 以上的坐标漂移部分点会落在楼顶、桥下甚至江面上不做地图匹配直接画热力图得到的是噪声热点而不是道路热点。基于 Python 的轨迹数据挖掘核心工作是把数据还原成车辆真实位移再在时间和空间两个维度上找出可运营的结论。对刚接触时空数据的人这套流程能少走半年弯路对已经在做网约车、物流调度的工程师也有参数和验证层面的参考价值。2. 武汉出租车轨迹的字段结构与数据清洗基线2.1 轨迹数据集的字段定义与质量陷阱武汉市出租车轨迹数据大多来自车载终端定时上报常见字段如下。字段类型含义常见脏数据vehicle_idstr车辆唯一编号前导零丢失、混入测试车timedatetime记录时间时区混用、时间戳乱序longitudefloatWGS84 经度0 值、超出武汉市域latitudefloatWGS84 纬度0 值、与经度不匹配passenger_statusint0 空车 / 1 载客状态跳变、全程为 1speedfloat瞬时速度 km/h负值、瞬时加减速异常directionint行驶方向角0-360 外取值落地第一步先确认坐标系。国内 GPS 设备多数直接记录 WGS84 经纬度但也有部分厂商做了 GCJ-02 偏移后入库。判断方法很简单取一个已知武汉地标坐标和车辆轨迹叠加到地图底图上看偏差。若偏移 300 米量级说明是加偏坐标后续要与高德、百度地图底图数据对齐时统一走对应坐标系的纠偏服务不要自己在经纬度上做线性缩放。这个细节不处理所有空间连接结果都会系统性偏移。质量陷阱集中在三处坐标落在境外、时间乱序、状态字段长期不变。坐标异常相对好处理按武汉地理范围筛掉即可时间乱序会在后续算速度、里程时制造负值必须先排序再处理状态字段长期不变说明车载计价器与 GPS 模块通信异常这类车要么剔除要么只保留其空间轨迹不参与载客行为分析。2.2 基于 pandas 的轨迹清洗最小实现清洗任务用 pandas 就够不需要一上来上 Spark。核心思路是统一列名、解析时间、按车排序、去重、范围过滤、算相邻点时间差与距离差。import pandas as pd import numpy as np df pd.read_csv( wuhan_taxi_sep01.csv, names[veh_id, tm, lng, lat, status, speed, direction], dtype{veh_id: str}, parse_dates[tm], ) # 统一排序同一辆车内按时间升序 df df.sort_values([veh_id, tm]).reset_index(dropTrue) # 同一辆车、同一时刻只留一条记录 df df.drop_duplicates(subset[veh_id, tm]).reset_index(dropTrue) # 按武汉市域大致范围过滤异常坐标 df df[ df[lng].between(113.7, 115.4) df[lat].between(29.9, 31.5) ].copy() # 速度合理性过滤 df df[(df[speed] 0) (df[speed] 150)].copy()read_csv里的names参数用于直接指定列名适用于原始文件没有表头的情况。parse_dates[tm]会尝试自动解析时间字符串如果格式是类似20230901143022的紧凑格式需要改成format%Y%m%d%H%M%S否则解析性能会很差且容易出错。排完序后去重才能保证后面用groupby.diff()算出的时间差是真实相邻点之间的值而不是被乱序样本污染的值。范围过滤的边界值按武汉市中心城区放宽到远城区不要卡得太死。武汉南北跨度约 160 公里东西也有近百公里过小的包围盒会误删天河机场、江夏、新洲方向的数据。2.3 载客状态识别与轨迹切分清洗后的数据仍然是一连串散点需要切分成“一段一段的行程”。这里有两个层级的概念一次载客行程通常由状态字段从 0 变为 1 开始到从 1 变为 0 结束一次连续上报指车辆熄火、信号丢失后重新出现的轨迹需要用时间间隔切分。实际操作中两步要同时做。# 相邻上报时间差 df[dt] df.groupby(veh_id)[tm].diff() # 空间距离用半正矢公式求相邻两点球面距离 def haversine(lng1, lat1, lng2, lat2): lng1, lat1, lng2, lat2 map(np.radians, [lng1, lat1, lng2, lat2]) a ( np.sin((lat2 - lat1) / 2) ** 2 np.cos(lat1) * np.cos(lat2) * np.sin((lng2 - lng1) / 2) ** 2 ) return 6371000 * 2 * np.arcsin(np.sqrt(a)) df[prev_lng] df.groupby(veh_id)[lng].shift() df[prev_lat] df.groupby(veh_id)[lat].shift() df[dist] haversine( df[prev_lng], df[prev_lat], df[lng], df[lat] ) # 时间间隔超过 5 分钟或距离超过 2 公里视为新轨迹段 df[seg_break] (df[dt] pd.Timedelta(minutes5)) | (df[dist] 2000) df[seg_id] df.groupby(veh_id)[seg_break].cumsum()groupby(veh_id)[seg_break].cumsum()之前要把每个车组的第一行seg_break视为 True否则一辆车最开始的记录会落到上一个车的段里。更稳妥的写法是直接加一个is_first df.groupby(veh_id).cumcount() 0把is_first也并进取seg_break。切完段之后记录状态翻转就是上下车事件df[prev_status] df.groupby(veh_id)[status].shift() # 状态从 0 变 1 是上车点从 1 变 0 是下车点 od_points df[df[status] ! df[prev_status]].copy() od_points[event] np.where( od_points[status] 1, pickup, dropoff )这里要警惕计价器状态本身切换延迟。很多定位终端在上车后几十秒才把状态置 1直接使用状态变化点做上下客坐标会有偏差。常见做法是把状态变化点往前回退一段距离或者与车速低点匹配载客开始前车辆通常处于怠速或缓行状态状态翻转点前面最近一个速度低于 5km/h 的点更接近真实上车位置。2.4 清洗结果校验先算指标再落盘清洗代码跑完不要直接进入绘图环节先打印一组质量指标把脏数据比例控制在可解释范围。total len(df) valid_ratio len(df) / total vehicle_count df[veh_id].nunique() sample_interval df[dt].dropna().median() zero_speed_ratio (df[speed] 0).mean() status_both ( df.groupby(veh_id)[status].nunique() 2 ).mean() print(f有效数据占比: {valid_ratio:.2%}) print(f车辆数: {vehicle_count}) print(f中位采样间隔: {sample_interval}) print(f零速占比: {zero_speed_ratio:.1%}) print(f出现过两种载客状态的车辆占比: {status_both:.1%})valid_ratio低于 90% 时优先检查坐标边界是否过紧sample_interval中位数一般在 10 到 60 秒之间如果大于 5 分钟说明数据本身就不是高频上报后续速度分析要降低时间分辨率zero_speed_ratio过高说明车辆长时间拥堵或终端在停车后仍持续上报聚类热点时要区分拥堵驻车点与真实上下客点。这一步还会暴露一个常见误用直接把所有点聚类而不区分载客状态。空车巡游点、载客行驶点、停车等待点是三种完全不同的行为混在一起做热点分析得到的结果无法指导调度。干净数据是后续所有分析的地基值得在这里多花时间。3. 轨迹地图匹配GPS 散点如何落回武汉路网3.1 为什么要先做地图匹配再谈分析清洗后的经纬度点仍然存在米级到十几米的误差在武汉这种高架、桥梁、隧道密集的城市误差会导致三类问题。第一车辆轨迹画出来后“飞”到楼宇或江面上视觉上不可信交付给业务方时解释成本极高第二计算路网速度时无法判断车辆在哪条路上汉口沿江大道和内部小路的车速会被混在一起拥堵识别失真第三OD 分析中上下客点落在道路外做“热点区域加周边道路”的关联时空间连接会错位。地图匹配就是给每一个 GPS 点找到它最可能所在的路段并给出一串连续的路段序列。它把“经纬度轨迹”变成“路网轨迹”之后的路径还原、行程时间预测、区域车速统计才具有道路语义。很多教程把地图匹配当成可选项实际做过真实城市数据的人都知道跳过这步等于把分析建立在错误的空间参照上。3.2 HMM 地图匹配的数学设定和关键参数常见的高精度方案是隐马尔可夫模型。GPS 观测点序列是显式的隐藏状态是车辆真实所在路段。模型要回答的问题是给定一条 GPS 轨迹最可能经过的路段序列是什么。观测概率描述 GPS 点与其候选路段之间的接近程度通常用点到路段的投影距离构造高斯分布。距离越近概率越高。转移概率描述车辆从上一时刻所在路段到达当前时刻候选路段的合理性由两条路段的道路网络最短路径距离与 GPS 点间直线距离的差异决定。两者差异越小转移概率越大。最终用 Viterbi 算法搜索全局最优路段序列。参数常见取值调参依据候选路段搜索半径120-200 米半径过小漏掉高架下层过大引入平行道路观测概率标准差15-30 米对应 GPS 典型误差水平转移距离惩罚系数0.3-0.5控制路径距离与直线距离差异的敏感度短时回退窗口3-5 个点用于处理定位丢星后的补偿常用在线算法取值更大高架桥是武汉地图匹配的最大难点。武汉长江大桥、二七长江大桥、天兴洲大桥上层与下层、城市快速路和地面辅路重叠GPS 点在垂直投影上可能同时靠近多条路。仅看单点距离必然选错HMM 的价值就在于通过前后位置关系纠正这种歧义。3.3 路网预处理与 Viterbi 解码实现路网数据一般用 OpenStreetMap 的武汉市 shp 或 pbf 文件。先做投影转换将经纬度坐标转到武汉市适用的 UTM 50N 投影之后所有距离计算都用平面欧氏距离避免每次算球面距离带来的性能开销。import geopandas as gpd road_gdf gpd.read_file(wuhan_road.shp) road_gdf road_gdf.to_crs(EPSG:32650) # 建立路网空间索引供候选路段搜索使用 road_gdf.sindex这是开源地图匹配库 Leuven.MapMatching 常见的输入形态。实际项目中我一般先用一个最近路段匹配作为快速基线def snap_points_to_road(points_gdf, road_gdf): points_utm points_gdf.to_crs(EPSG:32650) snapped [] for idx, point in points_utm.iterrows(): # 使用空间索引找最近的路段 nearest_idx road_gdf.sindex.nearest(point.geometry)[1][0] nearest_geom road_gdf.geometry.iloc[nearest_idx] snapped.append(nearest_geom.distance(point.geometry)) points_utm[dist_to_road] snapped return points_utm这个基线把点吸附到最近道路能解决大部分可视化问题但路径还原能力弱。要得到真正连续的路段序列需要走 HMM。核心解码逻辑不复杂可以直接实现def viterbi_decode(obs, states, emit_p, trans_p): obs: 每个时刻的 GPS 观测 states: 每个时刻的候选路段列表 emit_p[s][o]: 路段 s 产生观测 o 的概率 trans_p[s1][s2]: 从 s1 到 s2 的转移概率 T len(obs) S [len(candidates) for candidates in states] dp [{} for _ in range(T)] back [{} for _ in range(T)] for s in range(S[0]): dp[0][s] emit_p[0][s] back[0][s] -1 for t in range(1, T): for s in range(S[t]): candidates [ (dp[t - 1][prev] * trans_p[t][prev][s] * emit_p[t][s], prev) for prev in range(S[t - 1]) if trans_p[t][prev][s] 0 ] best_prob, best_prev max(candidates) dp[t][s] best_prob back[t][s] best_prev # 回溯最优路径 best_last max(range(S[T - 1]), keylambda s: dp[T - 1][s]) path [best_last] for t in range(T - 1, 0, -1): path.insert(0, back[t][path[0]]) return path代码中emit_p由 GPS 点到路段投影距离的高斯分布算出来trans_p由最短路径距离差算出来。candidates列表为空时要把对应的转移概率设为极小值而不是 0否则整条轨迹会断掉。实际工程中很少有人完整手写这个逻辑Leuven.MapMatching 和 GraphHopper 是常见选择但理解 Viterbi 回溯过程对调参数仍然十分必要。3.4 匹配失败看什么日志地图匹配跑完不要只看匹配率一个数。需要输出每段轨迹的候选路段数量、平均观测距离、路径平均长度这三个指标。候选路段数量为 0说明路网在该区域缺失优先补路网平均观测距离超过 50 米说明 GPS 点严重漂移需要回到清洗环节增强过滤匹配后路径长度比原始轨迹距离短三分之一以上往往是跨江隧道或高架桥的拓扑连接缺失。武汉特别容易出现隧道内丢星、过桥后跳变。遇到这种情况我一般把匹配结果按桥梁和隧道分段统计不能因为整体匹配率尚可就放弃细分场景的验证。地图匹配的输出要落成一张宽表包含轨迹 ID、路段 ID、进入时间、离开时间、覆盖长度这张宽表是后续所有速度分析与路径分析的数据源。4. 上下客热点聚类与运营指标的时空分析4.1 用 DBSCAN 提取载客热点参数与实现地图匹配完成后上下客事件点已经从原始坐标集里提取出来。接下来做需求热点分析DBSCAN 是轨迹数据挖掘里最常用的聚类算法因为它不需要预先指定簇数量也能自动识别噪声点。DBSCAN 两个核心参数是eps和min_samples。eps表示两个点被视为邻居的最大距离min_samples表示形成簇所需的最少点数。对于武汉市出租车上下客事件我一般把eps设在 200 到 300 米min_samples设在 5 到 10。250 米大致对应一个大型商场或地铁站的辐射范围过小会把同一商圈拆成多个碎片过大则把街道对面的不同需求点混为一谈。from sklearn.cluster import DBSCAN from sklearn.neighbors import NearestNeighbors import numpy as np # 将上下客点投到 UTM 投影后再聚类 coords_utm od_points[[x_utm, y_utm]].values # 用 K 距离图辅助选 eps k 8 nn NearestNeighbors(n_neighborsk).fit(coords_utm) distances, _ nn.kneighbors(coords_utm) k_dist np.sort(distances[:, -1]) # eps 取 K 距离图拐点处的值这里按实际数据调整 eps 250 model DBSCAN(epseps, min_samples8, metriceuclidean) od_points[cluster] model.fit_predict(coords_utm)K 距离图把每个点到第 8 近邻的距离排序后画出来曲线出现明显拐点的位置就是合理的eps。这比拍脑袋设参数可靠。聚类完成后cluster为 -1 的点是噪声代表零星需求不参与热点中心计算。热点中心用簇内点的几何中心即可如果要更精确可以用道路网络上的质心而不是经纬度平均值后者在道路弯曲处会偏移到路外。4.2 运营指标前置计算空驶率、带载率与有效里程聚类只能看出空间形态运营分析需要指标。三个最基础的指标分别是带载率、空驶率、平均载客里程。带载率是载客状态时间占总在线时间的比例反映司机实际有收入的时间占比空驶率是空车状态行驶里程占总行驶里程的比例反映巡游效率平均载客里程则反映单均需求距离与城市空间结构和计价规则直接相关。指标计算口径武汉市典型参考范围使用注意带载率载客时间 / 在线时间60%-75%不含停车等待时间会偏高空驶率空车里程 / 总里程30%-40%空车停车不计入里程平均载客里程载客总里程 / 载客次数5-8 公里受机场、火车站长单影响计算时注意空驶率的分母是空车行驶里程与载客行驶里程之和而不是在线时长。出租车在路边怠速等客时速度接近 0用时间算出的空驶率会偏低。里程指标由每段轨迹的距离加总得到已经是地图匹配后的路网距离比原始 GPS 点连线的直线距离多出 8% 到 15% 的路程这个差异在商业报告里必须说明清楚。trip_summary ( od_points.loc[od_points[event] dropoff] .groupby(seg_id) .agg( veh_id(veh_id, first), pickup_time(tm, first), dropoff_time(tm, last), trip_dist_km(trip_dist_m, sum), ) ) trip_summary[trip_duration_min] ( trip_summary[dropoff_time] - trip_summary[pickup_time] ) / pd.Timedelta(minutes1)这个聚合结果要回填到地图匹配输出的路段宽表上。因为匹配宽表已经有每段轨迹在每条路段上的行驶距离按轨迹 ID 汇总就能得到一条完整行程的实际路网里程远比用原始点经纬度连线的累计值可靠。4.3 分时段热力图输出轨迹数据的价值在时间维度上要有区分全天合并的热点图会把早高峰通勤和夜间娱乐需求混在一起建议按小时切片输出。常见的输出方式是 Folium 热力图直接叠加在 OpenStreetMap 底图上。import folium from folium.plugins import HeatMap # 取晚高峰 18-19 点的上车点 evening_mask (od_points[tm].dt.hour 18) (od_points[event] pickup) evening_pts od_points.loc[evening_mask, [lat, lng]].values m folium.Map(location[30.59, 114.30], zoom_start12) HeatMap(evening_pts, radius18, blur10, max_zoom15).add_to(m) m.save(wuhan_18h_pickup_heatmap.html)radius控制热力点的影响半径blur控制颜色过渡的柔和度这两个参数需要按数据量调整。上车点数量很大时radius可以缩小到 12 到 15避免热力块连成整片红色。生成 HTML 后用浏览器打开方便业务同事直接查看不需要 GIS 软件。分时段热点图与聚类结果要互为印证。如果 DBSCAN 在汉口火车站附近聚出一个簇但 18 点热力图显示该簇几乎不亮说明该热点主要由凌晨或午间需求贡献调度资源就不能堆在晚高峰。这种“聚类给出区域热力图给出时间”的组合分析是轨迹数据挖掘交付中比较实用的模式。5. 从热点到 OD 流用网格聚合验证调度价值热点回答了“哪里需求多”但没有回答“乘客从哪里来到哪里去”。对于出租车调度真正有决策意义的是 OD 流早高峰从住宅区集中流向办公区晚高峰反向流动夜间从商圈流向居住区。用网格聚合将武汉城区划分为约 1 公里见方的单元统计每个网格对之间的上下车转移量就能得到一张 OD 转移矩阵。# 按经纬度格网给上下车点打格 ID def cell_id(lng, lat, x0113.7, y029.9, size0.01): return f{int((lng - x0) // size)}_{int((lat - y0) // size)} od_points[cell] od_points.apply( lambda r: cell_id(r[lng], r[lat]), axis1 ) # 同一行程内上车格与下车格配对 trip_cells ( od_points.loc[od_points[event] pickup, [seg_id, cell]] .merge( od_points.loc[od_points[event] dropoff, [seg_id, cell]], onseg_id, suffixes(_from, _to), ) ) od_matrix pd.crosstab(trip_cells[cell_from], trip_cells[cell_to])cell_id中的size0.01约等于 1 公里这个粒度适合观察跨区流动。网格过大OD 流向被平均化看不出通勤廊道网格过小OD 矩阵稀疏度极高统计意义不足。武汉这种江水分隔的城市OD 流分析还要额外标注跨越长江、汉江的河段流向因为跨江行程的时间和费用显著高于同岸行程其需求弹性也不同。OD 矩阵做好后可以用一个简单指标验证调度价值返程空驶率。把 OD 矩阵拆成早高峰 7 点到 9 点与晚高峰 17 点到 19 点两个切片分别计算每个网格的净流入量净流入大的区域在对应时段必然需要更多空车从外部调入。若某个居住密度很高的网格晚高峰净流入小、但上车量极大说明这个区域是典型的“单向外送型”需求司机在高峰前主动空放过去反而是最优策略。# 记录每个行程的上车时间所在小时 trip_cells[pick_hour] od_points.loc[ od_points[event] pickup, tm ].dt.hour.values morning trip_cells[trip_cells[pick_hour].between(7, 9)] m_flow pd.crosstab(morning[cell_from], morning[cell_to]) morning_net_inflow m_flow.sum(axis0) - m_flow.sum(axis1) print(morning_net_inflow.sort_values(ascendingFalse).head(10))净流入前 10 的网格就是早高峰的车辆汇集目标。这套逻辑同样可以用于网约车平台的调度预先知道 OD 流方向就能减少司机的盲目巡游里程。最后把 OD 矩阵导出为 Parquet 并按小时分片后续要做周维度对比或接入可视化前端时读取性能会明显优于 CSV。本文还有配套的精品资源点击获取