森林既是碳汇又是碳源，在大气碳含量方面起着重要作用^[1]。森林生态系统在气候变化、水循环、生物多样性中发挥着重要作用，而森林地上生物量、碳储量等参数与森林生态系统结构和功能密切相关。碳储量是有机与无机环境间能量循环的衡量指标，且森林地上生物量与碳储量是探测气候变化重要指标^[2]。森林碳储量在森林生态系统生产、代谢、运转过程中起着重要作用，是研究森林固碳潜力、森林生产力、环境与气候变化^[1]的重要基础。此外，估测碳储量不仅有助于评估一些景观的生物能源，促进生态平衡与环境协调发展，并有助于监测森林资源的可持续性，可协助相关政府部门对森林可持续经营进行规划^[3]、制定不同的造林政策和有效实施森林管理战略。因此，以温带人工林为研究对象估测森林碳储量具有重要意义。

林分空间结构和森林参数之间存在着十分紧密的关系，而激光雷达对森林具有一定穿透力，具有精确描绘森林冠层三维结构的能力，该技术已成功应用于提取树高、冠幅、蓄积和生物量^[3-9]等森林参数以及广泛用于森林估测与制图^[10]。传统的森林地上生物量和碳储量以野外实地调查数据进行计算，这种地面调查方法可信度高，但较为繁琐，并对森林造成一定破坏性。传统调查在大区域范围难以实现对森林碳储量的动态监测，而激光雷达在国家尺度、州尺度及区域尺度^[11-13]进行森林资源监测已获得广泛运用，故许多学者利用无人机激光雷达估测森林参数。

激光雷达系统估算森林生物量与碳储量等森林参数方面有较多研究。早在前苏联时就有学者用激光雷达数据估测树高并获得较好的估测精度^[14]。高度变量已被证明与林分参数的估算密切相关^{[4, 15]}；Parker等^[16]运用激光雷达数据进行林分高度、郁闭度森林参数的估测。Lefsky等^[17]提取点云特征变量：最大值、中位数、平均值和二次平均值，用这4个高度变量建立地上生物量回归模型，Næsset等^[18]考虑百分位高度与郁闭度特征变量进行建模来估测挪威北部森林生物量。Fischer等^[19]研究人员通过水平和垂直结构的特征变量来分析森林结构与生物量关系，结果表明树冠顶部高度是描述水平森林结构的最佳预测指标，而枝叶剖面的标准偏差是森林垂直异质性的最佳预测指标。密度变量作为林分结构的重要指标，用于描述样地内水平分布情况，其中郁闭度可以很好的反应出林分枝叶的分布程度，森林垂直剖面描述了植被垂直结构分布信息^[20-21]。以往的研究中，以温带人工林为研究对象考虑垂直结构变量的研究甚少。

激光雷达估测森林参数一般以单木法和面积法两种方式进行：单木法以点云数据生成数字冠层模型或者直接以点云数据为基础^[22-23]，使用单木分割算法进行分割，得到单木树高、冠幅、材积和胸径等参数。面积法以样地面积为计算单位，统计点云特征参数与森林参数为基础，如平均高、郁闭度和垂直结构变量等参数。单木法需要较高的点云密度和精确分割算法，大区域估测数据量庞大和计算时间长^[24]。近年来面积法用于各类森林参数估算已有大量模型，包括参数模型和非参数模型^[25-26]。参数模型为多元线性模型和非线性模型，非参数模型包括人工神经网络、支持向量机、深度学习和回归树等^[27-28]机器学习方法。其中一部分研究，非参数模型优于参数模型，而另一部分参数模型估测精度更高。从林分三维结构出发，参数模型更好解释^[29-30]，而众多研究中多元线性回归模型是最常用的森林参数估测模型^{[5, 31]}。也有研究表明，对于大多数森林参数（如蓄积量、生物量等）估测而言，非线性模型的精度优于线性模型。因此，本文利用两种模型估测人工林碳储量，探究不同模型对森林参数估测的差异。

本研究以内蒙古自治区赤峰市喀喇沁旗旺业甸实验林场落叶松（Larix gmelinii）和油松（Pinus tabuliformis）针叶林为研究对象。从林分冠层的三维结构分析出发，通过具有明确生物物理和林学解析意义的激光雷达变量的组合。基于无人机激光雷达数据和样地实测数据估测森林碳储量，在以激光雷达高度变量、密度变量的同时，考虑加入垂直结构变量进行建模，探讨加入垂直结构变量对模型精度的影响。本文通过比较分析多元线性模型与多元乘幂模型估算人工林碳储量，并对精度进行评价。从而更精准地估测森林资源（碳储量、生物量和蓄积量等）参数，量化森林碳储量有利于了解森林生态系统的固碳能力，为碳循环和森林可持续经营提供数据支撑。

