一种深度学习网络的树木点云骨架重建方法

徐风， 张博， 袁星月， 安锋， 陈帮乾， 周立军， 云挺*

(1.南京林业大学信息科学技术学院， 南京 210037； 2.南京林业大学机械电子工程学院， 南京 210037； 3.中国热带农业科学院橡胶研究所， 海口 570000)

树木骨架空间模型的精准重建在研究林木资源信息化和表型结构特征反演中发挥着至关重要的作用[1-2]。近些年，针对树木模型重建的技术大致可以分为两类，分别是基于图像与激光点云的树木重建。这两种技术在数据获取上存在着差异性，但提取树木空间特征的方式上又具有一定的相似性。

基于图像的树木枝干重建方法通常采用双目视觉[3]、多目视觉[4]及深度相机的方法捕获树木的纹理与景深信息[5]，结合曲率约束实现树木枝干重建，或利用多目视觉的多张图像通过交互式编辑重建树木枝干[6]。基于图像的树干重建方法虽然简单易行，但提取的树木枝干停留在二维层面，很难直观反映树木的三维几何形态，且提取结果会被多样的环境背景、复杂的树干空间拓扑结构、立体匹配、误差、视角遮挡等多方面因素干扰[7]。

激光扫描具有高精度、稠密度等特性[8]，同时，可结合计算机图形学、机器视觉等理论算法开展林木三维点云的处理分析，例如，利用图论将采集到的点云数据存储为八叉树结构[9]，采用结合点法线的PROSAC算法建模提取林分树干[10]；或使用空间填补的三角面片构造枝干表面来建立树木冠层模型[11]；或通过建立方向场引导并结合图连通算法从激光点云重建树木骨架[12]；其他还有利用提取枝干特征点构建模型的方式，主要有：采用水平切片[13]、超体素分割[14]、基于矢量的枝干轴线提取模型[15]和基于拉普拉斯的骨架收缩方法[16]，从稠密的扫描点云中找寻局部的方向性[17]、结构性[18]、与不同的空间特征[16]也是树木枝干重建的要点。与基于二维图像的树木枝干重建方式相比[7]，基于激光点云数据的树木建模得到的树木几何精度高，普遍能够较好地呈现出树木的空间真实状态。但也存在一些问题，例如，树木具有多个一级二级分枝，冠层内部拓扑结构复杂，不能精准分类叶子和树木枝干点云分类精度不高等，这些因素对树木骨架后期重建会产生较多干扰。

近年来，面向点云处理分析的深度学习网络逐渐被学者们关注，如PointNet[19]、VoxelNetVOX[20]、PointCNN[21]、分层K-d树ContextNet[22]和SPLATNet[23]。前几种网络目前主要应用于工业中产品和规则几何物体的检测，以及车辆自动驾驶道路上目标识别，而应用于复杂拓扑结构的林木点云中的器官识别并不多见。同时，目前深度学习算法处理复杂点云还有很多未解决的问题，主要集中在3个方面。

(1)深度学习既需要考虑识别物体的全局特征，又需要融入物体的局部特征，这涉及局部邻域基团的特征计算与对象整体特征表达在深度学习网络中融合的问题。

(2)空间点云由排列散乱且无序的点构成，难以进行规则的卷积运算。如何面向空间点云开展有效的卷积运算并且提取特征，是待研究的问题。

(3)多尺度分解可以整合目标对象的整体与细节特征，如何有效地将多尺度分解应用到点云的处理分析中是研究的方向。

鉴于以上种种因素，现运用晶格投影策略的深度学习网络与计算机图形学算法，处理地面激光树木点云数据。该深度学习网络构造点云邻域基团，通过旋转不变模块与多尺度晶格投影得到该点在不同尺度下的二维特征，以便于卷积运算，并结合深度学习的池化与多层感知机操作，自动实现复杂林木点云的枝叶分类操作。以海南多类树木为研究对象，开展精准的枝叶分离并重建树木空间枝干的三维模型。同时，对不同级别枝干和相应的叶团簇进行参数提取，并与真实测量数据比对以验证算法的有效性。

1 激光点云数据获取与预处理

1.1 研究区与样本概况

研究区位于海南大学校园(20°3′N，110°20′E)与儋州市橡胶种植园(19°32′N，109°28′E)，试验地地理位置在低纬度热带北缘，属于热带海洋气候，春季温暖少雨多旱，夏季高温多雨，秋季多台风暴雨，全年日照时间长，辐射能量大，年平均气温23.8 ℃，最高平均气温28 ℃左右，最低平均气温18 ℃左右，年平均降水量1 664 mm。选取研究区内多株不同树种作为实验对象，如橡胶树(品种为PR107和CATAS 7-20-59)、香樟树、红枫树、樱花树等。

