基于LiDAR点云数据的树冠空隙度指数分析

郑莎莎，董品亮，王成，习晓环，吕毅斌

(1.昆明理工大学理学院,昆明 650500； 2.中国科学院遥感与数字地球研究所,北京 100094； 3.美国北德克萨斯大学地理系，丹顿 76203)

0 引言

树冠结构包括分枝结构、冠型以及树体的组成部分，能反映森林生态系统的生物物理过程，在林业方面有着重要意义[1]。国内外学者已对树冠结构进行了大量研究,如樊魏等[2]从泡桐树冠的生物量和叶面积等方面研究了其结构特征； 张海东等[3]通过对大青山油松人工林枝、叶分布和生长垂直变化的测定,研究了其树冠垂直结构特征。从研究方法上看，分形维数被许多学者用于描述树冠结构[4-7]； 但美国著名数学家Mandelbrot[8]指出分形维数存在一定局限性，分形维数相同的对象之间，空间结构表现形式会存在明显差异，故仅用分形维数这一指标不足以描述空间分布差异性，因此需要其他指标来区分不同对象。

1983年Mandelbrot首次提出空隙度“Lacunarity”的概念[8]，并将其用于描述物体结构的间歇或分形结构。目前，作为对象几何结构的间隙大小及其空间分布情况的技术指标被很多学者应用于衡量空间对象或空间形式描述大范围时空现象，它能直观地体现数据的空间分布特性。Plotnick等[9]指出空隙度是描述空间分布形式的一种多尺度分析方法，可用于二进制数据和三维以下的定量数据。由于激光雷达(LiDAR)点云数据可以获取地物高精度的垂直结构信息，能够表达复杂树冠的精细结构[10-16]，因此本文采用LiDAR点云数据来研究树冠结构。当前空隙度指数在林业上的应用多数针对林分冠层结构[17-18]，而利用该指标描述单个树冠的研究还很少，对于如何计算单个树冠点云的三维空隙度还鲜有研究，利用该指标的优势，探索识别树冠空间结构分布差异性的新方法是本文的目标，如何计算树冠的空隙度是本文的基础和关键。

空隙度的计算方法最早由Mandelbrot[8]提出，后来被Gefen等[19]做了改进； Allain等[20]提出估算确定分形对象和随机分形对象的空隙度算法，由于该方法原理简单、计算简便，得到了广泛应用。如Plotnick等[9]计算了5个平移不变性不同的一维点集空间分布形式的空隙度； Dong[21]提出了基于内像素值的空隙度估测方法； Myint等[17]计算了城市影像纹理的空隙度； Frazer等[22]对呈现随机、聚集及均一3种形式的5个不同模拟冠层高度模型进行了空隙度统计分析； Malhi等[18]采用空隙度分析法分析了亚马孙河热带森林冠IKONOS影像的纹理特征； Dong[23]用滑动盒方法计算了一维、二维和三维情况下点集空隙度。上述研究表明，空隙度指数不仅能应用于分形对象，还能应用于现实生活中的非分形对象，不局限于一维和二维数据集，对三维数据一样有效。遥感影像数据具有三维结构信息(坐标值x，y和亮度值z),利用三维滑动盒算法计算其空隙度，因更能体现研究对象的空间格局特征而被普遍应用。激光点云也具有三维结构信息(坐标值x，y和高程值z)，为此本文利用点云数据的空间信息优势，结合改进的三维滑动盒算法，提出了基于体素(体积元素)和三维凸包的树冠空隙度估算方法。实验表明，该方法能有效地区分不同形状的树冠，进而为LiDAR数据应用于森林植被分类提供算法支持，深化森林生态变化过程探索中对树冠结构三维空间异质性的研究。

1 数据与方法

1.1 数据

初步利用模拟数据研究树冠空隙度的估算方法。根据Dong[21，24]使用的3种简单几何模型(锥、半球和半椭球)和改进方法相关研究模拟所需的点云数据。模拟结果见图1。

图1 模拟的不同树冠点云及三维凸包

进一步利用RIEGLVZ-400在野外采集几棵不同树冠的点云，并做相应的预处理。图2中是4棵预处理后的实测树冠点云。

图2 4棵实测的树冠点云