1.   研究区概况与数据获取

1.1   研究区概况

研究区位于内蒙古自治区赤峰市喀喇沁旗旺业甸实验林场，试验区面积为566 km²，地理位置为 118°09′ ~ 118°30′E、41°21′ ~ 41°39′N，如图1所示。属于温带季风气候，冬季漫长寒冷，春季干旱多风，夏季雨热同期，秋季短促。研究区整体山势陡峭, 相对高差较大，平均海拔在1 300 ~ 1 900 m之间，坡度一般为15° ~ 35°。构成本区山地的基岩，主要是花岗岩、片麻岩等，土壤为山地草甸土、棕壤土和棕壤型粗骨土或石质土等。属于暖温带向寒带的过渡地区，冬季寒冷, 夏季温热。森林类型为暖温带针叶林和落叶阔叶林，优势树种主要为落叶松、油松和白桦（Betula platyphylla）蒙古栎（Quercus mongolica）和山杨（Populus davidiana）等。林下有少部分灌木和草本，以土庄绣线菊（Spiraea pubescens）、山玫瑰（Rosa davurica）、库页悬钩子（Rubus sachalinensis）、小红菊（Dendranthema chanetii）、紫鳞苔草（Carex angare）为主，植被类型为山地森林植被。

图  1  研究区地理位置

Figure  1.  Geographic location of the study area

下载: 全尺寸图片 幻灯片

1.2   无人机激光雷达数据采集和处理

本研究使用的无人机激光雷达数据于2017年9月在赤峰市喀喇沁旗旺业甸实验林场采集，以八旋翼无人机，搭载RIEGL VUX-1UAV激光扫描仪，飞行相对高度为300 m，速度为4.8 m/s，激光扫描角范围330°，脉冲重复频率380 kHz，最大有效测量速率500 000 点/s，点云密度约40点/m²。在旺业甸林场选取了3个样区，样区1位于南部，面积约2.6 km²；样区2位于西部，面积约2.1 km²；样区3位于东部，面积约0.6 km²；3个样区激光雷达数据的总覆盖面积约5.3 km²。样区1布设了样地11块，其中油松样地2块，落叶松样地9块；样区2布设了样地6块，其中油松样地3块，落叶松样地3块；样区3布设了样地9块，其中油松样地6块，落叶松样地3块。3个样区共布设了样地26块，其中油松样11块，落叶松样地15块；每个样地的尺寸为25 m × 25 m，取的点云格式为LAS标准格式。使用TerraSolid软件对激光点云进行分类处理得到空间分辨率为0.5 m的数字高程模型（digital elevation model，DEM），使用DEM对激光雷达点云数据进行归一化处理，得到归一化植被点云数据。

1.3   样地数据

本研究的地面调查数据是2017年9—10月根据不同树种、林龄选取了25 m × 25 m的样地26个，主要为油松和落叶松林。油松11个样地，落叶松15个样地，其中包含落叶松和白桦混交4块样地。样地中心点和角点坐标通过载波相位差分技术（real-time kinematic，RTK）定位，样地方向和边界长度通过RTK放样。在各个样地中，使用胸径尺进行每木检尺（DBH≥5 cm），测其胸径、树高、冠幅、枝下高和树冠垂直结构等单木因子。树高使用激光测高器测量，冠幅是使用皮尺在冠幅的两个主方向上测其投影长度。对于枯木标注出来不参与地上生物量的计算，根据实测的单木信息计算单木地上生物量，模型（1）、（2）和（3）所示^[32-34]，汇总得到样地地上生物量，在转换为单位面积地上生物量。生物量含碳率：油松0.521、落叶松0.521和白桦0.491，参数来源于国家发展与改革委员会2013年发布《碳汇造林方法学》。样地林分特征如表1所示。

表  1  样地林分特征

Table  1.  Stand characteristics of sample plots

林分特征 Stand characteristics	变化范围 Variation range	均值 Mean	标准差 Standard deviation
地上生物量/(t·hm−2) Aboveground biomass/(t·ha−1)	60.41 ~ 293.85	156.65	61.73
碳储量/(t·hm−2) Carbon storage/(t·ha−1)	31.47 ~ 153.10	81.48	32.18

下载: 导出CSV 

| 显示表格

注：Y_agb为油松生物量，L_agb为落叶松生物量，H_agb为白桦生物量，D为胸径，H为树高。

2.   研究方法

2.1   点云特征变量提取与筛选

