基于多时相光学和雷达遥感的太平湖生态保护区森林地上生物量反演

卢 佶，张国威，吴 昊

（国家林业和草原局 华东调查规划院，浙江 杭州 310019）

森林是陆地生态系统的主体，主导陆地碳循环的动态变化，在全球碳收支中发挥着重要作用[1-2]。森林碳储量的动态变化可以通过定期绘制森林地上生物量进行监测[3]。基于样地的森林调查方法可以提供精确的森林地上生物量。然而，在大尺度监测中存在局限性，例如调查周期较长、时间分辨率不高等问题[4]。因此，寻求大面积森林生物量和碳储量精准监测的方法显得尤为重要。

遥感技术为大面积反演森林地上生物量(AGB）或碳储量提供了全覆盖和长时间序列数据。目前，已有研究使用较高分辨率(＞30 m）的光学数据(Sentinel-2）[5]和合成孔径雷达(Sentinel-1）数据[6](25 m）进行大尺度森林地上生物量反演。但现有研究对地上生物量的估计往往出现相互矛盾和不一致的现象，可能是缺乏大量的地面数据参考；也可能是森林结构复杂存在空间异质性，遥感信号响应机制不清晰等所致[7-8]。尤其对于热带或者亚热带森林地区，高郁闭特征和异质性需要探讨光学和合成孔径雷达数据的互补机制，以提高热带森林地上生物量反演的准确性和一致性。

近年来，已有研究通过单个传感器或者联合光学雷达，来评估Sentinel 数据对不同生物群落的森林地上生物量反演潜力[9]。如PANDIT 等[10]采用Sentinel-2 数据及其衍生的光谱指数，发现与光谱波段相比，红边植被指数对地上生物量估计具有更高的贡献度。也有研究表明：Sentinel-1 数据在地上生物量和碳储量反演方面也有较大的优势[11]。周蔚等[12]对亚热带森林地上生物量建模中，发现Sentinel-2 的光谱指数比光谱波段贡献度更大，并且联合Sentinel-1 和Sentinel-2 数据反演精度最佳。CHEN 等[13]采用Sentinel-1 和Sentinel-2 及其衍生数据(如纹理、光谱指数、生物物理变量），系统比较了1 种参数化(地理图形加权回归）和3 种非参数化的机器学习算法(支持向量机、随机森林和人工神经网络）的地上生物量建模方法，结果表明：Sentinel-1 的纹理特性和Sentinel-2 的生物物理变量是地上生物量估计模型的最优变量。LIU 等[14]和CHEN 等[15]结合了Sentinel-1、Sentinel-2 和航天飞机雷达地形测绘任务(SRTM）的数字高程模型(DEM）数据，验证了Sentinel-1 和Sentinel-2 对森林地上生物量绘图的适用性，并且发现随机森林算法的精度最优。然而，以上研究大都是基于生长季时期的研究，鲜有采用不同季节的光学和雷达数据对森林地上生物量反演进行探究，尤其是亚热带地区夏季降水多，森林湿度较大，雷达对地表湿度比生物量更加敏感。鉴于此，本研究以太平湖生态保护区森林为例，采用随机森林算法评估不同季节Sentinel-1 和Sentinel-2 数据的光学和雷达遥感信息对地上生物量的响应特征，以期为亚热带森林碳储量精准评估提供科学参考。

1 研究地区与方法

1.1 研究区概况

本研究区位于黄山市太平湖生态保护区，该区地处亚热带季风湿润气候，年均降水量达1 617.0 mm，年均气温为15.5 ℃，年日照时数＞1 600.0 h，年均蒸发量＞1 100.0 mm。研究区森林覆盖面积约9 万hm2，平均海拔为220 m，高差为750 m，平均坡度为22°。生态保护区植物资源丰富，乔、灌树种900 多种。研究区为亚热带针阔混交林，以中幼林龄为主，常见主要树种包括马尾松Pinusmassoniana、栎树Quercus、木荷Schimasuperba、桦树Betula、楝Meliaazedarach、樟树Cinnamomumcamphora、杉木Cunninghamialanceolata、毛竹Phyllostachysedulis等。

1.2 样地数据

62 个样地数据来源于2019 年生态保护区一类清查数据，调查时间为2019 年8—9 月。每个样地大小为25.82 m×25.82 m，样地调查因子包括树种、株数密度、胸径、树高等。每个样地的地上生物量计算采用一元地上生物量方程：

式(1）中：Ma为地上生物量；D为胸径；a和b为不同树种系数。不同树种系数取值参考周蔚等[12]的研究。每个样地树木地上生物量累计值为样地尺度的地上生物量。样地参数统计如表1 所示。

表1 研究区样地参数统计Table 1 Statistics of forest plot parameters in the study area

1.3 遥感数据获取和预处理

考虑覆盖各个生长季节，选择2019 年1、3、6、9、12 月5 期Sentinel-1 和Sentinel-2 影像数据，其中Sentinel-2 影像云量均＜10%。卫星数据从欧空局Sentinel 数据中心(https：// scihub.copernicus.eu/dhus/#/home）下载。

1.3.1 Sentinel-1 数据和预处理 Sentinel-1 数据选择干涉宽幅模式(interferometric wide swath，IW），具有双极化(VV 和VH）的升轨数据。该产品类型是多视影像(GRD），空间分辨率为10 m，采集入射角为26°～40°。本研究共采集了5 景不同季节影像，并根据研究区范围进行裁剪。对每景Sentinel-1 影像数据预处理，采用欧空局官方公布的哨兵应用平台软件(SNAP）进行精确轨道校正、地形辐射校正(归一化）和斑点过滤(图1）。其中，采用SRTM 数字高程模型(DEM）用于地形校正，地形辐射校正后的图像采用Lee 斑点过滤算法以减少斑点。为了避免后向散射系数值较小对后续计算和分析造成的影响，将其转化成归一化值[16]。

式(2）中：dB为归一化后的后向散射系数，σ0为VH 或者VV 极化的后向散射系数。此外，根据LAURIN 等[16]的研究，采用dB的差值(VH-VV）与求和(VH+VV）有助于提高地上生物量的建模精度。

1.3.2 Sentinel-2 数据和预处理 本研究选择云量小于5%的Sentinel-2 影像，共采集10 张影像覆盖研究区域。对于每景影像，使用SNAP 中的Sen2cor 插件进行大气校正，将L1C (大气顶部TOA）产品转换为L2A (大气底部BOA）。大气校正过后的Sentinel-2 影像，将3 个可见光、3 个红边、2 个近红外、2 个短波红外(B11 和B12）共10 个波段数据的空间分辨率重采样为20 m，用于后续的预处理和分析。采用SNAP 中的“Biophysical Processor”提取叶面积指数(LAI）、植被覆盖率(fcover）和光合有效辐射吸收比(fapar） 3 个生物物理参数，描述植被状态和空间动态分布(图2）。这3 个生物物理参数是基于PROSAIL 辐射传输模型反演获得[17]。主要计算步骤包括：①标准化输入。将10 个波段的BOA 数据输入模型；②人工神经网络(ANN）算法反演。基于PROSAIL 模型构建查找表并采用ANN 训练输入参数和模拟光谱的非线性关系；③标准化输出。根据Sentinel-2 观测值输出选择的目标参数；④生成质量评估指标。评估反演值的不确定性。除3 个生物物理参数外，还基于Sentinel-1 和Sentinel-2 影像，提取了光谱波段、光谱指数、生物物理参数和后向散射系数用于地上生物量建模的遥感特征参数(表2）。

表2 地上生物量建模的Sentinel 特征参数Table 2 Sentinel feature parameters for above-ground biomass modelling

1.4 随机森林回归模型

对研究区亚热带森林地上生物量估算采用随机森林回归方法。随机森林算法通过生成大量的决策树，每一个决策树都是使用训练数据的唯一引导样本独立构建，对所有决策树预测器的平均值作为最终目标变量的估计值[18]。随机森林根据所有可用的预测器随机选择，减少决策树之间的相关性，从而减少噪声和提高预测的准确性[19]。随机森林算法对数据冗余和非线性处理具有很好的鲁棒性，能够处理具有不同属性和数值范围差异较大的变量，避免产生较大的过拟合现象[20]。

随机森林回归模型采用R 语言的“randomForest”包(适用于分类和回归训练）实现。该软件包通过给定的模型参数范围进行迭代，获得模型最佳性能的参数组合。决策树个数(ntree）和每个节点上随机选择的预测变量数量(mtry）采用网格搜索法优化设定。

1.5 模型构建和精度评价

为了比较时间序列上Sentinel-1 和Sentinel-2 数据对地上生物量的估算结果，建立了6 组回归模型(表3）。在进行每个回归模型训练前，采用方差膨胀因子(VIF＜10）排除共线性强的变量，减少模型过拟合问题。此外，采用十折交叉验证的方法进行模型训练和验证。最后，利用决定系数(R²）、均方根误差(ERMS）对每个模型的地上生物量估算精度进行评价。

表3 基于Sentinel-1 和Sentinel-2 的地上生物量模型构建Table 3 Above-ground biomass modelling based on Sentinel-1 and Sentinel-2

式(3）～(4）中：Oi为地上生物量观测值；Pi为地上生物量预测值；为地上生物量样地观测平均值；i为样地号；n为样地总数量。

2 结果与分析

2.1 森林冠层光学和雷达特征参数的季节变化特征

如图3 所示： 在近红外波段，LAI 对冠层反射率的贡献占主导地位，森林地上生物量越大LAI 越大。因此，对于相同季节，近红外波段的森林冠层反射率随森林地上生物量增大而增大。然而，在短波红外，森林冠层主要受水分影响。因此，在短波红外的森林冠层反射率随森林地上生物量的增大而减小。总体而言，6 和9 月的森林冠层反射率较大，并且不同地上生物量水平间的冠层反射率区分度较大；相反，1、3 和12 月不同地上生物量水平间的森林冠层反射率差异性减小，主要是由于在这些月份不同地上生物量水平的森林叶面积差异不显著。

图3 在不同等级地上生物量下森林冠层反射光谱曲线随季节的变化趋势Figure 3 Seasonal variation of forest canopy reflectance spectra under different above-ground biomass levels

从图4 可见：增强植被指数(EVI）和归一化植被指数(NDVI）在相同季节随着地上生物量等级升高而增大；相反，归一化植被水分指数(NDWI）在相同季节随着地上生物量等级升高而降低。此外，地上生物量在较低等级(0～150 Mg·hm-2）范围内，相同季节的EVI 和NDVI 指数变化难以区分；而在地上生物量较高等级(＞150 Mg·hm-2）范围内，NDVI 对地上生物量的敏感度要优于EVI 指数。相似地，在地上生物量较低等级范围内，植被水分敏感指数Ⅰ(STVI1）和植被水分敏感指数Ⅱ(STVI2）在相同季节性差异不大，但在地上生物量较高等级范围差异相对明显。与光谱指数相似，fapar、fcover、LAI 在不同季节变化趋势一致，呈现先增高后降低的趋势(图5）。当地上生物量在0～100 Mg·hm-2与＞100 Mg·hm-2范围内时，3 个生物物理参数在相同季节差异不明显。

图4 不同等级地上生物量的植被指数季节变化曲线Figure 4 Seasonal change curves of vegetation index in different levels of above-ground biomass

图5 不同等级地上生物量下生物物理参数的季节变化Figure 5 Seasonal variation of biophysical parameters in different levels of above-ground biomass

相对于Sentinel-2 的光学特征参数，Sentinel-1 的后向散射系数对研究区的针阔混交林在不同季节的差异性表现不明显(图6）。其中，对较高等级的地上生物量，VH、VV 和VH+VV 在1 和12 月的值要低于3—9 月，表现出旱季低雨季大的规律。但是，对于较低等级的地上生物量，这3 个指标在不同季节没有明显的差异。总体而言，相对于较高等级的地上生物量，Sentinel-1 的后向散射特征要比较低等级的地上生物量高。

图6 不同等级地上生物量下雷达特征参数的季节变化Figure 6 Seasonal variation of SAR characteristic parameters for different levels of above-ground biomass

2.2 基于不同季节Sentinel-1 和Sentinel-2 影像的地上生物量反演结果

如图7 所示：基于Sentinel-2 影像的模型(S2-all、S2-biophysical 和S2-spectral）要优于Sentinel-1 影像的模型(S1-all 和S1-backscatter）。对采用不同季节Sentinel-1 后向散射系数的S1-backscatter 模型，地上生物量反演精度最低，R2为0.36，ERMS为52.97 Mg·hm-2。同时，对较高等级的地上生物量(＞150 Mg·hm-2）同样存在较严重的低估现象。而对于S1-all 模型，加入了后向散射系数相关的雷达衍生指数(VH-VV 和VH+VV）后估计精度有所提高(R2=0.47，ERMS=49.11 Mg·hm-2），但仍然无法避免对较高等级地上生物量的低估现象。

图7 地上生物量实测值与估计值散点图Figure 7 Scatter plots of measured and estimated values of above-ground biomass

基于生物物理参数的S2-biophysical 模型对地上生物量估计精度(R2=0.66，ERMS=40.12 Mg·hm-2）优于基于光谱波段的S2-spectral 模型(R2=0.60，ERMS=42.95 Mg·hm-2）。并且，仅采用光谱波段的S2-spectral 模型对较高等级的地上生物量(＞150 Mg·hm-2）存在低估现象。S2-all 模型联合光谱波段、光谱指数及其衍生的生物物理参数后，地上生物量估计精度提高了5.62 Mg·hm-2。通过联合不同季节所有的Sentinel-1 和Sentinel-2 特征变量所构建的模型精度最高(R2=0.78，ERMS=31.56 Mg·hm-2），相比S2-all 模型，地上生物量均方根误差降低了8.56 Mg·hm-2。

2.3 地上生物量建模的最佳预测因子和影像采集季节

如图8 所示：在S1-backscatter (图8A）的预测变量中，排前5 位的变量分别为12 月_VH、6 月_VH、6 月_VV、1 月_VH 和12 月_VV，采集季节都来自夏季和冬季，其中贡献度最大的是12 月_VH，这一规律和S1-all (图8B）相同。在S1-backscatter 和S1-all 模型中，VH 极化的后向散射系数比VV 贡献度更大。

图8 地上生物量模型变量的重要性Figure 8 Variable importance of above-ground biomass modelling

S2-spectral 模型前15 个最重要预测变量中(图8C），有12 个变量来源于生长季节，再次证明了基于生长季的光学影像对常绿针阔混交林地上生物量估计的重要性。S2-biophysical 模型(图8D）变量中，3 个生物物理参数中LAI 贡献度最大，其次是fapar 和fcover。与光谱波段相比，S2-all 模型(图8E）中生物物理参数LAI、fapar 和STVI2 对地上生物量估计贡献度较高，并且这些参数采集的影像时期大部分来源于6 和9 月。结合所有变量的S1S2-all (图8F）模型中，前15 预测变量几乎全为Sentinel-2 变量，表明在对亚热带森林地上生物量反演方面，Sentinel-2 比Senttinel-1 更有优越性。

3 讨论

3.1 基于Sentinel-1 与Sentinel-2 数据的地上生物量建模比较

本研究发现：Sentinel-2 比Sentinel-1 更适合于亚热带太平湖生态保护区森林的地上生物量建模，这与VAFAEI 等[21]和ZHAO 等[22]的研究结果一致，但其他相关研究却获得了相反的结果[23]。本研究中Sentinel-1 地上生物量较差的估算结果可以归纳为2 个主要因素。第一个，Sentinel-1 的C-波段波长较短是主要限制因素。与波长较长的雷达波段(如L 和P）相比，C 波段雷达穿透植被冠层的能力有限，无法获取完整的森林结构信息，因此，采用较长波长(如L 和P 波段）的Sentinel-1 数据和Sentinel-2 联合，可以提高地上生物量反演的准确性[16]。第2 个潜在的因素是本研究没有考虑Sentinel-1 的纹理信息。有研究表明：Sentinel-1 的后向散射系数及其衍生参数(如：VH+VV 或VH-VV）的纹理信息对地上生物量制图具有较大的贡献度[23]。因此，在未来的研究中，应考虑将纹理信息结合起来，以改善地上生物量估计。

光学数据在绘制地上生物量估计方面容易存在饱和现象[16]，但对本研究森林的地上生物量反演影响不明显。一方面，Sentinel-2 具备较高的空间分辨率，有助于提高地上生物量的反演能力。另一方面，在本研究中，生物物理参数(fapar、fcover 和LAI）、植被指数和红边谱段对地上生物量估计具有重要作用。植被指数和光谱波段在以往的研究中已经被证实其有效性，但Sentinel-2 反演的生物物理参数对地上生物量建模的重要性和贡献性少有探究。在本研究中，仅使用生物物理参数(S2-biophysical，ERMS=37.69）比使用所有光谱波段(S2-spectral，ERMS= 42.95）得出的精度更高。并且，它在地上生物量估计中的表现优于所有Sentinel-1 数据，证明了生物物理参数对亚热带针阔混交林地上生物量估计的优势。

在所有的地上生物量模型中，结合Sentinel-1 和Sentinel-2 数据取得了最好的结果，其他相关研究也发现了类似的结论[24-25]。主要由于2 个系统的成像技术、数据特征和信息内容不同，两者优势互补能够有效提高地上生物量估计。在本研究中，使用年内的时间序列数据，包含不同季节进行组合，对Sentinel-1和Sentine-2 的互补性有了更深层次的探索。与周蔚等[12]仅采用1 个时间数据结果(ERMS=37.05 Mg·hm-2）相比，本研究结果(ERMS=31.56 Mg·hm-2）有一定程度的提高，这说明不同时间序列的Sentinel 数据越多，对地上生物量建模的准确性越好。

3.2 基于Sentinel 数据地上生物量建模的重要性变量和采集季节

本研究显示：生物物理参数LAI、fapar 和fcover 的贡献度比植被指数、光谱波段、雷达极化参数及其衍生指数大。CHEN 等[26]获得了类似的结果，发现LAI 和地上生物量之间的相关性更高。值得注意的是，不同生物物理变量的重要性与植被组成和特征关系紧密。本研究区森林比较郁闭，LAI 对地上生物量估计起到关键作用。然而，如果是稀疏的草原林地，林下灌草茂密对整个林分LAI 贡献较大，相反fapar 参数将起到较大的贡献作用。此外，对于仅考虑Sentinel-1 的模型，VH 极化和VH+VV 参数具有较高的重要性，这也验证了交叉极化(HV，VH）比同极化(HH，VV）更适用于地上生物量估计。

S2-all、S2-biophysical 和S2-spectral 模型的变量重要性图显示：生长季(6—9 月）对研究区的地上生物量建模很重要。这主要是由于大多数植被类型在该季节有较高的植被覆盖率。同时，涉及Sentinel-1数据的模型(S1S2-all、S1-all 和S1-backscatter）在旱季(10 月至翌年5 月）获取的图像对研究区地上生物量的建模也很重要。其他相关研究也表明：Sentinel-1 数据对植被冠层、水含量、土壤水分等季节性变化产生差异[27-28]。其主要原因是较短波长的X 波段和C 波段对地表水分比生物量更加敏感[29]。此外，Sentinel-1 对树冠穿透力有限，对稀疏的森林地上生物量估计是有利的，尤其是对落叶森林。因为在旱季(落叶期），树枝裸露容易被观测。通过结合不同季节(旱季和雨季），更有利于联合光学和雷达数据的互补性，提高地上生物量建模的准确性。

4 结论

Sentinel-2 模型对地上生物量反演的精度要优于Sentinel-1，联合两者的反演精度最大。其中，生物物理参数(LAI、fapar 和fcover 指数）贡献度最大，植被指数(STVI1 和STVI2）和后向散射系数(VH 极化）贡献度次之，光谱波段最小。再次证明了生长季对亚热带森林地上生物量反演的重要性。同时，也证明了联合生长季的Sentinel-2 光学影像与旱季的Sentinel-1 雷达数据对提高森林地上生物量具有较大潜力。