带Classification的LAS跑通PointNet++点云语义分割

简介：面向需要将 PointNet 或 PointNet++ 落地到自有 LAS 点云的研究者与工程人员，这份资料聚焦带 Classification 属性的 LAS 点云如何组织训练与预测，基于 yanx27 的 PyTorch 实现整理，可在两种网络结构间按需选择，适合已有一定深度学习基础、希望跑通自定义点云分类或分割流程的读者。压缩包共 19 个文件、约 26KB，主体为 15 个 Python 脚本，覆盖网络模型、数据加载、LAS 数据处理及训练测试入口；另附一个已分类 LAS 文件用于训练与验证、一个 LAS 文件用于预测，以及一份说明代码结构的 docx 文档，便于对照目录快速定位模块。已有 896 人学习。读者可据此理解从原始 LAS 到数据集构建、再到模型训练与推理的完整链路，并借助现成脚本减少环境配置与数据接口改造中的反复试错，快速形成可复现的自定义点云实验基线。

1. 从一份带 Classification 的 LAS 说起：PointNet++ 语义分割落地要跨过哪几道坎

手上拿到一份 LAS，每条记录里都带着 Classification 字段，看起来标签已经现成，直接丢进 PointNet / PointNet++ 就能跑点云语义分割。真动手就会发现中间隔着好几道坎：LAS 里的分类码是测绘体系的编码，地面是 2、建筑是 6、植被是 3/4/5，不是从 0 开始的连续类别；一整幅点云动辄几百万点，显存根本塞不下；点云没有网格结构，分块、采样、归一化任何一步训练和推理不一致，结果就会翻车。这篇笔记按 pytorch 点云语义分割的实际流程走一遍：从 LAS 读 Classification、做类别映射和分块采样、改 PointNet++ 的分割头、训练评估，到最后把预测结果写回 LAS。适合手里有带分类码的机载或车载点云、想拿 PointNet 系列做一个能跑通的 baseline 的从业者。

2. 拆开 LAS 的 Classification：用 laspy 把点云读成可训练的数组

2.1 LAS 的坐标、尺度与 Classification 字段语义

LAS 文件分成 header 和 point records 两块。header 里最该盯的是版本号、point format、scale/offset 和点数，point records 里才是每个点的坐标和属性。坐标在文件里存的是整数，读出来要按x = X_record * scale_x + offset_x还原成浮点。很多人第一次读 LAS 觉得坐标“怎么是几百万”，就是因为只看了整数部分没管 offset。laspy 的las.x已经替你做了这层换算，直接拿来用就行，但 header 里的 scales 和 offsets 建议打印出来看一眼，机载数据的 scale 常见是 0.001，offset 一般是大数，丢了它就等于把整幅点云平移到了别处。

Classification 在 LAS 1.0 到 1.3 里就是一个 uint8 字段，取值遵循 ASPRS 的分类码约定；LAS 1.4 的 point format 6 以上把分类标志位（synthetic、key-point、withheld、overlap）拆成了独立字节，Classification 本身仍然是 uint8。这意味着你拿到的分类码大概率是测绘软件写进去的语义码，而不是神经网络要的 0 起始连续标签。下表是机载点云里最常出现的几个码，做映射时基本绕不开：

提示：分类码 0 和 1 在多数开源数据集里都是“没标”的意思，把它们当正常类别训，模型会学到一个语义混杂的类，mIoU 看着不掉但可视化会很难看。

2.2 先统计再动手：读取与类别分布

不要一上来就写转换脚本。第一件事是把文件读进来，看点数和类别直方图，这一步能提前暴露 80% 的坑。

import numpy as np import laspy # 读入单幅 LAS，laspy 新版把子字段直接挂在对象上 las = laspy.read("tile_01.las") hdr = las.header print("version:", hdr.version, "point_format:", hdr.point_format.id) print("point_count:", hdr.point_count) print("scales:", hdr.scales, "offsets:", hdr.offsets) # x/y/z 已经按 scale/offset 还原为浮点，直接堆成 N×3 xyz = np.stack([las.x, las.y, las.z], axis=1).astype(np.float32) cls = np.asarray(las.classification).astype(np.uint8) print("xyz range:", xyz.min(0), xyz.max(0)) uniq, cnt = np.unique(cls, return_counts=True) for u, c in zip(uniq, cnt): print(f"class {u:3d}: {c:10d} {c / len(cls) * 100:6.2f}%")

