简介《北京市出租车GPS轨迹数据统计特征分析》论文PDF全文适合测绘、地理信息系统、交通数据分析等方向的师生及研究人员作为参考文献。内容以2015年5月11—17日北京市3万多辆出租车约3.85亿条GPS轨迹点数据为对象系统介绍了非载客段剔除、研究区范围筛选等预处理流程以及基于R语言的行驶距离、持续时长、分时段出行量、出行方向与平均速度五维统计特征分析并展示了指数分布拟合、距离衰减效应及工作日/周末日律性等结论对理解城市居民出行规律颇有帮助。资源为单份PDF文档大小2.56MB保留正式期刊论文版式图表与公式完整可直接引用或打印研读。该资料已有479人浏览学习适配轨迹数据挖掘、城市规划与交通预测等研究方向。1. 出租车GPS轨迹数据统计特征分析没做清洗直接统计结果全是噪声这份资源是一篇2020年发表在《测绘与空间地理信息》上的论文《北京市出租车GPS轨迹数据统计特征分析》研究对象是2015年5月11日至17日北京市3万多辆出租车产生的约3.85亿条GPS轨迹点数据。如果你的工作涉及轨迹数据挖掘、交通流分析或城市出行规律研究这篇论文最大的价值不是结论本身而是整套从原始轨迹到统计结论的处理链路非载客段剔除、区域过滤、速度异常清洗、特征量计算再到R语言的幂函数拟合与分布验证。它特别适合两类人——刚接触GPS轨迹数据、不知道从哪一步下手的初学者以及已经踩过“数据没洗干净就跑模型”坑的从业者。论文能用一周约7500万个有效载客点拟合出出行距离和持续时长服从的指数分布规律并从中读出距离衰减效应与工作日、周末的日律性差异。这正是轨迹数据统计分析的正确打开方式先清洗再计算最后用分布说话。2. 原始数据清洗从3.85亿点到7500万点三段过滤缺一不可2.1 字段理解与空载剔除VFLAG字段是第一个分水岭拿到出租车GPS数据第一件事不是写统计脚本而是搞清楚每条记录里有哪些字段、每个字段的业务含义是什么。论文的表1给出了四个关键字段SUID是每辆车的唯一标识码UTC是世界标准时间LAT和LON是经纬度而VFLAG描述空重车状态0表示没有载客。这个字段是整条分析链路能否成立的前提。原始数据里混着大量空车状态的数据点这些点在寻客或返回途中产生不携带乘客出行特征直接参与统计会把出行距离和时间的分布彻底带偏。论文的处理方式是直接根据VFLAG字段剔除非载客段。实际操作中需要注意不同数据源的VFLAG定义可能不一致有的用0/1有的用True/False有的用中文“重车/空车”第一步先做字段映射统一成二值变量。# 读取原始CSV假设字段名与论文表1一致 taxi_raw - read.csv(beijing_taxi_20150511_17.csv, stringsAsFactors FALSE) # VFLAG字段统一映射0空载1载客 # 有些源数据可能用TRUE/FALSE或“重车/空车”先转成标准0/1 taxi_raw$is_occupied - ifelse(taxi_raw$VFLAG 1 | taxi_raw$VFLAG 载客 | taxi_raw$VFLAG TRUE, 1, 0) # 剔除空载点只保留载客段 taxi_occupied - taxi_raw[taxi_raw$is_occupied 1, ] # 输出清洗前后行数对比确认剔除比例 cat(原始点数:, nrow(taxi_raw), \n) cat(载客点数:, nrow(taxi_occupied), \n) cat(空载占比:, round(1 - nrow(taxi_occupied) / nrow(taxi_raw), 4) * 100, %\n)这里的关键不是ifelse判断本身而是意识到空载数据对后续统计的污染程度。论文最终从3.85亿点清洗到7500万点剔除比例约80%大头就是空载段数据。如果你拿到的数据没有VFLAG字段也有替代方案通过轨迹段内相邻点的时间间隔和状态切换特征推断载客状态但精度会差很多优先找带空重车标识的数据源。2.2 空间范围裁剪范围设错会把长途订单全删光数据清洗的第二道关卡是空间范围。研究区是北京市论文给出的经纬度范围是经度115.7°E—117.4°E、纬度39.4°N—41.6°N。但这里有个细节值得注意直接按北京市范围裁剪会把跨省长途订单全部误删。论文的处理方式是考虑到部分车辆执行长途任务起点或终点在北京市内于是将经纬度范围扩大到河北省范围即经度113°E—120°E纬度范围原文写的是“113°—120°”结合上下文能判断这明显是OCR识别错误纬度不可能到120°N合理的范围应该是39°N—42°N量级。裁剪逻辑也有讲究不是删掉越界的那一个点而是如果当前轨迹段内至少有一个点超出范围就剔除整条轨迹段。这样保证了轨迹段的完整性不会出现一条订单轨迹被拦腰切断的情况。# 定义研究区扩展范围经度113-120E纬度39-42N lon_min - 113; lon_max - 120 lat_min - 39; lat_max - 42 # 标记越界点 taxi_occupied$out_of_bound - ifelse( taxi_occupied$LON lon_min | taxi_occupied$LON lon_max | taxi_occupied$LAT lat_min | taxi_occupied$LAT lat_max, 1, 0 ) # 按轨迹段分组只要段内有一个点越界整段剔除 library(dplyr) bad_segments - taxi_occupied %% group_by(SUID, trip_id) %% # trip_id需要根据时间差自行切分 summarise(has_outlier max(out_of_bound)) %% filter(has_outlier 1) taxi_clean - taxi_occupied %% anti_join(bad_segments, by c(SUID, trip_id))这段代码里出现了一个trip_id原始数据里不一定有现成的轨迹段编号需要根据SUID和时间戳自行切分。常见的做法是同一辆车相邻两个点时间间隔超过某个阈值比如300秒就认为是一条新的轨迹段。这一步在特征量计算之前必须完成否则后续的距离、时长统计全都算不对。2.3 速度异常剔除180km/h阈值不是拍脑袋定的GPS定位受卫星信号、隧道遮挡、设备故障等因素影响偶尔会冒出一些离谱的点位漂移。论文的处理方法很直接通过相邻轨迹点的行驶里程和时间计算出瞬时速度循环查找超过速度阈值的轨迹段并剔除该轨迹段终点。阈值设定为180km/h依据是车辆的行驶速度极限。为什么是180km/h而不是120km/h城市道路限速一般是60—80km/h但出租车在高速公路上也可能跑到120km/h左右如果阈值设到120会把正常高速行驶的订单误删。180km/h留出了足够余量只杀那些明显不可能的速度尖峰。需要注意的是这里剔除的是超速轨迹段的终点而不是整条轨迹段——因为超速点往往是GPS漂移造成的异常点把该点去掉后轨迹段的其余数据仍然有效。# 按轨迹段内相邻点计算速度用大圆距离/时间差 # 假设已有函数 haversine_dist(lon1, lat1, lon2, lat2) 返回米 speed_limit - 180 # km/h taxi_clean - taxi_clean %% arrange(SUID, trip_id, UTC) %% group_by(SUID, trip_id) %% mutate( prev_lon lag(LON), prev_lat lag(LAT), prev_time lag(UTC), seg_dist_m mapply(haversine_dist, prev_lon, prev_lat, LON, LAT), seg_time_h as.numeric(difftime(UTC, prev_time, units hours)), seg_speed_kmh ifelse(seg_time_h 0, seg_dist_m / 1000 / seg_time_h, NA) ) %% filter(is.na(seg_speed_kmh) | seg_speed_kmh speed_limit)写完这段要提醒一句直接filter删除超速点会破坏轨迹点的时间连续性导致后续计算相邻点距离时前后两个点的时间间隔拉大速度被低估。更稳妥的做法是先标记超速点检查超速点前后的速度是否连续如果只有单个点漂移删除后可以插值补点或者用前后两个正常点的直线距离代替漂移点。论文没有展开这一步但实际复现时这是最常见的翻车点。2.4 特征量计算大圆距离、速度与方向角数据清洗完成后进入特征量计算阶段。论文计算了四类特征时间、距离、速度和方向。时间方面把UTC时间转换为北京时间计算连续两点之间的持续时长。距离计算只考虑轨迹本身的几何特征不考虑路网信息直接计算轨迹点的大圆距离整条轨迹的总里程由所有相邻点的大圆距离累计而来。大圆距离公式用的是球面几何的标准形式论文式(1)里地球半径取6371.393km。# 大圆距离计算Haversine公式 haversine_dist - function(lon1, lat1, lon2, lat2) { R - 6371.393 # 地球半径单位km dlat - (lat2 - lat1) * pi / 180 dlon - (lon2 - lon1) * pi / 180 a - sin(dlat / 2)^2 cos(lat1 * pi / 180) * cos(lat2 * pi / 180) * sin(dlon / 2)^2 c - 2 * asin(min(1, sqrt(a))) return(R * c * 1000) # 返回米 }方向计算用的是线性方向平均值Linear Directional Mean把每条轨迹段看成由相邻两点构成的有向线段每段的方向角θi以正东方向为基准逆时针旋转得到。然后通过四个象限条件对合成方向角做调整得到最终的方向统计值。这一步的R语言实现可以直接调用circular包也可以手写公式。# 线性方向平均值计算弧度 theta - atan2(sum(sin(angle_vector)), sum(cos(angle_vector))) # 象限调整 if (sum(sin(angle_vector)) 0 sum(cos(angle_vector)) 0) { result - theta } else if (sum(sin(angle_vector)) 0 sum(cos(angle_vector)) 0) { result - 180 - theta } else if (sum(sin(angle_vector)) 0 sum(cos(angle_vector)) 0) { result - 180 theta } else { result - 360 - theta }速度计算更简单相邻点距离除以时间差即可。原始数据的采样间隔绝大部分为60秒这个间隔对城市出租车轨迹来说不算密意味着两个点之间的路径可能已经拐了好几个弯大圆距离会比实际行驶路径短。这是论文方法的一个固有边界它统计的是轨迹点间的直线移动特征不是道路实际行驶特征。理解这一点后面解释速度分布和工作日/周末差异时才能把握住粒度。3. 四维统计特征距离、时间、方向、速度怎么分别建模3.1 出行距离3km为界两段幂函数拟合载客段出行距离统计的第一步是划定分析范围。论文对全部载客段距离做分布统计发现0—50km范围覆盖了99%以上的数据因此聚焦这个区间以1km为间隔统计轨迹段数量。随后发现工作日与周末在峰值和出行总量上差异显著果断将工作日和周末分开讨论。工作日载客段出行距离在3km左右达到峰值随后逐渐减少呈现典型的距离衰减效应。周末的分布形状类似但出行量明显下降工作日峰值距离当天约有38000条轨迹段5月16日周六只有30000条左右5月17日周日进一步降到23000条左右。工作日出行量最高、周六次之、周日最低这个排序本身就构成一条出行意愿的量化结论。分段拟合的思路值得细看。论文以3km为界恰好对应北京市出租车起步价距离将出行距离分为1—3km短距离和3—50km长距离两段。短距离用式(7)p(d) a·d^(-β)·exp(-γd)拟合长距离用式(8)p(d) (d d0)^β·exp(-αd)拟合。这种分段处理的原因在于短距离和长距离的出行决策机制不同短距离受起步价、步行替代效应影响长距离更接近城市空间结构的重力模型。短距离拟合发现一个稳定结论一周内所有日期的拟合曲线都光滑且相似峰值出现在2.2km左右1.7—2.7km区间内的载客段数量占比约50%1.5—2.7km区间占比约80%。这说明3km以内的短途出行行为高度稳定受工作日/周末影响极小。长距离拟合则出现了重尾分布符合帕累托法则表4给出了7天内3—15km出行量占比稳定在78.3%—80.0%之间二八定律在城市出租车出行距离上得到了干净利落的验证。3.2 出行时间短时间幂律与长时间重尾时间维度的分析分了三层。第一层是分时段出行量统计——以小时为单位统计一周内不同时段的载客段数量得到的时间分布折线图能直观看出城市居民的生活节奏。核心发现是工作日出行趋势高度相似周末也有自己的相似性一周七天都在凌晨4点左右出现出行量最低值工作日有明显的波峰波谷周末没有明显峰值凌晨4点后上升到8—9点基本维持在17000左右波动。第二层是出行持续时间分布。统计发现乘客出行时间在10分钟左右达到峰值99%以上的出行持续时长在2小时以内因此分析范围锁定在0—120分钟。同样按10分钟为界分成短时间1—10min和长时间10—120min。短时间用式(9)p(t) a·t^β拟合长时间用式(10)p(t) a·t^(-β)·e^(-γt)拟合。短时间出行的统计特征是0—5min和5—10min各占50%左右整体随时间的增加出行量出现上升趋势。长时间出行则再次出现重尾分布10—34min以内的出行量占比稳定在76.5%—83.6%之间同样符合二八定律。第三层是长周期性时间分布——以一周为长周期统计每天每个小时内上车点和下车点数量。结果表明工作日普遍在9—10点、13—15点、21—22点达到峰值周末在9—10点达到峰值后基本在一定范围内波动周末出行量明显少于工作日。这组峰值时段是分析乘客出行意图的重要线索如果通勤是出租车出行的主要驱动早晚高峰的7—9点和17—19点应该出现峰值但实际峰值却落在9—10点、13—15点和21—22点。3.3 出行方向线性方向平均值揭示棋盘式路网方向分析是对轨迹段方向做分角度统计。论文采用线性方向平均值作为方向特征值统计一周七天载客段的方向分布后发现7天的方向分布高度相似但角度分布不均匀呈现明显的空间异质性。主线方向集中在0°—10°和270°—280°次线方向分布在180°—190°和90°—100°。这个结果跟北京市的道路格局高度吻合。作为棋盘式路网的代表城市北京主干道多为东西朝向或南北朝向东西朝向占主要部分次干道也基本保持这种朝向。东城区、西城区道路基本沿正东、正西、正南、正北四个方向朝阳区、海淀区、丰台区、石景山区的道路方向也集中在这四个方向。方向分布的主峰值对应东西向道路次峰值对应南北向道路而0°—10°和90°—100°之间的小角度偏移反映的是部分道路在正方向基础上存在的小范围扩散趋势。这里有一个容易被忽视的分析价值轨迹方向分布可以作为城市路网方向的侧面验证手段。如果某个城市的轨迹方向分布与已知路网方向明显矛盾说明轨迹数据可能存在系统性偏移比如GPS设备安装角度偏差或坐标转换参数错误。3.4 平均速度8点和18点的谷值为什么比出行量高峰更值得关注平均速度的分析直接反映城市道路拥堵状态。论文按小时统计每天的平均速度曲线后发现一天内平均速度从凌晨5—6点开始快速下降8点左右达到第一个谷值之后上升至12点左右形成小峰值随后逐渐下降18点左右达到全天最小值接着再迅速上升。8点和18点正是上下班高峰的典型拥堵时段。更有意思的是周末对比周末最小速度约20km/h高于工作日的约15km/h说明工作日道路普遍比周末更拥堵。这一点本身不意外但论文接下来的交叉对比才是亮点——把平均速度曲线和出行量曲线放在一起看发现出行量高峰9—10点、13—15点、21—22点与堵车高峰8点、18点并不一致。这个错位意味着北京市出租车出行的主体并不是早晚高峰通勤的上班族而是其他商务出行和事务性出行人群。如果你在做交通需求预测只看出租车订单量推断道路拥堵程度会在这个错位上栽跟头。4. 避坑与常见问题轨迹统计分析的五个实操陷阱4.1 空重车字段值不统一统计前没做映射导致载客段数量翻倍现象按VFLAG字段过滤后得到的载客段数量比预期多了近一倍分布曲线出现明显的双峰。原因数据源里VFLAG字段存在多种取值比如0/1、True/False、“重车/空车”混在一起。直接用数值过滤只保留了0和1但True/False被当成字符处理后全部保留导致空载数据混入统计。解决读取数据后先做字段去重统计table(taxi_raw$VFLAG)看所有取值统一映射成0/1二值变量再过滤。这个动作建议作为数据清洗的第一步在剔除空载之前完成。4.2 单点越界直接删点轨迹段被拦腰截断现象空间范围裁剪后统计出来的出行距离大量集中在1km以内的极短值轨迹段平均时长明显偏低。原因只删了越界的那一个点没有考虑轨迹段的完整性。一条从北京到河北的长途订单中间某几个点越界被删剩下的点被当成多条短轨迹参与统计。解决按论文的做法先逐点标记越界状态再按轨迹段分组段内只要有一个点越界就整段剔除。代价是损失一部分跨省长途订单但换来了统计口径的一致性和轨迹段的完整性。这一步必须在切分trip_id之后做否则没法按段分组。4.3 速度阈值设置过低正常高速订单被误删现象清洗后载客段数量比论文少了约15%长途订单几乎消失。原因速度阈值设置成120km/h出租车在机场高速、五环等路段行驶时瞬时速度能到110—130km/h加上60秒采样间隔下的大圆距离计算误差正常的120km/h高速行驶被误判为超速。解决按照论文的思路把阈值放宽到180km/h只剔除物理上不可能的速度尖峰。如果你的数据集采样间隔不是60秒而是更短比如10秒瞬时速度精度更高阈值可以适当收窄到150km/h左右但要先画出速度分布直方图确认尖峰的位置。4.4 直接把所有日期合在一起统计日韵律被平均淹没现象统计结果曲线平滑得看不出任何规律距离衰减效应消失工作日和周末差异也不明显。原因把一周7天的数据混在一起统计均值工作日和周末的分布差异被平均掉日律性特征被抹平。解决严格按工作日/周末分组统计必要时按星期逐日对比。论文就是发现工作日与周末在峰值和出行量上有较大差异后才把两者分开讨论后续所有的拟合和占比分析都分了两组。建议在探索阶段先把7天的分布曲线画在一张图上颜色区分星期几比直接算均值信息量大得多。4.5 出行量高峰≠拥堵高峰用订单量推断路况会翻车现象用出租车载客段数量推断道路拥堵时段发现推断结果和实际路况数据对不上。原因论文的交叉对比已经证明出行量高峰9—10点、13—15点、21—22点与堵车高峰8点、18点并不一致。出租车出行的主体不是通勤人群用出租车订单量表征全路网的交通需求存在系统性偏差。解决如果目标是分析城市拥堵应该直接使用速度分布数据如果只有出租车轨迹数据必须明确其代表的是出租车出行需求而非全量交通需求。在论文结论的解读和复用场景描述中始终把“出租车乘客”作为分析主体不要偷换概念。5. 从统计到可复现R语言幂函数拟合脚本与结果验证前四章把论文的完整分析链路拆完了这一章落地一套可以直接改参数复现的R语言拟合脚本。论文的核心结论——短距离出行指数衰减、长距离重尾分布、出行时间二八定律——全部由幂函数拟合得出拟合质量直接决定结论可信度。下表汇总了论文四个拟合结果的关键参数做验证时直接对着核对拟合对象拟合公式工作日参数周末参数短距离(1-3km)p(d) a·d^(-β)·exp(-γd)a0.076, β-1.749, γ0.778a0.074, β-1.655, γ0.737长距离(3-50km)p(d) (dd0)^β·exp(-αd)d04.548, β0.814, α0.085d03.723, β0.849, α0.081短时间(1-10min)p(t) a·t^βa0.028, β0.825a0.029, β0.797长时间(10-120min)p(t) a·t^(-β)·e^(-γt)a0.138, β-0.257, γ0.073a0.173, β-0.234, γ0.079用R语言做幂函数拟合时最直接的方案是对公式两侧取对数转换成线性模型后用lm()求解。以短距离拟合p(d) a·d^(-β)·exp(-γd)为例两侧取自然对数后得到ln(p) ln(a) - β·ln(d) - γ·d变成一个二元线性回归问题。注意这里β在原文公式里是负指数取对数后回归系数是负数写代码时符号别搞反。# 短距离拟合ln(p) ln(a) - β*ln(d) - γ*d # data: 数据框, 列名为 dist_km 和 count fit_short - lm(log(count) ~ log(dist_km) dist_km, data subset(data, dist_km 1 dist_km 3)) # 提取参数并还原 a - exp(coef(fit_short)[1]) beta - -coef(fit_short)[2] # 注意符号回归系数是-β取负还原 gamma - -coef(fit_short)[3] cat(sprintf(a%.3f, β%.3f, γ%.3f\n, a, beta, gamma))这段代码的核心是符号处理。取对数后原公式的-β·ln(d)项对应的回归系数是负值还原参数时要取负γ同理。如果数据量很小直接用nls()做非线性最小二乘拟合也可以但初值选择需要谨慎建议先用线性化方法算一组初值再用nls()精调比直接盲试初值稳定得多。结果验证环节论文用了两个可复现的占比指标不需要重新拟合也能快速检验你的清洗和统计流程是否正确出行距离在3—15km以内的载客段数量占比应稳定在78%—80%之间出行时间在10—34min以内的载客段数量占比应稳定在76.5%—83.6%之间。这两个占比指标对数据清洗质量非常敏感如果空载数据混入、轨迹段切分不当或者速度异常点没清干净占比立刻偏离区间。我在复现时先用这两个指标做自检占比落在区间内才继续做拟合省了很多排查时间。# 占比验证距离3-15km和时长10-34min dist_ratio - nrow(subset(trips, dist_km 3 dist_km 15)) / nrow(trips) time_ratio - nrow(subset(trips, duration_min 10 duration_min 34)) / nrow(trips) cat(sprintf(3-15km占比: %.1f%%\n, dist_ratio * 100)) cat(sprintf(10-34min占比: %.1f%%\n, time_ratio * 100))最后说一个我复现时踩过的细节原始数据的采样间隔不全是60秒部分车辆存在30秒或120秒的采样频率直接在整份数据上统一按60秒间隔计算速度会让部分轨迹段的速度系统性偏高或偏低。我的处理方式是按车辆分组后先检查采样间隔分布对非60秒采样的轨迹段单独标记计算速度时使用各自的实际时间差而非固定间隔。从那以后我每次拿轨迹数据都强制把“检查采样间隔分布”写进清洗流程的第一步这个习惯帮我避开了很多后续拟合阶段的诡异偏差。希望帮到你。本文还有配套的精品资源点击获取
