基于地面激光雷达点云的树木冠层变化定量分析

杨水荣，黄洪宇，唐丽玉，陈崇成

(福州大学空间数据挖掘与信息共享教育部重点实验室，福州大学地理空间信息技术国家地方联合工程研究中心，福建 福州 350108)

0 引言

树木空间结构及动态变化信息，对于树木经营管理、 生态环境建模和植物表型学的研究等具有重要意义. 森林尺度的树木变化检测以森林树木在物候期及人工干扰等情况下生物量及结构变化研究为主[1]. 目前对于单木变化检测大部分侧重于树木生长生理变化，借助微树芯、 生长锥[2]等方式进行单木对象的变化检测对树木有一定破坏性，且树木径向生长测量精度不高. 地面激光扫描仪(terrestrial laser scanner，TLS)具有获取单木尺度的树木形态结构信息的能力，能获取高精度的单木枝干、 冠层结构信息和对树木无破坏性获取等优点[3]. 目前，基于TLS点云数据的单木尺度树木冠层变化检测主要从两方面展开. 一方面，通过空间统计分析对不同阶段植物冠层点云的百分位高度进行比较，量化植物在特定状态或时间段冠层移动及节律性变化[4]. 该方式虽然简单，但对实验条件要求严格，需排除风扰动等外界影响因素，难以在常规室外环境中实施，无法体现树冠在表型结构上的变化. 另一方面，通过重建不同时期单木的定量结构模型(quantitative structure model，QSM)[5]或基于体素的枝干结构模型[6]，监测冠层结构参数变化情况，反映树木冠层扩展[7]及枝干体积量变化[8]. 该方式只适用于落叶或无叶树木，不适用于含有茂盛枝叶的树木对象. 有叶冠层可采用冠层体积描述变化，树木冠层体积包括了冠层轮廓以内所占的体积[9]. 目前已有多种方式[10]进行冠层体积估计. 由于树木冠层结构复杂，规则几何体无法描述真实轮廓，精准测定冠层体积较为困难[11]. 冠层内部结构复杂程度对于地面激光雷达采集数据质量影响较大，且冠层茂密存在遮挡，导致扫描获取的冠层点云数据存在一定的缺失. 目前，研究主要从冠层体积定义出发，以逼近冠层轮廓的体积作为树木冠层体积[12].

近年来，基于TLS点云数据的冠层变化检测研究较少，本研究利用体素概括性，采用基于体素的冠层分析方法从不同参考角度对冠层进行组分变化分析； 根据冠层体积定义，采用切片后点云逐层构成Alpha shape3d模型进行冠层体积估计，通过冠层外围轮廓与点云间的距离关系分析冠层变化情况. 以台风“玛利亚”过境前后园林树种芒果树为例，探究冠层变化定量分析方法的有效性.

1 研究区和数据来源

1.1 研究区概况

研究的案例树木为位于福州大学旗山校区校园内部道路两侧的芒果树，如图1所示.

图1 研究区及研究对象Fig.1 The research area and the research object

1.2 数据获取及预处理

使用Riegl VZ-400地面激光扫描仪，获取树木空间坐标信息及反射率属性. 通过同名点拼接多站扫描数据. 根据现场环境，台风前后扫描仪设站位置如图2所示，扫描设站数均为5站，扫描角分辨率均为0.04°，垂直扫描角范围为30°～130°，水平扫描角范围为0°～360°，多站点云拼接标准误差控制在0.004～0.005 m. 通过对两期树木点云构建KD-Tree数据结构，计算树木点云近邻点间间距，拼接后的原始树木点云间平均间距约为0.003 m. 采用半径滤波器(radius filter)进行噪声滤波，过滤后树木近邻点云平均间距为0.003 7 m(详见图3(a)). 由图3(b)窗口图亦可看出台风对树木冠层的破坏程度，冠层外围叶片点云少，果实点云消失(红点表示台风后点云； 黑点表示台风前点云).

图2 台风前后测站空间位置示意图Fig.2 Schematic diagram of scanner location before and after typhoon

图3 树木两期点云配准及冠层受损情况Fig.3 Two-phase tree point cloud registration and canopy damage

跟踪树木受台风影响时间短，成熟期树木枝干生长变化可忽略不计，分析冠层在台风期间变化情况，进行树木点云的枝叶分割. 由于地面激光扫描仪获取的树木点云数据仅含坐标及反射率信息，参考点云反射率信息，通过多人对台风前后树木对象进行分拣枝干点云与叶片点云并评估分拣结果，台风前后枝叶分类结果如图4所示.

图4 枝叶点云分拣结果Fig.4 Branches and leaves point cloud sorting results

2 研究方法

2.1 冠层组分分布变化分析方法

体素作为三维空间分割最小基本单位，对三维空间的概括性，在一定程度上不受噪声及遮挡现象的影响.

