简介这份资源是李健等人在《计算机应用研究》网络首发的论文PDF面向物联网、移动计算与轨迹数据挖掘方向的研究生、算法工程师与安全分析人员聚焦移动终端轨迹数据量大且分布不均时异常检测效率偏低的难题。论文提出双层聚类的用户轨迹异常检测思路先在特定空间距离与时间间隔下提取停留点集合经层次聚类划分停留区域并识别异常停留区域再对区域之间的运动轨迹段做二次层次聚类定位异常轨迹段实验显示该方法在检测全面性与速度上均优于传统算法。资源包共1个文件为约689KB的PDF单篇论文包含中英文摘要、关键词、引言、方法设计与实验分析等完整章节可直接用于论文精读、方法复现与研究综述引用。页面显示已有296人学习下载适合作为轨迹异常检测方向入门与方案对照的参考材料。1. 从一份网络首发 PDF 说起双层聚类到底能解决什么手里攒了几十万条从移动终端导出的定位记录想知道哪一段不对劲——这个需求听起来简单真动手第一步就卡住点太多。这份网络首发的 PDF 论文《面向物联网移动终端设备的用户轨迹异常检测研究》给了一套能照着搭的骨架先在时空双阈值下把稠密点压成停留点再对停留点做一次层次聚类找出异常停留区域最后对停留区域之间的运动轨迹段做第二次聚类找出异常轨迹段整条链路叫双层聚类DCTOD。它面向的是物联网移动终端设备上批量、分布不均的定位数据把检测目标拆成「停得不对」和「走得不对」两类。适合做位置数据取证、行为分析以及要给用户轨迹异常检测搭 baseline 的工程同学。2. 从原始定位点到停留点时空双阈值怎么卡2.1 为什么先做停留点而不是直接对点聚类论文里那组实验数据是单个用户的 248217 个 GPS 轨迹点。如果直接把这堆点丢进层次聚类簇间距离要反复算复杂度按点数三次方增长几十万点的量级在单机上基本没有跑通的可能。停留点的作用就是在这个环节把数据量砍下来一段连续的、空间上兜得住的轨迹被压成一个四元组(x, y, t_start, t_stop)其中 x、y 是段内所有点经纬度的均值t_start 取段首点时间t_stop 取段尾点时间。压缩比不是重点语义才是。睡觉、上班、在某个商场待一下午这些行为在数据上的表现就是长时间小范围停留压成停留点之后一个点就代表一次「有意义的行为」后面的聚类才有解释力。原始点序列和停留点序列的差别大致是这样维度原始轨迹点停留点字段(经度, 纬度, 时间戳)(经度均值, 纬度均值, t_start, t_stop, 点数)量级十万到百万几十到几百表达能力逐点位置一次停留行为下游用途清洗、切段聚类、异常判定停留点的判定条件在原文定义 2 里写得很清楚段内相邻轨迹点距离必须小于空间阈值 Wd段的结束时间与开始时间之差必须大于时间阈值 Wt同时段内点数超过 2。注意这里是「相邻点之间」的约束不是「每个点到段首点的距离」这意味着一段停留可以缓慢漂移成一条细长的带状区域。这在实车上很常见人坐在车里等红灯GPS 抖动会带着坐标一路挪真正判断一个停留段是否合格最好再补一条「段首点到段尾点距离 2~3 倍 Wd」的兜底条件否则聚类时会冒出一些形状诡异的簇。2.2 数据清洗速度阈值先剔掉瞬移点清洗这一步原文的思路很朴素按时间排序得到轨迹序列算相邻两点的平均速度超过 Vmax 就丢掉后一个点反复扫描直到没有可丢的点。真正动手时有两个细节必须先处理掉否则后面全是错的。第一经纬度不能直接相减当距离。原文公式 (2) 写的是(x_i - x_j)² (y_i - y_j)²开根号这是平面简化的写法落到代码里必须先把经纬度投影成米。第二速度单位。距离除以时间得到的是 m/s而阈值习惯按 km/h 给中间差一个 3.6 的换算系数。import math R 6371000.0 # 地球平均半径单位米 def to_local_xy(lon, lat, lon0, lat0): 把经纬度投影到以 (lon0, lat0) 为原点的局部平面单位米。 几十公里范围内等距圆柱近似足够用别拿经纬度差值直接当距离。 x math.radians(lon - lon0) * R * math.cos(math.radians(lat0)) y math.radians(lat - lat0) * R return x, y def clean_by_speed(points, vmax_kmh1000.0): points: [(lon, lat, ts)]ts 为秒级时间戳输入先按时间升序。 相邻点速度超过 vmax 就丢弃后一个点反复扫到序列稳定为止。 pts sorted(points, keylambda p: p[2]) changed True while changed and len(pts) 2: changed False keep [pts[0]] for p in pts[1:]: lon0, lat0, t0 keep[-1] x, y to_local_xy(p[0], p[1], lon0, lat0) dist math.hypot(x, y) dt p[2] - t0 if dt 0: # 同秒多点或时钟回拨直接丢 changed True continue if dist / dt * 3.6 vmax_kmh: # m/s 换算成 km/h 再比 changed True continue keep.append(p) pts keep return pts逻辑上是一个「删除后重扫」的循环因为丢掉一个点之后它前后两个点会变成新的相邻关系新的速度可能又超标所以要反复扫到稳定。实践中 1 到 3 轮就收敛时间复杂度最坏 O(n²)实际接近 O(n)空间 O(n)。参数vmax_kmh原文实验取 1000这个值非常宽松基本只能拦住 GPS 漂移造成的瞬移点。换成步行、骑行数据一般收到 30 到 60 之间比较合理车载数据可以放到 120 到 200。阈值定太紧的后果是正常数据被大量删掉肉眼看不出来但停留点的时长会整体缩水。2.3 停留点生成双指针写法与三个阈值原文算法 1 的伪代码有一个很坑的地方for循环里i内层while里又写了i直接照抄会出现跳点甚至死循环。我按双指针重写了一遍结构更清楚也方便加边界判断。def extract_stay_points(pts, wd100.0, wt1800.0, min_pts3): pts: 清洗后的 [(lon, lat, ts)]按时间升序。 wd: 空间阈值米wt: 时间阈值秒min_pts: 段内最少点数。 返回 [(lon_mean, lat_mean, t_start, t_stop, n)]。 stays [] i, n 0, len(pts) while i n: j i lon_sum, lat_sum pts[i][0], pts[i][1] cnt 1 while j 1 n: lon0, lat0, _ pts[j] x, y to_local_xy(pts[j 1][0], pts[j 1][1], lon0, lat0) if math.hypot(x, y) wd: # 相邻点跳出空间阈值段到此截断 break j 1 lon_sum pts[j][0] lat_sum pts[j][1] cnt 1 t_start, t_stop pts[i][2], pts[j][2] if cnt min_pts and (t_stop - t_start) wt: stays.append((lon_sum / cnt, lat_sum / cnt, t_start, t_stop, cnt)) i j 1 # 关键从断点后一个点重新开段不是 i 1 return staysi是段首指针j一路往前探直到相邻两点距离超过 Wd 才停。跳出内层循环后i j 1这个赋值是整段代码里最容易写错的地方——写成i 1会让同一个点被重复计入多个停留段停留点数量翻倍且互相重叠后面的聚类轮廓系数会一路往下掉。整个算法对轨迹序列只做一次线性遍历时间复杂度和空间复杂度都是 O(n)。三个参数的取值和影响整理如下原文实验用的是 Wd 100 m、Wt 30 min这个组合对应的是人的日常活动尺度比如在办公室待半天、在家过夜参数含义原文取值调大之后Wd相邻点空间阈值100 m相邻点更容易连成一段停留点更少更粗细碎停留被合并Wt停留时长下限1800 s30 min短暂驻足被过滤适合看长期行为模式min_pts段内最少点数3即 k 2高采样率数据可以设大一点稀疏采样必须设小Vmax清洗速度上限1000 km/h越大越宽松越大越可能漏掉真正的跳变点min_pts是最容易被忽略的一个。它隐含了采样间隔的假设如果终端每 5 分钟上报一个点30 分钟也就 6 个点min_pts 3勉强能出结果如果采样间隔拉到 15 分钟30 分钟只有 2 个点这时候所有停留段都会被cnt min_pts卡掉返回的列表是空的。遇到「一个停留点都没有」的情况第一反应应该是去数一下原数据的采样间隔而不是怀疑聚类。3. 第一层聚类层次聚类加轮廓系数定簇数3.1 为什么不用 DBSCAN原文给的理由很直接停留点在空间上分布不均匀用户常去的几个地方会形成密集簇偶尔去一次的偏远地点孤立得厉害而 DBSCAN 的邻域半径 eps 和最少点数 minPts 在这种「类间距离拉得很长」的场景下非常难定。eps 调大密集区域会被糊成一团eps 调小孤立点全被判成噪声。层次聚类的好处是不需要预设簇数跑完得到一棵聚类树然后拿轮廓系数在树上切一刀选评估最好的那一层。代价是复杂度高但这里已经做过停留点压缩输入规模从几十万降到几百O(n³) 在停留点这一层是完全可接受的——这正是双层结构里「第一层压缩、第二层聚类」的设计意图。维度DBSCAN层次聚类需要预设eps、minPts 两个参数无簇数由轮廓系数事后选对不均匀分布敏感参数难兼顾不敏感只依赖集合间距离输出结构扁平簇 噪声点完整聚类树可观察合并过程时间复杂度依赖索引接近 O(n log n)O(n³)适用规模大样本几百到几千的小样本3.2 集合间平均距离与轮廓系数的实现停留点之间的基础距离用欧氏距离但注意这里要用投影后的米制坐标不是经纬度差值。两个停留点集合之间的距离用平均连接average linkage计算d(Ci, Cj) 1 / (|Ci| × |Cj|) × ΣΣ d(s, t)也就是把两个集合里所有点对的距离加起来除以点对数。原文公式 (4) 表达的就是这个意思。轮廓系数的定义是每个点算一个 s(i) (b(i) - a(i)) / max{a(i), b(i)}其中 a(i) 是该点到同簇其他点的平均距离衡量簇内紧密度b(i) 是该点到最近的异簇所有点的平均距离衡量分离度。全部点取平均就是这一层聚类的 R 值范围 [-1, 1]越接近 1 越好。import numpy as np from itertools import combinations def pair_dist_matrix(xy): xy: (n,2) 米制坐标。返回 n x n 欧氏距离矩阵向量化算完。 diff xy[:, None, :] - xy[None, :, :] return np.sqrt((diff ** 2).sum(-1)) def avg_linkage_history(dist): 自底向上层次聚类平均连接返回每一层的簇标签。 用点索引数组表示簇直接取子矩阵求均值避免反复更新中心点。 n dist.shape[0] clusters [np.array([i]) for i in range(n)] history [list(range(n))] # 第 0 层每个点独立成簇 while len(clusters) 1: best, bi, bj None, -1, -1 for a, b in combinations(range(len(clusters)), 2): sub dist[np.ix_(clusters[a], clusters[b])] d sub.mean() # 平均连接 if best is None or d best - 1e-12: # 加 eps 保证并列时结果稳定 best, bi, bj d, a, b merged np.concatenate([clusters[bi], clusters[bj]]) clusters [c for k, c in enumerate(clusters) if k not in (bi, bj)] [merged] labels np.zeros(n, dtypeint) for lab, c in enumerate(clusters): labels[c] lab history.append(labels.tolist()) return history def silhouette(dist, labels): 平均轮廓系数。labels 长度与 dist 维度一致。 labels np.asarray(labels) uniq np.unique(labels) if len(uniq) 2: return -1.0 s [] for i in range(dist.shape[0]): same labels labels[i] same[i] False if same.sum() 0: s.append(0.0) continue a dist[i, same].mean() b min(dist[i, labels u].mean() for u in uniq if u ! labels[i]) s.append((b - a) / max(a, b) if max(a, b) 0 else 0.0) return float(np.mean(s))avg_linkage_history每一轮扫所有簇对找最小平均距离合并之后把新标签写进history于是「每一层的聚类结果」都被留了下来后面按层算轮廓系数取 R 最大那层的标签作为最终结果。pair_dist_matrix用广播做向量化几百个点的距离矩阵一次算完。有个实现细节值得说d best - 1e-12里那个 eps 是为了避免浮点并列。簇对距离相等时不同的遍历顺序会合出不同的树同一份数据两次运行结果不一致排查起来非常折磨人。加上这个容差之后至少在同一台机器上是可复现的。想省事也可以直接用scipy.cluster.hierarchy.linkage(xy, methodaverage)拿到连接矩阵再用fcluster切层速度更快。但轮廓系数还是得自己按公式算因为它依赖的是原始距离矩阵不是 scipy 返回的压缩形式。3.3 从聚类结果里挑出离群停留点原文算法 2 的最后一步写的是「C Min(Q) 获取到离群停留点集」这句话在复现时基本没法直接用Q 是轮廓系数最高那层的聚类结果Min 指的是哪一维度上的最小是簇内点数最少还是离最近簇最远原文没有展开这块得补一条判定规则。我一般的做法是两条同时成立才算离群该簇的停留点数量占全部停留点的比例低于阈值说明出现频次低并且该簇质心到其他簇质心的最近距离大于所有簇间最近距离的中位数说明位置也偏。阈值不写死成米用中位数做自适应避免换了数据集就跑不动。def outlier_clusters(xy, labels, min_ratio0.05): 挑离群簇样本占比过小 与最近簇质心距离偏大。 返回 (离群簇标签列表, 各簇质心, 各簇到最近簇的距离)。 labels np.asarray(labels) uniq np.unique(labels) centroids np.array([xy[labels u].mean(0) for u in uniq]) sizes np.array([(labels u).sum() for u in uniq], dtypefloat) ratio sizes / sizes.sum() cdist np.sqrt(((centroids[:, None, :] - centroids[None, :, :]) ** 2).sum(-1)) np.fill_diagonal(cdist, np.inf) nearest cdist.min(1) thr np.median(nearest) # 自适应阈值不硬编码距离 out [u for k, u in enumerate(uniq) if ratio[k] min_ratio and nearest[k] thr] return out, centroids, nearestmin_ratio是频次门槛取决于样本总量停留点只有几十个的时候5% 意味着少于 2 个点的簇会被判为稀疏样本上百时这个比例可以放到 3% 到 8% 之间慢慢试。nearest取中位数当阈值样本里的正常活动区域大多集中在城市范围内中位数天然把「城区内正常距离」和「远郊孤点」分开。输出直接就是异常停留区域的簇标签回填到地图上能一眼看出是哪些地方。4. 第二层聚类轨迹段距离的三要素加权4.1 轨迹段怎么切哪些段不参与聚类停留点序列按时间排好之后相邻两个停留点之间夹着的那段原始轨迹就是一条轨迹段整个集合记为 P {s1,s2, s2,s3, ..., s_{n-1},s_n}。切完之后段分成两类两端停留点落在同一个停留区域内的叫簇内段落在不同区域的叫簇间段。原文明确只对簇间段做异常检测理由写得很简短实际操作中确实如此。簇内段是用户在同一个地方不同时段之间的位移空间跨度小速度和方向都被 GPS 抖动主导区分度极低硬拿去聚类只会产出一堆假阳性。簇间段才是真正描述「从 A 地到 B 地怎么走」的运动模式绕路、异常换道、非典型路线都藏在这一层里。注意这里段的划分方式和经典的「按固定长度切分」不一样段与段之间不一定共享端点长度也不一致。正因为如此后面计算段间距离时不能用逐点对应的欧氏距离只能用集合对集合的度量也就是 Hausdorff 那套。4.2 速度、角度、位置三段距离的落地原文公式 (6) 把段间距离拆成三个分量加权速度差、角度差、位置差。三个分量的定义分别是速度距离取两条段平均速度之差的绝对值。段内平均速度 段内相邻点距离累加 / 段首到段尾的时间差单位 m/s。角度距离取两条段平均偏转角之差的绝对值。偏转角由连续三个轨迹点形成的夹角决定段内所有中间点各算一个取平均。位置距离用 Hausdorff 距离定义是 max{ max_{p∈P} min_{q∈Q} d(p,q), max_{q∈Q} min_{p∈P} d(p,q) }也就是两条段点集之间「最不匹配的那一对点」的距离。def seg_velocity(seg): 段的平均速度单位 m/s。seg: [(x, y, ts)]x/y 已是米制坐标。 dist sum(math.hypot(seg[k1][0]-seg[k][0], seg[k1][1]-seg[k][1]) for k in range(len(seg)-1)) dt seg[-1][2] - seg[0][2] return dist / dt if dt 0 else 0.0 def seg_angles(seg): 段内各中间点的偏转角弧度。前后向量长度接近 0 的点跳过。 ang [] for k in range(1, len(seg) - 1): ax, ay seg[k][0]-seg[k-1][0], seg[k][1]-seg[k-1][1] bx, by seg[k1][0]-seg[k][0], seg[k1][1]-seg[k][1] na, nb math.hypot(ax, ay), math.hypot(bx, by) if na 1e-6 or nb 1e-6: # 原地停留的点跳过别产生 0/0 continue cos max(-1.0, min(1.0, (ax*bx ay*by) / (na*nb))) ang.append(math.acos(cos)) return ang def hausdorff(p, q): 两条段点集之间的 Hausdorff 距离米双向取最大。 def one_way(a, b): return max(min(math.hypot(pi[0]-qi[0], pi[1]-qi[1]) for qi in b) for pi in a) return max(one_way(p, q), one_way(q, p)) def seg_distance(p, q, w(1/3, 1/3, 1/3), scale(10.0, math.pi, 500.0)): 三要素加权段距离。scale 把三个量纲拉到同一尺度 速度按 10 m/s、角度按 π、位置按 500 m 各算一个单位。 dv abs(seg_velocity(p) - seg_velocity(q)) / scale[0] ap, aq seg_angles(p), seg_angles(q) dth abs((sum(ap)/len(ap) if ap else 0.0) - (sum(aq)/len(aq) if aq else 0.0)) / scale[1] dl hausdorff(p, q) / scale[2] return w[0]*dv w[1]*dth w[2]*dl这里最需要强调的是scale那一步。原文公式 (6) 是三个距离直接乘权重相加实验里权重取 1/3。问题在于速度距离的量纲是 m/s、角度距离是弧度、位置距离是米数值上位置距离动辄几百速度和角度基本在个位数和 3.14 以内等权相加等于让位置距离单方面主导速度和角度几乎不起作用——那不如干脆只算 Hausdorff。加上归一化之后三个分量才真正在同一个尺度上竞争。scale的取值跟数据有关(10.0, math.pi, 500.0)是我在城区出行数据上常用的起点速度 10 m/s 差不多是电动车到低速汽车的区间500 m 对应一次跨街区移动。换成高速公路场景位置尺度要放到 5 km 以上否则所有段的 Hausdorff 都接近 1聚类退化。4.3 二次聚类与两层结果的对应关系第二层直接复用第 3 章的avg_linkage_history和silhouette输入从停留点的距离矩阵换成段的距离矩阵流程一模一样算段间距离矩阵、自底向上合并、每层算轮廓系数、取 R 最大那层、把占比小且距离远的簇标成异常轨迹段。两层的输入输出对应关系整理成表更清楚层级输入距离定义输出业务含义第一层停留点序列点间欧氏距离平均连接离群停留点簇去了不该去的地方第二层簇间轨迹段速度/角度/Hausdorff 加权离群轨迹段走了一条不寻常的路第二层的距离矩阵规模是段数的平方段数通常比停留点还少n 个停留点最多 n-1 条段所以这一层的计算开销比第一层更小。真正耗时的是 Hausdorff它对每条段是 O(|p| × |q|)段长几十个点无所谓但段长上千点时要把累加改成 numpy 广播否则 Python 的双层生成器能把 CPU 吃满。拿到异常轨迹段的标签之后最好再回到原始点序列去看一遍这几段的起止时间和地理位置确认它不是由数据缺失造成的「断头段」。论文里只做了运行时间的横向对比没有给异常段的语义解释这块得靠自己的数据来补。5. 避坑与排查复现这篇论文最容易翻车的五处5.1 数据与阈值层面现象停留点数量要么等于原始点数的一半要么一个都没有。原因几乎都是距离单位没统一。wd 100如果被当成「度」来比地球上任意两点的经纬度差值都远小于 100整个序列会被连成一段反过来如果没能把经纬度差值转成米而代码里又用1e-3之类的小阈值去卡就会一个段都成不了。解决方式是把投影函数做成唯一入口所有距离计算前必须先调to_local_xy代码里不允许出现「直接用经纬度差值比阈值」的分支。我当时图省事在清洗和停留点两处各写了一份距离函数后来改阈值只改了一处两处结果对不上查了半天。现象清洗后点数减了一半正常轨迹被删得七零八落。原因是 Vmax 的单位搞混了。dist / dt出来是 m/s如果你直接把 Vmax 设成 60 去比等价于把上限压到了 60 m/s还勉强设成 5 就变成 5 m/s等于把骑车的人全判成跳变。解决方式是把换算写在比较表达式里* 3.6 vmax_kmh并且在清洗前后各打印一次点数看删除比例是不是符合预期。正常情况下干净数据的速度过滤删掉的比例应该在 1% 以内删到 10% 以上就是阈值或者时间戳出了问题。现象轮廓系数曲线一路走低最好的那层也才 0.1 出头。这个时候先别怀疑算法去看停留点的分布图。如果停留点全都集中在一个很小的区域内簇内紧密度和簇间分离度都算不出差异轮廓系数自然上不去。原因通常是 Wd 设得太大把相邻的两个活动区域合并成了一个。把 Wd 从 100 m 降到 50 m 重跑往往能看到曲线恢复正常的「先升后降」形状。5.2 聚类与距离计算层面现象凭经验设的异常判定规则检出的全是出差、旅游这类偶发远行点。这是单维度位置密度聚类的通病。用户一年去两次外地在空间分布上就是典型的小样本远簇一定会被第一层挑出来但它对行为分析的价值很低。解决办法是加最小支持度把停留点按时间切桶比如按周统计每个簇在多少个桶里出现过只出现一两次的簇先打上「低频」标记而不是直接判异常同时把工作日和周末分开跑两套结果取交集能过滤掉大量节假日出行造成的噪声。现象一跑就吃满内存或者段距离里冒出 NaN。内存问题基本是把原始点直接喂进了距离矩阵24 万点的方阵谁也存不下。检查输入规模第一层的 n 应该是停留点数而不是原始点数几百以内才合理。NaN 有两个来源一是seg_angles里遇到重复点算出 0/0代码里已经用na 1e-6跳过了二是段长为 1 或 2 时没有中间点seg_angles返回空列表平均值退化成 0如果另一条段算出来是 0 而max(a,b)也为 0除法同样会炸。解决方式是在切段时加一条硬过滤长度小于 3 的段直接不参与第二层聚类它们本来就没什么形状可言。6. 进阶参数扫描、注入验证和一个顺手的小工具参数定不下来的时候别凭感觉试。我一般写个网格扫描把 Wd 和 Wt 两两组合跑一遍看停留点数量和平均段内点数怎么变找曲线拐点。import itertools def sweep(pts, wd_list, wt_list): 扫 Wd / Wt输出每个组合的停留点数量与平均每段点数用来选参数。 rows [] for wd, wt in itertools.product(wd_list, wt_list): stays extract_stay_points(pts, wdwd, wtwt) rows.append({ wd: wd, wt: wt, n_stay: len(stays), pts_per_stay: sum(s[4] for s in stays) / max(1, len(stays)), }) return rows扫描的重点看n_stay随 Wd 增大而下降的斜率斜率最陡的那一段说明参数正在粗暴合并真实停留拐点之后的一档就是比较稳的取值。pts_per_stay太低比如低于 5说明段太碎后续聚类没意义太高上百说明 Wd 太大把整天的活动糊成了一段。Wt 主要决定短时驻足要不要保留做长期行为分析就设大做实时提醒就设小。验证方式上原文只比了运行时间没有给出检测效果的量化指标这是一处空白。想自己补的话最省事的做法是注入异常在一份正常轨迹里人工插入两类样本一类是绕行在两点之间插一条明显偏离的路径一类是异地长时间停留然后看两层聚类各自能不能把它们捞出来算召回率和误报数。注入的异常最好按比例插入比如 5% 的段这样不至于把聚类结构带偏。这套流程跑通之后再拿真实数据看结果心里才有底。最后一个顺手的小习惯把每次实验的 Wd、Wt、Vmax、权重、scale 和结果统计一起写进一个 CSV文件名带上判断版本。调到第十几次组合的时候你一定会记不清哪组参数出来的结果长什么样。我从被一次参数对齐搞混的经历之后每次调整阈值都强制把参数和输出的停留点数、异常簇数一起落盘哪怕只是临时跑一下。希望帮到你。本文还有配套的精品资源点击获取
