一种基于UAV-LiDAR点云的多尺度单木分割

A Multi-Scale Method for 3D Individual Tree Extraction Using UAV-LiDAR

摘要: 小型UAV-LiDAR (Unmanned Aerial Vehicle-Light Detecting and Ranging)可以快速获取林区三维数据，已经被广泛运用于林业调查。并且通过单木分割可以在激光点云中快速获取单株植被，进而获取单木几何参数。本文利用UAV-LiDAR进行数据采集，针对现有单木分割算法邻域尺度不易确定、复杂场景下单一尺度无法满足等问题设计由粗到精多尺度的单木分割策略，避免对邻域大小等参数的依赖。此外以往区域合并方法主要基于分割点云的水平投影的曲率、大小等，模糊了高程方向信息，由此本文基于点云三维分布特征设计合并准则。实验表明，本文方法在更有效的分割树冠的同时极大的抑制过分割的产生，测区内F得分达到了89%，较原方法提高7%左右。

Abstract: The UAV-LiDAR (Unmanned Aerial Vehicle-Light Detecting and Ranging) can quickly produce 3D point clouds in forest industry, and has been widely used for forestry characterization to facilitate forest ecological and management studies. The geometric properties of individual tree can be easily got from the point clouds, which is handled by the method of individual tree extraction. In this paper, the original point cloud is acquired by low-cost UAV-LiDAR. Because of the problems of the existing method, such as the difficulty in determining the scale of neighborhood and the situation of complex scenes, which is hardly processed by the method with fixed scale, we design a coarse-to-fine multi-scale method. In addition, the methods of previous regional merging are based on the curvature and size of the horizontal projection of segmented point clouds, which blur the information in elevation direction. Therefore, we design a new one, based on the 3D distribution characteristics of point clouds. The experiments show that the framework identifies individual trees more precisely, and greatly inhibits the generation of over-segmentation. The result F-score reaches 89%, which is about 7% more than the original method.

1. 引言

森林与我们的生态环境息息相关，我国定期进行森林资源调查，实施林业保护和扩增等政策。而目前调查工作人工参与成分较高，自动化程度低，严重影响工作效率。急需满足精度要求适应各种复杂场景的遥感手段进行原始信息的自动采集与获取。其中单木结构参数信息，例如单木位置、树高、冠幅等，对于树木生长研究、森林资源调查等具有重要意义。相比于摄影测量的方式，激光雷达可以直接的获取地物的表面点位三维信息，进而更准确的获取单木结构参数，而小型UAV-LiDAR (Unmanned Aerial Vehicle-Light Detecting and Ranging)更是为数据的获取提供了便利。由此激光雷达技术已逐渐成为森林资源调查重要的遥感手段，而在无序、离散的点云数据中准确分割单木点云成为其关键问题，尤其在面对复杂场景，树木种类不一、形状大小各异，如何提高分割精度也一直是该问题的难点。

国内外在机载点云单木分割上已经做了很多研究，从方法上主要分为两大类：1、基于冠层高度模型(Canopy Height Model, CHM)的单木分割；2、基于点云的单木分割。冠层高度模型(CHM)一般通过地物表面模型(Digital Surface Model, DSM)与地面高程模型(Digital Elevation Model, DEM)做差，并用局部插值生成的栅格图像，其基于CHM的单木分割实质上采用的图像处理的技术。基于点云的单木分割，充分保留高程上的信息，更加直接与精确，但分割准则更加复杂。

由于影像处理的算法更为丰富，目前基于CHM的单木分割一直是工程应用的热点。Hyppa等最早采用局部最大探测CHM钟树冠顶点，再运用区域增长寻找边界 [1]，但局部最大窗口难以控制。为此，Papescu等对对树高与冠幅大小建立关系确定冠层大小，进而确定窗口范围 [2]。Chen等也用树高与冠幅之间的统计关系确定搜索局部最大值的窗口大小，再以此为注水点进行分水岭分割 [3]。但不同树种树高与冠幅关系并不单一，难以确定准确函数关系。Mongus等运用局部拟合曲面来探测树冠顶点，并结合树干探测保证更准确的树冠点，以此作为标记进行分水岭分割 [4]，但在复杂场景下，尤其是密集区域冠层对树干存在遮挡，机载激光雷达无法获取树干点。总体而言，CHM生成的插值过程带有一定的误差和不确定性，并且CHM只保留了表面高程信息，存在对低矮植被的覆盖，导致一些小树木无法准确提取。

