做GPS轨迹分析的时候最常遇到的问题就是一堆密密麻麻的轨迹点怎么把它们快速分成几类共享单车的骑行轨迹、外卖配送路径、早晚高峰的通勤路线全都叠在地图上根本看不清规律。这时候轨迹聚类就派上用场了。这篇我从实际项目出发讲讲怎么用Matlab写一套基于Kmeans的轨迹聚类包括预处理、相似度计算、聚类执行和结果评估核心代码可以直接拿去改着用。适合刚接触轨迹数据、或者已经在用Matlab做数据分析但不知道怎么处理“序列型数据”的朋友。1. 轨迹聚类的应用场景与方案选型1.1 轨迹聚类到底在解决什么问题轨迹聚类本质上是在做无监督学习没有标签不知道哪条轨迹属于哪一类只能靠轨迹本身的形态特征把它们归堆。应用场景比想象中广。拿共享出行来说平台每天产生几十万条骑行轨迹运营需要知道哪些路线是高频通勤走廊、哪些是休闲骑行路线才能决定在哪里投放车辆、设置电子围栏。用人工逐条看轨迹肯定不现实聚类可以直接把相似形状的轨迹聚合在一起。再比如交通管理部门的OD分析需要找出主要的出行通道物流调度系统要识别司机常用的配送路径模式。这些都是轨迹聚类的典型落地场景。业内常用的轨迹聚类方法不少基于密度的DBSCAN、基于层次的凝聚聚类、基于模型的GMM以及这篇的主角Kmeans。Kmeans能成为最常用的基线方案不是因为它最精确而是因为它简单、快、可解释性强。尤其当你处理的轨迹已经是某种固定长度的特征向量时Kmeans几乎是性价比最高的选择。1.2 为什么选择Kmeans作为基础算法Kmeans的目标很直白把N条轨迹分到K个簇让每个样本到其簇中心的平方距离之和最小。算法流程就是经典的迭代——初始化K个中心点计算每个样本到中心点的距离把它归到最近的中心再重新计算每个簇的中心重复直到中心不再变化。这个逻辑迁移到轨迹聚类上只有一个障碍轨迹不是高维空间里的点而是点的序列。所以用Kmeans做轨迹聚类的关键不是Kmeans本身而是“如何把轨迹转化成Kmeans能算距离的数据形式”。这也是这篇博文里我最想强调的部分——很多人直接把经纬度序列塞给Kmeans结果聚得一塌糊涂其实是没想清楚距离度量问题。Kmeans的另一个优势是扩展性好。2万条轨迹每条提取100维特征在普通笔记本上跑Matlab的kmeans函数也就几秒钟的事比DBSCAN这种需要建图算邻域的算法快得多。虽然它对非凸簇效果一般、对噪声敏感、需要事先指定K值但这些问题都有成熟的处理手段后面专门用一节来聊。2. 从原始轨迹到可聚类特征预处理与相似度度量2.1 原始轨迹为什么不能直接丢给Kmeans先看一个非常典型的场景。假设你有两条轨迹都是从同一条路从A点到B点但一条是步行每秒记录一个点另一条是开车每5秒记录一个点。两点之间直线距离一样轨迹形状几乎重合但点的数量差了好几倍。如果用原始的坐标点做距离计算这两条轨迹完全没法对齐计算出来的“距离”会大得离谱。这就是轨迹数据最核心的难点长度不固定、采样频率不一致、坐标存在漂移噪声。另外还有一个更隐蔽的问题——方向性。同一条路来回走两次如果单纯比对点的位置集合它们是相似的但如果比对点的顺序它们恰恰是相反的。你是关心“形状相似”还是关心“路径方向也一致”这个决策直接决定后面的方案选型。所以轨迹聚类的第一步不是选算法而是先确定“怎样算两条轨迹相似”。这个定义定不好后面全白搭。2.2 特征提取从“序列”到“可以算距离的向量”把轨迹从“序列”变成“向量”常见的有三条路。第一条路关键点抽稀。用Douglas-Peucker算法把一条轨迹压缩成少数几个特征点比如一条1000个点的通勤轨迹抽稀后可能只剩12个关键转角点。之后每个轨迹就变成12个点组成的序列长度差异缩小再做距离计算就稳定多了。这条路的缺点是点数量依然可能不一致还需要二次对齐。第二条路网格化编码。把地图划分成规则的网格比如32×32每条轨迹经过的格子记录下来变成一组网格编号。这条轨迹就从一个坐标序列变成了一个“网格集合”既可以算集合的Jaccard相似度也可以把网格映射成0/1向量。这也是我后面代码里采用的方法优点是对采样频率不敏感速度快。第三条路统计特征。直接算起点、终点、轨迹长度、平均速度、行驶方向分布等标量指标拼成一个向量。这个方法最粗糙但处理海量数据做初筛时效率极高。我见过一些工业项目就是用“起点网格终点网格轨迹长度”三个特征做Kmeans效果居然也不错——因为很多业务场景本身就只看OD起终点和距离。2.3 距离度量Hausdorff、DTW与网格Jaccard的区别定了特征提取方式接下来就是距离度量。这里我选三种最具代表性的距离对比一下。Hausdorff距离衡量的是两条轨迹点集之间的“最大偏离程度”。形象地说就是A轨迹上每个点到B轨迹上最近点的距离取所有距离中的最大值。它适合比较几何形状对轨迹点顺序不敏感。但缺点是对离群点极其敏感——只要有一个点漂移严重整个距离就被带偏。DTW动态时间规整是处理时序序列对齐问题的经典方法。它允许两个序列在时间轴上进行非线性伸缩找到最优的对齐路径。想象两个人在不同时间内唱同一首歌音符有快有慢但旋律轮廓一致DTW就能把它们齐步走地对齐。DTW很适合轨迹序列尤其是采样频率不一致的情况缺点是比较两条长度分别为M和N的轨迹复杂度是O(M×N)数据量大时计算成本很高。网格Jaccard距离则是先把轨迹转成经过的网格集合再算集合的交并比。它的优点是极其稳定不怕采样频率不一致也不怕小幅噪声缺点是会丢失顺序信息——一条直线穿过的网格和一条蜿蜒穿过同一批网格的轨迹Jaccard距离可能是一样的。三者的选择逻辑很简单如果是离线分析、几百条轨迹、精度要求高用DTW如果是在线处理、几万条轨迹优先网格化加Jaccard。我下面的完整代码用网格化方案因为最容易跑通、最适合演示Kmeans聚类流程。距离度量复杂度对噪声敏感性对顺序敏感性适用场景HausdorffO(M×N)极敏感不敏感几何形状比较DTWO(M×N)较敏感敏感采样率不一致的序列网格JaccardO(MN)稳定不敏感海量轨迹快速聚类3. Matlab完整实现与代码拆解3.1 生成仿真轨迹数据我先模拟三类轨迹直线通勤型、环线型、L型折线型加入适量高斯噪声模拟GPS漂移。这样聚类结果容易验证你也能直接看到不同簇的可视化效果。% 仿真轨迹数据生成 % 三类轨迹直线通勤、环线、L型折线 clear; clc; rng(42); % 固定随机种子保证可复现 numTraj 60; % 总共60条轨迹 trajCell cell(1, numTraj); for i 1:numTraj classType mod(i, 3); % 0/1/2 三类循环 switch classType case 0 % 直线通勤型从(0,0)到(100,50) n 40 randi(30); % 点数随机 t linspace(0, 1, n); x 100 * t randn(size(t)) * 1.5; y 50 * t randn(size(t)) * 1.5; case 1 % 环线型近似椭圆 n 50 randi(30); t linspace(0, 2*pi, n); x 30 25 * cos(t) randn(size(t)) * 1.0; y 30 15 * sin(t) randn(size(t)) * 1.0; case 2 % L型折线先水平再垂直 n randi([30, 50]); half round(n / 2); t1 linspace(0, 1, half); t2 linspace(0, 1, n - half); x [linspace(0, 50, half), linspace(50, 80, n-half)]; y [zeros(1, half), linspace(0, 60, n-half)]; x x randn(size(x)) * 1.2; y y randn(size(y)) * 1.2; end trajCell{i} [x(:), y(:)]; end % 绘制原始轨迹未着色便于观察三类形态 figure(Name, 原始轨迹); hold on; box on; for i 1:numTraj plot(trajCell{i}(:,1), trajCell{i}(:,2), Color, [0.6 0.6 0.6]); end title(仿真轨迹数据); xlabel(X坐标); ylabel(Y坐标);这段代码里有个关键细节每类轨迹的点数都做了随机化处理而不是统一的固定长度。这样更贴近真实场景也能检验后面的特征提取方案对长度差异是否稳定。3.2 轨迹网格编码与距离矩阵计算网格编码是整个流程的地基。我把地图范围定为xLim[-5,105]、yLim[-5,70]划分成32×32的网格。每条轨迹的每个坐标点都换算成对应的网格编号然后取唯一网格集合作为这条轨迹的“指纹”。% 网格化编码函数 function codeSeq trajectoryEncoding(traj, gridSize, xLim, yLim) % 将轨迹坐标序列映射为网格编号集合 % 输入traj - Nx2矩阵 [x, y] % gridSize - 网格边长上的格子数 % xLim - [xmin, xmax] % yLim - [ymin, ymax] % 输出codeSeq - 去重后的网格编号列向量 % 避免除零错误 if xLim(2) xLim(1) || yLim(2) yLim(1) error(网格范围不能是零区间); end n size(traj, 1); codeSeq zeros(n, 1); for j 1:n % 计算该点所在的网格行列号并映射到 [1, gridSize] gx floor((traj(j,1) - xLim(1)) / (xLim(2) - xLim(1)) * gridSize) 1; gy floor((traj(j,2) - yLim(1)) / (yLim(2) - yLim(1)) * gridSize) 1; gx max(1, min(gridSize, gx)); % 边界保护 gy max(1, min(gridSize, gy)); codeSeq(j) (gy - 1) * gridSize gx; % 二维索引转一维编号 end codeSeq unique(codeSeq); % 集合化去除重复经过的格子 end有了编码函数下一步就是计算两两轨迹之间的Jaccard距离矩阵。这一步在数据量大时是性能瓶颈我会在后面专门讲优化思路。% 计算所有轨迹两两之间的Jaccard距离矩阵 gridSize 32; xLim [-5, 105]; yLim [-5, 70]; % 对所有轨迹做编码预处理 cellCodes cell(1, numTraj); for i 1:numTraj cellCodes{i} trajectoryEncoding(trajCell{i}, gridSize, xLim, yLim); end % 距离矩阵计算 D zeros(numTraj, numTraj); for i 1:numTraj for j i1:numTraj inter length(intersect(cellCodes{i}, cellCodes{j})); % 交集大小 uni length(union(cellCodes{i}, cellCodes{j})); % 并集大小 D(i,j) 1 - inter / uni; % Jaccard距离 1 - 相似度 D(j,i) D(i,j); % 距离矩阵对称 end end3.3 Kmeans聚类执行与可视化距离矩阵算出来以后每条轨迹就变成了一行“到所有其他轨迹的距离向量”维度是numTraj。这个向量可以理解为这条轨迹在“轨迹距离空间”里的坐标描述。对这个矩阵直接跑Kmeans。这里有个值得说清楚的细节Matlab自带的kmeans默认用欧氏距离我传入的就是n×n的稠密距离矩阵集群中心是在这个“距离空间”里计算的。这种做法在轨迹聚类研究里很常见相当于先用相似度定义了轨迹的内积结构再在结构空间里聚簇。% Kmeans聚类 K 3; % 我们在仿真阶段明确知道有3类先固定下来 [idx, C] kmeans(D, K, Replicates, 5, Distance, sqEuclidean, MaxIter, 200); % 可视化聚类结果不同簇使用不同颜色 figure(Name, Kmeans轨迹聚类结果); hold on; box on; colors lines(K); for i 1:numTraj tr trajCell{i}; plot(tr(:,1), tr(:,2), Color, colors(idx(i), :), LineWidth, 1.2); end title([Kmeans轨迹聚类结果K num2str(K)]); xlabel(X坐标); ylabel(Y坐标);这里我特意用了Replicates, 5——这个参数让Kmeans算法从5组不同的随机初始中心开始迭代最终返回目标函数最小的那次结果。这是对抗Kmeans初值敏感最直接的手段。还有一点Distance选的是sqEuclidean即平方欧氏距离。因为Kmeans的优化目标是簇内误差平方和用平方欧氏距离和优化目标完全一致收敛更稳。3.4 聚类质量评估轮廓系数聚类完成后不能只看图“像不像”还得有个量化指标。轮廓系数Silhouette Coefficient是应用最广的聚类评估指标之一。对每个样本它计算该样本到同簇其他样本的平均距离a以及到最近其他簇样本的平均距离b得到轮廓值(b-a)/max(a,b)。取值范围在[-1,1]越接近1代表聚类效果越好。% 轮廓系数评估 sil silhouette(D, idx); meanSil mean(sil); fprintf(平均轮廓系数: %.4f\n, meanSil); % 绘制轮廓图 figure(Name, 轮廓系数); silhouette(D, idx); title(轮廓系数图);轮廓系数还能帮你发现异常情况。比如绝大多数样本轮廓值都很高但某个簇里有个样本值是负数说明这个样本大概率被分错了簇或者它本身就是一个离群轨迹。我实际跑过这组仿真数据平均轮廓系数一般在0.75以上——对轨迹聚类来说算是不错的结果。如果你在真实数据上测出来达不到这个水平多半是之前的地图范围、网格粒度或者K值选择出了问题而不是算法本身不行。4. 参数调优与常见的坑4.1 K值怎么选肘部法则还是轮廓系数真实项目里你不会像仿真数据这样提前知道K3。确定K有两种主流方式。肘部法则的思路是对不同K值分别运行Kmeans记录簇内误差平方和SSE画成折线图。随着K增大SSE单调下降但下降速度会在某个点突然变缓形似“手肘”这个点就是推荐K值。Matlab里可以用kmeans返回的sumd属性快速计算SSEKList 1:8; SSE zeros(size(KList)); for ii 1:length(KList) [~, ~, sumd] kmeans(D, KList(ii), Replicates, 3); SSE(ii) sum(sumd); end plot(KList, SSE, -o); xlabel(K值); ylabel(SSE);轮廓系数法更直观直接对每个K值算平均轮廓系数取最高点对应的K。两种方式各有缺陷肘部法则在数据分布不清晰时手肘点不明显轮廓系数计算量大一些。我的习惯是两者结合先用肘部法则圈定一个范围再用轮廓系数在里面精挑。4.2 网格粒度与坐标范围的影响我必须要提醒一个我踩过的坑网格粒度不能盲目取。网格太粗所有轨迹都变成“几个格子”区分度下降网格太细GPS噪声被放大同一条路反复走的轨迹可能被分到不同的格子集合里聚类稳定性变差。我做下来比较稳的经验值是32×32到64×64之间。如果轨迹覆盖的地图范围是几十公里可能还需要根据轨迹的平均长度自适应调整。一个判断标准是单条轨迹平均穿过的网格数应该在15到40之间。低于10个说明网格太粗高于60个说明太细聚类容易被噪声主导。坐标范围也必须固定且合理。代码里xLim和yLim要覆盖所有轨迹点的坐标范围并且留出少量边距。如果直接用轨迹点的min/max作为边界边缘轨迹可能被压到边界上的网格里造成失真。4.3 数据顺序与随机种子问题Kmeans对初始中心敏感这是它的天性。就算同一份数据连续跑两次结果都可能不同——第一次分出来“直线型环线型L型”三个簇第二次可能把某两条直线轨迹分开、合并了另外一类。解决方式有三种固定随机种子适合调试和复现结果用Replicates跑多轮取最优适合最终决策自己实现Kmeans初始化适合教学理解。我的建议是前两者组合使用调试阶段rng(42)固定种子保证每次结果一致出最终结果时删掉固定种子、加高Replicates让算法自己找最优解。5. 真实场景中的问题排查与避坑记录5.1 问题轨迹长度差距过大导致聚类效果差真实数据里一条轨迹可能只有5个点另一条可能有一千个点。这时候无论用什么距离度量短轨迹都容易成为“孤岛”被单独分到一个簇里。我的处理思路是先抽稀再编码。Matlab自带reducepoly函数可以基于Douglas-Peucker算法做轨迹简化。抽稀到统一的关键点数范围再进网格编码流程。另外在预处理阶段可以加一个过滤条件点数少于10的轨迹直接剔除。这种轨迹很可能只是GPS冷启动的噪声不具备聚类价值。5.2 问题聚类结果每次跑都不一样高频出现的问题。原因大概率是初始中心随机导致局部最优。看到一个奇怪的现象别慌先检查代码里是否设了Replicates。设了之后还是不稳定检查数据本身——如果某些簇之间距离太近Kmeans会在它们之间反复横跳本质是你的K值选大了。另一个容易忽略的点是MapReduce或并行环境下的随机数变量问题。如果有并行计算记得用parfor时在循环内部单独设置随机种子否则不同worker的随机序列可能相同。5.3 问题离群点拖垮聚类质量轨迹数据里经常有“幽灵轨迹”——GPS信号在高楼林立的区域来回反射生成一串毫无规律的折线。这些轨迹跟谁都不像却被Kmeans强行分配给某个簇还会把簇中心拉偏。经验做法是聚类前用轮廓系数做一轮离群点筛查先跑一次Kmeans计算每个样本的轮廓值把轮廓值低于0的样本剔除掉再重新聚类。这是最省事的离群点处理策略效果立竿见影代价只是多跑一轮聚类。5.4 问题数据量大导致内存爆炸距离矩阵是O(n²)的。1万条轨迹就要1亿个浮点数占用约800MB内存还能勉强扛住5万条轨迹就是25亿个浮点数20GB内存直接报警。我的解决方案是分块计算距离并只保存距离矩阵的上三角或者干脆用single精度D zeros(numTraj, numTraj, single); % 单精度省一半内存更进一步的做法是不存完整距离矩阵改用特征向量直接Kmeans。网格编码后的每条轨迹可以展开成一个gridSize*gridSize维的0/1向量用稀疏矩阵存储内存占用会下降一两个数量级。这是从“轨迹距离矩阵法”切换到“轨迹特征向量法”的重大收益。5.5 问题Matlab中文字符乱码不少同事在跑代码时遇到中文注释乱码原因通常是文件编码不是UTF-8。在Matlab R2021a之后推荐在“预设项-编辑器-语言”里把文件编码统一设成UTF-8。如果代码里有中文字符串常量比如title(平均轮廓系数)设置后需要重新打开脚本生效。这个坑不影响聚类结果但影响团队协作的体验提一笔。6. 扩展方向与个人实践体会6.1 从Kmeans升级到两步聚类法基于距离矩阵做Kmeans只是轨迹聚类的基操。我在真实项目里更爱用两步法第一步用DBSCAN基于网格编码密度剔除离群轨迹、找出大致的稠密簇第二步在簇内用Kmeans做精细切分。两步法的好处是既解决Kmeans对噪声敏感、又不至于像纯DBSCAN那样在高维空间里效率低下。6.2 把Kmeans换成K-Shape考虑时序形态如果轨迹数据本身带有时间戳而你希望聚类时把“时间节奏”也纳入考量可以了解下K-Shape算法。它专门针对时间序列的形态相似性设计用交叉相关作为相似度度量配合尺度不变的处理比KmeansDTW的思路快很多。Matlab社区里已经有人封装了K-Shape的实现搜一下就能找到。6.3 谈谈轨迹可视化的小技巧聚类出来之后单纯把所有轨迹画在一张图上会非常拥挤。我的习惯是每个簇单独画一张子图并在图例里标注簇内轨迹数量和平均轮廓系数。还有一个实用技巧取簇中心轨迹离簇中心最近的轨迹单独高亮它能最直观地代表这一个簇的典型路径形态。6.4 最后说一点个人感悟做轨迹聚类这两年我最深的一个体会是这类项目八成时间花在“定义什么是相似”上而不是花在跑算法上。距离度量一旦确定Kmeans本身只是最后一公里的工具。所以如果你正在做一个轨迹聚类的项目先别急着写kmeans函数多跟业务方确认清楚——“你要的相似是路径形状相似还是起终点相似还是时间节奏也一致”这个问题的答案直接决定了你的预处理方案选哪条路。我在初期做分享的时候也栽过跟头拿着一套DTW方案硬套一个需要实时处理几万条轨迹的系统算力完全扛不住。后来改成网格编码加Jaccard距离效果虽然损失了一点精度但性能和稳定性都上来了。算法选型永远是一个权衡没有绝对的好坏。这篇代码和思路算是一个比较平衡的起点希望你在真实数据上跑一跑、调一调形成自己的套路。