第15章 - 点云分割、聚类与特征提取

第 10 章把点云分成”地面/非地面”，那是二分类。真实项目很快会遇到更细的问题：一片建筑区里有多少栋独立建筑（聚类）、每栋楼的立面和屋顶怎么区分（几何特征）、小区范围怎么用一个多边形圈出来（边界提取）。这些问题不是”再加一个滤波器”能解决的，它们都需要先定义邻域，再在邻域上算某种量。

PDAL 在这块的 stage 可以分成三组：

本章的实测数据来自 PDAL 2.10.0：用 readers.faux 生成的合成点云，跑通 filters.cluster + filters.covariancefeatures + filters.hexbin，把真实列名、默认值和报错原文记录下来。

15.1 先理解邻域

15.1.1 knn 与 radius

几乎所有特征与聚类都归结为一个问题：对每个点，找它的邻居。PDAL 提供两种方式：

filters.eigenvalues 与 filters.normal 的 knn 默认 8，filters.covariancefeatures 的 knn 默认 10。半径模式要求指定 min_k（默认 3），邻居不足时跳过该点。

15.1.2 用最优邻域代替固定 knn

固定 knn 在不同密度区域表现不一致。filters.optimalneighborhood 为每个点搜索使特征最稳定的邻域，输出两个维度：

• OptimalKNN：最优邻居数；
• OptimalRadius：对应半径。

参数只有两个：k_min 默认 10，k_max 默认 14。之后 filters.covariancefeatures 设 optimized: true 就会读这两个维度，不再用固定 knn。代价是额外的搜索开销，但在地形起伏、植被与建筑混合的场景里，分类边界会明显更干净。

15.2 协方差特征

15.2.1 特征从哪来

对邻域内点的坐标求 3×3 协方差矩阵，做特征值分解得到 λ₀ ≥ λ₁ ≥ λ₂，再组合成各种描述子。filters.covariancefeatures 支持 11 个：

前四个（Linearity/Planarity/Scattering/Verticality）构成所谓 “Dimensionality” 预设，也是 feature_set 的默认值。

15.2.2 参数与算例

一个实测 pipeline（去掉 Density，理由见 15.2.3）：

[
  {
    "type": "readers.faux", "filename": "faux", "count": 500,
    "mode": "Normal", "mean_x": 0, "mean_y": 0, "mean_z": 0,
    "stdev_x": 5, "stdev_y": 5, "stdev_z": 1, "seed": 7
  },
  { "type": "filters.optimalneighborhood" },
  {
    "type": "filters.covariancefeatures",
    "knn": 8, "feature_set": "all", "optimized": true, "threads": 1
  },
  { "type": "writers.text", "filename": "clus.csv", "format": "csv" }
]

输出列名（PDAL 2.10.0 实测）：

"X","Y","Z","OffsetTime","OptimalKNN","OptimalRadius","Linearity","Planarity",
"Scattering","Verticality","Omnivariance","Anisotropy","Eigenentropy",
"EigenvalueSum","SurfaceVariation","DemantkeVerticality","Density"

注意 Density 这一列只在跑了 optimalneighborhood 之后才出现——这正是下一个坑。

15.2.3 坑：Density 需要前置 stage

直接把 feature_set 设成 "all" 而不跑 optimalneighborhood，会得到（实测原文）：

PDAL: filters.covariancefeatures: Density feature requires OptimalKNN and OptimalRadius
dimensions, which are missing in the input PointView.

两种修法：一是显式列出不依赖最优邻域的特征，"Linearity,Planarity,Scattering,Verticality" 或直接用默认 "Dimensionality"；二是先加 filters.optimalneighborhood 并把 optimized 设为 true。第二种是官方推荐路线，因为 optimalneighborhood 的输出本身就是后续分类的好特征。

15.2.4 eigenvalues 与 normal 的区别

filters.eigenvalues 输出三个原始特征值维度 Eigenvalue0/1/2（参数只有 knn 默认 8 和 normalize 默认 false），适合你自己组合特征；filters.normal 输出 NormalX/NormalY/NormalZ 与 Curvature，适合做表面朝向分析、法向量一致性检查。filters.normal 还有 always_up（默认 true，强制法向量朝上）和 refine（默认 false，把法向量传播到未计算点）。

15.3 聚类：四种算法怎么选

15.3.1 filters.cluster：欧氏距离聚类

先做欧氏聚类提取（Euclidean Cluster Extraction），再给每个点打标签。规则简单：两点距离 ≤ tolerance 就归为一类，再按点数过滤。

输出维度 ClusterID：从 1 开始编号，不属任何簇的点为 0。

[
  "input.las",
  { "type": "filters.cluster", "tolerance": 2.0, "min_points": 50 },
  { "type": "writers.las", "filename": "cluster.las", "extra_dims": "ClusterID" }
]

tolerance 是唯一需要调的量。点间距 0.1 m 的机载数据，建筑之间距离通常 5 m 以上，取 1.0 就够；地面激光扫描点间距 2 cm 且场景连续，取 1.0 会把整片区域连成一簇，这时要按实际空隙调到 0.1–0.3。

