基于Sentinel-1A和Landsat 8数据的区域森林生物量反演

许振宇，李盈昌，李明阳，李 超，汪 霖

（南京林业大学 南方现代林业协同创新中心，江苏 南京 210037）

森林生物量是评价森林生态系统固碳量和碳平衡能力的重要参数之一，准确测量森林生物量对研究大面积陆地生态系统碳循环具有重要意义[1-2]。传统的野外测量是较为准确估算生物量的方法，但是此方法费时、费力，不能提供大规模连续观测数据。遥感数据可以获得植被覆盖变化、植被光谱特征等信息，而且它们与森林生物量之间关系密切。因此，各种类型的遥感系统，包括光学传感器和主动传感器，已被广泛应用于估算森林生物量[3]。

2008年9月1日，美国地质调查局（USGS）免费开放了Landsat 存档遥感数据。由于Landsat具有覆盖面广、光谱信息充分、重访频率高、计算机处理程度高、空间分辨率与森林资源连续清查样地空间尺度相近等优势，已被广泛用于区域森林生物量遥感估测。但作为一种被动式光学遥感，Landsat 影像的穿透能力不强，获取的更多是森林冠层结构信息，而且容易受云雾天气影响[4]。

合成孔径雷达（SAR）作为一种主动的遥感技术，具有穿透树冠的能力，能够全天候工作。研究表明，合成孔径雷达获取的后向散射系数与森林的结构参数有关，可用于森林地上生物量估测，但不同波长的SAR 对生物量反演具有不同的饱和点。通常，波长越短，饱和点越低，对生物量变化的响应敏感度较低，因此长波雷达数据更适合生物量估算。由于大多数长波SAR 卫星都是商业卫星，大面积数据获取成本非常高，限制了其在区域森林生物量遥感估测中的推广应用。2013年，欧洲航天局的C 波段Sentinel-1 卫星开始为全球提供免费的较高分辨率的SAR 数据。由于C 波段的穿透力和饱和点较低，单独采用Sentinel-1 数据进行遥感估测会影响估测精度。研究表明，采取主被动式遥感的数据协同策略，可以提高森林生物量遥感估测精度[5-6]。

在区域尺度森林生物量遥感估测中，除了遥感数据源外，选择合适的算法来建立生物量估计模型也非常关键。传统的统计回归方法简单易行，但这种方法不能有效地描述森林生物量与遥感数据之间的复杂非线性关系，在地形复杂的丘陵山区遥感估测精度不高。为了提高生物量模型的非线性估计能力，多种机器学习方法（如决策树、k-最近邻、人工神经网络和支持向量机）开始应用于森林生物量估测[7-8]。

通过查阅国内外文献，尚未见到基于Sentinel-1A 和Landsat 8 OLI 数据这2 种不同性质数据源的区域森林生物量遥感估测比较研究成果。本文以湖南省重点林区桂东县为例，分别采用主动式遥感（Sentinel-1A 数据）、被动式遥感（Landsat 8 OLI 数据）、主被动相结合（Sentinel-1A 数据结合Landsat 8 OLI 数据）3 种数据集和多元线性回归、随机森林、人工神经网络、袋装算法等四种模型，进行森林地上生物量（简称生物量，下同）特征变量选取、参数建模、精度评价、空间制图，以期为Sentinel-1A C 波段极化雷达数据在我国区域森林生物量的遥感估测中的推广应用提供科学参考。

1 研究区概况

桂东县隶属于湖南省郴州市，位于湖南省东南部，地处113°37′～114°14′E，25°44′～26°13′N之间，东西宽61.2 km，南北长53.6 km，总面积1 452 km2。桂东县平均海拔900 m，地势高差悬殊，呈“九山半水半分田”的土地利用格局。桂东县地处亚热带湿润季风气候区，雨量充沛，年降水量1 742.4 mm，四季分明，年平均气温15.8 ℃。