基于归一化处理后的点云数据提取特征变量,点云特征参数变量描述见表2。根据森林参数的估测情况，本研究分别提取高度变量，密度变量，垂直结构变量。其中高度变量包括：点云平均高、高度百分位数、高度标准差和高度变动系数；密度变量包括：郁闭度和密度百分位数；垂直结构变量包括：叶面积密度均值、叶面积密度标准差和叶面积密度变动系数。研究表明森林地上生物量与蓄积量存在相关性^[35]，由于地上生物量是胸径和树高的函数，树冠与胸径间存在转换关系，即用郁闭度和平均高可以估测地上生物量。碳储量与生物量间以含碳率转换，在碳储量估测的模型中确定两个固定特征变量平均（H_m）、林分郁闭度（cc）后。根据皮尔逊相关系数分析LIDAR变量与预测森林参数变量之间的相关性，选出相关关系显著的变量，基于样地数据和激光雷达数据分别建立多元线性模型与多元非线性回归模型，通过两种建模方法的比较分析最终确定林分碳储量最优模型。初步筛选的3组变量，采用穷举法进行变量组合，x₁和x₂为选取的不同点云变量，分别得到线性模型和非线性模型16个。二参数模型1个，三参数模型3个和四参数模型12个。本研究使用SPSS26进行变量筛选，python3.8软件，用高斯牛顿迭代法进行拟合，并求解模型参数。

表  2  激光雷达点云特征变量

Table  2.  Characteristic variables of LiDAR point cloud

点云特征变量 Characteristic variable of point cloud	特征变量描述 Characteristic variable description
平均高 Mean height (Hm)	归一化高度的平均值 Average of normalized height
高度标准差 Height SD (Hs)	归一化高度的标准差 SD of the normalized height
高度变动系数 Height variation coefficient (Hv)	归一化高度的变异系数(标准差与平均数的比值) Variation coefficient of the normalized height (ratio of SD to the mean)
高度百分位数(hp10, hp15, …, hp95) Height percentile (hp10, hp15, …, hp95)	激光返回点的高度分布百分位数 Height distribution percentile of LiDAR return
郁闭度 Canopy closure (cc)	高于2 m的激光返回点所占的百分比 Percentage of LiDAR returns above 2 m to total returns
密度百分位数(dp10, dp15, …, dp95) Density percentile (dp10, dp15, …, dp95)	在各百分位高度等级以上的激光返回点在所有返回点中所占的百分比 Percentage of LiDAR returns above each percentile level to the total returns
叶面积密度均值 Mean leaf area density (Lm)	根据Beer-Lamber[20]法则计算其叶面积，计算平均值 Leaf area is calculated according to the Beer-Lamber[20] rule and the average value is calculated
叶面积密度标准差 SD of leaf area density (Ld)	叶面积密度的标准差 SD of leaf area density
叶面积密度变动系数 Variation coefficient of leaf area density (Lv)	叶面积密度的变异系数(标准差与平均值的比值) Variation coefficient of leaf area density (ratio of SD to the mean)

下载: 导出CSV 

| 显示表格

2.2   多元回归分析

2.2.1   回归模型

依据激光雷达变量优化与筛选得到的变量构建二变量、三变量和四变量的多元回归线性模型。以研究区地面实测所得到的林分地上生物量与碳储量数据作为因变量，激光雷达特征变量作为自变量，建立多元线性回归模型^[27]，模型结构式见公式（4 ~ 6）。非线性模型应用最广泛，其中相对生长模型最具有代表性，是所有模型中应用最普遍的一类模型。当前大多数学者普遍接受、符合生物学的地上生物量异速生长模型的基础上^[20-21]，异速生长方程可将点云特征变量与森林参数（如：地上生物量）相联系，本研究选择碳储量与点云特征变量的幂函数模型^[20-21]来估测森林参数，其表达式如（7 ~ 9）式所示。

式中：Y为碳储量（t/hm²），H_m为点云平均高，cc为郁闭度，a₀、a₁、a₂、a₃、a₄为模型待定回归系数， x₁为点云变量（L_m、L_d、L_v），x₂为点云变量（H_s、H_v、dp50、dp75）。

2.2.2   模型评价

采用留一法来验证模型的估测精度，确定系数R²和相对均方根误差rRMSE用于模型评价的指标。该方法是在模型已经确定的情况下，随机选取一个样地作为验证样地，而其他的样地数据进行建模，做n−1次验证。评价指标为n−1次随机试验结的值，筛选最优模型。公式如下（10 ~ 13）：

