张玥焜，余文杰，赵习之，吕艳星，冯伟桓，李峥嵘，胡少军

(1.西北农林科技大学信息工程学院，陕西 杨凌 712100；2.北京星闪世图科技有限公司，北京 100085)

树木建模在计算机图形学、虚拟现实、农林业等领域具有广泛地研究与应用价值。一方面，通过树木建模获取其生物量、叶面积指数和产量等表型参数，有助于掌握森林资源信息，并进行高效管理与决策；另一方面，计算机生成的真实感树木模型可增加虚拟现实场景的沉浸感，可直观向非专业人员展示农林业复杂模型，并为机器人喷药等提供重要的空间信息。

树木建模经过近50 年的研究，目前已形成基于点云[1-6]、图像[7-14]、规则或过程[15-19]及基于草图的方法[20-23]，其中最新的文献为杨垠晖和王锐[24]的“树木的真实感建模与绘制综述”。最近，国内学者在该领域的研究非常活跃，WANG 等[25]采用变分优化方法生成了真实感强的树模型；WANG 等[26]提出的基于插值混合模型可从2 棵给定模型中构建新的树模型；WANG 等[27]利用树形空间的统计建模方法，降低了树木建模复杂度，随后WANG 等[28]又提出了一种更加完善和精确的树木建模方法，即通过树形空间对树木模型进行统计分析，实现了从少量树木生成多棵树木的目的；张晓鹏团队研究了基于树的可见点云合成非可见点云数据并生成三维模型的新方法[6]，进一步提出在多视图像生成稠密树点云后采用基于过程的方法构建真实感植物模型的新方法[29]。

早期树木的建模主要采用基于规则的方法，如文献[15]开发的L 系统，随着激光雷达扫描技术的蓬勃发展，最近基于点云的建模方法成为了新的研究热点。文献[17]采用空间殖民算法从目标树冠体中随机分布的点中构建了真实感强的树结构模型。文献[19]和文献[23]拓展了文献[17]的方法，提出了更加高效的方法用以生成三维树模型。由于空间殖民算法提取的骨架点不是从原始点云中进行选择，难以保留树的分枝细节。文献[12]根据图像序列重建树点云，并通过复制可见分枝来生成被遮挡部分的分枝。文献[1]基于地面激光扫描点云生成树的分枝结构，首先利用最短路径方法构建子图，通过连接相邻簇质心的方式提取树的骨架结构，最后对点云缺失部分的分枝进行骨架合成。文献[2]提出了针对多棵树点云的自动全局优化算法，与文献[1]的算法相比，重建时间由分钟缩短至秒。随后，文献[3]又提出了一种新的基于叶簇的树模型生成方法。

上述基于点云的建模方法要求点云中树枝骨架清晰，因此大部分研究采用地面激光雷达获取的高密度树点云，其方法采集到的树点云密度低，难以有效实现重建。文献[4]提出了一种基于双约束的贪婪算法实现了稀疏激光雷达树点云的自动几何重建，但该方法建模鲁棒性较差，难以解决残缺树点云的建模问题。

在树木点云分割方面，可采用基于区域增长、基于分水岭、基于最小生成树和基于规范的分割方法实现[30]。文献[2]将三维点云投影到地面上，并将高密度点识别为根节点，然后利用Dijkstra 最短路径算法求解生成树，最后用零权值去除根边实现了单棵树的分割，但基于最小生成树的方法不适于树间距小且点云稀疏情况下的单木分割。文献[4]针对稀疏机载激光雷达(airborne lidar scanning，ALS)树点云，提出了一种基于规范分割的方法自动实现了大规模稀疏树点云的分割。但上述2 种自动分割方法难以避免过分割现象，无法满足小规模树点云局部分割和深度学习标注过程中对单棵树点云分割精度的要求。

本文主要针对ALS 采集的低密度树点云，面向单棵树点云的较高精度分割与鲁棒性建模应用领域，提出一种交互式分割和建模新方法，避免自动化分割和建模可能造成的过分割及缺失点云建模失败等问题。首先，基于高度图的交互式分割方法，分层从一片稀疏点云中提取出单棵树点云，然后基于改进的空间殖民算法，通过交互式调节相关建模参数实现单棵树的建模，使得原始点云细节特征得以保留，最终获得忠实于原始点云的树模型。

