前阵子接了个需求要把一批共享电单车的GPS轨迹按空间相似度归类。拿到数据的第一时间我脑子里冒出来的方案就是sklearn.cluster.KMeans几行代码就能出结果看起来轻轻松松。但真正动手之后才发现GPS轨迹这种数据根本没法和普通表格数据一样直接丢进KMeans。轨迹是长度不等的点序列每一条都有自己独特的顺序、采样间隔和噪声形态。这篇内容我会从“为什么不能直接KMeans”讲起把GPS轨迹预处理、空间相似度度量、K-medoids替代方案和完整Python实现一次讲透适合正在做轨迹挖掘、异常检测、路径归类的同学参考。1. 为什么轨迹聚类不能直接把坐标丢进KMeans1.1 轨迹数据比普通特征多出来的“序列属性”普通聚类任务里一行数据就是一个样本一列是一个特征样本之间可以一一对应。GPS轨迹不是这样一条轨迹是连续采集到的经纬度点序列每秒钟一个点、五分钟一个点、甚至路段中间掉线一段得到的序列长度完全不同。如果我天真地把一条轨迹的所有点拍平成一行第一条轨有200个点第二条只有50个点连特征维度都对不上模型根本无从下手。更麻烦的是GPS轨迹不仅有点坐标还有一个隐含的“顺序”。同一条路走一遍和走两遍再回来点集合可能类似但空间语义完全不同直行和绕路折线的形状不同相似度计算时也必须考虑顺序带来的差异。1.2 空间相似度到底在比什么很多人提到GPS空间相似度第一反应就是“计算两组经纬度点之间的距离”。这个想法对但不完整。两组轨迹要比较的不只是每个点的位置还包括整体形态、运动方向、以及路径在空间上的覆盖范围。举个例子外卖员从A点到B点一条轨迹是直行另一条先往东绕一圈再回到主干道。两条轨迹的起点、终点完全一样但空间形态差异很大。到底算相似还是不相似取决于你的业务目标。如果关心的是“货送到没有”起点终点足够如果关心的是“路线是否规范”就必须把中间形态也纳进来。这个问题决定了后面选哪种距离函数也会直接影响到聚类结果。我不建议一开始就追求复杂算法而是先明确你要聚类的是“位置范围”还是“走的路径形状”。2. 计算距离之前先把GPS轨迹整理成同类项2.1 坐标格式与坐标系WGS84和投影坐标的选择原始数据里的经纬度看起来是数字但藏了很多坑。有些设备导出的是度分秒格式比如116°2345.6有些是十进制小数如果不同来源混在一起必须统一成十进制度再计算距离。我自己的处理逻辑是度分秒转十进制度十进制度 度 分/60 秒/3600。统一用WGS84坐标系这是绝大多数GPS设备的默认坐标系。如果计算距离的范围在同一个城市距离小时直接用经纬度算Haversine距离即可。如果要把轨迹做成标准KMeans的向量就必须投影到米制平面坐标比如UTM否则直接用经纬度坐标算欧氏距离纬度方向和经度方向尺度不一致聚类结果会变形。2.2 噪声点与漂移点过滤按速度、跳变和时间戳判断GPS设备在城市里经常出现漂移点会突然跳到几百米外然后又跳回来。这种点如果不处理会在计算两个轨迹的相似度时重新引入节点不匹配项导致本应相似的轨迹被判成不相似。我的过滤规则比较简单对每条轨迹按时间排序计算相邻点之间的距离用距离(m) / 时间差(s)估算瞬时速度如果速度超过合理阈值比如骑行超过120km/h、机动车超过200km/h基本可以断定是漂移点或设备故障直接切除对静止不动的点做“停留点压缩”如果一个半径20米内的点持续超过5分钟只保留首尾两个点因为DTW计算时大量重复点会把相似度拉向“停留多不多”而不是“路径像不像”。2.3 轨迹抽稀和重采样给后面的计算减负实时采集的轨迹点非常密一条10公里骑行轨迹可能有几千个点。如果直接拿这几千个点去做两两DTW距离数据量一大计算会慢到怀疑人生。所以我会做两步预处理抽稀用Douglas-Peucker算法以20米为距离阈值把落在“折线包围盒”内部的冗余点去掉保留关键转弯点。重采样根据轨迹总长度等比例采样到固定点数比如每条轨迹统一取64个点。这样有两个好处一是减少DTW里的序列长度二是如果之后想做向量化KMeans也具备条件。需要说明的是重采样不是简单地从原序列里每隔N个取一个点而是按照累计距离在折线上插值。比如把整条轨迹按总长度分成63段每段末端取一个点。这样快慢、停不停对采样密度的影响被抹平了留下的只是空间形状本身。3. 空间相似度度量Haversine、Hausdorff、DTW怎么选3.1 点与点之间的距离Haversine公式不能省所有轨迹相似度的底层都是经纬度点之间的距离。在城市尺度下之前有同事用球面余弦公式算出来相差不大但为了稳妥我在所有核心代码里都用Haversine公式。它的逻辑是把地球看成半径6371公里的球把两点转换到球坐标后求弧长。公式写成代码很直接import numpy as np EARTH_RADIUS_KM 6371.0 def haversine(lat1, lon1, lat2, lon2): lat1, lon1, lat2, lon2 map(np.radians, [lat1, lon1, lat2, lon2]) dlat lat2 - lat1 dlon lon2 - lon1 a np.sin(dlat / 2) ** 2 np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2) ** 2 c 2 * np.arcsin(np.sqrt(a)) return EARTH_RADIUS_KM * c返回单位是公里。如果后续要算速度我会把公里转成米再除以时间差。3.2 Hausdorff距离衡量两条轨迹的“最大偏离”Hausdorff距离是集合之间的经典距离对轨迹A上的每一个点找它在轨迹B上最近的点的距离然后取最大值。反过来也要算一次最后取对称值中的最大值。这个指标很适合回答“A跑到了B完全没有覆盖过的区域吗”这类问题。但它对离群点非常敏感。GPS只要出现一个漂移点Hausdorff距离就可能瞬间变成几百米把两条明明同路的轨迹拉开。我在实跑中很少单独用它只有在做“是否经过同一片范围”这类粗粒度判断时才考虑。3.3 DTW距离解决长度不同和节奏不同的关键DTWDynamic Time Warping是我在项目里最常用的算法。它的优势在于它允许两条轨迹在时间轴上“错位对齐”。你走得快、我走得慢或者你中间停了一分钟这些都不会破坏整体形态的相似度。核心递推式是D[i][j] d(p_i, q_j) min(D[i-1][j], D[i][j-1], D[i-1][j-1])其中d(p_i, q_j)是Haversine距离。从左上角一直算到右下角最后的值就是两条轨迹的DTW距离。代码实现其实不复杂def dtw_distance(traj1, traj2, dist_fnhaversine): n, m len(traj1), len(traj2) dp np.full((n 1, m 1), np.inf) dp[0, 0] 0.0 for i in range(1, n 1): for j in range(1, m 1): p traj1[i - 1] q traj2[j - 1] cost dist_fn(p[0], p[1], q[0], q[1]) dp[i, j] cost min(dp[i - 1, j], dp[i, j - 1], dp[i - 1, j - 1]) return dp[n, m]复杂度是O(n*m)。如果两条轨迹重采样到64个点一次计算只有4096次Haversine调用完全能接受。但如果两条轨迹一个3000点、一个2500点那就是750万次基础距离计算跑起来非常慢。所以前面说抽稀和重采样不是为了好看而是让DTW能在合理时间内跑完。3.4 我的选型逻辑不同距离函数背后的语义差别很大。我整理了一张表方便各位按业务场景对照参考距离方式解决的核心问题主要缺点适合场景Haversine点对点计算两个点在地球表面的距离无法比较变长序列底层基础计算Hausdorff衡量两条轨迹覆盖空间的最大偏差对GPS漂移点极敏感粗粒度区域覆盖比较DTW解决长度不同、速度不同带来的错位复杂度较高需重采样同路不同速度的轨迹归类LCSS容忍点和点之间缺失、噪声参数多阈值难定设备丢点严重的数据在这套轨迹聚类项目里我的主距离函数选择是DTW。因为共享电单车轨迹最大的特点就是“同一条路骑行速度不同、停留时间不同”DTW能把这些干扰因素压到最低。Hausdorff只作为过滤异常轨迹的辅助手段如果两条轨迹的Hausdorff距离已经小到阈值以下那么DTW大概率也很近可以跳过复杂计算。4. K-means遇到轨迹质心失效与K-medoids替代4.1 标准K-means为什么要算“质心”KMeans的精髓是反复迭代先把样本分给最近的质心再重新计算每个簇的质心。在普通特征空间里质心就是该簇所有样本的坐标平均值计算简单、含义明确。但轨迹数据的质心是什么如果每条轨迹都被重采样成64个固定点理论上我可以把64个点取平均得到另一条64个点的“平均轨迹”。实际跑一次你会发现很多点在平均后落到了两条路之间变成一段没有任何车辆走过的“虚拟道路”。尤其是十字路口、双向车道这类地方平均轨迹会在两条道路中间画出一条离谁都近、又谁都不像的线。这个问题在路网匹配场景里尤其严重。项目里几条轨迹分别走一条双向道路的南北两侧平均后形成的质心轨迹落在了隔离带中间后续再用这条质心轨迹做路线推荐业务上完全不可解释。4.2 K-medoids用真实轨迹当质心既然平均轨迹不可靠一个直接的想法就是不用平均直接从每个簇里选一条最有代表性的真实轨迹当“质心”。这就是K-medoids的基本思想。它和KMeans的唯一区别在于“质心的更新方式”KMeans质心是簇内所有样本特征的均值K-medoids质心是簇内与其他样本距离之和最小的那一个真实样本。用DTW距离矩阵做K-medoids非常自然因为整个过程不需要计算任何平均轨迹只需要在已有的距离矩阵里做查询。得到聚类结果后每个簇的代表路径就是一条真实存在的GPS轨迹可以直接可视化直接给业务方解释。目标函数也明显变成了轨迹版的经典形式对所有簇C_k找medoidm_k最小化sum(D(i, m_k))其中i属于该簇。4.3 固定长度向量化KMeans的适用场景我不是说KMeans完全不能用。有一种场景可以凑合把轨迹重采样成固定点数然后展平成[lat1, lon1, lat2, lon2, ...]的二维向量丢进sklearn的KMeans。这个方案胜在速度快、依赖少原型验证时很有用。但使用前必须做投影转换把经纬度换算成米制坐标否则欧氏距离在南北方向和东西方向尺度不一致。另外它对轨迹速度动态变化仍然敏感同一条路正着走、反着走或者只是采样起点偏了一点展开后的向量就可能离得很远聚类结果不稳定。所以我把这种方案定位为“快速粗聚类”真正交付业务结论时还是用DTW配合K-medoids。5. Python实现从GPS点到聚类簇的完整代码5.1 数据结构约定与基础工具函数我在项目里把每条GPS轨迹统一成这样的结构一个二维NumPy数组第一列是纬度第二列是经度行数不固定。import numpy as np # 每条轨迹示例 traj_example np.array([ [31.2304, 121.4737], [31.2305, 121.4739], [31.2309, 121.4742], # ... [31.2456, 121.4890] ])底层的距离函数就用前面写的haversine。再补一个按累计距离重采样的函数def resample_trajectory(traj, num_points64): seg_lens np.zeros(len(traj) - 1) for i in range(len(traj) - 1): seg_lens[i] haversine(traj[i][0], traj[i][1], traj[i 1][0], traj[i 1][1]) cum_len np.concatenate([[0], np.cumsum(seg_lens)]) total_len cum_len[-1] target_dist np.linspace(0, total_len, num_points) resampled [] j 0 for dist in target_dist: while j len(cum_len) - 2 and cum_len[j 1] dist: j 1 seg_ratio (dist - cum_len[j]) / max(seg_lens[j], 1e-12) lat traj[j][0] (traj[j 1][0] - traj[j][0]) * seg_ratio lon traj[j][1] (traj[j 1][1] - traj[j][1]) * seg_ratio resampled.append([lat, lon]) return np.array(resampled)这个函数可以放在数据预处理阶段把每条轨迹都压缩成64个均匀分布的点。5.2 构建轨迹距离矩阵拿到所有重采样后的轨迹下一步就是构建一个n x n的距离矩阵矩阵第(i, j)个元素是轨迹i和轨迹j的DTW距离。def build_dtw_matrix(trajs): n len(trajs) D np.zeros((n, n)) for i in range(n): for j in range(i 1, n): d dtw_distance(trajs[i], trajs[j]) D[i, j] d D[j, i] d return D这里有个性能优化点矩阵是对称的所以只算上三角省一半时间。如果n超过1000建议再加并行比如使用多进程每个进程负责一部分行。5.3 用K-medoids迭代找聚类簇K-medoids的迭代逻辑很简单先随机选K条轨迹作为初始medoid然后重复“分配簇”和“更新medoid”两步。class TrajectoryKMedoids: def __init__(self, n_clusters3, max_iter100, random_state42): self.n_clusters n_clusters self.max_iter max_iter self.random_state random_state self.medoid_indices_ None self.labels_ None def fit(self, D): n D.shape[0] rng np.random.default_rng(self.random_state) medoids rng.choice(n, sizeself.n_clusters, replaceFalse) for _ in range(self.max_iter): labels np.argmin(D[:, medoids], axis1) new_medoids medoids.copy() for k in range(self.n_clusters): cluster_idx np.where(labels k)[0] sub_dist D[np.ix_(cluster_idx, cluster_idx)] cost sub_dist.sum(axis1) best_local cluster_idx[np.argmin(cost)] new_medoids[k] best_local if np.array_equal(new_medoids, medoids): break medoids new_medoids self.medoid_indices_ medoids self.labels_ np.argmin(D[:, medoids], axis1) return self注意我更新medoid时用的是“簇内到其他点总距离最小”的标准而不是平方距离最小。因为DTW距离本身不满足欧氏空间的代数性质用平方和会把个别长距离轨迹的影响放大。业务上更关心的是“代表性轨迹和绝大多数轨迹靠得近”因此用和为优化目标是合理的。5.4 输出聚类结果和典型路径聚类结束后把标签挂回原始轨迹DataFrame同时取出medoid对应的轨迹作为该簇的“典型路径”。我通常还会保存一个JSON里面存每条轨迹的ID、簇编号、以及到medoid的距离。这样业务方可以直接看到某类轨迹的代表路径是哪一条距离中位数是多少有没有离群异常样本。6. K值、轮廓系数和人肉检验6.1 先看簇内总距离的衰减曲线K值怎么定是所有聚类任务都躲不开的问题。这里不能用分类准确率只能靠“组内紧凑程度”。对K-medoids来说我定义每个簇的组内距离为所有样本到本簇medoid的距离之和整体指标就是所有簇的和。def total_within_cluster_distance(D, labels, medoids): total 0.0 for i in range(len(labels)): total D[i, medoids[labels[i]]] return total把K从2遍历到10画出组内总距离随K变化的折线。理论上曲线会有一个“拐点”在拐点之前增加K带来的距离下降非常快过了拐点后下降幅度变小。这个拐点可以作为K值参考。不过拐点法在实际轨迹数据里经常不清晰。我的习惯是结合作业目的来定比如业务方希望“区分出几条不同的主流路线”那么K等于主流路线数量最合理此时不会刻意追求极小化组内距离。6.2 轮廓系数的理解和局限另一个常用的验证方法是轮廓系数。对每条轨迹i如果和它同簇的其他轨迹之间的平均距离记为a(i)和它最近的那个不同簇中所有轨迹的平均距离记为b(i)那么轮廓系数就是s(i) (b(i) - a(i)) / max(a(i), b(i))s(i)接近1说明这条轨迹的同类聚合度高距离其他簇远接近-1说明它可能被分错了簇。所有样本取平均就是全局轮廓系数。但轨迹聚类里轮廓系数高不代表结果一定好。原因在于轨迹相似度本身可能“欺骗”指标比如大多数样本是短途轨迹短途之间距离天然就小轮廓系数会偏高。所以我每次都会把轮廓系数和可视化的dealed轨迹放在一起看先抽几个簇让业务方人眼判断再决定是否调K。6.3 可视化最容易骗人GPS轨迹聚类结果落到地图上时因为线路多、交叉方向多地图上看起来很容易显得很“乱”。用matplotlib直接画经纬度散点会把所有轨迹叠在一起颜色是靠标签区分的交叉路口处各色混在一起不好分辨。我一般是这样化的每个簇随机抽10到20条轨迹用同一颜色画在独立子图里再把medoid轨迹用粗线叠上去。这样既能看到该簇内部的离散程度也能看到代表轨迹的形态。簇与簇之间不要画在一张图上否则交叉点太多人眼根本看不过来。7. 真实跑聚类时遇到并解决的几个坑7.1 时间戳没排序导致轨迹倒着画有一批数据来自不同采集端有一个采集端传过来的轨迹时间字段排序异常。聚类前如果不按时间戳排序直接重采样折线会来回穿插空间形态完全失真聚类结果几乎没有可解释性。处理办法很笨但有效先把每条轨迹按时间戳排序过滤掉时间差为负的脏点再检查首尾点是否在目标区域范围内。这一步不会让结果变“聪明”但能保证后续计算不离谱。7.2 停留点不压缩DTW结果被污染有些轨迹在商圈附近停留很久一分钟内产生了上百个点。DTW会把“停留”识别成路径上的一个重要形态。结果就是两条真正同路的轨迹因为一条停过、一条没停距离被显著拉大。解决方案就是前面提到的停留点压缩。我设置规则连续5分钟内点群半径小于15米就把这群点压成首尾两个点。这个操作跑完之后聚类结果里“同一路线不同停留情况”的轨迹被稳定归到同一簇了。7.3 距离矩阵内存和计算量失控两千条轨迹不算多但两两DTW后的距离矩阵也有约200万个浮点数内存8MB左右看起来不大真正吓人的是计算时间。如果每条轨迹不重采样平均1000个点一次DTW就要计算100万次距离两千条轨迹两两匹配总计算量非常大。我实际采用的办法是先对轨迹做抽稀将每个轨迹控制在64个点再计算DTW矩阵时用并行worker。把两千条轨迹的完整矩阵控制在十几分钟内跑完。如果数据量到万条以上我会先抽样一部分样本来选K和参数再用全量数据跑最后一轮。7.4 空间相似度定义必须跟业务对齐最后也是最容易被忽略的坑技术上的“空间相似度”不等于业务上的“同类”。比如咖啡馆配送场景业务方认为只要送到同一栋楼就算同类不管你是从南门进还是北门进但直接用DTW聚类可能会把南北门两个方向分成两簇。这时候不能猛调算法参数而是要在预处理阶段把“只保留第一公里和最后两百米”作为特征区域或者把起终点单独加入距离函数。相似度算法只是工具定义清楚“哪部分空间特征重要”才是项目成功的前提。我自己跑过好几轮之后最大的体会是轨迹聚类里距离函数和质心选择比聚类算法本身更值得花时间。标准KMeans不是不能用但要建立在固定长度、投影坐标、欧氏距离的良好前提下一旦距离函数变成DTW这种非欧度量就别硬着头皮算平均质心换成K-medoids让真实的代表性轨迹来说话。如果你正在处理类似GPS轨迹归类的问题建议先把预处理和距离度量敲定再回头调K值你的聚类结果会比直接套模型稳定得多。