简介这是一份用Python实现的ST-DBScan时空密度聚类算法代码包面向需要处理空间聚类任务的开发者和算法学习者。ST-DBScan属于基于密度的聚类方法无需预先指定簇数量仅需半径和最小邻居数两个参数即可在不同密度的数据中找出形状各异的密集区域同时抑制噪声干扰可用于森林砍伐范围划定、医学影像中异常区域识别等场景。压缩包共8个文件整体约660KB包含3个Python源码、1个CSV示例数据、1张聚类效果图、1个README说明、1个License授权文件及Git忽略配置结构清晰便于直接阅读、运行和二次修改。目前已有707人学习下载源码附带示例数据与可视化输出能够帮助读者快速搭建实验环境并验证算法效果适合用于课程设计、算法实验或实际项目中的时空聚类模块开发是理解ST-DBScan从理论到落地的一份实用参考。1. 为什么需要 ST-DBSCAN从 DBSCAN 到时空聚类的边界有些场景簇的形状不是圆形密度也不均匀比如森林砍伐斑块和肿瘤区域。K-Means 这类质心聚类很难处理因为簇数未知且形状不规则。ST-DBSCAN 从密度连通分量出发不需要指定簇数只需给出邻域半径和最小邻居数就能把高密度区域连成簇并标记低密度点作为噪声。这份 Python 源码基于 [Sander et al., 1998] 的思路实现适合地理空间分析、遥感图像处理和医学辅助诊断场景。我拿到代码后先跑通了一个小样例再把核心循环拆开看发现细节都在邻域查询和边界判断上下面逐一展开。2. 算法原理与参数选型半径、最小邻域与密度可达2.1 核心定义ε 邻域、核心点、边界点与噪声ST-DBSCAN 的全部逻辑建立在一个前提上高密度区域被低密度区域分隔。算法把每个点分成三类。给定邻域半径 ε 和最小邻居数 MinPts若某点 ε 邻域内的样本数大于等于 MinPts它就是核心点若自身不是核心点但落在某个核心点的 ε 邻域内它就是边界点两者都不是的样本就是噪声点。簇定义为密度相连的核心点集合以及依附在这些核心点上的边界点。这里的「密度可达」是通过核心点之间直接或间接相邻来定义的。两个核心点距离小于 ε就称彼此直接密度可达一组核心点通过链式连接构成同一个簇的骨架。边界点不参与扩展但会被分配给第一个覆盖它的核心点所在的簇。这种设计让算法能够识别任意形状的簇同时对噪声有天然的隔离能力也正是它适合森林砍伐区域划定的原因——砍伐区域不会自动呈球形往往是不规则的多边形。需要特别留意的是ε 邻域内统计的样本数通常包含样本自身。项目中如果region_query返回的邻居列表包含查询点自己那么判断核心点的阈值 MinPts 就按包含自身计算否则要加 1。这种细微差别会让同一组参数产生完全不同的聚类结果。建议在阅读源码时先确认这一点。提示region_query返回的邻居数量是否包含查询点直接决定了 MinPts 的语义建议先看代码再定参数。2.1.1 时空数据下的密度定义当数据增加时间戳或属性维度后欧氏距离不再适合直接计算邻域。ST-DBSCAN 的做法是把距离扩展为混合度量。常见做法是空间部分仍用欧氏距离时间部分单独计算再加权合成。例如距离函数def st_distance(p, q, alpha1.0): spatial ((p[0]-q[0])**2 (p[1]-q[1])**2) ** 0.5 temporal abs(p[2]-q[2]) return spatial alpha * temporal在这段代码里p和q是由[x, y, t]构成的一维数组spatial表示欧氏距离temporal是时间差的绝对值。alpha权重越大时间差在邻域判断中占的比重越高适合对时间敏感的场景。如果你要处理多维度属性把它看作加权曼哈顿距离即可。实际地理位置坐标经纬度建议先投影到平面坐标系比如 UTM否则直接用经纬度的角度差计算距离没有物理意义。2.2 与普通 DBSCAN 的差异ST-DBSCAN 的处理维度普通 DBSCAN 只处理一个特征空间所有维度都视为同质距离。ST-DBSCAN 的「ST」指 space-time它把空间和时间/属性分开建模。时间维度和空间维度计量单位不同合在一起会产生量纲失配距离 1 公里的差距和 1 秒的差距不能直接相加。常见做法分两种一种是上面提到的加权线性组合另一种是用两个独立的阈值分别判断。维度普通 DBSCANST-DBSCAN距离函数统一欧氏距离空间距离 时间加权距离参数ε, MinPtsspatial_eps, temporal_eps/alpha, MinPts适用场景单一特征空间聚类周期事件识别、轨迹分析、区域演化噪声处理全局同阈值可分别控制空间和时间容忍度如果数据不止时间和空间例如加上污染浓度、道路等级等属性ST-DBSCAN 同样可以把这些属性作为附加距离项。只是每个属性都需要一个权重或归一化步骤否则高量纲特征会主导邻域判断。这个特性让算法不仅能找出空间聚集簇还能找出生理指标相近的病灶区域。2.3 参数选择策略ε 和 MinPts 怎么定ε 的选择通常借助 k-距离图把每个点到第 k 近邻的距离排序并绘制折线在曲线拐点处取 ε。MinPts 取 2 倍数据维度只是经验起点空间聚类一般取 4-10带时间维度时取 10-30 更稳妥。MinPts 太小会让少量散点形成核心点导致噪声被并入簇太大则会抹掉小簇。调参时应当先固定 MinPts再观察 k-距离图确定 ε两者交替微调。时间维度的 alpha 可以直接设置为空间距离尺度的比值。比如空间单位是公里时间单位是小时若认为 10 小时与 1 公里等价则 alpha 0.1。更稳健的办法是先对每一维做标准化再统一使用同一个 ε这样最省事。3. Python 实现源码拆解从距离计算到聚类输出3.1 项目结构与主流程解压后的项目包含py-st-dbscan-master目录下面是.gitignore、README.md、LICENSE、map.png和src源码目录。src里就是 Python 实现的全部代码。这类小项目一般不会采用复杂分层主流程通常是读入数据 → 初始化标签数组 → 逐个点判断是否访问过 → 邻域查询 → 扩展簇。map.png是示例数据在地图上的可视化结果用来确认聚类坐标系是否为经纬度或投影坐标。这种小规模聚类项目源码通常不超过几百行阅读顺序应该是main→cluster→distance。我解压后第一件事就是grep -n def列出所有函数。这样做能快速判断作者把距离度量和簇扩展分开了没有。我先看了src中的入口文件。入口负责解析命令行参数通常包括数据路径、空间半径、时间半径、最小邻居数。为了让代码可复现我建议用argparse代替硬编码参数这也是多数开源项目选择的写法。下面是一个可运行的数据读取片段import argparse import pandas as pd def load_data(path, cols): df pd.read_csv(path) return df[cols].values if __name__ __main__: parser argparse.ArgumentParser() parser.add_argument(-i, --input, requiredTrue) parser.add_argument(-e, --eps, typefloat, default100.0) parser.add_argument(-m, --min_pts, typeint, default5) args parser.parse_args()参数说明-i指定 CSV 路径-e是邻域半径-m是 MinPts。这里把args.eps直接传给聚类函数即可。注意pd.read_csv只接收文件路径如果数据带表头需要先确认列名或者用headerNone跳过。加载后的values是 NumPy 数组后续距离计算都基于此。3.2 距离计算与邻域查询实现ST-DBSCAN 性能瓶颈在邻域查询。最直接的实现是双重循环每个点与其他所有点计算距离。代码如下import numpy as np def naive_region_query(data, idx, eps): neighbors [] for i in range(len(data)): diff data[i] - data[idx] if np.sqrt(np.dot(diff, diff)) eps: neighbors.append(i) return neighbors这段代码用np.dot(diff, diff)计算向量内积避免np.linalg.norm的额外开销。对于几千个点双重循环还能接受到几万点时就非常慢。因此源码如果提供了scipy.spatial.cKDTree版本应优先使用。邻域查询返回的是索引列表不是点的坐标这样在标记标签时可以直接用索引操作。如果一个点已经被访问过则跳过邻域查询避免重复计算。3.3 簇扩展与噪声标记核心的扩展函数维护一个队列。初始种子是当前核心点的邻居从中取出一个点如果没访问过就标记为当前簇并检查它是否也是核心点若是则把它的邻居加入队列。这个队列可以用collections.deque实现出队复杂度 O(1)而list.pop(0)是 O(n)会拖慢大规模数据。from collections import deque def expand_cluster(data, labels, point_idx, cluster_id, eps, min_pts): seeds deque(naive_region_query(data, point_idx, eps)) if len(seeds) min_pts: labels[point_idx] -1 # 噪声 return False labels[point_idx] cluster_id while seeds: current seeds.popleft() if labels[current] -1: labels[current] cluster_id # 边界点归入 if labels[current] ! 0: continue labels[current] cluster_id neighbor naive_region_query(data, current, eps) if len(neighbor) min_pts: seeds.extend(neighbor) return True这里labels初始全 0-1表示噪声正整数表示簇编号。expand_cluster返回是否真的创建了新簇。注意边界点的处理如果当前点是之前某个簇标记的边界点-1 但尚未分配簇这里直接覆盖为当前簇这是 DBSCAN 的常见行为。如果一个点已经被访问过但不是噪声标签非 0continue跳过。这样避免重复扩展。这个函数里seeds.extend(neighbor)会引入大量重复点实际生产中需要再判断标签是否为 0或者用visited数组否则队列可能膨胀。可以从这里入手优化。3.4 结果输出与地图叠加聚类结束后标签数组需要写回文件并可视化。简单做法是直接把labels作为新列合并到 DataFrame。地图叠加时可以用matplotlib绘制散点图按簇编号着色。import matplotlib.pyplot as plt def plot_clusters(data, labels, output): plt.figure(figsize(8, 6)) for cid in np.unique(labels): mask labels cid if cid -1: plt.scatter(data[mask, 0], data[mask, 1], cgray, s5, labelnoise) else: plt.scatter(data[mask, 0], data[mask, 1], s5, labelfcluster {cid}) plt.legend() plt.savefig(output, dpi150)plot_clusters把同一簇的点画成同一种颜色噪声点用灰色区分。如果数据是经纬度底图可以用contextily或直接叠加到现成地图项目里map.png应该就是这么生成的。输出图片主要用来确认簇在空间上是否连续是否有异常孤立点。4. 实战在真实数据集上运行与参数调优4.1 数据准备CSV 格式要求与坐标预处理运行算法前数据至少要包含空间坐标列时间列可选。CSV 如果是下面这种格式可以直接用pandas.read_csv读入x,y,t 120.123,30.456,2024-01-01 08:00:00但日期字符串无法参与距离计算必须先转成数值。常见做法是把时间转换成 Unix 时间戳或者从起始时间开始按小时计数。更实际的做法是用项目自带的示例数据跑通流程再替换成自己的数据。地理坐标需要注意投影如果 x 和 y 都是经纬度用欧氏距离计算时 0.01 度的变化约等于 1 公里但不同纬度下这个比例会变化。建议先投影到 UTM 坐标系再计算距离否则高纬度地区会失真。由于 ST-DBSCAN 对量纲敏感我把数值预处理封装成一个函数from sklearn.preprocessing import StandardScaler def preprocess(df, time_col): df[t_num] (pd.to_datetime(df[time_col]) - pd.Timestamp(1970-01-01)).dt.total_seconds() scaler StandardScaler() df[[x_scaled, y_scaled, t_scaled]] scaler.fit_transform( df[[x, y, t_num]] ) return df, scalerStandardScaler会把空间和时间统一到均值为 0、方差为 1 的尺度之后所有维度共用同一个 eps 就合理了。这里有一个容易被忽略的点标准化后eps 的单位不再是原始米或秒而是标准差。所以调参时要从 0.1 到 1.0 之间去试不要用原来的 100 米之类的直觉值。4.2 命令行运行与参数解析项目源码如果支持命令行运行时应该类似python src/main.py -i data.csv -e 0.5 -m 10 --spatial 1 --temporal 1-e是标准化后的半径-m是最小邻居数--spatial和--temporal是权重标志。如果源码没有封装 CLI可以写一个 Python 包装脚本把参数传给 ST-DBSCAN 类。核心参数建议放在一个字典里集中管理方便批量实验。输出是新的 CSV 文件每行附带cluster列噪声点的值通常是-1实际业务里可以单独导出做异常分析。输出建议使用带 BOM 的 UTF-8 编码否则 Windows 环境下用 Excel 打开中文路径可能乱码。另外聚类结果的簇编号不是按面积或重要性排序的如果业务需要对簇排序可以在输出后按簇内点数降序排列把最大簇编号固定为 1。4.3 调参流程与效果评估先从 k-距离图确定 eps 候选值。把数据的第 k 近邻距离排序画出折线拐点对应的距离就是合适的 eps。代码如下from sklearn.neighbors import NearestNeighbors def knn_distance_plot(X, k): nn NearestNeighbors(n_neighborsk) nn.fit(X) dist, _ nn.kneighbors(X) dist np.sort(dist[:, -1]) plt.plot(dist) plt.xlabel(points sorted by distance) plt.ylabel(f{k}-th NN distance) plt.show()k可以取 MinPts-1因为要排除自身。曲线如果是陡峭上升的 L 形说明数据簇结构明显拐点前的平坦区长度对应簇内点的密度。有了候选 eps 后再用轮廓系数对比参数组合。对带噪声的密度聚类轮廓系数会把噪声视作差簇因此要同时报告噪声比例。评估指标计算方式适用场景轮廓系数(b-a)/max(a,b)簇紧密度与分离度噪声比例噪声点数 / 总点数数据质量与参数敏感度DBI簇内离散度与簇间距离比多簇形状相似时建议同时记录三个指标不要只看轮廓系数。我调参时通常固定 MinPts10把 eps 从 0.2 到 1.0 按 0.1 步进跑一遍选择噪声比例不超过 20% 且轮廓系数最高的组合。这样做的好处是能快速排除明显过小或过大的 eps。5. 进阶技巧索引加速、边界噪声与结果验证5.1 用 KD-Tree 把邻域查询从 O(n²) 降到 O(n log n)大规模时空数据下双重循环是主要瓶颈。我一般会用scipy.spatial.cKDTree替换naive_region_query只针对空间维度建树时间维度在过滤时再判断from scipy.spatial import cKDTree def indexed_region_query(data, spatial_idx, eps, temporal_epsNone): tree cKDTree(data[:, :2]) neighbors tree.query_ball_point(data[spatial_idx, :2], reps) if temporal_eps is not None: t0 data[spatial_idx, 2] neighbors [i for i in neighbors if abs(data[i, 2] - t0) temporal_eps] return neighborscKDTree.query_ball_point返回候选索引二次过滤时间维度后得到最终邻域。注意eps的单位与空间坐标一致temporal_eps独立控制时间差别。这个优化对几十万点也非常有效聚类本身仍是串行但邻域计算不再成为瓶颈。如果还会更慢考虑用joblib对点分块做邻域查询再合并结果但要注意避免跨块重复计算。5.2 边界点的归属与簇间重叠DBSCAN 的边界点可能同时落入多个核心点邻域最终归属取决于遍历顺序。ST-DBSCAN 加入时间维后这类情况更常见。解决思路是把核心点先找全用并查集把所有密度相连的核心点合并最后再遍历一次每个边界点选择距离最近的核心点簇。距离度量要和聚类时一致。这样结果不再依赖迭代顺序稳定性更好。5.3 结果验证的三种方式第一种是用合成数据验证正确性。生成几个已知形状的密集区域加上均匀噪声跑完聚类后与真实标签对比计算调整兰德指数ARI。第二种是检查实际场景的合理性比如森林砍伐区域聚类后每个簇是否对应一条完整的砍伐带是否存在跨越明显空隙的簇。如果簇面积异常大或有长条形穿透说明 eps 偏大。第三种是稳定性验证对同一数据重排行顺序或只取 80% 数据重复聚类比较两次结果的 ARI。若 ARI 低于 0.9说明参数处于临界位置边界点分配不稳定需要调小 eps 或增大 MinPts。上述验证可以在项目map.png基础上叠加生成结果用视觉判断后再用 ARI 量化。一次完整验证流程就封装成一个脚本python scripts/evaluate.py -i data.csv -e 0.5 -m 10 --truth labels.csv脚本内部计算 ARI、噪声比例和簇数输出到终端。这样每次调参都有量化反馈而不是只靠眼睛看散点图。本文还有配套的精品资源点击获取 时空聚类算法ST-DBSCAN解析：Python实现与调参实战
简介这是一份用Python实现的ST-DBScan时空密度聚类算法代码包面向需要处理空间聚类任务的开发者和算法学习者。ST-DBScan属于基于密度的聚类方法无需预先指定簇数量仅需半径和最小邻居数两个参数即可在不同密度的数据中找出形状各异的密集区域同时抑制噪声干扰可用于森林砍伐范围划定、医学影像中异常区域识别等场景。压缩包共8个文件整体约660KB包含3个Python源码、1个CSV示例数据、1张聚类效果图、1个README说明、1个License授权文件及Git忽略配置结构清晰便于直接阅读、运行和二次修改。目前已有707人学习下载源码附带示例数据与可视化输出能够帮助读者快速搭建实验环境并验证算法效果适合用于课程设计、算法实验或实际项目中的时空聚类模块开发是理解ST-DBScan从理论到落地的一份实用参考。1. 为什么需要 ST-DBSCAN从 DBSCAN 到时空聚类的边界有些场景簇的形状不是圆形密度也不均匀比如森林砍伐斑块和肿瘤区域。K-Means 这类质心聚类很难处理因为簇数未知且形状不规则。ST-DBSCAN 从密度连通分量出发不需要指定簇数只需给出邻域半径和最小邻居数就能把高密度区域连成簇并标记低密度点作为噪声。这份 Python 源码基于 [Sander et al., 1998] 的思路实现适合地理空间分析、遥感图像处理和医学辅助诊断场景。我拿到代码后先跑通了一个小样例再把核心循环拆开看发现细节都在邻域查询和边界判断上下面逐一展开。2. 算法原理与参数选型半径、最小邻域与密度可达2.1 核心定义ε 邻域、核心点、边界点与噪声ST-DBSCAN 的全部逻辑建立在一个前提上高密度区域被低密度区域分隔。算法把每个点分成三类。给定邻域半径 ε 和最小邻居数 MinPts若某点 ε 邻域内的样本数大于等于 MinPts它就是核心点若自身不是核心点但落在某个核心点的 ε 邻域内它就是边界点两者都不是的样本就是噪声点。簇定义为密度相连的核心点集合以及依附在这些核心点上的边界点。这里的「密度可达」是通过核心点之间直接或间接相邻来定义的。两个核心点距离小于 ε就称彼此直接密度可达一组核心点通过链式连接构成同一个簇的骨架。边界点不参与扩展但会被分配给第一个覆盖它的核心点所在的簇。这种设计让算法能够识别任意形状的簇同时对噪声有天然的隔离能力也正是它适合森林砍伐区域划定的原因——砍伐区域不会自动呈球形往往是不规则的多边形。需要特别留意的是ε 邻域内统计的样本数通常包含样本自身。项目中如果region_query返回的邻居列表包含查询点自己那么判断核心点的阈值 MinPts 就按包含自身计算否则要加 1。这种细微差别会让同一组参数产生完全不同的聚类结果。建议在阅读源码时先确认这一点。提示region_query返回的邻居数量是否包含查询点直接决定了 MinPts 的语义建议先看代码再定参数。2.1.1 时空数据下的密度定义当数据增加时间戳或属性维度后欧氏距离不再适合直接计算邻域。ST-DBSCAN 的做法是把距离扩展为混合度量。常见做法是空间部分仍用欧氏距离时间部分单独计算再加权合成。例如距离函数def st_distance(p, q, alpha1.0): spatial ((p[0]-q[0])**2 (p[1]-q[1])**2) ** 0.5 temporal abs(p[2]-q[2]) return spatial alpha * temporal在这段代码里p和q是由[x, y, t]构成的一维数组spatial表示欧氏距离temporal是时间差的绝对值。alpha权重越大时间差在邻域判断中占的比重越高适合对时间敏感的场景。如果你要处理多维度属性把它看作加权曼哈顿距离即可。实际地理位置坐标经纬度建议先投影到平面坐标系比如 UTM否则直接用经纬度的角度差计算距离没有物理意义。2.2 与普通 DBSCAN 的差异ST-DBSCAN 的处理维度普通 DBSCAN 只处理一个特征空间所有维度都视为同质距离。ST-DBSCAN 的「ST」指 space-time它把空间和时间/属性分开建模。时间维度和空间维度计量单位不同合在一起会产生量纲失配距离 1 公里的差距和 1 秒的差距不能直接相加。常见做法分两种一种是上面提到的加权线性组合另一种是用两个独立的阈值分别判断。维度普通 DBSCANST-DBSCAN距离函数统一欧氏距离空间距离 时间加权距离参数ε, MinPtsspatial_eps, temporal_eps/alpha, MinPts适用场景单一特征空间聚类周期事件识别、轨迹分析、区域演化噪声处理全局同阈值可分别控制空间和时间容忍度如果数据不止时间和空间例如加上污染浓度、道路等级等属性ST-DBSCAN 同样可以把这些属性作为附加距离项。只是每个属性都需要一个权重或归一化步骤否则高量纲特征会主导邻域判断。这个特性让算法不仅能找出空间聚集簇还能找出生理指标相近的病灶区域。2.3 参数选择策略ε 和 MinPts 怎么定ε 的选择通常借助 k-距离图把每个点到第 k 近邻的距离排序并绘制折线在曲线拐点处取 ε。MinPts 取 2 倍数据维度只是经验起点空间聚类一般取 4-10带时间维度时取 10-30 更稳妥。MinPts 太小会让少量散点形成核心点导致噪声被并入簇太大则会抹掉小簇。调参时应当先固定 MinPts再观察 k-距离图确定 ε两者交替微调。时间维度的 alpha 可以直接设置为空间距离尺度的比值。比如空间单位是公里时间单位是小时若认为 10 小时与 1 公里等价则 alpha 0.1。更稳健的办法是先对每一维做标准化再统一使用同一个 ε这样最省事。3. Python 实现源码拆解从距离计算到聚类输出3.1 项目结构与主流程解压后的项目包含py-st-dbscan-master目录下面是.gitignore、README.md、LICENSE、map.png和src源码目录。src里就是 Python 实现的全部代码。这类小项目一般不会采用复杂分层主流程通常是读入数据 → 初始化标签数组 → 逐个点判断是否访问过 → 邻域查询 → 扩展簇。map.png是示例数据在地图上的可视化结果用来确认聚类坐标系是否为经纬度或投影坐标。这种小规模聚类项目源码通常不超过几百行阅读顺序应该是main→cluster→distance。我解压后第一件事就是grep -n def列出所有函数。这样做能快速判断作者把距离度量和簇扩展分开了没有。我先看了src中的入口文件。入口负责解析命令行参数通常包括数据路径、空间半径、时间半径、最小邻居数。为了让代码可复现我建议用argparse代替硬编码参数这也是多数开源项目选择的写法。下面是一个可运行的数据读取片段import argparse import pandas as pd def load_data(path, cols): df pd.read_csv(path) return df[cols].values if __name__ __main__: parser argparse.ArgumentParser() parser.add_argument(-i, --input, requiredTrue) parser.add_argument(-e, --eps, typefloat, default100.0) parser.add_argument(-m, --min_pts, typeint, default5) args parser.parse_args()参数说明-i指定 CSV 路径-e是邻域半径-m是 MinPts。这里把args.eps直接传给聚类函数即可。注意pd.read_csv只接收文件路径如果数据带表头需要先确认列名或者用headerNone跳过。加载后的values是 NumPy 数组后续距离计算都基于此。3.2 距离计算与邻域查询实现ST-DBSCAN 性能瓶颈在邻域查询。最直接的实现是双重循环每个点与其他所有点计算距离。代码如下import numpy as np def naive_region_query(data, idx, eps): neighbors [] for i in range(len(data)): diff data[i] - data[idx] if np.sqrt(np.dot(diff, diff)) eps: neighbors.append(i) return neighbors这段代码用np.dot(diff, diff)计算向量内积避免np.linalg.norm的额外开销。对于几千个点双重循环还能接受到几万点时就非常慢。因此源码如果提供了scipy.spatial.cKDTree版本应优先使用。邻域查询返回的是索引列表不是点的坐标这样在标记标签时可以直接用索引操作。如果一个点已经被访问过则跳过邻域查询避免重复计算。3.3 簇扩展与噪声标记核心的扩展函数维护一个队列。初始种子是当前核心点的邻居从中取出一个点如果没访问过就标记为当前簇并检查它是否也是核心点若是则把它的邻居加入队列。这个队列可以用collections.deque实现出队复杂度 O(1)而list.pop(0)是 O(n)会拖慢大规模数据。from collections import deque def expand_cluster(data, labels, point_idx, cluster_id, eps, min_pts): seeds deque(naive_region_query(data, point_idx, eps)) if len(seeds) min_pts: labels[point_idx] -1 # 噪声 return False labels[point_idx] cluster_id while seeds: current seeds.popleft() if labels[current] -1: labels[current] cluster_id # 边界点归入 if labels[current] ! 0: continue labels[current] cluster_id neighbor naive_region_query(data, current, eps) if len(neighbor) min_pts: seeds.extend(neighbor) return True这里labels初始全 0-1表示噪声正整数表示簇编号。expand_cluster返回是否真的创建了新簇。注意边界点的处理如果当前点是之前某个簇标记的边界点-1 但尚未分配簇这里直接覆盖为当前簇这是 DBSCAN 的常见行为。如果一个点已经被访问过但不是噪声标签非 0continue跳过。这样避免重复扩展。这个函数里seeds.extend(neighbor)会引入大量重复点实际生产中需要再判断标签是否为 0或者用visited数组否则队列可能膨胀。可以从这里入手优化。3.4 结果输出与地图叠加聚类结束后标签数组需要写回文件并可视化。简单做法是直接把labels作为新列合并到 DataFrame。地图叠加时可以用matplotlib绘制散点图按簇编号着色。import matplotlib.pyplot as plt def plot_clusters(data, labels, output): plt.figure(figsize(8, 6)) for cid in np.unique(labels): mask labels cid if cid -1: plt.scatter(data[mask, 0], data[mask, 1], cgray, s5, labelnoise) else: plt.scatter(data[mask, 0], data[mask, 1], s5, labelfcluster {cid}) plt.legend() plt.savefig(output, dpi150)plot_clusters把同一簇的点画成同一种颜色噪声点用灰色区分。如果数据是经纬度底图可以用contextily或直接叠加到现成地图项目里map.png应该就是这么生成的。输出图片主要用来确认簇在空间上是否连续是否有异常孤立点。4. 实战在真实数据集上运行与参数调优4.1 数据准备CSV 格式要求与坐标预处理运行算法前数据至少要包含空间坐标列时间列可选。CSV 如果是下面这种格式可以直接用pandas.read_csv读入x,y,t 120.123,30.456,2024-01-01 08:00:00但日期字符串无法参与距离计算必须先转成数值。常见做法是把时间转换成 Unix 时间戳或者从起始时间开始按小时计数。更实际的做法是用项目自带的示例数据跑通流程再替换成自己的数据。地理坐标需要注意投影如果 x 和 y 都是经纬度用欧氏距离计算时 0.01 度的变化约等于 1 公里但不同纬度下这个比例会变化。建议先投影到 UTM 坐标系再计算距离否则高纬度地区会失真。由于 ST-DBSCAN 对量纲敏感我把数值预处理封装成一个函数from sklearn.preprocessing import StandardScaler def preprocess(df, time_col): df[t_num] (pd.to_datetime(df[time_col]) - pd.Timestamp(1970-01-01)).dt.total_seconds() scaler StandardScaler() df[[x_scaled, y_scaled, t_scaled]] scaler.fit_transform( df[[x, y, t_num]] ) return df, scalerStandardScaler会把空间和时间统一到均值为 0、方差为 1 的尺度之后所有维度共用同一个 eps 就合理了。这里有一个容易被忽略的点标准化后eps 的单位不再是原始米或秒而是标准差。所以调参时要从 0.1 到 1.0 之间去试不要用原来的 100 米之类的直觉值。4.2 命令行运行与参数解析项目源码如果支持命令行运行时应该类似python src/main.py -i data.csv -e 0.5 -m 10 --spatial 1 --temporal 1-e是标准化后的半径-m是最小邻居数--spatial和--temporal是权重标志。如果源码没有封装 CLI可以写一个 Python 包装脚本把参数传给 ST-DBSCAN 类。核心参数建议放在一个字典里集中管理方便批量实验。输出是新的 CSV 文件每行附带cluster列噪声点的值通常是-1实际业务里可以单独导出做异常分析。输出建议使用带 BOM 的 UTF-8 编码否则 Windows 环境下用 Excel 打开中文路径可能乱码。另外聚类结果的簇编号不是按面积或重要性排序的如果业务需要对簇排序可以在输出后按簇内点数降序排列把最大簇编号固定为 1。4.3 调参流程与效果评估先从 k-距离图确定 eps 候选值。把数据的第 k 近邻距离排序画出折线拐点对应的距离就是合适的 eps。代码如下from sklearn.neighbors import NearestNeighbors def knn_distance_plot(X, k): nn NearestNeighbors(n_neighborsk) nn.fit(X) dist, _ nn.kneighbors(X) dist np.sort(dist[:, -1]) plt.plot(dist) plt.xlabel(points sorted by distance) plt.ylabel(f{k}-th NN distance) plt.show()k可以取 MinPts-1因为要排除自身。曲线如果是陡峭上升的 L 形说明数据簇结构明显拐点前的平坦区长度对应簇内点的密度。有了候选 eps 后再用轮廓系数对比参数组合。对带噪声的密度聚类轮廓系数会把噪声视作差簇因此要同时报告噪声比例。评估指标计算方式适用场景轮廓系数(b-a)/max(a,b)簇紧密度与分离度噪声比例噪声点数 / 总点数数据质量与参数敏感度DBI簇内离散度与簇间距离比多簇形状相似时建议同时记录三个指标不要只看轮廓系数。我调参时通常固定 MinPts10把 eps 从 0.2 到 1.0 按 0.1 步进跑一遍选择噪声比例不超过 20% 且轮廓系数最高的组合。这样做的好处是能快速排除明显过小或过大的 eps。5. 进阶技巧索引加速、边界噪声与结果验证5.1 用 KD-Tree 把邻域查询从 O(n²) 降到 O(n log n)大规模时空数据下双重循环是主要瓶颈。我一般会用scipy.spatial.cKDTree替换naive_region_query只针对空间维度建树时间维度在过滤时再判断from scipy.spatial import cKDTree def indexed_region_query(data, spatial_idx, eps, temporal_epsNone): tree cKDTree(data[:, :2]) neighbors tree.query_ball_point(data[spatial_idx, :2], reps) if temporal_eps is not None: t0 data[spatial_idx, 2] neighbors [i for i in neighbors if abs(data[i, 2] - t0) temporal_eps] return neighborscKDTree.query_ball_point返回候选索引二次过滤时间维度后得到最终邻域。注意eps的单位与空间坐标一致temporal_eps独立控制时间差别。这个优化对几十万点也非常有效聚类本身仍是串行但邻域计算不再成为瓶颈。如果还会更慢考虑用joblib对点分块做邻域查询再合并结果但要注意避免跨块重复计算。5.2 边界点的归属与簇间重叠DBSCAN 的边界点可能同时落入多个核心点邻域最终归属取决于遍历顺序。ST-DBSCAN 加入时间维后这类情况更常见。解决思路是把核心点先找全用并查集把所有密度相连的核心点合并最后再遍历一次每个边界点选择距离最近的核心点簇。距离度量要和聚类时一致。这样结果不再依赖迭代顺序稳定性更好。5.3 结果验证的三种方式第一种是用合成数据验证正确性。生成几个已知形状的密集区域加上均匀噪声跑完聚类后与真实标签对比计算调整兰德指数ARI。第二种是检查实际场景的合理性比如森林砍伐区域聚类后每个簇是否对应一条完整的砍伐带是否存在跨越明显空隙的簇。如果簇面积异常大或有长条形穿透说明 eps 偏大。第三种是稳定性验证对同一数据重排行顺序或只取 80% 数据重复聚类比较两次结果的 ARI。若 ARI 低于 0.9说明参数处于临界位置边界点分配不稳定需要调小 eps 或增大 MinPts。上述验证可以在项目map.png基础上叠加生成结果用视觉判断后再用 ARI 量化。一次完整验证流程就封装成一个脚本python scripts/evaluate.py -i data.csv -e 0.5 -m 10 --truth labels.csv脚本内部计算 ARI、噪声比例和簇数输出到终端。这样每次调参都有量化反馈而不是只靠眼睛看散点图。本文还有配套的精品资源点击获取