杨晓慧 吴金卓 刘浩然 钟 浩 林文树

（1. 东北林业大学机电工程学院 哈尔滨 150040；2. 东北林业大学土木与交通学院 哈尔滨 150040）

作为森林生态系统的重要组成部分，人工林不仅改善了我国森林资源短缺的问题，而且也为国民经济建设提供了大量木材和林副产品，在发挥经济效益的同时对维持生态系统平衡、减缓全球气候变暖等起着重要作用（银彬吾等，2019）。 冠层是森林与外界环境相互作用最直接和最活跃的界面层，是植物进行呼吸作用、蒸腾作用和光合作用等生理过程的主要场所，本身承载了森林生物多样性的主体部分。郁闭度是森林冠层的主要参数，是反映森林结构和森林环境的重要因子，在水土流失、水源涵养、林分质量评价、森林景观建设等方面有着广泛应用。

林分郁闭度的传统测定方法主要有目测法、树冠投影法、样点法、鱼眼照片法和光学遥感法等，其中，地面实测方法获取的数据通常只在小范围内具有代表性且测量过程耗时耗力（孙钊等，2020），光学遥感因穿透性低一定程度上影响森林垂直结构参数估测。激光雷达（light detection and ranging，LiDAR）是一种新兴的主动遥感技术，能够有效穿透森林，在多时空尺度上获取森林生态系统高分辨率的三维地形、植被结构参数、叶面积指数等参数。Lefsky 等（2002）提出的地基激光雷达技术（terrestrial laser scanning，TLS）可用于测量单木三维结构和植被冠层结构；机载激光雷达技术（airborne laser scanning，ALS）的出现使得提取高精度森林结构参数成为可能（瞿帅等，2018）；星载激光雷达覆盖范围广，可用于大区域林业资源调查研究。无人机激光雷达相比上述3 种平台具有灵活性高、成本低、效率高、受地形影响小、影像重叠率大、分辨率高、姿态角大、相幅小、不受云雾干扰等特点（王娟等，2020；Salamíet al., 2014），近年来在森林结构参数提取中得到越来越多的应用（Davieset al., 2014）。

随着遥感技术的不断发展，国内外已有利用激光雷达数据进行林分郁闭度估测的研究报道。李丹（2012）基于配准后的多站地基激光雷达数据，设定不同高度阈值提取乔木样地郁闭度，估测精度达83.91%。赵勋等（2020）以机载 LiDAR 点云数据为基础，采用首次回波归一化后 2 m 以上点云密度与所有点云密度之比估测郁闭度，估测值与实测值之间的R2达0.478。Arumäe 等（2018）利用机载激光扫描（ALS）数据，设定不同高度阈值和回波计算郁闭度，并与基于数字半球照片（digital hemispherical photographs，DHP）和Crown模型（CCRCrown）的估测值进行了比较。García 等（2012）采用地球科学测高系统（geoscience laser altimetry system，GLAS）数据中的冠层基高度（canopy base height，CBH）估算林分郁闭度，精度达89%。吴项乾（2018）结合单木分割方法，先由UAV-LiDAR 点云生成冠层高度模型（canopy height model，CHM），再将基于点云提取的特征变量逐步引入郁闭度估测模型，R2max达0.78。汪霖等（2019）以地面实测数据、无人机高分影像、激光雷达点云数据为主要信息源，采用局部最大值法和种子点分割法提取林分郁闭度并对其精度进行了检验。Aicardi 等（2017）将TLS 点云与UAV影像点云数据结合，建立了一个具有高分辨率和高精度的森林数据集提取高郁闭度森林区域关键参数。

综上可知，目前林分郁闭度估测大多基于地基或机载雷达数据进行，关于东北地区人工林林分郁闭度的研究较少，特别是林分高度、强度、冠层特征变量对针、阔叶林郁闭度反演过程贡献率方面的研究鲜见报道。鉴于此，本研究选取我国东北地区具有代表性的兴安落叶松（Larix gmelinii）、 樟子松（Pinus sylvestrisvar.mongolica）、 黑 皮 油 松 （Pinus tabuliformis）、 蒙古栎（Quercus mongolica）、 白桦（Betula platyphylla）为研究对象，基于UAV-LiDAR 获取的高密度点云数据，结合植物冠层分析仪获取的40 个样本郁闭度实测值，采用30 个样本应用统计模型法进行针叶林和阔叶林林分郁闭度建模，采用10个样本验证其精度，比较影响针、阔叶林郁闭度的雷达特征变量并进行郁闭度制图，以期为快速准确估测人工林林分郁闭度提供基础数据和技术参考。