通过地面三维激光扫描仪获取的点云密度很高，一棵树可达上百万个点，为方便后续计算，对点云数据进行简化(滤波)。利用奥地利Riegl公司生产的Riegl VZ-400获取点云数据，使用配套软件Riscan Pro对采集的数据进行预处理，为后续研究做准备。

1.2 方法

滑动盒算法[20](gliding-box algorithm，GBA)是一种计算空隙度的典型算法。三维情况下，实施对象是一个立方体区域，现实中树冠边界组成的立方

体是非规则的，要对其采用滑动盒算法，需找出树冠的边界，再进行二值化处理。

1.2.1 三维数组

首先，构建点云数据的最小包围盒，对包围盒内点云数据进行体素化[25]。通过设定体素单元(简称体元)的长L、宽W、高H,包围盒划分为行数r×列数c×层数l个体元，r，c，l分别指沿X,Y,Z轴方向划分的体素个数。体元所在的空间称为体素空间，包围盒内的体素空间用数组A表示。点云(xi,yi,zi)与体元(Ii,Ji,Ki)关系[26]为

,

(1)

(2)

式中：x0，y0，z0分别表示点云在X，Y，Z轴上的最小值；Lx/n1，Ly/n2，Lz/n3分别指体元的长、宽和高； ‘⎣」’表示向下取整；I，J，K的取值范围分别为I∈[1，r]，J∈[1,c]，K∈[1，l]，体素集合记为V； 取A的行数、列数、层数分别为r，c，l，首先生成一个全0数组A。

其次，找出各点云转化到体素空间中的位置，并将其在A上的值设为1。树冠各点与相应体素呈一一对应关系，A(Ii，Ji，Ki)=1,i∈[1，Min],Min是V的点数。

最后，三维滑动盒算法计算树冠空隙度时，只针对树冠表面及其内部进行，利用三维凸包找出树冠边界体素，将边界外部体素设为2。

1.2.2 三维滑动盒算法

当盒子尺寸为R时，表示边长为R个体素个数，将边长为R的立方体盒子放到A上滑动，即先将边长为R的立方体盒子放在A的左上角。立方体盒子覆盖的体元中，研究对象占据体素的个数记作S(称边长为R的立方体盒子的质量)。然后将立方体盒子按先行、后列、再层的方式移动，每次移动1个体元位置，直到整个A上每个位置都被滑动盒覆盖过为止。此时，获得一组随滑动盒移动而获得的质量等于S的盒子数n(S,R)。设ρ(R)是尺寸为R的盒子的最大质量,ρ(R)=R3，S∈[0，ρ(R)]，且当S≠0时，S为整数。将A中可覆盖边长为R的立方体盒子总个数记作NR,则

NR=(r-R+1)(c-R+1)(l-R+1)。

(3)

Q(S,R)是n(S，R)概率分布函数，盒子的尺寸为R时，计算Q(S，R)的一阶矩Z(1)和二阶矩Z(2)，空隙度Λ定义为

(4)

以上称为三维滑动盒法3D GBA[22]。

利用上述三维滑动盒法计算树冠空隙度时，A中体素2不参与运算。为使计算更精确，将盒子分为5类： ①全为体素2，记为Ⅰ； ②含体素0和2，记为Ⅱ； ③含体素1和2，记为Ⅲ； ④含体素0，1和2，记为Ⅳ； ⑤不含体素2，记为Ⅴ。

Ⅰ不参与计算，直接剔除； Ⅴ在边界内，必须参与运算； Ⅱ，Ⅲ，Ⅳ类盒子都处于边界处，需找到某个合适的阈值适当地丢弃部分盒子，使其对空隙度值的影响达到最小。

1.2.3 阈值的选取

阈值的选择考虑以下3个原则： ①含足够体素1； ②含足够体素0； ③含较少体素2。其公式为

Nlacu=Ncan+Nin，

(5)

式中：Nlacu为参与运算的盒子总数；Ncan为边界上参与计算的盒子数；Nin为树冠内部参与运算的盒子数。为确定阈值，利用模拟数据做表1所示的寻找阈值的研究。

表1 边界上3类盒子分布

①NⅡ，NⅢ，NⅣ分别表示Ⅱ，Ⅲ，Ⅳ类盒子的个数；Nborder为边界上盒子数。

