1. 为什么Kmeans对原始轨迹坐标“失效”了三条最容易犯的错大概两年前我接过一个共享电单车的骑行轨迹分析需求客户只说了一句“我们用Kmeans聚一下就行”。我当时心里就咯噔一下因为轨迹聚类这个事最容易翻车的地方恰恰就在“直接聚”这三个字上。后来果不其然他们对原始GPS坐标点跑了Kmeans出来的图完全看不出任何骑行模式每条轨迹都像被剪刀剪碎了一样东一块西一块。这件事之后我养成了一个习惯拿到轨迹数据的第一反应不是写Kmeans而是先想清楚“我要聚的对象到底是什么”。Kmeans本身只是一个在欧氏空间里找簇心的算法它天生处理的是“点”和“点”之间的距离。但轨迹是一条带顺序和时间的序列同样是从A点到B点有人走直线有人绕路有人在某个路口等了两分钟红灯有人中途停下去买东西。这些行为差异如果直接反映在坐标点上Kmeans根本区分不出“这是一条绕路的轨迹”和“这是两个不同地点的散点”。所以这篇文章要解决的问题就是怎样把Kmeans正确用在轨迹数据上并且给出可以直接跑的Matlab代码。我把新手最容易犯的错总结成三条先看明白这三条后面再谈算法路线和代码思路会顺很多。1.1 轨迹长度不一致特征矩阵根本构造不起来这是最朴素也最致命的问题。Kmeans的输入是一个二维矩阵每一行是一个样本每一列是一个特征。但轨迹数据天生长短不一A用户骑了5分钟GPS采了30个点B用户骑了40分钟采了300个点。你没法直接把这两条轨迹放进行矩阵因为维度都对不上。大多数人的第一反应是“补零”或“截断到最短”。补零的问题在于零本身会成为一个有意义的特征Kmeans会莫名其妙多分出一类专门吸收那些被补了超多零的短轨迹截断更粗暴直接把轨迹后半段的形态扔掉了比如那条绕路去充电桩换电的轨迹后半段恰恰是它和其他轨迹最大的区别。正确做法是先对轨迹做弧长重采样把所有轨迹统一成固定数量的坐标点。这一步做好了后续的Kmeans才有“矩阵”可吃。1.2 把所有坐标点堆在一起聚类出来的不是轨迹簇而是热力点这是“假聚类”里面最多人踩的坑。有些人知道轨迹长度不一但嫌重采样麻烦索性把100条轨迹的上万个点全部倒进同一个矩阵然后跑Kmeans。跑完之后聚类结果确实很干净簇与簇之间边界明显但每个簇的含义是“这些坐标点在空间上挨得比较近”不是“这些轨迹在行为模式上相似”。举个很简单的例子三条轨迹都从小区门口出发前500米完全重叠然后分道扬镳一条去地铁站一条去菜市场一条去公园。把所有点混在一起聚类前500米的点一定被划成同一簇因为它们在空间上紧密抱团而后面的点按照目的地被切成两三簇。最后你得到的结果名为“轨迹簇”实为“地理热点”。如果你想分析的是“早上从小区出发的人分别去哪”这个结果当然可用但如果你想分析的是“三条完整轨迹的走法结构”这种聚类方式从根上就错了。1.3 采样率不一致同名点对之间的“假距离”让Kmeans无从下手就算你统一了轨迹长度还有一个隐蔽问题采样率不一致。假设两条轨迹都从A地到B地A用户的车是10Hz高频率上报B用户的车是1Hz低频上报。两者重采样到同样50个点之后理论上看起来可比了但如果中间有转弯高频轨迹记录下了细腻的弧线低频轨迹则是一条折线对应点之间的空间距离会非常大。Kmeans算欧氏距离时这种误差会被放大结果就是同一条路线被分成两簇仅仅因为采样率不同。我的判断标准很简单如果你的聚类目标是“整条轨迹的形态或行为”Kmeans就不能直接吃原始坐标点如果目标是“找出经常被经过的区域”那直接用坐标点聚类没问题。想清楚这一点你才敢往下走。2. 轨迹聚类的三条技术路线选路线比写代码更重要前面说的都是“不能怎么做”接下来聊“该怎么做”。轨迹聚类没有唯一标准答案我在实际项目里一般会根据数据条件和业务诉求在三条技术路线里选一条。这三条路线在工程里都非常常见选对了能省掉大量返工时间。2.1 路线A弧长重采样加展平特征向量这是本文Matlab代码采用的主路线也是我给大多数项目做“第一版结果”时的首选。核心思路分成三步对每条轨迹做弧长重采样统一成N个坐标点把N个点的横纵坐标拼接成一个1行2N列的向量作为这条轨迹的特征对这个特征矩阵跑标准Kmeans。这个方案最大的好处是简单粗暴且可解释每条轨迹的每个重采样点都有固定含义特征向量里的第i个位置就是“整条轨迹在里程第i段处的坐标”。Kmeans聚类后簇心向量可以重新折叠成一条轨迹这就是这个簇的“平均轨迹”画出来非常直观。缺点是对轨迹的形变和错位比较敏感两条形态相似但起点终点偏移较大的轨迹可能因为末端坐标差过大被拆开。2.2 路线B距离矩阵加MDS或谱聚类如果你觉得轨迹形态差异很明显起点终点都不固定用路线A效果不好那我建议用路线B。先算任意两条轨迹之间的相似度距离得到一个距离矩阵D再通过多维缩放或者谱聚类把轨迹分组。这个路线不要求轨迹等长你可以用DTW距离、Hausdorff距离、最长公共子序列距离甚至自定义业务距离非常灵活。缺点是计算复杂度和存储开销都不小。60条轨迹的距离矩阵是60×60600条轨迹就是600×600如果每条轨迹的DTW计算都要几十毫秒整体就是几分钟起步。所以这个路线更适合离线分析不适合在线实时聚类。2.3 路线C手工特征压缩轨迹再跑Kmeans如果业务上面要求“每个簇必须讲得出商业含义”那路线C往往是最明智的选择。不要拿整条轨迹去聚类而是先从轨迹里提取出有业务含义的标量特征全程平均速度、最大瞬时速度、总里程、转弯次数、平均航向角变化率、起终点直线距离与里程的比例等。把这些特征拼成一个向量再做Kmeans。这种做法的好处是聚类结果可以直接落到业务动作上。“簇1是速度高、转弯少的长途直行轨迹簇2是速度低、转弯多的短途巡游轨迹”这句话讲给运营听对方一下就懂了。缺点是特征怎么设计很考验经验特征选得不好聚类结果等于随机分组。2.4 三条路线怎么选一张表讲清楚场景特征推荐路线理由轨迹长度差不远采样率经过预处理路线A实现成本最低结果最稳定轨迹形态差异大起点终点各不相同路线B距离度量更灵活能捕捉形变业务要求解释性要直接给运营用路线C特征向量语义清晰可落行动数据量在几百条以下离线分析路线A或B计算量可控Debug方便数据量上万条在线聚类路线A特征固定后可直接用增量式Kmeans3. 预处理是命门弧长重采样、数据清洗、坐标归一化的细节很多人在Kmeans上报错“矩阵维度不一致”就开始换框架、换语言其实问题出在前面处理原始轨迹时没做好标准化的动作。预处理不是可有可无的步骤它对聚类结果的影响比Kmeans参数大一个数量级。3.1 弧长重采样为什么必须按里程等间隔取点弧长重采样的原理说白了就是按照轨迹的总长度把整条轨迹等分成若干段在等分点处取坐标。这样做的好处是彻底消除采样率和速度差异的影响。拿两条从公司到家的轨迹举例一条是早高峰堵车车速慢GPS点位密集另一条是半夜畅通车速快GPS点位稀疏。如果不重采样点位密度直接反映车速Kmeans会误把“堵车时段”当成一个聚类特征重采样之后每条轨迹都用等量的点描述相同的地理距离比较的就只剩下“路线形态”了。Matlab里实现弧长重采样很直接我自己常用的函数长这样function trajN resampleTrajectory(traj, N) % traj: n x 2 的轨迹点列例如经投影后的平面坐标 % N: 重采样点数也就是统一后每条轨迹的点数 % trajN: N x 2 的重采样结果按全程弧长等间隔取点 if size(traj, 1) 2 error(轨迹点数太少无法重采样); end segLen sqrt(sum(diff(traj, 1, 1).^2, 2)); cumLen [0; cumsum(segLen)]; totalLen cumLen(end); if totalLen 0 trajN repmat(traj(1,:), N, 1); return; end tq linspace(0, totalLen, N); % interp1 会对矩阵的每一列分别插值 trajN interp1(cumLen, traj, tq, linear); end关于重采样点数N怎么选我的经验是看轨迹里最短的那条。最短轨迹如果有20个点那你重采样到30~50个点是安全的信息不会丢太多如果最短轨迹只有5个点那这数据本身质量就存疑硬重采样到50只会把噪声插值得漂漂亮亮没意义。3.2 清洗跳变点和NaN先删脏数据再谈聚类GPS轨迹常见的脏数据有两大类一类是NaN丢星另一类是跳变点。跳变点的典型表现是相邻两个采样点之间距离超大比如1秒内“瞬移”了500米这通常是GPS漂移造成的。清洗跳变点的方法也比较朴素先计算每相邻两点的速度设定一个阈值比如超过整条轨迹速度中位数的5倍就判定为跳变把跳变点删除再用前后有效点线性插值补上。Matlab里判断跳变的核心代码段我贴出来dt diff(t); % t 是时间序列单位秒 segSpeed sqrt(sum(diff(xy, 1, 1).^2, 2)) ./ dt; medSpeed median(segSpeed); outlierIdx find(segSpeed 5 * medSpeed); % 设定5倍中位数为跳变阈值注意不要直接用绝对速度阈值因为不同场景的速度分布差太远了。步行轨迹的中位数速度只有1.5m/s左右机动车轨迹的中位数速度轻松上10m/s用绝对阈值会误杀。3.3 坐标归一化比你想的更微妙如果你已经在轨迹数据上做了重采样那么最后的归一化步骤直接影响Kmeans的聚类形状。这里有个细节经常被忽略到底是对每一列单独zscore还是对全部坐标做统一缩放我的建议是如果你要保留轨迹的形状比例就做统一缩放不要对x坐标和y坐标分别标准化。分别标准化等于把横轴和纵轴的尺度强行拉成一样圆形轨迹会被压成椭圆直角三角形会被压成斜边变短的畸变三角形。正确的做法是先求出所有轨迹在x和y方向上的整体标准差然后用同一个尺度因子缩放。实际操作里我更多是直接对展平后的特征矩阵做一次统一标准化X reshape(trajRes, numTraj, N * 2); scaleFactor std(X(:)); X X ./ scaleFactor;这样每个特征维度的相对比值还保留着原始的物理意义同时又让量级进入了Kmeans友好的范围。4. 相似度度量RMSD、DTW与Hausdorff的取舍逻辑Kmeans的实质是基于距离的迭代优化所以距离度量才是那个真正决定聚类结果的东西。很多人的代码跑得很熟练但从不回头看一眼自己用的是哪种距离出了奇怪结果也不知道从哪查起。4.1 RMSD路线A背后的默认度量在重采样之后的轨迹空间里最自然的距离是逐点距离的平均值也就是RMSDRoot Mean Square Deviation。两条重采样轨迹都是N×2的矩阵相减之后求每个点的欧氏距离然后取平均。d sqrt(sum((trajA - trajB).^2, 2)); dRMS mean(d);RMSD的好处是计算快而且对小幅噪声有天然平滑效果。但它对时间轴错位非常敏感如果一条轨迹在某个弯道比另一条早了10米开始转弯逐点对应时转弯前的直线段会被强行错位比较距离值虚高。这也是为什么路线A对“形变”不友好。4.2 DTW路线B最常用的距离代码不长但逻辑要懂DTW动态时间规整解决的就是RMSD对时间轴错位的毛病。它的思想是允许轨迹点按顺序“错位对齐”代价最小的对齐方式就是最优匹配。比如两条轨迹同样是绕个S弯一条弯得早、一条弯得晚DTW能够正确地把两边的弯心匹配上距离自然就小了。Matlab实现版本很多我自己用的简洁版本如下function d dtwDist(tr1, tr2) % tr1, tr2: n1 x 2, n2 x 2 的两条轨迹点列 n1 size(tr1, 1); n2 size(tr2, 1); D pdist2(tr1, tr2); % 点对点距离矩阵 C inf(n1 1, n2 1); C(1,1) 0; for i 1:n1 for j 1:n2 % 只能向前对齐不能回头 C(i1,j1) D(i,j) min([C(i,j1), C(i1,j), C(i,j)]); end end d C(n11, n21); end这个双重循环实现直观但坏消息是复杂度是O(n1×n2)轨迹一长、样本一多就跑不动。实际项目里我一般只在轨迹数量不超过200条时用DTW全量距离矩阵。如果想提速可以限制扭曲窗口宽度比如只允许点对在前后20%长度范围内对齐效果几乎不变速度能快好几倍。4.3 Hausdorff距离和Frechet距离一句话讲清楚差异在路线B里还可以用Hausdorff距离或Frechet距离。Hausdorff距离的定义是“轨迹A上每个点到轨迹B的最近距离的最大值”它的特点是极其关注两条轨迹的最大偏离程度。这个性质很双刃它擅长找出“某一段明显背离”的轨迹但也很容易被单个漂移点带偏。Frechet距离则更接近“人在遛狗狗绳拉直时最短能多短”的模型它同时约束了点位的顺序性和连续性比Hausdorff更符合轨迹形态比较的直觉。但Frechet的计算比DTW更复杂Matlab没有内置函数需要自己实现离散版本工程成本高一般我不推荐第一个版本就上它。4.4 实际该用哪个结合你的数据特性拍板如果你的轨迹已经做了弧长重采样并且采样率差异不大直接用RMSD就行Kmeans的收敛速度快结果也好解释。如果轨迹在空间上有明显的“先经过A再经过B”的顺序结构但时间节奏不同用DTW。如果只想粗筛“有没有某一段轨迹严重偏离”用Hausdorff。如果业务上把轨迹当成“一条绳”要求整条路径的形状逼真度再考虑Frechet。5. 完整Matlab实现从模拟轨迹到Kmeans聚类的一站式脚本接下来是这篇文章的重头戏一份可以直接运行的Matlab脚本。为了让你不依赖外部数据就能看到完整效果我先生成三组形态不同的模拟轨迹然后做弧长重采样、特征展平、Kmeans聚类和可视化。你把自己真实的轨迹数据套进对应的步骤就行。5.1 主脚本全流程一口气跑通%% 基于Kmeans的轨迹聚类完整示例 clear; clc; close all; rng(42); %% Step 1: 模拟生成三种形态的轨迹 numTraj 90; % 总轨迹条数 trajCell cell(numTraj, 1); for i 1:numTraj if i 30 % 第一类平缓直线型 nPts 18 randi(10); t linspace(0, 1, nPts); x 0.8*t 0.04*randn(nPts,1) 0.1*randn; y -0.3*t 0.04*randn(nPts,1) 0.1*randn; elseif i 60 % 第二类上凸曲线型 nPts 16 randi(12); t linspace(0, 1, nPts); x t.^2 0.05*randn(nPts,1); y sin(1.5*pi*t) 0.06*randn(nPts,1) 0.1*randn; else % 第三类先直后折的折线型 nPts 24; t linspace(0, 1, nPts); x min(0.6*t, 0.6) 0.03*randn(nPts,1); y 0.8*t 0.03*randn(nPts,1); y(t 0.6) 0.48 1.2*(t(t0.6)-0.6) 0.03*randn(sum(t0.6),1); end trajCell{i} [x, y]; end %% Step 2: 弧长重采样到统一点数 N 50; % 统一后的轨迹点数 numTraj length(trajCell); trajRes zeros(numTraj, N, 2); for i 1:numTraj trajRes(i,:,:) resampleTrajectory(trajCell{i}, N); end %% Step 3: 展平成特征向量并做统一缩放 X reshape(trajRes, numTraj, N*2); scaleFactor std(X(:)); X X ./ scaleFactor; %% Step 4: 用轮廓系数扫描K值K从2到6 Krange 2:6; silScores zeros(1, length(Krange)); for k 1:length(Krange) idxTmp kmeans(X, Krange(k), Replicates, 10); silTmp silhouette(X, idxTmp); silScores(k) mean(silTmp); end [~, bestKPos] max(silScores); K Krange(bestKPos); %% Step 5: 用最优K值跑最终Kmeans [idx, C] kmeans(X, K, Replicates, 20); %% Step 6: 分别计算每个簇的平均轨迹用于可视化簇中心 avgTrajCell cell(K, 1); for k 1:K members find(idx k); avgTrajCell{k} squeeze(mean(trajRes(members, :, :), 1)); end %% Step 7: 可视化 figure(Position, [100 100 1200 420]); colors lines(K); subplot(1,3,1); hold on; grid on; for i 1:numTraj traj_i squeeze(trajRes(i,:,:)); plot(traj_i(:,1), traj_i(:,2), Color, [0.65 0.65 0.65]); end title(重采样后的所有轨迹); axis equal; subplot(1,3,2); hold on; grid on; for i 1:numTraj traj_i squeeze(trajRes(i,:,:)); plot(traj_i(:,1), traj_i(:,2), Color, colors(idx(i), :), LineWidth, 1.0); end title([Kmeans聚类结果K, num2str(K)]); axis equal; subplot(1,3,3); hold on; grid on; for k 1:K plot(avgTrajCell{k}(:,1), avgTrajCell{k}(:,2), ... Color, colors(k,:), LineWidth, 2.5); end title(每个簇的平均轨迹); axis equal;5.2 一个容易被可视化带歪的细节簇中心来自标准化空间第一次用Kmeans做轨迹聚类的人很容易直接把kmeans返回的簇心C拿来画“中心轨迹”画出来的图往往乱七八糟。原因在于我们喂给kmeans的特征矩阵X是经过统一缩放的簇心C也是缩放空间里的坐标直接reshape回轨迹形状时横纵坐标已经脱离了原始数据的物理尺度。所以在Step 6里我没有用C去画图而是重新取出每个簇内的原始重采样轨迹按簇内成员做平均。这样得到的平均轨迹才是“这个簇的代表性路线”画出来和原始坐标对得上业务人员看了也不会懵。这个细节我强烈建议你保留因为它直接影响你给同事或客户汇报时的可信度。5.3 拿到聚类结果之后至少看一眼簇内的平均轨迹代码跑通只是开始。我的习惯是每次聚类结束第一件事不是看轮廓系数而是把每个簇的平均轨迹和两三条典型单条轨迹叠在一张图上快速确认“这个簇的代表性轨迹是否符合直觉”。如果平均轨迹杂乱无章说明K值选大了或者预处理出了问题如果平均轨迹平滑清晰才值得继续往下做业务分析。6. K值选择与聚类质量评估轮廓系数、稳定性与业务校验Kmeans有个绕不开的宿命K要你自己给。很多教程讲到这里就扔给你一个轮廓系数让人背公式但我更愿意把K值选择理解成“聚类结果的复现性校验”。因为轨迹数据往往没有标准答案同一个K在不同初始点下跑出来的稳定程度才是更靠谱的评估指标。6.1 轮廓系数Matlab一行搞定但别只取一个值轮廓系数的计算方法你自己写也不难对每个样本算它到同簇其他样本的平均距离a再算它到最近其他簇所有样本的平均距离b轮廓系数就是(b-a)/max(a,b)。取值范围从-1到1越接近1说明样本离自己簇越近、离其他簇越远。Matlab里直接调用内置函数sil silhouette(X, idx); meanSil mean(sil);但我不建议只扫描一次就拍板K。脚本里我循环了K2到6然后取了平均轮廓系数最大的K。这里有个隐形问题轮廓系数容易偏袒K小的情况因为簇少的时候簇间边界天然更清晰。所以我通常会把轮廓系数排名前两三名都跑一遍再对比可视化结果选那个业务上最有意义同时轮廓系数也不差的K。6.2 肘部法则看簇内距离平方和随K的变化肘部法则是另一种常见的K值选择方法看的是簇内距离平方和WSS随K增大的下降趋势。下降速度骤减的那个拐点就是“肘部”。Matlab计算WSS可以手动实现for k Krange idxTmp kmeans(X, k, Replicates, 10); wss(k-1) 0; for kk 1:k members X(idxTmp kk, :); center mean(members, 1); wss(k-1) wss(k-1) sum(sum((members - center).^2)); end end plot(Krange, wss, o-);不过说实话真实轨迹数据里WSS肘部经常不是那么明显曲线平滑下降很难说哪个点才是肘。所以我的策略是肘部法用来圈定候选范围轮廓系数用来进一步筛选最后用稳定性拍板。6.3 稳定性校验同一K在不同初始化下是否总是给出相似分组这个指标很多教程不提但在业务项目里非常实用对同一个K用不同的随机种子跑10次Kmeans然后比较10次结果的相似度。相似度高的说明这个K和这组数据是“匹配”的相似度低说明数据在这个K下本来就分不清任何结论都不可靠。Matlab里实现也不复杂用RandStream控制随机种子跑10次聚类之后逐样本比较标签的一致性。如果标签之间的平均互信息AMI低于某个阈值我就直接建议用户换距离度量或者换路线而不是继续调K。6.4 别让数学指标凌驾于业务之上做轨迹聚类到最后我极少遇到“数学上的最优K”和“业务上的最优K”完全一致的情况。有一次共享单车项目的数据轮廓系数最高指向K6但运营那边实际只有三种可执行的调度策略K6意味着每种策略要被拆成两簇操作上根本无法落地。后来我折中选了K4在多簇之间加了细分标签轮廓系数从0.55掉到0.47但每个簇都能对应到明确的运营动作。这件事给我的教训是K值选择归根结底是服务业务目标的决策数学指标是工具不是判决书。7. 真实轨迹数据上的踩坑记录采样率、停留点与离群点的处理最后这部分是这些年在真实项目里积累下来的一些经验。模拟数据再怎么完美到了真实轨迹场景该遇到的坑一个都少不了。我把最典型的几个记录下来希望你不用重新踩一遍。7.1 采样率差异的威力以出租车GPS数据为例出租车GPS数据的采样间隔并不是恒定的空载省电模式下可能是30秒一条载客接单后变成5秒一条。两条从机场到市中心的轨迹如果一条是空载状态录的一条是载客状态录的重采样之前算距离数值可以差到几公里。弧长重采样能解决点位数量不均的问题但没法彻底解决“重要拐弯处点数少”的问题。我通常会额外做一步先按时间插值到1秒间隔再做弧长重采样。这样等于先补密再统一拐弯处的形态特征保留得更好。7.2 停留点会扭曲轨迹形状必须识别并剁掉轨迹里最坑人的是停留行为。用户在某地停了20分钟GPS以1Hz频率一直在上报这段轨迹在空间上表现为一个点团。如果你把整条轨迹拿去重采样这段点团会占掉重采样点里很大比例于是这条轨迹的“形状”被这个停留点完全主导。聚类时不管它真正的行驶路线是什么样子都会被分到“停留时间长的轨迹”那一簇去。我的处理方式是先算每个相邻点对的速度把速度低于0.5m/s的连续片段标记为停留段然后将停留段内的点做降采样比如每20个点保留1个或者干脆把停留段的中心点作为单点保留。具体保留策略要看业务如果分析的是路径规划停留点直接删掉如果分析的是出行行为停留时长本身反而是一个特征那就把它提取出来不要混在坐标里。7.3 离群点Kmeans最怕那种“又长又怪”的轨迹离群轨迹对Kmeans的影响比离群点对大得多。比如100条正常轨迹里混进来一条把整个城市的对角线都跑了一遍的轨迹它的长度和形状都极端重采样之后这个样本在特征空间里距离其他样本极远。这会导致两个后果要么它单独成一簇K被它浪费掉一簇要么它强行拽动某个簇心把正常轨迹也带偏。我处理离群轨迹的顺序是先用简单的长度和里程阈值过滤掉明显异常的轨迹再用DBSCAN在展平特征空间上做一次粗聚类把落在任何簇外或者簇很稀疏的样本标记为离群。最后才把这些干净样本送入Kmeans。这样不仅结果稳还能在报告里多出一个“离群轨迹”类别很多情况下这个类别反而能发现异常驾驶或设备故障比正常聚类结果更有业务价值。7.4 在线跑Kmeans时的坑中心漂移与增量更新最后提醒一下想做在线轨迹聚类的朋友。Kmeans的原始版本是批处理的每来一批新轨迹就要重新跑一遍全量数据随着数据量增长越来越慢。实际工程中我多数时候改用MiniBatchKMeans或者在线式增量更新新样本进来后先算它到现有簇心的距离归入最近簇同时按一定学习率更新簇心。Matlab里没有内置MiniBatchKMeans但你可以自己写一个几十行的更新循环核心就是那个簇心更新公式。如果你的轨迹流是实时上报的这个方向会比反复调用全量kmeans靠谱得多。……关于代码最后再补充一点如果你自己只有经纬度坐标建议先把经纬度转换成平面坐标再跑聚类。Matlab的Mapping Toolbox有deg2km或geodetic2enu这类函数没有工具箱的话用等距投影的近似公式也够用。这个步骤不做你在真实地理数据上计算欧氏距离误差在低纬度地区还勉强能接受到了高纬度地区单位经度对应的地面距离变化非常大聚类结果会失真到没法看。我个人的习惯是拿到轨迹后的第一件事永远是画图先把所有轨迹叠一张图看一遍再决定用什么距离、什么路线、什么K。这个习惯帮我避开了至少一半的无效调试希望你也能用上。