出版日期: 2016-07-25
点击次数：
下载次数：
DOI: 10.11834/jrs.20165061
2016 | Volumn20 | Number 4

技术方法

地基激光雷达的玉兰林冠层叶面积密度反演

王洪蜀^1,2 , 李世华¹ , 郭加伟³ , 梁祖琴¹

1. 电子科技大学资源与环境学院, 四川成都 611731

2. 四川水利职业技术学院测绘工程系, 四川崇州 611231

3. 西南交通大学地球科学与环境工程学院, 四川成都 611756

收稿日期: 2015-03-27; 修改日期: 2016-02-03; 优先数字出版日期: 2016-07-25

基金项目: 国家重点基础研究发展计划（973计划）（编号：2013CB733405）；国家自然科学基金项目（编号：41471294）

第一作者简介: 王洪蜀(1990-),女,硕士研究生。研究方向为激光雷达遥感。E-mail:wanghongshu5013@sina.com

通讯作者简介: 李世华(1977-),男,副教授。主要从事定量遥感、激光雷达遥感等方面的研究。E-mail:lishihua@uestc.edu.cn

中图分类号: TN958.98;S771.8

文献标识码: A

文章编号: 1007-4619(2016)04-0570-9

摘要

叶面积密度LAD（Leaf Area Density）是表征冠层内部叶面积垂直分布的重要参数，其分布廓线的准确反演对研究植被碳氮循环、初级生产力和生物量估算等具有重要意义。本文在电子科技大学校内建立实验样区，利用地基激光雷达Leica ScanStation C10和数码相机获取玉兰林高分辨率3维激光点云数据和真彩色影像。利用监督分类将真彩色影像中枝干等非光合组织与叶片分离，再将像素分类信息映射给点云数据，从而提取叶片点云。通过点云数据体元化，并引入2维凸包算法确定垂直方向分层树冠边界，获取激光接触冠层的频率；随机选择不同高度的多个叶片，利用特征值法进行叶片平面拟合，估算出叶倾角，并结合天顶角估算叶倾角校正因子；最后基于体元的冠层分析VCP（Voxel-based Canopy Profiling）方法实现树林冠层LAD反演。结果表明体元化的叶片点云数据能准确确定树林冠层边界和统计接触频率实现LAD反演；反演的LAD变化走势与区域林木冠层叶片垂直分布相吻合，在冠层中下部随着高度的增加叶面积密度也随之增加，在4 m高度处达到最大值1 m²/m³，之后随着高度的增加叶面积密度逐渐降低。根据LAD计算得到的累积叶面积指数LAI为3.20 m²/m²，与LAI-2200实测的叶面积指数相比，相对误差为1.26%。

关键词

地基激光雷达 , 叶面积密度 , VCP方法 , 阔叶林冠层

Retrieval of the leaf area density of Magnolia woody canopy with terrestrial Laser-scanning data

WANG Hongshu^1,2 , LI Shihua¹ , GUO Jiawei³ , LIANG Zuqin¹

1.School of Resources and Environment, University of Electronic Science and Technology of China, Chengdu 611731, China

2.Department of Surveying and Mapping Engineering, SiChuan Water Conservancy Vocational College, Chongzhou 611231, China

3.Faculty of Geosciences and Environmental Engineering, Southwest Jiaotong University, Chengdu 611756, China

Abstract

As a canopy vertical structure parameter, the Leaf Area Density (LAD) of horizontal layers is generally applied to quantify the leaf area within a given volume. It is important in characterizing the vertical distribution of leaf chemical components and canopy biomass. Light Detection And Ranging (LiDAR) has been widely applied to obtain the three dimensional structural properties of vegetation in the past ten years. Terrestrial LiDAR data can be utilized to estimate LAD. The objective of this studyisto develop a new method for woody canopy LAD estimation through the use of terrestrial LiDAR data. High-resolution point clouds and true-color images of Magnolia woody were acquired with Leica Scan Station C10 and a digital camera device, respectively.The study sites arelocated in the University of Electronic Science and Technology of China. Leaves were labeled through supervised classification of true-color images based on the registration of point clouds and true-color images.The laser-scanning points over leaves were accurately extracted with the classification acquired from the corresponding image pixels. Through the voxelization of point clouds, the two dimensional convex hull algorithm was introduced to determine the boundaries of each vertical layer of tree crowns and the contact frequency of laser canopy. With the Eigenvalue method, leaves at different tree heights were randomly selected to fit leaf planes and estimate leaf angles.Together with the zenith angle, the adjustment factor of leaf angles was identified. Finally, the vertical profiles of LAD were retrieved with the Voxel-based Canopy Profiling (VCP) method. The voxelization of point cloud data can be used to accurately determine forest canopy boundaries and the statistical contact frequency of LAD inversion. The LAD profile is in agreement with the vertical distribution of canopy leaves. The LAD values increased with the height of trees in the middle and lower canopies, and the maximum LAD was 1 m²/m³ at a height of 4 m.The LAD valuesdecreased in the upper canopy as the height increased. The Leaf Area Index(LAI) calculated from LAD was 3.20 m²/m². Comparison with the field measurements acquired by LAI-2200 suggested that the relative error was 1.26%. By combining true-color image classification and canopy point cloud registration, point clouds over leaves can be accurately extracted from original LiDAR data. The LAD of woody canopy can be estimated from point clouds over leaves through the VCP method;in particular, the blade of the point cloud voxel can accurately determine the forest canopy boundaries and statistical contact frequency and can also be used in LAD estimation. However, the accuracy of LAD estimation cannot be evaluated directly because LAI does not reflect the vertical layer of LAD. LAD values need to be investigated in future studies.

