三峡库区森林叶面积指数多模型遥感估算

董立新

(1.中国遥感卫星辐射测量与定标重点开放实验室，北京 100081； 2.国家卫星气象中心，北京 100081； 3.中国气象局卫星应用联合研究中心，北京 100081)

0 引言

叶面积指数(leaf area index，LAI)是定量研究森林生态系统能量交换的一个重要结构参数，定义为单位地表面积上的叶子表面积总和的一半[1]。LAI地面测量可通过收获法及异速生长模型估测，或利用基于冠层辐射传输模型的观测仪器直接获得[2]。但地面测量较难获取大范围LAI。因此，遥感估算区域尺度LAI是十分有意义的工作[3]。遥感估算主要有2类： 一是植被指数法，即利用植被指数与LAI的关系； 另一类是模型反演，主要是通过植被冠层模型，如辐射传输模型、几何光学模型和计算机模拟模型等反演[4]。

植被指数法是利用植被指数与实测LAI之间的关系，实现区域尺度的LAI估算。如Chen等[5]利用LAI与增强植被指数(enhanced vegetation index,EVI)和归一化植被指数(normalized difference vegetation index，NDVI)的关系实现了针叶林LAI估算； 朱高龙等[6]利用Landsat5 TM数据计算植被指数与LAI之间的关系，建立了区域森林LAI估算模型； 刘婧怡等[7]提取10个植被指数，在与实测LAI相关性基础上进行了区域LAI估算； 韩婷婷等[8]在2014年利用TM数据反演了西双版纳地区森林LAI； 刘振波等[9]在2015年基于植被指数与实测LAI构建小兴安岭冠层LAI模型； 姚雄等[10]在2017年通过计算12种遥感植被指数，分析了实测LAI数据和植被指数的相关性，构建了林地LAI估算模型。还有许多研究者都利用植被指数与LAI的关系，计算了不同植被类型的LAI[10-17]。

反射率受到大气、土壤和植被类型等因素的影响，所以通过植被指数计算LAI没有通用的关系式。植被指数法虽然有一些局限性，但易于得到大范围的LAI，有很强的实用性，且精度较高[3]。针对植被指数法的局限性，提出了混合像元分解法(spectral mixture analysis, SMA)，解算端元组分在像元中占的百分比，由此减少LAI计算误差。陈丽等[18]基于SMA模型进行了森林LAI反演。

辐射传输模型通过迭代方法使损失函数值达到最小，即采用优化技术使得模型结果与实际相符。最早用于植被遥感的模型是SUIT模型[19]。在此基础上，SAIL模型考虑了LAI、平均叶倾角及冠层组分的透射、反射和土壤反射等参数，是一个具有代表性的模型。Kuusk[20]对SAIL模型又进行了改进，考虑了热点效应，由此产生了SAIL模型的另一个版本SAILH，其在计算单次散射贡献对二向反射率的贡献时，考虑了叶片的尺寸以及相应的阴影影响。

几何光学模型解决的是植被几何结构和空间分布模型化等问题，用结构参数(株密度、树冠大小、高度等)表达几何结构。该模型考虑了地物的太阳承照面和阴影，计算光照植被、阴影植被、光照地面和阴影地面4个分量。其中最有代表性的为Li-Strahler几何光学模型[21-22]和Jupp等[23]提出的模型。为了能够使模型适用于浓密的林地，Li等[24]进一步提出了考虑树冠相互阴影效应的GOMS模型； 其后，Li等[25-26]在间隙率模型的基础上，进一步提出了几何光学—辐射传输混合模型，实现了几何光学模型向辐射传输模型的逼近。

辐射传输模型和几何光学模型等都具有一定的物理基础，且独立于植被类型的影响，有一定的普适性。但有些反函数不收敛，会造成LAI的估算错误，同时计算复杂耗时。计算机模拟法采用蒙特卡罗方法(Monte Carlo)进行光子跟踪，但所需参数太多，而且计算量较大，不适用于大面积的LAI反演。本文以三峡库区为研究区，利用植被指数法针对不同森林类型建立多个LAI反演模型，以期估算区域尺度LAI。

1 研究区概况及数据源

1.1 研究区概况

三峡库区是一个特定的区域概念(图1)，泛指175 m水位方案淹没涉及的20个县、市，位于长江上游低位地区，北靠大巴山，南邻云贵高原。地理位置位于E106°00′～111°50′，N29°16′～31°25′之间，面积约5.8×104km2。开始于湖北省宜昌市，向西至重庆市江津区，著名的三峡大坝位于宜昌市西部边界。

