基于MODIS时间序列数据的竹林地上生物量估算

杨绍钦，王 翔，许 澄，商天其

（1.浙江省森林资源监测中心，浙江 杭州 310020；2.浙江省林业勘测规划设计有限公司，浙江 杭州 310020）

森林地上生物量(AGB)是表征森林生产力和陆地生态碳循环的重要指标，对研究全球气候变化、人类经营活动等对森林生态系统的影响具有重要意义[1]。精确估算森林地上生物量对森林可持续发展、增加森林碳汇能力及减缓全球气候变暖具有关键作用[2]。样地调查是精确获取森林地上生物量的主要手段[3-5]，但该方法耗时费力[6]，难以实现大范围时空动态监测。近年来，众多研究者利用遥感变量信息(如光谱特征、植被指数、冠层结构参数、纹理等)，结合样地调查构建遥感估算模型，实现了对地上生物量的时空估算。森林地上生物量遥感估算有参数和非参数模型[7-12]，其中，以支持向量回归(SVR)模型为代表的非参数模型被广泛应用于森林地上生物量时空估算研究[13-15]。LI等[16]比较了不同模型(如最小二乘法模型、广义加性模型、随机森林、SVR等)在估算地上生物量上的效果后发现，SVR模型估算效果最好；FENG等[14]发现随机森林、SVR模型等非参数模型估算的地上生物量精度要高于参数模型(线性模型和乘数模型)，且SVR模型估算精度最高。

竹林具有“爆发式生长”的特点。在生长和抽枝展叶期，竹林叶面积指数(LAI)逐渐增加，地上生物量快速积累[16-17]。LAI时间序列上的变化趋势能够反映森林的季节性变化，并常作为森林生态系统碳水循环和地上生物量估算等研究的重要参数之一[18-20]。相较于归一化植被指数(NDVI)，增强型植被指数(EVI)对高生物量区域的敏感性更高，受土壤背景和大气效应的影响更小[21-22]。比值指数(RVI)是绿色植被的灵敏指示参数，能够反映叶干生物量和叶绿素含量的动态变化。目前，对竹林地上生物量的遥感估算研究多是基于单一时相遥感植被指数、光谱信息、纹理等特征进行，估算精度较低[23-24]，且不能表征竹林快速生长过程中地上生物量快速积累的动态变化[6]。本研究以浙江省竹林资源为研究对象，以2014年浙江省中分辨率成像光谱仪(MODIS)的LAI、EVI、RVI时间序列产品为数据源，结合2014年竹林地上生物量调查数据，构建SVR模型，估算浙江省竹林地上生物量，以期为竹林碳汇遥感监测提供参考。

1 研究方法

1.1 研究区概况

浙江省 (27°06′～31°11′N，118°01′～123°10′E)位于中国东南沿海、长江三角洲南翼，是中国竹林面积最大的省份之一，竹资源十分丰富，竹产业发达，素有“世界竹子看中国，中国竹子看浙江”的美誉。目前，浙江省竹林面积为94万 hm2，其中毛竹Phyllostachys edulis林83万 hm2，竹林面积占全省森林面积约15%。

1.2 研究数据

1.2.1 MODIS 数据及预处理 MODIS LAI (MOD15A2H)、地表反射率 (MOD09A1)和 MODIS NDVI(MOD13Q1)数据产品在美国国家航空航天局(NASA)官网下载。MODIS LAI和地表反射率的时间和空间分辨率分别为8 d和500 m，MODIS NDVI的时空分辨率分别为16 d和250 m。数据利用重投影软件MRT (MODIS Reprojection Tools)投影，采用最近邻域法将空间分辨率统一重采样为 500 m×500 m。基于MODIS NDVI时间序列产品和地表反射率数据提取浙江省竹林丰度信息。首先，利用NDVI数据和最大似然法提取林地专题信息；其次，掩膜提取林地的地表反射率，并进行最小噪声变换；然后，利用影像端元法提取竹林、针叶林和阔叶林端元光谱曲线；最后，利用完全约束最小二乘混合像元分解法，得到浙江省竹林丰度图。竹林丰度提取方法参考LI等[25]和MAO等[26]。

1.2.2 地面调查数据 调查和收集182个不同立竹密度、立地质量和经营条件的标准竹林样地(30 m×30 m)。调查因子包括胸径(DBH)和年龄(0～1年生的竹子是1度竹，2～3年生是2度竹，4～5年生是3度竹)，记录样地中心坐标。基于样地竹林胸径和年龄，分别计算毛竹和雷竹Ph.praecox单竹的地上生物量[19]。

其中：M毛和M雷分别表示毛竹和雷竹单竹地上生物量；DBH和A分别为胸径和竹龄。为与MODIS遥感数据匹配，本研究将样地竹林面积线性扩大为500 m×500 m[11]。

1.3 遥感变量设置

