基于Sentinel-1 和Sentinel-2A 数据的森林蓄积量估算

雷令婷，高金萍，张晓丽，高显连

(1.北京林业大学 林学院，北京 100083；2.国家林业和草原局 调查规划设计院，北京 100714)

森林蓄积量反映一个国家森林资源规模，是中国森林资源调查的重要内容[1-2]，也是衡量森林质量、森林固碳作用[3-4]以及评价森林经营效果的关键性指标[5].遥感技术对于高效准确地获取森林蓄积量、进行森林资源调查和监测等具有重要意义.如何利用遥感技术结合少量地面实测数据进行蓄积量估算[6]，减少野外调查工作量，提高调查的精度和时效性，是目前森林资源调查的重点问题[7].

光学遥感数据凭借多时相、多分辨率、光谱信息丰富、获取方便等优势[8]，已被广泛应用于森林蓄积量估算[9-11].基于光学遥感数据的蓄积量估算方法主要通过遥感因子结合地面实测数据，利用多元线性回归[12]、随机森林[13]、k-近邻算法（k-Nearest Neighbor,k-NN）[14]等方法建立反演模型.光学影像提供丰富的冠层信息，但因波长较短、对森林空间结构信息不敏感[15]、无法穿透树冠且光学信号易饱和[16-17]等局限性影响蓄积量估算精度.合成孔径雷达（Synthetic Aperture Radar,SAR）等获取微波遥感数据的方法森林穿透性较好，更好地刻画森林垂直结构特征、森林密度等信息[18-19]，且不受天气影响[20]等，在森林蓄积量估算方面具有一定优势.近年来，学者们利用SAR 数据开展了较多森林蓄积量估算研究[21].然而，SAR 数据无法提供与植被生理相关的光谱信息，总体的估算精度并不高.此外，由于SAR 数据还受地形因素影响[22]，一定程度上影响了蓄积量估算效果.仅利用光学影像提取的森林水平结构信息或SAR 垂直结构信息进行蓄积量估算均具有局限性.若将光学和SAR 数据结合估算森林蓄积量，可以优势互补，具有提高估算精度的潜力.因此，研究者们开展大量研究挖掘协同森林水平与垂直结构信息进行森林蓄积量估算的方法[23].总之，协同主动、被动遥感数据进行蓄积量估算，精度具有一定程度的提高.然而，目前国内外利用光学遥感影像与SAR 数据结合的研究多集中于生物量估算和树高估算中，对森林蓄积量估算研究相对较少.

在吉林省临江市西小山林场内，本文以微波遥感数据的Sentinel-1 和光学数据的Sentinel-2A 为遥感数据源，通过多元逐步回归和随机森林方法构建森林蓄积量估算模型，并利用地面实测数据进行精度验证，确定森林蓄积量估算最佳模型，同时分析各变量对森林蓄积量估算模型的贡献，最终获取西小山林场森林蓄积量的空间分布特征，可以实现森林蓄积量高效精准调查.

1 研究区与数据

1.1 研究区概况研究区位于吉林省临江市西小山林场（127°18′8″～127°32′52″N，41°45′46″～41°36′25″E），林场总面积为218.869 km2，地处长白山腹地，鸭绿江畔.该区属于温带大陆性季风气候，年平均气温为2～4 ℃，年平均降水量为750～1 000 mm.西小山林场以人工林为主，树种结构复杂，包括针叶林和阔叶林，其中针叶树由云杉（Picea asperata Mast）、红松（Pseudotsuga brevifolia Cheng et L.K.Fu）、臭松（Abies sibirica）、樟子松（Pinus sylvestris）等组成；阔叶树包括白桦（Betula platyphylla）、大青杨（Populus ussuriensis Kom）和榆树（Ulmus pumila L）等.森林覆盖率达66.7%，成熟林面积达133.220 km2，约占林场总面积的61.11%，平均树高约为13.3 m.

1.2 数据收集与处理

1.2.1 地面调查数据 地面调查数据主要为样地数据和小班数据，均来源于2016—2018 年国家林业和草原局开展的东北内蒙古重点国有林区二类调查.其中，研究区按照林场面积计算该林场内布设样地的数量，从而确定抽样网格密度，每个网格的中心作为样地的位置.最终确定共布设56 块样地（东西为30 m，南北为20 m）（图1），记录了环境因子数据（坡位、坡向和海拔等）、每木检尺数据（胸径、树高等）以及郁闭度、林龄和株数等信息.