1 点云数据获取

ALS 树点云采用无人机Swiss Drones’ Dragon (载重35 kg，最大飞行高度2 000 m，满载续航时间1 h)搭载激光雷达扫描设备Riegl Vux-Sys (激光脉冲频率最高550 kHz，扫描速度10～200 line/s，位置精度0.05～0.30 m)获取，采集的树点云密度为40～50 点/m2。为便于后期交互式分割处理，首先采用lastools (https://lastools.github.io)进行地面点云和树木点云的提取。如图1(a)所示，ALS 系统具有一次性扫描大规模点云的优势，但采集的树点云内部出现枝干稀疏或缺失等现象(图1(b))，从而影响后期分割与重建。

图1 机载激光雷达扫描系统获取的树点云((a)扫描系统获取的山脉树点云片断；(b)分割出的单棵树点云) Fig.1 Tree point clouds acquired by ALS system ((a) A slice of mountain point clouds with thousands of trees scanned by ALS system; (b) Segmented single tree point cloud)

2 交互式树点云层次分割

由于ALS 树点云具有密度低、枝叶交叉错综复杂等问题，采用自动分割方法会造成过分割，难以满足未来基于深度学习树点云分割研究中对人工标注方面的精度需求。因此，本文提出一种交互式分割方法，可利用高度图和点云切片技术，通过交互式调整每个切片生成精度较高的单棵树点云。主要包括生成高度图、交互式粗选单棵树区域和交互式精细调整切片凸包3 个步骤。

2.1 生成高度图

ALS 系统采集的树点云中树冠部分较多，而一般树冠形状呈现树干中心位置高、周围低且从中心到四周高度呈递减趋势，基于这一特点，可将点云的最高位置投影到地平面上，颜色用对应高度值表示形成高度图，从而便于分辨出单棵树点云。

图2 选取一片ALS 树林点云，点云中任意一点的位置和地面高度分别记为Pi=(Pix,Piy,Piz)和Gi(i=1,2,…,n)，Pmin和Pmax分别为点云包围盒最小和最大位置点，Gmin和Gmax分别为所有树木点云相对地面的高度最小值和最大值，Δs为分辨率，将点云离散化投影的二维高度图长宽分别为W=(Pmax…x-Pmin…x)/Δs和H=(Pmax…z-Pmin…z)/Δs，则对应高度图生成算法为：

图2 基于机载雷达点云树点云生成的点云高度图 Fig.2 Height map generated from ALS tree point clouds

根据算法1 生成的高度图(图2)，可清晰发现每棵树的大致中心位置和覆盖区域，为下一步进行单棵树的粗选提供了依据。

2.2 交互式粗选单棵树区域

在生成的高度图中可直观发现单棵树的大致区域，基于本文开发的交互式分割系统，鼠标单击高度图区域选择树冠中心(a,b)，然后在树的边缘位置再次单击鼠标可确定树冠的大致区域半径r，对应到三维空间系统将联动生成一个圆柱形区域，其中圆柱上下底中心位置为(aΔs+Pmax…x,Pmax…y,bΔs+Pmax…z)和(aΔs+Pmin…x,Pmin…y,bΔs+Pmin…z)，圆柱半径为rΔs。计算落入圆柱体内的点并记为分别为圆柱体内点云包围盒最小和最大坐标值点。设定Δh为高度方向切片分辨率，可将落入圆柱体内的点分为个切片，落入第i个切片的点Qi(i=1,2,…,N)对应高度范围为针对切片i内的点Qi，采用单调链凸包算法求得对应凸包，步骤如下：

(1) 将给定Qi按x坐标升序排序，若x坐标值相同则按y坐标升序排序；

(2) 创建凸包上部。将排序后的点按照x坐标从小到大加入凸包S，若新加入点使得S不再是凸多边形，则逆序删除之前加入S的点，直到S重新成为凸多边形为止；

(3) 创建凸包下部。将排序后的点按照x坐标从大到小加入凸包V，若新加入点使得V不再是凸多边形，则逆序删除之前加入V的点，直到V重新成为凸多边形为止；

(4) 合并集合S=S∪V形成最终Qi对应凸包。

图3 右上角蓝色线段形成区域即为切片i中的点形成的凸包，观察可见，粗选点云分割精度低，出现了不相关2 簇点云被分到同一切片的问题。

图3 根据高度图交互式粗选单棵树点云柱形区域生成对应点云切片与凸包 Fig.3 Corresponding 2D convex envelope and slice of a point cloud generated from a rough single tree point cloud in a cylinder space chosen from a height map

2.3 交互式精细调整切片凸包

为解决图3 中点云切片误分割问题，针对凸包，系统增加了鼠标交互式调整或删除凸包多边形顶点的操作，并动态判断哪些点落入了新多边形区域。设定多边形顶点为Si(i=1,2,…,N)，判断点p是否在多边形内部，见算法2。

在调整凸包顶点的过程中，需动态判断Qi中所有点是否在新的多边形内，如图4 所示，若在多边形内部，则全部选定为当前单棵树切片中的点云，至上而下依次调整所有切片如图5(a)所示，并依次连接相邻切片形成如图5(b)所示的三维包络，导出其点云即形成分割后的单棵树点云。

图4 交互式调整切片凸包确定单棵树点云区域 Fig.4 Refine single tree point cloud in a selected slice

图5 连接各层调整后的切片点云形成单棵树点云包络 ((a)调整好的所有点云切片多边形区域覆盖的点即为目标点；(b)对应切片生成的三维包络) Fig.5 Generate 3D hull of a single tree point cloud by connecting adjusted slices ((a) Slices of a tree point cloud; (b) 3D hull generated from the slices)

3 交互式树木建模

由于ALS 树点云密度低、枝干结构不清晰且常在树干部分存在缺失信息，传统的用于高密度树点云的重建算法不能直接用于稀疏ALS 树点云的鲁棒重建，空间殖民算法具有算法简单、鲁棒性好的优点[17]，后面经过改进已成功在基于草图的建模中得到有效应用[19,23]，因此本文采用该算法用于生成树枝骨架结构。

3.1 改进的空间殖民算法

空间殖民算法基于树木生长空间竞争的原理，能较好地提取被冠层遮挡的枝干部分骨架信息。假设p为初始骨架点，对应可生长空白空间点集合为S(p)。空间殖民算法首先需计算骨架点到其受影响的每一空白空间点q之间的方向向量，然后对方向向量求和并标准化为向量n，向量n用于确定新树枝的生长方向，设定骨架点间距为d，则新的骨架点位置可表示为p+dn，然后设定删除阈值dk，若空间待搜索点与新骨架点的距离小于删除阈值，则将待搜索点删除。此外，默认的空白空间点搜索区域是一个球域，容易生成角度过大的树枝分杈且偏离点云位置，因此需改进该算法[5,19]，将搜索空间约束为半顶角角度θp、影响半径为rp的锥形区域，即

其中，vpʹp和vpq分别为上一级树枝骨架和当前骨架的生长方向。改进的空间殖民算法3 流程如下：

3.2 树枝粗度计算与几何表示

空间殖民算法类只能提取生长空间点或树点云对应的三维骨架，为获取骨架对应的树木几何模型，本文基于异速生长模型计算树枝粗度[1]，然后采用广义圆柱体进行树枝几何重建，从而实现树枝的几何表示，树叶根据叶序规则添加[4]，最终得到树木的三维模型。

假设树干最大粗度为droot，当前树枝长度为Lp，其上第i个子树枝长度为Lp,i，则基于异速生长模型计算的当前树枝粗度dp为

其中，Subtree(p)为树枝p上的所有子树枝集合，λ为值在1～2 之间的粗度调节因子，由树种决定。设置树干最大粗度后，采用深度优先搜索算法可计算所有树枝粗度。

3.3 交互式建模

自动化建模鲁棒性较差，默认的参数组合难以解决残缺树点云和自然树木的真实感建模问题。因此本文设计了如图6 所示的用户界面，通过交互式调节约束角度θp、删除阈值dk和影响半径rp等3 个参数，用户可快速生成不同参数组合下形状各异的重建树几何模型。此外用户还可在图6 对话框中交互式设置树点云重采样密度、树干最大粗度droot、树叶类型和叶序规则[7]，通过少量参数设置即可实时生成不同类型的树木模型，为建模人员提供了灵活选择。

图6 交互式树木建模用户界面(框区域内为可交互式调节的约束角度、删除阈值和影响半径参数) Fig.6 User interface of interactive tree modeling (The interactively adjustable constraint angle,kill distance and influence radius are shown in the red rectangle)

4 实验结果分析与讨论

本试验平台为Windows 7 操作系统，计算机配置为CPU i5-7200U 2.50 GHz，内存8.0 G，显卡 NVIDIA GeForce 940MX，所有实验树点云对象均采用第2 节介绍的ALS 系统采集。交互式树点云分割和建模软件演示视频参见链接:https://cie.nwsuaf.edu.cn/docs/2020-07/08679ab69a1d49a3aebf95a6a7f 79214.zip.

4.1 交互式树点云分割结果

基于本文算法，测试了如图7 所示难以自动分割、从单一角度观察难以识别的小片ALS 树点云，实验结果表明，交互式分割能够有效实现单棵树点云的分割。

图7 交互式分割测试结果((a)调整好的点云切片；(b)形成的三维包络；(c)提取的分割后单棵树点云) Fig.7 Result of interactive segmentation ((a) Adjusted slices of a tree point cloud;(b) 3D hull generated from the slices; (c) Extracted single tree point cloud from the 3D hull)

此外，以2 棵相邻激光雷达树点云为测试对象，与文献[2]提出的最小生成树分割(图8(a))和文献[4]提出的规范割分割(图8(b))相比，本文交互式分割方法(图8(c))较好地避免了误分割。

图8 不同方法分割结果比较((a)最小生成树分割；(b)规范化分割；(c)交互式分割) Fig.8 Comparison of various segmentation results from different methods ((a) Minimum spanning tree segmentation; (b) Normalized-cut;(c) Interactive segmentation)

4.2 交互式树木建模结果

图9(a)为一棵ALS 树点云，实验中固定交互式调节约束角度、删除阈值和影响半径中2 个参数，仅改变其中一个参数，重建结果表明：①在删除阈值和影响半径不变的情况下，约束角度从整体上抑制了树枝的走向，在无角度约束时，部分枝条会出现向下生长和枝条开角不合理的情况，约束角度θp越大，枝条开角随之增大，生成的枝条数量越多(图9(b,c))，不合理角度的枝条被修剪，但约束角度不能无限增大，当约束角度大于其实际所有枝条中的最大角度，约束便不再有意义，故找调整树木对应的最佳生长约束角度对树木的几何重建非常重要，经大量实验证明，本系统测试的树点云最佳生长约束角度在80°～100°之间；②在约束角度和影响半径不变的情况下，删除阈值dk主要影响细枝的生成数量，如图9(c,d)所示。

图9 改进的空间殖民算法中约束角度θp、删除阈值dk 和影响半径rp 对建模结果的影响((a)输入点云；(b) θp=85°，dk=2d，rp=20d；(c) θp=95°，dk=2d，rp=20d；(d) θp=95°，dk=3d，rp=20d；(e) θp=85°，dk=2d，rp=10d)Fig.9 The impact of the different combinations of constraint angle,kill distance and influence radius in the improved Space Colonization Algorithm (SCA) ((a) Input point cloud;(b) θp=85°,dk=2d,rp=20d;(c) θp=95°,dk=2d,rp=20d; (d) θp=95°,dk=3d,rp=20d;(e) θp=85°,dk=2d,rp=10d)

为检验采用交互式操作是否能找到较好保留原始树枝点云特征的参数组合，以一棵裸眼能观察到部分树枝结构的树点云为例(图10(a)～(c))，从蓝色区域重建结果可知，通过设置较小的删除阈值和影响半径，可重建出与原始树枝点云更加吻合的树枝几何模型。

图10 不同参数组合下重建结果与原始树枝点云对比((a)原始点云；(b) θp=90°，dk=2d，rp=20d；(c) θp=90°，dk=1.8d，rp=12d)Fig.10 Reconstructed results of branch details with different combinations of parameters in SCA ((a) Input point cloud; (b) θp=90°,dk=2d,rp=20d;(c) θp=90°,dk=1.8d,rp=12d)

为验证改进算法的鲁棒性，测试了如图11 所示的2 组残缺树点云，重建结果表明，虽然缺失大量点云信息，改进的空间殖民算法仍然能跨过缺失点云区域，稳定生成较自然的树枝形态。进一步研究发现，即使对于明显无树形特征的“Stanford Bunny”及“Happy Buddha”点云，通过交互式调整合适参数，亦能重建出兼具点云特征和树形特征的树模型(图12)。

图11 残缺点云重建结果 Fig.11 Reconstructed results of two incomplete point clouds with different parameters ((a) θp=90°,dk=2d,rp=17d;(b) θp=90°,dk=2d,rp=20d)

图12 “Stanford Bunny”及“Happy Buddha”重建结果 Fig.12 Reconstructed results of“Stanford Bunny”and“Happy Buddha”point clouds with different parameters ((a) θp=90°,dk=1.2d,rp=18d;(b) θp=90°,dk=1.2d,rp=7d)

图13(a)和(b)为航拍真实树木照片，图13(c)～ (f)为对应ALS 树木点云，采用本文交互式重建算法，可生成与对应点云和真实照片接近的树木几何模型(图13(e)和(g))。但是，树木几何模型重建质量与树点云的质量紧密相关，因为采集的树点云密度低，所以重建结果与真实照片相比在细节上依然有较大差异。

图13 重建结果与真实航拍树木图片对比((a)航拍区域及对应(b)单棵树木照片；(c)与区域(a)对应的机载激光雷达扫描点云；(d)单棵树点云；(e)采用本文算法对树点云(d)重建的树枝结构俯视图；(f)树点云(d)对应的侧视图；(g)重建结果) Fig.13 Comparison of reconstructed tree model and aerial photograph and point cloud ((a) The area of photograph with (b) A single tree;(c) The corresponding point clouds of the area (a);(d) The top view of the single tree point cloud;(e) Top view of the reconstructed tree model;(f) Side view of the single tree point cloud (d);(g) The reconstructed tree model)

图14(上)展示了6 棵形状不同、点云数量不同(点云数量对应为6 534，4 809，8 407，10 234，2 163和4 462)、密度不同的树点云，设置相同的重建参数(θp=90°，dk=2d，rp=20d)，采用改进的空间殖民算法重建树枝模型并添加树叶后的效果如图14(下)所示，结果表明，改进的空间殖民算法能够较好地重建出树的自然形态。

图14 6 种不同形态(上)树点云对应的(下)重建结果展示 Fig.14 (Top) Six tree point clouds with various shapes and their corresponding (Bottom) reconstructed tree models

在重建效率方面，对图11 和图14 中对应的8棵树点云数量及建模时间进行统计(图15)，可以发现，当点云数量小于5 000 时，重建时间约在1.2 s之内，当点云数量大于等于5 000 时，重建时间明显增长，因此后期需进一步对算法优化，降低时间复杂度。

图15 点云数量与重建时间关系 Fig.15 Relationship between the number of points and the reconstructed time

5 结束语

本文针对ALS 树点云密度低、点云缺失导致的自动分割与建模困难等问题，提出了一种基于高度图和切片技术相结合的交互式分割方法，与最小生成树分割及规范自动分割方法对比发现，手动交互式分割有效减少了树木重叠部分对分割结果的影响，避免过分割问题，且对于单一角度难以识别的复杂点云，可借助多角度调整实现精细分割。

在完成单棵树点云的分割之后，基于改进的空间殖民算法，通过交互式调节约束角度、删除阈值和影响半径等3 个参数的不同组合，可实现树枝茂密程度、光滑度以及和原始点云吻合度的调整，从而最终决定树木的生长形态，此外对于残缺的稀疏树点云，该算法亦能较好地重建出与点云吻合度高的自然树模型。

由于内在生物的多样性和外在生长环境的差异性，自然界树的形态千差万别，很难找到一种通用的方法精确描述自然界中所有的树种，本文通过实验发现，几乎没有任何一组参数，可以让所有树木点云都能得到完美重建。但对于同一树种，其外在的差异性相对较小，今后，一方面将基于学习的方法研究不同树种的参数取值空间，例如根据树种的不同自动选取相应的参数组合，提高重建自动化效果；另一方面，将结合kd树和并行计算优化空间殖民算法，提高重建效率。最后，交互式分割存在费时、费力等问题，今后拟结合基于区域增长或规范的分割方法，在二维切片上先进行自动预分割，然后在预分割基础上进行交互式调整，实现半自动分割，提高操作效率。