ST-DBSCAN时空聚类实战：Python实现与参数调优

简介：ST-DBSCAN算法Python实现代码包，面向具备一定Python基础的数据分析与机器学习开发者，用于处理带噪声的空间点数据聚类任务。该算法最大特点是不需预先指定簇的数量，而是依据半径与最小邻居数两个参数，自动识别高密度连通分量，从而得出不同形状、密度相近的簇；因此尤其适合森林砍伐范围划定、医学影像中肿瘤影响区域识别等空间数据场景。资源共8个文件，压缩包仅660KB，内含3个Python源码文件（核心算法与调用示例）、1个CSV示例数据集、1张PNG示意图、1份Markdown说明文档与LICENSE许可证，并配备.gitignore工程文件，目录结构简洁，便于对照文档快速上手。目前已有707人学习下载，这份代码可帮助读者深入理解ST-DBSCAN的分簇逻辑与参数影响，也能作为工程化改造的基础模板；无论用于课程设计，还是在此基础上扩展时空大数据聚类应用，都具有较好的参考价值。

1. ST-DBSCAN 要解决什么问题：时空聚类的边界与演化

处理一周的外卖骑手轨迹数据时，会出现一个反直觉现象：把全部点直接丢给 DBSCAN，结果市中心那个簇直径只有 800 米，时间跨度却有 14 个小时——凌晨 2 点的夜宵单和下午 6 点的晚高峰单被归进同一个簇。问题不在 DBSCAN，而在于输入里时间被当成了普通属性而不是约束条件。ST-DBSCAN（Spatio-Temporal DBSCAN）把邻域定义拆成两个独立半径：空间阈值 Eps1 和时间阈值 Eps2，一个点要被归入某簇，必须同时满足“你在它附近”且“你们发生在同一段语义时间里”。这是共享单车潮汐识别、轨迹停留点检测、城市事件聚集分析里最常用的密度聚类方案，比给 DBSCAN 多加一列 timestamp 要严谨得多。下面按“原理 → Python 实现 → 参数标定 → 簇演化分析”的顺序推进，代码可以直接下载改造成你自己的时空数据分析工具。

2. ST-DBSCAN 邻域定义与密度扩展：两个半径如何处理时间维度

2.1 从 DBSCAN 到 ST-DBSCAN 的邻域变化

原始 DBSCAN 的邻域是纯空间的：给定半径 Eps，点 p 的邻域 N(p) = {q | dist(p,q) <= Eps}。ST-DBSCAN 把 N(p) 切成两个独立条件：

N(p) = {q | dist(p,q) <= Eps1 AND |t_p - t_q| <= Eps2}

Eps1 仍然叫空间邻域半径，单位是米或经纬度度数；Eps2 叫时间邻域半径，单位是秒、分钟或小时。这里的关键是 AND：两个条件都成立才构成时空邻居。如果只满足空间近（比如都在同一个路口）但时间相隔两小时，那么这两个点彼此不可达。

这个 AND 有直接的数学后果：ST-DBSCAN 的邻域大小受时间维度约束，聚类结果是对空间密度的时间切分。同样的空间热点，如果一天里出现了三个不同的聚集高峰，ST-DBSCAN 会得到三个不同的簇，而不是一个时间上坍缩的大簇。这是它与带时间属性的 DBSCAN 的本质区别：DBSCAN 不会因为表里多放一列 timestamp 就把时间变成约束，它默认 timestamp 参与了欧氏距离计算，结果受单位影响极大——经纬度差 0.001 和时间差 600 秒根本没有可比性。

2.2 密度可达的传递性与时间窗口的截断效应

DBSCAN 里密度可达的定义是存在一条链 p=p_0, p_1, ..., p_n=q，每一步都是从核心点出发的邻域内的点。ST-DBSCAN 沿用这个定义，但因为每一步的邻域判断都加入了 |t_{p_i} - t_{p_{i+1}}| <= Eps2，链式扩展的时间跨度会被逐步累积。

