DCTOD双层聚类轨迹异常检测：停留点与轨迹段工程复现

简介《面向物联网移动终端设备的用户轨迹异常检测研究》是一篇面向物联网轨迹数据挖掘、移动终端用户行为分析与异常检测方向研究者的网络首发论文作者为李健、付雄、王俊昌刊载于《计算机应用研究》。内容针对海量且分布不均的移动轨迹数据提出双层聚类检测思路先在特定空间距离与时间间隔下提取停留点经层次聚类划分停留区域并识别异常停留区域再对停留区域间的运动轨迹段二次聚类以发现异常轨迹段。文中结合实验验证该方法在检测全面性与速度上的优势并讨论其用于用户行为分析、安全保障与个性化服务的实践价值。资源包为1个PDF文件约689KB结构完整、便于检索阅读已有295人学习下载适合相关方向的研究生、教师与工程技术人员参考。1. 轨迹异常检测最容易翻车的地方不是算法不够花哨做移动终端取证或者用户行为分析的人迟早会撞上同一个问题手里有几万个带时间戳的经纬度点想知道这个人哪天哪段路走得不对劲。大多数现成方案会把停留、移动、区域三件事塞进一个距离公式里算结果就是参数调到崩溃异常名单里一半是正常通勤。这篇《面向物联网移动终端设备的用户轨迹异常检测研究》给的思路是把问题拆成两层先把原始轨迹压缩成停留点对停留点做层次聚类从簇的分布里挑出异常停留区域再对停留区域之间的运动轨迹段做第二次层次聚类挑出异常轨迹段。作者团队把方法缩写为 DCTODdouble-layer clustering trajectory outlier detection双层聚类异常轨迹检测。它解决的不是识别异常有多准这种泛问题而是在终端设备算力有限、数据分布极不均匀的前提下怎么把异常检测的规模压下来。适合两类人做位置数据取证分析的、做用户行为建模需要清洗轨迹的。下面按我复现时的顺序把阈值、距离公式和踩过的坑一条条摊开。2. 从原始轨迹点到停留点三个阈值决定后面全部结果停留在这一层的判定直接决定后面所有聚类的输入质量。这一章要处理的是把(x, y, t)三元组序列压成带起止时间的四元组同时把 GPS 漂移点清掉。压得粗异常段被吞掉压得细停留点数量爆炸第二层聚类的距离矩阵直接吃满内存。2.1 数据清洗按瞬时速度剔除跳变点原文的做法很朴素按时间排序生成轨迹序列 G比较相邻两点的空间距离与时间差算出平均速度 v若 v 大于阈值 Vmax 就丢掉后一个点然后重新遍历一遍直到序列里所有相邻点对都满足速度约束。原文实验里 Vmax 取 1000单位是公里每小时。这个 1000 不是正常行驶速度上限而是一条 GPS 跳变的兜底线。真实驾车很难持续超过 200 km/h取值 1000 意味着只清掉那种一秒钟跨了几十公里的信号漂移。我一般会往下压到 300500取决于采集设备的采样频率采样间隔越短正常位移对应的瞬时速度越容易被噪声放大阈值反而要放宽采样间隔在 15 分钟量级的手机数据300 已经足够。实现上有两个必须自己补的细节。一是距离不能用经纬度差直接算得先做球面距离或者等距圆柱投影否则在高纬度地区误差会到 30% 以上。二是原始数据里时间戳单位不统一秒、毫秒、字符串都见过不做统一会在算时间差时得到 0 或者天文数字速度判定就彻底失效。import numpy as np import pandas as pd R 6371000.0 # 地球平均半径单位米 def haversine(lon1, lat1, lon2, lat2): 球面距离输入为十进制度输出为米 lon1, lat1, lon2, lat2 map(np.radians, [lon1, lat1, lon2, lat2]) dlon lon2 - lon1 dlat lat2 - lat1 a np.sin(dlat / 2.0) ** 2 np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2.0) ** 2 return 2 * R * np.arcsin(np.sqrt(a)) def clean_trajectory(df, vmax_kmh500.0, max_round10): df 需含列: lon, lat, ts(秒级时间戳) 返回按时间排序、剔除速度异常后的序列 g df.sort_values(ts).reset_index(dropTrue).copy() for _ in range(max_round): if len(g) 2: break d haversine(g[lon].values[:-1], g[lat].values[:-1], g[lon].values[1:], g[lat].values[1:]) dt np.diff(g[ts].values).astype(float) dt[dt 0] 1e-6 # 防除零同时暴露时间戳重复问题 v (d / dt) * 3.6 # m/s - km/h bad np.where(v vmax_kmh)[0] 1 # 标记后一个点 if len(bad) 0: break g g.drop(indexbad).reset_index(dropTrue) return gvmax_kmh控制清洗强度max_round是迭代上限——每剔除一轮就可能产生新的相邻关系理论上要循环到不动点实际数据里 23 轮就收敛了设 10 只是防死循环。代码里dt 0被压成极小值而不是直接删点是因为时间戳重复往往意味着采集端时钟回拨属于另一类问题我倾向于让它在这里以速度爆炸的形式暴露出来而不是静默吞掉。提示清洗后的点数一定要打印出来和原始条数对比。如果掉了 30% 以上先别急着往下走多半是时间戳单位搞错了而不是数据质量差。2.2 停留点生成一次扫描的四元组压缩原文对停留点的定义是一段连续轨迹段段内相邻点距离小于空间阈值 Wd且段的首尾时间差大于时间阈值 Wt。满足条件时这整段被压成一个四元组s (x, y, t_start, t_stop)其中 x、y 取段内所有点经纬度的算术平均t_start取段首点时间t_stop取段尾点时间。对应算法是单次线性扫描从轨迹序列头部开始只要当前点和下一个点的距离还在 Wd 内就继续并入当前段一旦超出就以该点为界重新开段。整个算法时间复杂度 O(n)空间复杂度 O(n)。原文的伪代码里还有一个隐含条件——段内点数必须大于 2也就是至少三个点才算停留点两个点的停留更可能是采样噪声。def extract_stay_points(g, wd100.0, wt1800.0, min_pts3): g: clean 之后的 DataFrame(lon, lat, ts) wd: 空间阈值米 wt: 时间阈值秒1800 30 分钟 min_pts: 停留点最少包含的原始点数 stay [] lon g[lon].values lat g[lat].values ts g[ts].values n len(g) i 0 while i n - 1: j i 1 # 向内扩张当前停留段 while j n and haversine(lon[j-1], lat[j-1], lon[j], lat[j]) wd: j 1 seg_len j - i duration ts[j-1] - ts[i] if seg_len min_pts and duration wt: stay.append({ lon: lon[i:j].mean(), lat: lat[i:j].mean(), t_start: ts[i], t_stop: ts[j-1], n_pts: seg_len, }) i j # 停留结束点作为下一段起点 else: i i 1 # 不满足则滑动一个点重试 return pd.DataFrame(stay)wd和wt是这套方法里唯一需要人工干预的两个参数原文实验取Wd 100米、Wt 30分钟。这两个数不是拍脑袋100 米大致对应一个建筑群或一个街区30 分钟大致对应坐下来吃饭、开会、办事这类真实停留行为。往下压到 50 米 / 10 分钟会把等红灯、堵车也算成停留往上升到 300 米 / 2 小时家和工作地这种真正的高频停留反而会被切开——因为家庭住址的 GPS 抖动半径在城市楼宇间经常超过 150 米。i j和i i 1这个分支值得单独说一句。前者表示这段确实是停留直接跳到段尾继续后者表示这段不够格只往后挪一个点。如果两个分支都写成i j一旦采样间隔大导致每个点都不满足 Wd整条轨迹会被一个点一个点地吃掉最后返回空列表。2.3 参数之间的相互牵制三个阈值不是独立的。Wd 调小会缩短每段的跨度间接让能通过 Wt 的段变少Wt 调小会让停留点数量上升直接把第二层聚类的规模推上去。原文自己也提到实验是建立在已有算法对时空阈值和轨迹距离权重的规定之上这意味着这套参数更像一个可用起点而不是最优解。复现时我一般按这个顺序定先看数据采样间隔间隔中位数决定 Wt 的下限Wt 至少要覆盖 510 个采样点否则永远凑不够 min_pts再根据区域建筑密度定 Wd城区 80120 米、郊区 150300 米最后核对输出的停留点数量和原始轨迹点数量的比例健康区间大概在 1:50 到 1:500 之间偏离这个量级就说明参数跑偏了。3. 第一层层次聚类把停留点折成停留区域停留点生成完数据就变成了一张用户常去哪儿的稀疏点图。第一层聚类要回答的问题是这些点里哪些是用户反复出现的地方形成停留区域哪些是孤零零冒出来一次的可疑停留。3.1 为什么原文绕开了 DBSCANDBSCAN 在这类任务里是最常见的默认选项但原文明确说了不用它理由有两个类与类之间距离过长时邻域半径 eps 和邻域内最小点数 minPts 极难确定而停留点的聚类往往是小数据量问题层次聚类的开销可以接受。这个判断我认同。DBSCAN 的核心假设是同一簇内密度相近但用户轨迹的密度分布是天然不均匀的——家和工作地可能各有一两千个点偶尔去的医院、餐厅只有十几个点。用同一组 eps/minPts 去套要么把稀疏的正常停留判成噪声要么把密集区域的边缘点吸进核心簇。层次聚类不需要预设簇数也不需要密度均匀代价是要自己决定在哪一层截断。原文用自底向上的凝聚式层次聚类。两个停留点之间的距离用欧式距离因为停留点坐标已经是段内均值可以近似当平面点处理两个停留点集合之间的距离用集合内所有点对距离的平均值也就是平均连接average linkage。3.2 凝聚过程与 O(n³) 的现实代价算法的推进方式很直接初始每个停留点自成一簇每轮找出距离最近的两个簇合并合并后重算距离矩阵直到所有点并成一簇。原文给出的复杂度是时间 O(n³)、空间 O(n³)这个量级在小样本下没问题但停留点一旦过千纯 Python 实现会直接卡死。实际复现建议用成熟库的距离矩阵接口把距离计算下沉到 C 层同时只保留压缩上三角矩阵condensed matrix空间能从 n² 降到 n(n-1)/2。原文实验用的用户 002 有 248217 个原始轨迹点压缩成停留点后规模已经小了一个数量级这是双层设计里最实在的一步优化先用阈值把点数砍掉再做复杂度高的聚类。3.3 轮廓系数定簇数R 最大不等于业务最优原文没有让人手工指定簇数而是每合并一次就算一遍轮廓系数 R最后取 R 最大的那次聚类结果。轮廓系数的定义是R (1/N) · Σ [(b(i) - a(i)) / max{a(i), b(i)}]其中 a(i) 是停留点 s_i 到同簇其他点的平均距离b(i) 是 s_i 到最近的其他簇中所有点的平均距离的最小值。R 落在 [-1, 1]越大说明簇内越紧凑、簇间越分离。import numpy as np from scipy.cluster.hierarchy import linkage, fcluster from scipy.spatial.distance import pdist, squareform def best_clustering(stay_df, metriceuclidean): 在凝聚层次聚类上扫描轮廓系数返回最优簇划分 # 停留点是经纬度先投影到米制平面避免度为单位造成的各向异性 lon np.radians(stay_df[lon].values) lat np.radians(stay_df[lat].values) x R * lon * np.cos(lat.mean()) y R * lat pts np.column_stack([x, y]) Z linkage(pts, methodaverage, metricmetric) best_k, best_r, best_lab 2, -2.0, None for k in range(2, min(len(pts) - 1, 30)): lab fcluster(Z, tk, criterionmaxclust) if len(np.unique(lab)) 2: continue r silhouette(pts, lab) if r best_r: best_k, best_r, best_lab k, r, lab return best_k, best_r, best_lab def silhouette(pts, labels): 标准轮廓系数距离用欧式 D squareform(pdist(pts)) n len(pts) s np.zeros(n) for i in range(n): same (labels labels[i]) same[i] False a D[i, same].mean() if same.any() else 0.0 b min(D[i, labels c].mean() for c in np.unique(labels) if c ! labels[i]) s[i] 0.0 if max(a, b) 0 else (b - a) / max(a, b) return s.mean()先把经纬度投影到米制平面这一步不能省。直接拿度做欧式距离经度方向 1 度对应的实际距离随纬度变化同一份数据在不同城市会得到不同的聚类结构这是典型的玄学复现失败来源。代码里用cos(平均纬度)做等距圆柱近似跨城市尺度的数据够用跨省数据建议换 UTM 分带。k的扫描范围卡在2到min(n-1, 30)上限 30 是我自己的经验值——用户停留区域的合理数量基本不会超过几十个扫到 30 已经能覆盖全部有意义的结构再往上纯属浪费。注意轮廓系数对非凸簇和密度差异大的数据是偏悲观的。原文用它选簇数选出来的结果里簇内紧凑的往往是家和工作地这种高频区域稀疏但真实的停留会被推到边缘。这正是原文下一步把最小簇直接视为异常停留点的前提也是它最大的争议点。原文最后一步是C Min(Q)在 R 最大的那次聚类结果里取样本数最少的那个簇作为异常停留点集。这个规则的隐含假设是用户极少去的地方就是异常在取证场景下成立但在个性化推荐场景下完全不成立——用户第一次去的城市会被判成异常。要落地到别的场景这里得换成相对频次低于阈值或者最近 N 天首次出现。4. 第二层聚类轨迹段的三元距离与异常段判定第一层给的是用户在哪些区域停留第二层要回答用户在区域之间怎么移动。后者才是异常轨迹的主体——大部分异常行为不体现为去了一个奇怪的地方而体现为用一条奇怪的路到达一个正常的地方。4.1 轨迹段的切分只对区域间段做二次聚类原文按停留点把整条轨迹切成段形成集合 P {s1,s2, s2,s3, ..., s_{n-1},s_n}每条段以两个停留点为端点中间夹着原始轨迹点。切完之后分两类停留点聚类簇内部两个停留点之间的段以及不同停留区域之间的段。只有第二类进入二次聚类。理由是簇内段本质上还是在同一个区域里转反映的是区域内不同时段的停留速度慢、距离短、方向随机做聚类只会把噪声放大。这个取舍很关键直接决定了异常检测的召回率把区域内段也塞进去正常的小幅移动会把距离分布整体拉平真正的异常段就被淹没了。4.2 三个距离分量速度、角度、位置原文的段间距离是三个分量的加权和d(p_i, p_j) w1·d_v(p_i, p_j) w2·d_θ(p_i, p_j) w3·d_L(p_i, p_j)速度距离 d_v 是两段平均速度之差的绝对值段平均速度等于段内所有相邻点距离之和除以总时长。角度距离 d_θ 基于偏转角对段内每个点 g_{i1}用 g_i、g_{i1}、g_{i2} 三点连线形成的夹角定义偏转角 β段的距离取两段偏转角差值的累加。位置距离 d_L 用 Hausdorff 距离即一个段内所有点到另一段最近点距离的最大值再对两个方向取 max。原文实验里三个权重都取 1/3。def seg_features(seg): seg: DataFrame(lon, lat, ts)返回该轨迹段的三个基础特征 d haversine(seg[lon].values[:-1], seg[lat].values[:-1], seg[lon].values[1:], seg[lat].values[1:]) dt np.diff(seg[ts].values).astype(float) dt[dt 0] 1e-6 v float(d.sum() / dt.sum()) # 平均速度 m/s ang [] lon, lat seg[lon].values, seg[lat].values for k in range(1, len(seg) - 1): # 相邻两点连线的方位角 b1 np.arctan2(lat[k] - lat[k-1], lon[k] - lon[k-1]) b2 np.arctan2(lat[k1] - lat[k], lon[k1] - lon[k]) diff abs(b2 - b1) ang.append(min(diff, 2 * np.pi - diff)) # 折叠到 [0, pi] mean_ang float(np.mean(ang)) if ang else 0.0 # 平均偏转角弧度 return v, mean_ang, seg[[lon, lat]].values def segment_distance(f1, f2, w(1/3, 1/3, 1/3), scale(3.0, np.pi, 300.0)): f1, f2: seg_features 的输出 scale: 各分量的归一化尺度 (速度 m/s, 角度 rad, 位置 m) v1, a1, p1 f1 v2, a2, p2 f2 d_v abs(v1 - v2) / scale[0] d_a abs(a1 - a2) / scale[1] d_p hausdorff(p1, p2) / scale[2] return w[0] * d_v w[1] * d_a w[2] * d_p这段代码里scale是我自己加的原文没有。原因见下一节。hausdorff用标准双向 Hausdorff两个方向各取 max 再取 max实现上直接对两个点集算全对距离矩阵行最小值取最大、列最小值取最大然后取二者较大值。4.3 权重 1/3 是个陷阱先归一化再谈权重原文三个权重都取 1/3看起来是三者同等重要但直接照抄会翻车。因为三个分量的量纲完全不同速度差可能是 030m/s角度差在 03.14弧度位置差在 0几千米。不归一化的情况下位置距离一项就能吃掉整个加权和的 99%速度和角度形同虚设。这解释了为什么有些复现出来的结果异常段全是最长的那些段——模型其实只看到了长度。我的做法是先按业务上多大算显著给每个分量定一个尺度归一化到 01 再加权。速度尺度取 3 m/s步行和骑车的分界角度尺度取 π完全反向位置尺度取 300 米一个街区。这套尺度在城区轨迹上表现稳定换到城际数据要相应放大。第二层聚类的过程和第一层完全一样用段间距离矩阵构建层次聚类扫描轮廓系数取最优划分最小簇即异常轨迹段。原文没有重复伪代码只说将算法 2 的输入改为轨迹段集、距离方法改为轨迹段间距离这个复用设计是这套方法里比较优雅的部分——同一个聚类框架换两次输入和距离函数就覆盖了两种粒度。要注意一个实现细节段间距离矩阵的构建成本是 O(m²·n)m 是段数、n 是段内平均点数Hausdorff 每一对都要算全对距离。停留区域多的用户段数可能上百这一步是整套流程里最慢的。我一般先把段内点做等间隔降采样到 20 个点以内Hausdorff 的结果几乎不变耗时能降一个数量级。5. 复现踩坑五个让 DCTOD 跑不出原文结果的问题下面五条都是我在复现同一份逻辑时真实遇到过的按现象、原因、解决三段写。5.1 停留点层三个把结果带偏的参数坑现象一停留点数量为 0或者比预期少一个数量级。原因通常是 Wd 用错了单位。原文公式里的距离 d 是两点间的欧式距离但经纬度直接相减得到的是度1 度约 111 公里拿 100 当阈值意味着要求两点距离小于 100 度等于永远满足反过来若先转成米又忘了纬度方向的经度缩放阈值会偏小。解决办法是统一用球面距离或等距投影后的米制坐标且把 Wd 的单位在函数签名里显式写出来我在extract_stay_points里直接命名为wd并在注释标明米就是为了避免这个。现象二停留点覆盖了整条轨迹几乎每个点都变成停留点。通常是 Wt 设得太小或者时间戳单位是毫秒却按秒传进去了。原始数据里 1 分钟停留变成 60000直接跨过 1800 的门槛。解决办法是在入口处做一次单位断言取时间差的中位数如果大于 10000 就认定是毫秒统一除以 1000。现象三坐标均值算出来落在停留段外面。原文伪代码里有个小混乱——sumLatitude累加的是gi.x而 x 定义是经度。虽然均值的数学形式不变但复现时按字面抄会把经纬度对调最后画出来的图横竖颠倒。我用 DataFrame 的lon/lat命名强制区分比在代码里靠注释保证更可靠。5.2 聚类层轮廓系数与量纲的两个隐雷现象四轮廓系数扫描结果全为负或者最优簇数永远是 2。原因有两个一是停留点少于 3 个时轮廓系数无定义代码里要显式跳过二是平均连接average linkage在样本高度不均衡时会倾向于把大簇拆成两个去提升平均分离度。解决办法是把扫描下限提到 3同时把最小簇样本数小于 2的划分直接过滤掉——一个点自成一簇时 a(i) 恒为 0轮廓系数会被虚高。现象五二次聚类把正常通勤段全判成异常。这是量纲问题前面已经说过。另一个容易忽略的原因是轨迹段数量太少区域间段只有十几条时层次聚类的第一层合并就相当于把一半的段划成一类最小簇必然异常。解决办法是给二次聚类设一个最小样本数门槛比如区域间段少于 10 条时直接跳过改用规则判据平均速度偏离全局中位数 3 倍以上。现象六能跑通但跑不完距离矩阵吃满内存。原文的 O(n³) 复杂度在小数据下无感停留点到几百个就开始卡。解决办法是用压缩距离矩阵替代全矩阵把 Hausdorff 换成降采样后的近似以及最有效的——在每个停留区域内部先做一次点级去重同一区域重复出现的停留点只保留代表性的一条。原文实验选的是单个用户的数据量级本来就不大这一点在原文里没有被强调但工程落地上绕不过去。6. 从单用户到批量跑我常用的验证与调参习惯原文实验用的是某公开 GPS 轨迹数据集182 名志愿者、2007 到 2012 年其中被选中的用户有 248217 个轨迹点实验环境是四核 2.5 GHz、8 GB 内存、JDK 1.8。单个用户跑通不代表整套流程可用真正要落地得先把验证方式定下来。我的做法是注入合成异常而不是等标注。具体操作在清洗后的轨迹里人为构造三类异常——把某一段的坐标整体平移 2 公里空间异常、把中间若干点的时间戳压缩到 1/10速度异常、把某段轨迹反向角度异常然后看二次聚类能不能把这三段捞出来。召回率低于 60% 就回头查距离归一化的尺度是不是设错了。这个方法的好处是不依赖任何人工标注而且能分别暴露三个分量各自的敏感度。参数敏感性上我整理过一张常用的起点表参数作用城区起点郊区起点调整信号Vmax清洗跳变点300 km/h500 km/h清洗后掉点超 10% 则放宽Wd停留空间阈值80100 m150300 m停留点占比低于 1:500 则放大Wt停留时间阈值30 min20 min停留点少于 30 个则缩小段数门槛二次聚类最小样本1010低于门槛改用规则判据速度尺度归一化基准3 m/s8 m/s异常段全是长途则放大调参顺序我固定为先定 Wd 让停留点数量和原始点数的比例落在健康区间再扫 Wt 让停留区域数量落在 520 个之间最后调二次聚类的尺度。反过来先调聚类参数会一直在错误的输入上打转。还有一个容易忽略的收尾动作把每一层的中间结果都落盘。停留点表、段表、两个距离矩阵、每次合并的轮廓系数序列全部存成独立文件。我吃过一次亏——某次结果和预期差很远怀疑是聚类写错了回头排查两天最后发现是停留点生成阶段把 i 的滑动写成了i j空转了一整条轨迹。从那以后我每次跑这套流程都会在停留点生成后强制走一遍点数守恒检查所有停留点包含的原始点数之和加上被清洗掉的点数必须等于原始总数。对不上就立刻停不往下走。这套双层结构真正的价值不在聚类算法本身而在于它把停在哪和怎么走拆成了两个可以独立验证的环节——哪一层出问题中间结果一看就知道。希望帮到你。本文还有配套的精品资源点击获取