图1 研究区位置及野外样点分布示意图

研究区属于典型的亚热带季风气候区，冬季温暖，夏季干热，秋季湿度较大，多云雾。年平均气温为15～19 ℃，平均年降水量为1 140～1 450 mm。库区河谷地与平坦地区占4.3%，丘陵占21.7%，山地占74%。由于地处我国亚热带常绿阔叶林区，主要的植被类型包括常绿针叶林、落叶针叶林、常绿阔叶林、混交林、灌木林和农作物，在林内通常有一至数个优势种，常分为2个乔木层。区内人烟稠密，农业与林业活动历史悠久。现有的森林中人工林、次生林较多，由于自然演替与人为干扰的原因，使演替的系列、阶段与速度的变化较大，加之树种较多，因而较之我国东北针阔叶混交林与华北暖温带落叶阔叶林区域情况更为复杂。

1.2 数据源及其预处理

1.2.1 野外数据

所用野外数据是在“三峡工程生态与环境遥感动态与实时监测”项目支持下，参照《中国生态系统研究网络观测与分析标准方法——陆地生物群落调查观测与分析》以及《测树学》[27]的调查方法，在2003—2007年间对土地利用及LAI进行多次广泛调查(图1)，主要包括龙门河自然保护区(2003年4—6月)和三峡库区(2006年8—9月； 2007年9—10月)。

选取林相比较完整且有代表性的森林设置样地[28]，样地的分配兼顾不同的森林类型和地形条件，在龙门河自然保护区分别布设样地40块，样地大小为100 m×100 m，每样地内按系统抽样法设立5个20 m×20 m样方，对样地内胸径大于5 cm的乔木进行每木检尺，内容包括树木名称、胸径、树高、冠幅、枝下高、覆盖度和LAI，同时记录了当地的环境条件； 在三峡库区调查样地共28个，在距林分中心60～100 m的范围内设置4～5个样地(成熟林为100 m，其他为60 m)，样地大小为20 m×20 m，每样地再分成4个样方进行乔木层调查。

1.2.2 LAI野外测量

测量仪器为Li-Cor公司的LAI-2000 Plant Canopy Analyzer。实际测量中，选取有代表性的样地，样地内的测量采用随机抽样布点法。与光学遥感影像的空间分辨率(30 m×30 m)对应。布设考虑了仪器的最大覆盖范围，避免测量范围的重叠。

LAI-2000要求测量尽量选择在阴天,以避免晴天条件下的直射光影响。实际操作中，当天空和植物冠层的条件不理想的情况下需要调整测量的方式，如通常需要测量多个数据进行均值处理得到LAI的值； 同时，太阳和操作员不能在传感器的视角里； 对于树叶很浓而又有大的空隙时需要传感器使用窄的视角，可以把树叶和空隙结合起来考虑。

1.2.3 遥感数据与预处理

光学遥感数据主要为2002年覆盖三峡库区的7景Landsat TM影像数据，成像日期分别是2002年8月29日、2002年8月30日、2002年9月1日和2002年10月2日。由于光学遥感传感器受到太阳高度、地形和大气等因素的影响，要得到地物真实反射率必须先经过数据预处理，主要包括： 辐射校正、大气校正、几何纠正和地形光谱校正。

本文采用ERDAS9.0软件中的ATCOR2模块进行TM影像的大气校正； 并采用控制点校正方式，以三峡地区1∶ 50 000比例尺的地形图作为参考，对TM数据进行了几何纠正。本文采用最近邻法进行重采样，以保证光谱信息不变，重采样后的文件为高斯—克吕格投影。坐标转换时的精度均控制在0.8个像元以内。

考虑到三峡地区地形起伏较大，对反射率有很大影响，特别是在陡峭地区还会形成太阳光在坡面上的多次反射，因此，对TM数据采用Minnaert校正法[29]进行地形光谱校正。

1.2.4 辅助数据

辅助数据主要为土地利用图，是利用30 m空间分辨率Landsat TM影像制作的2002年三峡地区1∶ 50 000土地利用图[30]，经过2003年8月进行野外验证，2002年分类结果总体精度为84.31%，Kappa系数为0.75[31]。本文从中提取出针叶林、阔叶林和混交林等地物类型。

2 植被指数法LAI估算方法