图1 研究区位置图Fig.1 Location of study area

样地蓄积量计算是根据每木检尺结果，计算样地林木株数、胸径断面积，采用立木材积表[24]分树种计算各林木材积，各树种材积之和得到样地材积，样地材积除以样地面积得到每公顷蓄积量.

小班数据主要包括小班面积、平均胸径、平均蓄积量和树种类型等，主要用于生成林场小班蓄积量分布图.

1.2.2 影像数据 本研究影像采用欧洲“哥白尼计划”中的Sentinel 系列卫星数据，分别是Sentinel-1 和Sentinel-2.其中Sentinel-1 SAR 是一个基于C波段的成像系统，分辨率最高5 m，幅宽达到400 km[25].Sentinel-2 卫星携载MSI 多光谱成像仪，覆盖13 个波段[26]，观测幅宽达到290 km，空间分辨率分别为10、20 m 和60 m，重返周期为5 d[27].通过美国ASF 网站（https://vertex.daac.asf.alaska.edu/）下载2017 年9 月15 日sentinel-1 干涉宽模式（Interferometric Wide swath,IW）S1B GRDH 影像，极化方式为 VV、VH，空间分辨率为 5 m×20 m；在USGS 网站（https://earthexplorer.usgs.gov/）下载2017年9 月23 日Sentinel-2A 的 Level-1C 数据.由于SAR 影像数据存在因地形、斑点噪声引起的几何和辐射畸变，为了能够实现地面实测数据和SAR影像精确匹配，对数据进行轨道校正、辐射定标、多视化、地形校正、分贝化和有效散射面积校正以及重采样（采样大小为10 m）等处理，得到后向散射系数.由于Sentinel-2A Level-1C 数据是经过几何精校正的正射影像，并没有进行辐射定标和大气校正，因此利用Sen2Cor 大气校正处理器（2.5.5 版本）进行大气校正得到Level-2A 大气底层反射率数据，同时采用C 模型进行地形校正，并将所有波段影响重采样到10 m.

2 研究方法

2.1 纹理变量纹理信息可以反映森林冠层信息，避免由于树木相互遮挡造成的影像空间色调变化，在一定程度上反映森林蓄积量信息[3].本研究利用灰度共生矩阵提取Sentinel-2A 波段的纹理信息，窗口大小设置为3×3.具体见表1.

表1 Sentinel-1 纹理信息列表Tab.1 The list of texture information from Sentinel-1

2.2 光学变量由于光学影像的波段光谱信息及其相关植被指数能够反映植被生长状况，可以作为反映森林蓄积量的直接变量.本研究从Sentinel-2A 影像中提取10 个波段的光谱信息和14 个植被指数，具体见表2.

表2 Sentinel-2A 特征变量列表Tab.2 The list of variables derived from Sentinel-2A

2.3 SAR 的投影角校正的后向散射系数合成孔径雷达（SAR）数据对于森林垂直结构信息较为敏感，在一定程度上可以反映森林蓄积量状况.然而，由于SAR 传感器侧视成像的特点，后向散射特征受到局部地形起伏的影响.因此，我们对SAR 数据进行投影角校正，并提取VV 极化方式下的变量VH 极化方式下的变量和VH/VV 的变量

2.4 森林蓄积量模型本研究采用多元逐步回归方法，分别以提取的纹理信息、光学变量和后向散射系数以及不同类型变量组合为自变量，以实测样地蓄积量为因变量，构建森林蓄积量估算模型.建模过程中，每个自变量随着对回归方程贡献的变化，将自变量引入或剔除出回归方程，最终在回归方程中的自变量均为显著的变量.此外，利用随机森林方法构建森林蓄积量估算模型.随机森林的流程主要为：①输入样本集；②随机选择训练集；③从M个特征变量中选择m个特征变量（m＜M）；④基于以上步骤，构建n个决策树并求取平均值从而得到最终预测结果.该方法可以克服回归方程可能造成的过拟合问题，其建立多个决策树，并将其融合从而得到一个较为准确的模型[1].