基于点云的单木分割，直接对原始点云进行处理，保留了更多高度信息。经典流程是先运用局部最大探测树冠顶点，再进行K-Means聚类 [5]，但没有从根本上解决基于CHM单木分割所存在的问题。为了减少高植被对低矮植被的影响，王濮等在局部最大探测树冠顶点基础上，运用图割算法分割点云，最后用全局最优代替局部最优实现对漏检单木的检测 [6]。郭庆华等总结了植被生长的规律，基于规则从顶至下聚类分割单木 [7]，但对于不同树种，其规则不尽适用。也有自底向上，先运用圆柱拟合探测树干，再运用图优化或图割的方式判断周围点云的归属，实现单木分割 [8] [9]，但其对点云质量要求较高，要包含底部完整的树干点云。Ferraz等分析植被树冠机载点云分布特点，其树冠中心不仅在高程方向上最高，在水平方向点密度也最大，由此最早运用空间特征设计均值漂移算法(MeanShift)进行单木分割，之后又采用不同高度不同带宽和高斯加权进行单木分割 [10] [11]，并应用于热带雨林复杂场景 [12]。Yan等人直接从数据本身出发，为每一个树冠点自适应确定带宽大小 [13]。Xiao等人将核函数形状改为Pollock核用参数具体确定核函数形状 [14]。代文霞等又将MeanShift改进增添光谱维度，进一步提高分割精度 [15]。最后Xiao等又对现有的MeanShift算法归纳总结并应用不同场景对比分析，探索不同情况下的适用性 [16]。相对而言，运用均值漂移进行单木分割，规则简单、参数控制较少，且适用场景多样，但对邻域带宽要求较高。

总之，现有单木分割算法基于植被生长特性和点云的分布特征，一般都是针对树种单一区域效果很好，而面对复杂场景，特别是树种冠层、高度、形状等差异较大的地方效果不佳。并且大多数算法对窗口、邻域等有一定的要求，不同环境，需要设置不同参数，单一阈值更是无法满足复杂场景的处理。随后一些自适应、多尺度方法应运而生，在一定程度上弥补单一尺度的不足，但其设计规则上同样基于一定参数，多数情况需要人工干预。针对以上问题，本文基于均值漂移算法，设计由粗到精多尺度分割策略，主要贡献有：1、减少对核函数带宽参数的依赖；2、弥补单一尺度无法实现复杂场景中多种树木的准确分割；3、基于三维点云分布特征设计过分割点云的合并原则。

2. 方法

2.1. 基于超体素分割的均值漂移

相对于其他算法，均值漂移算法形式更为简洁，没有过多的参数控制，仅依赖机载点云单木中心在水平方向密度最大、高程方向最高的特点。结合单木机载点云分布特点，均值漂移将三维空间拆分为两个方向：水平方向和垂直方向。核函数采用高斯核，依赖距离进行加权。在三维点集 $P_{3D}=\left\{X_i=\left(x_i,y_i,z_i\right),i=1,2,3,\cdots ,n\right\}$ 中，根据式(1)得任意点 $X_A=\left(x_A,y_A,z_A\right)$ 的均值漂移向量。其中s和r分别代表了水平方向和垂直方向， $g^s$ 和 $g^r$ 分别代表它们的核函数， $h^s$ 和 $h^r$ 为其带宽。 $X_A^s=\left(x_A,y_A\right)$， $X_A^r=\left(z_A\right)$， $X_i^s=\left(x_i,y_i\right)$， $X_i^r=\left(z_i\right)$， $X_i$ 由 $X_A$ 与核函数形状与带宽大小确定。其中水平核函数 $g^s$ 为式(2)，垂直核函数为式(3) (4)。