本文首先利用TM数据计算7种常用植被指数和5个自定义植被指数，并结合野外观测的LAI，通过多元逐步回归与主成分分析(principal component analysis,PCA)等多个回归模型筛选不同森林类型的LAI估算模型，最终建立了不同森林类型的LAI估算模型，实现区域尺度森林LAI多模型估算。

2.1 植被指数计算

植被指数是指由多光谱遥感数据经线性和非线性组合而构成的数值[32-33]，是反映地表植被生长、覆盖、生物量和种植特征的间接指标。本文计算的植被指数包括7个常用指数与5个自定义指数。常用指数包括： NDVI[34]、垂直植被指数[35](perpendicular vegetation index，PVI)、土壤调节植被指数[36](soil-adjusted vegetation index，SAVI)、修正的土壤调节植被指数[37](modified soil-adjusted vegetation index，MSAVI)、EVI[38]、大气阻抗植被指数[39](atmospherically resistant vegetation index，ARVI)和土壤调节大气耐抗植被指数[39](soil adjusted and atmospherically resistant vegetation index，SARVI)。

植被指数与所用的波段组合方式有很大关系，为了充分利用其他波段，以突出植被信息为原则，在反复试验的基础上，本文自定义了5个指数，计算公式分别为

(1)

(2)

(3)

(4)

(5)

式中ρbi为Landsat TM数据第i波段的反射率。

2.2 PCA

在实际问题中，多变量问题中不同变量之间存在一定相关性，由于变量较多，势必增加问题的复杂性。PCA将原变量重新组合成新的互相独立的几个综合指标，同时可根据需要从中选取较少指标尽可能地反映原变量的信息。PCA也是一种降维处理，抓住了问题的主要矛盾，并提取了新的信息，有利于实际问题分析处理。PCA可以单独处理一些问题外，也可以与回归分析结合起来，即主成分回归。

设原始矩阵X的p个向量作线性组合，综合指标向量简写为

Fi=a1iX1+a2iX2+…+apiXp，a1i2+a2i2+…+api2=1。

(6)

根据协方差矩阵求出特征值λi、主成分分量贡献率和累计方差贡献率，确定主成分分量个数，并将特征值λi按大小顺序排列。特征值是各主成分分量的方差，反映了各个主成分分量的影响力。根据选取主成分分量个数的原则，选取累计贡献率达80%～95%的特征值λ1，λ2，…，λm，其中整数m即为主成分分量的个数。主成分分量贡献率β和累积贡献率B的计算方式分别为

(7)

(8)

然后，建立初始变量载荷矩阵，解释主成分分量。主成分分量载荷P(yi,xj)实际上也是主成分分量yi与各变量xj之间的相关系数R，即

(9)

式中Lij为单位特征向量。

最后，建立主成分回归模型进行计算分析。

3 结果与验证

3.1 植被指数与实测LAI相关性分析

本文计算了植被指数与LAI的R值(表1)。由表1可见，VI1，NDVI，SAVI，ARVI，MSAVI及SARVI等植被指数与LAI的相关性相对较高。但不同的森林类型，其相关性稍有不同。由此依据，可进一步建立植被指数与LAI的回归关系模型。

表1 各种植被指数与LAI的R值

注： 显著水平P≤0.05。

3.2 基于植被指数的LAI估算回归模型筛选

3.2.1 一元回归分析

根据表1初步筛选植被指数，针叶林选用VI1，NDVI，SAVI，ARVI，MSAVI和SARVI共6个植被指数； 阔叶林则选用VI1，NDVI，MSAVI和SAVI共4个植被指数，虽然VI3的R值也达到了0.34，但在实际回归计算中效果不太理想，故未做一元回归分析； 混交林选用NDVI，VI1和VI2共3个植被指数。同时，进一步建立了一元线性回归分析模型，并统计了R、判定系数(R2)等变量(表2)。

表2 各种植被指数与LAI的一元线性回归分析

注： 显著水平P≤0.05。

在一元线性回归分析后，尝试一元非线性回归拟合。选用LAI与NDVI的多个拟合方程进行比较，确定最能充分描述数据的模型。主要有： 线性模型、二次多项式、复合模型、生长模型、对数模型、三次多项式、S曲线、指数模型、双曲线模型和幂指数模型。针叶林的拟合结果如图2和表3所示，图2中NDVI值域拉伸至[0，255]。通过拟合模型的R2对各种模型的优劣进行比较，由表3可以发现，NDVI与LAI的二次多项式模型拟合的R2为0.518，远大于其他模型结果，似乎应采用二次多项式对不同森林类型的LAI进行非线性拟合。但从图2可以发现，在NDVI值的低位(179)，LAI的值竟然达到了2以上，这与现实不相符。

