基于地面激光点云的树木胸径快速量测方法

刘静云

（上海工程技术大学 城市轨道交通学院，上海 201620）

0 引 言

在林业资源调查中，树木的胸径（Diameter at Breast Height，DBH）是一项关键基础数据，其对树木生长趋势、生物量估算及森林的碳循环研究有着重要意义。树木胸径通常指树高1.3 m处的直径数值。传统的胸径量测方式常采用卡尺直接量测或使用皮尺量测圆周大小，并以此计算树木胸径数值。然而这种量测方式时间长，且容易产生由人工目视带来的量测误差。随着测量技术的发展，全站仪、经纬仪等量测设备不断应用于林业调查中。但是，由于林下环境复杂，相关森林参数难以被精确获取，这些设备的应用仍存在一定的局限性。

地面激光扫描（Terrestrial Laser Scanning，TLS）技术近些年发展迅速，具有扫描范围大、采集速度快、数据质量优、自动化程度高等优点，非常适合于获取精细的林下结构数据，现已成为森林调查、林业管理和生态系统研究的重要方法。

近年来，许多学者针对三维激光点云展开了树木胸径量测的相关研究。大多数方法的量测思路是：先对点云数据高程归一化处理，再取胸径高度数据将其投影到二维平面中，借助圆或椭圆拟合方法，进一步量测出圆的直径、椭圆的长短半轴数值作为胸径值。常见拟合方法有Hough变换、最小二乘拟合、基于随机采样一致性（Random Sample Consensus，RANSAC）的圆拟合等，但这种二维的方式存在倾斜树木投影时很难被检测与拟合的缺陷。部分方法还采用圆柱拟合的方式完成树干的检测，得到圆柱模型的轴向与半径参数，以圆柱的直径作为树木胸径值。此种方式相较于二维的方法能够避免投影所带来的问题，鲁棒地拟合出树木模型，但其拟合的精度略有降低。此外，多数研究关注于树木胸径的量测精度，却很少兼顾方法的计算效率，而实际应用中不仅需要考虑算法结果的精度，其运行效率也应当被纳入评价体系中。

针对上述问题，本文提出一种自动化的树木胸径检测与量测方法。该方法先通过布料模拟滤波（Cloth Simulation Filter，CSF）算法，得到点云数据的地面模型，并对原始数据进行高程归一化；再通过点云切片与聚类自动获取潜在的树干点云簇，最后基于RANSAC的圆柱拟合方法，检测点云中的树干并计算树干胸径。

1 方法与原理

本文提出的方法基于地面激光点云实现对树木胸径的自动化量测，是一套端到端的方法，整个量测流程可概括为地面模型生成、高程归一化、点云切片、点云聚类、圆柱拟合等5个步骤完成。

1.1 点云数据预处理

树木胸径参数量测计算时，与地面点云的数字高程模型（Digital Elevation Model，DEM）密切相关，地面的起伏不可避免地会对胸径量测高度的起点产生影响。为了统一点云数据中树木的量测起点，本文对原始点云进行高程归一化，以消除地形带来的影响。

1.1.1 地面模型生成

点云数据的高程归一化关键在于获取原始点云地面数据的DEM。本文选取CSF算法提取样地点云地面DEM。该算法是一种提取点云地面模型的算法，其原理是先将点云数据反转，再模拟出一块布料，将布料置于点云上方，待其自然落下，布料的形状即为网格模型。由于该算法不直接处理非地面点的数据，因此运算速率较高，适合用来生成点云数据的DEM；CSF的运行时间受输入数据的点数影响较小，可直接将原始数据作为输入。此外，该算法参数设置简单，只需要通过调整布料柔软度、运算迭代次数、生成栅格大小等参数即可生成不同的地面模型。为了提取更加准确的DEM，本文设置CSF的参数为：布料柔软程度为2，进行陡坡处理，迭代次数为500次，栅格大小为0.05 m。

1.1.2 高程归一化

通过CSF算法获取点云的DEM后，先将原始点云投影到DEM上，获取每个点投影后的高度，再对原始点云中所有点减去其投影后的高度，重新生成点云数据，即高程归一化后的点云数据。高程归一化前后的点云效果如图1所示。

图1 高程归一化前后的点云Fig.1 Point cloud before and after elevation normalization

1.2 自动化胸径量测