Key words

Terrestrial LiDAR , Voxel-based Canopy Profiling (VCP)method , Leaf Area Density (LAD) , forest canopy

1 引言

森林冠层中不同高度的叶片对总光合作用和碳储量的贡献不同，叶面积垂直分布状况是能量、水和养分的重要决定因素，在确定生物栖息地适宜性、物种数量和多样性中发挥着重要作用(Swatantran等，2011)。LAD是表征植被垂直方向上不同分层内叶面积差异的重要参数，其定义为单位地面面积单位高度范围内的单面叶面积总和(Weiss等，2004)，其沿垂直方向的积分为叶面积指数LAI(Leaf Area Index)(Hosoi和Omasa，2006)。LAD可用来表述不同高度处叶面积分布情况，其取决于物种、生长阶段和环境干扰等因素，是生态过程模型在全球碳循环中模拟陆地植被角色的一个重要输入参数(Treuhaft等，2002)。近年来，激光雷达LiDAR(Light Detection and Ranging)作为一种主动遥感技术，发展十分迅速，机载LiDAR和地基LiDAR在森林参数的定量测量和反演上取得了成功的应用(庞勇等，2005；李旺等，2013；刘鲁霞等，2014)。尤其地基LiDAR分辨率高、光斑尺寸小，对植被冠层具有一定的穿透力，能够获得植被冠层内外详细的3维结构信息，是进行植被冠层LAD估测的重要遥感手段。

目前，利用小光斑地基LiDAR主要通过计算激光光束与冠层的接触频率或基于间隙率理论进行LAD估算(赵静等，2013)。Hosoi和Omasa(2006)基于激光光束与冠层的接触频率这一思想提出基于VCP方法用于单木LAD估算，该方法首先配准多个站点扫描的地基LiDAR数据；然后建立立方体元模型，估算每个水平层叶片与激光的接触频率，并构建叶倾角校正因子，从而求得LAD值。为检验VCP方法的适用性，Hosoi等人利用地基LiDAR数据先后成功地提取了榉树(Japanese zelkova，Hosoi和Omasa，2007)，枫树(Japanese maple)和山茶树(Camellia，Hosoi和Omasa，2011)树冠垂直叶面积密度廓线，监测了榉树(Japanese zelkova)冠层LAD和叶倾角的季节变化，结果表明从春季到夏季，冠层上部LAD值增加，从夏季到秋季，冠层中下部LAD值降低。Hosoi和Omasa(2007)评价了基于地基LiDAR数据利用VCP方法对树木冠层LAD反演精度的影响因素，结果表明反演精度受到非光合作用组分、叶倾角取向、入射激光束数量以及激光光束入射角度等因素影响。Béland等人(2011)基于体元与接触频率的思想，利用地基LiDAR数据对草原上的树木进行了叶面积分布估算，通过校正阴影、叶倾角、非光合作用组分以及激光束大小，得到的估计值与直接测量值平均相差14%。Takeda等人(2008)通过将森林冠层的LiDAR点云数据进行体元化，然后基于间隙率模型估算了落叶松(Japanese larch)植物面积密度PAD(Plant Area Density)，与实测数据相比较，结果显示由于遮挡冠层上部有低估的现象。Zhu和Zhang(2011)通过构建3维体元模型，利用间隙率模型估算了单棵树不同高度层的叶面积指数，与实际测量值相比，在植株上部和下部估算值偏高，而在中部，估算值偏低。由此可见，基于体元的思想和接触频率或间隙率理论方法，已广泛地运用于基于地基LiDAR数据的单木和区域林木叶面积密度反演，但目前利用VCP方法进行LAD反演主要集中于单棵树，对区域林木快速进行叶面积密度估算并验证有待进一步的发展与完善。

