基于局部区间极大值的激光扫描点云树木自动提取方法与流程

本发明涉及一种基于局部区间极大值的激光扫描点云树木自动提取方法。

背景技术：

利用激光扫描技术进行森林调查，可以安全高效地获得扫描样地内的完整三维点云数据。该点云数据反映了扫描样地内植被，树木，地形等对象完整的形态结构，为树木参数调查提供了可靠的数据支持。在三维点云数据中进行单株木提取是进行树木参数调查和森林资源调查的基础和前提，其结果直接决定了后续调查及计算的正确性。然而点云林区场景的数据量极大，且场景的可视性受森林冠层的干扰严重，不利于从整体上直接进行观测。此外，采用交互式手动割取单株木的方法工作量大，需要在三维空间上进行切割，操作难度大，分割效果粗糙，易带入噪声。而基于点云数据的树木自动提取技术将各株树木所属数据点单独划分出来，可以为自动森林资源调查提供有效的解决途径。

在三维点云数据中进行树木自动检测和提取，是将属于同一株树木(包括树干和树冠)的数据点集从原始数据中提取出来。一方面，提取的树木需要保留完整的结构，即提取完整的树干以及树冠数据，从而满足单株胸径和树高等参数的测量需求。另一方面，提取的树木数据应该尽可能排除杂草以及地面的数据点，为后续计算减少负担。

林区通常环境复杂，不仅地形坡度大，且局部地形起伏不平，对自动提取方法带来考验。此外，各林区树种存在差异，而不同的树种形态、高度、大小各异，变化多样，要适应不同树种的有效提取，对提取方法的鲁棒性产生极大的要求。而点云数据本身受扫描仪位置影响，其数据密度分布不均、数据遮挡缺失等问题也容易对提取结果的稳定性造成负面影响。从海量点云中快速、自动地提取树木虽然难度大，要求高，但该技术具有重要的经济和实用价值，一直是国内外的研究热点。

目前，已有的面向三维激光点云的单木自动提取方法主要分为面向机载点云提取方法和面向地面点云提取方法两类：面向机载激光点云数据的提取方法主要考虑的是树冠的提取，这是由于机载点云俯视扫描的特点，使树木信息更多的表现在树冠结构上，常用的方法为基于局部邻域主成分分析(PCA)的方法，基于分水岭的方法，基于几何分割规则的方法等，其缺陷在于：利用树冠进行处理难以克服树冠顶端被冠层遮挡的情况，也难以保证树干的完整提取。而面向地面激光的单木提取方法主要考虑了树干和树冠的完整保留，对于林区场景，主要方法是在固定高度(如1.3米)采用截面拟合圆判定的方法，对于城区场景，则通常采用去除地面后，欧式聚类或超体素分割的方法，然后通过计算PCA维数，RGB颜色，FPFH等特征，从而识别和提取树木目标，其缺陷主要在于：

1、采用固定高度截面拟合圆判定的方法在数据部分遮挡缺失的情况下容易遇到困难，森林植被茂密，数据获取时远处的树木因为前面物体的遮挡会导致部分树干缺失，从而导致无法拟合圆形，最终丢失树木；

2、采用先去除地面，再利用欧式聚类构建目标对象的方法，首先地面滤除在起伏较大的林区场景中容易留下地面点噪声，此外，欧式聚类难以克服重叠目标的分割，林区中树冠容易交错，聚类后易形成大片树木相连的情况，森林层次复杂，杂草的生长和特异的树木形态干扰形态特征的计算，而颜色信息也会受到林地郁闭度的影响(光照)而不稳定。

技术实现要素：

本发明的目的在于克服上述现有技术的不足，提供一种基于局部区间极大值的激光扫描点云树木自动提取方法。

为实现上述目的，本发明采用以下技术方案：

基于局部区间极大值的激光扫描点云树木自动提取方法，包括以下步骤：

(1)对获得的三维点云数据进行网格划分，建立非空网格索引，将空网格排除；

(2)遍历每一个非空网格，计算非空网格中高程区间的累积值；

(3)获取树木的中心点集:

a、设置能量粗分隔T_E对网格进行阈值粗分割，保留累积值E≥T_E的网格，排除累积值E＜T_E的网格并将其置为空网格；

b、计算阈值粗分割后保留下来的每个网格的水平中心c＝[c_x,c_y]，所有保留网格水平中心c＝[c_x,c_y]构成了嫌疑树干中心点集C；

c、为嫌疑树干中心点集C构建Kd-tree索引结构，设置树木半径r，遍历嫌疑树干中心点集C中的点c_i，以该点为查询点，搜寻其邻域半径内的邻域中心点c_j(c_j∈C)，若为空集(除查询点本身)则保留，若非空则依次查找各个邻域中心点c_j，比较该邻域中心点c_j与查询点c_i的高程区间累积值，若查询点的累积值为该邻域内所有中心点中的最大值则标记为树木中心点，并将其他点标记为非树木中心点，若查询点不是最大值，则标记为非树木中心点，在遍历完所有点后将该点排除，最终得到树木中心点集为

