前言

在处理森林场景的激光雷达点云数据时，我过去往往不进行去噪处理。因为在我看来，这一步似乎没有太大必要——真的如此吗？
从流程上看，点云去噪通常被视为点云预处理的标准步骤之一。但在实际操作中，是否有必要、该如何去噪与降采样，其实取决于点云的采集方式、质量以及我们最终的分析目的。

什么是点云中的“噪点”？

点云数据中的噪点，并没有一个明确的、统一的定义。但从实用角度看：

凡是影响你后续数据处理或建模精度的异常点，都可以视为噪点。

举个常见的例子：在机载激光雷达（ALS）扫描过程中，若下方恰好飞过一只鸟，它会被误记录为一个高于林冠的异常点。这类噪点通常具有明显的空间孤立性，可通过以下方式剔除：

• 高度阈值法，设定某个高度阈值进行去除；

• 裁剪法，在软件上直接手动裁剪进行去除；

• 孤立点搜索法，即在每个点的邻域（如3米半径）内统计邻点数量，若低于某一阈值（如2个），则将其视为孤立点。

这类“孤立点噪声”在森林点云中相对较少，而且易于识别。因此，我的处理方式往往是直接裁剪，速度快，效果也够用。

去噪实践：孤立点搜索法

我曾尝试使用孤立点搜索法对样地点云数据进行去噪，思路很简单：判断每个点在3米内是否有至少2个邻点，若不足则视为噪点删除。但由于此方法需逐点搜索邻域，计算量非常大，我试过了几种算法，最终Open3D 的半径滤波比较快，文末附代码。

不那么明显的“噪声”

还有一类“噪声”不那么显眼，却常常被忽视：由于设备本身的扫描精度限制，激光束在扫描同一表面时可能出现微小的空间漂移，导致同一位置上记录了多个偏移点，如下图所示。这在可视化中常表现为“重影”或“毛刺”，不仅影响美观，更重要的是显著增加了点的数量，影响处理效率。

但是这种情况无法使用现有的去噪算法进行去除，需要使用降采样的算法，去除冗余数据。但是降采样我们需要注意的是保持单木点云特征，因为我们需要提取胸肌、树高等关键参数，保证参数提取的准确性是很重要的。

点云降采样的必要性

目前我们常用的地面手持扫描样地，点密度大大，但一个样地数据在3 GB 左右，机载点云比较稀疏，扫描完一个林分差不多也有几GB。未来，随着无人机林下激光扫描、林下机器人等新型采集技术的发展，点云数据的分辨率和覆盖范围都将大幅提升，数据量将显著增长。

为保证在不牺牲胸径、树高等精度的前提下提升处理效率，合理降采样将变得越来越重要。所以我就在想，怎样保持树木的点云特征，还能删除冗余点呢，我查了下资料。点云降采样方法很多，包括：

• 随机抽样；

• 均匀抽样；

• 体素网格降采样；

• 法向空间均匀抽样；

• 基于微分几何的曲率加权抽样。

其中，法向空间降采样和基于微分几何的曲率加权抽样据介绍是能够保持点云特征的。如法向空间降采样，通过在法向量空间中进行均匀抽样，使得所保留点之间的法线方向尽可能分散。结果表现为：几何特征变化大的区域保留点较多，变化小的区域则点较少，是一种特征保持效果较好的方法。

降采样实践：

• 法向空间均匀抽样，其是将x/y/z的法向量范围分区间，比如每个分六个区间就有6的3次方组合方式，得到216个箱子，计算每个点的法向量，将每个点分到每个箱子内，在去取每个箱子内的代表点。但是这种特征保持的算法并不适用于森林场景，一颗树木的法向量几乎代表整个样地所有点的法向量，这会导致点云很破碎。文末附代码

• 基于微分几何的曲率加权抽样。我在网上找了一个代码，运行速度慢，而且单木主干的部分点云还被过采样了，我没找到更好的算法了。

• 随机抽样：最不起眼的随机抽样，到是运行速度快，且结构特征保持不错。如果你觉得不放心，还可以对点云分箱在降采样

结论与建议