$m_{h^s{h}^r}\left(X_A\right)=\frac{{\displaystyle {\sum }_{i=1}^n{X}_i{g}^s\left({\Vert \frac{X_A^s-X_i^s}{h^s}\Vert }^2\right)g^r\left({\Vert \frac{X_A^r-X_i^r}{h^r}\Vert }^2\right)}}{{\displaystyle {\sum }_{i=1}^n{g}^s\left({\Vert \frac{X_A^s-X_i^s}{h^s}\Vert }^2\right)g^r\left({\Vert \frac{X_A^r-X_i^r}{h^r}\Vert }^2\right)}}-X_A$ (1)

$g^s\left({\Vert \frac{X_A^s-X_i^s}{h^s}\Vert }^2\right)=\{\begin{array}{l}{\text{e}}^{-\frac{1}{2}{\Vert \frac{X_A^s-X_i^s}{h^s}\Vert }^2}\Vert \frac{X_A^s-X_i^s}{h^s}\Vert \le 1\\ 0其他\end{array}$ (2)

$g^r\left({\Vert \frac{X_A^r-X_i^r}{h^r}\Vert }^2\right)=\{\begin{array}{l}1-{\Vert 1-dist\left(X_A^r,X_i^r\right)\Vert }^{2}mask\left(X_A^r,X_i^r\right)=1\\ 0其他\end{array}$ (3)

$dist\left(X_A^r,X_i^r\right)=\{\begin{array}{l}\min \left(\Vert \frac{\left(X_A^r-\frac{h^r}{4}\right)-X_i^r}{\frac{3h^r}{8}}\Vert ,\Vert \frac{\left(X_A^r+\frac{h^r}{2}\right)-X_i^r}{\frac{3h^r}{8}}\Vert \right)mask\left(X_A^r,X_i^r\right)=1\\ 0其他\end{array}$ (4)

均值漂移算法是一个迭代的过程，对每个种子点进行MeanShift向量计算，即该点的偏移均值，并将其移动到偏移后位置。更新点位后继续迭代进行，达到稳定位置停止(偏移向量模长小于ε = 0.1 m)，对不稳定的点设置最大的迭代次数。然后将相邻距离较小的点进行合并聚类，最终得到单木树冠点位和单木分割结果，其中最终迭代点代表了树冠顶点或树冠中心点，并称之为模点。

在初始种子点选取上，假如将原始点云全部参与迭代计算，其计算量过于庞大，并存在大量的冗余，而算法本身代表了区域性(核函数范围内)的点云分布规律，由此可以采用体素点来代替体素内所有点参与计算。而规则体素在一定程度上会将相邻的不同类点包含在一起，导致分割结果边界不准确。在此采用超体素 [17] 分割代替体素规则分割，获得边界明确，体素内点类一致的不规则超体素，并用超体素中心点作为初始种子点，提高算法的运算效率。

此外，原算法逐次对每一个种子点进行迭代计算，并最后相邻合并。考虑到在跌代偏移的过程中，种子点存在近邻情况，并且由于算法的区域性特点，它们此后的迭代偏移位置一定相近，并最终稳定到相近结果。为此，我们采用边迭代边合并的策略，对每一步迭代后的种子点进行相邻合并操作。

2.2. 由粗到精多尺度单木分割

在一般情况下，均值漂移核函数带宽 $h^s$ 设置为树冠半径， $h^r$ 设置为植被高度，能达到较优的分割结果。而对于区域内树种冠层、高度、形状等参差不齐，很难通过固定带宽大小对整个区域进行单木分割，往往在带宽大小偏小时过分割明显，偏大则相邻植被区分性变差，欠分割明显。由此设计由粗到精多尺度的分割策略，总体流程框架如图1所示。其中对每一次分割结果通过机器学习的方式进行分割状态(过分割、欠分割、正常)判别。对欠分割点云减小尺度再分割，对过分割点云进行判断与合并操作。通过迭代进行多层次的分割与合并处理，实现点云由粗到精多尺度分割。具体算法流程如下：

输入数据：原始点云。

参数：核函数初始核函数带宽 $h^s,h^r$ 、步长 ${\lambda }^s,{\lambda }^r$。

1) 原始点云预处理，将地面点投影并内插为DEM；

2) 对目标点云运用基于超体素分割的均值漂移算法(带宽 $h^s,h^r$ )进行单木分割；

3) 统计分割点云特征，并逐个分类。对于欠分割点云令 $h^s=h^s-{\lambda }^s$， $h^r=h^r-{\lambda }^r$ 返回步骤2)，对于过分割点云执行步骤4)；