图2 针叶林LAI-NDVI不同模型拟合图

表3 植被指数与LAI的非线性回归模型R2(针叶林)

而线性模型拟合与双曲线模型比较，具有较高的R2，但双曲线模型更能表达LAI随NDVI的变化趋势，这与以往的研究结果一致[40]。所以对于针叶林，采用双曲线模型进行分析(表4)，分析可见，其R2最大为0.634，不太令人满意，故尝试利用多元逐步回归分析。

表4 植被指数与LAI的双曲线分析模型

注： 显著水平P≤0.01。

3.2.2 多元逐步回归分析

由以上结果可知，一元回归分析模型在对LAI预测中的结果不尽如意，所以，考虑对其进行多元逐步回归分析。分析发现，多元逐步回归分析模型比一元回归分析模型具有更好的R值和预测能力(表5)。经多元逐步回归分析，针叶林LAI对植被指数VI1，VI2与VI4较为敏感，R达到0.803，R2判定为0.644，经t检验，VI1与VI4的检验值均达0.000，按0.05水平均有显著性意义。其他类型比针叶林结果较差，原因可能是各指标之间存在波段相关性，故尝试采用PCA法进行建模。

表5 各种植被指数与LAI的多元逐步回归分析

注： 显著水平P≤0.05。

3.3 PCA模型

利用PCA方法，对针叶林、阔叶林和混交林分别进行了第1到第4主成分分量的提取(表6)，分析了各主成分分量贡献率(表7)，从而建立了基于各植被指数的LAI估算的主成分回归模型(表8)。表中分量为不同森林类型进行PCA而导出的主成分分量Ci(i=1,2,3,4)。

表6 不同森林类型的PCA结果

表7 PCA各分量贡献率

表8 各种植被指数与LAI的主成分回归模型

注： 显著水平P≤0.05；Ci(i=1,2,3,4)由主成分分量载荷确定。

由表8可见，与多元逐步回归分析相比，针叶林的主成分回归模型R反而有所下降，但是阔叶林和混交林的R值有所上升。根据以上分析，发现对于不同的森林类型，可采用不同的模型进行估算。根据R2的大小，在研究区不同森林类型LAI预测中，针叶林采用多元回归模型，阔叶林与混交林都采用主成分回归模型。

3.4 结果及验证

将多模型应用到TM影像中，计算得到了研究区LAI结果(图3)。可以看出，研究区大部分LAI值范围在1～9之间，整体上呈“中北低、东南高”的分布格局； 东部林地LAI值相对较高，大部分区域的LAI值高于4，主要原因在于该区分布龙门河自然保护区，林分结构以针叶林和阔叶林为主； 南部林地LAI值也相对较高，该地区地形复杂，森林覆盖率非常好； 中北部地区林地LAI值相对较低，原因在于该区植被稀疏，过度开发使得林下水土流失严重。可见，基于多模型估算的LAI值较好地还原了三峡库区LAI值的空间分布状态。

图3 三峡库区LAI结果

对不同的森林类型，其预测精度可用总均方根误差(root mean square error, RMSE)来评价(图4)。LAI的精度通过与样本中预留的观测样地(分针叶林、阔叶林和混交林)进行检验，其RMSE分别为0.829 4,1.111 5和1.790 9。实测值与预测值比较如图4所示，R2也达到了0.77以上，结果比较满意。

(a) 针叶林(b) 阔叶林 (c) 混交林

4 结论

区域尺度LAI遥感反演结果是森林生态系统和全球碳循环研究的基础参数，本文通过多个模型估算得到三峡库区区域尺度森林LAI分布结果，主要结论有：

1)利用Landsat TM计算的7种常用植被指数和5个自定义植被指数，通过模型筛选建立了不同森林类型的LAI估算模型。其中，针叶林采用多元逐步回归模型，阔叶林与混交林采用PCA模型，最终通过多个模型估算得到三峡库区区域尺度的森林LAI。

2)通过与样地实测LAI数据比较进行精度验证，针叶林、阔叶林和混交林LAI的RMSE分别为0.829 4,1.111 5和1.790 9，R2均达到了0.77以上。

3)研究发现，森林LAI容易出现饱和现象，且利用遥感数据进行建模时，对不同季节的数据，模型也有所不同，因此进行模型筛选是非常必要的，这正是植被指数法LAI遥感估算的一大难点和重要工作。