1.2 数据获取与预处理

点云数据的获取使用Leica Scanstation C10 型号三维激光扫描仪，其扫描角度为360°× 270°，扫描速率为5万点/s，各类树木的数据采集时间为2019年8月15日。考虑到树木的冠层内部结构复杂、树叶相互遮挡[24-25]，为了获取完整的树木三维点云数据，以每棵树冠为中心，分别对其进行对称的2站扫描，每站扫描仪距扫描树冠的中心3 m远，扫描精度设置为中等精度，具体扫描参数为：角度精度(水平/垂直)：中等0.057°/0.057°、视场角(水平/垂直)：360°/270°、扫描点间距：距离3 m时最小<3 mm、标靶扫描精度：3～5 mm。接着采用 Cyclone 软件结合标靶数据对原始两站点云进行了手动拼接，从而得到完整的目标树叶子与枝干点云数据，扫描的点云密度平均为12 000 点/m2。由于获取到的树木三维点云数据中含有噪点，采用文献[26]的方法对采集的点云进行了去噪处理，接着把标记好枝叶信息的点云数据带入晶格投影的深度学习网络中训练。

2 实验方法

2.1 晶格投影的深度学习网络实现树木枝叶分离

图1 枝叶分类的晶格投影的深度学习网络Fig.1 Deep learning network for the classification of wood and leaf point clouds

2.1.1 旋转不变性模块

构造的旋转不变性模块的具体如下：首先将点云数据看作一个整体，分步骤旋转，运用三个旋转矩阵：

(1)

(2)

(3)

式中：RX、RY、RZ分别为绕X、Y、Z轴的旋转矩阵；α、β、γ分别为绕各自轴的旋转角度，输入以pi为中心点，直径为0.1 m呈立方体领域的点云基团，该点云的坐标为n×3的矩阵pj(xj，yj，zj)，j=1，2，…，n，矩阵pj经旋转变换后得到新坐标矩阵：

p′j=pjRXRYRZ

(4)

旋转变换结果如图1(b)所示，一个不规则放置的局部枝干和叶子经旋转变换后，枝干的轴向几乎与Z轴平行(即直立)，且叶子平行于XOZ平面。

2.1.2 晶格投影与重心插值

将旋转后的树木局部点云p′j作为输入置于一个立方体之中。该立方体的XOY平面、XOZ平面和YOZ平面上都排布有晶格顶点，单个平面上的顶点有共σ2个，呈σ×σ二维矩阵状规则排列。这里σ为晶格尺度，它是一个常数，且三层晶格投影模块的σ不相同，是逐步减小的。定义第一层的σ为10，那么第二层的σ为「10/2⎤，第三层的σ为「10/4⎤，如图1 (a)所示。「⎤指向上取整。晶格的大小是通过晶格尺度σ来控制的，σ的变化会改变晶格顶点之间的距离，可以类比为对晶格进行等比例放大的操作。σ越大，晶格尺寸越大，点云细节表现更多。

(5)

式(5)中：mod为求余操作；{}内左边为被除数；右边为除数；「⎤指向上取整；⎣」指向下取整。

晶格投影完成后，使用重心插值将对投影生成的点进行溅射操作，将一个点的数据散开至最近的三个晶格顶点处。

(6)

式(6)中：只有bj，0、bj，1和bj，2三个未知数，解得插值系数为

(7)

以一个局部枝干和一片局部树叶为例，上述操作结束后能够得到如图2和图3所示结果。图2(d)与图3(d)是使用到的树叶局部点云图与枝干局部点云图，图2(a)～图2(c)是树叶点云在不同晶格尺度与不同平面的晶格投影与重心插值结果，图3(a)～图3(c)是枝干点云在不同晶格尺度与不同平面的晶格投影与重心插值结果。各图像对应的晶格尺度σ已标注在分图题中，对应的投影平面已标注在各图像左上方。三条颜色条分别对应10、5和3这三种晶格尺度下，晶格投影与重心插值结束后，晶格点所累加的插值系数大小。

将所有点在同一个平面上生成的重心插值系数存储在同一个矩阵内，因为投影到三个平面，且平面数量与存储插值系数的矩阵数量相同。所以存储插值系数的矩阵编号k=0、1、2，共有三个存储插值系数的矩阵，每个矩阵记为

图2 树叶点云的晶格投影图Fig.2 Schematic diagram showing the lattice projection for leaf point clouds

图3 枝干点云的晶格投影图Fig.3 Schematic diagram showing the lattice projection for branch point clouds

(8)

2.1.3 激活函数与数据集