(4)树木的提取

将原始的三维点云数据集P投影于XOY平面上，从而获得投影点并为该点集建立Kd-tree索引，然后遍历每一个树木中心点c_tr(c_tr∈C_tr)，每遍历一个树木中心点，以该中心点作为查询点，搜索投影点集中距离小于r的数据点并记录该索引值i，构成单株木的点云索引序列，最后在原始点云中将属于该索引序列的点集提取出来作为单株树木点云，完成提取树木操作。

进一步地，步骤(1)具体包括以下步骤：

a、对获得的整个三维点云数据集P，按照XOY平面进行网格划分，其中网格为正方形网格，设置网格边长为l，将每个点p_i＝[x_i,y_i,z_i]，p_i∈P，其中i＝1,2,...,s，s为点的个数，按每个点的[p_x,p_y]值分别将点划入对应所属的网格中，并建立相应的索引关系；

b、记录非空网格的索引值，将空网格排除。

进一步地，步骤(2)具体包括：

a、计算每一个非空网格的区间的最大值及最小值，建立高程区间[z_min,z_max]，利用高程区间分辨率z_res对高程区间进行分段统计，分段操作为

bin(n)＝[z_min+(n-1)z_res,z_min+nz_res)

式中bin(n)为高程区间，该高程区间包含区间内所有的三维数据点，n为高程区间的序号；

b、通过高程分段统计后，计算每个分段区间是否有数据占据该空间，因此对于高程区间的判断有二值函数

其中#表示统计当前高程区间内的点个数，n为当前高程区间的序号，N为个数阈值，

则高程区间的累积值定义为：

其中m为高程区间数量,从而每一个网格获得对应的累积值E，并记录。

进一步地，在步骤(2)与步骤(3)之间还有对网格进行重叠平移划分的步骤，具体为：采用部分重叠网格划分法对已划分的网格顶点V＝[x_v,y_v]重新计算新的顶点位置V′＝[x_v',y_v']，并将每个原始点p_i＝[x_i,y_i,z_i]重新划入对应的网格，再次建立新的网格与原始点的索引关系，在新的网格上重复进行非空网格索引建立和非空网格高程区间累积计算操作，获取新划分网格坐标下每个网格对应的累积值E′。

进一步地，所述部分重叠网格划分法的重叠部分取决于设置的重叠比例n_offset，则新的网格顶点位置V′＝[x_v',y_v']＝[x_v+l/n_offset,y_v+l/n_offset]。

进一步地，在步骤(3)中利用能量粗分隔T_E分别对两次网格划分后的结果进行分割，获得两次的网格中心点集C₁及C₂，则所述嫌疑树干中心点集C＝C₁∪C₂。

采用上述技术方案后，本发明与背景技术相比，具有如下优点：

1、本发明通过计算高程区间累积值来描述各个网格内的点集空间分布情况，利用树干在垂直空间相对于地面和杂草对象显著性强的特点，通过寻找高程累积值较大的网格来获取潜在的树干位置。该方法考虑的是空间点集累积的结果，对密度和遮挡不敏感，对树干任意局部缺失具备抵抗力，克服了点云数据密度远近(距离扫描仪)分布不均的缺点，使得结果受采集设备影响小，更具稳定性；

2、本发明通过对网格的半重叠平移，进行多次区间累积计算，从而精确获得局部潜在的高程空间累积极值，避免了树木定位不精确的提取结果，也间接避免了在提取计算过程中树边杂草对树干定位的干扰；

3、本发明采用了非极大值抑制的方法对中心点集进行判定，从而精确定位潜在树干的位置，可以更好地区别地面与杂草的干扰，避免了形态特征在不同树种间计算的不稳定，对多种树木的自动提取都具有有效性，采用将潜在树干附近的中心点(即树冠位置)消除的方法，并在最后以树干中心点为基点，将范围内的点集一并提取，达到整树合并的效果，避免了同一株木的反复提取，使结果具有可靠性；

4、本发明采用网格块为单元进行计算，计算速度快，可以快速，自动提取林区场景三维点云数据中的单株树木点集，提取的单株树木包含树干和树冠结构，更适用于大规模点云场景，具有实际意义和应用价值。

附图说明

图1为本发明算法的流程示意图；

图2a为本发明实施例原始点云数据；

图2b为本发明树木提取结果示意图；

图3为实施例中不同树种的单株木提取结果示意图，由左往右依次为马尾松、杉木、阔叶林。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

实施例

如图1所示，基于局部区间极大值的激光扫描点云树木自动提取方法包括以下步骤：

(1)对获得的三维点云数据进行网格划分，建立非空网格索引，具体为：

a、输入整个三维点云数据集P，计算点集P在三个正交坐标轴方向上的最大最小值，然后设置正方形网格，每个网格的边长设置为l(本实施例中l＝0.2米)，根据计算获得的最大最小值在XOY平面对原始点云进行网格划分，将每个点p_i＝[x_i,y_i,z_i](p_i∈P，i＝1,2,...,s，s为点的个数)按每个点的[p_x,p_y]值将该点划入对应所属的网格中，并建立相应的索引关系；

