联合LiDAR和多光谱数据森林地上生物量反演研究

巨一琳，姬永杰，黄继茂，张王菲

(1.西南林业大学林学院，云南 昆明 650224;2.西南林业大学地理与生态旅游学院，云南 昆明 650224)

森林是陆地生态系统的主体，在生物圈的生物地球化学过程中发挥重要的“缓冲器”和“调节器”的功能，对改善生态环境、维护生态平衡有着重要的作用[1]。森林地上生物量(above ground biomass, AGB)作为森林生产力的重要评价指标，是评估森林碳收支的重要参数，是系统发挥其他生态功能的物质基础，也是陆地生态系统进行碳循环的重要内容[2-3]。森林地上生物量的传统调查方法是通过现地测量森林的胸径及树高等因子，使用异速生长方程或蓄积量-生物量模型计算森林地上生物量。基于资源清查的人工调查方法可准确估算森林地上生物量，但其需要花费大量时间、人力和物力，还会对森林造成一定的干扰与损伤，且只适于较小区域的森林AGB估算。不同于传统资源调查方法，使用遥感技术立足不同传感器获取影像，可快速准确、连续、动态、长期地估算森林AGB，同时可在大区域范围内实现森林AGB的估算。基于遥感技术的森林AGB反演有光学遥感、激光雷达、微波遥感等方式，在不造成森林生态损伤的前提下可极大提高对森林AGB反演的效率和精度[4]。

多光谱遥感具有丰富的光谱特征及植被指数，但其穿透性差，难以获取森林垂直结构特征，并且估算森林生物量时容易出现易饱和、敏感性差的现象[5-6]。激光雷达(light delection and ranging, LiDAR)具有穿透植被叶冠，并能直接获取目标的三维坐标，在地形、林木高度和植被空间结构的探测方面具有极大优越性[7]。将多光谱遥感与LiDAR数据相结合，进行森林AGB反演是当前研究的热点之一。Popescu等[8]融合机载LiDAR数据和多光谱数据，结合地面调查数据建立森林生物量反演模型，结果显示两种数据的融合比单一数据反演AGB精度更高。徐婷等[9]借助LiDAR和Landsat8 OLI多光谱数据采用线性多元逐步回归法分别建立森林生物量反演模型，结果表明两种数据建立的综合模型反演精度更高。卜帆[10]结合机载LiDAR和高光谱数据进行生物量反演，结果表明两种数据结合时采用支持向量机建立森林生物量反演模型精度较高，决定系数(R2)达到0.7。

利用遥感数据进行森林AGB反演通常采用参数模型和非参数模型两类方法，多元线性逐步回归为参数模型方法中反演森林AGB既常用又有效的算法。非参数模型中的KNN-FIFS(k-nearest neighbor with fast iterative features selection, KNN-FIFS)算法是韩宗涛[11]2017年基于马氏距离K-NN算法，提出以迭代方式优选遥感特征，并将特征因子优化组合的反演模型。韩宗涛等[12]以Landsat8 OLI和SAR数据为数据源，采用KNN-FIFS算法对大兴安岭根河研究区森林AGB进行估测，结果表明KNN-FIFS算法精度显著优于多元线性逐步回归法(SMLR)及K-NN算法，R2达到0.77。张少伟等[13]采用KNN-FIFS算法结合2009年和2014年两期主、被动遥感数据，进行特征组合及快速优化构建估测模型，结果表明KNN-FIFS算法可大幅度提升高维度、多模式遥感特征优选效率。

本研究以LiDAR数据和Landsat8 OLI数据为数据源，结合地面实测数据，进行数据源的不同组合，并筛选优化特征，再以多元线性回归及KNN-FIFS两种算法建模，对研究区典型的寒温带落叶松林进行森林AGB反演，最后，对各自模型的适宜性及精度及准确性进行对比分析，从而为推动大兴安岭地区典型森林AGB动态监测提供参考。

1 材料与方法

1.1 研究区概况