由表1可以看出，边界上的盒子包含Ⅱ，Ⅲ，Ⅳ这3类盒子，随着R增大,NR逐渐减小，同时出现如下规律： ①R越大，边界上盒子的比例就越大； ②R越大，Ⅳ所占的比例就越大，当R一定时，NR=NⅣ; ③Ⅱ和Ⅲ的个数逐渐减少，当R达到某值后,Nborder=NⅣ; ④Ⅱ与Ⅲ相比，当R在某个值范围内时(表1中R=14),R越小，NⅡ与NⅢ的差值就越大; ⑤当NR=Nborder时，Nin=0,此时只含Ⅳ; ⑥当Nin=0时，阈值必须慎重选取，否则对空隙度会有很大影响; ⑦当NⅡ>0时，从盒子中总体的分布情况来看，体素0的个数都小于体素2的个数，此时须选出含足够体素0、较少体素2的盒子，使丢弃的盒子对空隙度的影响降到最小。当R取不同值(13，18，27，30)时，Ⅳ中体素0，1，2的分布情况如图3所示。

图3 边界上IV类盒子中体素0，1，2的分布情况

从图3可以看出： 当R逐渐增大时，每个Ⅳ类盒子含体素1的个数逐渐稳定(绿色)。体素0(蓝色)和2(红色)2条线之间有交点，随着R的增大，交点逐渐沿X轴左移。图中当R=27后，没有了交点，在此变化过程中，体素2的个数占单个Ⅳ盒子中体素的比例逐渐增大，体素0的比例逐渐减小； 当R值接近A的最小维数时，每个Ⅳ盒子含有3类体素的情况逐渐相似。利用相同方法可以分析Ⅱ和Ⅲ类盒子的体素分布情况。

据上述分析，按如下方法选取边界处盒子：

1)Ⅳ中选出体素0和1的总数大于2，体素0的个数大于2或体素1的个数大于2的盒子，当NR=NⅣ时，只存在Ⅳ。为不影响计算结果，必须选出盒子参与运算，由于此时盒子中体素1的比例相互接近，计算体素0和2个数的差值，找出差值小于总差值均值的盒子，选出含足够体素0、较少体素2的盒子，将其参与计算。

2)对于Ⅲ，体素个数较少且含有体素1，为减少信息丢失，全部参与计算。

3)Ⅱ的选取同Ⅳ。

1.3 算法流程

首先获取点云数据，然后对预处理后的数据按图4的算法流程进行操作。滑动盒算法采用重叠式的移动方式[22]。

图4 空隙度计算流程

2 结果分析与验证

2.1 结果分析

模拟树冠参数为冠直径10 m，点数1 178个，冠高分别为14.868 5 m,5.702 3 m,10.858 9 m。体素划分尺寸分别为0.329 7，0.328 1，0.517 7； 0.328 9，0.328 3，0.192 2； 0.327 7，0.331 8，0.353 4。盒子选用2体素×2体素×2体素，即R=2，此时lnR=0.693 1。空隙度结果见表2。

表2 3种冠型空隙度

如表2所示，体素划分尺度不同，3种冠型的空隙度就不同。当R=2时，锥型<半球型<半椭球型。此时，含0个体素1的盒子个数最多，8个体素1(即全1或有一个0或2)的盒子数最少甚至没有。理论上认为，本文模拟点云分布均匀，在点云密度适中的情况下，整个盒子中体素完全属于树冠的概率很小。当S=1时，Q(1,1)=1×Q(1,1)=P,Z(2)/[Z(1)]2=P/P2，Λ(1)=1/P，即当R=1时，P表示树冠上点占据整个树冠空间位置的比例，当P趋于+∞时，空隙度值趋近于0。表2中Q(1，1)分别为0.239,0.251,0.237，代表这3棵树冠体素化后树冠点所占比例情况。

图5表示3种树冠空隙度在不同体素划分尺度下的对比。这3种树冠的点云数、体素划分、冠直径等与表2的情况一致。

图 5 不同体素划分尺度下3种树冠空隙度分布

由图5可知，锥型树冠的空隙度分布情况相对较低，半球型树冠和半椭球型树冠的空隙度曲线的变化趋势相似，空隙度接近。

图6呈现了在相同体素划分尺度下3种树冠空隙度曲线分布。