桂东县属中亚热带常绿阔叶林区，森林植被茂盛，全县森林覆盖率高达85%，是湖南省重点林业县。境内植被类型多样，且具有代表性，为野生动物提供了良好的觅食和休憩场所。此外，桂东县还是众多河流的发源地，境内有大小河流133 条，分属湘江和赣江两大水系，水资源十分丰富。

2 数据来源及预处理

论文的主要数据包括：1）来自于美国地质调查局网站（https://eartxplorer.usgs.gov/）的湖南省郴州市桂东县2014年Landsat 8 OLI 遥感影像。影像含云量为2%，空间分辨率为30 m×30 m。由于研究对象为森林生物量，所以选择的遥感影像时相为森林植被生长处于成熟期的10月；2）来自欧空局（https://scihub.copernicus.eu/）的Sentine-1A 数据，获取时间为2014年10月，空间分辨率为5 m×20 m，极化模式为VV、VH；3）地理空间数据云网站（http://www.gscloud.cn/）下载的数字地形模型（DEM），空间分辨率为30 m×30 m；4）桂东县2014年43 块森林资源连续清查固定样地数据，样地的属性表包括土地利用与覆盖、森林资源状况、生态状况三大类共60 多个调查因子。

2.1 数据预处理

2.1.1 样地森林生物量计算

论文中样地单木的地上生物量计算采用曾伟生的基于木材密度的树种组一元生物量模型[9]。单木的地上生物量方程为：

式中：Ma为地上生物量（kg）；D为林木胸径（cm）；a为参数（a=0.3p）；p为木材基本密度（g/cm3）。在计算出每块样地单木地上生物量后，再换算为每公顷森林地上生物量（t/hm2）。通过计算，43 块样地的平均生物量为52.08 t/hm2，不同级别生物量样地空间分布情况如图1所示。

图1 研究固定样地生物量空间分布Fig.1 Spatial distribution of permanent plots with biomass class in study area

2.1.2 SAR 影像及处理方法

论文所使用的合成孔径雷达数据是从Sentinel-1A 卫星C 波段传感器上获得的，中心频率为5.405 GHz。研究所用的数据为Sentinel-1A 干涉宽幅模式下的GRD 数据，采用SNAP（Sentinels application platform）软件对数据进行预处理。预处理的过程包括几何校正、辐射定标、地理编码、斑点滤波、地形校正等。通过重采样将SAR 图像空间分辨率调整到固定样地同样大小（25.82 m×25.82 m）。纹理图像使用灰度共生矩（GLCM）计算，窗口大小为3×3，灰度量化等级为64。

2.1.3 Landsat 8 OLI 影像及处理方法

Landsat 8 OLI 影像的预处理软件为ENVI。首先对图像进行辐射校正，消除和校正辐射误差和畸变，然后进行FLAASH 大气校正[10]。由于桂东县为海拔高差大、地形复杂的山区，地形校正能有效地消除地形的影响，消除地形起伏引起的阳坡和阴坡光谱特征的差异，更好地反映物体的真实光谱特征。地形校正采用基于余弦校正模型的C 校正算法。然后，对图像进行重采样，计算的纹理图像的参数与Sentinel-1A 相同。此外，论文提取了6 种常见的植被指数。

通过图像预处理，生成基于Sentinel-1A 数据的预测变量18 个，基于Landsat 8 数据的预测变量60 个。预测变量的类型、名称、个数、描述见表1。

2.1.4 遥感估测模型

论文采用多元线性回归（Multiple linear regression）、随机森林（Random forest）、神经网络（Neural network）、袋装算法（Bagging）共4 种模型进行研究区生物量遥感估测。

多元回归方法用来进行森林生物量遥感估测，是以原始波段信息、植被指数和纹理特征等为作为自变量，通过建立遥感光谱数据与样地实测森林生物量之间的相关关系，估算森林生物量。因为此方法直观易懂，且对遥感数据的处理技术要求相对较低，被众多的研究者所采用。多元回归可以定量地描述变量之间的相关性和显著性，多元回归假定响应和一组解释变量之间存在线性关系，可以表示为模型：