研究区位于内蒙古根河市大兴安岭森林生态系统国家野外科学观测研究站(121°30′～121°31′E，50°49′～50°51′N)，面积为102 km2，高程范围为810～1 116 m，为我国目前纬度最高的森林生态系统野外科学观测站。研究区属寒温带大陆季风气候，年均气温-5.3 ℃，气温年较差47.4 ℃，为高纬度多年冻土、寒温带森林生态系统的典型区域。该生态站森林植被覆盖度达75%以上，优势树种有兴安落叶松(Larixgmelinii)、樟子松(Pinussylvestrisvar.mongolica)、阔叶树白桦(Betulaplatyphylla)和山杨(Populusdavidiana)等。研究区地理位置见图1。

1.2 研究区数据获取

1.2.1 机载LiDAR数据获取

本研究使用的机载LiDAR数据以“运-5”为机载平台，载有Leica机载雷达系统，获取时间为2012年8—9月，平均飞行高度2 700 m，共32个航带，激光脉冲发射频率为100～200 kHz，扫射角度为±35°。传感器记录了首次和末次激光返回信息，记录点云数据的三维坐标信息、高程值、强度信息以及返回类型等信息[14]。

1.2.2 Landsat8 OLI数据获取

本研究使用的多光谱数据为Landsat8 OLI数据中的B1—B7波段，从地理空间数据云官方平台获取。考虑到2013年8—9月研究区光学遥感影像云层覆盖比较严重，另选取成像时间为2013年10月19日、影像无云层覆盖、航带号为123/24的一景数据，具体参数见表1。

表1 Landsat8 OLI陆地成像仪参数

1.2.3 样地调查数据获取

研究区样地调查数据分别获取于2012年37块30 m×30 m固定样地和2013年18块45 m×45 m临时样地，共55块。样地位置采用差分GPS定位，精度控制在1 m以内。样地进行每木调查，因子为胸径、树高、枝下高、冠幅、相对坐标，同时记录样地林分特征，灌木、草本的种类以及高度等参数。样地优势树种包括兴安落叶松和白桦两种，样地内两类树单木生物量的计算采用式(1)和式(2)得到，样地生物量采用陈传国等[15]的幂指数异数生长方程式(3)得到。

(1)

(2)

W样地=∑W1+W2+……+Wn。

(3)

式中：W为森林生物量，DBH为实测胸径。

1.3 遥感数据处理

1.3.1 LiDAR数据预处理

LiDAR数据处理前期部分由数据供应商完成转换、拼接及质量检测；后期部分主要借助于数字绿土LiDAR360软件。主要包括以下步骤：点云数据定义坐标系、点云去噪、提取地面点、提取植被点、提取DEM以及点云归一化。归一化后得到的冠层高度即去除地形高度影响得到的真实树高，最后将点云特征提取到样地尺度上。

1.3.2 Landsat8 OLI数据预处理

本研究获取的Landsat8 OLI数据已完成地形的几何校正，因此在预处理影像时不用做几何校正。使用ENVI5.3软件进行辐射定标和大气校正，通过辐射定标可以实现DN值与辐射亮度、反射率值等物理量之间的转化；大气校正消除大气和光照等因素对地物反射的影响从而获取其真实的地表反射率。

1.4 遥感影像特征提取

机载LiDAR第1次回波信息鲁棒性较强，将第1次脉冲回波形成的参数应用在森林AGB反演结果中效果较好[4,16]。

Landsat8 OLI影像提取4类变量因子，分别为单波段因子、植被指数、主成分分析因子、纹理因子。单波段因子包括B1—B7，植被指数包括归一化植被指数(INDVI)、差值植被指数(IDVI)、比值植被指数(IRVI)、土壤调节植被指数(ISAVI)、增强植被指数(IEVI)、有效叶面积指数(ISLAVI)[4]、地表反照率(RAlbedo)、B4/RAlbedo、B547、B65、B74等变量。主成分分析因子提取了前3个波段(PCa1—PCa3)；纹理因子提取了B1—B7波段的8种纹理信息，包括均值(Me)、方差(variance，Var)、对比度(contrast，Con)、信息熵(entropy，En)、协同性(homogeneity、Homo)、相异性(dissimilarity,Dis)、二阶矩(second order moment，Sec)、相关性(correlation，Cor)[17]。LiDAR和Landsat8 OLI变量含义及公式见表2。