树木点云体素化前，首先提取树木根部地面以上5 cm切片点云，利用最小二乘法拟合圆，求取圆心，并以根节点为空间坐标原点，从不同参考角度分析树木整体变化情况. 点云体素化公式为：

(1)

其中：res指的是体素分辨率； coord(x, y, z)为点云三维坐标； round()函数为四舍五入函数；V(x, y, z)为体素中心坐标.

针对体素分布情况，从不同参考轴及角度分析，首先计算体素相对距离：

(2)

式中：v为点云数量是n的体素集； (xvi,yvi,zvi)为体素集内体素中点坐标；p(xp,yp,zp)为定义参考点坐标. 体素相对坐标轴距离的计算公式如下式所示：

(3)

这里，a为定义参考坐标轴，当a为x轴，xvi=0； 当a为y或z，yvi或zvi=0. 根据体素向量与坐标轴向量的内积公式，体素相对坐标轴角度的计算公式如下式所示：

(4)

其中：β1为体素向量；β2为选定参考坐标轴向量，本次实验取|β2|=100，当a为x轴时，β2为x轴向量，当a为y或z轴时，β2为y或z轴向量.

冠层体素投影图可以反映冠层部分受生长、 环境等因素的影响情况. 两期的点云处理后，对得到的叶片点云进行体素化，然后进行平面投影，每个投影像素包含点云个数Np和体素个数Nv，分析其体素投影密度分布、 冠幅和投影面积变化情况，以下为冠幅计算公式：

(5)

其中：ΔXp和ΔYp分别为投影体素中横、 纵坐标的极值差. 投影面积计算公式为：

P=Np×res2

(6)

这里，Np为投影面体素个数.

2.2 冠层体积变化分析方法

各种冠层体积估计方法有各自的优缺点，适用于不同的冠型. 根据冠层体积的定义，目前常用的冠层体积估计方法是基于凸壳或逐层叠加台体的方法. 图5(a)为常用冠层体积估计模型凸壳模型; 图5(b)～(c)分别为α=1的Alpha shape 3D模型及其切片形式. 由Alpha shape 3D构建原理可知，当参数α=∞时，网格模型将为凸壳模型. 取冠层一定厚度切片点云，根据凸壳模型与Alpha shape 3D算法原理[13]形成图6(a)～(b)两种方法线框图以及图6(c)的两种方法模型叠加图，结合树的形态结构特点的不规则性，Alpha Shape 3D切片模型更加贴近于切片外围轮廓. 相关研究说明通过切片形式能够提高冠层体积估计精度[9, 14]，采用切片化Alpha shape 3D模型估算冠层体积，根据冠型特点及点云分散情况将冠层点云按一定厚度水平切片后，逐层构建Alpha shape 3D模型，并累加模型体积作为冠层体积估计值.

图5 三种冠层体积估计模型Fig.5 Three canopy volume estimation methods

图6 Alpha shape 3D与凸壳法方法对比Fig.6 The comparison between convex hull method and Alpha shape 3D

进一步分析冠层变化点云，分析的基本流程如图7所示. 首先，使用R语言图形库对台风后非枝干冠层点云构建Alpha shape 3D网格模型，并输出该模型. 由于输出模型三角面片未带面片法向量信息，使用开源点云处理软件cloudcompare对网格模型进行表面重采样，产生模型表面高密度点云，并重新计算点云法向量. 计算台风前点云与重建网格模型间的距离(cloud to mesh distance，C2M distance)[15]，通过直通滤波器对台风前冠层点云C2M distance属性进行滤波，获取网格外冠层点云，即冠层变化点云. 基于提取的冠层变化点云，分析变化体积量及变化点云分布.

图7 基于网格模型的冠层变化分析基本流程Fig.7 Basic processing steps of canopy changes analysis based on mesh model

3 结果与分析

根据公式(3)，结合高斯核密度函数统计台风前后体素相对z轴的分布变化情况(具体详见图8). 由图8可知，台风前后树木整体体素密度分布趋势相近，相对z轴0.3～1.5 m区间内，台风后体素密度略高于台风前; 相对于z轴2.1 m外，台风前体素密度大于台风后体素密度. 由此可知，台风过后，远离主干的冠层外部组分受到台风的影响较靠近主干的冠层内部组分要大，冠层内部组分变得更加紧致，而外冠层对应的数据点减少. 根据体素相对轴角度公式，设方位角0°为测站的x轴朝向，并将空间坐标原点移至树木根节点，台风前后体素变化详见图9. 由图9可知，方位角55°～165°区间内，树木体素密度增加，该方位角现场区域相邻为另一株树木，方位角165°～280°区间内台风后密度降低. 这表明同一棵树在不同方位(水平和垂直方向)受到的风力作用的大小有差异，并体现了风场与环境(地形、 建筑、 树木分布)相互作用的复杂性.