1 研究区概况

东北林业大学城市林业示范基地（127°35′—127°39′E，45°42′—45°44′N，海拔 136~140 m）地处黑龙江省哈尔滨市中心城区，占地面积 44 hm2，属中温带大陆性季风气候，年降雨量400~800 mm，年均气温3.5 ℃（刘浩然等，2021）。基地内栽种有兴安落叶松、樟子松、黑皮油松、胡桃楸（Juglans mandshurica）、蒙古栎、水曲柳（Fraxinus mandshurica）、白桦等针叶和阔叶树种，目前已形成小面积纯林块状混交。基地地形坡度平缓，林木较为密集，郁闭度较高（范伟伟等，2020）。研究区地理位置及样地空间分布见图1。

图1 研究区地理位置及样地空间分布Fig. 1 Geographical location of the study area and spatial distribution of the sample plots

2 研究方法

2.1 数据选择与预处理

2.1.1 地面数据采集 2021 年5 月，在研究区选取4块针叶林样地，包括2 块兴安落叶松样地、1 块黑皮油松样地、1 块樟子松样地；2 块阔叶林样地，包括1 块白桦样地、1 块蒙古栎样地（表1）。采用WinSCANOPY 2010a 植物冠层分析仪获取林木冠层半球图像，得到各样地实测郁闭度。具体做法是在针叶林和阔叶林样地内分别随机建立20 个半径为10 m 的圆形样方，利用 GPS 记录各样方地理位置，坐标参考WGS-84。在光线均匀、避免直射的天空条件下于各样方内随机选取3 个拍摄点，通过WinSCANOPY 2010a 植物冠层分析仪配有三星NV3 型号相机的180°鱼眼镜头拍摄林分冠层获取其半球图像，鱼眼镜头竖直向上并与水平面垂直。为消除灌木和杂草的影响，鱼眼镜头离地高度设为1.5 m。

表1 研究样地内5 种人工林概况Tab. 1 Basic information of five kinds of plantation in the sample plots

2.1.2 激光雷达数据采集 2021 年5 月，使用四轴大疆经纬M300 RTK 无人机搭载禅思L1 激光扫描仪获取激光点云数据。无人机载激光雷达数据获取参数如下：飞行高度100 m，飞行速度5 m·s-1，飞行间隔60 m，旁向重叠50%，主航线角度25°，平均点密度约197 pts·m-2，采样频率160 kHz，采取实时真彩色点云上色模式。

2.1.3 数据预处理 首先，对LiDAR 原始点云数据进行去噪、裁剪处理，通过CloudCompare 软件中的布料模拟滤波（cloth simulation filter, CSF）方法获得地面文件，并将分类的地面点插值以生成数字高程模型（digital elevation model，DEM）、数字表面模型（digital surface model，DSM），利用DEM 进行高度（z）归一化获得高度归一化点云数据（艾萨迪拉·玉苏甫等，2020；尤号田等，2020）。然后，应用Matlab 平台构建冠层体积模型（canopy volume model，CVM）计算冠层体积分布，采用PCM 点云数据处理软件中基于CHM 的方法进行单木分割，分割后的点云导入CloudCompare 软件中进行树冠点云提取，以计算冠层密度、冠层回波的偏度和峭度等。对采集的树冠球形真彩色图像，将其转换为灰度图像并绘制灰度图像直方图，结合树干剔除效果，调整分割阈值以确定合适阈值点（图2），取3个拍摄点的平均值作为该样方的郁闭度实测值。各样方郁闭度测量值见表2。

表2 各样方郁闭度测量值Tab. 2 Measured value of canopy closure of each quadrate

图2 图像二值化Fig. 2 Image binarization

2.2 技术流程

首先基于点云三维坐标获得冠层高度模型CHM和树冠点云，由CHM、回波能量值、树冠点云计算高度、强度、冠层3 类特征变量，然后采用主成分分析法和逐步回归法将处理后的特征变量与植物冠层分析仪获取的郁闭度进行建模，最后评价模型精度和适应性并进行郁闭度制图（图3）。

图3 技术流程Fig. 3 Technical flow chart