表2 LiDAR和Landsat8 OLI特征变量汇总

1.5 变量筛选和建模分析

1.5.1 偏最小二乘法及变量投影重要性分析

偏最小二乘回归是对多元线性回归模型的一种扩展，是在普通多元回归的基础上增加主成分分析、典型相关分析的思想[17], 很好地解决了只用一个线性模型来描述独立变量Y与预测变量组X之间的关系。

变量投影重要性(VIP)是基于偏最小二乘法基础上进行变量筛选的过程，用来判断单个自变量在解释因变量的重要性，如果自变量对因变量的重要性越大，则VIP值越大；如果各自变量对y的解释作用一样，所有的VIP值均为1；有学者建议用0.8作为临界值来区分重要变量与不重要变量[18]，其计算公式为：

(4)

式中：k为自变量个数，ch为相关自变量提取的主成分，r(y,ch)为因变量和主成分的相关系数，表示主成分对y的解释能力，whj为自变量在主成分上的权重。

1.5.2 线性多元逐步回归

线性多元逐步回归是一种线性多元回归模型进行变量筛选的算法[19]，其回归模型为：

Y=β0+β1x1+β2x2+……+βmxm+ε。

(5)

式中：Y是生物量的值；x1、x2……xm是预测变量，β0为常数项，β1、β2……βm为回归方程的系数；ε是随机误差值。

1.5.3 KNN-FIFS方法

本研究采用韩宗涛提出的KNN-FIFS方法，是在K-NN方法的基础上进行特征因子优化组合，提高森林参数反演效率[11]。KNN-FIFS组合优化基本原理如下：

1)设样地数为n，特征数为m，由样地数据和遥感特征提取训练数据，即样地对应特征所在像元的值；

2)初始化最优特征子集为空集，因此最优模型均方根误差理论上为最大值；基于K-NN法, 依次利用特征建立森林AGB反演模型，得到最优特征子集个数的K-NN反演模型及每个模型对应的均方根误差(RMSE)；

3)得到的最优RMSE, 即RMSE最小值, 设该值为σ(RMSEb), 将研究区RMSE最大值设为RMSE0[σ(RMSE)]。若σ(RMSEb)<σ(RMSE0)则将σ(RMSEb)赋给σ(RMSE0), 并将RMSEb对应的特征子集赋给特征子集，反之迭代结束。

KNN-FIFS算法在迭代运行中，当距离度量标准确定时,K值将影响反演结果，而K值则受样地信息以及参数等因素影响[12]。其中K值默认为1～11，光谱特征信息提取窗口为1～11。

1.5.4 精度验证评价方法

反演精度评价采用决定系数(R2)和均方根误差[RMSE，式中记为σ(RMSE)]，以及相对均方根误差[RMSEr，式中记为σ(RMSEr)]来检验，R2越接近1，代表着模型精度越高；精度验证采用留一交叉验证方法。RMSE值越低，说明回归模型更准确[20]。计算公式如下：

(6)

(7)

(8)

2 结果与分析

本研究先采用偏最小二乘法对不同的遥感数据源组合方式提取的参数进行信息重要性排序，再依据其重要性排序确定最优参数组合，最后使用线性多元逐步回归和KNN-FIFS算法构建森林AGB反演模型。

2.1 遥感影像参数优选

采用偏最小二乘法回归算法进行参数优选，依据变量投影的重要性(VIP)做变量筛选，选择0.8以上的变量后续分析建模。将LiDAR和Landsat8 OLI所提取的参数作为自变量输入SIMCA 14.1中进行偏最小二乘法回归算法的建模和VIP值排序，筛选VIP>0.8的影像参数进行后续森林AGB反演模型构建。单一LiDAR遥感数据时筛选得出18个参数；单一Landsat8 OLI遥感数据时筛选得出38个参数；LiDAR、Landsat8 OLI两种遥感数据组合时筛选得出33个参数(图2)。