举个例子：设 Eps2 = 10 分钟，采样点 p0 在 09:00，p1 在 09:08，p2 在 09:16。p0 和 p1 满足时间邻域，p1 和 p2 满足时间邻域，但 p0 和 p2 相差 16 分钟，不满足直接邻域。如果 p0 是核心点，从 p0 出发扩展时，p2 首先要通过 p0 的邻域检查，时间差超限，扩展失败。但如果 p1 先被扩展，而 p1 在扩展前尚未标记为其它簇的成员，那么 p2 会被划入 p1 所在簇，这一刻簇的时间跨度已经超过 Eps2。也就是说：时间窗口是逐跳约束，不是簇全局约束，簇内首尾两点的时间差可能超过 Eps2。

这个特性常被忽略，但它是理解结果里“簇的时间宽度”的钥匙。做轨迹停留点检测时，如果希望簇的总时间跨度严格不超过某个阈值，需要另加一个后置过滤条件（见第 5 章），ST-DBSCAN 本身不做全局时间闭合。

2.3 核心概念与数据预处理要求

2.3.1 核心点、边界点与噪声点的判定

沿用 DBSCAN 的三类点定义，判定条件是时空邻域内的点数 >= MinPts。边界点的判定要特别注意：边界点不要求自身是核心点，但它必须位于某个核心点的时空邻域内。一个常见的实现误区是，边界点如果同时落在两个簇的邻域里，会被先访问的那个簇吸收，ST-DBSCAN 原版不解决“边界归属歧义”。对大多数事件聚集场景这个歧义影响不大，但如果你要做严格的簇关系图谱，需要在聚类后按空间重叠对簇做二次合并。

2.3.2 预处理：时间换算与数据排序

时间维度在实现前要先统一成数值。常见做法是把时间戳统一成 Unix 秒，或者把一天内的时间换算成当天秒数（0~86400）。前者适合连续多日轨迹，后者适合只关心日内模式的场景。还需要按时间排序或建立时间索引，这能让邻域查询提前剪枝：查询一个点的空间邻域后，用 Eps2 对候选时间做一次区间过滤。百万级数据量下，按时间排序后可以用二分查找定位时间窗口边界，再与空间候选取交集，能省掉大量无效的距离计算。

排序还有一个隐藏收益：簇的编号会带上“时间先后”的属性，便于后续分析中按时间回放簇的出生与消亡。

2.4 ST-DBSCAN 与 DBSCAN、OPTICS、层次聚类的选型对照

这个对照表的结论很直接：如果你的目标不是“找形状”，而是“找在同一时间段内反复出现的空间聚集”，ST-DBSCAN 是密度聚类里对现有代码改动最小的选择。层次聚类虽然也能表达时空语义，但距离矩阵的 O(n²) 成本在小样本之外基本不可用。

3. Python 实现 ST-DBSCAN：从暴力邻域到 cKDTree 加速

3.1 数据表结构与时空字段解析

实现前先约定输入格式。ST-DBSCAN 最小输入只需要三列：longitude（经度）、latitude（纬度）、timestamp（时间戳）。我用 pandas 读入后直接转成 numpy 数组，避免逐行遍历 DataFrame 的开销。

import numpy as np import pandas as pd df = pd.read_csv("geolife_sample.csv", parse_dates=["timestamp"]) # timestamp 统一转成 Unix 秒，后续时间差计算才是纯数值运算 df["ts"] = df["timestamp"].astype("int64") // 10**9 coords = df[["longitude", "latitude"]].values.astype(np.float64) times = df["ts"].values.astype(np.float64)

参数说明：astype 转换时间的目的是让 |t_p - t_q| <= Eps2 的判断退化成一次浮点减法。很多初版实现用 datetime 对象直接做差会慢一个数量级，因为 datetime 减法每次都会构造 timedelta 对象。实际项目中如果 timestamp 已经是 int 类型，可以直接跳过 parse_dates 这一步。

3.2 空间距离：Haversine 还是平面近似

500 米以内的聚类（比如共享单车潮汐点），直接用经纬度做平面欧氏会有偏差：纬度 1 度约 111 公里，经度 1 度在赤道约 111 公里，但在北纬 40 度约 85 公里。如果不做投影，Eps1 的物理含义会随纬度漂移。常见做法是：小范围数据（城市级）先用等距圆柱投影把经纬度转成米制坐标，再做欧氏距离；或者直接用 Haversine 计算球面距离。这里给出投影转换的最小实现：