在进行样地调查时，由于有4 个样地中蓄积量为零，即在这4 个样地中没有树木予以剔除，以排除蓄积量为零对建模产生影响，从而导致蓄积量估算模型不准确.

为了验证Sentinel-1 以及Sentinel-2A 影像在森林蓄积量估算方面的潜力并分析模型的估算精度，利用留一交叉验证方法对5 个模型结果分别进行精度评价，并通过决定系数（coefficient of determination，R2）、留一交叉验证均方根误差(Root Mean Square Error of leave-one-out Cross-Validation,RMSECV)和相对留一交叉验证均方根误差(relative Root Mean Square Error of leave-one-out Cross-Validation，rRMSECV)验证模型的预测能力.

3 结果与分析

3.1 特征变量筛选将基于Sentinel-1 和Sentinel-2A 所提取的因子分成4 组建立森林蓄积量估算模型：①仅有纹理信息；②仅有光学特征变量；③仅有后向散射系数；④3 组变量组合；对每组均采用逐步回归方法筛选出森林蓄积量估算的最佳变量（表3），进而完成对森林蓄积量的估算.

4 组模型特征变量筛选结果如表3 所示，R2介于0.115～0.513，不同类型变量建立的森林蓄积量估算效果差异较大.其中，在纹理变量中筛选出CorB4、VarB5和VarB11与森林蓄积量相关性较高；在光学特征变量中筛选出B2、B11、PSSRa和REIP等4 个变量相关性较高；在后向散射系数中，等3 个变量相关性均较好.利用3 类变量的模型（d）具有较好的表现，R2为0.513，仅使用单独纹理信息、光学特征变量以及后向散射系数等分别建立的模型，R2较低.

表3 特征变量筛选结果Tab.3 The results of feature variable screening

3.2 森林蓄积量建模本文将筛选的变量分别利用多元逐步回归和随机森林方法构建模型，如表4所示.结合3 种类型变量建立森林蓄积量模型（d）表现出明显的优势，其中模型（d）的RMSECV为49.70 m3·hm-2，rRMSECV为0.26.而仅利用单一类型的变量进行森林蓄积量估算时，表现均低于模型（d），RMSECV相对较高为53.91 m3·hm-2，其中基于纹理信息模型表现相对较好.基于后向散射系数模型表现最差，RMSECV达到57.78 m3·hm-2，无法进行森林蓄积量估算.利用随机森林的模型5 表现相对于线性模型较差，RMSECV达67.45 m3·hm-2.

表4 模型精度评价表Tab.4 The estimation accuracy of different models

利用留一交叉验证对5 个模型估算性能进行了验证，建立预测蓄积量与实测蓄积量的散点图（图2）.从图2 可以看出，模型（d）表现相对较好，散点均匀分布在两侧，表明这个模型具有一定的适用性.而单一类型变量构建模型均相对较差，RMSECV均较高.此外，利用SAR 数据提取变量建立的模型（c）（图2（c））由于变量信息较少，无法较好估算森林蓄积量，RMSECV相对较高.综上说明，模型（d）相较于单一变量模型和随机森林模型更适合研究森林蓄积量估算与区域制图.

图2 不同模型预测与实测蓄积量散点图Fig.2 Scatter plot of predicted and measured stock volume from different models.

3.3 森林蓄积量空间分布比较分析5 种森林蓄积量模型可知，结合纹理变量、光学特征变量和后向散射系数构建的模型（d）在研究区森林蓄积量估算精度最高.因此，本研究以小班为单位，利用模型（d）对林场森林蓄积量进行反演，结果如图3（a）所示，并以小班蓄积量实测图进行对比（图3（b））.从图3 可以看出，西小山林场蓄积量呈现中南部较低，东部和中部偏西蓄积量较高的空间格局.小班实测平均蓄积量为147.14 m3·hm-2，预测平均蓄积量为188.31 m3·hm-2，其中，蓄积量在100～200 m3·hm-2分布最多，约占总面积的41.96%，蓄积量0～100 m3·hm-2和300～400 m3·hm-2均较少，分别占总面积的10.44%和7.55%，在预测过程中出现蓄积量值高估的现象，高估的原因可能是建模样本数量不够多.