4) 通过合并准则判断并寻找满足点云 $C_i$ 合并条件的点云 $C_k$，并将 $C_i$ 与 $C_k$ 合并，返回步骤3)；

5) 输出点云单木分割结果。

2.2.1. 点云预处理

预处理主要目的是改善原始点云的质量，在一定程度上排除噪点的干扰，并且将地面点与非地面点进行分离，提取目标点云。由于布料模拟模型(Filtering Method Based on Cloth Simulation, CSF) [18] 符合多种自然状态的地表情况，通过模拟布料自由下落状态并平铺向地面，将贴近布料的点设为地面点，其余为非地面点。由此选用CSF算法进行地面滤波，滤波结果如图2所示，黄色为地面点，绿色为非地面点，其中非地面点为主要研究区域。对于地面点，为了方便后续树高等特征的计算，将其投影并重采样获得该区域的DEM。对于非地面点，即目标点云，删除1 m邻域点数小于5的噪点，并剔除地面高度小于1 m的低矮地物。

Figure 1. The overall framework

图1. 总体框架

Figure 2. After filtering, the ground points and the non-ground points

图2. 滤波后地面点和非地面点

2.2.2. 点云状态分类

对于单一邻域带宽分割结果存在一定的欠分割和过分割现象，本文采用机器学习的方法对分割结果进行分类。其关键在于通过特征设计和特征组合等方式，设计可供区分不同状态的特征，已实现对过分割、正常、欠分割三种状态的准确判断。对于点云，应尽可能的充分应用其分布的几何特征，从不同层次，不同程度表达与描述。如图3所示，结合单木点云的分布特点和不同状态下的具体情况设计如表1所示9个特征，其中包含了树冠面积、高度、水平投影曲率、树冠层次分布、点云主方向等特点，从不同角度刻画点云特性，以便将其区分，并采用随机森林分类器进行分类，通过预训练该分类器得到点云分割状态的预测模型。

Figure 3. Diagram of feature calculation of crown clusters

图3. 单木点云特征计算示意图

Table 1. Features extracted from crown clusters

表1. 单木点云特征提取

2.2.3. 合并原则

以往合并主要基于二维平面，借助影像处理的合并准则，只考虑欠分割目标在水平投影的形状、曲率、邻接关系等，模糊z方向的信息。例如由于高植被对周边低矮植被的部分覆盖，影响其在水平方向的判断。由此本文基于三维点云特征对欠分割点云是否需要合并和与哪个点云合并进行判断。合并判别式如式(5)。

$F_{ij}=L_{ij}+p_{{\sigma }_d}\left(d\right)+p_{{\sigma }_{overlap}}\left(1-overlap\right)$ (5)

$L_{ij}=\left|1-l_{ij}\right|l_{ij}=\{\begin{array}{l}0过分割\\ 1正常\\ 2欠分割\end{array}$ (6)

$p_{\sigma }\left(x\right)=\{\begin{array}{l}{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}\left|x\right|\le limit\\ 0其他\end{array}$ (7)

其中 $F_{ij}$ 表示点云 $C_i$ 与点云 $C_j$ 合并系数， $l_{ij}$ 代表 $C_i$ 与 $C_j$ 合并后运用分类器模型对其标签预测结果， $p\left(x\right)$ 采用高斯函数确定依赖变量的可能性大小。 $p\left(d\right)$ 通过距离判断两点云合并可能性，d为点云模点之间的水平距离，r为点云 $C_j$ 树冠半径，可通过分割点云的最小外包框的长轴确定，并且当点云 $C_i$ 在点云 $C_j$ 边界位置时合并可能性为50%，由此可得：

${\sigma }_d=0.85r$ (8)

由于点云重叠区域几乎为一个平面，而平面占立体的大小几乎为零，由此应先定义如何计算重叠率才能更好的表达点云的重叠关系。我们将点云的临界区域占点云 $C_i$ 水平投影的比值定位overlap，比值越大代表将 $C_i$ 并入 $C_j$ 的可能性越大。假使仅通过水平投影的重叠度判断，那么在高植被与低矮植被近邻情况下，存在三维空间不重叠而二维水平投影重叠的现象。为了更好的保留点云的三维空间分布又方便表达两片点云的重叠程度，设计如下计算方式：

1) 寻找点云 $C_i$ 中与点云 $C_j$ 邻近的点(点到 $C_j$ 的距离小于 $d_t$ ) $P_i$，并且 $P_i\subseteq C_i$ ；

2) 计算 $P_i$ 水平投影面积占 $C_i$ 的比例即为overlap。