式中：${x}_{i}$为地面实测的林分特征数据，${\hat x_i}$为模型估测的林分特征数据，$\overline x$为${x}_{i}$平均值，n为样本数量，k为变量个数。

3.   结果与分析

3.1   模型拟合效果

通过决定系数与相对均方根误差筛选得到的不同特征变量的最优模型（见表3）。碳储量固定特征变量（H_m，cc）模型的拟合效果决定系数都高于0.70。证明了在估算碳储量森林参数中，郁闭度与平均高对模型拟合具有很大的贡献。就不同模型引入垂直结构特征变量而言，线性模型在未引入垂直结构特征变量时（${R_{\rm{adj}}^2}$ = 0.69，rRMSE = 21.39%），选入垂直结构特征变量（L_v）的模型（${R_{{\rm{adj}}}^2}$ = 0.67，rRMSE = 21.31%）。而非线性模型选入叶面积密度均值（ladmean）后模型拟合效果提高，由（${R_{\rm{adj}}^2}$ = 0.74，rRMSE = 19.39%）提高到（${R_{{\rm{adj}}}^2}$ = 0.82，rRMSE = 16.18%）。线性模型中四参数拟合效果最优（${R_{{\rm{adj}}}^2}$ = 0.88，rRMSE =12.90%），非线性模型中四参数拟合效果最优（${R_{{\rm{adj}}}^2}$ = 0.90，rRMSE =11.23%）。可见，点云特征变量的高度变量对模型的拟合效果也具有较大影响。

表  3  多元回归预测模型

Table  3.  Prediction models of multiple regression

模型 Model	碳储量预测模型 Prediction model of carbon storage	${R_{{\rm{adj}}}^2}$	rRMSE/%
(4_1)	$y=-40.16+9.79H_{\rm{m}}+2.89{\rm{cc}}$	0.69	21.39
(5_1)	$y=-66.53+10.03H_{\rm{m}}+9.95{\rm{cc}}+11.84L_{\rm{v}}$	0.67	21.31
(6_1)	$y=-135.98+9.78H_{\rm{m}}+4.84{\rm{cc}}+7.77H_{\rm{s}}+72.09L_{\rm{m}}$	0.88	12.90
(7_1)	$y=0.72{H_{\rm{m}}}^{1.88}{{\rm{cc}}}^{0.15}$	0.74	19.39
(8_1)	$y=0.67{H_{\rm{m}}}^{1.89}{{\rm{cc}}}^{0.14}{L_{\rm{m}}}^{0.86}$	0.82	16.18
(9_1)	$y=0.76{H_{\rm{m}}}^{1.98}{{\rm{cc}}}^{0.23}{H_{\rm{v}}}^{0.18}{L_{\rm{v}}}^{0.05}$	0.90	11.23

下载: 导出CSV 

| 显示表格

3.2   适应性检验

通过对比线性回归方法与非线性回归方法基于森林参数构建的估测模型，得出非线性回归模型拟合效果较线性回归模型拟合效果更好。其中非线性回归模型（9_1）检验效果最好（R²为0.86，rRMSE为9.92%），线性回归模型（6_1）（R²为0.85，rRMSE为12.38%）。线性与非线性模型R²相差0.01，而rRMSE相差2.46%。模型（5_1）与（4_1）相比精度降低，相对均方根误差增大。模型（7_1）与（8_1）相比精度提高，相对均方根误差减小。当选入参数相同的情况下，非线性回归模型拟合效果与稳定性明显优于线性回归模型。碳储量（图2）实测值与估测值的散点图分布趋势可以看出，森林参数的预估值大多分布于1∶1线的两侧，趋势与拟合结果相符，说明总体上模型预测效果较好。

图  2  碳储量参考值与预测值的散点关系

实线为1∶1验证线。The solid line is 1∶1 verification line.

Figure  2.  Scatter relationship between reference values and predicted values of carbon storage

下载: 全尺寸图片 幻灯片

4.   讨　　论

在本研究中，我们以面积为566 km²的内蒙古自治区赤峰市喀喇沁旗旺业甸实验林场为研究区，通过人工林的地面样地数据和无人机LiDAR数据，采用多元线性回归模型、多元乘幂模型，对地上碳储量进行估测，并对不同模型和变量间的估测影响进行了分析，就研究范围、样地量和研究方法而言，我们的研究结论对人工林的碳储量估测具有一定意义。

4.1   垂直结构变量分析