设置141个遥感变量，其中LAI、EVI和RVI时间序列数据各46个，各序列年均值各1个。由于MODIS产品易受大气、气溶胶、雪覆盖等因素影响，数据存在误差大、不连续等问题[20]，本研究利用粒子滤波同化算法提高MODIS LAI时间序列数据的精度，主要步骤为：①利用3次样条帽盖平滑算法对MODIS LAI产品进行平滑处理[27]；②将平滑的MODIS LAI输入到LAI动态模型，模拟叶面积指数；③将MODIS地表反射率的红光和近红外反射率、PROSAIL模型模拟冠层反射率、模拟LAI值输入到粒子滤波同化算法，得到叶面积指数同化结果。粒子滤波同化算法详细参考LI等[17]和MAO等[26]。增强型植被指数和比值指数利用MODIS地表反射率的第1 (red)、第2 (nir)和第3 (blue)波段分别计算。

其中：IEVI、IRVI分别表示增强型植被指数(EVI)和比值指数(RVI)，ρnir为近红外波段反射率，ρred为红光波段反射率，ρblue为蓝光波段反射率。利用mSG平滑算法对EVI和RVI进行处理，以减少噪声影响[28]。

1.4 变量筛选及 SVR 模型构建

变量数量的增加可能会导致“维数灾难”，造成运算复杂，处理速度下降[19]。为充分利用模型输入变量的特征，需要选择最优变量特征构建SVR模型。随机森林算法预测能力强，鲁棒性较好，适合处理高维数据集，计算自变量对因变量的相对重要性[29]。因此，可采用随机森林算法选择SVR模型的输入变量，即通过计算均方差增量百分数来度量和评估各变量特征对竹林地上生物量的相对重要性。随机森林参数的最优ntree值为2 000，mtry采用模型默认值，即变量个数的1/3[25]。

SVR通过选择合适的核函数将样本数据从低维空间转换到高维空间，使得低维空间内的非线性问题可在高维空间中线性处理，并保证其良好的泛化能力。本研究设置radial、linear、polynomial和sigmoid等4种核函数[30]，比较分析后选择最优核函数构建竹林地上生物量估算模型。采用随机不重复抽样法从182个竹林样本中抽取50%作为训练样本，剩余的50%作为验证样本。

1.5 模型精度评价

利用决定系数(R)、均方根误差(RMSE)和P(t检验)评价SVR模型模拟地上生物量的精度、误差和显著性水平，R较高且RMSE较低，说明模型估算结果较好。计算公式如下：

其中：N为样本个数，i=1, 2 ,···,N；ym表示模型模拟值，yo表示实测值；为N个实测值的平均值。

2 结果与分析

2.1 变量筛选

基于随机森林算法得到EVI、RVI、LAI等141个输入变量对竹林AGB影响的重要性得分，由图1可知：不同时间的叶面积指数、增强型植被指数和比值指数对竹林地上生物量影响存在较大的差异，其中LAI时间序列变量占比最大，RVI变量占比最小。在前50个影响较大的变量中，有43个对竹林地上生物量影响的重要性得分大于0.10%，LAI年均值(LAI_Ann)重要性得分最大(2.29%)，LAI时间序列变量中，LAI_153(第153天LAI值)最大，为1.44%，LAI_113最小，为0.10%。因此选取变量重要性得分大于0.10%的43个变量作为SVR模型输入变量，构建竹林地上生物量估算模型。

图1 输入变量对竹林地上生物量影响的重要性得分Figure 1 Importance score of input variables on aboveground biomass of bamboo forest

2.2 SVR 模型构建及分析

利用筛选的变量，分别以radial、linear、polynomial和sigmoid等4种核函数构建竹林地上生物量估算模型。由图2可知：基于radial函数的SVR模型的训练和测试精度(R)均最高，RMSE最低，估算地上生物量和实测地上生物量的精度(R)分别为0.76和0.72，RMSE分别为5.15和8.03 Mg·hm-2；其次是linear函数，训练和测试精度(R)分别为0.64和0.60，RMSE分别为5.95和8.38 Mg·hm-2；sigmoid函数训练和测试精度(R)最低，均＜0.15，RMSE最大，均＞9 Mg·hm-2。polynomial函数构建的SVR模型在估算竹林地上生物量的训练和测试精度均较高，但测试样本中估算地上生物量时出现负值，因此本研究选择基于radial核函数构建的SVR模型估算浙江省竹林地上生物量。

基于radial核函数构建SVR模型，分别利用单一变量时间序列数据和植被指数组合(EVI+RVI)测试模型精度。由图3可见：利用单一变量(EVI或RVI)时间序列数据构建的SVR模型估算竹林地上生物量，精度较低(R＜0.65)，RMSE较大；利用单一变量LAI时间序列数据和植被指数组合(EVI+RVI)模拟的AGB的精度(R＞0.65)较高，RMSE较小。结合图2A可知：基于radial核函数的SVR模型，利用LAI、EVI和RVI时间序列数据估算竹林地上生物量，训练精度和测试精度均最高，误差均最低。即利用LAI、EVI、RVI时间序列数据，基于radial核函数构建SVR模型，可有效模拟竹林地上生物量的时空分布。