def lonlat_to_meters(lon, lat, ref_lat=None): # 等距圆柱投影，ref_lat 取数据集中位数纬度，减小投影变形 if ref_lat is None: ref_lat = np.median(lat) k = 111320.0 * np.cos(np.deg2rad(ref_lat)) x = lon * k y = lat * 111320.0 return np.column_stack([x, y])

参数说明：111320 是每纬度对应的米数近似；乘 cos(ref_lat) 是对经度方向的尺度压缩。投影后 Eps1 的单位从“度”变成“米”，聚类半径的语义变得可解释。城市级数据（半径 1~2 公里内）这个近似的误差小于 1%，不需要引入 pyproj 的重型依赖。

提示：投影的参考纬度取数据集中位数即可。如果你的数据横跨多个城市，按城市分块做聚类，不要用同一个 ref_lat 强行投影全量数据。

3.3 数据结构与索引设计

ST-DBSCAN 的邻域查询实际是两个查询的交集：

1. 空间近邻：在半径 Eps1 内找到所有候选点；
2. 时间过滤：在候选中筛出 |dt| <= Eps2 的点。

为了避免在第二步里遍历全部候选，我会对时间数组建排序索引，用 np.searchsorted 做二分查找；空间近邻用 scipy.spatial.cKDTree。cKDTree 在二维数据的查询复杂度接近 O(log n)，替代了暴力 O(n) 的遍历。注意 cKDTree 只能处理欧氏距离，所以经纬度必须先投影成米制坐标，否则查询条件里 Eps1 要换算成度数，而度数随纬度变化，索引查询会出错。

3.4 完整可运行的 ST-DBSCAN 代码

下面的代码实现一个时空密度聚类器，保存为 stdbscan.py 后可直接 import 使用。为了可读性，我把查询函数和聚类主流程分开封装。

import numpy as np from scipy.spatial import cKDTree class STDBSCAN: def __init__(self, eps1=150.0, eps2=600.0, min_pts=5, projection=None): """ eps1: 空间半径，单位米（投影后） eps2: 时间半径，单位秒 min_pts: 时空邻域内最少点数 projection: 可选的 (lon, lat) -> (x, y) 函数 """ self.eps1 = eps1 self.eps2 = eps2 self.min_pts = min_pts self.projection = projection self.labels_ = None def _neighborhood(self, idx, tree, times): # 1. cKDTree query_ball_point 得到空间候选 candidates = np.asarray( tree.query_ball_point(self.pts[idx], self.eps1), dtype=np.int64 ) # 2. 时间过滤：向量化做差值比较，避免 Python 循环逐个判断 t = times[idx] dt = np.abs(times[candidates] - t) return candidates[dt <= self.eps2] def fit(self, lon, lat, timestamps): n = len(lon) if self.projection is not None: self.pts = self.projection(lon, lat) else: self.pts = np.column_stack([lon, lat]) times = np.asarray(timestamps, dtype=np.float64) tree = cKDTree(self.pts) visited = np.zeros(n, dtype=bool) # label 为 -1 表示噪声，0..k 表示簇编号 labels = np.full(n, -1, dtype=np.int64) cluster_id = 0 for i in range(n): if visited[i]: continue visited[i] = True nbrs = self._neighborhood(i, tree, times) if len(nbrs) < self.min_pts: continue # 暂标记为噪声候选，后面可能被核心点吸收 labels[i] = cluster_id seeds = set(nbrs) # 用集合做去重，避免同一个边界点被重复压栈 while seeds: q = seeds.pop() if not visited[q]: visited[q] = True q_nbrs = self._neighborhood(q, tree, times) if len(q_nbrs) >= self.min_pts: # q 是核心点，其时空邻域并入扩张集合 seeds.update(q_nbrs) if labels[q] == -1: # 边界点归入当前簇；核心点早前已被标记 labels[q] = cluster_id cluster_id += 1 self.labels_ = labels return labels def fit_predict(self, lon, lat, timestamps): return self.fit(lon, lat, timestamps)

逻辑说明：流程与原始 DBSCAN 的扩张式一致，区别只在 _neighborhood 里多了时间过滤这一行。注意核心点的标记时机：我在将邻域种子入栈前就把 i 标记为 cluster_id，在弹出种子点后，若该点是核心点则保持已标记状态，只有边界点才会走到if labels[q] == -1分支。这里也处理了一种边界情形：某点在主循环里曾因邻域不足被当成噪声候选，后来落在了一个核心点的时空邻域内，它会被重新吸收为边界点，这正是 DBSCAN 的语义。