由图4(a)计算两个临近的树冠最大重叠度为：

$overla{p}_{down}=\frac{S_{P_i}}{S_{C_i}}\approx 0.14$ (9)

由图4(b)计算需要合并的最小重叠度为：

$overla{p}_{up}=\frac{S_{P_i}}{S_{C_i}}\approx \frac{4d_t}{\pi r}$ (10)

$overla{p}_{mid}=\frac{overla{p}_{up}+overla{p}_{down}}{2}$ (11)

由此假设当重叠度为 $overla{p}_{mid}$ 时合并可能性为50%，得

${\sigma }_{overlap}\approx \frac{1-overla{p}_{mid}}{1.18}$ (12)

(a) (b)

Figure 4. Adjacent canopy clusters distribution, (a) Maximum overlap in case where merging isn’t required. (b) Minimum overlap in case where merging is required

图4. 相邻点云分布，(a) 代表不合并两临近树冠最大重叠情况；(b) 代表需合并两点云重叠最小情况

由以上关系式，确定合并判别式的各项参数，而各项参数只由数据本身影响，不受人为控制。最后通过合并原则判断的具体流程如下：

1) 找出与过分割点云 $C_i$ 距离(模点距离)小于5 m的所有点云 $C_j$， $j=1,2,3,\cdots ,n$ ；

2) 分别计算式(5) $F_{ij}\left(C_i,C_j\right)$， $j=1,2,3,\cdots ,n$ ；

3) 选出最大 $F_{ik}=\max \left(F_{i1},F_{i2},\cdots ,F_{in}\right)$ ；

4) 判断 $F_{ik}>1$, 满足则点云 $C_i$ 与 $C_k$ 合并。

3. 实验分析

3.1. 实验数据

我们采用轻小型低成本无人机激光扫描系统 [19] 进行数据采集与实验，选区为湖北省武汉市马鞍山森林公园某区域，研究区内地势相对平坦，海拔19 m~32 m，坡度0˚~15˚，主要树种为阔叶林，采集时间为2018年夏末，平均点密度24 pts/m²。采集原始点云如图5所示，其中点云经高程赋色显示。

Figure 5. The UAV laser scanning system collects raw data

图5. 无人机机载激光扫描系统采集原始数据

测区内有边界明显稀疏区域，也有连接紧密的密集区，主要有低矮灌木、常绿乔木、落叶乔木等，依据树种更是有5种以上植被类型，并且形状、大小、高度差异较大。测区内人工与自然林交杂，人工林占主体。

3.2. 实验结果

3.2.1. 点云状态分类

如图6表现了各个特征在欠分割、正常和过分割状态下的统计信息。为了方便比较，统一每类样本总量。从实际样本的统计特性来看，正常分割结果树冠面积大部分分布在9 m²左右，而相对欠分割样本更大，过分割偏小。C-1与C-2表现了树冠的曲率特征，其中正常分割结果的值更接近于1。Ratio-1表现了主干的承受能力，正常分割结果更为适中，而过分割结果树冠面积偏小，其计算结果偏小，欠分割与其相反。其他四个偏置项，正常样本的分布规律明显，并且数值偏小(Bias-3代表了PCA主方向与z方向夹角的余弦，由此偏大)。以上几个特征从不同方面刻画点云状态，增强其可分性。

随机森林设计单棵树所能达到的最大深度不超过10，树节点持续分裂的最小样本数量为10，树的每个节点随机选择变量的数量控制在4个左右。总计标住801样本，训练样本434，测试样本367。对样本进行训练并测试，其训练精度97%，测试精度可达85.8%。