15.3.2 filters.dbscan：密度聚类

DBSCAN 不要求簇的形状，能在噪声里找出任意形状的簇。参数：

与 filters.cluster 的区别：DBSCAN 有”噪声点”概念，密度不足的点不进任何簇；欧氏聚类只要连通就算一类。数据里有大量孤立噪声时，DBSCAN 的结果更干净。

15.3.3 filters.lloydkmeans：按 K 均值切

Lloyd K-means 适合”已知要分成几类”的场景（如按高度分带、按强度分级）。参数只有 clusters（默认 10）、max_iterations（默认 10）和 dimensions（默认 X,Y,Z）。它不做几何连通性判断，簇可能在空间上不连续——这是特性不是缺陷，但用它做建筑单体提取会得到奇怪的结果。

15.3.4 supervoxel 与 litree

filters.supervoxel 做超体素分割，输出 ClusterID，参数 seed_resolution 默认 1.0、spatial_importance、normal_importance 等，适合先过分割再合并的分割流水线。filters.litree 是树木单株分割：minimum_points 默认 10、minimum_height 默认 3.0、search_radius 默认 100.0，输入要求有 HeightAboveGround 维度（先用 filters.hag_nn/filters.hag_dem 算）。

15.4 边界与密度：hexbin

15.4.1 它解决什么问题

给一个测区点云计算”数据覆盖范围”的多边形，最直接的做法是算凸包，但点云常有凹口和空洞。filters.hexbin 的思路是：把区域切成六边形格网，统计每格点数，保留点数超过阈值的格子，再把这些格子并起来求边界。

15.4.2 实测输出

[
  {
    "type": "readers.faux", "filename": "faux", "count": 50000,
    "mode": "Uniform", "bounds": "([0,100],[0,100],[0,5])", "seed": 11
  },
  {
    "type": "filters.hexbin",
    "density": "hex_density.geojson",
    "boundary": "hex_boundary.geojson",
    "edge_length": 10.0, "threshold": 15
  }
]

100×100 的均匀点云、边长 10 的六边形，产出 56 个密度格与 1 个边界多边形；每个密度 feature 的属性只有两个字段：

{ "ID": 8589934593, "COUNT": 1292 }

ID 是 H3 风格的整数格 ID（h3_grid 关闭时也仍用这个字段名），COUNT 是落在该格内的点数。boundary 输出的是 1 个多边形 feature，可直接给 filters.crop 做后续裁剪，或作为成果范围写进元数据。

阈值 threshold 直接决定边界形状：调大，边缘零散点被排除，边界收紧；调小，边界外扩。用默认 15 时，低密度边缘区域通常会被削掉一圈，这是符合预期的。

15.5 一个完整的建筑提取流程

把上面的 stage 串起来，从原始点云到建筑单体：

[
  "input.las",
  { "type": "filters.expression", "expression": "Classification == 6" },
  { "type": "filters.optimalneighborhood", "k_min": 10, "k_max": 14 },
  {
    "type": "filters.covariancefeatures",
    "optimized": true, "feature_set": "Linearity,Planarity,Scattering,Verticality",
    "threads": 4
  },
  { "type": "filters.expression", "expression": "Planarity > 0.6 && Verticality < 0.4" },
  { "type": "filters.cluster", "tolerance": 1.5, "min_points": 200 },
  { "type": "filters.expression", "expression": "ClusterID > 0" },
  {
    "type": "writers.las", "filename": "buildings.las",
    "extra_dims": "ClusterID,Planarity,Verticality", "forward": "header"
  }
]

每一步的作用：先只用已标为建筑的点（LAS 分类 6）；用最优邻域算特征；用 Planarity/Verticality 剔除植被和车辆；聚类；丢掉噪声簇；写出带 ClusterID 的结果。实际项目里 Planarity/Verticality 的阈值要按数据源标定，0.6/0.4 只是起点。

15.6 排错清单

15.7 生产建议

先定尺度，再选参数。 knn、tolerance、edge_length 都是长度量纲，必须由数据点间距推导。入库前先跑 pdal info --stats 拿到 X/Y/Z 的 count 与范围，估算平均点间距，再设参数——照着别人的参数表抄是最常见的错。

特征计算放在分类之前。 协方差特征是”局部形状”的客观描述，可以先用它们做无监督筛查，再决定分类策略。把特征计算放在 pipeline 前段，也能让后续多个分支复用同一份结果。

聚类结果一定要保留 ClusterID。 写成 LAS 时用 extra_dims 带上，否则下游无法按簇聚合。对 BPF/Parquet 这类支持任意列的格式更省事。

边界数据与点云分开管理。 filters.hexbin 产出的 GeoJSON 多边形本质上是”数据覆盖范围”，应该跟着元数据走（入库 PostGIS、写进 STAC item 的 geometry），而不是每次分析重新算。

15.8 本章小结