2.2.1 LiDAR 特征变量提取与标准化 离散点云归一化后，计算样地尺度的激光雷达特征变量，包括树木高度和强度特征变量；应用Matlab 平台构建冠层体积模型，计算冠层体积分布；基于单木分割提取样地冠层信息，计算冠层密度、冠层回波的偏度和峭度等。冠层体积分布（canopy volume distribution, CVD）是描述三维点云形态的重要指标，通过对每块样地的点云集合进行三维分析，获取基于三维体的冠层体积分布点云变量，步骤如下：按照三维坐标，将点云分为若干个5 m×5 m×0.5 m 的小三维体，若小三维体含有点云，则该小三维体属于“非空三维体”，反之为“空三维体”；将位于冠层以上的“空三维体”定义为“孔隙体”（Open），反之定义为“封闭体”（Close）；将位于上65%区域的“非空三维体”定义为“透光体”（Euphotic），反之定义为“遮光体”（Oligophotic）；汇总所有三维体，形成冠层体积分布图并汇总计算CVD 内4 种类别三维体分别占有的体积比例（以某针叶样本为例，见图4）。考虑到地面存在干扰点，滤除0.3 m 以内的地面点（Voseet al., 1995）。本研究计算高度、强度、冠层3 类特征变量（表3） 。利用IBM SPSS Statistics 软件对提取完成的3 类特征变量进行标准化处理。

表3 LiDAR 特征变量及其描述①Tab. 3 LiDAR characteristic variables and descriptions

图4 冠层体积模型Fig. 4 Canopy volume model

2.2.2 主成分分析 本研究共计算30 组×26 个特征变量，在每组特征变量内部，变量之间可能存在相关性，若对每个指标分别进行分析，则不能完全利用数据中的信息，因此采用IBM SPSS Statistics 软件，在保证主成分分析有效的前提下（KMO 检验系数>0.5，巴特利球体检验P<0.05），分组对高度、强度、冠层特征变量进行分析，减少分析指标的同时，尽可能多地保留原始变量信息且使变量彼此不相关，作为新的综合指标。主成分选取上，每组保留特征值大于1 的主成分。主成分载荷矩阵Ui与因子载荷矩阵Ai和特征值λi之间的关系为：

将Ui与每组特征变量的标准化值Zxi相乘即得主成分F的表达式：

2.2.3 逐步回归分析 针对针叶林和阔叶林，利用提取的多组LiDAR 特征变量结合多元逐步回归方法构建郁闭度估测模型，形式如下（吴项乾等，2020）：

式中：y为因变量；q为自变量；b为常数；ε 为误差。

以地面实测林分郁闭度为因变量，以高度、强度和冠层特征变量为自变量，应用逐步进入法（stepwise）将激光雷达特征变量按照与实测林分郁闭度相关性由高到低的顺序引入方程，在保证每一步引入变量检验值F达到显著水平（P<0.05）的前提下，筛选出高相关性变量，剔除低相关性变量，以得到最优模型。

2.2.4 回归模型精度评价 选取精度和真实程度高的变量组合作为反演模型。本研究采用反映自变量与因变量之间相关性的判定系数R2、调整判定系数AdjR2、均方根误差RMSE（root mean square error）和相对均方根误差rRMSE（relative root mean square error）评价回归模型精度（洪奕丰等，2019）。R2越高，则二者之间相关性越高；但样本容量一定时，增加解释变量必定使自由度减少，AdjR2可在模型的复杂程度和衡量模型的优良程度上取一个平衡，使模型趋于简单；RMSE 反映实测值与预测值之间标准误差大小，与评价对象本身的数值关系很大，RMSE 越小，则模型预测效果越好；rRMSE 为 RMSE 与预测值算数平均值的比值，反映模型总体预测精度，rRMSE 越小，表明模型预测精度越高（刘浩等，2018）。各评价指标计算公式如下：

式中：n为样地数量；X为鱼眼镜头测得的林分郁闭度；为实测郁闭度的平均值；Xˆ为模型估测的林分郁闭度；k为自变量个数。

3 结果与分析

3.1 主成分分析