2.2 森林AGB反演模型构建

将在单一LiDAR、单一Landsat8 OLI及LiDAR & Landsat8 OLI联合遥感数据3种组合形式下的优选参数及地面实测样地数据，分别作为自变量和因变量输入线性多元逐步回归和KNN-FIFS模型，使用R2和RMSE分析各模型反演的结果，并分析各模型精度和适宜性。

2.2.1 线性多元逐步回归模型建立

在多源遥感数据分别是单一LiDAR、单一Landat8 OLI、LiDAR & Landat8 OLI 3种组合时，将使用偏最小二乘法筛选而得的18、38、33个参数输入线性多元逐步回归模型，最后进一步筛选得出以下特征和形成最优模型。模型结果及最优参数见表3和图3。

表3 多元逐步回归参数方程

2.2.2 KNN-FIFS模型建立

在两种遥感数据源及3种组合方式下，将以偏最小二乘算法筛选出的参数输入KNN-FIFS模型。KNN-FIFS算法的最佳模型参数见表4，进一步特征筛选结果见图4。在单一Landat8 OLI、LiDAR & Landat8 OLI两种数据组合方式下，均筛选得出B65参数,说明该参数对于森林AGB相关性较大且较为敏感。这是因为植被在红光波段有着很强的吸收特性,同时在近红外波段有着很强的反射特性，故选择近红外和红光波段进行植被遥感监测效果最佳。另在两种数据组合最后筛选参数中分别包含了冠层密度和高度分位参数，说明这两种体现森林水平和垂直结构因子的参数，在反演森林AGB中发挥着关键性的作用。

表4 KNN-FIFS特征组合

2.3 森林地上生物量反演结果对比

由模型评价结果(表5)可知模型反演精度如下：①单一LiDAR数据的线性多元逐步回归模型和KNN-FIFS模型反演精度相差不大，R2、RMSE以及RMSEr仅相差0.02、2.05 t/hm2和1.57%，总体上线性多元回归模型的精度更好，对于数据的适宜性更好。②单一Landsat8 OLI数据的线性多元逐步回归模型和KNN-FIFS模型反演精度相差较大，R2、RMSE以及RMSEr已相差0.36、9.64 t/hm2以及15.42%，KNN-FIFS模型明显具有较高的精度及对Landsat8 OLI数据的适宜性。③LiDAR & Landsat8 OLI数据组合的两种模型反演精度均较高，分别为线性多元逐步回归模型和KNN-FIFS模型在各种数据组合方式的最高。两种模型的反演精度也相差较小，R2、RMSE和RMSEr相差只有0.04、2.99 t/hm2和2.36%。

表5 模型评价结果

通过以上反演结果分析表明：①采用LiDAR和Landsat8 OLI两种数据相结合，构建线性多元逐步回归和KNN-FIFS模型均可达到最高的反演精度，机载激光雷达和多光谱遥感数据的多源数据结合可发挥各自数据优势，能表征森林水平和垂直结构，最大程度反演森林的三维结构，进而提高森林AGB反演精度和准确度。②采用参数优选技术可有效提高森林AGB反演精度。首先采用偏最小二乘算法对两种数据源3种组合方式时总体参数筛选，后又在两种模型构建的同时进行了特征优选，两次的参数优选过程去除了和森林AGB相关性低、敏感性差的参数。③在单一LiDAR和LiDAR&Landsat8 OLI两种组合下，针对该研究区典型寒温带落叶松林森林AGB，使用两种模型的反演精度几乎相差不大，说明模型在该数据组合及森林类型模型适宜性相似。而在单一Landsat8 OLI时，KNN-FIFS反演精度明显高于线性多元回归模型，KNN-FIFS模型的适宜性更强。

