Multi-DBSCAN多密度聚类算法解析与Python实践

做空间聚类的人大概都遇到过这种尴尬数据里有几坨挤得很紧的点簇又有几片散得跟芝麻似的稀疏簇。经典的 DBSCAN 只给一套eps和min_samples半径定小了稀疏簇被当成噪声半径定大了密集簇又粘成一团。我最近在几个混合密度的空间轨迹数据集上反复折腾最终把 Multi-DBSCAN 在 Python 里完整跑通并且顺手把它用到了传感器时间序列的行为分段上。这篇就讲清楚它的核心思路、可直接复用的实现代码以及我调参时踩过的几个坑适合正在做一般聚类、时序聚类的朋友参考。1. 为什么固定 eps 搞不定多密度场景的核心矛盾1.1 经典 DBSCAN 的隐含假设DBSCAN 的确是个好算法它对簇的形状不敏感也不需要预先指定簇数。但它有个隐含假设数据里各个簇的密度是“同一个量级”的。eps决定了一个点的邻域半径min_samples决定核心点要求邻域内最少有多少个点。这两个参数本质上是在定义一个密度阈值密度低于这个阈值的区域一律归为噪声。这在很多玩具数据集上没问题但真实数据很少这么听话。比如一份城市出行轨迹数据CBD 区域的点密到每秒好几个到了郊区可能几分钟才有一个。用一个固定阈值去覆盖这两种区域结果一定是顾此失彼。1.2 多密度数据的三类典型表现我在实际项目里总结出三种最常见的失败形式稀疏簇被误杀eps按密集簇的尺度来定稀疏簇内部点间距离远大于eps里面的点全部变成噪声点。密集簇被错误合并eps按稀疏簇的尺度来定密集簇之间原本有明显空隙但因为半径太大边界被互相穿透几个簇糊成一个。同一簇被拆成碎片有些簇密度分布不均匀内部也有稀密变化固定半径会把边缘部分切出去形成大量碎小的伪簇。这三种情况我用过一个很直观的说法DBSCAN 的eps是一把固定刻度的尺子可数据里的物体大小不一样。Multi-DBSCAN 的思路就是给每个点发一把适合自己的尺子让它只和“看得顺眼”的邻居抱团。注意标准库里的sklearn.cluster.DBSCAN只支持全局单一eps如果数据是从多个密度来源混合出来的一定要先在eps上做文章而不是盲目调参强行“碰运气”。2. Multi-DBSCAN 的关键思路让每个点有自己的邻域半径2.1 局部密度估计k 近邻距离的合理用法Multi-DBSCAN 的核心是给每个点计算一个局部半径eps_i。最常用的做法是看它的 k 近邻距离对每个点找到离它最近的 k 个点统计这 k 个距离的中位数或均值这个值就代表了这个点周围的“尺度”。如果一个点周围很密集它的 k 近邻距离就小eps_i就小如果周围很稀疏k 近邻距离就大eps_i就大。这样一来不管是密集簇里的一个点还是稀疏簇里的一个点都能拿到和自己密度匹配的半径。实际操作里我一般取 k min_samples - 1这和 DBSCAN 里核心点判定的逻辑保持一致性。如果你用min_samples5那 k 取 4 到 8 都比较稳取中位数比取均值更抗离群点干扰。2.2 点对之间的“互认邻居”规则有了每个点的局部半径接下来要想清楚一个问题点p的邻居是dist(p, q) eps_p还是dist(p, q) eps_q这两个条件都不对称可能出现p认为q是邻居、q却不认为p是邻居的情况。聚类最终产生的是无向的簇结构所以必须让邻居关系对称。我用的是互认半径eps_pq (eps_p eps_q) / 2 dist(p, q) eps_pq 时p 和 q 互为邻居这个规则比 max 温和比 min 宽松吃掉了“高密度点向低密度点伸一伸、低密度点也够一够”的中间地带。实测下来它能让密度有明显梯度变化的区域平滑过渡而不是出现一条生硬的切割线。2.3 两种工程路线多阈值合并 vs 每点自适应实现 Multi-DBSCAN 大体有两类路线路线做法优点缺点多阈值合并从 k-距离曲线里挑出多个显著拐点分别用几个eps跑 DBSCAN再把结果按重叠关系合并逻辑简单可以利用成熟的 DBSCAN 实现拐点数量不好定合并规则写起来麻烦每点自适应每个点用自己的局部半径做邻域直接走一遍类似 DBSCAN 的扩张流程不需要手动挑阈值效果更细腻无法直接用现成库得自己写核心循环我重点用的是第二种因为它更接近“多密度自适应”这个本质。第一种路线我把它作为辅助手段用来观察数据里大致有哪几个密度层级顺带验证参数合理性。3. Python 完整实现从 k-距离曲线到自适应聚类3.1 数据准备与 k 参数选择先准备一份模拟数据方便验证效果。这里生成三个高斯簇密度明显不一样import numpy as np rng np.random.default_rng(42) g1 rng.normal(0, 0.3, (300, 2)) # 很密集 g2 rng.normal(3, 1.0, (200, 2)) # 很稀疏 g3 rng.normal([6, 3], 0.15, (150, 2))# 极密集 X np.vstack([g1, g2, g3])在跑 Multi-DBSCAN 之前我先画 k-距离曲线。这个曲线能告诉我数据里面密度差别有多大也能帮我定eps_scale的初始值。from sklearn.neighbors import NearestNeighbors def k_distance_curve(X, k8): nn NearestNeighbors(n_neighborsk).fit(X) dist, _ nn.kneighbors(X) # 每个点到第 k 个近邻的距离排个序 k_dist np.sort(dist[:, -1]) return k_dist k_dist k_distance_curve(X, k8)如果有 R 型曲线呈现多个明显“膝盖”也就是斜率突然变大的位置基本可以判断数据是多密度的。这种情况下固定eps的 DBSCAN 一定不会省心。3.2 核心代码每点自适应半径的 Multi-DBSCAN下面是完整实现可以直接复制用。代码分两个输入模式普通特征矩阵或者预先算好的两两距离矩阵。import numpy as np from sklearn.neighbors import NearestNeighbors from sklearn.metrics import pairwise_distances def multi_dbscan(X, k8, min_samples5, eps_scale1.0, precomputedFalse, clip_quantile(0.05, 0.95)): n X.shape[0] # 第一步估计每个点的局部半径 if precomputed: kth_dist np.partition(X, k 1, axis1)[:, 1:k 1] else: nn NearestNeighbors(n_neighborsk 1).fit(X) dist, _ nn.kneighbors(X) kth_dist dist[:, 1:k 1] # 去掉自身那列 local_eps np.median(kth_dist, axis1) * eps_scale lo, hi np.percentile(local_eps, [clip_quantile[0] * 100, clip_quantile[1] * 100]) local_eps np.clip(local_eps, lo, hi) # 第二步计算两两距离大数据请见后面的优化说明 # 注意precomputed 模式下 X 本身就是距离矩阵不能再用 pairwise_distances D X if precomputed else pairwise_distances(X) labels np.full(n, -1, dtypeint) visited np.zeros(n, dtypebool) cluster 0 # 第三步扩张聚类和 DBSCAN 的流程一致 for i in range(n): if visited[i]: continue visited[i] True neighbors np.where(D[i] (local_eps[i] local_eps) / 2)[0] neighbors neighbors[neighbors ! i] if len(neighbors) min_samples: labels[i] -1 # 暂时标为噪声 continue labels[i] cluster stack list(neighbors) while stack: p stack.pop() if not visited[p]: visited[p] True pn np.where(D[p] (local_eps[p] local_eps) / 2)[0] pn pn[pn ! p] if len(pn) min_samples: stack.extend(pn) if labels[p] -1: labels[p] cluster cluster 1 return labels几个实现细节说明一下local_eps用了中位数而不是均值受单个异常邻居影响更小。clip_quantile把极端大、极端小的局部半径截断防止个别稀疏离群点把半径拉到离谱程度。互认邻居条件D[i] (local_eps[i] local_eps) / 2是向量化一次性算完的速度比双循环快很多。扩张时对“已经被其他簇标记过的点”不做覆盖保持了边界点和簇归属的稳定性。3.3 用 k-距离曲线辅助找参考 eps如果你不想完全放弃多阈值路线可以用下面这个函数快速找到几个候选eps拿来检查结果或者做对比实验from scipy.signal import argrelextrema def detect_eps_list(X, k8, top_n3): nn NearestNeighbors(n_neighborsk).fit(X) dist, _ nn.kneighbors(X) k_dist np.sort(dist[:, -1]) diff np.diff(k_dist) # 一阶差分大说明曲线在这里斜率突变 max_idx argrelextrema(diff, np.greater)[0] if len(max_idx) top_n: max_idx np.argsort(diff)[-top_n:] eps_list k_dist[max_idx 1] return np.sort(eps_list)[::-1]这段代码找的是 k-距离曲线上的局部跳变点每一个跳变点都对应一个密度层级的边界。拿到的eps_list可以直接喂给标准 DBSCAN 做对比也可以反过来帮你在eps_scale上找到一个合理的基线值。用前面生成的模拟数据跑一下效果非常明显labels multi_dbscan(X, k8, min_samples5, eps_scale1.0) print(np.unique(labels, return_countsTrue))标准 DBSCAN 在这种数据上eps调小会丢掉稀疏簇调大又会合并密集簇。Multi-DBSCAN 几乎不用试参三坨密度差异很大的数据都能被正确分开。4. 在时序聚类里的落地特征窗口与行为分段4.1 时序数据如何变成聚类可用的样本时序聚类有两种常见思路一种是直接用原始序列做距离度量比如 DTW另一种是把时间序列切成窗口每个窗口提取统计特征变成普通向量再聚类。Multi-DBSCAN 对两种思路都适用。我推荐先用滑动窗口提取特征。下面这段代码把一个一维传感器序列变成一组特征向量def window_features(ts, window20, step5): X, index [], [] for start in range(0, len(ts) - window, step): seg ts[start:start window] X.append([ np.mean(seg), np.std(seg), (seg[-1] - seg[0]) / window, np.median(np.abs(seg - np.mean(seg))), ]) index.append(start window // 2) return np.array(X), np.array(index)然后对特征做标准化再丢进 Multi-DBSCANfrom sklearn.preprocessing import StandardScaler feat, idx window_features(ts, window20, step5) feat StandardScaler().fit_transform(feat) seg_labels multi_dbscan(feat, k6, min_samples3, eps_scale1.0)每个聚出来的簇对应一段行为模式平稳段、剧烈波动段、缓慢漂移段……因为不同行为段的波动幅度差异可能非常大比如静止时噪声几乎为零运动时波动剧烈这种场景天然就是多密度聚类该管的事。4.2 基于 DTW 距离的时序聚类变体如果时序形状才是关键而不是窗口统计量那就先算 DTW 距离矩阵再把距离矩阵传给multi_dbscanfrom tslearn.metrics import cdist_dtw D cdist_dtw(series_list) # 得到 n x n 距离矩阵 seg_labels multi_dbscan(D, k8, min_samples5, eps_scale1.0, precomputedTrue)这时候precomputedTrue会让函数跳过NearestNeighbors拟合直接从距离矩阵的行向量里取 k 近邻距离。要注意cdist_dtw返回的矩阵必须自己检查一下对角线是否为 0、是否对称否则局部半径估计会出问题。4.3 参数与效果观察时序分段场景下我一般把window取成业务上能解释的时长单位比如 20 个采样点约等于 2 秒那每个段就有明确的物理含义。step控制输出标签的密度step越小分段边界越精细但计算量也越大。调参时优先看两个现象分段太碎说明eps_scale偏小分段太长、把不同状态混在一起说明eps_scale偏大。时序数据和二维坐标点的区别在于特征维度里的“单位”可能并不一致所以标准化这步千万不要省否则均值和标准差的量纲差异会直接污染 k 近邻距离。5. 调参经验、坑点与工具对比5.1 参数组合速查表参数作用常见取值备注k计算局部半径时参考的近邻数6~12可以取min_samples - 1min_samples核心点判定阈值3~10太小会出大量碎片簇eps_scale全局松紧缩放系数0.8~1.51.0 时直接采用局部中位数clip_quantile局部半径截断范围(0.05, 0.95)防止异常点干扰5.2 常见问题实录与排查我把调试中遇到的实际问题整理成了一张速查表现象可能原因解决办法噪声点特别多eps_scale太小或k太小调大eps_scale增大k原本分开的簇被粘上互认半径过大边界点太多调小eps_scale降低k同一个簇内部被切碎局部半径抖动剧烈收紧clip_quantile分位数增大min_samples高维特征距离基本失去意义维度灾难先降维PCA/UMAP 都行再做聚类运行内存直接爆掉pairwise_distances存了 n x n 矩阵改用 KD-Tree 半径查询或用多阈值合并路线第 5 条值得多说两句。当n超过五万pairwise_distances的 n x n 矩阵就很难受了。这时候我一般会把数据分块或者改用“最大半径查询”的近似方案先算出所有局部半径的上限max_eps local_eps.max()用 KD-Tree 做半径搜索得到候选邻居再用互认半径条件过滤。这样虽然会多查一些候选点但内存占用从 O(n^2) 降到接近 O(n)。5.3 Multi-DBSCAN、OPTICS、HDBSCAN 怎么选不少朋友会问既然有现成的 OPTICS 和 HDBSCAN为什么还要自己写 Multi-DBSCAN。我的看法是OPTICS 本质上是用可达距离把不同密度层级的结构都保存下来再通过提取分簇来使用。它很强大但输出不是直接标签需要额外处理而且参数敏感度不高不代表不需要理解。HDBSCAN 在大多数情况下参数更省心对多密度数据非常友好但它内部基于最小生成树调参逻辑相对黑盒出问题以后很难手工干预。Multi-DBSCAN 的优势在于每一步都可解释局部半径是怎么来的邻居是怎么互认的谁被标成噪声都能逐个检查。对于要写进算法描述文档的项目这种透明性很值钱。如果项目看重自动化且没时间调参直接上 HDBSCAN 通常更省事如果要解决特定业务数据、还要跟团队解释每个簇的边界Multi-DBSCAN 这种“自己人写的透明算法”反而更好用。最后分享一个经验Multi-DBSCAN 不是万能的它对“密度连续渐变”的数据效果最好比如从市中心到郊区的扩散型分布。如果是多个完全独立的密度层级也就是簇与簇之间有明确空白区那标准的 DBSCAN 配合一次 k-距离曲线调参就够了不用自我折磨。我后来一直保留着这套自适应聚类代码主要就是看重它能少做一轮参数搜索并且每次调试时每个点都“心里有数”。