基于三维激光点云的树冠体积计算方法研究

李宏星，刘玉春

(1.湖北生态工程职业技术学院园林与建工学院，湖北 武汉 430200； 2.武汉市测绘研究院，湖北 武汉 430010)

1 引 言

三维激光扫描技术作业速度快，测量精度高，能够获取丰富的树木点云数据。高精度的树木点云数据信息量大，对林业、生态环境等研究具有重要的衍生价值，在城市森林调查、园林绿化、树木管理等工作中得到了许多成功应用，例如，绿化数据的采集[1，2]、城市绿地建库[3]、园林绿化工程竣工测量[4]、城市树木特征提取[5]、三维绿量估计[6]等。

树冠体积是评价树木生长状况、生态效益的一个重要因子。树冠是不规则的非实心体，且内部充满空隙，导致树冠体积难以准确定义和计算。常用做法是将包裹整个树冠的最小包络体体积作为树冠的体积。传统方法将树冠近似看作规则几何体来估算体积，精度较差[7]。随着三维激光扫描技术的引入，基于三维点云计算树冠体积成为一种效率高、结果可靠的方法，得到了广泛应用[8～11]。为此，国内外学者提出了许多种树冠体积计算方法。除了前述的简单几何体近似法，还有体元法和计算几何法。体元法将三维空间划分成一个个的小立方体，统计出内部含有数据点的小立方体数量，再乘上小立方体的体积得到树冠体积[8]。这种方法的优点是算法简单，可以排除树冠内部的空洞，缺点是计算量较大，无法顾及因遮挡产生的数据空缺现象。此外，计算的体积与体元尺度相关，体元越大，得到的树冠体积越大。计算计算几何法直接根据三维点云构建不规则形体，并将其体积作为树冠体积，例如凸包法[12]、α-shape法[11]。为了提高这类方法的计算精度，通常将树冠进行分层，每层看作是一个台体，各层体积总和即为树冠体积。分层方法的优点是可较好地避免体积高估，但是层数的划分对估计结果有较大影响。分层法的一个关键是计算每层切面的面积，常用的凸包多边形算法仍会导致高估[13，14]。为克服凸包算法估计结果偏大的问题，笔者提出了一种基于扇形分割的算法[15，16]。由于树冠体积真值未知，无法评价各种体积算法的绝对精度。实际应用中只能根据主观因素(如研究目的和期望的准确性)和客观因素(如树冠特征、时间和成本因素)来选择适当的树冠体积计算方法。这需要对各种算法的特性有清楚的认识，为此，本文分析几种常用的树冠体积计算方法的优缺点，并通过实例对各种方法的计算结果进行比较。

2 树冠体积计算方法

2.1 简单几何体法

简单几何体法将树冠看作一个规则的三维几何实体，如图1(b)，其体积可以几何体的体积公式计算。常用的几何体有圆锥体、半球体、椭球体、圆柱体和抛物体等，它们的体积公式如下：

(1)

式中，V为树冠体积，D和H分别为冠径和冠高。

这种方法只需要测量冠幅和冠高作为输入，计算简单，是传统测树学最常用的树冠体积估计方法。但是，由于树冠形状不规则，近似几何体的选择取决于树种、树木的生物特征，以及研究人员的主观判断。不同树种的树冠形态各异，没有哪一种几何体适合逼近所有树种的树冠。因此，简单几何体法通常比较适用于计算形状比较规则的树冠体积。对于形状复杂的树冠，该方法的估计结果一般比实际体积严重偏大。

2.2 体元法

体元法将树冠所在的空间范围划分成边长为d的小立方体，遍历这些小立方体，判断每个立方体是否包含有点云数据，有则将其视为有效体元并保留，否则视为无效体元进行删除，如图1(c)所示。把最后剩下的有效体元体积之和作为树冠体积，即：

V=N·d3

(2)

其中N是有效体元的数量。