图6 相同体素划分尺度下3种树冠空隙度分布

由图6可知，体元边长0.5 m。当R=2时，锥型和半椭球型的空隙度几乎相同，半球型相对较低； 当R>2时，锥型开始远离半椭球型，靠近半球型。总体上，半球型的空隙度较低，半椭球型较高，锥型介于二者之间。

与图5相比，图6除不同体素划分，其他条件完全相同。从总体上看，3种冠型空隙度分布有共同点： 盒子尺寸R不同，空隙度不同。随滑动盒尺寸R的不断增大，空隙度会逐渐减小。从模拟的冠型来看，半球型的点分布相对较密，故其空隙度值较低。与锥型相比，半椭球型外表面总体变化减缓，而锥型表面收缩较快，相对于半椭球型来说，分布较密，故其空隙度值较半椭球型小。综上所述，利用空隙度能有效分析树冠空间结构。图6显示相同体素划分尺度情况下的比较更为合理，同时也说明了空隙度指标具有尺度效应。图5和图6表明观测尺度(滑动盒)对树冠的空隙度有一定的影响，尺度越小，树冠空隙度越大，反之亦然。

2.2 验证分析

为体现空隙度在体素划分上的尺度依赖性，将体素空间划分为0.5 m和0.2 m这2种尺度，实测结果分别为图7(a)(b)。据树冠参数，将4棵树冠分别近似为： 半球型、锥型、半球型和半椭球型。本文选用2个半球型树冠，虽冠型相同，但树种不同，这主要是为了说明冠型和树种都会对空隙度曲线产生影响，这进一步说明空隙度描述空间分布的有效性。

(a) 体素空间0.5 m(b) 体素空间0.2 m

由图7可知，8条曲线随着R增大，空隙度均不断减小。从定义上看，此结论与Plotnick等[9]提到的空隙度的另一种定义相符合。

从空间的占据能力上看，半椭球型树冠的P值最小，所以其空间占据能力最小，意味着该尺度下，相应空隙量最多。当R较小时，它们的空隙度值差异明显。锥型树冠的空隙度值介于半球型树冠和半椭球型树冠之间，半球型树冠的值相对较小。

从体素划分尺度上看，与图7(a)相比，图7(b)的体素划分尺度更细，曲线值的范围更大。在2个划分尺度下，4棵树冠空隙度曲线分布均相似，半椭球型≥锥型≥半球型。半球型树冠空隙度值在整个R取值范围上几乎小于锥型树冠和半椭球型树冠。此结果与图6模拟数据结果接近。

由以上分析可知，空隙度值依赖于尺度(体素化分尺度和滑动盒尺度)，然而，这2种尺度依赖于数据本身和前期的预处理操作，其中主要的影响因素包括： 树冠原始点云数据的密度、数据滤波处理等。文中采用0.05 m的尺度对于冗余程度较大的点云数据进行滤波。数据采用多大尺度滤波关系到空隙度指数能否有效区分不同的树冠，特别是不同树种的树冠，其叶片和枝条分布的疏密程度会有很大差异，尺度选取不当，可能会出现无法区分或区分困难的问题。利用空隙度指数能否有效地区分树冠依赖于树冠枝条或叶片分布的密度，滤波后点云数据的密度应不大于叶片和树枝分布的密度，否则会造成空隙度值小，空间异质性程度低，空间结构完全不同的树冠的异质性程度差异会非常微弱。

有学者指出具有自相似性特点的对象，空隙度曲线近似于直线[7]。图7(a)中锥型树冠和半椭球型树冠曲线，图7(b)中锥型树冠曲线在一定尺度范围内都近似于一条直线，故可推断这些树冠在一定尺度范围内具有自相似性。图7(a)中半球型树冠曲线，图7(b)中半椭球型树冠曲线在某个R值之后变化相对急剧，当R大于某个值之后，变化又趋于平缓，其值接近甚至相等。这一现象表明当盒子尺寸R达到空间中团状结构尺寸后，空隙度曲线便会迅速下降[7]。文中的实测数据与模拟数据虽然相似，但存在一定的区别，其原因有2点： ①与模拟数据(分布均匀)和实测数据(非均匀)本身特点有关； ②模拟数据具有规则的外轮廓，而实测数据蕴含生态过程中的随机性因素居多。