逻辑很直白：先把 header 里的版本、点格式、缩放参数打出来，再用las.classification拿到分类数组做直方图。参数上唯一要注意的是 dtype，las.classification读出来是 uint8，转 numpy 时别再降精度；xyz必须转 float32，float64 会让后面的 Dataset 内存翻倍。跑完这段你会得到两个关键信息：类别码都有哪些、地面类占比多少。地面占到 70% 以上是常态，这直接决定了后面损失函数和采样策略怎么设。

2.3 类别码映射：用 LUT 把测绘编码压成 0..N-1

网络输出的是 0 到 N-1 的连续类别，LAS 给的是稀疏的分类码，中间必须有一张映射表。我一般写成字典，再展开成 256 长度的查找表，比np.vectorize快一个数量级，而且反查也方便。

import numpy as np # 左边是 LAS 分类码，右边是网络类别 id；没列出的码统一落到 ignore CLASS_MAP = {2: 0, 6: 1, 3: 2, 4: 2, 5: 2, 9: 3} IGNORE = 255 lut = np.full(256, IGNORE, dtype=np.uint8) for src, dst in CLASS_MAP.items(): lut[src] = dst label = lut[cls] # 映射后的标签，取值 0..3 或 255 print(np.unique(label, return_counts=True)) # 反查表：把网络预测的类别 id 还原成 LAS 分类码 inv_lut = np.zeros(256, dtype=np.uint8) for src, dst in CLASS_MAP.items(): inv_lut[dst] = src

这段的核心是lut和inv_lut两张表。lut用分类码直接索引，256 项覆盖 uint8 全空间，没配的码默认 255，天然就是ignore_index，不需要额外写 if。inv_lut是写回 LAS 时的后悔药——网络输出 1，你得知道它对应 LAS 的 6 而不是原样写 1。参数上要留意一个细节：多个源码映射到同一个目标类（3/4/5 都映射到 2）时，反查只能还原成其中一个码，写回前想清楚这是否可接受；如果下游要区分植被高度，就得在映射阶段就拆开。

2.4 落盘：npy 存坐标、标签和原始索引

转换脚本跑一次要几分钟到几十分钟，别每次都重读 LAS。我习惯把一幅图拆成三个数组存：坐标、标签、原始点序号。第三个最容易被忽略，但写回预测结果时全靠它。

import numpy as np import laspy las = laspy.read("tile_01.las") xyz = np.stack([las.x, las.y, las.z], axis=1).astype(np.float32) cls = np.asarray(las.classification).astype(np.uint8) lut = np.full(256, 255, dtype=np.uint8) for src, dst in {2: 0, 6: 1, 3: 2, 4: 2, 5: 2, 9: 3}.items(): lut[src] = dst np.save("tile_01_xyz.npy", xyz) np.save("tile_01_label.npy", lut[cls]) np.save("tile_01_idx.npy", np.arange(len(xyz), dtype=np.int64))

idx就是 0 到 N-1 的自然序列，看着没用，但分块采样会打乱点序，推理阶段要靠它把每个点的预测投回原始位置。坐标用 float32、标签用 uint8、索引用 int64，是内存和精度之间比较划算的组合。文件大到内存放不下时，换成np.memmap或者按 tile 分文件存，不要试图一次性读完整测区。

3. 把整幅点云切成训练样本：分块、采样与标签对齐

3.1 为什么必须分块：显存、感受野与边界效应

PointNet++ 的 set abstraction 层做的是球查询加最远点采样，计算量随点数呈超线性增长。一幅机载点云几百万点，直接喂进去不是显存爆就是跑一晚上一个 epoch。分块是刚需，但块大小不是随便定的：块太小，单块里可能只有地面，模型学不到建筑和植被的上下文；块太大，显存吃不住，而且最远点采样在超大点集上会漏掉小目标。工程上常用的起点是水平块 2 到 5 米、每块采 4096 个点，机载数据可以放宽到 8192，车载或室内数据往 1 米、2048 点收。