本研究提取高度、强度和冠层特征变量并将3 类特征变量标准化。经主成分分析，针叶林样本的高度特征变量得到3 个主成分F1c、F2c、F3c，强度特征变量得到4 个主成分F4c、F5c、F6c、F7c，冠层特征变量得到2 个主成分F8c、F9c；阔叶林样本的高度、强度和冠层特征变量分别得到3、3 和2 个主成分，分别记为F1b、F2b…F8b。由式（2）得出各组特征变量的主成分得分及排名，以针叶林的高度特征变量得到的3 个主成分为例，X、Z轴分别代表第一和第二主成分，Y轴代表第三主成分，颜色不同的3 类点在空间中的位置表示兴安落叶松、黑皮油松、樟子松3 类针叶树种的第一、第二和第三主成分得分，结果见图5。

图5 针叶林高度特征变量主成分得分Fig. 5 Principal component score of height characteristic variables for coniferous forests

3.2 逐步回归分析

经逐步回归法分析，得到针叶林（式8）、阔叶林（式9）林分郁闭度估测的回归方程：

回归方程的精度如表4 所示。在针叶林郁闭度估测中，回归方程的AdjR2=0.722，模型拟合优度较高，显著水平P=0.000 6<0.05，模型显著性较高，变量F8c前的系数绝对值最大，表明冠层特征变量在针叶林郁闭度预测中贡献最大；在阔叶林郁闭度估测中，回归模型的AdjR2=0.725，模型拟合优度高，P=0.000 6<0.05，模型显著，强度特征变量对阔叶林郁闭度预测精度影响最显著。

表4 回归方程的精度Tab. 4 Accuracy of regression equation

3.3 精度评价

将提取的特征变量进行主成分分析，代入回归方程可得郁闭度估测值。针、阔叶林各选择5 个样方作为验证样本，通过回归方程预测的林分郁闭度与WinSCANOPY 2010a 植物冠层分析仪获得的郁闭度的交叉验证结果见图6。针叶林与地面实测郁闭度建立的预测模型精度为R2=0.72，阔叶林与地面实测郁闭度建立的预测模型精度为R2=0.77，筛选10 个检验点的郁闭度与实测郁闭度显示出较高相关性（r=0.859）。通过反演模型计算得到的林分郁闭度与实测郁闭度较均匀地分布在拟合线两侧，表明模型可用来反演研究区典型针、阔叶林林分郁闭度，反演结果精度较高。

图6 林分郁闭度预测模型交叉验证结果Fig. 6 Cross validation results of forest canopy closure prediction model

利用ArcGIS 的反距离权重插值法将激光雷达估测郁闭度与地面实测郁闭度可视化。由图7 可知，郁闭度范围在0.81~0.87 之间，样地郁闭度呈东北西南高、西北东南低的趋势，与每块矩形样地的走向基本一致。无人机激光雷达数据在预测林分郁闭度方面优势明显，多变量模型能充分发掘森冠层信息，结合多组UAV-LiDAR 特征变量估测林分郁闭度能够取得较好效果。

图7 实测与估测森林郁闭度Fig. 7 Measurement and estimation of forest canopy closure

4 讨论

本研究基于无人机激光雷达点云数据和WinSCANOPY 2010a 植物冠层分析仪样地实测数据，采用逐步回归法构建人工林林分郁闭度估测模型，基本满足郁闭度反演需要，可实现基于无人机激光雷达点云数据的林分郁闭度快速估测。谷金英等（2014）提取遥感光谱因子和地形因子，基于多元回归分析法对林地郁闭度进行反演，精度达81.6%；吴飏等（2012）利用Spot5 数据结合96 个纹理特征估测林分郁闭度，主成分分析与逐步回归分析结果显示，近红外波段的光谱和纹理特征与郁闭度相关性最显著，加入纹理特征后估测精度可达84.32%；李擎等（2019）基于GF-2号遥感影像，采用主成分法进行天山云杉（Picea schrenkiana）林郁闭度估测研究，以纹理特征+光谱信息+地形因子为自变量构建的估测模型拟合度为0.823。本研究在点云高度、冠层特征变量基础上引入回波强度特征变量，针叶林回归模型R2为0.782，阔叶林回归模型R2为0.784。这是因为针叶树冠形为圆锥形，无人机激光雷达在发射激光束的过程中易错失冠顶，而阔叶树冠顶呈椭球形，截获点云概率较高，相对于针叶树提取冠顶的概率更大，高度特征变量的估测精度更高；此外，激光点云穿透针叶林冠层以及通过冠层间的空隙进入林内形成地面回波点的概率更高，在植被点与全部回波点的比值运算中更易造成低估现象。因此，阔叶林郁闭度反演精度高于针叶林。由于本研究将对郁闭度影响较大的CC2m归入冠层特征变量，故估测方程中冠层特征变量超过高度特征变量，成为对郁闭度预测精度影响更为显著的因素。

