从GPS轨迹到热点与OD：Python出租车数据挖掘全流程解析

简介基于Python实现的武汉市出租车轨迹数据挖掘与分析项目代码已经过测试可运行适合计算机相关专业学生用于毕业设计、课程设计或项目实训。压缩包共36个文件大小约11.58MB以10个Python脚本为核心覆盖轨迹读取排序、插值平滑与路网匹配、上下车点可视化与热区分析、轨迹聚类、OD热点空间交互网络构建、目的地预测及异常轨迹分析等步骤另含武汉行政区划的shapefile空间数据和少量说明文档方便在真实地理数据上进行空间统计。通过阅读代码能完整掌握出租车GPS轨迹从原始清洗到业务建模的流程例如如何识别客流热点、构建OD网络和检测异常轨迹。用户下载后可按目录顺序演练也可二次修改用于自己的实验或毕业设计。目前已有213人学习下载整体体量适中是快速上手轨迹数据挖掘与空间分析的良好参考。1. 为什么要专门写一套出租车轨迹挖掘而不是调库跑通如果只是把出租车GPS数据读进DataFrame再画几张散点图两个小时就够了。但真实场景里武汉这样千万级人口城市的出租车轨迹数据动辄上亿条单日记录就可能超过2000万行叠加信号漂移、停车等待、跨江大桥信号遮挡、司机交接班等因素后原始数据基本不能用。所谓数据挖掘不是跑一个聚类算法看个热闹而是先把轨迹还原成可用信息载客点在哪、空驶去了哪、哪条路几点最堵、司机交接班对载客点分布的影响有多大。这篇文章从数据结构出发完整走一遍清洗、停留识别、轨迹压缩、特征工程、热点挖掘和OD分析全部用Python实现代码可以直接在本地Jupyter或VS Code里跑起来。适合刚接触时空数据挖掘的Python使用者也适合有数据分析经验但没处理过轨迹数据的人——你真正缺的不是算法而是轨迹数据特有的预处理和验证手段这是普通DataFrame教程里不会讲的。2. 轨迹数据的原始形态与三个跑通全流程的Python环境准备2.1 出租车轨迹记录里到底有哪些字段可以用武汉出租车轨迹数据最常见的来源是车辆GPS终端按固定频率上报位置典型字段如下。不同数据源字段名略有差异但核心信息一致。import pandas as pd # 武汉出租车轨迹典型字段按真实数据常见顺序排列 sample_cols [ vehicle_id, # 车辆唯一标识 time, # 记录时间格式一般为 YYYY-MM-DD HH:MM:SS lng, # WGS-84 或 GCJ-02 坐标系下的经度 lat, # 纬度 speed, # 瞬时速度单位 km/h direction, # 行驶方向角0-360度 passenger, # 载客状态0空车 1载客 mileage, # 累计里程部分数据源有 ] df pd.read_csv(wuhan_taxi.csv, namessample_cols, dtype{ vehicle_id: str, passenger: int8, }, parse_dates[time]) print(df.shape) print(df.dtypes)这里有个关键点vehicle_id必须读成字符串否则前导零会丢失time直接用parse_dates转换避免后面反复调用pd.to_datetime拖慢速度。passenger用int8足够因为取值只有0和1。读入后先别急着清洗花30秒看一下数据范围。print(df[time].min(), df[time].max()) print(df[lng].describe()) print(df[lat].describe())从经度范围可以直接判断坐标是否合理。武汉市中心大约在东经114.2°到114.5°、北纬30.4°到30.7°之间如果数据里出现经度小于73或大于135、纬度小于3或大于54的记录说明混入了脏数据或火星坐标偏移未处理需要在下游处理中排除。2.2 处理亿级轨迹的Python环境要怎么组织处理亿级数据Pandas单机跑得动但要用对方式。三个要点直接决定效率第一用Parquet或Feather格式落地中间结果不要反复读写CSV第二清洗和特征工程尽量用向量化操作避免逐行iterrows()第三安装GeoPandas用于地理空间索引安装scikit-learn用于聚类安装OSMnx或Shapely用于路网匹配。下面是推荐的环境组织。conda create -n trajectory python3.10 -y conda activate trajectory conda install -c conda-forge geopandas scikit-learn shapely pyproj -y pip install pandas pyarrow osmnx matplotli提示轨迹数据涉及经纬度投影计算pyproj用于坐标系转换GeoPandas在空间过滤和栅格聚合时性能远超纯Pandas循环。如果只需要最简功能可以暂时不装OSMnx用GeoPandas自带的sjoin也够。环境准备好后先把原始数据按天切分落盘这是一种常见做法能显著降低后续多次全表扫描的代价。df[date] df[time].dt.date for d, group in df.groupby(date): group.drop(columns[date]).to_parquet(fdata/{d}.parquet, indexFalse)按天切分的好处有两个一是单日数据量大约百万到千万行内存可控二是后面做OD分析、热点挖掘时天然按天聚合更符合业务习惯。切分后删掉原始DataFrame释放内存import gc del df gc.collect()2.3 坐标系不一致会让聚类结果偏离几百米先统一再挖掘轨迹数据最容易踩的坑是坐标系混用。国内GPS设备输出的原始坐标是WGS-84但部分数据源为了合规会转成GCJ-02火星坐标有的地图厂商SDK直接输出BD-09百度坐标。三者偏差大约在几十米到几百米之间这在城市道路级别的挖掘中足以让一个路口的热点被分到旁边街区。判断当前坐标系的简单方法是取一个已知地标对照比如武汉长江大桥中心约在东经114.2848°北纬30.5463°如果数据在该点位附近偏移约0.0001到0.0005度大概率是GCJ-02未转换。用pyproj统一转换。from pyproj import Transformer # 假设原始数据是GCJ-02统一转成WGS-84 trans Transformer.from_crs(gcj02, epsg:4326, always_xyTrue) df[lng], df[lat] trans.transform(df[lng].values, df[lat].values)如果数据本身就是WGS-84此步骤可跳过。但要注意如果后续路网匹配使用OSMnx或公开路网数据路网通常也是WGS-84必须保证轨迹和路网在同一坐标系下。3. 轨迹预处理的四个必要环节清洗、停留识别、轨迹切分、压缩3.1 清洗不是删空值那么简单要看漂移点和速度突变轨迹数据的异常分为三类超出城市范围的离群点、由GPS漂移造成的瞬时跳变、由设备故障导致的重复或乱序时间戳。第一类直接用经纬度边界过滤第三类按下述方式删除。import numpy as np # 基于经纬度边界过滤武汉城区范围 mask ( (df[lng] 113.8) (df[lng] 114.8) (df[lat] 30.2) (df[lat] 30.9) ) df df[mask] # 删除同一车辆同一时间戳重复记录 df df.drop_duplicates(subset[vehicle_id, time]) # 按车辆和时间排序保证轨迹数据时间顺序一致 df df.sort_values([vehicle_id, time]).reset_index(dropTrue)速度突变需要记住一个关键点不要单看绝对速度。城市出租车在桥梁上和隧道口速度差距极大但GPS在隧道内可能产生几十秒的定位漂移瞬时速度会跳到120km/h以上。过滤规则用「相邻点距离除以时间差得到的平均速度」更合理这种速度叫载体速度能反映两个点间的真实行驶快慢。# 计算相邻点位移和方向角排除因漂移产生的虚假高速 df[prev_lng] df.groupby(vehicle_id)[lng].shift(1) df[prev_lat] df.groupby(vehicle_id)[lat].shift(1) df[prev_time] df.groupby(vehicle_id)[time].shift(1) # 用Haversine公式计算距离此处用简化球面距离 R 6371000.0 dlat np.radians(df[lat] - df[prev_lat]) dlng np.radians(df[lng] - df[prev_lng]) a np.sin(dlat / 2) ** 2 np.cos(np.radians(df[prev_lat])) * np.cos(np.radians(df[lat])) * np.sin(dlng / 2) ** 2 df[dist] 2 * R * np.arcsin(np.sqrt(a)) # 时间差转小时速度单位 km/h df[dt_hour] (df[time] - df[prev_time]).dt.total_seconds() / 3600.0 df[moving_speed] df[dist] / 1000.0 / df[dt_hour].replace(0, np.nan) # 载体速度超过 200km/h 视为漂移删除该段连接 df.loc[df[moving_speed] 200, [prev_lng, prev_lat, prev_time]] np.nan删除漂移点不是直接把速度异常的行删掉因为前后两段轨迹可能都是正常的错在当前点连接了不该相连的两点。处理方式是置空prev字段这样后续按空值切分轨迹时轨迹就会在这个位置断开而不是把两段不相干路线强行连在一起。3.2 载客状态切换是天然的轨迹切分点出租车的载客状态passenger字段是轨迹切分的天然标注。一次行程的起点是passenger从0变1的记录终点是1变0的记录。先找到所有状态切换点再切出单次行程这种做法的好处是行程定义和业务口径一致后面做的起讫点OD分析、行程时长统计都不需要重新定义。# 找到每个车辆的载客状态变化点 df[pkg_prev] df.groupby(vehicle_id)[passenger].shift(1) df[trip_start] (df[pkg_prev] 0) (df[passenger] 1) df[trip_end] (df[pkg_prev] 1) (df[passenger] 0) # 对每条行程分配唯一ID df[trip_id] df[trip_start].cumsum() # 行程结束时下一段行程开始前的记录都属于当前行程 df.loc[df[trip_end], trip_id] df[trip_start].cumsum() df.loc[df[passenger] 0, trip_id] np.nan注意一个容易出错的地方如果原始数据没有passenger字段比如只有空驶记录怎么办。这时退而求其次用停留检测切分轨迹。停留意味着车辆长时间低速且位移很小一般是等人或交接班将停留点前后的轨迹视为两次行程。停留检测做法# 基于时间窗和距离阈值识别停留点 df[pt_delta] df[dist].fillna(0) df[pt_time] (df[time] - df[prev_time]).dt.total_seconds().fillna(0) # 滑动窗口内累计位移小于50米且持续超过180秒标记为停留 def mark_stay(g): stay np.zeros(len(g), dtypebool) i 0 while i len(g): j i while j len(g) and g[pt_time].iloc[i:j1].sum() 180: j 1 seg g.iloc[i:j] if seg[pt_delta].sum() 50 and seg[pt_time].sum() 180: stay[i:j] True i max(j, i1) return stay df[is_stay] df.groupby(vehicle_id, group_keysFalse).apply(mark_stay)这段代码的核心逻辑是滑动累积时间不是固定窗口。固定窗口会漏掉长时间低速缓行的情况而这里只要窗口中累计位移很小且持续3分钟以上就判为停留。停留段落在后续轨迹压缩中要保留原始点因为它们包含载客点位置的语义信息。3.3 DTW和道格拉斯-普克压缩轨迹挖掘里的取舍原始轨迹点太密直接用原始点做聚类或路网匹配计算量巨大且噪声点多。常见的轨迹压缩算法有道格拉斯-普克Douglas-Peucker和基于时间阈值的均匀抽稀。道格拉斯-普克保留形状特征适合几何形态分析时间抽稀则适合速度、拥堵类分析——因为你需要等间隔时间戳而不是等距离点。以下用递归实现道格拉斯-普克先算每点到首尾连线的垂直距离。from typing import List, Tuple def douglas_peucker(points: List[Tuple[float, float]], epsilon: float) - List[Tuple[float, float]]: if len(points) 2: return points # 首尾连线 x1, y1 points[0] x2, y2 points[-1] # 垂直距离向量化计算 dx x2 - x1 dy y2 - y1 if dx 0 and dy 0: dists np.sqrt((np.array(points)[:, 0] - x1)**2 (np.array(points)[:, 1] - y1)**2) else: dists np.abs((dy * (np.array(points)[:, 0] - x1) - dx * (np.array(points)[:, 1] - y1)) / np.sqrt(dx*dx dy*dy)) idx_max np.argmax(dists) if dists[idx_max] epsilon: left douglas_peucker(points[:idx_max1], epsilon) right douglas_peucker(points[idx_max:], epsilon) return left[:-1] right else: return [points[0], points[-1]]实际使用时不要对整条轨迹一次压缩而是先按轨迹ID或行程ID分组逐一压缩。GROUP_EPS 0.0005 # 约50米车辆GPS精度在此量级可接受 def compress_group(group): pts list(zip(group[lng], group[lat])) pts_c douglas_peucker(pts, GROUP_EPS) return pts_c # 仅对载客轨迹压缩 df_trip df[df[passenger] 1].copy() compressed df_trip.groupby(trip_id).apply(compress_group)压缩比一般在80%到95%之间即原始1000个点最后只剩50到200个点但路线形状和载客点位置保持不变。压缩后的轨迹才适合做全局聚类和OD矩阵否则1亿点跑DBSCAN的内存消耗直接让机器卡死。4. 轨迹特征工程与两个高价值挖掘任务热点区域聚类和OD时空分析4.1 轨迹点的高维特征不只是经纬度和速度原始轨迹只有坐标、时间、速度和方向但真正能用于预测或分析的标签性特征需要手工构造。下表是轨迹挖掘任务里最常用的一批衍生特征特征名计算方式适用场景时间槽按小时切分(0-23)时段热度分析方向变化率相邻点方向角差值的绝对值识别绕路、急转弯单位时间载客次数按车辆一小时内的载客状态切换次数司机运营效率载客率载客时间/总工作时长司机收入分析速度分位数单段轨迹速度的P50/P90拥堵判别、路线规划乘客上车点经度纬度行程起点坐标热点区域、候车点选址行程直线距离起点到终点Haversine距离OD分析、运力调度方向变化率的计算要小心角度跨越360度的边界比如从350度变到10度差值应该是20度而不是340度。下面这段代码处理了这个问题。# 当前方向角与下一方向角差值 df[dir_diff] np.abs(df[direction].diff()) df[dir_diff] np.where(df[dir_diff] 180, 360 - df[dir_diff], df[dir_diff]) # 按行程ID聚合得到急转弯次数 turn_features df_trip.groupby(trip_id)[dir_diff].agg([mean, max, lambda x: (x 30).sum()]) turn_features.columns [turn_mean, turn_max, turn_count]注意方向角在静止停放时可能剧烈跳动因此要先剔除速度接近0的点再算方向变化率否则噪声直接淹没信号。4.2 用DBSCAN找出热点区域参数要按武汉的街道尺度设热点区域挖掘是出租车轨迹数据挖掘最经典的任务之一输出的是一批中心点经纬度和覆盖半径可以直接用于调度、广告投放选址、候车区规划。主流算法是DBSCAN聚类它不需要预先指定簇数量能识别噪声点这对轨迹数据非常友好。武汉热点区域聚类前先选特征只用上车点经纬度然后用DBSCAN。from sklearn.cluster import DBSCAN # 提取上车点 pickup_points df_trip[df_trip[trip_start] True][[lng, lat]].values # 用 haversine 距离度量替代欧氏距离 from sklearn.metrics import pairwise_distances import math def haversine_dist(X, Y): R 6371000 X np.radians(X) Y np.radians(Y) dlat Y[:, 0][:, None] - X[:, 0] dlng Y[:, 1][:, None] - X[:, 1] a np.sin(dlat/2)**2 np.cos(X[:, 0])[None, :] * np.cos(Y[:, 0])[:, None] * np.sin(dlng/2)**2 return 2 * R * np.arcsin(np.sqrt(a)) D haversine_dist(pickup_points, pickup_points) db DBSCAN(eps150, min_samples20, metricprecomputed) labels db.fit_predict(D) # 聚类中心用簇内点均值 hotspots [] for lab in set(labels): if lab -1: continue pts pickup_points[labels lab] hotspots.append({ cluster_id: lab, lng: np.mean(pts[:, 0]), lat: np.mean(pts[:, 1]), count: len(pts) }) hotspot_df pd.DataFrame(hotspots).sort_values(count, ascendingFalse)这里两个参数必须说明。eps150表示聚类半径150米武汉市区出租车上车点密度高150米约等于一个街区的尺度能精确到路口级别如果设成500米会把多个路口合并成一个热点失去调度参考价值。min_samples20是最小簇点数全天数据下20个点可以过滤偶发上下客但如果你分析的是晚高峰1小时的数据20个点要求过高热点会过少。一个常见做法是分时段跑DBSCAN而不是全量跑一次早高峰用min_samples10平峰用min_samples20晚高峰用min_samples30。时段划分代码df_trip[hour] df_trip[time].dt.hour for hour, grp in df_trip.groupby(hour): pts grp[grp[trip_start] True][[lng, lat]].values if len(pts) 50: continue D haversine_dist(pts, pts) db DBSCAN(eps150, min_samples15 if hour in [7,8,9,17,18,19] else 20, metricprecomputed) labels db.fit_predict(D) # 打印每个时段簇数量与点数 print(hour, len(set(labels)) - (1 if -1 in labels else 0), (labels ! -1).sum())聚类完成后需要验证结果。常见做法是回到地图上随机抽样50个簇人工查看是否落在商圈、地铁站或大型社区门口。如果大量簇落在立交桥或江河中央先回去检查坐标转换和漂移点清洗而不是调聚类参数。4.3 OD矩阵和跨江需求武汉特有的空间约束分析武汉被长江和汉江分割成三镇跨江出行是极其重要的需求形态。OD矩阵分析需要把上车点和下车点映射到空间网格再统计网格间的流量。网格大小采用0.01度×0.01度约1公里见方比较合适过小则矩阵稀疏过大则丢失街道级信息。用GeoPandas做网格映射import geopandas as gpd from shapely.geometry import Point, Polygon import numpy as np # 生成武汉城区0.01度网格 lng_min, lng_max 113.8, 114.8 lat_min, lat_max 30.2, 30.9 grids [] for i in np.arange(lng_min, lng_max, 0.01): for j in np.arange(lat_min, lat_max, 0.01): grids.append((i, j, i0.01, j0.01)) grid_df gpd.GeoDataFrame({ grid_id: range(len(grids)), geometry: [Polygon([(a,b),(c,b),(c,d),(a,d)]) for a,b,c,d in grids] }, crsEPSG:4326) # 将行程起终点映射到网格 df_trip[geo_o] gpd.points_from_xy(df_trip[o_lng], df_trip[o_lat]) df_trip[geo_d] gpd.points_from_xy(df_trip[d_lng], df_trip[d_lat]) o_join gpd.sjoin(df_trip[[trip_id, geo_o]], grid_df, howleft) d_join gpd.sjoin(df_trip[[trip_id, geo_d]], grid_df, howleft) od_matrix pd.crosstab(o_join[grid_id], d_join[grid_id])pd.crosstab输出的是一个N×N稀疏矩阵N是网格数。武汉城区大约0.1度×0.7度的范围网格数在700个左右OD矩阵就是700×700完全可以直接操作。跨江需求分析不要只看OD矩阵绝对值要做期望流量修正——因为网格OD矩阵天然和网格内人口密度、上车点总数相关。一种常见做法是计算每个网格对间的流量占比即流量除以该网格总上车数消除规模差异后跨江OD的实际强度才可比。提示OD矩阵的结果不要直接读成热点结论。如果发现长江两岸网格间流量很高先确认是否是轮渡和地铁上下客点在出租车轨迹附近因为出租车上车点可能在地铁站出口周围流量会叠加。5. 从轨迹速度分布反推路况拥堵验证用对照实验确认挖掘结果可信5.1 把轨迹速度聚合到路段上的交叉验证热点聚类和OD分析做完如何确定结果可信拿官方路况数据对不上的情况很常见尤其是公开轨迹数据往往只有部分车辆覆盖不全。一个不依赖外部数据的验证方法用同一份轨迹数据分别挖掘出拥堵路段然后与你聚类结果里的载客热点做时间相关性对照。拥堵路段识别从轨迹速度做起。# 每辆车每5分钟的平均速度用于路况时段分析 df_trip[time_bin] df_trip[time].dt.floor(5min) speed_bin df_trip.groupby([time_bin, vehicle_id])[speed].mean().reset_index() # 按5分钟窗口聚合区域平均速度 area_speed speed_bin.groupby(time_bin)[speed].agg([mean, median, count]) area_speed[time] area_speed.index area_speed area_speed.set_index(time) # 计算相对拥堵指数当前速度/全天平均速度 area_speed[congest_idx] area_speed[median] / area_speed[mean].mean()拥堵指数小于0.8可以认为该时段出现拥堵状态。将这个拥堵指数和对应时段的载客热点簇数量做相关性分析。逻辑是晚高峰交通拥堵时路上车速慢、乘客等待时间更长但热点簇的总数不应急剧下降因为需求仍然存在。如果某一时段拥堵指数上升伴随热点簇数量显著下降说明这个聚类受噪声影响需要检查DBSCAN参数。from scipy.stats import pearsonr # 按小时聚合拥堵指数和热点簇数量 hourly_congest area_speed.resample(1h).median() hotspot_count_per_hour df_trip[df_trip[trip_start] True].groupby(df_trip[time].dt.hour).apply( lambda g: len(set(DBSCAN(eps150, min_samples20, metricprecomputed).fit_predict( pairwise_distances(g[[lng, lat]].values, metrichaversine_dist) ))) - 1 ) print(pearsonr(hourly_congest[congest_idx], hotspot_count_per_hour))相关系数绝对值大于0.5时说明拥堵时段和热点分布有同步变化聚类结果在时间维度上自洽。如果系数接近0优先检查轨迹时间戳是否有时区偏移以及载客状态是否在真实交接班时段产生了错误标记。5.2 用剪枝后的特质数据和Baseline对比验证另一种验证方法更直接不调整算法而是更换数据样本构造对照实验。把轨迹数据按车辆类型分为两类——武汉的出租车中有传统的双班制车辆和少量新能源单班车它们的载客状态切换频率不同。分别跑热点聚类对比热点中心点经纬度的偏移距离。# 假设已识别双班车与单班车分别聚类 for label in [double_shift, single_shift]: grp df_trip[df_trip[vehicle_type] label] pts grp[grp[trip_start] True][[lng, lat]].values D pairwise_distances(pts, metrichaversine_dist) labels DBSCAN(eps150, min_samples15, metricprecomputed).fit_predict(D) centers [] for lab in set(labels): if lab -1: continue centers.append(np.mean(pts[labels lab], axis0)) print(label, clusters:, len(centers))如果两类的热点中心点坐标差超过200米或者簇数量差异超过30%说明你的聚类结果受车辆类型分布影响过大热点不代表真实需求而更可能反映了抽样偏差。此时不要急着发表结论需要按车辆类型分层抽样再做一次聚类才能得到稳健结果。5.3 轨迹压缩验证恢复原轨迹的误差控制道格拉斯-普克压缩本身带有信息损失验证压缩对下游影响的一种做法是将压缩后的轨迹点还原成折线然后计算与原始轨迹的平均Hausdorff距离。这个指标若小于30米说明压缩没有改变道路属性。from scipy.spatial.distance import directed_hausdorff def hausdorff_dist(orig_pts, comp_pts): orig_arr np.array(orig_pts) comp_arr np.array(comp_pts) h1 directed_hausdorff(orig_arr, comp_arr)[0] h2 directed_hausdorff(comp_arr, orig_arr)[0] return max(h1, h2) # 取一条行程对比 sample_trip df_trip[df_trip[trip_id] 1] orig list(zip(sample_trip[lng], sample_trip[lat])) comp compressed[1] # 该行程压缩后的点序列 print(Hausdorff distance:, hausdorff_dist(orig, comp), meters)如果Hausdorff距离过大把epsilon从0.0005调整到0.0002即约20米精度。注意不要调得太小否则压缩退化失去意义也不要只验证一条轨迹抽样100条取平均和P90只有P90小于30米才算全量可接受。本文还有配套的精品资源点击获取