另一个容易被忽视的是边界效应。滑窗切块时，落在块边缘的建筑会被切断，模型在边缘学到的形状是残缺的。训练阶段我一般用不重叠切块（stride 等于 block），推理阶段才用 50% 重叠加投票，两者分开处理，别用同一套参数硬套。

3.2 滑窗分块与点数补齐的实现

分块函数要解决三件事：按平面网格取点、点数不足时补齐、点数过多时降采样。补齐用有放回采样比补零好，补零会让网络学到一堆位于原点的假点。

import numpy as np def make_blocks(xyz, label, idx, block=2.0, stride=2.0, num_point=4096, min_pts=512, seed=0): rng = np.random.default_rng(seed) xmin, ymin = xyz[:, 0].min(), xyz[:, 1].min() xmax, ymax = xyz[:, 0].max(), xyz[:, 1].max() xs = np.arange(xmin, xmax, stride) ys = np.arange(ymin, ymax, stride) blocks = [] for bx in xs: for by in ys: m = ((xyz[:, 0] >= bx) & (xyz[:, 0] < bx + block) & (xyz[:, 1] >= by) & (xyz[:, 1] < by + block)) n = int(m.sum()) if n < min_pts: # 点太少，放弃这一块 continue sel = np.nonzero(m)[0] if n >= num_point: sel = rng.choice(sel, num_point, replace=False) # 降采样 else: extra = rng.choice(sel, num_point - n, replace=True) # 有放回补齐 sel = np.concatenate([sel, extra]) blocks.append((xyz[sel], label[sel], idx[sel])) return blocks

参数含义逐个说清楚：block是水平块边长，单位跟坐标系一致，投影坐标下就是米；stride是滑窗步长，训练时设成和block相等做不重叠切块，推理时设成block/2换取重叠；num_point决定单样本点数，直接对应显存；min_pts是块的最小点数阈值，低于它就跳过，避免大量只有几十个点的垃圾块进入训练集。seed固定住能保证每次转换结果一致，调参排查时很省事。

注意：这里只按 x/y 切块，z 方向不切。机载数据 z 跨度大，按三维切会产生大量空块；但如果你的数据是高层建筑立面加地面，块内高差可能超过 100 米，归一化时要留意。

3.3 归一化：绝对坐标是训练里的隐形炸弹

投影坐标的数值动辄几十万，直接送进网络，第一层 MLP 的激活值就炸了。归一化必须做，而且训练和推理必须用同一个函数，这是血泪经验。

import numpy as np def normalize(xyz): # 先平移到块质心，再缩放到单位球 center = xyz.mean(axis=0, keepdims=True) p = xyz - center scale = np.max(np.linalg.norm(p, axis=1)) return (p / max(scale, 1e-6)).astype(np.float32)

这段做的是“去质心 + 球归一化”。先用块的质心做平移，消掉绝对坐标；再除以块内点到质心的最大距离，把所有点压进单位球。为什么不逐轴 min-max？因为那会改变点云的几何比例，一个狭长的建筑块会被拉成方块，法向和曲率全变了。为什么不除以固定常数（比如 50 米）？因为不同块的尺度差异很大，固定常数在稀疏块上会把点压成一团。

有个边界情况必须提：如果块里只有地面，点近似共面，scale可能很小，除完会放大噪声点的扰动。max(scale, 1e-6)只是防除零，真正的解法是在分块阶段就把纯地面块按比例丢弃，或者给scale设一个下限。

3.4 训练/验证划分：按区域切，不要按块随机切

按块随机划分是最常见的错误之一。同一幅点云相邻块之间有大量重叠点，随机切会让验证集里的点出现在训练集中，指标虚高得离谱，上线就现原形。正确做法是按 tile 或按空间区域划分：拿几幅独立的 LAS 做验证，或者在同一幅图里按坐标范围切出验证区。

类别平衡同样要处理。地面占 70% 时，模型把所有点预测成地面就能拿到很高的 OA，但 IoU 惨不忍睹。两个手段配合用：一是在分块阶段对含稀有类的块做上采样（建筑块重复采样两次），二是损失函数加类别权重，权重按1 / log(1.02 + freq)这类公式算，别用 1/freq，容易把稀有类权重拉到几百倍导致训练不稳。

4. 模型侧改造：把 PointNet/PointNet++ 的类别数对齐到你的标注体系