图6中使用的数据的冠直径和点数均相等，只有冠型不同。图2中4棵真实树冠因树种本身的不同，树冠高度、冠直径、树枝分布格局、叶片大小和叶片分布格局均存在差异，且它们自身所处的环境也会对其本身的树枝和叶片分布存在一定影响，所以图7中相应的4条曲线存在明显差异，体现了前文所说的空间分布不同而造成的影响。图7中半球型树冠在2种划分尺度下线型很相似，区别是分布位置的高低，这可能与所属树种本身有关(它们属于不同树种)。虽冠型类似但空间结构不同，所以空隙度分布不同。通过实测数据的验证分析表明，利用空隙度能有效地分析树冠的结构特征。

3 结论

1)本文以模拟数据为基础，论述了基于体素和三维凸包的树冠空隙度计算方法，通过该算法分析了3种树冠冠型的空隙度。同时，利用实测点云数据做相应验证分析。研究表明，从分析各个树冠的空间格局入手，将空隙度作为描述树冠空间分布状态的一个指标，应用在判断树冠空间异质性程度上，可为森林生态变化过程探索中采用LiDAR数据研究不同树冠的空间分布特征和进行森林分类提供方法支持。文中介绍的方法也为未来计算非规则对象的空隙度指数提供了技术思路。

2)经对比分析，空隙度值的差异除冠型外，树冠内部的空间分布在一定程度上影响着空隙度的大小和空隙度曲线的变化趋势。冠型相似而空间结构不同，空隙度分布也会有差异。但差异的大小和影响程度问题今后需做进一步研究。

3)基于空隙度的尺度依赖性，不同尺度条件下体素划分对空隙度差异的影响需进一步研究(本文仅对比了2个尺度)。

4)采集LiDAR点云数据的角度分辨率应该适当。角度分辨率小虽然有利于获取目标更精细的结构信息，但数据会存在大量冗余，影响计算效率。后期对数据进行的滤波处理是去除冗余数据，提高算法效率的重要关键。本文采用的八叉树滤波，间隔尺度的选取标准是保证利用空隙度指数能有效的区分树冠。对这一问题今后还要继续深入研究。

5)本文方法的关键在于准确区分边界，即如何选择边界处参与空隙度计算的盒子，才能使计算结果达到更高的精度。虽然本文采用的办法已取得较好效果，但对于盒子的选取今后还可探讨其他更好的方式。

参考文献(References)：

[1] 王政权,王庆成,李哈滨.红松老龄林主要树种的空间异质性特征与比较的定量研究[J].植物生学报,2000,24(6):718-723.

Wang Z Q,Wang Q C,Li H B.Characteristics and comparison of spatial heterogeneity of the main species of korean pine in old growth forests[J].Acta Phytoecologica Sinica,2000,24(6):718-723.

[2] 樊魏,赵东,王奇瑞,等.泡桐农田防护林带单株树冠结构特征[J].中国农学通报,2012,28(7):85-88.

Fan W,Zhao D,Wang Q R,et al.Individual tree crown structure of paulownia fortunei(seem)hemsl in Shelterbelts[J].Chinese Agricultural Science Bulletin,2012,28(7):85-88.

[3] 张海东,田有亮,何炎红,等.大青山油松人工林树冠垂直结构特征的研究[J].安徽农业科学,2009,37(18):8739-8742.

Zhang H D,Tian Y L,He Y H,et al.Study on the crown vertical structure characteristic of chinese pine plantation in Daqing Mountain[J].Journal of Anhui Agri Sci,2009,37(18):8739-8742.

[4] Zeide B,Pfeifer P.A method for estimation of fractal dimension of tree crowns[J].Forest Science,1991,137(5):1253-1265.

[5] 陈军,李春平,关文彬,等.林带小钻杨树冠的分维结构[J].林业科学,2006,42(12):6-12.

Chen J,Li C P,Guan W B,et al.Fractal characteristics of tree crown of populus Xiaozhuanica in Shelterbelts[J].Scientia Silvae Sinicae,2006,42(12):6-12.

[6] 刘兆刚,刘继明,李凤日,等.樟子松人工林树冠结构的分形分析[J].植物研究,2005,25(4):465-470.

