5.2 空间数据算法
文档摘要
5.2 空间数据算法 本节摘要：scipy.spatial 处理一切与"点的位置关系"有关的计算。KDTree 用空间分割把最近邻查询从逐点扫描提速到对数级别，是大规模点云应用的加速器；Delaunay 三角剖分与 Voronoi 图给出点集的拓扑结构；pdist 与 cdist 提供欧氏、曼哈顿、余弦等距离度量；ConvexHull 找出点集的边界轮廓。本节四件工具一次讲清，并给出选型对照。 本节导航 阅读完本节，你应当能够： 用 KDTree 完成最近邻查询与半径查询，并解释它为什么比逐点扫描快； 说清 cKDTree 与 KDTree 的实现差异与选型建议； 生成 Delaunay 三角剖分与 Voronoi 图，并说明二者的对偶关系；
5.2 空间数据算法
本节摘要：scipy.spatial 处理一切与"点的位置关系"有关的计算。KDTree 用空间分割把最近邻查询从逐点扫描提速到对数级别，是大规模点云应用的加速器；Delaunay 三角剖分与 Voronoi 图给出点集的拓扑结构；pdist 与 cdist 提供欧氏、曼哈顿、余弦等距离度量；ConvexHull 找出点集的边界轮廓。本节四件工具一次讲清，并给出选型对照。
本节导航
阅读完本节，你应当能够：
- 用 KDTree 完成最近邻查询与半径查询，并解释它为什么比逐点扫描快；
- 说清 cKDTree 与 KDTree 的实现差异与选型建议；
- 生成 Delaunay 三角剖分与 Voronoi 图，并说明二者的对偶关系；
- 在 pdist 与 cdist 的多种距离度量中选出适合当前数据的那一个；
- 用 ConvexHull 计算点集的凸包，并识别退化情形。
一、问题与直觉
先感受一下暴力的代价。给一万个二维点，要找每个点的最近邻，朴素做法是两两算距离，一亿次运算，Python 里跑下来要几十秒。点云规模到十万、百万，这个数字直接不可接受——而激光雷达一帧就是几十万点。
空间数据算法的核心思想就一句话：别把空间当成一维列表，利用点的位置做索引。KDTree 把空间递归切成小块，查询时只搜邻近的小块，大部分点根本不用碰。这一节从最近邻查询讲起，然后是几何结构（三角剖分、Voronoi、凸包），最后是距离度量。四件工具合起来，就是 scipy.spatial 的全貌。
二、核心原理
2.1 KDTree：把空间切成小块
KDTree 的建树过程像切西瓜：取所有点，沿某个坐标轴（第一层 x 轴，第二层 y 轴，交替进行）找中位数，把点集分成左右两半，再对每半递归。查询时从根节点往下走，每次只进一侧子树，找到候选后再回溯检查另一侧——大部分分支根本不用访问。
最小例子：
import numpy as np from scipy.spatial import KDTree rng = np.random.default_rng(0) pts = rng.random((10000, 2)) # 一万个二维点 tree = KDTree(pts) q = np.array([0.5, 0.5]) dist, idx = tree.query(q, k=3) # 最近三个点的距离与下标 print(dist, idx) ball = tree.query_ball_point(q, r=0.05) # 半径内的所有点 print(len(ball))
query 返回最近 k 个点的距离与下标，query_ball_point 返回半径内的点列表。十万点规模下，单次查询是微秒级，比暴力扫描快两三个数量级。为什么这么快？因为每次比较只花一次坐标差运算，就能排除掉一半的点。
2.2 cKDTree 与 KDTree：同一个接口，两种实现
cKDTree 是 KDTree 的 Cython 重写，接口基本一致，速度通常快数倍到数十倍。选型建议很直接：生产环境用 cKDTree，教学和调试用 KDTree——后者的纯 Python 实现更容易读，报错信息也友好些。不过要注意，两个类在部分高级参数上略有出入，老项目里混用前先查版本说明。
图 5-2 KDTree 近邻搜索示意
KDTree 的代价在建树：一次性 O(n log n)。如果你的查询次数很少、点数也少，暴力法反而省事；点多了、查询多了，KDTree 的收益就体现出来了。这是一个典型的"预处理换查询"的取舍。
2.3 Delaunay 三角剖分与 Voronoi 图
点集连成三角形网格，怎么连最合理？Delaunay 剖分的准则是：每个三角形的外接圆内不包含其他点。这个准则保证了三角形尽量饱满，不会出现又扁又长的畸形三角形，对插值、有限元、地形建模都重要。
from scipy.spatial import Delaunay, Voronoi tri = Delaunay(pts[:200]) print(tri.simplices.shape) # 三角形顶点下标矩阵 vor = Voronoi(pts[:200]) print(len(vor.regions)) # Voronoi 区域个数
Delaunay 与 Voronoi 是一对对偶结构：Voronoi 图把平面划分成每个点附近的区域——区域里任意位置到该点的距离比到其他点都近；而 Delaunay 三角形的一条边，恰好连接两个 Voronoi 区域相邻的点。看懂了这对关系，就理解了"谁和谁相邻"的两种等价表达。Voronoi 图在地理上叫势力范围图：基站覆盖、商圈划分、动物领地，都是同一个模型。
2.4 distance 模块：距离度量怎么选
pdist 算一组点两两之间的距离矩阵，cdist 算两组点之间的交叉距离矩阵：
from scipy.spatial.distance import pdist, cdist pair = pdist(pts[:50], metric='euclidean') cross = cdist(pts[:10], pts[:5], metric='manhattan')
metric 参数决定距离的定义。选错度量的后果比选错参数严重得多——曼哈顿距离和欧氏距离在坐标轴上差不了太多，但余弦距离和欧氏距离可能给出完全相反的排序。
| 度量 | 直观含义 | 适用场景 |
|---|---|---|
| euclidean | 直线距离 | 普通物理坐标 |
| manhattan | 沿坐标轴走的距离 | 城市街区、网格路径 |
| chebyshev | 各轴差的最大值 | 棋盘移动、切比雪夫距离 |
| minkowski | 参数化一族，p 值可调 | 需要统一调参时 |
| cosine | 方向夹角，忽略长度 | 文本向量、用户画像 |
| hamming | 对应位不同的比例 | 二进制编码、基因序列 |
欧氏距离在大多数场景够用；数据维度很高时，欧氏距离会失去区分度（维度灾难），此时余弦距离往往更稳。文本与推荐系统里，cosine 是默认选择。
三、工程实践要点
3.1 ConvexHull：点集的边界
凸包是能包住所有点的最小凸多边形。找凸包在图像轮廓、碰撞检测、聚类可视化里常用：
from scipy.spatial import ConvexHull hull = ConvexHull(pts) print(len(hull.vertices)) # 凸包顶点个数
一个需要注意的点：如果所有点几乎共线，凸包会退化——面积趋近于零，后续计算（比如面积、质心）会不稳定。真实数据里这种退化不罕见，处理前先检查 hull.volume 或 hull.area 是否接近零。
3.2 工具选型对照
| 需求 | 工具 | 复杂度 | 备注 |
|---|---|---|---|
| 找最近邻点 | KDTree.query | 对数级 | 点多次查询收益最大 |
| 找半径内所有点 | query_ball_point | 与半径相关 | 返回下标列表 |
| 点集三角网格 | Delaunay | 近线性 | 插值与有限元前置 |
| 势力范围划分 | Voronoi | 近线性 | 与 Delaunay 对偶 |
| 两两距离矩阵 | pdist | 平方级 | 点数大时内存爆炸 |
| 两组交叉距离 | cdist | 平方级 | 常用于匹配问题 |
| 边界轮廓 | ConvexHull | 近线性 | 注意退化情形 |
pdist 的输出是压缩格式，只存上三角，取用时要配 squareform 转换；点数上万时，距离矩阵本身就能吃掉几百兆内存，这时候该想的是降采样或者换 KDTree 思路，而不是硬撑。
⚠️ 常见坑：KDTree 默认用欧氏距离，换度量要用 metric 参数指定，但部分度量下查询效率会下降。另外树在点集变动后必须重建，别指望原地更新——动态点集场景要自己维护重建策略。
💡 关键直觉：KDTree、Delaunay、Voronoi、ConvexHull 是同一件事的四个侧面——把"点的空间关系"结构化。先想清楚要回答什么问题，再选工具，而不是拿着工具找问题。
3.3 cKDTree 与 KDTree：实测说了算
老资料里几乎都写着"cKDTree 比 KDTree 快几十倍"。这个说法在旧版本里没错——那时的 KDTree 是纯 Python 实现，慢是实打实的。但新版本把 KDTree 本身重写成了 Cython 实现，cKDTree 变成了兼容别名。别信传说，跑一下：十万个三维点建树，两个类都在五十毫秒上下；一千次近邻查询，都是四毫秒左右，差距在误差范围内。
import time from scipy.spatial import KDTree, cKDTree rng = np.random.default_rng(0) pts = rng.random((100000, 3)) q = rng.random((1000, 3)) for cls in (KDTree, cKDTree): t0 = time.perf_counter(); tree = cls(pts); build = time.perf_counter() - t0 t0 = time.perf_counter(); tree.query(q, k=5); query = time.perf_counter() - t0 print(cls.__name__, '建树', round(build, 3), '查询', round(query, 3))
整条链路从建树到查询画出来，就是下面这张图：
查询快的原因在上半段：每次比较一次坐标差，就能排除掉一半候选。这份实测还说明一件更普遍的事——网上流传的性能结论有保质期，算法选型前自己跑一遍基准，比背结论可靠。
3.4 Delaunay 剖分的应用：插值与网格生成
Delaunay 剖分最大的工程用途不是画三角形好看，而是给插值和网格生成当地基。scipy.interpolate 的 griddata 在二维以上默认就走这条路线：先查查询点落在哪个三角形里，再用三角形三个顶点的值做加权平均。这套机制对应地形建模里的不规则三角网——测了一批点的高程，用 Delaunay 连成网，任意位置的标高就能估出来。有限元网格生成同样离不开它：剖分质量直接决定求解精度，而 Delaunay 准则恰好保证三角形饱满，避免又扁又长的单元拉垮数值解。想手动确认一个点落在哪个三角形，用 find_simplex：
tri = Delaunay(pts[:200]) simplex = tri.find_simplex(np.array([[0.5, 0.5]])) print(simplex) # 查询点所在的三角形编号
气象插值是最典型的落点：稀疏站点测到的温度、降雨，先做 Delaunay 剖分，再逐三角形插值，整片区域的场就出来了；第 6 章的气象案例里，站点相关性分析也会用到同一套邻近结构。有限元前处理里，剖分质量检查是标准流程——三角形最小角过小、长宽比过大，求解器会直接报警，而 Delaunay 准则在源头就压住了这类问题。
3.5 Voronoi 图：最近服务点的分配模型
Voronoi 图回答的问题是：空间里任意位置，最近的种子点是谁。这个模型在工程里到处都是：基站覆盖范围划分、外卖配送区域、连锁店选址评估、机器人路径规划里的障碍物扩张，本质都是 Voronoi。物理里它叫维格纳赛茨原胞，把晶体切成每个原子周围的一块区域。实现上有个细节要当心：边界点的区域可能是无界的，vor.regions 里出现索引 -1 就表示这个区域延伸到无穷远，画图时要么裁掉，要么把边界框设大。
进阶用法是质心 Voronoi：把种子点迭代地挪到各自区域的质心，多跑几轮，点集就会均匀铺满区域。想做均匀采样或网格优化，这比随机撒点收敛快得多。
3.6 距离度量的更多选择
2.4 的表覆盖了常用度量，工程里还会撞见几个。seuclidean 是标准化欧氏距离，先按各维方差归一化再算，量纲混乱的数据先想到它；mahalanobis 考虑变量间的协方差结构，适合强相关的数据，但要把协方差矩阵的逆传进去，样本少时这个逆矩阵本身就不稳；correlation 把向量中心化后再算余弦，适合有基线漂移的场景；braycurtis 常用于生态学的物种计数。选度量的顺序建议是：先想数据语义——是坐标、是计数、还是方向；再想量纲——各维单位不同就先标准化；最后才翻度量表。高维数据里欧氏距离区分度下降，转余弦；稀疏计数数据用专门度量，别硬套欧氏。
要点速记
- 要点一：KDTree 用递归空间分割把最近邻查询从线性扫描提速到对数级别，适合点多、查询多的场景。
- 要点二：生产环境用 cKDTree，学习调试用 KDTree，接口基本一致。
- 要点三：Delaunay 剖分保证外接圆内无点，三角形饱满，是插值与网格生成的地基。
- 要点四：Voronoi 图与 Delaunay 剖分对偶，前者管势力范围，后者管相邻关系。
- 要点五：pdist 与 cdist 的度量选择要看数据语义，高维数据优先考虑余弦距离。
- 要点六：ConvexHull 找边界，但共线退化会让面积与质心不可靠。
- 要点七：pdist 的距离矩阵是平方级内存，大规模数据要换思路。
下一节换个画风：从空间数据跳到公式里的函数。伽马、误差、贝塞尔，这些名字吓人的函数，scipy.special 一行代码就能算。