简介这是一份基于Python语言实现的ST-DBScan算法源代码资源适合机器学习、数据挖掘方向的学习者及需要处理空间密度聚类的开发者使用。ST-DBScan作为经典空间聚类算法能够在含噪声的数据集中自动识别高密度连通分量无需预先指定簇数量仅需设置扫描半径与最小邻居数两个参数常用于森林砍伐区域划定、肿瘤区域识别等场景。压缩包共8个文件大小约660KB其中包含3个Python脚本用于核心算法与调用流程1个CSV示例数据便于直接运行测试另有图片文件辅助展示聚类效果以及LICENSE、README等说明文档结构简洁清晰。目前已有708人学习下载适合希望阅读源码、理解算法细节并在此基础上扩展应用的读者。通过这份资源可以快速掌握ST-DBScan的工程实现思路了解参数选择对聚类结果的影响并依据自带示例数据验证算法效果为后续二次开发或论文实验提供参考。1. 用 Python 实现 ST-DBScan 前先看它解决了什么你手上有一批带时间戳和经纬度的点可能是共享单车订单、出租车轨迹、维修工单也可能是犯罪记录。想找的是“在空间上聚在一起时间上也聚在一起”的簇比如某个路口连续一个月每晚都有单车淤积。直接用 DBSCAN 只比空间距离会把时间相隔半年的点全归成一类得到的是静态热点图而不是事件簇。ST-DBScan算法 正是在 DBSCAN 的邻域判定里增加时间维度空间上在 Eps1 范围内时间上在 Eps2 范围内才算邻居。这篇文章用纯 Python 写一版可运行的 ST-DBScan覆盖数据预处理、聚类代码、调参和踩坑。适合想搞懂算法内部逻辑、需要本地修改而不是调库交差的工程师。2. ST-DBScan 的核心逻辑与数据结构空间、时间两个半径怎么共同决定邻居2.1 从 DBSCAN 到 ST-DBScan多出来的不是坐标轴而是邻域判定条件DBSCAN 里“邻居”的定义只有一条距离小于 Eps。两个点只要空间够近就算邻居不管它们是不是同一天出现的。这个假设在静态空间聚类里没问题在时空数据里就会翻车。一个热门景点常年有人聚集静态 DBSCAN 会把它识别成一个全年无休的大簇但实际上白天和深夜的人群是两类事件需要分开治理。ST-DBSCAN 的改动很小把原来的邻域判定从“距离小于 Eps”改成“空间距离小于 Eps1 且 时间距离小于 Eps2”。注意这里是 AND不是 OR。两个点哪怕经纬度完全相同只要时间差超过 Eps2就不是邻居反之时间完全一致但相隔十万八千里同样不是邻居。这样聚出来的每个簇既能保证空间上紧凑又能保证时间上集中。为什么不直接把时间戳作为第三维塞进距离公式因为量纲问题。经纬度差是度时间差是秒两者直接相加没有任何物理意义Eps 也无法解释成空间半径或时间窗口。即便把时间归一化到 [0,1]不同应用的归一化尺度完全不一样调参没有可迁移性。ST-DBSCAN 拆成两个半径好处是每个参数都有直观含义可以分别用空间分布和时间分布去估计这也是它在时空聚类任务里远比“三维 DBSCAN”实用的原因。另一个容易忽略的点是“密度”的定义。在时空场景里一个点的时间邻居可能分布在完全不同的空间位置空间邻居也可能来自完全不同的时刻。ST-DBSCAN 把两者同时纳入邻域核心点的含义就变成了“在一个时空圆柱体里有足够多邻居的点”。这个圆柱体的底面半径是 Eps1高度是 2 倍 Eps2方向是时间轴。理解了这张图后面调 Eps1 和 Eps2 的思路就会清晰很多Eps1 决定圆柱底面在空间上圈多大Eps2 决定圆柱在时间轴上多高。2.2 输入数据长什么样时间戳、经纬度、对象ID的预处理约定我一般把输入整理成一个三列的 DataFrame每行是一个事件点或轨迹点longitude、latitude、timestamp。如果数据里带对象ID比如车辆ID或用户ID可以先按对象分组逐组做 ST-DBScan也可以全部放一起做取决于业务上是否需要“同一个对象的点才能相连”。轨迹停留点提取通常按对象分组事件热点聚类通常不分组。预处理第一件事是时间戳统一。CSV 里可能是“2024-05-01 08:30:00”也可能是 Unix 毫秒如果不统一成同一种数值类型后面算时间差全是灾难。我习惯全部转成 Unix 秒import pandas as pd df pd.read_csv(events.csv) df[timestamp] pd.to_datetime(df[timestamp]).astype(int64) // 10**9 print(df.dtypes)这里pd.to_datetime会把字符串解析成纳秒级时间戳astype(int64)取的是纳秒数整除10**9转成秒。统一成秒后Eps2 可以直接用 1800 表示半小时不用在每次循环里反复解析字符串。注意如果原始时间戳已经是毫秒级就不要整除10**9要整除10**6否则所有时间差都会缩小一千倍聚类结果会变得一团糟。这个坑在第五节详细说。经纬度处理上如果数据覆盖范围很小比如一个城市内的订单直接用经纬度做欧氏距离误差在低纬度地区不大但纬度越高越离谱稳妥做法是用 haversine 公式算球面距离或者把经纬度投影成平面坐标。我的建议是范围在几十公里内用 haversine 算出距离矩阵范围更大按数据所在区域选投影坐标系。核心点是 ST-DBSCAN 的 Eps1 需要和距离单位一致用米就都用米用度就都用度。时间戳这一列我不建议再保留成 Python 的 datetime 对象参与运算。datetime 之间虽然能做差但每次减法的成本比数值差高一个数量级数据量到几万时差距非常明显。统一成秒之后时间差计算就是一次整数减法后面所有邻域查询都能用向量化掩码一次算完。2.3 核心参数 Eps1、Eps2、MinPts、Delta 的分工与初始值估算ST-DBSCAN 有四个参数比 DBSCAN 多一个 Delta后两个有一点联动但整体是各管各的。Eps1 是空间邻域半径单位米或经纬度。它回答“多近算空间邻近”。Eps2 是时间邻域半径单位秒回答“多近算时间邻近”。MinPts 是密度阈值回答“邻域里至少几个点才算核心”。Delta 是 ST-DBSCAN 独有的回答“一个簇的时间跨度相对点数足够紧凑时是否允许吸收密度略低的边界点”。实际工程里我看到的落地实现中Delta 经常被弱化甚至省略。它最早在论文里用来处理簇边界的过渡噪声点判断依据是簇内“平均时间间隔”是否小于 Delta。不同代码实现对这个参数的定义差异较大有的做成后处理有的做进聚类主循环。为了可复现这篇采用后处理方式先跑完标准时空邻域聚类再对每个噪声点找最近簇如果该簇的时间跨度与点数之比小于 Delta就把它吸收进簇。这样不影响主循环逻辑出问题时也好排查。初始值估算有个固定套路。Eps1 用 K 距离曲线计算每个点到第 MinPts 近邻居的距离排序后画曲线取曲线拐点。Eps2 用时间差直方图把排序后的时间戳做一阶差分画直方图取主体区间的边界。MinPts 在没有先验知识时取 4和多数密度聚类实践一致。Delta 先不设等基础聚类跑出来观察簇的时间分布再定。要注意的是这四个参数不是独立的Eps1 调太大会让空间上松散的簇合并Eps2 调太大会把不同时段的事件揉在一起而 MinPts 越小的参数组合对噪声越敏感。3. 用 Python 手写 ST-DBScan 聚类核心代码与三个关键参数3.1 先搭骨架数据读取与距离计算函数用 vscode 配好 python 环境之后新建一个st_dbscan.py把下面这段骨架代码放进去。我习惯先写距离计算因为后面对邻域查询的调用都会复用这两个函数。import numpy as np import pandas as pd def haversine_distance(lon1, lat1, lon2, lat2): 计算两个经纬度坐标的球面距离单位米 R 6371000.0 lon1, lat1, lon2, lat2 map(np.radians, [lon1, lat1, lon2, lat2]) dlon lon2 - lon1 dlat lat2 - lat1 a np.sin(dlat / 2) ** 2 np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2) ** 2 return 2 * R * np.arcsin(np.sqrt(a)) def build_spatio_temporal_matrix(coords, timestamps): 构建空间距离矩阵和时间差矩阵只保留上三角再镜像对称节省一半计算量 n len(coords) dist_mat np.zeros((n, n)) time_mat np.zeros((n, n)) for i in range(n): for j in range(i 1, n): dist_mat[i, j] dist_mat[j, i] haversine_distance( coords[i, 0], coords[i, 1], coords[j, 0], coords[j, 1] ) time_mat[i, j] time_mat[j, i] abs(timestamps[i] - timestamps[j]) return dist_mat, time_mathaversine_distance把经纬度弧度化后算球面距离返回值是米。build_spatio_temporal_matrix一次性把空间距离和时间差都算成矩阵后面所有邻域查询都在这两个矩阵上用布尔掩码完成。注意这里用abs(timestamps[i] - timestamps[j])保证时间差非负排序方向不影响聚类结果。这个实现只适合数据量在几千量级的场景。n 个点会生成 n 乘 n 的两个矩阵一个一万点的数据集光距离矩阵就是 10000 乘 10000 的 float64占 800MB。对于更大数据量应在 3.2 节里用 KDTree 替代矩阵后面会专门讲。3.2 核心聚类主循环邻域查询与簇扩张ST-DBSCAN 的主循环和 DBSCAN 几乎一样差别只在邻域查询函数里多加了时间条件。下面这段是完整聚类函数def query_neighbors(i, dist_mat, time_mat, eps1, eps2): 返回点 i 的时空邻居距离小于 eps1 且时间差小于 eps2 的点索引 return np.where((dist_mat[i] eps1) (time_mat[i] eps2))[0] def st_dbscan(dist_mat, time_mat, eps1, eps2, min_pts): n dist_mat.shape[0] labels np.full(n, -1, dtypeint) # -1 表示噪声等 Delta 后处理再决定去留 cluster_id 0 for i in range(n): if labels[i] ! -1: continue neighbors query_neighbors(i, dist_mat, time_mat, eps1, eps2) if len(neighbors) min_pts: continue # 先保持噪声不急着分配簇 # 新建一个簇通过种子队列扩张 labels[i] cluster_id seeds set(neighbors) seeds.discard(i) while seeds: q seeds.pop() if labels[q] ! -1: continue labels[q] cluster_id q_neighbors query_neighbors(q, dist_mat, time_mat, eps1, eps2) if len(q_neighbors) min_pts: seeds.update(q_neighbors) cluster_id 1 return labels逻辑拆开看外层循环找没被标记的点如果它的时空邻域里至少有min_pts个点就把它当核心点开启一个新簇。seeds是该核心点的邻居集合每次取出一个点如果没被标记过就归入当前簇如果它自己也是核心点就把它的邻居也加进待处理集合。这样所有密度相连的点都会被同一个簇吸收。seeds.discard(i)是为了避免把当前点自己再处理一遍。参数说明eps1单位米eps2单位秒min_pts是整数。这个函数没有返回到底哪些点是核心点只返回每个点的簇标签。如果后续要做 Delta 吸收labels 里值为 -1 的就是候选噪声点。注意query_neighbors用的是而不是|这是 ST-DBSCAN 和普通 DBSCAN 的唯一本质区别漏掉时间条件就退化成了空间聚类。3.3 把 Delta 参数接进去噪声点与边界点的再判定Delta 的后处理逻辑很简单对每个噪声点找到它时空邻域里出现次数最多的簇如果那个簇本身时间足够紧凑就把这个噪声点吸收进去。紧凑度定义为簇内时间跨度除以簇内点数值越小代表单位时间的密度越高。def absorb_noise_by_delta(labels, timestamps, dist_mat, time_mat, eps1, eps2, delta): n len(labels) cluster_idx {} for i, lbl in enumerate(labels): if lbl ! -1: cluster_idx.setdefault(lbl, []).append(i) cluster_compactness {} for lbl, idxs in cluster_idx.items(): t_min min(timestamps[i] for i in idxs) t_max max(timestamps[i] for i in idxs) cluster_compactness[lbl] (t_max - t_min) / len(idxs) for i in range(n): if labels[i] ! -1: continue neighbors query_neighbors(i, dist_mat, time_mat, eps1, eps2) label_counts {} for nb in neighbors: if labels[nb] ! -1: label_counts[labels[nb]] label_counts.get(labels[nb], 0) 1 if not label_counts: continue best_label max(label_counts, keylabel_counts.get) if cluster_compactness[best_label] delta: labels[i] best_label return labelscluster_compactness是簇内时间跨度与点数的比值。Delta 设得越小越严格的簇才有资格吸收噪声点设得越大空间上靠近簇、时间上也比较集中的边界点都会被吸收但可能把一些孤立噪声也拉进来。Delta 的初始值我一般设成 Eps2 的两倍除以 MinPts。比如 Eps2 取 1800 秒、MinPts 取 4那 Delta 初始值就是 900。跑完看结果再微调。如果发现某个簇吸收了过多相隔很远的事件点说明 Delta 太大降一半重跑就行。注意 Delta 后处理只影响噪声点不影响已形成的簇结构所以可以反复调整而不必重跑主循环。4. 调参顺序与效果验证从 K 距离曲线到时空轮廓系数4.1 先用 K 距离图定 Eps1再用时间直方图定 Eps2调参顺序不能乱。我见过不少同事上来就凭感觉填 Eps1 和 Eps2结果聚类结果像撒胡椒面然后又去调 MinPts越调越乱。正确顺序是先定空间半径再定时间半径最后微调 MinPts 和 Delta。Eps1 用 K 距离曲线。数据量不大时可以直接用已经算好的dist_matimport matplotlib.pyplot as plt # 用 min_pts 作为第 K 个近邻的索引 k 4 k_distances np.sort(dist_mat, axis1)[:, k] k_distances_sorted np.sort(k_distances) plt.figure(figsize(8, 5)) plt.plot(k_distances_sorted) plt.xlabel(point index sorted by distance) plt.ylabel(distance to {}-th neighbor (m).format(k)) plt.grid(True) plt.show()曲线通常有一个明显拐点拐点左侧是簇内点右侧是离群点。Eps1 就取拐点对应的纵坐标。注意dist_mat对角线是 0所以要对每一行排序后取第 k 列k 对应的是 MinPts 的值。如果数据量超过一万dist_mat已经很大这一步可以抽样 2000 点算 K 距离趋势不会差太多。Eps2 不用 K 距离而是看时间间隔分布。把时间戳排序后做一阶差分画直方图time_steps np.sort(df[timestamp].values) time_diffs np.diff(time_steps) plt.figure(figsize(8, 5)) plt.hist(time_diffs, bins50, edgecolorblack) plt.xlabel(time gap between consecutive events (s)) plt.ylabel(count) plt.title(Time gap distribution) plt.show()如果事件在时间上均匀分布这个直方图会有一条长尾取主体区间右侧边界作为 Eps2。比如大部分时间差集中在 600 秒以内Eps2 取 600 或稍微放大到 900。如果直方图有几个峰说明事件天然分成几个时段Eps2 取第一个峰后的谷底比较合适。4.2 MinPts 与 Delta 的联动怎么避免“一个大簇吞掉所有点”MinPts 直接影响 Eps1 的取值因为 K 距离曲线的第 k 列就是按 MinPts 取的。MinPts 设得越小K 距离越小Eps1 也越小聚类结果会出现大量碎片簇设得越大K 距离越大Eps1 越大小簇被合并甚至出现一个簇把全城事件吞掉的情况。经验值MinPts 从 4 开始跑一次聚类统计簇数量和噪声比例。如果噪声超过 50%先把 MinPts 降到 3 或者把 Eps1 调大 20%如果簇数量超过预期数量的一倍说明 Eps1 太小或 MinPts 太大反向调整。Delta 的联动在于MinPts 调大后边界点更容易被判定为噪声Delta 后处理可以把其中一部分拉回来但不会改变簇的核心结构。所以正确姿势是先固定 MinPts调 Eps1 和 Eps2 到簇形态基本合理再开 Delta 做边缘吸收。千万不要一上来就开 Delta否则你根本分不清是主聚类参数的问题还是 Delta 的问题。4.3 验证聚类结果可视化着色与轮廓系数调参后要做两件事可视化着色和定量指标。可视化是最直观的验证plt.figure(figsize(10, 8)) scatter plt.scatter(df[longitude], df[latitude], clabels, cmaptab20, s8) plt.colorbar(scatter, labelcluster id) plt.xlabel(longitude) plt.ylabel(latitude) plt.title(ST-DBScan clustering result) plt.show()噪声点会被tab20色表里的某一色覆盖建议在绘图前把噪声点单独处理成灰色避免混淆。定性验证看两点簇是否在空间上紧凑同一个簇的点是否落在相近的时间段。如果同一个簇里既有凌晨又有中午说明 Eps2 偏大或者时间差矩阵没有正确参与邻域查询。定量验证用改进的轮廓系数。sklearn 的silhouette_score需要传入距离矩阵但 ST-DBSCAN 的邻域是时空双条件不能只用空间距离需要把空间和时间两个距离归一化后合并成一个综合距离矩阵from sklearn.metrics import silhouette_score max_spatial dist_mat.max() 1e-9 max_temporal time_mat.max() 1e-9 combined_dist dist_mat / max_spatial time_mat / max_temporal # 只对非噪声点计算排除孤立噪声对指标的干扰 mask labels 0 if len(set(labels[mask])) 1: score silhouette_score(combined_dist[mask][:, mask], labels[mask], metricprecomputed) print(silhouette score (spatio-temporal):, score) else: print(only one cluster, silhouette score is not meaningful)combined_dist把空间距离和时间差分别压到 [0,1] 再相加两个维度权重默认相等。如果你的业务里时间重要性更高可以把时间项乘 1.5 或 2也可以在归一化后用加权和。轮廓系数在 0.5 以上说明簇间分离度不错低于 0.2 说明参数不合适回到 4.1 的曲线重新调 Eps1 和 Eps2。5. 避坑指南ST-DBScan 落地的 5 个翻车现场5.1 现象不同时间的事件被聚成一类聚类结果和普通 DBSCAN 没区别原因邻域查询里漏掉了时间条件或者把写成了|。我用 Python 做python 代码调试时最常遇到的就是这个query_neighbors函数里只写了一个 eps1时间矩阵压根没参与。还有人图省事直接用 sklearn 的 DBSCAN给它喂了经纬度加时间戳三列结果把时间戳当第三维坐标聚出来全是空间主导。解决检查query_neighbors的返回条件确认同时包含dist_mat[i] eps1和time_mat[i] eps2。再用一个简单测试验证取两个经纬度相同的点一个时间为 0一个时间为 3600设 Eps2 为 1800这两个点必须被分到不同簇或其中一个被标为噪声。把这个测试写进单元测试里以后改代码不怕回归。5.2 现象时间戳格式统一后聚类结果完全变了之前跑得好好的结果突然不对劲原因原始时间戳有的是秒有的是毫秒统一成秒时整除的系数不对。CSV 里的时间是 “2024-05-01 08:30:00” 这种字符串还好用pd.to_datetime转出来是纳秒级。但如果你从数据库里导出的字段已经是 Unix 毫秒直接astype(int64) // 10**9会把所有时间差压到接近 0Eps2 形同虚设。解决先看时间戳的数值范围。计算max(timestamp) - min(timestamp)如果结果是十位数的量级基本是秒如果是十三位数是毫秒。统一转换逻辑秒级原样用毫秒级整除10**6如果是从纳秒字符串转的整除10**9。写一个函数处理不要在每个脚本里各写一遍转换逻辑。5.3 现象数据量到一两万点聚类跑了十几分钟还没结束内存先崩了原因build_spatio_temporal_matrix生成的 n 乘 n 矩阵空间和时间两个 float64 矩阵加起来 1.6GBn10000。再加上中间变量内存很容易爆。主循环里的set和np.where每次都要扫一整行计算量也是 O(n²) 级别。解决数据量大时不要构建全矩阵改用空间索引加时间过滤。常见做法是用 scipy 的cKDTree做空间近邻查询先按 Eps1 找到空间候选点再在候选点里过滤时间差小于 Eps2 的点。下面这段可以直接替换query_neighborsfrom scipy.spatial import cKDTree # 把经纬度投影成平面坐标单位米这里用简单的等距投影 def project_xy(lon, lat): R 6371000.0 x R * np.radians(lon) * np.cos(np.radians(lat.mean())) y R * np.radians(lat) return np.column_stack([x, y]) coords_xy project_xy(df[longitude].values, df[latitude].values) tree cKDTree(coords_xy) def query_neighbors_fast(i, timestamps, eps1, eps2): spatial_idx tree.query_ball_point(coords_xy[i], eps1) time_diff np.abs(timestamps[spatial_idx] - timestamps[i]) return np.array(spatial_idx)[time_diff eps2]注意project_xy里用了lat.mean()做等距投影的基准纬度这只在数据范围较小、纬度跨度不大时才准。严格的做法是按数据中心所在 UTM 分带投影。主循环改成邻域查询返回数组后其他逻辑不用动一万点可以在几分钟内跑完。5.4 现象MinPts 设成 2结果所有点都被吞进同一个簇找不到噪声原因MinPts2 时任意两个相邻点都能互相成为核心点簇会沿着时间上的链条不断扩展早上 8 点的点连接 9 点9 点连接 10 点最后一整天的点连成一条长链。这在时间聚类里尤其严重因为时间差是可以无限累加的两个相隔很远的事件可以通过中间点“传导”成密度相连。解决MinPts 至少 3推荐 4 或 5。如果你确实希望簇比较细碎可以保持 MinPts3但要把 Eps2 调小切断时间长链。更系统的做法是固定其他参数扫描 MinPts 取 2、3、4、5、6画簇数量和噪声比例曲线找到曲线变平的位置。5.5 现象低纬度数据聚类结果正常高纬度城市的数据偏移严重簇形状被拉长原因直接拿经纬度算欧氏距离在纬度 60 度的地方1 度经度的实际地面距离只有赤道的一半但代码里 1 度经度还是当 1 度算。聚类结果会在经度方向上被“放大”簇看起来横向拉长。解决所有距离计算一律用haversine_distance或者先做投影再进 KDTree。不要用np.linalg.norm(x1 - x2)直接算经纬度坐标差。如果项目里已经有经纬度欧氏距离的老代码建议用一个小数据集做对比测试分别用两种距离跑聚类统计簇数量和成员差异大概率会发现高纬度簇形状差异明显。6. 进阶把 ST-DBScan 用成轨迹停留点提取器ST-DBScan 最常见的落地场景之一是轨迹停留点提取。车辆或行人的轨迹点每隔几秒采一个停在一个地点时会有一堆空间接近、时间集中的点。用 ST-DBScan 聚类后每个簇就是一次停留簇中心的经纬度是停留位置簇内时间跨度是停留时长。具体做法先按对象ID分组对每个对象的轨迹点单独跑 ST-DBScan。参数上Eps1 取停留半径比如 50 米Eps2 取最小停留时长比如 5 分钟就是 300 秒MinPts 取 3保证至少 3 个连续点确认一次停留。聚类完成后把每个簇的点取经纬度中位数作为停留点坐标输出一行记录results [] for obj_id, group in df.groupby(object_id): coords group[[longitude, latitude]].values times group[timestamp].values dist_mat, time_mat build_spatio_temporal_matrix(coords, times) labels st_dbscan(dist_mat, time_mat, eps150, eps2300, min_pts3) for lbl in set(labels): if lbl -1: continue idxs np.where(labels lbl)[0] results.append({ object_id: obj_id, lon: np.median(coords[idxs, 0]), lat: np.median(coords[idxs, 1]), start_time: times[idxs].min(), end_time: times[idxs].max(), duration: times[idxs].max() - times[idxs].min(), point_count: len(idxs), }) stay_df pd.DataFrame(results)参数验证上我习惯做一个二维扫描固定 MinPts让 Eps1 从 20 米到 200 米、Eps2 从 60 秒到 600 秒变化画出停留点数量变化曲线。曲线会有一个明显平台区平台区内的参数组合对结果影响不大选中间值最稳。如果平台区很窄说明数据本身停留特征模糊这时候不要硬调参应该回去看原始轨迹是否有周期性抖动。我早期在这个项目上踩过最深的一次是 Eps2 设太大把三次短暂停车连成了一次长停留后来改成先做轨迹点时间下采样再聚类问题才解决。参数这东西跑出来的结果能解释业务现象才算数。希望帮到你。本文还有配套的精品资源点击获取