森林参数与森林冠层的三维结构密切相关，激光雷达高度变量和密度变量精确地刻画了森林冠层的三维结构，垂直结构变量（叶面积密度均值、标准差和变动系数）反映了亚优势木的分布，因此，由高度变量和密度变量的基础上，加入垂直结构变量，更准确地刻画林分冠层的三维结构。研究结果发现乘幂回归模型在高度变量和密度变量估测时的相对均方根误差为19.39%，而加入垂直结构变量后相对均方根误差为16.18%。本文中4个变量的线性模型与非线性模型相比于2个变量和3个变量的模型而言，其模型拟合效果与预测精度更好，表明加入垂直结构变量对估测人工林碳储量的重要性。这一结果与刘浩等^[8]研究结果一致，他们研究基于北亚热带沿海平原人工林研究结果表明利用冠层垂直结构剖面参数估测林分参数可以提高预估精度。出现这一结果一方面是森林在垂直方向上具有的分层现象，垂直结构变量有助于获取植被分层的变化^[9]。枝叶剖面的标准偏差是森林垂直异质性的最佳预测指标^[20]，垂直结构的特征变量可以用来分析森林结构与生物量关系。另一方面研究区落叶松和油松为不同龄林，激光雷达点云垂直结构特征变量能够较好的代表林分整体状态。故本文以温带针叶林和落叶阔叶林为对象，在点云高度变量与密度变量的基础上，加入垂直结构变量估测人工林碳储量。

4.2   模型分析

模型变量的优化与筛选通常采用主成分分析方法、逐步回归法、贝叶斯模型平均法^[36]等。主成分分析法在变量降维过程中，变量解释意义很模糊；逐步回归法、贝叶斯模型平均法虽然基于原始变量，但从林分结构描述而言，筛选的变量同样缺乏明确的生物物理意义和林学解析性。本文确定两个固定变量：H_m、cc，通过穷举法结合其他变量进行组合，考虑森林三维结构获取最优的模型。建模的方法有很多，估测森林参数的建模方法大致分为两种：一种是基于统计分析的参数方法，由激光雷达数据提取林分高度、密度等特征变量，构建线性或非线性回归模型估算森林参数；另一种是采用机器学习的非参数方法来构建森林参数估测模型，如：支持向量回归（SVR）、K-近邻算法（KNN）、随机森林（RF）、人工网络神经（BP）、多元自适应回归（MARS）、堆叠稀疏自动编码（SSAE）^{[5, 37]}等。多元线性回归模型建立了参数变量之间关系，但却可能忽略非线性的因果关系和交互效应，机器学习一定程度上使得模型拟合精度提高，但一系列的暗箱操作却未考虑模型机理性，不同的方法各有其优缺点。本研究使用无人机激光雷达估测碳储量结果与Jayathunga等^[38]一致，Jayathunga等使用了点云平均高、郁闭度和高度标准差来预测阔叶混交林碳储量，相对均方根误差为16.7%。相对均方根误差出现差异的原因可能是森林类型或者模型变量不一样。例如Takagi等^[39]、Puliti等^[40]和Bohlin等^[41]的结果得出相对均方根误差分别为23.3%、15.0%和25.5%。结果之间的细微差异可能是由于研究之间使用的数据源、数据采集参数、森林的结构复杂性和模型变量的不一样。本文从模型机理出发对线性模型与非线性模型方法进行比较，无论是2变量、3变量、还是4变量模型，其非线性模型都优于线性模型。从学科角度上来看，在复杂的森林生态系统中，多是非线性关系，故估测碳储量所建立的非线性模型表现更好。

众多学者致力于机载激光雷达估测森林碳储量的研究，精确估测森林碳储量不仅可以推动碳循环和碳汇研究工作的进展，而且可为森林资源监测和经营管理提供决策服务^[3]。然而关于北方人工林的相关研究较少，本文以落叶松和油松针叶人工林为研究对象，采用非线性模型与线性模型估测碳储量，揭示垂直结构变量对估测碳储量的影响，但由于研究样地数量有限，难以保障估测碳储量在大区域范围的精确性，故今后的工作在确保样本量的基础上，探索大样本量情况下模型精度的变化，以进一步提高人工林碳储量估测的精度。

5.   结　　论

本研究使用无人机激光雷达数据通过线性模型和非线性模型对内蒙古赤峰市喀喇沁旗旺业甸人工林的森林碳储量估测。研究结果表明：加入垂直结构变量可提高森林参数模型的拟合效果，碳储量最佳的模型中使用了两个高度变量，一个密度变量和一个垂直结构变量，R²为0.86，rRMSE为9.91%。非线性回归模型的拟合效果和稳定性要优于线性回归模型。该研究具有一定的局限性，通过线性和非线性模型分析不同变量组合，其预估精度较好，需要选取不同类型的人工林进一步分析，评价模型迁移性和林业应用潜力。