由散点图(图5)分析可知，两种模型在两种数据源的3种组合方式下的估测值均出现不同程度的高估和低估现象，其中单一LiDAR数据源时，两种模型反演的结果均出现明显的低估现象，而LiDAR & Landsat8 OLI组合下高估与低估现象却不明显。关于估测值与实测值的线性拟合性上，森林AGB值约80 t/hm2以下时，两种模型在两种数据3种组合方式下总体上均表现出较好的线性拟合性，实测值与估测值之间相关性较高；但随着森林AGB数值的增加，拟合性表现得较为离散。这是因为机载LiDAR对森林的穿透能力，能够提供准确的树高和直观的冠层垂直结构信息，结合反应森林水平结构的光谱信息，表现得在单一机载LiDAR数据和两种数据结合反演森林AGB中精度为较高，线性拟合性较好。而在单一Landsat8 OLI数据时，整体拟合性均低于LiDAR拟合结果，尤其在线性多元逐步回归模型反演时表现最差，低于KNN-FIFS反演效果。这是因为在KNN-FIFS拟合方法中，估测参数只与相邻的k个样本有关，因此样地点的加权求值可以减少样本不平衡问题，能更好地描述各参数与影像之间的非线性关系[12]。对于整体拟合结果较差，这与多光谱数据具有丰富的光谱特征及植被指数，但穿透性差，难以获取对森林垂直结构特征参数，对森林AGB反演具有较大的结构性局限性有关。

2.4 森林地上生物量空间分布

基于3种数据源采用线性多元逐步回归算法和KNN-FIFS算法进行单一反演以及联合反演得到航飞区森林地上生物量制图(图6)。

从图 6a—6c可以看出，基于多元逐步回归算法的LiDAR数据的制图结果中，生物量分布<40、≥40～70以及≥160～200 t/hm2的森林区域各占30%，其余占10%；采用Landsat8 OLI数据的制图结果中生物量≥40～70 t/hm2的区域约占研究区的40%，≥100～130 t/hm2的区域约占研究区的40%，<40 t/hm2的区域约占研究区的20%；采用融合LiDAR & Landsat8 OLI数据的制图结果中生物量≥40～70 t/hm2的区域约占研究区的50%，≥130～160 t/hm2的区域约占研究区的30%，<40 t/hm2的区域约占研究区的20%。另外结合实际土地利用现状(图6b)可发现Landsat8 OLI数据的森林地上生物量结果中的高值区出现在了林区主干道附近，有明显的错估现象。

从图6d—6f中可以看出，基于KNN-FIFS算法的LiDAR数据的制图结果中，生物量分布<40 t/hm2、≥40～70 t/hm2的区域各占研究区的40%，其余占20%；采用Landsat8 OLI数据的制图结果中生物量≥40～70 t/hm2的区域占50%，<40 t/hm2的区域约占研究区的20%，≥130～160 t/hm2约占30%；采用融合数据的制图结果中生物量≥40～70 t/hm2约占研究区的50%，<40 t/hm2约占研究区的30%，其余占20%。结合地面调查信息来看，基于KNN-FIFS算法反演出的生物量制图与调查信息基本相符，<40 t/hm2的区域和≥40～70 t/hm2的区域占比较大。另外结合根河生态站航飞区森林地上生物量高值区主要分布在航飞区北部以及山脊线附近，低值区主要分布在海拔较低的林区主干道。航飞区域森林地上生物量的空间分布趋势与实际的地貌特征相符。航飞区北部处于山坡的阴坡处，可以较好地存留水分，因此森林地上生物量较高；在航飞区域的林区主干道附近，由于修建主干道这种工程可以改变土层结构以及条件，因此类似的人为干扰破坏植被活动会导致森林地上生物量偏低。

3 讨 论

1)森林地上生物量的监测对区域生态植被保护及生态系统平衡都有着重要的意义，相比采用传统人工地面调查森林AGB方法，使用遥感技术可快速获取大区域尺度、实时及经济的森林AGB分布情况。本研究以LiDAR和Landsat8 OLI影像为数据源结合地面实测数据，对不同组合形式下的遥感参数进行筛选，并采用线性多元逐步回归和KNN-FIFS算法，对航飞区森林AGB进行反演及分析研究。