参数说明：eps2=600 表示时间窗口为 10 分钟，适合订单类数据；min_pts=5 表示至少需要 5 条记录在同样的时空窗口中才算聚集。projection 参数直接传第 3.2 节的 lonlat_to_meters 函数即可。如果数据量只有几千条且不想装 scipy，可以删除 cKDTree 相关代码，直接对全量点算距离矩阵。

3.5 大数据的两个优化点与失败时的排查

3.5.1 减少重复邻域查询

cKDTree 构建本身是 O(n log n)，百万点规模大约几秒到十几秒。真正的瓶颈在 _neighborhood 会被调用多次。第一个优化：只在主循环里对未访问点执行邻域查询，跳过已聚类的点；第二个优化：如果某个簇扩展特别慢，把 seeds 换成双向链表结构的自定义集合，避免 set 的 pop 顺序随机导致 CPU cache miss。对绝大多数分析任务，前一个优化已经足够。

3.5.2 参数错误导致的典型失败特征

如果所有点都聚成 1 个簇，通常是 eps1 或 eps2 设置过大，或者 min_pts 设成了 1。如果几乎全是噪声（标签全为 -1），先检查投影函数是否用了全局参考纬度，再把 eps1 临时放大 10 倍观察簇数量是否下降，用来判断是不是空间半径太小。时间维度的失败特征是：空间上明显成团、但聚类结果零散——这时优先调大 eps2，而不是 eps1。

4. 参数标定：Eps1、Eps2 与 MinPts 的经验区间和自动估计

4.1 三个参数的业务含义与初始值

参数定标不需要一开始就追求最优，先按表里的业务含义给一组可解释的初值，跑通全流程后再做网格搜索。要注意这三个参数不是独立的：Eps2 决定了一个簇在时间轴上“活多久”，MinPts 决定了在这个存活窗口内要凑够多少人。

4.2 k-距离曲线：Eps1 的客观估计方法

k-距离的基本逻辑：对每个点，计算它到第 k 近邻的距离，按距离升序排列并画曲线，曲率最大处对应的距离就是候选 Eps1。这个 k 一般取 MinPts - 1 或 MinPts。注意这个曲线只考虑空间距离，不做时间过滤，因为它估计的是“空间上集聚的尺度”，而不是“同时出现的尺度”。

from scipy.spatial import cKDTree def find_eps1(pts_xy, k=5): # pts_xy 是投影后的米制坐标，不是经纬度 tree = cKDTree(pts_xy) # k+1 是因为 query 会把自身当第一近邻 dist, _ = tree.query(pts_xy, k=k + 1) kth = np.sort(dist[:, -1]) return kth

输出后，按距离为 x 轴、kth 为 y 轴画折线，找斜率突变点。斜率突变点的含义是：距离小于该值时点数增长极快，说明这些点都嵌在密集区；距离大于该值后曲线变平，说明进入了渐远邻域。这是密度聚类里最经典的启发式，比拍脑袋定半径可靠得多。

4.3 Eps2 的两种标定路径：采样间隔倍数与事件语义

路径一是按数据采样率：如果轨迹是每秒采样一次，Eps2 取 30~60 秒；如果是订单数据，最小语义间隔是 1 分钟，Eps2 取 5~15 分钟。路径二是按目标事件：检测“红绿灯前滞留”用 60~120 秒，检测“早高峰办公区聚集”用 30 分钟。Eps2 过小的症状是：同一空间簇被切成时间碎片；过大的症状是：不同通勤波次被合并成一个簇。

4.4 联合调参：网格搜索加一个可解释的评估指标

ST-DBSCAN 没有类别标签时，常用轮廓系数做内部评估，但纯空间轮廓系数无法反映时间正确性。常见做法是构造一个时空版本的轮廓系数：a_i 是点 i 到同簇内其他点的时空距离均值，b_i 是到最近异簇的时空距离均值，时空距离定义为两个维度的加权欧氏距离。权重靠人工给定，不参与搜索。

from sklearn.metrics import silhouette_samples def spatiotemporal_silhouette(pts_xy, times, labels, time_weight=0.5): # time_weight 表示时间秒数折算成米数的比例，例如 0.5 = 0.5 米/秒 scaled = np.column_stack([pts_xy[:, 0], pts_xy[:, 1], times * time_weight]) return silhouette_samples(scaled, labels).mean()