生成的输出数据根据使用到的晶格尺度的不同，对应使用不同的多层感知器进行卷积。具体来说：晶格尺度为10时，得到的3×n×3的输出数据全部输入卷积模板为MLP(3， 16)的多层感知机中，卷积后得到3×n×16的数据。晶格尺度为5时，得到的3×n×3的输出数据全部输入卷积模板为MLP(3， 32)的多层感知机中，卷积后得到3×n×32的数据。晶格尺度为3时，得到的3×n×3的输出数据全部输入卷积模板为MLP(3， 64)的多层感知机中，卷积后得到3×n×64的数据。接着链接三组输出数据，得到一个3×n×112的特征，进行累加，得到n×336个特征。再使用全局池化操作，形成一个336的向量，最后使用多层感知机MLP(64，32，2)与激活函数，得到2×1的枝和叶置信度值。

枝叶分类问题是一个二分类问题，为了处理该问题，深度学习网络在最后使用sigmoid作为激活函数，计算公式为

(9)

sigmoid函数将网络输出的2×1的枝和叶置信度值转化至0～1，并取较大值代表对应输入点云的类别。

数据集由4种树木组成，分别是橡胶树(CATAS 7-20-59)、橡胶树(PR 107)、香樟树和红枫树。橡胶树是热带树种，具有很大的经济价值， CATAS 7-20-59与PR 107相比，较为低矮。香樟树是亚热带树种，树形高大，多用作绿化与美化景观。红枫是亚热带树种，树形优美，多用作观赏植物。数据集树木参数如表1所示。

使用机器学习[27]，辅助人工标注对数据集进行了枝叶分离，给每个点标好标签以便深度学习网络训练时使用。部分枝叶分离后的训练集如图4所示。

表1 深度学习数据集树木参数Table 1 The tree parameters of the deep learning dataset

棕色为分割后的枝干点云; 绿色为分割后的叶子点云图4 部分深度学习训练集 Fig.4 Partial training sets for our deep learning network

2.2 树木骨架重建

2.2.1 树干高度分层与中心点求取

图5 部分橡胶树干点云及每层聚类中心点提取Fig.5 Part of the scanned rubber tree branch point cloud and the corresponding extracted central points of each layer

2.2.2 枝干链表确立及圆柱体拟合

图6 树木骨架边缘点与分枝节点提取Fig.6 Extraction of the bifurcation and edge points from each tree skeleton

直到该链表遇到边缘结点时，从根结点开始的链的中心点被分类为主枝干结点。同时，在主枝干链上延伸的分枝结点称为一级分枝结点，进而定位整个树干的其他一级分枝。

(10)

(11)

式(11)中：q1、q2分别为拟合直线L1上的任意2个点，‖‖2为2范数。

重复上述过程，即从当前枝干段上随机采样两个点，并不断优化得到d的最小值，直至d不再变化。

图7 生长角度变化最小准则确立主枝干链表Fig.7 Tree trunk chain determination using the rule of minimal variation in the growth angle

图8 自适应半径计算的圆柱体分段枝干点云拟合Fig.8 Adaptive radius calculation of cylindrical subsection branch point cloud fitting

2.2.3 一级枝干对应叶团簇分割

最后，以分离的主枝干和各一级分枝扫描点云为聚类中心，根据空间分水岭聚类算法[28]，实现不同枝干所对应的叶子点云分类，并运用Alphashape算法[30]计算单株树不同分类点云对应的空间叶团簇的体积。

3 结果与分析

3.1 枝叶分离与骨架重建

深度学习网络运算在英特尔i7-7700 CPU@2.80 GHz处理器和16 GB内存的Windows 10 64位PC上执行，在PyCharm和MATLAB软件平台上执行了晶格投影策略的深度学习网络和树木重建程序，同时使用NVIDIA RTX 2080Ti GPU代替CPU来减少训练时间。在晶格投影策略的深度学习网络模型中，学习率为0.000 1，批量大小设置为16，迭代次数为300。设置的总训练时间约为80 h。

具体训练过程中的准确率和损失函数值如图9所示，随着学习迭代次数的不断增加，训练样本的分类准确率呈上升趋势，而损失函数呈下降趋势，这表明所使用的深度学习模型参数满足全局优化且收敛。同时，训练期间神经元网络在每个批次中遇到了一些复杂的样本，如点云存在局部遮挡和数据缺失、一些叶片扫描不全呈现枝干形态而被当前的深度网络权值进行误判等，进而导致回归损失函数值的局部波动。

在经历了100个Epoch之后，准确性和损失函数值分别到达0.90和0.05，表明模型参数较优。通过晶格投影的深度学习网络的分割结果如图10所示。