4.1 分割头与 num_classes、ignore_index 的对应关系

PointNet 系列做分割时，分类分支和分割分支共享特征提取，分割头的输出通道数必须等于你的类别数。原始实现里num_classes通常是 13 或 16（对应公开数据集），换成自己的数据必须改，改完还要检查两处：最后一层 conv 的 out_channels，以及损失函数里的ignore_index。

import torch import torch.nn as nn NUM_CLASSES = 4 # 对应 CLASS_MAP 里的 0..3 IGNORE = 255 # 分割头：输入是 (B, C, N)，输出 (B, num_classes, N) seg_head = nn.Sequential( nn.Conv1d(128, 128, 1), nn.BatchNorm1d(128), nn.ReLU(inplace=True), nn.Dropout(0.5), nn.Conv1d(128, NUM_CLASSES, 1), ) criterion = nn.CrossEntropyLoss(ignore_index=IGNORE)

NUM_CLASSES必须和CLASS_MAP的值域上限加一严格一致，写成 4 就意味着映射后的标签只能是 0/1/2/3/255。ignore_index=255让被标为 ignore 的点不参与梯度，这是处理未分类点的标准做法。Dropout(0.5)放在分割头前是原始实现的配置，小数据集上可以保留，数据量大时降到 0.3 更稳。

提示：PyTorch 的CrossEntropyLoss默认期望输入是(B, C, d1, d2, ...)，点云模型输出通常是(B, N, C)，调用前记得logits.transpose(1, 2)，这个顺序错了不会报错，但训练结果就是玄学。

4.2 Dataset 与 DataLoader：把归一化锁在一个函数里

Dataset 只做两件事：取块、归一化。所有跟坐标变换相关的逻辑都收在这里，训练和推理共用同一个类，杜绝“训练用 A 方法、推理用 B 方法”的低级错误。

import numpy as np import torch from torch.utils.data import Dataset def normalize(xyz): center = xyz.mean(axis=0, keepdims=True) p = xyz - center scale = np.max(np.linalg.norm(p, axis=1)) return (p / max(scale, 1e-6)).astype(np.float32) class BlockDataset(Dataset): def __init__(self, blocks, augment=False): self.blocks = blocks self.augment = augment def __len__(self): return len(self.blocks) def __getitem__(self, i): xyz, label, idx = self.blocks[i] if self.augment: # 绕 z 轴随机旋转 + 小幅抖动，是最便宜的增广 theta = np.random.uniform(0, 2 * np.pi) c, s = np.cos(theta), np.sin(theta) rot = np.array([[c, -s, 0], [s, c, 0], [0, 0, 1]], dtype=np.float32) xyz = xyz @ rot.T xyz = normalize(xyz) return (torch.from_numpy(xyz), torch.from_numpy(label.astype(np.int64)), torch.from_numpy(idx.astype(np.int64)))

坐标、标签、索引三件套一起返回，idx在训练时用不上，但推理阶段必须靠它把预测投回原始点序。增广放在归一化之前做旋转，是因为旋转不改变尺度，顺序反过来结果一样但读起来更顺。label转 int64 是CrossEntropyLoss的硬性要求，忘了转会在第一个 batch 就报 dtype 错误。DataLoader 的num_workers设成 CPU 核数的 1/4 到 1/2，设太大反而因为进程切换变慢。

4.3 训练循环：lr、batch、epoch 的实际取值

超参没有万能值，但有一组能跑通的起点：Adam，初始 lr 1e-3，weight decay 1e-4，batch 8 到 16，epoch 80 到 120，学习率用多项式衰减 power 0.9。

import torch from torch.optim.lr_scheduler import PolynomialLR device = torch.device("cuda" if torch.cuda.is_available() else "cpu") model = PointNet2Seg(num_classes=4).to(device) # 换成你实际的模型类 opt = torch.optim.Adam(model.parameters(), lr=1e-3, weight_decay=1e-4) crit = torch.nn.CrossEntropyLoss(ignore_index=255) sched = PolynomialLR(opt, total_iters=100, power=0.9) for epoch in range(100): model.train() for xyz, label, _ in train_loader: xyz, label = xyz.to(device), label.to(device) logits = model(xyz.transpose(1, 2)) # (B, N, C) -> (B, C, N) loss = crit(logits, label) opt.zero_grad() loss.backward() torch.nn.utils.clip_grad_norm_(model.parameters(), 1.0) opt.step() sched.step()