1.2.1 点云切片

由于树木胸径是1.3 m处的树干直径，为了获取有效的点云数据，本文通过切片的方式提取归一化后点云数据中地面高度1.2～1.4 m处的所有点。通过点云切片，可大大减少点云数据的处理量，切片区域的点云数据如图2所示。

1.2.2 点云聚类

为了提高点云聚类时的运算效率，首先对切片后的点云数据使用分辨率为0.5 cm的体素网格进行体素化处理；再对点云进行欧式距离聚类，获得若干个点云簇；欧式聚类的距离阈值设置为0.1 m。在欧式距离聚类的过程中，构建八叉树来加速聚类。通过八叉树的结构，为每个点搜索临近点，可提升计算效率。由于聚类后仍会存在少数杂点，可通过设置点数阈值，剔除点数小于100的点云簇，保留点数多的点云簇作为潜在的树木对象。

图2 点云切片Fig.2 Point cloud slicing

1.2.3 圆柱拟合

水平圆与圆柱拟合是常见的树干模型拟合方式，但由于树木并不总是垂直于地面，所以投影得到的结果无法很好地拟合树干。如果将圆投影到与树干轴线垂直的方向来计算树木胸径，则会增加冗余计算，而用圆柱体估计模型，则可以同时获取树干轴线方向与半径。因此，本文采用圆柱拟合的方式以估计树木胸径，并对每个点云簇使用基于RANSAC的圆柱拟合方法来获取圆柱的参数。点云簇拟合圆柱模型如图3所示。

图3 点云簇拟合圆柱模型Fig.3 Point cloud cluster cylinder model fitting

基于RANSAC的圆柱拟合首先计算输入点云中每个点的法线，再选用任意两点及其法线计算拟合圆柱模型，最后将点云中到该圆柱模型的距离阈值以内的数据点设定为内点，记录内点数量，完成第一次迭代。通过不断迭代，将保留最大内点数的模型作为最佳圆柱模型。由于基于RANSAC的模型估计方法仅能估计得到一种模型结果，当数据存在两个模型时则无法迭代收敛，因此本文采用先聚类得到点云簇，再将点云簇分别进行RANSAC圆柱拟合的方式来避免这一缺陷。此外，如果按顺序依次对点云簇进行圆柱拟合，并不高效；而分别对点云簇执行圆柱拟合过程，相互之间并不会产生影响，可以采用多线程的并行计算加速这一过程，提高圆柱拟合的效率。本文通过这种方式，快速检测出每个点云簇中的类圆柱体，将圆柱模型半径的2倍视为树木胸径估测值。

2 实验与分析

本文基于PCL开源库使用C++语言编程实现整个算法，并且基于OpenMP方法实现并行加速运算过程。实验环境为一台Intel Core i7-10750H CPU＠2.60 GHz、32 GB运行内存的计算机。

2.1 实验数据

本文实验数据为FGI开源基准数据集中的两块样地点云数据（图4），其点数均超过1亿点。该数据集收集于芬兰埃沃附近的森林，植被以苏格兰松、挪威云杉为主。每块样地面积为32×32 m，其中包含2014年5～7月所有DBH大于5cm的树木量测数据。

图4 实验数据Fig.4 Experiment data

该数据集使用Leica HDS6100地面激光扫描仪进行扫描，其量测视野为360°×310°，距离扫描仪25 m处的距离测量精度为±2 mm。该数据集从样地中的5个位置进行扫描：中心扫描一次，中心扫描的东北、东南、西南、西北方向扫描4次。边缘扫描与中心扫描之间的距离约12 m。通过标志物，将5站扫描数据进行配准，使其中心区域的树木点云更加完整，极大地减少了树木相互遮挡造成的数据缺失现象。

2.2 实验设计与精度评定

本文将实验数据中的两块样地作为数据输入，经过算法处理后得到树木胸径预估值后，与实测值进行比较分析。此外，将本文算法与Computree软件的测算结果进行了对比分析。Computree是一款面向林业点云数据处理的软件，其中包含类似的树木胸径量测功能。因此，本文使用Computree软件在两块样地中进行树木胸径值的估测，并与实测值进行比较分析。由于该软件无法直接处理点数较大的原始数据，因此在使用过程中先对样地数据执行下采样过程，将数据下采样至空间点间距为1cm，输入点数约2000万，软件运算时均使用其默认参数值。