表1 基于Sentinel-1A 和Landsat 8 生物量预测变量†Table 1 List of biomass prediction variables based on Sentinel-1A and Landsat 8

式中，Y是生物量的值，X1，X2，…，Xn是预测变量，α0是常数，α1，α2，…，αn是与相应变量相关的回归系数，n是预测变量的个数，ε是误差项。

随机森林是Breiman[11]在2001年提出的一种基于决策树的分类回归算法。由于随机森林生成决策树的样本是随机选择的，所以随机森林模型不会出现过拟合问题，运算速率和精度都很高，对训练数据的适应能力也很强[12]。随机森林通过Bootstrapping 算法从原始样本数据集随机收集新数据集，从原始样本数据集中随机且有放回地抽取2/3 样本。随机森林用于回归的时候,将所有决策树预测结果的平均值作为最终的预测结果。随机森林模型用R 语言random Forest 包执行，它需要调整的参数为建立的决策树的数量（ntree）和决策树分裂时抽取的变量个数（mtry），本研究采用random Forest 包内嵌的tuneRF 函数来优化ntree和mtry。ntree 设置为2 500，mtry 为26。

人工神经网络应用于森林生物量遥感估测，是通过选取遥感数据有关波段的像元灰度值、植被指数等因子作为模型的输入变量，以像元所对应地理位置的森林生物量作为输出变量，利用样地数据进行训练来构建模型[13]。人工神经网络模型的优点是具有较高的精度，缺点是其模拟过程是“黑箱”操作，不能很好地解释模型的内在机理。人工神经网络选择BP 网络，是一种按误差逆传播算法训练的多层前馈网络，是目前应用最广泛的神经网络模型之一。BP 神经网络模型拓扑结构包括输入层、隐含层和输出层。BP 神经网络的学习过程包括由输入信号的正向传播和误差的逆向传播。输入信号的正向传播是指样本信号从输入层输入，再经过网络的阈值、神经元和权重的转移函数等作用后,从输出层输出的过程。若期望值与输出值之间的误差大于规定量,则进行修正。误差的逆向传播是指将误差按照梯度下降的原则分摊给各层的神经元，从而从各层神经元获得误差信号，以此作为修改权重的依据。反复进行以上两个过程，对权重进行不断地修改，这就是网络的训练过程，直到网络的输出误差满足要求或者训练次数达到设定的次数。神经网络模型用R 语言nnet 包执行，本研究采用常规的由一个输入层、一个输出层和一个隐含层组成的常规3 层BP 神经网络。本研究将隐含层单元数确定为10 个，权重调整速度为0.1，最大迭代次数为5 000 次。

Bagging（Bootstrap aggregating）[14]是集成不同的个体学习器成为一个学习器的一种算法，Bagging 算法基础是通过可重复取样得到不同的样本子集，使不同样本子集上训练得到的学习器具有较大的差异度和较高的泛化性能。袋装算法采用均匀采样从样地样本集中抽取出多个训练集，在每个训练集中重复取样，对这部分样本进行实验，分别建立决策树，再将这些决策树结合[15]。本研究采用R语言中adabag 包中的bagging 函数对样本进行建模，需要调整的参数为mfinal，表示迭代的次数，本研究将mfinal 设为20。

2.1.5 模型特征变量选择

在区域森林生物量遥感估测中，特征变量的选择是影响估测精度的一个重要因素，尤其对于随机森林、人工神经网络等机器学习算法而言。特征变量选择方法多种多样，论文主要采用逐步回归法和随机森林法两种方法进行变量选择。多元线性回归模型的变量选择采用逐步回归法，而随机森林、人工神经网络、袋装算法3 种机器学习方法模型变量的选择采用随机森林法。