transpose(1, 2)那行是整段最容易翻车的地方，前面提过。clip_grad_norm_设 1.0 是防 PointNet++ 在稀疏块上梯度爆炸的保险，去掉也能跑，但偶尔会碰到 loss 变 NaN。batch 大小主要受num_point约束：4096 点时 batch 16 大概吃 8 到 12 GB 显存，显存不够先降 batch 再降点数，别反过来，点数降太狠会破坏局部结构。

4.4 评估：混淆矩阵算 mIoU 与 OA

语义分割不能只看 loss。每个 epoch 结束在验证集上累计混淆矩阵，最后一次性算 IoU 和 OA，比逐 batch 平均靠谱。

import numpy as np def update_cm(cm, pred, gt, num_classes, ignore=255): m = gt != ignore idx = gt[m] * num_classes + pred[m] cm += np.bincount(idx, minlength=num_classes ** 2).reshape(num_classes, num_classes) return cm def compute_metrics(cm): tp = np.diag(cm).astype(np.float64) fp = cm.sum(axis=0) - tp fn = cm.sum(axis=1) - tp iou = tp / np.maximum(tp + fp + fn, 1e-9) oa = tp.sum() / np.maximum(cm.sum(), 1e-9) return iou, iou.mean(), oa

update_cm把预测和真值压成一个整数索引再bincount，比双重循环快几个数量级，768 类以内都不会有性能问题。compute_metrics返回逐类 IoU、mIoU 和整体准确率。看结果时别只盯 mIoU，逐类 IoU 里出现 0 的类要单独查——通常是这一类在验证集里样本太少，或者映射表里根本没配。

5. 避坑清单：LAS 标签训练里最容易翻车的五件事

5.1 现象：loss 卡在某个值不降，预测全是同一类

原因通常不是模型不行，而是标签被大面积 ignore 掉了。CLASS_MAP里只写了几个码，剩下的全落进 255，如果数据来自不同测绘软件，分类码体系可能完全不一样，实际有效标签只剩 5%。另一种情况是稀有类占比太低，梯度被地面类淹没。

解决办法：转换前先跑 2.2 的直方图，把类别分布打出来；训练第一个 epoch 后统计 batch 里非 ignore 点的比例，低于 60% 就回去改映射表。稀有类用分块上采样加类别权重双管齐下，权重上限压到 20 倍以内。

5.2 现象：写回 LAS 后点云错位，或者整片变成一个颜色

原因是分块采样打乱了点序，写回时按预测数组的顺序往原始文件里塞。分块时有放回补齐还会产生重复点，一个原始点可能对应多条预测。

解决办法：分块时把idx一起带出来，推理时按idx把概率累加回原始数组，最后对每个原始点取 argmax。这一步没有捷径，idx是唯一的后悔药，转换阶段就必须存下来。

5.3 现象：验证集 mIoU 高得离谱，换一块数据就崩

原因是训练集和验证集存在空间泄漏。随机按块划分时，相邻块共享大量边界点，模型等于在背答案。

解决办法：按 tile 或按坐标范围划分验证区，保证验证区的点从未在训练中出现。划分完抽查一下两边的 x/y 范围有没有交叠，交叠了就是切错了。

5.4 现象：训练正常，推理时同一条航带的数据预测全乱

原因是归一化方式不一致。训练时在块内做去质心加球归一化，推理时图省事用了全局 min-max，同一栋建筑在两种变换下的尺度差了好几倍。

解决办法：把normalize抽成独立模块，Dataset 和推理脚本 import 同一个函数，禁止各写一份。推理时依然按块处理，每块单独归一化，不要试图对整幅图归一化后再切。

5.5 现象：预测结果里出现了根本没定义的类别号

原因是映射表和反查表方向搞反了，或者写回时直接把网络输出的 0..N-1 当成 LAS 分类码写进去了。LAS 里 0 是“未分类”、1 是“未指定”，网络类别 1 写进去就变成了未指定。