图2 不同核函数构建的SVR模型估算竹林地上生物量的精度Figure 2 Estimating accuracy of 4 kernel functions in SVR model for aboveground biomass of bamboo forest

图3 基于 SVR 模型的 LAI、EVI、RVI及 EVI+RVI时间序列数据估算竹林地上生物量Figure 3 Model accuracy of LAI, EVI, RVI and EVI+RVI time series data for estimating aboveground biomass of bamboo forest using SVR model

2.3 竹林地上生物量空间分布估算

基于radial核函数构建的SVR模型，结合竹林丰度信息，估算得到浙江省竹林地上生物量统计变化值，由表1可知：浙江省全省竹林地上生物量均值为7.85 Mg·hm-2，总地上生物量为3.31×107Mg；浙江省竹林地上生物量在各市具有明显的差异性，其中：湖州市、杭州市、金华市、绍兴市和宁波市的竹林地上生物量均值均大于全省均值；湖州市竹林地上生物量均值最大为13.56 Mg·hm-2，总地上生物量为2.13×106Mg，占全省总地上生物量的6.44%，但其竹林面积仅占3.73%；杭州市地上生物量均值为9.86 Mg·hm-2，总地上生物量为8.48×106Mg，占全省总地上生物量的25.42%，竹林面积占20.24%；舟山市地上生物量均值最低，为5.72 Mg·hm-2，总地上生物量为6.89×104Mg，占全省总地上生物量的0.21%，竹林面积占0.29%。

表1 基于SVR模型模拟的浙江省竹林地上生物量统计值Table 1 Statistical values of aboveground biomass of bamboo forest in Zhejiang Province based on SVR model

3 讨论

本研究基于叶面积指数、增强型植被指数和比值指数等多变量，构建并优选SVR模型，实现浙江省竹林地上生物量较高精度的估算。SVR模型耦合LAI、EVI、RVI时间序列数据估算竹林地上生物量与实测值间仍存在一定的误差。当实测值较小时(＜20 Mg·hm-2)，估算值高于实测值，即低值高估；而实测值较高时(＞30 Mg·hm-2)，估算的地上生物量低于实测数据，即高值低估。可见本研究构建的SVR模型并不能完全消除地上生物量估算的饱和现象。

本研究发现：基于单一变量的SVR模型模拟竹林地上生物量效果较差，而多变量的参与尤其是LAI的参与有效提高了模型预测精度。本研究中，植被指数(EVI+RVI)构建的SVR模型模拟竹林地上生物量的精度分别为0.68和0.71，RMSE分别为5.79和8.13 Mg·hm-2，而加入LAI后，R分别提高了11.76%和1.41%，而RMSE分别降低了11.05%和1.23%，模型预测能力明显增加。一方面，叶面积指数是反映植被固碳能力的重要参数，能够体现竹林观测结构及时间动态；另一方面，本研究采用粒子滤波同化算法对MODIS LAI时间序列数据进行同化[17,26]，使得叶面积指数年际动态变化趋势与竹林生长特征相吻合[31]，减弱了信号饱和对地上生物量估算的影响。

已有研究[7,11]认为：调查样地与遥感影像尺度匹配的差异会加大SVR模型的估算误差。为与MODIS数据匹配，本研究将30 m×30 m竹林样地进行了线性转换，但由于没有考虑混合像元的影响，样地尺度转换出现误差，并由此加大了模型估算误差。相关研究[6,25]表明：在大尺度森林地上生物量估算中利用森林丰度数据能够减少估算的不确定性。因此，本研究采用混合像元分解方法提取了竹林丰度信息，并乘以模型模拟的竹林地上生物量，进而得到浙江省竹林地上生物量时空分布，尽可能地降低了混合像元对生物量空间估算的影响。

4 结论

本研究基于MODIS LAI、EVI、RVI时间序列数据，采用支持向量回归模型结合随机森林算法，在变量筛选和模型优选的基础上，实现了浙江省竹林地上生物量较高精度的估算。结果表明：①基于随机森林模型筛选出43个对竹林地上生物量影响最大的变量，其中LAI时间序列数据的贡献最大；②基于radial核函数耦合LAI、EVI、RVI时间序列数据，构建的SVR模型训练和测试精度均最高(R分别为0.76和0.72)，均方根误差均最低(分别为5.15和8.03 Mg·hm-2)；③由SVR模型估算得到浙江省竹林地上生物量，全省竹林地上生物量均值为7.85 Mg·hm-2，总地上生物量为3.31×107Mg；浙江省竹林地上生物量在各市具有明显的差异性，其中，湖州市、杭州市、金华市、绍兴市和宁波市的竹林地上生物量均值均大于全省均值，表明其竹林经营效果相对较好。研究结果为大范围竹林碳汇遥感精准监测提供了较好的方法。