b、在遍历每个点的过程中，记录每个网格内存在点的个数，建立非空网格索引，并在后续的处理中将空网格排除；

(2)遍历每一个非空网格，计算非空网格中非空高程区间的累积值E，具体为：

遍历每一个非空网格,每遍历一个非空网格时，先计算该区间的最大最小高程值，建立高程区间[z_min,z_max]，设置高程区间分辨率z_res(如z_res＝0.1米)，则利用高程区间的范围，进行高程分段统计，分段操作为

bin(n)＝[z_min+(n-1)z_res,z_min+nz_res)

通过高程分段统计后，计算每个分段区间是否有数据占据该空间，因此对于高程区间的判断有二值函数f(bin(n))有

式中，bin(n)为高程区间，该高程区间包含区间内所有的三维数据点，#表示统计当前高程区间内的点个数，n为当前高程区间的序号，N为个数阈值。遍历每一个区间段进行统计，然后定义该网格的高程能量为该网格内非空高程区间的累积值E，因此该能量累积函数定义为：

式中，m为高程区间数量，即m＝(z_max-z_min)/z_res,从而每一个网格获得对应的累积值E，并记录。

完成步骤(2)后进行网格重叠平移划分，具体如下：

采用部分重叠网格划分法对已划分的网格顶点V＝[x_v,y_v]重新计算新的顶点位置V′＝[x_v',y_v']，所述部分重叠网格划分法的重叠部分取决于设置的重叠比例n_offset，则新的网格顶点位置V′＝[x_v',y_v']＝[x_v+l/n_offset,y_v+l/n_offset]。本实施例中采用半重叠网格划分法，则对于已划分的网格顶点V＝[x_v,y_v]重新计算新的顶点位置V′＝[x_v+l/2,y_v+l/2]，并将每个原始点p_i＝[x_i,y_i,z_i]重新划入对应的网格，再次建立新的网格与原始点的索引关系。在新的网格上重复进行非空网格索引建立和非空网格高程区间累积计算操作，获取新划分网格坐标下，每个网格对应的累积值E′；

(3)以非极大值抑制法获取树木的中心点集:

首先设置能量粗分割阈值T_E，T_E为无量纲变量，例如利用T_E＝10。粗分割阈值T_E将两次网格划分后的结果(设为E，E′)进行分割，保留大于等于该阈值的网格(E≥T_E或E′≥T_E)，排除小于该阈值的网格，将网格置为空网格；

计算阈值粗分割后保留下来的每个网格的水平中心c＝[c_x,c_y]，就是计算该网格内点集的水平中心，所有保留的网格中心构成了中心点集C_n作为嫌疑树干中心点集，将两次网格的中心点集C₁,C₂进行合并操作，则最终的嫌疑树干中心点集即为C＝C₁∪C₂，c∈C。

为嫌疑树干中心点集C构建Kd-tree索引结构，首先设置树木半径r，遍历每个中心点集中的点c_i，c_i∈C，以该点为查询点，搜寻其邻域半径内的邻域中心点c_j，c_j∈C，若为空集(除查询点本身)，则保留；若非空，则依次查找各个邻域中心点c_j，比较该邻域中心点c_j与查选点c_i之间的高程区间累积值E_c，若查询点的累积值为该邻域内所有中心点的最大值则标记为树木中心点，并将其他点标记为非树木中心点，该步骤即获得k为该中心点在中心点集中的索引值；若查询点不是最大值，则标记为非树木中心点，在遍历完所有点后，将该点排除。因此，在经过极大值抑制后得到的树木中心点集为

5.树木提取

将原始点集P投影于XOY平面上，从而获得投影点集并为该点集建立Kd-tree索引。然后遍历每一个树木中心点c_tr(c_tr∈C_tr)，每遍历一个树木中心点，以该中心点作为查询点，搜索投影点集中距离小于r的数据点并记录该索引值i，构成单株木的点云索引序列。最后在原始点云中将属于该索引序列的点集提取出来作为单株树木点云，完成提取树木操作。

如图2a所示的是本实施例原始点云数据,如图2b为本实施例树木提取结果示意图，图2b中灰色为背景点，黄色为植被提取点。如图3为实施例中不同树种的单株木提取结果示意图，由左往右依次为马尾松、杉木、阔叶林的提取结果。

综上所述，本发明可以快速，自动提取林区场景三维点云数据中的单株树木点集，提取的单株树木包含树干和树冠结构。与现有技术相比，本发明采用高程区间累积值进行树木定位，在树干数据部分丢失的情况下依然可以有效提取，由于采用了累积量作进行计算，克服了点云数据密度远近(距离扫描仪)分布不均的缺点，使得结果受采集设备影响小，更具稳定性。此外，由于采用局部极大值方法进行树木检测，可以更好地区别地面与杂草的干扰，避免了形态特征在不同树种间计算的不稳定，对多种树木的自动提取都具有有效性。本发明采用网格块为单元进行计算，计算速度快，更适用于大规模点云场景，具有实际意义和应用价值。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。