逐步回归是一种线性回归模型自变量选择方法，其基本思想是将变量一个一个引入，引入的条件是其偏回归平方和检验是显著的。每引入一个新变量后都要进行F检验，对已入选回归模型的变量逐个进行t检验，将经检验认为不显著的变量删除，以保证所得自变量子集中每一个变量都是显著的。此过程经过若干步直到不能再引入新变量为止。这时回归模型中所有变量对因变量都是显著的。逐步回归在SPSS 软件中执行，步进标准使用F概率，进入和删除阈值分别设置为0.05和0.1。

利用随机森林算法，进行特征重要性度量，选择重要性高的特征。随机森林算法有两个度量变量重要性的值，可用于对变量进行排序。第一个度量值（从包外数据的排列中计算）是每棵树的预测的均方误差的百分比增加值（%IncMSE），第二个度量值是从对所有树的平均变量的拆分中节点杂质（IncNodePurity）的总减少值，用残差平方和来衡量。%IncMSE 和IncNodePurity 值越高，表明预测变量越重要。

最优变量子集的获取是一个连续的搜索过程，通常包括四个步骤：

1）子集生成：根据一定的搜索策略生成候选变量子集。本研究采用广义序贯逆向选择方法。搜索的起点是原始的全变量集。将数据集输入随机森林模型，根据度量值分别得到变量的重要性和降序。然后，移除一定数量（10%）最不重要的变量，生成变量子集。

2）子集评估：通过评估函数评估变量子集的预测性能。生成的子集被输入随机森林模型，并使用确定系数（R2）评估预测精度。

3）停止准则：确定变量搜索算法何时停止。子集求值后，应确定停止准则。如果没有停止条件，则无法停止搜索过程。本研究设置了两个停止准则：一是当子集的变量个数不大于集合个数时，对于不同的森林类型，该集合个数等于逐步回归选择的变量个数；二是当子集的预测R2连续三轮没有提高时。

4）子集验证：用于验证所选变量子集的有效性。本研究采用10 折交叉验证的方法对每一轮的可变子集的性能进行评估，因此子集验证不是独立的步骤。

本研究用R语言实现了的随机森林建模和变量选择。变量选择的工作流程如图2所示。

2.2 模型精度评价

遥感估测模型评价中经常使用到的评价指标有很多，如决定系数（R2）、均方根误差（RMSE）、相对误差（RE）、平均绝对误差（MAE）和预估精度（P）等。论文采用RMSE 和R2这2 个指标对模型预测精度进行评价：1）R2值越大，则模拟值和实际值之间的相关性越强；2）模拟值和实际值之间的RMSE 越小，则表明模型预测的效果越好。R2、RMSE 的计算方法，详见参考文献[16]，模型精度验证采用10 折交叉验证法。

3 结果与分析

3.1 特征变量选择结果

表2为分别采用逐步回归法、随机森林法进行特征变量筛选的结果。其中，多元线性回归模型的特征变量按照进入的顺序排列，而随机森林模型的特征变量按照重要性从大到小的顺序排列。

图2 随机森林模型基于变量重要性的变量选择流程Fig.2 Workflow of the variable selection based on variable importance for RF models

表2 各模型变量筛选结果†Table 2 Screening results of model variables

从表2可以看出，对于Sentinel-1A，交叉极化（VH）变量对生物量模型更重要，主要由于交叉极化对森林的垂直结构信息更为敏感[17]。在使用Landsat 8 和Landsat 结合Sentinel 时，变量筛选结果中都有红光波段（B4）和近红外波段（B5）及其纹理信息（特别是均值和相关性信息），表明这两个波段在森林生物量估算中具有十分重要的作用，这主要是红光波段处于叶绿素吸收区，对植被差异尤为敏感，而近红外波段位于水汽吸收带之间，受两个吸收到影响，对植被和土壤含水量敏感。此外，Landsat 8 生成的植被指数信息同样在森林估算中扮演重要角色。