绿色代表叶子点云；棕色代表枝干点云图10 不同树木枝叶分割结果图Fig.10 Wood-leaf classification results for scanned points of different tree species

深度学习方法与传统机器视觉算法枝叶分类性能对比如表2所示，分别在分类精度、IOU与分类时间三个参量上阐述算法的优劣。其中IOU定义为

(12)

面向不同树种枝干数据，根据每层高度的中心点与枝干扫描点云进行圆柱体拟合以构造树木骨架模型，部分结果如图11所示。

表2 机器学习与深度学习的点云枝叶分类性能表Table 2 Performance evaluation of the wood-leaf point classification using a deep learning approach versus the machine learning method

3.2 不同分枝与叶团簇提取

面向不同树种，运用本文方法提取的不同分枝及对应的叶团簇如图12和13所示。提取的不同一级枝干显示为不同颜色，而用不同颜色的最小凸包来展示相应分割的叶团簇。

3.3 树木参数反演结果比较

表3罗列了单株红枫树和香樟树以及不同品种的多棵橡胶树的生长参数实际测量值和本文方法的计算值，其中包括树高、冠积、叶团簇体积、胸径、一级枝干、二级枝干直径和主枝干与一级枝干之间夹角，同时，把实地测量值与本算法计算值进行对比验证，定量化表达本文方法的有效性。

图14给出了运用本文方法与实地验证的橡胶树比较结果。图14(a)显示了4棵PR107基因型和4棵CATAS 7-20-59基因型橡胶树计算的一级枝干直径与实测值的比较结果，具体为PR107 (决定系数R2=0.93， 均方根误差RMSE=0.77 cm， 相对均方根误差rRMSE=5.75%)，CATAS 7-20-59 (R2=0.92， RMSE=0.51 cm， rRMSE=7.21%)。 图14(b)显示PR107基因型和CATAS 7-20-59基因型橡胶树主枝干与一级枝干夹角的计算与测量比对结果，可以看出PR107型橡胶树具有较大的分枝角，与实测值比对为(R2=0.92， RMSE=4.86°， rRMSE=7.56%)；而CATAS 7-20-59型橡胶树的分枝角相对较小，反演结果为(R2=0.91， RMSE=1.77°， rRMSE=8.06%)。实验整体结果表明算法在树木枝干直径与分枝角度估算上具有较高的精度。图14(c)显示了两种橡胶树一级枝干直径和对应的叶团簇体积呈现正相关性，其中PR107型橡胶树具有更多的一级枝干和对应的叶团簇，整体树冠结构为发散形，而CATAS 7-20-59具有较少的一级枝干，树冠结构为倒花瓶形。同时，图14(c)表明较粗的枝干可以支撑起较大的叶团簇，这与更多叶花果等器官需要更粗的枝干进行营养传输和重力支撑的原理相吻合。

不同颜色代表不同枝干分段圆柱体拟合结果图11 圆柱体分段拟合的树木骨架重建结果图Fig.11 Skeleton reconstruction results of tree species using many cylinder models to fit each segment of tree branch points

不同颜色对应不同一级分枝图12 不同树木一级枝干提取结果图Fig.12 Every first-order branch extraction results using our approach for different tree species

棕色代表枝干点云其他颜色代表分割后的叶团簇图13 一级分枝对应的叶团簇分割结果图Fig.13 Foliage clump borne on the first-order branch of the same tree was segmented

表3 本文方法获得的林分参数与实地测量值对比Table 3 Tree property retrieval using our method versus field measurements

图14 算法结果与实测数据的比对分析Fig.14 The retrieved growth properties by our method versus the field measurements

4 结论

树木骨架的精准重建对于分析树木的表型结构、树体自身的物理特性、地理环境因子的影响都起到了信息支撑的作用。本文研究面向地面激光雷达数据开展多类阔叶树骨架形态的建模，通过面向枝叶点云分类的深度学习模型设计、空间链表结构的主枝干与一级枝干查找、参数自适应调整的圆柱体拟合、空间分水岭分类的叶团簇提取，精准且智能的重建了林木的枝干骨架结构与参数获取。

本文研究设计了晶格投影策略的深度学习网络，并融入了点云局部基团特征、旋转不变模块与多尺度晶格投影等模块，解决了局部特征提取、空间点云卷积和多尺度特征融合等问题。

虽然深度学习网络枝叶分类准确率达到91.31%，但导致准确率降低的因素主要有：实际扫描中枝叶互相遮挡，造成数据缺失；叶子和枝干点云点云紧密贴合，特征混淆而难以分离；点云训练集与测试集中存在噪声，影响了深度学习网络最终分割性能，这些都是下一步研究主要提升的方向。