Liu Z G,Liu J M,Li F R,et al.Fractal analysis of crown structure in pinus sylvestris mongolica plantaiton[J].Bulietin of Botanical Research,2005,25(4):465-470.

[7] 马克明,祖元刚.兴安落叶松分枝格局的分形特征[J].植物研究,2000,20(2):234-241.

Ma K M,Zu Y G.Fractal characteristics of the branching pattern of larix gmalini[J].Bulletin of Botanical Research,2000,20(2):234-241.

[8] Mandelbrot B B.The Fractal Geometry of Nature[M].New York:W H Freeman and Company,1983:310.

[9] Plotnick R E,Gardner R H,Hargrove W W,et al.Lacunarity analysis:A general technique for the analysis of spatial patterns[J].Physical Review E,1996,53(5):5461-5468.

[11]Lee J H,Biging G S,Radke J D,et al.An improved topographic mapping technique from airborne LiDAR:Application in a forested hillside[J].International Journal of Remote Sensing,2013,34(20):7293-7311.

[12]Smart L S,Swenson J J,Christensen N L,et al.Three-dimensional characterization of pine forest type and red-cockaded woodpacker habitat by small-footprint,discrete-return LiDAR[J].Forest Ecology and Management,2012,281:100-110.

[13]Varhola A,Frazer G W,Teti P,et al.Estimation of forest structure metrics relevant to hydrologic modeling using coordinate transformation of airborne laser scanning data[J].Hydrol Earth Syst Sci,2012,16(10):3749-3766.

[14]Vogeler J C,Hudak A T,Vierling L A,et al.LiDAR-derived canopy architecture predicts brown creeper occupancy of two western coniferous forests[J].The Condor,2013,115(3):614-622.

[15]Yi L,Eetu P,Juha H.Investigation of tree spectral reflectance characteristics using a mobile terrestrial line spectrometer and laser scanner[J].Sensors:Basel,2013,13(7):9305-9320.

[16]梁欣廉,张继贤,李海涛,等.激光雷达数据特点[J].遥感信息,2005(3):71-76.

Liang X L,Zhang J X,Li H T,et al.The characteristics of LiDAR data[J].Remote Sensing Information,2005(3):71-76.

[17]Myint S W,Lam N.Examining lacunarity approaches in comparison with fractal and spatial autocorrelation techniques for urban mapping[J].Photogrammetric Engineering and Remote Sensing,2005,71(8):927-937.

[18]Malhi Y,Román-Cuesta R M.Analysis of lacunarity and scales of spatial homogeneity in IKONOS images of Amazonian tropical forest canopies[J].Remote Sensing of Environment,2008,112(5):2074-2087.

[19]Gefen Y,Aharony A,Mandelbrot B B.Phase transitions on fractals:III.Infinitely ramified lattices[J].Journal Physics A:Mathematical and General,1984,17(6):1277-1289.

[20]Allain C,Cloitre M.Characterizing the lacunarity of random and deterministic fractal sets[J].Physical Review A,1991,44(6):3552-3558.

[21]Dong P.Test of a new acunarity estimation method for image texture analysis[J].International Journal of Remote Sensing,2000,21(17):3369-3373.

[22]Frazer G W,Wulder M A,Niemann K O.Simulation and quantification of the fine-scale spatial pattern and heterogeneity of forest canopy structure:A lacunarity-based method designed for analysis of continuous canopy heights[J].Forest Ecology and Management,2005,214:65-90.

[23]Dong P.Lacunarity analysis of raster datasets and 1D,2D,and 3D point patterns[J].Computers and Geosciences,2009,35(10):2100-2110.

[24]Dong P.Characterization of individual tree crowns using three-dimensional shape signatures derived from LiDAR data[J].International Journal of Remote Sensing,2009,30(24):6621-6628.

[25]Zheng G,Moskal L M.Computational-geometry-based retrieval of effective leaf area index using terrestrial laser scanning[J].IEEE Transactions on Geoscience and Remote Sensing,2012,50(10):3958-3969.

[26]邱华.三维体数据生成及三维缓冲区分析[D].长沙:中南大学,2011:1-67.

Qiu H.Three-dimensional volume data generation and three-dimensional buffer analysis[D].Changsha:Central South University,2011:1-67.