图8 台风前后体素密度变化情况(相对z轴距离)Fig.8 Changes of canopy voxel density before and after typhoon (relative z axis distance)

图9 台风前后体素密度变化情况(相对x轴角度)Fig.9 Changes of canopy voxel density before and after typhoon (relative x axis angle)

图10 台风前后冠层变化情况(x-y平面)Fig.10 Changes of crown projection area and crown width before and after the typhoon(x-y plane)

根据树木冠层组分的变化情况，利用枝叶分割提取非枝干点云，体素化后投影至x-y平面，以每个投影像素中含有的体素个数作为其像素密度值，将冠层沿垂直高度三等分后，对冠层上、 中和下部冠层分别取台风前后冠层2.787、 4.186和5.687 m处20 cm厚度的点云切片，通过对台风前后的切片体素化后进行差分，分析其变化情况, 如图10所示. 其中x代表体素栅格距x轴距离，y代表体素栅格距y轴距离(单位：m). 将冠层体素投影像素分为五种变化状况：密度未变化体素，即台风前后同一位置投影像素密度值未发生变化； 密度减少像素，即台风后同一位置投影像素密度减少； 密度增大像素，即台风后同一位置投影像素密度增大，可由台风扰动造成的冠层组分位移造成； 台风后形成像素，即由于受台风影响造成冠层组分位移，这些像素来源可由台风后消失像素及密度变化像素产生； 台风后消失像素，即台风后同一位置投影像素消失. 由图10可知，切片01损失像素个数最多，且随高度增加切片层消失体素个数逐渐降低，而台风后形成像素个数随高度增加而增加，从台风后形成像素个数与消失体素个数可以看出，由于冠层下部冠层枝叶更为茂盛，台风对冠层下部损坏较严重，而冠层上部受风扰动造成冠层位移较大. 密度减少像素个数变化趋势与消失像素个数趋势一致.

表1 冠幅及冠层投影面积

利用公式(5)～(6)分析台风前后树木冠层的投影面积及冠幅变化情况, 详见表1，从中可以看出台风过境后树木冠幅和冠层投影面积均减少.

树冠结构参数主要包括树木冠层枝条数量、 冠层体积和叶面积等. 由于树木对象的冠层稠密度对地面激光扫描仪获取树木点云数据质量的影响，冠层枝条数量除现场验证外，并无更准确的方式. 本研究主要分析台风前后树木冠层体积变化情况.

对于冠层树木的冠层体积估计方法，基于切片形式的Alpha shape 3D冠层模型相较于其他模型更接近冠层外围轮廓，因此本实验采用基于切片形式的Alpha shape 3D冠层估计方法对台风前后树木冠层进行体积估计，冠层体积分别为40.678 3、 39.010 9 m3，冠层体积减少量为1.667 4 m3.

基于网格模型提取台风前后变化点云，实验获取台风前后冠层变化点云如图11所示.

图11 变化点云提取及分布情况Fig.11 Extraction and distribution of change point clouds

通过计算点云与主干平面位置夹角，分析变化点云分布情况(详见图11(b))，冠层变化点云大部分处于距离冠层主干2.0～2.5 m范围内. 考虑变化点云分散程度，选用基于体素的变化体积量估计方式，体素尺寸选择0.05 m，计算得变化体积量为0.786 5 m3. 图11(c)为变化点云体积剖面图，由图可知，台风过后，冠层中下位置体积量变化较大.

4 结语

树木冠层结构及组分分布的变化影响着林木生长发育，而树木冠层的复杂性给冠层变化检测带来一定的困难. 本研究从树木冠层点云本身特点入手，提出从点云再到体素两个尺度上的树木冠层变化定量分析方案. 基于体素具有空间概括性的特点，利用体素化方法，从体素相对空间位置、 体素相对轴距离、 体素相对轴夹角和体素平面投影的角度，分析冠层组分分布变化情况. 鉴于目前未出现一种普适性、 鲁棒性较强的冠层体积估计方案，提出基于切片形式的Alpha shape 3D冠层体积估计方法. 该方法计算原理更接近于冠层实际轮廓. 在构建Alpha shape 3D模型基础上，结合冠层网格模型与点云间距离关系，提取冠层变化点云，并从方位角和垂直方向分析变化点云分布情况. 以受台风影响的芒果景观树冠层变化分析为例，分析冠层在台风影响下冠层组分分布、 冠层体积变化情况，并分析受台风影响变化点云分布情况. 该定量化分析方案可应用于分析植物生长过程的冠层变化情况，有利于探索冠层内各组分生长发育的空间模式，目前仅针对单木对象，未来将进一步扩大研究尺度. 由于不同树木冠层复杂程度不同，采用的冠层体积估计方法有所差异. 鉴于冠层体积精准测量仍是研究难点，故本研究所提出的切片形式Alpha shape 3D冠层估计方法在不同树木冠层适用性有待验证，未来将结合光谱信息等属性进行精度评价.