Estimating bamboo forest aboveground biomass based on MODIS time series data

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

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]。

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

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

设置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)波段分别计算。

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

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

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

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

其中：N为样本个数，i=1, 2$, \;\cdots , $ N；y_m表示模型模拟值，y_o表示实测值；$ {\bar{y}}_{{\rm{o}}} $为N个实测值的平均值。

基于随机森林算法得到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模型输入变量，构建竹林地上生物量估算模型。

Figure 1.  Importance score of input variables on aboveground biomass of bamboo forest

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

Figure 2.  Estimating accuracy of 4 kernel functions in SVR model for aboveground biomass of bamboo forest

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

Figure 3.  Model accuracy of LAI, EVI, RVI and EVI+RVI time series data for estimating aboveground biomass of bamboo forest using SVR model

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

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

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

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

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