此方法的优点在于适用于任何不规则的树冠，算法简单易于实现，仅需要判断体元的有效性，有效体元为树冠体积部分，无效体元为树冠之外或树冠内部空隙部分。通过删除无效体元可以剔除树冠外围的凹陷部分和内部空隙，避免过大估计树冠体积。但是，体元法的缺点也同样明显。首先是通过点云数据无法区分内部空隙是否是真的不存在植被元素，还是因为被遮蔽导致没有激光反射点。例如枝干内部空间就不可能形成反射点，如果体元设置太小，一些大的枝干内部空间将不被计入树冠体积。另一个问题是算法虽然简单但是比较耗时，需要遍历数量巨大的体元，而且对每个体元都要遍历一次点云数据以判断是否为有效体元。更重要的是体元法计算的体积与采用的体元大小高度相关。大的体元难以刻画树冠结构细节，但可以补偿遮蔽区。小体元能精细描绘树冠结构，但是无法补偿遮蔽区体积。如何确定合适的体元大小是关键，也是一个难点。韦雪花等[8]利用圆柏进行实验指出体元边长的最佳取值为D/10，因其采用的样本平均冠径为 2.39 m，故采用的体元边长为 0.2 m。

2.3 凸包法

凸包法可分为三维凸包法和二维凸包法。三维凸包法构建包含整个树冠三维点集的最小凸包体，以这个三维凸包体的体积作为树冠体积，如图1(d)所示。也可以首先将树冠进行分层，对每层的点云数据构建三维凸包，并计算它们的体积之和。二维凸包法是先将树冠等间隔分层切片，如图1(e)采用二维凸包算法对每层切片构建最小凸多边形并计算其面积，树冠体积采用下列公式计算[9，14]：

图1 树冠体积算法示意图

(3)

式中，h是分层间距，n是切片层数，Si为第i层切片面积。

常用的凸包算法有增量算法、包裹法、Graham扫描法等。这些传统凸包算法的不足是会将树冠边缘的凹陷空间计入树冠体积，导致体积估值偏大，三维凸包算法尤为严重。为了改进二维凸包法，董亚涵等提出了迭代渐进的凸包算法[14]，但改进算法提取得到的外轮廓实际不是凸多边形。对于分层算法，不同的分层间距值将对计算结果产生影响。理论上，这些间距取值越小，体积计算越准确。但是，由于树冠是离散采样，划分过细不仅降低计算效率，还可能得出不合理的结果。林松探讨了不同分层间距对树冠体积计算结果的影响，不同形状的树冠影响方式不同[9]。

2.4 扇形分割法

如图1(f)，扇形分割法基于分层切片的思想，先将树冠分层，再把每层切面分为若干小扇区。取切面中心为极点，扇区内最大极距为扇区半径，对应的数据点作为轮廓点得到切面的外轮廓，切面面积就是这些扇形面积的总和，可表示为[15]：

(4)

式中，m是每层划分的扇区数量，△θ是扇区的圆心角，Si，j为第i层切面第j个扇区的面积，ri，j是对应的最大极距。将式(3)代入式(4)可得体积计算式为：

(5)

除了分层间距h外，扇形分割法中还有一个关键参数是扇区圆心角△θ。当△θ=2π时切面面积最大，为每层数据点的外接圆面积。随着△θ取值减小，扇区划分越多，得到的外轮廓趋近于树冠外缘，得到的树冠体积越精确。然而，与分层间距类似，并非扇区划分得越小就越好，需要根据点云数据具体特征选择合适的圆心角大小。

3 实验结果与分析

3.1 实验数据