图3 西小山林场蓄积量反演图和小班蓄积量实测图Fig.3 The inversion map of stock volume and the measure map of stock volume in subcompartment of Xixiaoshan Forest Farm

4 讨论与结论

4.1 讨论本研究将Sentinel-1 和Sentinel-2A 提取的纹理信息、光学特征变量和后向散射系数作为特征变量，通过逐步回归方法最终筛选出VarB5、VarB11、B11、PSSRa、REIP、作为建模变量，利用多元逐步回归和随机森林方法分别对森林蓄积量进行估算.在纹理变量中方差（VarB5和VarB11）对森林蓄积量的贡献较大，主要反映给定窗口内光谱异质程度和像元间的相似程度，是影像灰度统计和空间特征信息的有效反映，可以有效地减少影像中阴影对森林参数提取的影响，同时可以提高森林与其他地物的差别[28]，从而使得纹理变量与森林蓄积量的相关性提高；短波红外的B11 波段为蓄积量估算模型的第一贡献者，主要原因可能为短波红外对生物量较为敏感，而生物量与蓄积量存在一定关系[29]；PSSRa 和REIP 对于蓄积量估算的贡献较大，主要原因可能是Sentinel-2A 数据中包含3 个与植被各类理化参数相关的红边波段[30]，而PSSRa和REIP 两个参数是由红边波段构成[31]，可以反映森林蓄积量的变化；SAR 数据后向散射系数主要由于SAR 数据对于树高等因子较为敏感[32-33]，而树高与森林蓄积量关系密切[34-35]，因此也作为森林蓄积量建模的变量之一.

筛选特征变量时，利用纹理信息、光学特征变量、投影角[36]的后向散射系数单一特征变量建立模型R2较低，而基于3 类变量建立的模型（d）精度较高，这是由于森林蓄积量对森林水平结构和垂直结构信息均具有较强的依赖性，仅使用单独结构的变量不足以充分描绘.此外，不同森林结构的特征变量为森林蓄积量估算提供不同权重信息，因此需要进一步明确各个结构特征变量对森林蓄积量的估算能力的影响，从而提高森林蓄积量估算精度.

对比5 个模型发现，模型（a）筛选的纹理变量能反映影像中森林纹理信息的分布和变化情况；模型（b）筛选出的变量能够较为全面地代表植被信息，起到图像增强的作用；模型（c）变量仅反映森林垂直信息.总之，以上模型筛选出的变量较为单一.利用不同类型的变量建立的模型（d），在一定程度上可以较为完整地反映森林蓄积量信息，同时较好地抑制了仅利用光学遥感和微波遥感所造成的蓄积量估算不准确的问题.基于光学和微波遥感数据利用随机森林方法进行森林蓄积量估算效果相比于多元线性回归方法较差.尽管目前利用机器学习等方法构建的非线性模型应用较多，但是机器学习给人的感觉就像一个黑盒子，我们无法控制模型内部的运行，只能在不同的参数和随机种子之间进行尝试.同时对因变量的解释性较低，即无法直接解释测量变量与蓄积量估算模型的关系以及各变量对模型的贡献[37].

值得注意的是，在本研究中存在蓄积量高估的情况，主要原因可能是样地数据相对较少，且56 个样地实测蓄积量值多集中于100 m3·hm-2以上，在森林蓄积量建模时未考虑蓄积量较小值，影响了森林蓄积量精度.目前，广义神经网络和深度学习在森林蓄积量估算中的应用取得较好效果，在今后的研究中将尝试利用此类方法并与线性模型对比，以实现快速准确的森林蓄积量动态监测和定量评价.

4.2 结论

（1）通过Sentine-1 和Sentinel-2A 为数据源提取的80 个纹理变量、24 个光学特征变量和3 个后向散射系数3 种类型变量，CorB4、VarB5、VarB11、B11、PSSRa、REIP、等7 个变量对蓄积量估算模型贡献较大.

（2）5 组森林蓄积量估算模型中，结合不同类型变量构建的模型（d）效果最好，RMSECV为49.70 m3·hm-2，森林蓄积量估算精度达74%.

（3）西小山林场森林蓄积量整体呈现中南部较低，东部和中部偏西蓄积量较高的空间格局，预测平均蓄积量达到188.31 m3·hm-2.