研究结果中采用单一LiDAR数据的线性多元回归和KNN-FIFS模型反演结果的R2为0.76和0.74。庞勇等[21]将机载LiDAR数据划分为针叶林、阔叶林和针阔混交林3个类型进行森林AGB反演，决定系数可达0.82以上，优于笔者的建模结果。本研究采用的机载LiDAR数据虽能准确获取森林垂直结构信息，但仍受真实地形起伏的影响，使归一化的点云高度产生误差，又未考虑森林类型的影响，在一定程度上影响了本研究在单一LiDAR数据情况下的精度。

Landsat8 OLI数据的两种模型的反演精度R2为0.24和0.60，反演精度结果相比LiDAR数据结果较低。Foody等[22]研究表明体现光学特征的光谱信息在不同的地物之间存在较大差异，地形存在起伏较大区域或郁闭度较低区域，光谱信息受到地形等影响较大导致不能准确表述植被特征，使得森林AGB数值产生10～15 t/hm2的反演误差。罗洪斌等[4]同样结合Landsat8 OLI和LiDAR两种数据采用偏最小二乘法对橡胶林地上生物量进行遥感反演中，采用单一Landsat8 OLI数据建模精度为0.56，同样比单一LiDAR数据时建模精度低。上述学者在采用单一Landsat8 OLI数据进行森林AGB反演时精度与本研究结果大致相当，多光谱数据反演精度低的原因多与光谱信息不能反映森林垂直结构、受地形影响较大及采用数据分辨率低等密切相关。

本研究融合LiDAR & Landsat8 OLI两种数据源采用两种模型反演森林AGB，反演精度R2分别0.84和0.80。胡凯龙等[23]利用光学数据的纹理特征和机载激光雷达的点云特征进行森林地上生物量反演分析，在不区分森林类型的情况下，建模精度低于本研究的建模结果，分别为0.73和0.79。徐婷等[9]基于机载LiDAR和光学数据进行建模，分析在不同森林类型的情况下估算综合森林地上生物量，得到的R2可以达到0.88，证实了将典型森林类型划分为不同森林类型会提高建模精度。本研究在未考虑森林类型的情况下，联合两种数据源建模总体上提高了反演的精度，这说明LiDAR数据和多光谱数据联合、参数组及特征优选对我国寒温带森林AGB反演精度提高是行之有效的，另采用两种模型也是适宜的。

2)本研究通过对LiDAR数据和Landsat8 OLI数据提取出的特征变量，采用偏最小二乘法进行特征优选，再通过线性多元逐步回归算法和KNN-FIFS算法反演森林地上生物量形成结论如下：LiDAR和Landsat8 OLI两种数据源的融合相比单一数据源对森林AGB的反演精度更高；在两种数据3种组合方式下进行特征的优选，采用偏最小二乘法算法进行筛选对反演精度的提高是有效的；线性多元回归及KNN-FIFS两个模型在不同组合下，对本研究区的森林类型反演是适宜的。

在未来的研究中可将光谱分辨率高的高光谱数据纳入，加大不同数据的优势互补，深刻反馈森林三维结构的细化特征；还可尝试参数及非参数其他模型的算法研究，寻找最优模型。另在加大特征提高的同时，应深入其他筛选算法，去除冗余特征降维提取相关性最大特征；对于热带森林等不同森林类型也应加大树种、模型适宜性研究及精度和饱和点提高的研究。

根据本研究结果看，采用偏最小二乘法进行特征变量的筛选，可以为后续的森林地上生物量建模提高精度。从特征变量VIP值排序上看，LiDAR百分位高度变量特征在排序中较靠前，反映了其在森林地上生物量有一定优势。此外，在纹理因子的提取中，可以考虑设置窗口的大小以便更好地进行特征选择。