上个月帮同事处理一批货车GPS轨迹需求是把5000多条行驶记录按行车模式归类。我打开Matlab第一反应是轨迹聚类嘛直接调kmeans不就行了结果报错——因为轨迹是变长的二维坐标序列kmeans吃的是定长特征矩阵。把坐标强行补零对齐后聚类结果又是一团浆糊。后来我把思路拆成先特征化再Kmeans两步链路才真正跑通。这篇文章就把这条完整链路写清楚轨迹预处理、特征提取、Kmeans聚类、结果评估Matlab代码全部贴出粘贴就能跑。适合正在做GPS轨迹分类、物流路径分析、运动模式挖掘或者课程作业里正好卡在轨迹聚类上的朋友参考。不用一次看懂所有原理先跑一遍再回头对照文字效率高很多。1. 先想清楚Kmeans到底能聚类什么1.1 轨迹数据和普通表格数据的三个差别新手最容易犯的错就是把轨迹数据当成普通数据表直接丢给算法。普通数据是一行一个样本一列一个特征的定长矩阵而轨迹数据本质上是变长序列两者有本质区别。第一个差别样本不是向量而是序列。一条轨迹由一串带顺序的GPS点组成比如[ (x1,y1), (x2,y2), ..., (xm,ym) ]m就是这条轨迹的点数。有的轨迹300个点有的轨迹只有15个点Kmeans要求输入是n×d的数值矩阵每一行是一个样本d必须固定。轨迹长度不一致矩阵都拼不到一起。第二个差别轨迹的特征藏在关系里不是某个坐标值里。一辆车是否在高速行驶取决于相邻点之间的距离和时间差不取决于它当时在路网哪个位置。如果直接把坐标当成特征聚类结果基本等于按空间位置划堆完全没有运动模式的语义。第三个差别轨迹之间的距离不能用普通欧氏距离衡量。两条长度一样的轨迹就算第k个点的坐标完全相同也可能走了完全相反的路顺序颠倒。轨迹相似度必须考虑形状、走向、时序这些都不是简单向量空间能表达的。想明白这三点就知道基于Kmeans的轨迹聚类绝不是把轨迹塞进kmeans这么简单而要想办法把轨迹转换成Kmeans能理解的形式。1.2 三条落地路线各有各的适用场景把轨迹喂给Kmeans业内常见的有三条路线我列个表对比一下路线核心做法优点缺点适用场景A. 特征向量化 Kmeans每条轨迹提取固定维度的特征速度、长度、弯曲度等组成矩阵后直接调kmeans快、稳定、可解释性强特征靠人工设计可能丢掉部分形态信息想知道轨迹分成哪几类运动模式B. 距离矩阵 Kmedoids先算两两轨迹的Hausdorff/DTW距离再用Kmedoids聚类不需要特征设计能刻画形状相似性距离矩阵O(n²)计算贵样本一大就跑不动轨迹形态差异明显样本量几百条以内C. 轨迹点聚类 投票把所有GPS点拿出来做Kmeans看每条轨迹落在哪个簇的点最多能直接发现空间热点区域本质是点聚类长轨迹横跨多个簇时会乱共享单车潮汐点、打卡热点挖掘这篇文章的主代码走路线A因为它是Kmeans这个名字最正统的用法也最好上手。路线B在第6章给出升级思路路线C跟轨迹聚类严格说不是一回事这里就不展开了。2. 轨迹预处理聚类好坏七成在数据清洗2.1 重采样把不等长、不等间隔变成统一节奏预处理这一步很多人会跳过但实际项目中它往往决定聚类成败。先说重采样。轨迹数据经常来自不同设备有的GPS上报频率是1秒一次有的是5秒一次还有的是行驶时上报、停车不上报。这些轨迹的采样间隔五花八门如果直接计算相邻点距离当速度等于拿不同尺度的数据硬比。解决办法是按累计路程重采样到固定点数function trajResampled resampleTrajectory(traj, N) % 按累计距离把一条轨迹重采样到N个点路径长度信息保留 dxy diff(traj); segDist sqrt(sum(dxy.^2, 2)); cumDist [0; cumsum(segDist)]; totalDist cumDist(end); if totalDist 1e-6 % 轨迹完全静止直接按原始点复制避免除零 idx round(linspace(1, size(traj,1), N)); trajResampled traj(idx, :); return; end queryDist linspace(0, totalDist, N); trajResampled [interp1(cumDist, traj(:,1), queryDist, linear), ... interp1(cumDist, traj(:,2), queryDist, linear)]; endinterp1是Matlab自带的一维插值函数按累计距离在原始轨迹上取N个点。cumulative distance的好处是不管原来的点疏密如何重采样后轨迹在路程维度上是均匀的速度特征就有了可比性。调用方式trajCell{1} resampleTrajectory(trajCell{1}, 50);2.2 GPS漂移点剔除别让乱飘的点统治聚类结果GPS数据最常见的脏点是漂移点——设备在隧道、高楼旁边突然跳出去几十米甚至几百米。这种点如果不处理会让轨迹的最大速度特征变成天文数字弯曲度也被拉高聚类结果直接被少数异常点带偏。判断漂移点其实就一条逻辑相邻两个GPS点之间的距离除以时间间隔超过物理运动极限。比如货车的物理极限可以认为不超过120km/h也就是33.3m/s如果某一段的瞬时速度远超这个值基本可以断定是漂移。function trajClean removeDriftPoints(traj, speedThreshold) % speedThreshold 单位是单位采样间隔的欧氏距离按实际数据标定 dxy diff(traj); segDist sqrt(sum(dxy.^2, 2)); keepIdx true(size(traj,1), 1); % 前一个点 keepIdx(1:end-1) segDist speedThreshold; % 后一个点也要满足否则从中间切断 keepIdx(2:end) keepIdx(2:end) (segDist speedThreshold); trajClean traj(keepIdx, :); end注意我同时检查了段的前点和后点因为一段距离异常通常是前后两个点里至少一个是坏的只去掉一个还可能在轨迹里留下跳变。阈值怎么定如果采样间隔是1秒货车阈值建议设40米如果是步行数据10米就很宽松了。2.3 经纬度坐标先做投影别直接算欧氏距离手里轨迹如果用的是经纬度还有个隐藏坑1度纬度的地面距离约111km1度经度的地面距离随纬度变化在赤道附近也是111km但在北纬40°大约只有85km。直接用经纬度做欧氏距离x方向和y方向的比例尺不一样会扭曲轨迹形状。最简单的处理办法如果研究区域不大用cos(mean(lat))对经度做个缩放lat0 mean(allLat); lonScaled allLon * cos(lat0 * pi / 180);更正规的做法是投影到UTM坐标Matlab的map工具箱里有现成函数geodetic2utm或者用projfwd。对聚类项目来说只要保证x、y单位一致且是米制精度完全够用。3. 完整Matlab实现从模拟数据到聚类可视化3.1 构造三类有辨识度的模拟轨迹为了让你能完整复现链路我先生成三类轨迹来模拟真实场景直线快速移动车辆、弯曲慢速移动行人、原地徘徊静止停留。每类20条共60条。这样有个好处——我们自己知道真实标签后面可以验证聚类效果。clc; clear; close all; rng(42); % 固定随机种子确保结果可复现 numTraj 60; trajCell cell(numTraj, 1); trueLabel zeros(numTraj, 1); % --- 第1类直线快速轨迹20条轨迹总长约80~120米20个采样点 --- for k 1:20 startPt [5, 5] rand(1,2) * 30; angle rand * 2 * pi; dist 80 rand * 40; n 20; t linspace(0, 1, n); x startPt(1) dist * cos(angle) * t randn(n,1) * 0.4; y startPt(2) dist * sin(angle) * t randn(n,1) * 0.4; trajCell{k} [x, y]; trueLabel(k) 1; end % --- 第2类弯曲慢速轨迹20条路径带正弦摆动40个采样点 --- for k 21:40 startPt [5, 5] rand(1,2) * 30; angle rand * 2 * pi; dist 15 rand * 10; n 40; t linspace(0, 1, n); x startPt(1) dist * cos(angle) * t 3 * sin(4*pi*t) randn(n,1) * 0.3; y startPt(2) dist * sin(angle) * t 3 * cos(4*pi*t) randn(n,1) * 0.3; trajCell{k} [x, y]; trueLabel(k) 2; end % --- 第3类原地徘徊轨迹20条随机游走但整体位移很小 --- for k 41:60 center [5, 5] rand(1,2) * 30; n 30; dxy randn(n,2) * 0.3; traj cumsum(dxy) - dxy(1,:) center; trajCell{k} traj; trueLabel(k) 3; end这段代码模拟了三种典型的运动模式。说明一点为了演示方便这里假设所有轨迹采样间隔相同默认1秒所以相邻点距离可以直接代表速度真实数据如果采样间隔不一致必须先用2.1节的重采样处理。3.2 轨迹特征提取把一条轨迹浓缩成一个向量这是整条链路的胜负手。特征选得好Kmeans就是锦上添花选得不好调什么参数都没用。我常用的轨迹特征如下表序号特征变量物理含义对聚类的作用1起点坐标startPt轨迹出发位置区分空间起点不同的运动2终点坐标endPt轨迹到达位置区分空间终点不同的运动3平均步长meanSpeed相邻点平均距离采样间隔一致时等价于平均速度区分快慢运动4最大步长maxSpeed相邻点最大距离区分稳定运动与突发加速5步长标准差speedStd速度波动程度区分匀速与变速运动6总路径长度totalDist轨迹累计长度区分长距离与短距离7直线位移directDist起点到终点的直线距离区分位移大与小的轨迹8弯曲度sinuosity总路径长度 / 直线位移区分直线轨迹与蜿蜒轨迹9累计转向角turnAngle所有轨迹点方向变化的绝对值之和区分直线、缓弯、乱绕对应函数function feat extractTrajFeatures(traj) % 输入m×2的轨迹坐标 % 输出1×11的行向量特征 dxy diff(traj); segDist sqrt(sum(dxy.^2, 2)); totalDist sum(segDist); directDist norm(traj(end,:) - traj(1,:)); meanSpeed totalDist / (size(traj,1) - 1); maxSpeed max(segDist); speedStd std(segDist); angles atan2(dxy(:,2), dxy(:,1)); turnAngle sum(abs(diff(angles))); sinuosity totalDist / (directDist eps); startPt traj(1,:); endPt traj(end,:); feat [startPt, endPt, meanSpeed, maxSpeed, speedStd, ... totalDist, directDist, sinuosity, turnAngle]; end为什么选这些特征因为不同运动模式在速度分布、弯曲程度、方向变化上差异足够大。举例直线快速轨迹的弯曲度接近1累计转向角很小慢速弯曲轨迹的弯曲度在1.5~3之间转向角明显更大原地徘徊轨迹的位移很小但方向乱变转向角容易被噪声抬得很高。三者一进特征空间Kmeans能轻松切分。3.3 特征标准化与Kmeans聚类特征提取完之后大部分人的第一反应是直接调kmeans。这里有个必须做的动作标准化。比如我的模拟数据里坐标值在0~50之间而累计转向角可能在0~10之间速度特征在0~8之间。如果不标准化Kmeans的欧氏距离几乎完全被坐标值主导运动模式的特征根本排不上号。我自己手写了标准化顺便处理了零方差列的坑% 组装特征矩阵 featMat zeros(numTraj, 11); for i 1:numTraj featMat(i,:) extractTrajFeatures(trajCell{i}); end % 手写标准化比zscore稳方差为0的列不除零 mu mean(featMat); sigma std(featMat); sigma(sigma 0) 1; % 防止零方差列产生NaN featNorm (featMat - mu) ./ sigma; % Kmeans聚类 K 3; [idx, C, sumd] kmeans(featNorm, K, Replicates, 10);kmeans是Statistics and Machine Learning Toolbox里的函数。Replicates参数表示重复10次、取组内距离平方和最小的那一次作为最终结果能有效缓解Kmeans对初始中心敏感的问题。默认使用的是kmeans初始化比纯随机初始化稳得多。3.4 可视化把聚类结果画回轨迹上聚类完必须画图否则无法直观判断效果。我的做法是按簇标签给每条轨迹上色figure; hold on; box on; colors lines(K); for i 1:numTraj plot(trajCell{i}(:,1), trajCell{i}(:,2), Color, colors(idx(i),:), LineWidth, 1.2); end title(Kmeans轨迹聚类结果按簇着色); xlabel(X / m); ylabel(Y / m);模拟数据跑出来的效果会非常理想三堆轨迹颜色分明基本不存在混色。这时可以用交叉表对比真实标签和聚类标签T crosstab(idx, trueLabel); disp(聚类标签 vs 真实标签); disp(T);crosstab输出的矩阵第i行第j列表示聚类为第i类、真实为第j类的轨迹条数。如果每一行只有一个大数说明聚类结果和真实类别几乎一一对应。需要注意Kmeans给的簇编号1、2、3是随机的可能聚类第1类对应真实第3类这没关系看对角线重排后的趋势即可。4. K值怎么定肘部法则、轮廓系数与可复现性4.1 肘部法则看组内距离平方和的拐点聚类之前第一个问题就是K选几。拍脑袋定K3在演示数据里没问题真实数据绝不能用猜的。最经典的方法是肘部法则对不同的K值分别计算组内距离平方和SSE看曲线在哪里出现肘部。Klist 2:8; SSE zeros(size(Klist)); for j 1:length(Klist) [~, ~, sumd] kmeans(featNorm, Klist(j), Replicates, 10); SSE(j) sum(sumd); end figure; plot(Klist, SSE, o-, LineWidth, 1.5); xlabel(K); ylabel(组内距离平方和 SSE); title(肘部法则选择K);sumd是kmeans返回的每个簇内点到中心的距离平方和向量sum(sumd)就是整个聚类的组内平方和。K越大SSE自然越小但下降幅度会越来越平缓。出现骤降转为缓降的拐点就是合理的K。真实数据的SSE曲线往往没有模拟数据那么漂亮的拐点会平滑下降。这时候不能只靠肘部法则要配合轮廓系数一起判断。4.2 轮廓系数量化簇内紧密、簇间分离轮廓系数是Matlab一句话能跑的评估工具s silhouette(featNorm, idx); mean_s mean(s); fprintf(平均轮廓系数%.3f\n, mean_s);轮廓系数的取值范围是-1到1。大于0说明该样本和自己的簇更近、离其他簇更远负数说明样本可能分错簇了。我做了很多项目之后的体感模拟数据均值在0.7以上很轻松真实轨迹数据能到0.4~0.6已经算不错的聚类低于0.3就说明特征选得不好或K值不对需要回头调整。silhouette会直接画出一张轮廓图横轴是轮廓系数一条条横线是样本。如果出现大量负值横线这些样本就是嫌疑点有必要单独拎出来看原轨迹到底长什么样。4.3 随机种子与Replicates为什么两次聚类结果不一样Kmeans自带随机性尤其初始化中心随机。同一份数据跑两次簇编号可能完全对调第一次的簇1变成第二次的簇2。样本的分组其实一致但如果你要出图、写报告、对比实验不固定随机种子会痛苦到崩溃。解决方案就一行命令rng(42);放在kmeans调用之前。Replicates10已经让结果稳定不少但固定随机种子才是论文可复现的正道。这个项目里我先后踩过几次随机性的坑后来凡是发布结果脚本第一行永远是rng固定值。5. 实测踩坑记录这些问题我不希望你再遇到5.1 特征不归一化位置特征把运动特征压死第一次拿同事的真实数据跑聚类出来的三堆轨迹勉强能用但仔细看分类依据几乎全是起点在城西还是城东运动模式完全没体现。我查了半天最后发现就是归一化问题那份数据的坐标范围是0~100km而速度特征是0~25m/s转向角是0~10弧度坐标的方差比其他特征大了几个数量级欧氏距离被位置完全主导。所以标准化不是可有可无的步骤。手写标准化那几行代码是我在这个项目里最值的20行。5.2 零方差列导致NaN聚类结果变成一团NaN另一个真实数据里的坑某批轨迹大部分是原地不动的直接位移这一列特征几乎全为0标准差也是0。用Matlab自带的zscore一算0除以0直接出NaNKmeans拿到NaN后聚类结果极度不稳定每跑一次都不一样。后来我换成手写标准化sigma(sigma0)1直接跳过零方差列这个问题彻底消失。这个细节在官方文档里根本不会告诉你。5.3 停留段把弯曲度和转向角带偏物流轨迹里有个非常普遍的噪点源车在服务区停10分钟GPS接收器一直轻微飘动点与点之间距离在1米以内徘徊。这段轨迹如果不处理总路径长度被噪声叠加得虚高方向变化特征更是被抖动抬到爆炸最后一条停着不动的轨迹反而被聚类成乱走的一类。处理办法很简单预处理阶段用速度阈值找出低速停留段把连续低速点合并成一个点或直接删除。判断标准可以定为连续5个以上相邻点速度都低于0.3m/s就算停留段。5.4 轨迹点数差异大特征计算口径不一致有的轨迹有1000个点有的只有15个点。少的那条可能点密度很低中间隔了几百米才有一个GPS点算出来的平均步长虚高跟高频采样轨迹的平均步长完全不是一个物理口径。如果不做重采样就提取特征速度特征基本失真。这也是我把重采样放在预处理第2位的原因。遇到结构性差异大的轨迹数据先统一到同一个点数N再谈特征。5.5 没有时间戳的轨迹速度只能叫空间步长真实业务数据经常只有经纬度和点号没有时间字段。没有时间就没法算真实速度只能算相邻点空间距离。这时候聚类结论只能表达空间形态相似不能表达运动速度相似写报告时务必注明否则业务方拿着结果当速度分类用容易出问题。能补时间就尽量补结合路网、红绿灯位置估算不能补就把速度类特征从特征表里去掉只用形状类特征。6. 进阶路线当轨迹真的没有均值时走向Kmedoids6.1 为什么说轨迹没有平均值Kmeans的均值是把一个簇内所有样本逐维求平均。但两条轨迹逐点相加再平均得到的是什么如果两条轨迹长度不一样均值轨迹在几何上没有任何意义就算长度一样第k个点在两条轨迹上也不一定对应同一个物理时刻或同一个路网位置。这种平均轨迹没法代表任何一类轨迹的典型形态。这时候更稳的做法是Kmedoids让簇中心必须是簇内一条真实存在的轨迹而不是计算出来的虚像。Kmedoids聚类的核心是距离矩阵所以要先定义轨迹之间的距离。6.2 Hausdorff距离最简单的轨迹形状距离Hausdorff距离的思路很直观两条轨迹之间的最近距离最大值衡量的是其中一条轨迹上离另一条轨迹最远的点离得有多远。Matlab实现很简洁function d hausdorffDist(traj1, traj2) % 双向Hausdorff距离 D pdist2(traj1, traj2); h12 max(min(D, [], 2)); % traj1的点到traj2最近距离的最大值 h21 max(min(D, [], 1)); d max(h12, h21); end注意Hausdorff距离只比较几何形状忽略了轨迹走向的时间顺序。更严格的替代品是Fréchet距离想象人牵着狗沿两条轨迹走绳子允许伸缩绳长的最小值它保留了时间顺序语义更贴合轨迹走法是否相似的判断但计算复杂度更高。论文和比赛里直接用Hausdorff起步是比较务实的做法。6.3 手写Kmedoids骨架代码Matlab的统计工具箱也有kmedoids函数但为了让你彻底理解我给出一个不依赖工具箱调用的PAM算法骨架function [idx, medIdx] myKmedoids(D, K, maxIter) % D: n×n距离矩阵 % idx: 每个样本的簇标签 % medIdx: 每个簇的代表轨迹下标 n size(D, 1); rng(42); medIdx randsample(n, K); % 随机选K条轨迹做中心 for iter 1:maxIter % 分配每个样本选最近的medoid [~, idx] min(D(:, medIdx), [], 2); % 更新在每个簇内部找离其他成员总距离最小的点作为新中心 for j 1:K members find(idx j); if isempty(members) continue; % 空簇用上一次中心占位复杂场景可重新初始化 end subD D(members, members); [~, pos] min(sum(subD, 2)); medIdx(j) members(pos); end end end调用方式n length(trajCell); D zeros(n); for i 1:n for j i1:n D(i,j) hausdorffDist(trajCell{i}, trajCell{j}); D(j,i) D(i,j); end end [idxMed, medIdx] myKmedoids(D, 3, 30);这套方案的瓶颈是距离矩阵n条轨迹要算n(n-1)/2次Hausdorff距离每条轨迹几十个点还好上千条轨迹就很吃力。那时要么对轨迹先抽稀要么改用Mini-Batch的思路分批算。样本量在几百条以内时Kmedoids是比特征Kmeans更正宗的轨迹聚类。最后说一点个人体会代码前后加起来不到200行真正花时间的全在把轨迹讲清楚这件事上——重采样、除漂移、特征设计、归一化。你如果拿着真实轨迹数据跑先别急着调kmeans参数把数据的采样间隔、坐标单位、停留段处理干净聚类结果会自己好起来。这也是基于Kmeans的轨迹聚类这个题目里最容易被忽略却最值钱的部分。