本文采用均方根误差（）与回归方程决定系数（），评定两种方法对树木胸径量测结果的精度。计算公式如下：

用来衡量实测值与估测值之间的偏差，该数值越小，表明预估值越准确，精度越高。表示估测值与实测值的接近程度，该数值越接近1，则估测值越接近实测值。

2.3 实验结果分析

2.3.1 胸径量测数量结果

两种方法量测与漏检的树木胸径数量，以及样地实测树木胸径的数量以及树木胸径数量的检测率见表1。从表1可以看出，Computree软件与本文方法量测胸径有效值均超过80%，都做到了树木胸径的有效量测。其中，本文方法在样地2中胸径检测率最高，Computree方法量测率均低于本文方法。经过检查，实测数据样地中存在部分林下灌木，扫描仪无法扫描到树干部位，仅能扫描其树冠处点云，从而造成缺少树干胸径点云无法量测。通过分析还发现虽然两块样地均是5站数据配准后的结果，但是边缘处仍存在部分树木胸径处点云呈半月状残缺，导致漏检情况出现。此外，Computree软件计算结果中样地1中一棵树木、样地2中两棵树木与实测值相差过大，将其视为粗差归为漏检树木胸径值中。

表1 两种方法树木胸径估测数量分析Tab.1 Analysis of the estimation number of trees DBH by two methods

2.3.2 胸径量测结果

将Computree软件与本文算法计算的树木胸径估测值分别按上文评价标准进行精度评定，两种方法在不同样地中的回归分析结果如图5所示。从图5中可以发现，两种方法在两块样地中的计算结果均接近于1，均小于3 cm，有较好的预估效果。本文算法在样地2中的估测结果更好，Computree软件在样地1上估测结果最差。此外本文算法的结果相较于Computree软件计算结果的更加接近于1，数值更小。

2.3.3 胸径量测误差分析

本文将两种算法估测值与实测值相减再取绝对值的结果视为单棵树木的量测误差，图6展示了两块样地胸径量测误差的箱型图。从图6中可以看出，两种方法的估测结果误差主体部分均位于0～2 cm之间，少数异常值分布于3～8 cm。Computree软件在两块样地中的最大误差均大于本文方法。其中，样地1中的误差均值与本文相近，样地2的误差均值大于本文方法。

2.3.4 时间性能分析

除精度外，算法的运行时间也是一个关键的衡量指标。由于本文算法采用了多线程加速计算，在样地1中进行了单线程与多线程的时间测试，统计了1～6线程算法运行5次所消耗时间的平均值，如图7所示。从图7中可以看出，使用单线程计算时，算法执行时间最慢（17.36 s）；而多线程运算均产生了时间效率的提升，其中6线程计算耗时最短（10.64 s），比单线程缩短了接近7 s的时间，由此证明并行运算大幅提升了时间效率。

图5 估测DBH与实测DBH的散点图及其关系Fig.5 Scatter plot and relationship between the estimated DBH and true DBH

图6 DBH误差箱型图Fig.6 Box plot of the DBH error

图7 运行时间与线程数的关系Fig.7 The relationship between execution time and the number of threads

此外，本文对两种方法的运行时间进行了统计，两块样地中消耗时间如图8所示。即便Computree软件的输入数据是下采样后的点云数据，但在两块样地中仍消耗时间超过2500 s。本文方法采用6线程并行计算，其量测树木胸径平均耗时10 s左右，在树木胸径量测的时间效率方面展现了一定的潜力。

图8 两种方法的运行时间Fig.8 Execution time of two methods

3 结束语

目前，三维激光扫描技术在森林调查工作中的应用越来越多，基于地面激光点云展开对树木胸径的量测工作具有较强的实际价值与研究意义。本文提出了一套端到端的自动化树木胸径提取算法，能够以多站配准点云数据作为输入，快速精确地获取树木胸径值。通过树木胸径量测对比实验，表明本文算法不仅在精度上优于Computree软件中的算法，又具有明显的时间性能优势。针对多站配准的、数据量较大的林业数据，该方法很好地解决了树木胸径量测的时间成本问题。然而，本文算法也存在许多需要改进的方面，目前该算法面向的林业场景仅适用于树干明显的数据，对于树干严重遮挡的森林还需要进一步探索。