解决办法：写回前强制走一遍inv_lut，写完再读回来做一次np.unique检查，确保所有值都在CLASS_MAP的键集合里。这一步我一般写成断言，不在集合里直接抛异常，别让它静默通过。

6. 从预测回到 LAS：重叠推理、投票与写回自查

训练跑通只是半程，真正交付的是带预测分类码的 LAS 文件。推理阶段和训练阶段有两个关键差异：块要重叠、结果要投票。用 50% 重叠滑窗，每个点会被预测 2 到 4 次，把这些预测的概率累加起来再取 argmax，比单次预测稳得多，尤其在块的边界处。

import numpy as np import torch def predict_tile(model, xyz, block=2.0, stride=1.0, num_point=4096, num_classes=4, device="cuda"): model.eval() prob_sum = np.zeros((len(xyz), num_classes), dtype=np.float32) hit = np.zeros(len(xyz), dtype=np.int32) xs = np.arange(xyz[:, 0].min(), xyz[:, 0].max(), stride) ys = np.arange(xyz[:, 1].min(), xyz[:, 1].max(), stride) for bx in xs: for by in ys: m = ((xyz[:, 0] >= bx) & (xyz[:, 0] < bx + block) & (xyz[:, 1] >= by) & (xyz[:, 1] < by + block)) if m.sum() < 64: continue sel = np.nonzero(m)[0] if len(sel) > num_point: sel = np.random.choice(sel, num_point, replace=False) pts = xyz[sel] center = pts.mean(axis=0, keepdims=True) p = pts - center p = p / max(np.max(np.linalg.norm(p, axis=1)), 1e-6) t = torch.from_numpy(p.astype(np.float32)).unsqueeze(0).to(device) with torch.no_grad(): logits = model(t.transpose(1, 2)) prob = torch.softmax(logits, dim=1)[0].cpu().numpy().T # (n, C) # 只把参与预测的点累加回原始位置 prob_sum[np.nonzero(m)[0][:len(sel)] if len(sel) == m.sum() else sel] += prob hit[np.nonzero(m)[0][:len(sel)] if len(sel) == m.sum() else sel] += 1 hit[hit == 0] = 1 return (prob_sum / hit[:, None]).argmax(axis=1).astype(np.uint8)

上面这段里有一处需要按实际降采样逻辑对齐的写法：降采样后sel和块内全部点不再一一对应，累加时要严格按sel索引回原始数组，不能按m的顺序。stride取block/2是重叠推理的关键，重叠率越高边界越稳，代价是推理时间线性增长。hit记录每个点被预测了几次，没被任何块覆盖的点（通常是研究区边缘）会保持 0，除之前要兜一下底。

拿到逐点类别后，写回这一步反而简单，但两个检查不能省：

import numpy as np import laspy inv_lut = np.zeros(256, dtype=np.uint8) for src, dst in {2: 0, 6: 1, 3: 2, 4: 2, 5: 2, 9: 3}.items(): inv_lut[dst] = src pred = predict_tile(model, xyz) # 网络类别 0..3 las_code = inv_lut[pred] # 还原成 LAS 分类码 assert set(np.unique(las_code)).issubset({2, 6, 3, 4, 5, 9}), "出现未定义分类码" out = laspy.read("tile_01.las") assert len(out.points) == len(las_code), "点数必须和原文件一致" out.classification = las_code.astype(np.uint8) out.write("tile_01_pred.las")

断言是给自己看的。第一处检查反查后有没有冒出没定义过的码，第二处检查点数有没有对上——predict_tile如果中途丢过点，这里会立刻报出来，而不是等到打开查看器才发现少了一片。写回时保持原 header 不动，scale/offset 沿用原文件，别新建 header，否则坐标会漂。

一个能省很多时间的自查习惯：挑一块既有地面又有建筑的区域，把真值和预测并排渲染成两幅图对比，重点看建筑边缘和植被交界。mIoU 涨两个点可能只是地面拟合更好了，边缘才是真正决定可用性的地方。

从那以后我每次做 LAS 标签训练，都先把类别直方图打出来，把CLASS_MAP和inv_lut一起存进实验目录，跑完第一个 epoch 就抽一块可视化核对，宁可 early 阶段多花二十分钟，也不要在写回阶段对着一整片单色点云找问题。希望帮到你。

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