实验样本采集自圭亚那和印度尼西亚的热带森林，由荷兰瓦赫宁根大学土地利用和碳排放数据网站(http：//lucid.wur.nl/datasets/terrestrial-lidar-of-tropical-forests)提供[17]。树木点云数据采集使用的仪器是RIEGL VZ-400 3D地面激光扫描仪，该仪器使用的激光波长为 1 550 nm，角分辨率为0.06°，可实现360°全方位扫描，最大天顶角可达100°。表1列出了样树的基本信息，GUY表示圭亚那，IND表示印度尼西亚，共20棵树木，包括11个树种。为了获得高质量的点云数据，每棵数设置了13个扫描测站，每站扫描仪沿垂直轴和水平轴旋转进行了2次扫描。获得的高密度点云如图2所示，其中最小的点云也超过24万个采样点。由图2可以看到这些树木的形态差异较大，树高最大值、最小值分别为 38.3 m、21.2 m，平均值为 29.75 m。利用卷尺实测的冠径最大值为 21.82 m，最小值为 5.5 m，均值为 14.8 m。

表1 样树的基本信息

图2 20棵样树的点云

3.2 结果分析

分别采用圆锥体法、半球体法、体元法、三维凸包法、二维凸包法和扇形分割法计算了20棵树木的树冠体积。由于点云数据量大，顾及计算精度和计算效率，体元法的体元边长取为 0.4 m。对于二维凸包法和扇形分割法，为了确定合适的分层间距，以 0.01 m的步长改变分层间距，得到树冠体积估值的相对变化量，取相对变化量小于给定阈值时对应的分层间距作为最佳值。图3是20棵树木树冠体积相对变化量均值随分层间距的变化曲线，取相对变化量的阈值为0.2%，对应的分层间距为 0.5 m。利用类似的方法我们确定了扇形分割法中△θ的最优取值为2°。各种方法估算的树冠体积如图4所示，表2给出了不同方法结果的线性回归方程、决定系数(R2)和相关系数。

图3 不同分层间距对应的树冠体积相对变化量(步长0.01 m)

图4表明在选用的6种方法中三维凸包法估算的树冠体积严重偏大，对于大型树冠尤为严重，例如第1、7～9、16和20号树木。三维凸包算法将树冠外围凹陷空间计入树冠体积，树冠越大外围的凹陷空间可能就越多，因此，三维凸包法对于大型树冠或外围存在较多凹陷的树冠容易高估其体积。两种简单几何体法(圆锥体和半球体)的估值也偏大，但是两种方法的结果还算比较接近，相关系数达到0.907。第7～9号树木的树冠较为扁平，冠径比冠高大得多，导致两种几何体方法的估值相差较大。对两种凸包法，二维凸包法得到的结果明显小于三维算法的结果，但它们之间的相关系数高达0.97，显著高于其他方法与三维凸包法的相关性。体元法和扇形分割法估算的体积较为接近，在全部方法它们的结果相对偏小。

图4 不同方法估算的树冠体积

由表2可以看出，不同方法计算得到的树冠体积在统计上具有较强相关性，两两之间的相关系数都在0.8以上。表2中的回归模型也展示了不同方法结果之间良好的线性关系，线性拟合的R2值大多都超过0.8。其中，只有体元法和简单几何体法之间的R2值相对较小(<0.8)。原因主要是样本树冠形态较为不规则，所用的圆锥体或半球体不能很好地刻画树冠形状，而其他几种方法基本能反映树冠的轮廓形状。扇形分割法与二维凸包法之间的相关性最强，相关系数为0.989，这可能是由于两种方法均基于分层切片的缘故。采用扇形分割法计算的树冠体积要小于二维凸包法的计算结果，说明扇形分割法能够有效克服传统凸包算法将树冠边缘的凹陷空间计入树冠体积的缺陷。

表2 不同方法树冠体积计算结果的回归方程、决定系数(R2)和相关系数

对于高密度的完整树冠点云，采用较小体元的体元法树冠体积估值可以近似看作真值[9]。从样树点云的采集方式和表1给出的点云数量可知样树的点云密度极高，因此，可以将体元法得到的树冠体积为真值作为参考，来评价其他5种方法的效果。对比5种方法，扇形分割法与体元法的相关系数最大为0.958，说明扇形分割法结果与真值的相关性最强。从图5给出的各方法体积估计值与体元法结果的线性回归关系可以看出，扇形分割法与体元法的结果最为接近。而且它们的回归模型R2值为0.918，也是5种方法中最大的，说明线性拟合度最优。由此可见，扇形分割法是5种方法中最为精确的树冠体积计算方法。

图5 以体元法为参考的线性回归关系

4 结 论

结果分析表明，虽然不同树冠体积计算方法得到的数值结果有较大差异，但是在统计上有显著的相关性。简单几何体法和凸包法计算的树冠体积相对偏大，而扇形分割法和体元法的结果相对偏小。若以体元法结果作为参考，扇形分割法的准确性最高。相比于体元法的时间复杂度，扇形分割法更具优势，是一种能够兼顾精度和计算效率的树冠体积计算方法。