参数说明：time_weight 取多少决定了时间维度在评估中的权重。0.5 表示每秒记 0.5 米，适合低速聚集场景；骑行场景可以取 5.0，即每秒最多移动 5 米。这个缩放只用于评估，不用于聚类本身。网格搜索时以这个指标为目标，每次跑 ST-DBSCAN 后计算分数，取最大值的参数组合。几千条数据下搜索成本可接受；百万条数据时建议先抽样一版做搜索。

4.5 MinPts 与 Eps2 的数量联动

MinPts 的一个可靠公式：MinPts >= 2 ×（Eps2 窗口内平均单点出现次数）。如果采样周期是 5 秒、Eps2 = 60 秒，一个点最多出现 12 次，MinPts 取 5~6 就能识别人群聚集；如果点只是一天一次的到访打卡，MinPts 取 2~3 才有意义，此时 ST-DBSCAN 退化成“共同出现”检测器。要注意 MinPts=1 会让每个点都自成一簇，因为每个点的邻域至少包含它自己，这是实现里要严防的。

5. 从静态簇到生命周期：用时间桶输出簇的演化轨迹

5.1 簇生命周期的时间桶设计

聚类完成后，一个簇在时间轴上可能跨越几分钟到几天。做事件检测时，最关心的是簇的新增、消失、合并、分裂。标准做法是把时间切成等宽时间桶（桶宽建议取 Eps2 的 1~2 倍），统计每个簇在每个桶内的成员集合，再用 Jaccard 相似度判断相邻桶的同一性。桶宽取 Eps2 的整数倍，是为了让一个“语义时间窗口”刚好落在 1~2 个桶里，既不丢细节，也不把单点抖动放大成事件。

def cluster_timeline(labels, times, bucket_size): t0 = times.min() num_buckets = int((times.max() - t0) / bucket_size) + 1 timeline = [[] for _ in range(num_buckets)] for label, t, i in zip(labels, times, range(len(labels))): if label == -1: continue # 噪声点不参与演化分析 b = int((t - t0) / bucket_size) timeline[b].append((label, i)) return timeline

这个函数返回每个时间桶内的簇成员索引。后续检查第 b 桶和第 b+1 桶之间是否有簇编号交集：交集非空则簇延续，为空则检查是否有新簇编号从桶内第一次出现。

5.2 事件检测：新增、消失、合并与分裂的判定规则

相邻桶之间同一个簇标号的成员集合用 Jaccard 系数算重叠率，重叠率 >= 0.5 认为簇延续。两个簇标号在相邻桶重叠率都低于 0.5、但又都与下一桶的同一标号重叠，则判定为合并；反之为分裂。这个阈值不需要调得很精确，因为 ST-DBSCAN 的簇成员变化是连续位移的，边界情形极少。

• 新增簇：簇标号 c 在第 b 桶首次出现，且与前驱桶没有任何一个标号的重叠率超过 0.5。
• 消失簇：c 在第 b 桶最后一次出现。
• 合并：两个及以上标号在相邻桶中同时与同一标号存在重叠。
• 分裂：一个标号在下一桶与两个及以上标号存在重叠。

常见误区是直接把簇标号当成稳定标识。每次聚类运行结果里的编号是随机的，进化分析必须基于成员集合的重叠关系，而不是编号本身。

5.3 输出事件表并接到可视化链路

把 cluster_timeline 的输出转换成一张事件表，接上原数据即可支持后续的可视化回放或业务告警。

def build_event_frame(df, labels, timeline, bucket_size): rows = [] t0 = df["timestamp"].min() for b, items in enumerate(timeline): bucket_start = t0 + b * bucket_size for label, idx in items: rows.append({ "bucket": b, "start_time": bucket_start, "label": label, "count": len(items), "lon": df.loc[idx, "longitude"], "lat": df.loc[idx, "latitude"], }) return pd.DataFrame(rows)

生成的 DataFrame 可以直接喂给 Plotly 做时间滑块的散点回放，或用df.groupby(["label", "bucket"]).size()画簇大小热力图。驻留时间分析、热点波动告警、潮汐调度策略验证，都在这张事件表上继续做。

本文还有配套的精品资源，点击获取