本文在假设从多个位置获取的激光点云能够最大限度地减少遮挡效应，充分地表达整个林木冠层的基础上进行。首先通过图像分类，实现叶片点云的提取；然后构建体元模型，获得区域林木树冠的体元数据，基于2维凸包算法，确定树冠边界；最后结合基于体元的冠层分析方法VCP实现区域林木叶面积密度反演。

2 地基LiDAR数据获取及预处理

2.1 实验数据获取

本文研究区域位于电子科技大学校园内，在该区域建立8 m×12 m大小的样方，样方内地势平坦，有22棵玉兰树，长势均匀，平均树高为6.6 m。实验使用的地基LiDAR系统是Leica ScanStation C10，该仪器拥有较高的扫描速度及测量级精度，并整合了高分辨率数码相机，能快速、非接触地提供目标场景的3维点云数据及彩色影像，技术参数见表 1。仪器配套的cyclone软件可实现点云及影像数据的预处理工作。

表 1 Leica ScanStation C10 技术参数
Table 1 Technical parameters of Leica ScanStation C10

下载CSV

参数类型	参数值
激光波长	绿色，波长为532 nm
扫描速率	50000 点/s
测量范围	300 m(反射率为90%的物体)； 124 m(反射率为18%的物体)
单点测距精度	±6 mm
扫描视场范围	水平：0°—360°；垂直：-45°-90°
扫描分辨率	有超高、高、中、低等分辨率，点间距最小间隔< 1 mm

为最大限度地降低枝叶遮挡，获取区域林木完整实景，本次实验共设置4个观测站点，测站分布如图 1所示。为了将这4个独立观测站获取的激光点云数据统一在一个坐标系下，布设了5个固定的黑白标靶，并保证这5个标靶的高度不一致；样区中心测站P1能够看到所有标靶并能进行精确扫描，P2、P3、P4测站至少精确扫描3个标靶。在P1测站，进行360°×270°的全景扫描，其他各测站扫描区域水平视场以测站与测站之间有部分重叠为准，垂直视场为-45°到90°。扫描分辨率为高分辨率即100m处点与点的垂直与水平间距为0.05m。激光扫描是在晴空无风的条件下进行的，每个测站扫描完毕后，进行仪器一体化彩色数码相机拍摄，单帧图像为17°×17°，400万像素。

图 1 样区扫描站点及LAI观测点分布图

Fig. 1 The distribution of scan position and LAI observation position in the study site

采用LAI-2200植被冠层分析仪测定样方的LAI，为避免操作者和过强的太阳辐射造成的误差，实验选用180°的遮盖帽。具体设置为ABBBBBBB，即在空旷地方取一个冠上的测量值A，再在样方中，以样方中心为中心点，沿东西南北方向共选择7个冠下观测点B，尽可能准确地表达该样方的叶面积指数。观测点以小圆点表示，其在样区的分布示意图见图 1。为避免偶然误差，进行了3次测量，测得的LAI值分别为3.15m²/m²、3.20m²/m²、3.13m²/m²，平均值为3.16m²/m²。

2.2 地基LiDAR数据预处理

地基LiDAR数据的预处理工作主要包括点云去噪和点云配准。点云去噪目的是剔除两类噪声，一类是与扫描目标无关的点云，如电线杆、建筑物、样方外的其他树木等；二是由于风和仪器信噪比等导致的噪声点云。

对于第一类噪声，可以利用cyclone软件中Fence命令框选择噪声点云，即可剔除；对于第二类噪声，可通过一些去噪算法来删除。本文采用Rusu等人(2008)提出的邻域统计分析技术进行第二类离群噪声点的去除。4个测站点云数据的配准选择基于公共标靶的配准方法。即对场景扫描后，仪器会对标靶进行精确扫描，标靶的中心相当于是两站数据的公共点。通过自动添加约束条件(Auto-add Constraints)，便可实现基于标靶的配准。一般情况下，利用标靶配准就可以得到良好的配准结果，如果偏差大，可采用精配准的方式再次配准点云。本次实验标靶配准误差最大为3 mm，最小为0 mm，满足研究拼接精度需求。

3 树林冠层LAD反演方法

地基LiDAR系统获取的激光点云数据能够真实地呈现目标对象的3维结构，但点云数据杂乱、无章可循且数据量偏大，本文将从叶片点云的提取、冠层叶片点云3维体元模型构建和基于体元的冠层分析方法3个方面建立叶面积密度反演模型，技术流程见图 2。

图 2 叶面积密度反演技术流程

Fig. 2 The workflow of LAD retrieval

3.1 叶片点云提取

叶片点云的提取在于剔除影响叶面积密度估算精度的枝干等非光合作用组分。地基LiDAR系统内置的高分辨率数码相机与扫描仪激光同轴，可利用相机和扫描仪的固定位置关系将相机拍摄的图像映射到点云中去，该步骤可在cyclone软件中实现。因而扫描、拍照后的树木不仅能得到点的3维坐标信息，还能获取点的色彩信息。