3.2 模型精度评价

将表2中筛选的特征变量分别代入多元线性回归、随机森林、人工神经网络、袋装算法模型中，建立森林生物量遥感估测模型，并采用10 折交叉验证法进行验证。不同数据源、不同模型遥感预测精度评价结果见图3。

图3 不同回归模型预测结果的散点Fig.3 Scatter plots of predicted results of different regression models

从图3中可以看出，在4 种遥感估测模型中，无论是单一数据源还是二者结合，随机森林算法预测精度最高，人工神经网络、袋装算法次之，多元线性回归精度最低。随机森林具有较高的预测准确率，对噪声和异常值具有较高的容忍度，且很少出现过拟合。随机森林是从所有的特征中，随机选取部分特征，用来构建决策树，而袋装算法则是用所有特征来构建分类器。一般来说，随机森林的效果要好于袋装算法，偏差相对更小，方差下降幅度更大。神经网络相比传统的机器学习算法，需要的样本数量要更多，由于本研究的样地数据只有43 块，造成人工神经网络的精度低于预期。多元线性回归模型由于不能够描述生物量与遥感变量之间复杂的非线性关系，因此精度最差。从图3可以看出，这4 种模型都存在对低生物量森林过高估计，对高生物量森林过低估计的问题，这在以前的研究中同样存在，因此是一个亟须解决的问题。

从图3也可以看出，3 种不同数据源遥感估测模型，Landsat 8 与Sentinel-1A 两种数据协同估算的精度最高，Landsat 8 精度次之，Sentinel-1A 最低。由于Sentinel-1A 为C 波段SAR，波长较短，穿透能力较弱，只能反映冠层结构，而森林生物量主要分布在树干、树枝等垂直结构上；而Landsat 8 OLI 是光学被动式遥感，具有可见光、近红外、远红外波段，与植被信息紧密联系，本研究结果也证明Landsat 8 在森林生物量估算中具有至关重要的作用。Landsat 8 结合Sentinel 数据，既可以得到植被光谱信息，又可以获得植被冠层结构信息，因此两者结合在一定程度上可以提高森林生物量估算精度。

3.3 生物量空间制图

通过以上建模方式的对比，最终选择性能最优的基于主被动式遥感数据的随机森林模型，进行研究区森林生物量估算（图4）。结果显示，桂东县森林生物量平均值为53.68 t/hm2，生物量的方差为33.42 t/hm2。作为湖南省重点林业县，桂东县森林生物量平均值高于湖南省森林生物量平均值（41.274 t/hm2），但却大大低于全国森林生物量的平均值（85.64 t/hm2）[18]。

图4 2014年桂东县森林生物量空间分布Fig.4 Spatial distribution of forest biomass in Guidong county in 2014

从图4可以看出，桂东县森林生物量高值区域主要分布在海拔较高、坡度较陡的东南、西南部，而生物量较低的林分则主要集中在海拔较低、坡度平缓的低丘、河谷地带。这种森林生物量的空间分布趋势与研究区的地势地貌特征、社会经济状况相符。在桂东县，海拔高、坡度陡的地方，人口密度低、交通落后、经济不发达，森林采伐、林地被占用、森林火灾等人为干扰活动较少，森林生物量较高；而低丘、河谷地区，人口密度大、经济发达，城镇扩张、乱砍滥伐、森林病虫害等人为干扰活动频繁，导致森林生物量较低。

为进一步统计研究区生物量分布情况，将森林生物量划分为低（＜20 t/hm2）、较低（20-55 t/hm2）、较高（55～90 t/hm2）和高（＞90 t/hm2）四个等级。结果表明，桂东县森林生物量以较低（31.97%）、较高（35.02%）为主，低生物量面积比例为16.68%，而高生物量等级面积比例最低，只有16.03%。造成研究区森林生物量平均值不高、高等级面积比例较低的原因主要有2 个：一是研究区人工林比重高，中幼林面积比例大，森林单位面积蓄积量低；二是研究区系南方集体林区，作为经营主体的林农重采伐轻抚育，主动开展森林经营意愿不强、积极性不高。