与基于机载雷达点云数据进行郁闭度估测的研究相比，本研究数据获取方式更灵活，时效性更强，精度更高，如张瑞英等（2016）结合ALS 点云数据与Landsat ETM+数据，采用Cubist 模型估测温带森林郁闭度，R2为0.722；与基于无人机激光雷达点云提取的特征变量构建统计模型进行郁闭度估测的研究相比，本研究估测精度有一定提高，如吴项乾（2018）采用雷达高度特征变量和冠层特征变量建模估测阔叶林郁闭度的R2为0.69（rRMSE=6.5%），而当通过高度特征变量和密度特征变量联合进行郁闭度估测时R2为0.78（rRMSE=5.9%）；与使用某一高度以上点云所占比例直接代表郁闭度的研究相比，本研究构建模型在不同森林类型郁闭度估测中稳定性更好，估测精度也有所提高，如穆喜云等（2015）采用高于2 m 的植被点回波与样地内全部点云回波的比值代表郁闭度，全部60 块样地的R2为0.721，阔叶林估测精度为56.62%。因此，本研究基于无人机激光雷达数据提取高度、冠层和回波强度特征变量构建的人工林林分郁闭度估测模型具有较高精度，可满足森林资源调查中的精度要求。

传统人工调查方法实现森林资源动态监测费时费力，针对样地或区域等尺度的林木调查，利用无人机LiDAR 扫描是一种灵活、高效和低成本的地面数据获取方法，能够满足快速的林分郁闭度估测需求。本研究提出的郁闭度估测方法只需获取目标区域较高密度和较强穿透力的LiDAR 点云及少量地面实测验证数据，提取相关特征变量就能建立适配研究区森林特征的反演模型，进而得到精度较高的林分郁闭度估测值。本研究建立的人工针/阔叶林林分郁闭度模型，预估精度能达到森林资源调查相关技术规定要求，可在实践中推广应用。

需要说明的是，尽管无人机激光雷达在森林结构参数提取和估测中具有诸多优势，但本研究还存在一些不足之处：因飞行数据无法同时保证点云的高密度性和强穿透性，本研究仅选取林分结构单一的人工林为研究对象，未对林间情况更复杂的天然林进行郁闭度估测；由于鱼眼镜头使用条件限制，所选取的样本点分布不够均匀；另外，由于UAV-LiDAR 技术处于发展阶段，历史数据缺乏，现阶段对比分析不同时期的森林监测数据难度较大。在未来的针对天然林郁闭度估测中，有待选取大量具有代表性的天然林样本，开发更精确、更简洁的算法增强模型的鲁棒性，如采用多元线性回归（multiple linear regression，MLR）、偏最小二乘回归法和随机森林（random forest，RF）等多种方法构建反演模型，讨论评价不同方法对估测精度的影响；另外，联合GF、Landsat-TM、ETM+数据，对激光雷达计算的森林参数进行准确性和可靠性验证，可为人工林精细化管理和森林资源监测提供技术支持。

5 结论

1） 本研究构建的基于无人机激光雷达数据的人工林郁闭度估测模型具有较高精度，且人工阔叶林郁闭度估测模型的判定系数、调整判定系数和均方根误差等精度指标均优于人工针叶林。

2） 引入多类特征变量估测林分郁闭度效果较好，采用逐步回归法建立统计模型反演的郁闭度与鱼眼照片提取的郁闭度之间有良好的关系。影响针叶林和阔叶林林分郁闭度反演的特征变量不同，冠层特征变量在针叶林郁闭度反演中所占权重更大。

3） 利用回归预测模型和空间插值法得到整个样地的郁闭度，检验点郁闭度与实测郁闭度显示出较高相关性（r=0.859），表明利用该预测模型和反距离权重插值法反演制图预测林分郁闭度是可行的且优势明显。本研究对东北地区常见的针叶林兴安落叶松、黑皮油松、樟子松，阔叶林白桦、蒙古栎进行林分郁闭度预测对比研究，取得较好预测效果，表明UAVLiDAR 点云在估测东北地区针叶林和阔叶林林分郁闭度上具有较好潜力。