简介:面向点云单木分割初学者与需要快速验证欧氏距离聚类流程的开发者,这份资料聚焦PCL中欧氏距离聚类方法的API调用与工程实现,可用于林木点云聚类实验、算法对比和课程示例。压缩包共3个文件,包含1个cpp实现源码、1个xyz点云样例和1个txt聚类输出,约93KB,体积轻便,便于直接查看核心代码、了解输入数据形态与结果组织方式。已有813人学习下载,说明其在点云聚类入门场景中具有一定参考价值。读者可从中获得可运行的欧氏距离聚类实现思路、点云样例与聚类结果对照,并据此理解参数设置、类别划分和输出结构,适合作为单木分割实验的起点或二次开发的基础素材。
1. 欧氏距离聚类:为什么单木分割的第一刀常常砍在这里
“用欧氏距离实现从点云数据中单木分割”这个说法容易造成误导,它对“树”没有做任何建模:既不看树冠形状,也不看分枝结构,更不知道什么叫树干。欧氏距离聚类做的事只有一件——在点云上找 ε-邻域图的连通分量。可它偏偏是几乎所有单木分割流程的第一步,原因不是它准,而是它错得可解释:每一块错误的形状,都能顺着参数一路追溯到具体原因。tolerance 给小了,一棵树被劈成三块;给大了,三棵树粘成一块。这种“参数和结果之间关系明确”的性质,在实测数据上比任何精度指标都珍贵。搞懂它,等于拿到了一把可以校准的尺子。
2. 从点云到簇:欧氏距离聚类到底在算什么
2.1 距离阈值不是树的半径,而是点之间的连通判据
绝大多数人第一次调 tolerance 时的直觉是错的。看到树冠直径 3 米,就把阈值设成 1.5 米,结果稀碎。因为 tolerance 在算法里从来不代表几何尺寸,它代表一个判定:两个点之间的距离小于等于 ε,就在它们之间连一条边。整个点云被这张边连成一张图,算法求的是这张图的连通分量。
这个视角一换,很多现象立刻说得通。一棵树内部的点如果因为遮挡、穿透而稀疏到间距 1.2 米,而 ε 只有 1.0,那么这棵树在图上就是断开的两个分量,聚类必然劈开。反过来,相邻两棵树的树冠如果在下层有少量交叠点,只要存在一条“点间距都小于 ε”的链,两棵树在图上就是同一个连通分量,聚类必然粘连。注意“链”这个词:连通性是有传递性的,A 连 B、B 连 C,A 和 C 距离两米也没关系。这就是欧氏聚类最典型的失效模式,链式效应。
所以 ε 的物理含义不是“树有多大”,而是“点多密”。它必须和点间距同一个量纲去谈,而不是和冠幅同一个量纲。这一点决定了后面所有参数标定的思路。
2.2 队列扩张:一次 BFS 如何把一个簇泡出来
机制层面,欧氏聚类就是一次广度优先搜索加一次半径查询。取一个没访问过的点当种子,把它压进队列,然后不断弹出队首、查询它 ε 邻域内的所有点、把没访问过的邻居压进队列并标记为已访问,直到队列空。这一轮走完,一个连通分量就被“泡”出来了。
下面是一份和主流点云库内部逻辑等价的 Python 实现,写清楚了每个数组的职责:
import numpy as np from scipy.spatial import cKDTree from collections import deque def euclidean_cluster(points, tolerance, min_size=10, max_size=10**9): pts = np.ascontiguousarray(points, dtype=np.float64) n = pts.shape[0] tree = cKDTree(pts) # O(N log N) 建空间索引,只建一次 labels = np.full(n, -1, dtype=np.int32) # -1 表示噪声/被过滤 visited = np.zeros(n, dtype=bool) # 全局访问标记,防止重复入队 cid = 0 for seed in range(n): if visited[seed]: continue visited[seed] = True queue = deque([seed]) members = [] while queue: i = queue.popleft() members.append(i) # 半径查询:返回所有满足 ||p_j - p_i||_2 <= tolerance 的下标 for j in tree.query_ball_point(pts[i], r=tolerance): if not visited[j]: visited[j] = True # 入队即标记,避免同一层内重复压栈 queue.append(j) # 簇大小过滤发生在扩张之后,是后置过滤而非分割条件 if min_size <= len(members) <= max_size: labels[members] = cid cid += 1 return labels, cid代码里有三个参数值得单独说。tolerance是半径查询的 ε,单位与输入点云一致,直接决定连通判据的松紧。min_size和max_size是整型点数,作用在扩张完成之后——小于 min_size 的簇被判定为噪声(标签保持 -1),大于 max_size 的簇被整块丢弃。最后一个细节是visited[j] = True写在入队的位置而不是出队的位置,这一行决定了同一个点不会被重复压进队列,是复杂度的关键。
query_ball_point返回的是半径内的全部邻居,平均返回数量记为 k。整个算法的总代价是 O(N log N + Σk)。极端情况下点全挤在一起,Σk 会退化到 O(N²),此时聚类不再是瓶颈,点云本身的信息密度才是。
2.3 和 DBSCAN 的差别:没有核心点判定导致链式效应
很多人把欧氏聚类当成 DBSCAN 的特例,这个说法只在“都用 ε 邻域”这一层成立,机制上差了一个关键零件。
DBSCAN 有两套点:核心点和边界点。一个点是不是核心点,看它的 ε 邻域内点数是否达到 minPts,只有核心点才能向外扩张。欧氏聚类没有这个概念——任何一个未被访问的点都能当种子,任何一个邻居都会被无条件吸收。区别体现在稀疏区:DBSCAN 会让稀疏区里的孤立点保持噪声身份,欧氏聚类会因为一颗种子顺着链一路爬过去,把整条稀疏链吞进同一个簇。
这解释了单木分割里一个很常见的现象:林分边缘、林窗附近那些点稀疏的树,经常被并到旁边的大簇里,而不是被丢掉。想抑制链式效应,标准做法不是换算法,而是在距离上加各向异性权重,或者先做密度均衡(下采样补齐),把“点间距”这个变量先摁住。
另一个常被忽略的性质是结果的确定性。对固定的 ε,图的连通分量划分是唯一的,点和点的输入顺序只会影响簇的编号顺序,不影响谁和谁在一组。这一点在需要两期数据对比、或者要做结果复现的场景里非常重要——区域生长这类算法对种子顺序敏感,换个遍历顺序结果就变,而聚类不会。
3. 单木分割的前置工序:不做归一化,聚类结果一定翻车
3.1 地面点:把整片林分焊成一个簇的“胶水”
原始点云里最危险的不是树,是地面。激光打到裸地、灌木、枯枝层会返回大量点,这些点在地表上几乎是连续分布的,间距远小于任何合理的 ε。一旦它们进了聚类,就会形成一张覆盖全图的连通网,每一棵树的树干都通过根部连到这张网上,最终结果是一个包含整片林分的巨型簇,而 max_size 一过滤,这个簇被整块丢掉,输出一片空白。
所以地面滤波不是可选项,是前置硬约束。常见的实现思路是先把地面点分离出来(形态学开运算、渐进三角网加密、或者布料模拟这类思路都有人用),再从原始点云里减去,然后计算归一化高度 H = Z − DTM(x, y)。DTM 用地面点做薄板样条或反距离加权插值得到。
这里有个容易混淆的点:纯 3D 欧氏聚类本身对地形起伏并不敏感,因为距离是三维算的,坡地上同一棵树的点对距离不会因为海拔变化而变大。归一化真正的价值在后面——高度切片要用绝对高度,树高统计要用绝对高度,而按高度截断(只保留 H > 2 m 的点)是抑制灌木桥接最便宜的一招。把这两件事分开看,就不会在“到底要不要归一化”上纠结。
3.2 体素下采样:让 tolerance 从“点数”回到“米”
机载和背包平台的点云密度差异极大,同一片林分里,树冠顶部的点间距可能只有 2 厘米,树冠内部因为遮挡能到 30 厘米。ε 是个全局常数,遇到这种非均匀密度,怎么设都是错的:按顶部设,内部处处断裂;按内部设,顶部密不透风地粘连。
体素下采样(voxel grid)是把这个变量摁住的标准手段。把空间切成边长 v 的立方体,每个立方体里所有点用质心代替,输出点云的局部密度就被拉平到大约 1/v³ 的量级,点间距稳定在 v 附近。此时 ε 的取值有了明确锚点:ε 取 2v 到 3v 之间,保证同一棵树的相邻体素质心能连上,又不会跨过树间空隙。
def voxel_downsample(points, voxel): """按体素质心下采样,把非均匀密度拉平到 ~1/voxel^3""" keys = np.floor(points / voxel).astype(np.int64) # 每个点归属的体素索引 uniq, inv, cnt = np.unique(keys, axis=0, return_inverse=True, return_counts=True) out = np.zeros((uniq.shape[0], points.shape[1]), dtype=np.float64) np.add.at(out, inv, points) # 体素内累加 out /= cnt[:, None] # 除以点数得质心 return outvoxel是唯一参数,单位是米,直接对应输出点间距。工程上一般取 0.05 到 0.2 米,取值依据是树冠半径和目标单木的点数预算:voxel 太小,一棵树几万个点,聚类慢且粘连;voxel 太大,小树只剩几十个点,会被 min_size 当噪声清掉。np.add.at这一行是必须的,用普通索引赋值会在同一体素多次写入时丢点,这是血泪经验级的坑。
顺带说一句,如果原始点云带强度或回波次数,体素质心会把这些属性抹平。需要保留的话,得在np.add.at之后按属性再算一次加权平均,或者干脆取体素内离质心最近的那个原始点作为代表。
3.3 高度切片与各向异性缩放:两种把 3D 问题变简单的做法
3D 聚类在单木分割里有个结构性困难:树冠上部枝条稀疏、间隙大,下部树冠互相交叠、间隙小。同一个 ε 不可能同时适配这两层。工程上有两种绕开的办法。
第一种是高度切片。按 1 米左右的层厚把点云切成水平层,每层只取 (x, y) 做 2D 聚类,再把相邻层的簇按水平重叠度或质心距离关联起来,拼成完整的单木。切片的好处是把“树上树下密度不同”变成“每层内部密度相近”,上层切片因为树冠分得开,往往是分割质量最好的一层,可以作为下层合并的锚点。代价是切片边界上的簇容易断,需要一层额外的跨层合并逻辑,而且层厚本身又成了新参数。
第二种更取巧:不改算法,改坐标。欧氏聚类用的是各向同性的球邻域,如果把 Z 轴先乘上一个缩放因子 s,那么在缩放后的空间里做球邻域查询,等价于在原始空间里做椭球邻域查询——水平半径 ε,垂直半径 ε/s。
def anisotropic_cluster(points, height, tol_xy=0.5, tol_z=0.6): """水平容忍 tol_xy、垂直容忍 tol_z 的椭球邻域,复用同一套聚类内核""" s = tol_xy / tol_z # Z 轴缩放因子 scaled = np.column_stack([points[:, 0], points[:, 1], height * s]) return euclidean_cluster(scaled, tolerance=tol_xy)tol_z大于tol_xy时,垂直方向更宽容,同一棵树的上下层更容易连上,而水平方向保持收紧以抑制邻树粘连——这正是单木分割想要的方向性。这个技巧的价值在于它不修改聚类内核,只是喂进去一份被拉伸过的坐标,所有现成的实现都能直接复用。注意缩放只作用于距离计算,输出的簇成员索引仍然对应原始点,标签数组不用回映射。
4. 参数与 API 契约:tolerance、min/max size 怎么定
4.1 tolerance 从点间距反推,而不是从冠幅猜
既然 ε 的物理含义是“点多密”,标定它就应该从点间距出发。做法很直接:对下采样后的点云做 k 近邻查询,取每个点到第 k 个邻居的距离,看这个距离分布的分位数,再乘一个安全系数。
from scipy.spatial import cKDTree import numpy as np def suggest_tolerance(points, k=8, q=0.90, scale=2.5): """按 kNN 距离的高分位数给出初始 tolerance,避免被离群点带偏""" tree = cKDTree(points) d, _ = tree.query(points, k=k + 1) # 第 0 列是点自身,距离为 0 knn = d[:, 1:] # 丢掉自身 return float(np.quantile(knn.mean(axis=1), q)) * scale参数含义:k取 8 到 16 之间,太小会被单个离群点干扰,太大在树冠边缘会被跨树的邻居拉高;q用 0.9 而不是中位数,是因为连通性由最稀疏的那部分点决定,用中位数会低估;scale是安全系数,2 到 3 之间比较稳。这个函数给的是起点而不是终点,真正的取值还得在切片后的上层点云上验证——树冠分开的那一层,簇的数量应该接近目视估计的株数。
如果整片林分密度差异很大,单一 ε 一定不够用,此时可以按空间分块(比如 20 米见方的网格)分别标定 ε,再在块边界上做合并。这个做法会增加复杂度,但在异龄林里往往比反复调全局参数省时间。
4.2 min/max cluster size 是过滤器,不是分割器
这两个参数的定位必须说清楚:它们不参与扩张,只在簇形成之后决定保留还是丢弃。所以它们解决不了“树粘在一起”或“树被劈开”的问题,只能做规模筛选。
估算方法是从点数预算倒推。下采样到体素 v 之后,单位面积的点数密度约为 1/v²,一棵冠幅半径 r 的树,投影面积 πr²,所以点数大致是 πr²/v²。冠幅 2 米、v 取 0.1 米,就是约 1250 个点。min_size取这个估算值的 1/10 左右,能滤掉灌木和噪声;max_size取最大冠幅对应点数的 3 到 5 倍,用来拦住那个“整片林分焊成一块”的巨型簇。
要提醒的是 max_size 一旦触发,那个簇是整块被丢弃,标签保持 -1,而不是被拆分。所以看到输出里出现大片无标签区域,第一反应应该是去看 max_size 是不是被撞了,而不是怀疑聚类本身。
4.3 API 契约:标签语义、确定性、并行边界
把聚类包成 API 时,有几条契约最好在文档里写死,否则调用方一定踩坑。
输入约定为 N×3 的连续内存数组,float64 或 float32 都可以,但坐标单位必须和 tolerance 一致——传进来是毫米、tolerance 传的是米,这个错在跨团队协作里出现频率高得离谱。输出是一个长度 N 的 int32 标签数组,-1表示噪声或未通过规模过滤,其余非负整数是簇编号,编号从 0 开始连续。
确定性这件事值得单独强调:对固定的 ε 和点集,连通分量划分唯一,但簇编号取决于种子的遍历顺序。如果调用方需要跨次运行对齐编号(比如两期数据对比树高生长),就得让编号有稳定的排序依据,最简单的做法是按簇内最小原始索引排序后重新编号,而不是按发现顺序。
并行有明确的边界。扩张阶段本质上是串行的——BFS 的访问标记是全局共享的,多线程直接改visited会出现同一批点被两个线程同时认领的问题,结果虽然还是连通分量,但标签写回时可能互相覆盖。可行的并行化路线是换算法形态:先并行地对每个点做半径查询生成边集,再用并查集(union-find)把这些边合并成簇。这个版本没有全局访问标记的竞争,只在最后归并时加锁,代价是边集可能很大,得按块流式处理。
5. 常见问题排查:五个把单木分割搞坏的坑
5.1 整片林分被聚成一个簇,或者输出全空
现象:聚类返回的簇数量是 1,成员数等于全部点;或者簇数量是 0,输出标签全是 -1。
原因:前者几乎一定是地面点没滤干净,地表那些密集点把每棵树的根部连成了网。后者通常是 max_size 被撞了——那个巨型簇整块被丢弃,于是什么都剩不下。两个现象其实是同一个病的两种表现,取决于 max_size 设得够不够大。
解决:在聚类之前做两件事,一是地面滤波后按归一化高度截断(保留 H > 2 米),二是把 max_size 临时设成无穷大跑一遍,看看最大簇的成员数是多少。这个数字直接告诉你地面点残留的严重程度,比任何日志都直观。
5.2 一棵树被劈成好几块
现象:目视估计 60 株,输出 130 多个簇,且这些簇在空间上互相贴近、高度范围重叠。
原因:局部点间距超过了 ε。树冠内部的遮挡、激光穿透不足、或者下采样过狠(voxel 取得比 tolerance/3 还大)都会造成这问题。链式效应在这里帮了倒忙——原本应该传递连通的路径,因为中间缺了几个点而断掉。
解决:先确认 voxel 和 tolerance 的比例关系,voxel 应该明显小于 tolerance;再检查这棵树所在位置的点密度分位数,如果明显低于全局中位数,说明单一 ε 不适用,改用各向异性缩放(把 tol_z 放大)或者按空间分块标定。不要急着把全局 ε 调大,那会让邻树粘连的坑同时爆发。
5.3 两棵相邻的树粘在一起
现象:输出的簇形状像个哑铃,或者冠幅估算值异常大,接近两棵树的和。
原因:树冠交叠处存在点间距小于 ε 的连通路径,或者下层灌木、藤本把两棵树的树干连了起来。前者是数据本身的分辨率限制,后者是前置滤波不彻底。
解决:优先级最高的是垂直分层,先只对上层(比如 H 大于 8 米的部分)做聚类,把树冠分开的结果作为锚点,再向下扩展;其次是收紧水平方向、放宽垂直方向的各向异性缩放;最后才考虑减小全局 ε,因为那会同时加剧 5.2 的劈裂问题。这里没有银弹,只能按“上层优先、水平收紧”的顺序试。
5.4 两次运行结果不一样,或者上了多线程就乱
现象:同样的输入和参数,两次调用返回的簇编号对不上;或者改成多线程后,某些点的标签在两次运行间跳变。
原因:编号对不上是正常的,连通分量划分唯一但发现顺序依赖种子遍历。标签跳变则是真 bug——多线程并发写 labels 或 visited 造成了竞争。
解决:把“划分唯一”和“编号稳定”当成两件事处理,编号在最后统一按簇内最小索引重排。并行版本彻底放弃共享 visited 的设计,改成边集加并查集。验证方法很土但有效:把输入点云随机打乱顺序跑两遍,比较两次的划分结果(按集合比较,不比编号),如果不一致,说明实现里混进了顺序依赖。
5.5 边缘和下层的小树全被当噪声丢掉
现象:林分中心区域分割正常,边缘和下层更新层的树完全消失,标签为 -1。
原因:min_size 是按主林层大树的点数预算设的,边缘树因为扫描角度、遮挡、或者本身冠幅小,点数只有大树的十分之一,直接被规模过滤清掉。另外,如果做了体素下采样,小树可能只剩几十个点,同样过不了门槛。
解决:min_size 不要用全局常数,改成和局部密度挂钩——按空间分块统计点密度,每块的 min_size 按该块密度重新估算。另一个办法是把过滤拆成两步:先按绝对下限(比如 20 个点)清掉纯噪声,再把剩余簇按点数和冠幅的比值做二次筛选,避免误伤。
6. 设计取舍:欧氏聚类该在哪一步退场
欧氏聚类的边界非常清晰:它只回答“这些点在几何上连不连”,不回答“这是不是一棵树”。凡是需要“树”这个概念的地方——判断两个簇是不是同一棵树的不同部分、判断一个簇是不是被压在下层的更新木——它都无能为力。认清这条边界,比调参重要得多。
| 方法 | 依赖的先验 | 结果确定性 | 主要失效模式 | 适合的位置 |
|---|---|---|---|---|
| 欧氏距离聚类 | 点间距、ε | 划分唯一,与输入顺序无关 | 密度不均导致断裂或粘连 | 流程第一步,粗分割 |
| 区域生长 | 种子点、相似性准则 | 依赖种子顺序 | 种子选错则整片污染 | 有先验锚点时替代聚类 |
| 冠层高度模型局部最大值 | 冠层表面、窗口尺寸 | 确定 | 平顶树冠漏检、下层树不可见 | 主林层株数估计 |
| 有监督分割 | 标注样本 | 依赖模型与权重 | 跨区域泛化差 | 有标注预算时的精修 |
表格里最值得琢磨的是“确定性”这一列。欧氏聚类是唯一一个结果与输入顺序无关的方法,这个性质在需要长期监测、两期对比的场景里价值极高——你不需要解释“为什么这次多出三棵树”,因为划分是唯一的,多出来的三棵一定对应数据或参数的真实变化。
另一个容易忽略的取舍是“什么时候该加复杂度”。各向异性缩放几乎零成本,只是把 Z 轴乘一个因子,就能显著缓解上下层密度差异,性价比极高,建议默认开启。高度切片加跨层合并能提升质量,但引入了层厚、重叠度阈值、合并策略三个新参数,只有当单层聚类的质量明显不够时才值得上。至于换用带监督的模型,那是标注预算的问题,不是算法优劣的问题。
我自己的习惯是每次换数据源都先做一遍“裸聚类”——不做切片、不做缩放、只用建议的 tolerance 跑一遍,把结果按高度着色看一眼。这一步花不到五分钟,但能立刻暴露出地面残留、密度断层、坐标系单位错误这三类问题。有一次拿到一份坐标单位是毫米的数据,tolerance 按米填的 0.5,聚类返回零个簇,排查了半小时才发现是单位的事。从那以后我每次接入新数据都强制先打印点云的包围盒尺寸和 kNN 距离中位数,两个数字一对,量纲错了立刻现形。希望帮到你。
本文还有配套的精品资源,点击获取