4 结论与讨论

森林生物量基础研究资料不全、地面实测人力物力消耗较大，一直是进行区域尺度上森林生物量遥感估测面临的最大难题。从1978年开始，我国在全国范围内建立了以省、自治区、直辖市为总体的森林资源连续清查体系，很多县市积累了丰富的固定样地调查数据。2008年后，随着航空航天技术的发展，包括Landsat 8、Sentinel-1A在内的各种免费遥感信息源日趋丰富。因此，以少量森林资源连续清查固定样地数据、同期主被动遥感数据为主要信息源，采用多种模型进行森林生物量遥感估测，有可能为区域尺度上森林生物量动态监测探索出一条科学适用的道路。

研究表明，Landsat 红光波段（B4）、近红外波段（B5）地表反射率及纹理特征、归一化植被指数（NDVI），Sentinel-1A 交叉极化（VH）后向散射系数及其纹理特征，在森林生物量反演中具有重要作用。在4 种遥感估测模型中，随机森林算法预测精度最高，人工神经网络、袋装算法次之，多元线性回归最低。因此，在海拔相差悬殊、地形复杂、森林垂直分布明显的南方亚热带林区，可以选择随机森林作为最佳遥感估测模型。

3 种不同数据集遥感估测精度比较分析表明，基于Sentinel-1A 的估测精度最低，基于Landsat 8 OLI 的反演精度第二，两种数据协同反演的精度最高。因此，作为一种免费的主动式遥感数据，由于波长较短、林冠穿透能力弱，Sentinel-1A 极化雷达数据不宜单独用来进行区域森林生物量遥感估测，只有与被动式遥感数据结合，才能发挥其作用。

本研究中采用的C 波段Sentinel-1A 合成孔径雷达数据，穿透能力较弱，使之在森林生物量反演中作用较小。因此，在今后的研究中，我们将尝试结合穿透能力更强的P、L 波段合成孔径雷达数据进行森林生物量反演。Landsat 8 多光谱数据空间分辨率为30 m，因此其图像混合像元现象较为严重，因此，在森林生物量较低的区域，Landsat 8 像元值很容易受到灌木、草本和裸土的影响，使像元值不能正确的反映该像元实际的生物量信息；而在森林生物量较高的区域，则由于Landsat 8 生物量估算的饱和问题，同样不能正确的反映该像元实际的生物量信息。在今后的研究中，我们将尝试采用混合像元分解方法、高空间分辨率或高光谱遥感影像解决混合像元和光谱饱和问题。此外，尽管国家森林资源连续清查数据是目前可用的质量最高的区域性森林调查数据，但其空间间距较大，导致县域尺度空间密度较小。因此，可以尝试在国家森林资源连续清查的基础上，适当补测调查样地，提高样地空间密度。

尽管在预测变量的选取中综合考虑了波段反射率、植被指数、纹理特征，然而，无论采用何种数据源，4 种估测模型的预测精度均低于拟合精度。以预测精度最高的随机森林为例，在3 种数据源中，拟合精度R2均超过0.9，但是预测精度R2均小于0.8，说明在森林参数反演中普遍存在着过拟合现象。这种现象的产生可能与研究区复杂的地形地势条件、亚热带常绿阔叶林复杂的林种、树种结构有关，而导致遥感估算森林生物量存在饱和现象，这也是在今后的研究中需要重点解决的问题。

此外，在预测变量的选择上，本研究没有考虑地形、气候等因素的影响。如何结合地形因素、气候等因素，分树种（组）、郁闭度或样地生物量分段建模，以提高森林生物量遥感估测的精度，也是今后研究的重要内容。