由于拍摄视角范围小，每张真彩色影像中叶片与枝干是清晰可见的，可利用监督分类方法将枝干等非光合组织与叶片分离，再将相片像素分类信息映射给点云数据，从而提取叶片点云。

3.2 3维体元模型构建

3维立体图像中的体元与2维平面图像中的像元相对应，将区域林木点云数据体元化，有利于信息表达，减少数据量。本文根据预处理后的点云数据区域范围，以区域坐标X、Y、Z的最小值(X_min、Y_min、Z_min)为起始点，以体元大小为步长，将区域划分为有限个体元，并确定点云在体元坐标系中，对应的体元坐标值及体元值。体元大小由体元的长(l)、宽(w)、高(h)决定，整个点云区域被划分为N_l×N_w×N_h个体元，其中，N_l=(X_max- X_min)/l，N_w=(Y_max-Y_min)/w，N_h=(Z_max- Z_min)/h(Zheng和Moskal，2012)。点云体元化后的坐标值由式(1)得到

$ \left\{ {\begin{array}{*{20}{c}} {i = {X_{{\rm{min}}}}{\rm{ + }}\left({{\rm{int}}\left({X - {X_{{\rm{min}}}}} \right){\rm{/}}l} \right){\rm{ \times }}l}\\ {j = {Y_{{\rm{min}}}}{\rm{ + }}\left({{\rm{int}}\left({Y - {Y_{{\rm{min}}}}} \right){\rm{/}}w} \right){\rm{ \times }}w}\\ {k = {Z_{{\rm{min}}}}{\rm{ + }}\left({{\rm{int}}\left({Z - {Z_{{\rm{min}}}}} \right){\rm{/}}h} \right){\rm{ \times }}h} \end{array}} \right. $

(1)

式中，int是取整符，直接取出小数前面的整数部分，(i，j，k)是点云数据(X，Y，Z)对应的体元坐标。本文中采用的体元大小l×w×h与扫描采用的点间距一致。

体元的体素值通过判断体元内包含的激光点个数来确定，如果激光点个数大于等于1，代表激光束被拦截，体元体素值赋为1，否则体素赋值为0。

3.3 叶面积密度计算模型

选用基于体元的冠层分析方法VCP(Hosoi和Omasa，2006)，进行叶面积密度计算，其公式如下

$ {\rm{LAD}}\left({h，\Delta H} \right)= \frac{{\cos \left(\theta\right)}}{{G\left(\theta\right)}} \cdot \frac{1}{{\Delta H}}\sum\limits_{k = {m_h}}^{{m_h} + \Delta H} {\frac{{{n_{\rm{I}}}\left(k \right)}}{{{n_{\rm{I}}}\left(k \right)+ {n_{\rm{P}}}\left(k \right)}}} $

(2)

式中，ΔH为水平厚度层，θ为水平厚度层ΔH范围的所有入射激光光束天顶角平均值；n_I(k)是第k层拦截激光的体元个数，即体素值为1的个数；n_P(k)是第k层被激光穿透的体元数，即体素值为0的个数；n_I(k)/(n_I(k)+n_p(k))是用于表示激光接触冠层的频率；cos(θ)/G(θ)是叶倾角和激光光束方向的校正因子；G(θ)是在假定叶片水平方位角分布均匀，单位叶面积在垂直于激光束方向平面上的平均投影，其计算公式可定义为

$ \begin{split} G\left(\theta\right)=& \frac{1}{{2\pi }}\textstyle\int_0^{2\pi }\!\!\! {\int_0^{\pi /2} {\rm{g}} } \left({{\displaystyle\theta _{\rm{L}}}} \right)\left| {\cos \left({{n_{\rm{B}}}，{n_{\rm{L}}}} \right)} \right|{\rm{d}}{\theta _{\rm{L}}}{\rm{d}}{\varphi _{\rm{L}}} = \\ & \textstyle\int_0^{\pi /2} {\rm{g}} \left({{\displaystyle\theta _{\rm{L}}}} \right){\rm{S}}\left({\theta，{\theta _{\rm{L}}}} \right){\rm{d}}{\theta _{\rm{L}}} \end{split} $

(3)

$ {\rm{S}}\left({\theta，{\theta _{\rm{L}}}} \right)\!=\! \left\{ \begin{array}{l} \cos \theta \cos {\theta _{\rm{L}}}，\quad\quad\quad\quad\quad\quad\quad\quad\! \text{当}\theta\! \leqslant \!\displaystyle\frac{{\rm{\pi }}}{2} - {\theta _{\rm{L}}}\\ \cos \theta \cos {\theta _{\rm{L}}}\left[ {1 \! + \! \displaystyle\frac{{2\left({\tan x - x} \right)}}{{\rm{\pi }}}} \right]，\;\; \text{当}\theta \! > \!\displaystyle\frac{{\rm{\pi }}}{2}\! - \!{\theta _{\rm{L}}} \end{array} \right. $

(4)

式中，x=cos^-1(cotθcotθ_L)，θ是激光束入射天顶角，θ_L是叶倾角，φ_L是叶表面法线方位角，n_B，n_L分别是激光束方向和叶表面法线的单位矢量。天顶角和叶倾角可从激光点云数据中获取。

为通过获取的叶倾角计算G(θ)，其公式可等价为

$ G\left(\theta\right){\rm{ = }}\sum\limits_{{{q}} = 1}^{Nq} {{\rm{g}}\left({{q}} \right)} {{S}}\left({\theta {\rm{，}}{\theta _{\rm{L}}}\left({{q}} \right)} \right) $

(5)

式中，q为不同叶倾角的类别，Nq为总的类别数，如果以5°为间距，Nq为18，θ_L(q)为q类中间叶倾角值，g(q)为q类叶倾角的分布，是q类的叶面积与总的叶面积的比率。

3.3.1 接触频率计算

激光接触冠层的频率可以通过分别统计每个水平厚度层内，拦截激光的体元个数n_l(k)和被激光穿透的体元个数n_P(k)计算得到。由于冠层结构的不规则性，根据区域点云最大最小值范围来构建3维体元模型，必然存在冠层之外的无效体元。这些冠层之外的无效体元，不能视为是被激光光束穿透的体元，因而确定冠层边界，剔除无效体元，是计算接触频率中的一个重要环节。点集P的凸包(convex hull)是指一个最小凸多边形，满足P中的点或者在多边形边上或者在其内(汪嘉业等，2011)。本文选用简单高效的2维凸包算法Graham Scan 确定每个水平层的冠层外轮廓范围。当点集P内的点数大于等于3时，如图 3，Graham 扫描凸包算法的过程如下(Graham，1972)：

图 3 Graham scan算法示意图

Fig. 3 Schematic diagram of Graham scanalgorithm

(1) 找出Y轴中坐标最小的点，如果存在多个点的Y坐标都为最小值，则选择X坐标最小的一点，作为基点，当作凸包顶点的起点p₀。

(2) 对剩下的其他点p_i，按照其与基准点p₀构成的向量与X轴的夹角大小进行排序。

(3) 以夹角按从小至大排序为例，从基准点开始，依次扫描排序后的点，如果这些点都在凸包多边形上，则每3个相继的点p_i-1，p_i，p_i+1应满足以下性质：p_i+1在向量< p_i-1，p_i> 的左侧。如果不满足上述性质，则p_i点一定不是凸包上的顶点，应删除。

获取每个水平厚度层体素值为1的体元位置，采用Graham scan 算法获取每层区域冠层外轮廓。图 4是激光束在某一水平厚度层被区域内树木叶片拦截的情况，反映了该水平厚度层树木叶片的覆盖情况，黑色的小点代表拦截激光束的体元位置，红色大圆点是这些点集的凸包多边形顶点，即该水平层区域冠层边界顶点。

图 4 区域林木冠层边界轮廓示意图

Fig. 4 Outline of one layer of crown

3.3.2 叶倾角校正

叶倾角的校正cos(θ)/G(θ)意为叶倾角和激光光束方向的校正因子，可通过式(4)、(5)来实现，需要对叶倾角和激光光束天顶角进行计算。叶倾角是叶片表面法线方向与天顶方向的夹角，叶倾角分布直接决定着植被冠层对入射辐射的截获量。以地基LiDAR高精度扫描模式获取的树木点云，其叶片点云是清晰可见的(图 5)。平均每片叶片由600多个激光点组成，能够很好地再现叶片的3维结构，利于叶倾角的估算。可随机选择不同高度的多个叶片，利用特征值法对每个叶片点云进行平面拟合，估计出平面的法线，法线方向与天顶方向的夹角即为叶倾角(官云兰等，2008)。

图 5 单叶片点云

Fig. 5 Leaf point clouds

激光光束天顶角是入射光束与天顶方向的夹角，可通过激光直角坐标和极坐标转换获取，如式(6)(Abbas等，2013)。

$ \left\{ \begin{aligned} r & {\rm{ = }}\sqrt {{x^{\rm{2}}}{\rm{ + }}{y^{\rm{2}}}{\rm{ + }}{z^{\rm{2}}}} \\ \beta & {\rm{ = ta}}{{\rm{n}}^{{\rm{ - 1}}}}\left({\frac{z}{{\sqrt {{x^{\rm{2}}}{\rm{ + }}{y^{\rm{2}}}} }}} \right)\\ \alpha & {\rm{ = ta}}{{\rm{n}}^{{\rm{ - 1}}}}\left({\frac{x}{y}} \right) \end{aligned} \right. $

(6)

式中，(x，y，z)为点的直角坐标，(r，α，β)为点的极坐标，β为仪器扫描方向与水平线的夹角，仰角为正，俯角为负，可换算成天顶角θ。

4 结果与分析

4.1 叶片点云

对高分辨率数码相机拍摄的2维照片进行监督分类后，枝干和叶片分别赋予不同的颜色，根据其与点云之间的位置校准信息，将照片颜色映射到对应的点云。林木点云着色后的全局图如图 6(a)所示，红色表示枝干，绿色表示叶片。从内部细节图 6(b)可以看出，绝大部分的叶片点云和枝干点云的颜色映射是正确的，但是仍然存在一些误差。可到看到少许叶片边缘，枝干边缘处存在一些颜色映射错分点，这与边缘处的噪点有一定的关系。在枝叶密集处，受光谱信息的影响，内部一些小枝干并未与叶片区分开。因此，根据点云颜色提取叶片点云，并无法完全去除枝干所贡献的点云，可选择通过人眼识别，进一步去除一些明显错误的点，得到最后的叶片点云结果如图 6(c)所示。

图 6 林木点云着色图

Fig. 6 Colorized trees point clouds

4.2 叶面积密度反演

4.2.1 叶倾角估算

从区域冠层中不同高度处随机选择200个叶片点云，利用特征值法估计平面法线，计算叶倾角。估算的叶倾角角度均分布在[0°，90°]，将[0°，90°]以 10°为间隔分为9个区间，统计各个区间的角度概率，得到叶倾角角度概率分布图，如图 7所示。

图 7 叶倾角概率分布图

Fig. 7 Distribution of leaf inclination anglesprobability

通过概率分布图可看出该玉兰树林叶倾角分布均匀，在给定的每个角度区间范围都有分布，集中分布在小于60°的区间内，平均叶倾角角度为35.5°。

4.2.2 接触频率和天顶角估算

通过凸包算法确定出冠层边界轮廓，接触频率可定义为植被冠层内激光光束被叶片挡住的概率，即体元体素值为1的个数和冠层内总体素个数的比值。经激光点直角坐标和极坐标的转换后发现天顶角范围为55°—89°，分布范围较大，这是因为区域林木高大，地基LiDAR进行扫描时，距离较短。如果直接利用整个冠层的天顶角平均值来计算叶倾角与激光光束的校正因子cos(θ)/G(θ)将产生较大误差，因而选择通过对冠层进行水平分层，利用每个水平层的平均天顶角计算校正因子。表 2为计算得到的不同水平厚度层的平均天顶角、接触频率和校正因子。

表 2 接触频率和校正因子
Table 2 The contact frequency and correction factor

下载CSV

	平均天顶角/(°)	接触频率	校正因子
第1层	88.34	0.13	0.99
第2层	80.11	0.12	0.99
第3层	78.79	0.25	0.99
第4层	75.97	0.36	1.00
第5层	71.63	0.39	1.05
第6层	68.70	0.42	1.05
第7层	65.03	0.47	1.05
第8层	63.16	0.52	1.11
第9层	60.62	0.48	1.11
第10层	59.24	0.32	1.11
第11层	59.39	0.15	1.11
第12层	59.99	0.14	1.11
第13层	58.30	0.10	1.11

4.2.3 叶面积密度反演

根据激光点云点与点之间的间距及区域树木的分布范围，考虑以0.1×0.1×0.1 m为体元大小构建区域林木冠层3维体元模型，并给体元赋予属性。以0.5 m为垂直间隔，将区域冠层划分成若干水平厚度层，分别利用VCP模型计算不同水平厚度层的叶面积密度，其垂直分布廓线如图 8所示。

图 8 叶面积密度廓线图

Fig. 8 LAD profile

从图 8可以看出，叶面积密度廓线与区域树叶片垂直分布基本一致，每个水平层内枝叶越茂盛该层叶面积密度值LAD就越大。在冠层中下部分随着高度的增加，叶面积密度增加，在4 m高度处，达到最大为1 m²/m³；之后随着高度的增加，叶面积密度逐渐降低。

4.3 叶面积密度实验验证

叶面积密度验证方法有直接和间接两种方法。直接方法是通过分层采集叶片，然后测量每层叶片单叶面积，这种方法工作量大且具有破坏性；对于区域范围的高大树木，实施该方法非常困难。间接方法是利用叶面积指数仪器测量树林冠层的叶面积指数，然后与叶面积密度计算的累积叶面积指数进行对比验证，这种方法容易实现但无法反映每个水平厚度层叶面积密度的反演精度。

本文利用LAI-2200冠层分析仪多次测量的区域平均叶面积指数3.16 m²/m²进行对比验证。LAI-2200冠层分析仪通过若干同心环传感器接收冠层上方和下方的光辐射并计算其比例，来推算光线透过冠层时被削弱的程度，从而得出冠层空隙度并计算LAI(梁顺林等，2013)。此计算结果是基于叶片在冠层内的空间分布是随机的，不考虑集聚效应的假设理论，计算结果为包括树干在内的总的叶面积指数，即有效叶面积指数，而非实际叶面积指数。集聚效应对实际叶面积指数的影响很大，即使相同的有效叶面积指数，随集聚指数的变化，实际叶面积指数变化也很大，一般情况下，实际叶面积指数值较有效叶面积指数值偏大。因此，实际叶面积指数较有效叶面积指数的稳定性差，不易测量；本文测量的有效叶面积指数只能作为一个参考数据进行对比验证。

利用每个水平厚度层的LAD计算的区域累积真实叶面积指数为3.20 m²/m²，比观测的有效叶面积指数值偏高0.04 m²/m²，相对误差为1.26%。考虑到LAI观测的叶面积指数有枝叶的影响，因此本文计算的LAD与实际相比有稍许偏高。

本文采用的VCP模型无需对冠层组分空间分布、大小和形状做任何假设，利用地基LiDAR系统发射的激光束对冠层进行充分采样后，通过去除树干等非叶片组分后计算其与叶片的接触频率，从而求得LAD。由于去除了非叶片组分的影响，通过分层处理，地基LiDAR又能较准确地发现那些被忽视掉的叶片，因而该方法计算出的累积叶面积指数与实际叶面积指数更接近。

5 结论

本文利用地基LiDAR设备ScanStationC10获取了阔叶林(玉兰树)点云数据，然后通过彩色照片监督分类快速提取叶片点云、构建体元模型和基于VCP进行了区域林木叶面积密度反演。本文的主要结论如下：

(1) 将LiDAR系统获取的彩色照片分类后的色彩信息映射到点云上，根据分类结果能较满意地分离出叶片点云，降低了非光合作用组分对叶面积密度反演的影响。

(2) 利用高分辨率的叶片点云数据进行叶倾角的估算，得到了区域样地真实的叶倾角概率密度分布，平均叶倾角为35.5°。分层估算的平均天顶角更适合计算叶倾角与激光光束方向的校正因子。

(3) 体元化的叶片点云数据更容易确定冠层边界和统计接触频率，得到的叶面积密度曲线变化走势与区域林木冠层叶片垂直分布相吻合，根据叶面积密度廓线计算得到的累积叶面积指数非常接近于LAI-2200实测的有效叶面积指数，仅偏高0.04 m²/m²。

本文基于LiDAR从多个位置获取的激光点云能够最大限度地减少遮挡效应这一假设进行研究。但在实际中，地基LiDAR即使从多站扫描目标，还是存在遮挡，导致有些叶片的点云获取不完整，从而影响最后的反演结果，如何填补和修复这些不完整的部分，有待于进一步研究。同时，树林冠层LAD直接测量非常困难，本文仅使用LAI-2200冠层分析仪测量的LAI对LAD累积量进行了验证，无法反映LAD分层反演结果精度，如何进行树冠分层LAD精度验证也是下一步需要重点考虑的问题。

参考文献(References)

Abbas M A, Setan H, Majid Z, Chong A K, Idris K M and AspuriA. 2013. Calibration and accuracy assessment of Leica ScanStation C10 terrestrial laser scanner//Developments in Multidimensional Spatial Data Models:LectureNotes in Geoinformation and Cartography.Berlin Heidelberg:Springer:33-47[DOI:10.1007/978-3-642-36379-5_3]

Béland M, Widlowski J L, Fournier R A, Côté J F andVerstraete MM.2011.Estimating leaf area distribution in savanna trees from terrestrial LiDARmeasurements. Agricultural and Forest Meteorology, 151 (9) : 1252–1266 . [DOI:10.1016/j.agrformet.2011.05.004]

Graham R L.1972.An efficient algorith for determining the convex hull of a finite planar set. InformationProcessing Letters, 1 (4) : 132–133 . [DOI:10.1016/0020-0190(72)90045-2]

Guan Y L, Cheng X J, Shi G G.2008.A robust method for fitting a plane to point clouds. Journal of Tongji University(Natural Science), 36 (7) : 981–984 .

( 官云兰, 程效军, 施贵刚. 2008. 一种稳健的点云数据平面拟合方法. 同济大学学报(自然科学版), 36 (7) : 981–984. )

Hosoi F, Omasa K.2006.Voxel-based 3-D modeling of individual trees for estimating leaf area density using high-resolution portable scanning lidar. IEEE Transactions on Geoscience and Remote Sensing, 44 (12) : 3610–3618 .

Hosoi F, Omasa K.2007.Factors contributing to accuracy in the estimation of the woody canopy leaf area density profile using 3D portable lidar imaging. Journal of Experimental Botany, 58 (12) : 3463–3473 . [DOI:10.1093/jxb/erm203]

Hosoi F, Omasa K.2011.Estimation of leaf area density profiles of Japanese maple and camellia woody canopies using portable scanning lidars. Eco-Engineering, 23 (4) : 105–109 . [DOI:10.11450/seitaikogaku.23.105]

Li W, Niu Z, Gao S, Qin Y C.2013.Analyzing and retrieving structural information of Piceacrassifolia based on airborne light detection and ranging data. Journal of Remote Sensing, 17 (6) : 1612–1626 .

( 李旺, 牛铮, 高帅, 覃驭楚. 2013. 机载激光雷达数据分析与反演青海云杉林结构信息. 遥感学报, 17 (6) : 1612–1626. [DOI:10.11834/jrs.20132330] )

Liang S L, Li X W, Wang J D. Quantitative Remote Sensing:Ideas and Algorithms. Beijing: Science Press 2013 .

( 梁顺林, 李小文, 王锦地. 2013. 定量遥感:理念与算法. 北京: 科学出版社 . )

Liu L X, Pang Y, Li Z Y, Xu G C, Li D, Zheng G.2014.Retrieving structural parameters of individual tree through terrestrial laser scanning data. Journal of Remote Sensing, 18 (2) : 365–377 .

( 刘鲁霞, 庞勇, 李增元, 徐光彩, 李丹, 郑光. 2014. 用地基激光雷达提取单木结构参数——以白皮松为例. 遥感学报, 18 (2) : 365–377. [DOI:10.11834/jrs.20143091] )

Pang Y, Li Z Y, Chen E X, Sun G Q.2005.Lidar remote sensing technology and its application in forestry. ScientiaSilvaeSinicae, 41 (3) : 129–136 .

( 庞勇, 李增元, 陈尔学, 孙国清. 2005. 激光雷达技术及其在林业上的应用. 林业科学, 41 (3) : 129–136. )

Rusu R B, MartonZ C, Blodow N, Dolha M, Beetz M.2008.Towards 3D point cloud based object maps for household environments. Robotics and Autonomous Systems, 56 (11) : 927–941 . [DOI:10.1016/j.robot.2008.08.005]

Swatantran A, Dubayah R, Roberts D, Hofton M, Blair J B.2011.Mapping biomass and stress in the Sierra Nevada using lidar and hyperspectral data fusion. Remote Sensing of Environment, 115 (11) : 2917–2930 . [DOI:10.1016/j.rse.2010.08.027]

Takeda T, Oguma H, Sano T, Yone Y, Fujinuma Y. Estimating the plant area density of a Japanese larch(LarixkaempferiSarg. )plantation using a ground-based laser scanner. Agricultural and Forest Meteorology. 2008 : 428 -438.

Treuhaft R N, Asner G P, Law B E, Van Tuyl S. Forest leaf area density profiles from the quantitative fusion of radar and hyperspectral data. Journal of Geophysical Research:Atmospheres(1984-2012). 2002 : ACL-7-1 -ACL 7-13.

Wang J Y, Wang W P, Tu C H, Yang C L. Computational Geometry and its Applications. Beijing: Science Press 2011 .

( 汪嘉业, 王文平, 屠长河, 杨承磊. 2011. 计算几何及应用. 北京: 科学出版社 . )

Weiss M, Baret F, Smith G J, Jonckheere I, Coppin P.2004.Review of methods for in situ leaf area index(LAI)determination:Part Ⅱ. Estimation of LAI, errors and sampling. Agricultural and Forest Meteorology, 121 (1/2) : 37–53 . [DOI:10.1016/j.agrformet.2003.08.001]

Zhao J, Li J, Liu Q H.2013.Review of forest vertical structure parameter inversion based on remote sensing technology. Journal of Remote Sensing, 17 (4) : 697–716 .

( 赵静, 李静, 柳钦火. 2013. 森林垂直结构参数遥感反演综述. 遥感学报, 17 (4) : 697–716. [DOI:10.11834/jrs.20132183] )

Zheng Gand Moskal L M.2012.Computational-geometry-based retrieval of effective leaf area index using terrestrial laser scanning. IEEE Transactions on Geoscience and Remote Sensing, 50 (10) : 3958–3969 . [DOI:10.1109/TGRS.2012.2187907]

Zhu Z and Zhang W M. 2011.Estimating the LAI of a single tree from terrestrial laser scanner data//Proceedings of the2011 International Conference on Remote Sensing, Environment and Transportation Engineering(RSETE). Nanjing, China:IEEE:4947-4951[DOI:10.1109/RSETE.2011.5965422]

摘要

关键词