3.2.2. 单木分割结果

通过反复实验选取均值漂移算法较优得带宽邻域进行比较分析。图7中显示了单尺度下各种过分割、欠分割和混杂存在的情况，图7(a)中由于单株多峰或带宽较小等情况使得单木被分割为多块，通过合并准则判断，其距离相近，重叠范围较大，很容易将其合并。图7(b)显示了由于树木连接紧密，带宽选取较大而导致主树干(较高)与周边树木无法区分，在多尺度分割过程中，总会通过较小带宽将其区分。图7(c)显示了在复杂场景数据质量不佳的情况下，欠分割与过分割往往是交错出现，通过多次迭代，即多次的分割与合并操作，以实现该情况下较优分割结果。

Figure 6. Feature statistical chart

图6. 特征统计图

Figure 7. Comparison between original method and multi-scale segmentation

图7. 原算法与多尺度分割前后对比

此外我们选取区域性结果进行对比分析，如图8显示了区域单木分割结果，并与原均值漂移算法、基于距离聚类方法 [7] (基于点云分割)进行对比分析。分类精度采用以下3个指标进行衡量 [20] [21]。

$r=\frac{N_{match}}{N_{match}+N_{om}}$ (13)

$p=\frac{N_{match}}{N_{match}+N_{com}}$ (14)

$F=\frac{2\left(r\times p\right)}{r+p}$ (15)

式中， $N_{match}$ 、 $N_{om}$ 、 $N_{com}$ 分别表示正确分割、漏分割和过度分割的数目；r代表冠幅探测率，p代表探测出树冠的准确率，并通过F评分整体评价。

Figure 8. The results of individual tree extraction in different areas

图8. 不同区域单木分割结果

表2统计了不同区域不同算法的分割精度。对比原始MeanShift算法，核函数带宽过大漏检过多，欠分割显著，带宽减小错检增多，即过分割现象严重，采用本文方法兼顾了不同带宽大小对不同大小树冠的分割能力，又一定程度上抑制了过分割现象，综合F评分提高7%左右。相比与另一种基于点云的单木分割，也有更准确的数木探测能力。

Table 2. Accuracy evaluation and methods comparative analysis

表2. 实验结果精度评定与对比

在plot1与plot2内过分割却并不明显，而在全区域内出现较多的过分割现象。由此在原始数据中查找过分割严重区域情况如图9所示。通过与原场景对比发现，一般在边界并且植被高大繁密的情况下，由于无人机搭载激光雷达测距有限，飞行高度较低，造成这些地方数据采集不完整，如图9区域1仅为一棵树，由于侧面和内部数据不完整导致被分割成若干段，且每段之间的邻接程度有限，无法满足合并要求，由此在单木点云数据区域性缺失较大时，很难将其分割为完整单木。

Figure 9. Over-segmented region

图9. 过分割明显区域

4. 结论

以往点云单木分割算法对窗口、邻域等十分敏感，只能处理植被密度，树种种类相对单一的场景，面对树种较多且分布不均等情况，单一尺度很难满足分割要求。随后一些自适应的方法也是通过一定情况下的统计信息，建立树高与冠幅的关系，很难普遍适用于纷繁复杂的现实情况，此外它们往往还需要一些额外的经验参数控制。本文基于以上问题，设计多尺度的分割策略，由粗到精的对整个测区进行单木分割，更全面的对场景内不同形状、不同大小单木进行提取。此外在小尺度下往往产生严重的过分割现象，由此设计一种基于三维点云分布特征的合并准则，对过分割点云进行判断与合并操作，极大抑制过分割的产生。最后本文运用低成本无人机机载雷达扫描系统进行实地数据采集和算法验证，实验表明，本文方法综合F评分达89%，较原算法提高7%左右，极大减少了过分割和欠分割等现象的影响，较其他基于点云的单木分割算法精度更高。但是对于数据采集不完整的区域，点云数据较为破碎，分割点云之间的重叠度较小，很难满足合并要求。由此在数据采集过程中，应尽量覆盖完整的测区范围，避免边界区域因遮挡导致数据缺失。