在森林点云处理任务中，我的实践经验如下：去噪手动裁剪，降采样从目前的数据量来说可以不做，或者随机降采样。

• 去噪：存在的明显孤立点（如飞鸟、探头干扰等），直接通过目视裁剪方式清理，简单高效；

• 降采样：若数据量较大，有必要降采样，随机降采样反而是性价比最高的选择，运行快，效果也能接受。适用于森林场景的特征保持的点云降采样算法还得研究。
  写这段时我在想点云数据很多不需要的字段或许可以清空，以降低数据量的大小
  点云去噪代码

import laspy
import numpy as np
import open3d as o3d

#1. 读取点云数据
las = laspy.read(r"C:\Users\YSYmc\Desktop\近期内容\lidar_data\test.las")
points = np.vstack((las.x, las.y, las.z)).T
## 转为 Open3D 点云
pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(points)

#2. 使用 Open3D 的半径滤波，高效过滤（底层为 C++ 实现）
radius = 3           # 搜索半径（米）
min_neighbors = 2      # 半径内至少邻居点数

filtered_pcd, ind = pcd.remove_radius_outlier(nb_points=min_neighbors, radius=radius)
filtered_indices = np.array(ind)

#3.保留原 LAS 属性字段
new_las = laspy.create(point_format=las.header.point_format, file_version=las.header.version)
new_las.header = las.header

for dim in las.point_format.dimension_names:
    if hasattr(las, dim):
        setattr(new_las, dim, getattr(las, dim)[filtered_indices])

#保存
new_las.write(r"C:\Users\YSYmc\Desktop\近期内容\lidar_data\test_deno.las")

法向空间均匀抽样代码

import laspy
import numpy as np
import open3d as o3d
import time
'''
normal_space_sampling可以这么理解：
x,y,z都分成几个区间，比如这里是6个区间，然后再组合，这有6的3次方组合方式，就得到216箱子。
根据每个点计算的法向量，将各点分到对应的箱子里
最后在每个箱子内选择一个代表点保留，如果不足我们最终要的保留点数10000，则随机补充.
该代码借鉴：https://blog.csdn.net/limeigui/article/details/144891669
'''

def normal_space_sampling(pcd, num_bins=6, num_samples=10000):
    #pcd是点云对象，num_bins是x/y/z要分的箱数，num_samples表示最终保留的点数
    # 估算法向量
    pcd.estimate_normals(search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=10, max_nn=30))
    #radius表示半径为10m，最多点数为30
    normals = np.asarray(pcd.normals)#此为得到各点的x/y/z的向量，[N, 3]

    # 法线方向分箱
    bins = np.linspace(-1, 1, num_bins)
    normal_bins = np.digitize(normals, bins)
    unique_bins = np.unique(normal_bins, axis=0)

    #取每个箱子的内的代表点
    sampled_indices = []
    for b in unique_bins:
        indices = np.all(normal_bins == b, axis=1)
        bin_points = np.where(indices)[0]
        if bin_points.size > 0:
            sampled_indices.append(np.random.choice(bin_points))

    # 不足则随机补充点
    while len(sampled_indices) < num_samples:
        sampled_indices.append(np.random.randint(0, len(pcd.points)))

    return np.array(sampled_indices)


#读取点云数据
input_file = r"C:\Users\YSYmc\Desktop\近期内容\lidar_data\test.las"
output_file = r"C:\Users\YSYmc\Desktop\近期内容\lidar_data\test_sample.las"

las = laspy.read(input_file)
points = np.vstack((las.x, las.y, las.z)).T

pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(points)

#2.法线空间采样
start = time.time()
sampled_indices = normal_space_sampling(pcd, num_bins=6, num_samples=10000)
end = time.time()
print(f"耗时 {end - start:.2f} 秒")

#3. 构建新 las 对象，保留原字段
new_las = laspy.create(point_format=las.header.point_format, file_version=las.header.version)
new_las.header = las.header

for dim in las.point_format.dimension_names:
    if hasattr(las, dim):
        setattr(new_las, dim, getattr(las, dim)[sampled_indices])

# 保存数据
new_las.write(output_file)