依据光学遥感特征优选的森林地上生物量反演1)

王熙媛 张王菲 李云 杨仙保

(西南林业大学，昆明，650224)

森林对维持全球碳循环平衡和减缓CO2引起的温室效应起到了至关重要的作用，能调节区域气候和保护水源，因此掌握森林资源信息意义重大[1]。森林生物量指森林生态系统中植被有机体的干物质总量，是进行森林生态系统碳循环与碳储量研究的基础[2]。而森林地上生物量(AGB)作为森林生物量遥感估测的主要森林参数受到了广泛的关注[3-4]，其测定和估计成为森林研究及林业生产应用的热点问题[5]。

随着遥感技术的快速发展，合成孔径雷达(SAR)、激光雷达(LiDAR)和光学影像已经成为获取大范围森林AGB和碳储量空间分布的基础遥感影像。SAR虽然受云层的影响相对较小，且对树冠有一定的穿透能力[6]，但信号受地形起伏影响较大，在森林植被结构复杂、地形起伏较大的云贵高原山地地区有一定的局限性[7]。激光雷达为光斑激光传感器，无法达到无缝覆盖，很难在点云密度、数据分析精度和分析效率之间取得平衡，因此在大尺度应用上存在限制[8-9]。

基于此，光学遥感数据能获取到水平方向连续的区域性数据，其获取途径更为广泛，且部分数据免费对用户开放，易于获取。此外，高分辨率光学遥感影像包含了大量的纹理、形状等空间几何信息，在对不同森林类型以及地形起伏较大区域的森林AGB进行反演时，其提取的纹理和细节信息更加丰富[7，10]。与其他国外免费获取的光学遥感数据相比，例如Landsat 8 OLI、MODIS、哨兵2号(Sentinel-2A)等数据，我国高分系列产品在数据可得性方面同样更具优势，发展高分系列民用卫星对监测森林动态变化分析具有重大意义。

综上所述，研究以昆明市宜良县小哨林场为研究区，云南松林为研究对象，结合Landsat 8 OLI、高分1号(GF-1)、哨兵2号(Sentinel-2A)数据的植被指数、纹理特征和样地调查数据，根据韩宗涛等[11]提出的快速迭代特征选择的K最近邻法(KNN-FIFS)，从高维遥感数据中快速选取相关特征进行森林AGB优化建模，从而对森林AGB进行反演研究，对比分析光学数据反演森林AGB的潜力及探索GF-1在森林AGB反演上的潜力以及可替代性。

1 研究区概况

试验区位于云南省中部宜良县的花园林场(24°30′36″～25°17′2″N，102°58′22″～103°28′75″E，图1)。研究区整体地势北高南低，属于北亚热带季风气候，降水主要集中在夏秋季节，冬春季节干旱少雨。研究区内优势树种为云南松、华山松、高山栎，全县森林覆盖率达46.4%[12]。

图1 研究区地理位置及样地分布

2 研究方法

2.1 光学卫星数据及预处理

Landsat 8 OLI Landsat 8 OLI数据由2013年发射携带OLI陆地成像仪获得。研究获取的影像成像时间为2019年12月，云量小于10%，分辨率为30 m，波段为2～7波段(0.45～2.35 μm)。本文采用ENVI对获取的Landsat 8 OLI影像进行了辐射定标和FLAASH大气校正，预处理后的影像见图1(a)。

Sentinel-2A Sentinel-2A是高分辨率多光谱成像卫星，有13个光谱波段，空间分辨率分别为10、20和60 m。在光学数据中，Sentinel-2A数据是唯一在红边范围含有3个波段的数据，对监测植被变化情况相对敏感，同时提供了更多波段选择[13-14]。本文获取的Sentinel-2A影像成像时间为2019年12月，为L2A级影像，使用的波段为2～8波段(0.49～0.84 μm，其中5～7为红边范围波段)和11～12波段(1.61～2.19 μm)共9个波段。本文采用SNAP软件对获取的Sentinel-2A影像进行了预处理，并将各波段数据重采样至10 m分辨率，镶嵌裁剪得到覆盖研究区影像图，见图1(b)。

GF-1 GF-1卫星是我国在高分辨率对地观测卫星系统重大专项中的第一颗卫星，其空间分辨率最高可达2 m。研究获取了覆盖研究区的一景GF-1影像，成像时间为2019年4月，空间分辨率为8 m，使用的4个波段波长为0.45～0.89 μm。本文采用ENVI插件对获取的GF-1数据进行了辐射定标、大气校正及正射校正，预处理后的影像见图1(c)。

2.2 样地数据获取及预处理

研究区样地地面调查数据分别于2019年8月和2020年8月分两次开展。2019年的样地数据采用手持GPS定位，定位误差在5 m以内，共获得89个云南松纯林角规样地点；2020年采用载波相位差分技术(RTK)定位，定位误差在1 m以内，共获得云南松纯林37个样地点，样地点布见图1。

样地生物量的计算采取蓄积—生物量的转换方法进行[12]。首先利用角规控制检尺法获得单位面积林木蓄积量，再利用云南松林蓄积量—生物量的转换公式完成云南松生物量的计算。研究中采用黄从德等[15]提出的蓄积量—生物量转换公式：

W=0.856 9×V0.856 4。

(1)

式(1)中V为云南松林分每公顷蓄积量；W为相对应的地上生物量。

2.3 特征提取

植被指数是森林参数估测中最常用的特征变量，它是由多光谱数据经过线性或非线性变换构成的能反映植被生长状态和分布的指数，构造基础是绿色植物对可见光和近红外波段截然不同的吸收和反射特征[16]。本文中提取的植被指数包括：归一化植被指数(NDVI)、简单植被指数(SR)、可见光大气阻抗植被指数(ARVI)、差值植被指数(DVI)、垂直植被指数(PVI)、比值植被指数(RVI)、B4和B5归一化差值植被指数(NDVI45)、绿波归一化差值植被指数(GNDVI)、反红边叶绿素指数(IRECI)、“哨兵2号”红边位置指数(S2REP)[17]共10种植被指数，其中后4种为哨兵2号数据特有，相关计算公式见文献[17]。

纹理是遥感图像的重要特征之一，它不仅可以反映粗糙程度，还可揭示影像中地物的结构信息及它们与周围环境的关系，是反映地表覆盖类型空间变化的重要信息[18]。本研究基于灰度共生矩阵(GLCM)提取了8种纹理信息，具体包括均值(Mea)、方差(Var)、均匀性(Hom)、对比度(Con)、相异性(Dis)、熵(Ent)、二阶矩(Sec)和相关性(Cor)[19]。

2.4 特征重要性排序

在提取的特征数据集中，众多特征全部参与建模会造成信息冗余，从而导致反演精度降低，为了得到高精度森林生物量反演模型，需要筛选出和森林AGB相关性较强的特征参与建模。通过R语言实现随机森林(RF)算法进行特征的重要性排序[19-20]，再选择重要性排名在前20的特征进行后续建模。

2.5 KNN-FIFS算法建模

K-NN是一种监督学习算法，其既可以解决分类又可以用于回归，即对于给定的测试样本，根据距离度量找出与其最靠近的K个训练样本，然后基于K个相邻训练样本的信息来预测未知量参数的值。本文中K个样本的属性值为待选样地的森林地上生物量，带估测像元的生物量由马氏距离下最靠近的K个训练样地的森林地上生物量通过距离加权获得。KNN算法不依赖于特定的函数分布，采用留一法(LOO)进行训练和验证，因此对训练和检验样本依赖性低，当输入特征向量较少而训练样本较多，其预测结果较稳健。然而，在特征维数较高时(比如在森林AGB遥感定量反演的实际应用中)，会造成特征组合优选效率低下，降低模型预测效率[21-22]。

为了克服KNN森林AGB估测中，特征组合效率低、模型预测精度低的缺陷，本文采用韩宗涛等[11]提出的KNN-FIFS算法估测森林AGB，该算法目前在多个森林AGB、森林郁闭度遥感定量反演中精度相对最高[19，21，23]。采用该算法进行森林AGB反演时，首先从样地数据和重要性排名靠前的遥感特征数据中提取一组训练数据，并将其带入KNN-FIFS算法，得到初始化均方根误差(RMSE-0)，将RMSE-0设置为理论上的极大值，然后利用训练数据建立森林AGB估测模型，从中优选出RMSE-b，如果RMSE-b

2.6 精度评价方法

KNN-FIFS方法估测AGB的精度由R2和均方根误差(RMSE)来衡量，R衡量模型点的准确性，常用R2表示，其越接近于1，表示模型精度越高；RMSE采用留一交叉验证计算，RMSE值越低，回归模型更加准确。R2和RMSE的计算公式如下[19]：

(2)

(3)

3 结果与分析

本文首先使用RF算法对3种光学数据提取出的参数特征信息进行重要性排序，并选择重要性排序结果中前20个遥感参数特征作为KNN-FIFS算法的输入特征，反演研究区森林AGB，并分别对3种数据估测结果及3种数据组合估测结果进行精度评价。研究结果中，阐述了Landsat 8数据、Sentinel-2A数据、GF-1数据以及结合3种数据分别采用KNN-FIFS算法最终选中的最优特征参数组合；图2(a)为基于Landsat 8森林AGB反演制图结果；图3(a)为基于Landsat 8森林AGB反演结果与地面调查结果的相关性分析。图2(b)和图3(b)为基于Sentinel-2A数据的结果；图2(c)和图3(c)为基于GF-1数据的结果；图2(d)和图3(d)为3种光学数据组合进行森林AGB反演的结果。

3.1 Landsat 8数据地上生物量反演结果

根据李春梅[16]的研究结果，在KNN算法中，当参数为20个时，反演结果的RMSE值最小，因此本文选用RF排序结果中的前20个参数作为KNN-FIFS算法的初始输入特征，并用“%IncMSE”和“IncNodePurity”做为重要性评价指标，其值越大表示该变量的重要性越大[24]。

Landsat 8数据提取的参数中，重要性排序前20的特征为：WL_G3_mea、WLS1_7_mea(WL为纹理，G3_mea为绿光波段窗口3×3的均值纹理信息；S1_7_mea为短波红外1波段窗口7×7的均值纹理信息，其余以此类推)、WLS2_3_mea、WLS1_9_mea、WLS2_9_mea、WLS2_7_mea、WLS1_3_mea、WLS2_9_var、WLS1_3_var、WL_G9_mea、WLS1_9_sec、WLR7_sec、WLR9_var、WLS2_3_var、WL_G9_var、WL_R9_mea、WLR7_mea、WLR9_hom、WLR3_mea、WL_G9_ent，这些选中的特征的“%IncMSE”和“IncNodePurity”值的范围分别在6.84～14.36和1 208.57～7 442.80之间。由此可以看出，Landsat 8数据提取的参数中，红色和红外波段的均值纹理占比较大，说明重要性高。

本文K值和窗口大小值的设置范围在1～11之间。采用KNN-FIFS取得研究区最优的森林AGB估测结果时，Landsat 8数据基于马氏距离的KNN-FIFS生物量反演模型获得的最优特征组合为归一化植被指数，VARSWIR2_3*3，DISSWIR2_7*7，CONgreen_7*7,VARR_7*7，MEAR_7*7(VARSWIR2_3*3为短波红外波段SWIR2的3×3窗口的方差纹理信息；DISSWIR2_7*7为短波红外波段SWIR2的7×7窗口的相异性纹理信息，以此类推)，当K=2，窗口大小为11×11时，估测结果与地面验证数据之间的R2为0.47，RMSE为23.29 t/hm2。采用Landsat 8数据的制图结果见图2(a)，可以看出生物量分布不足30 t/hm2的森林区域约占10%，30～90 t/hm2的森林区域约占80%，其余区域约占10%。

3.2 Sentinel-2A数据地上生物量反演结果

Sentinel-2A数据提取的参数中，重要性排序前20的特征为：s_SW2mea9、s_R9mea(s为Sentinel-2A简写，SW2mea9为短波红外2窗口9×9的均值纹理信息；s_R9mea为红波段窗口9×9的值纹理信息，其余以此类推)、s_R9con、s_SW1ent7、s_SW2ent9、s_SW2sec9、s_SW1mea3、s_R9dis、s_SW1ent9、s_R9sec、sen_R7mea、s_SW2var9、s_R9hom、s_SW2hom9、s_R9ent、s_R9var、s_SW2con9、s_SW1mea7、s_SW1sec3、s_R7var，其“%IncMSE”和“IncNodePurity”的值范围在5.55～18.32和802.79～5 471.42，与Landsat 8数据类似，Sentinel-2A数据提取的参数中，短波红外和红色波段的纹理信息占比较大，重要性较高。

Sentinel-2A数据基于马氏距离的KNN-FIFS生物量反演模型获得的最优特征组合为SECSWIR1_9*9，MEASWIR1_3*3，VARR_9*9，DISR_9*9，当K为4，窗口大小为5×5时，研究区森林AGB估测结果最优，其估测精度为R2=0.60，RMSE=21.40 t/hm2，采用哨兵2A数据的制图结果中(见图2(b))，生物量分布不足30 t/hm2的森林区域约占42%，30～90 t/hm2的森林区域约占55%，其余区域约占3%。

3.3 GF-1数据地上生物量反演结果

GF-1数据提取的参数中，重要性排序前20的特征为：gfB9_mea、gf_SR(gfB9_mea表示GF-1数据蓝波段9×9窗口的均值纹理信息，gf_SR即表示GF-1简单植被指数，其余以此类推)、gfB7_mea、gf_NDVI、gfG9_mea、gfR9_mea、gf_RVI、gfN7_dis、gfR3_mea、gfR7_mea、gfG3_mea、gfG3_cor、gfR9_mea、gfG3_sec、gfR9_mea、gfG7_mea、gfG3_ent、gfDVI、gfG7_sec、gfN9_sec、gfR7_ent、gfN9_ent，其“%IncMSE”和“IncNodePurity”的值范围在7.45～15.34和1 295.92～2 885.30，可以看出植被指数NDVI、RVI、SR以及绿波段的均值纹理占比较大，重要性较高。

图2 宜良县森林AGB分布图

基于GF-1数据、马氏距离的KNN-FIFS生物量反演模型获得的最优特征组合为MEAB_9*9，DISN_7*7，当K=9，窗口大小为1×1时，研究区森林AGB估测结果最优，其估测精度为R2=0.59，RMSE=22.11 t/hm2，采用GF-1数据的制图结果中(见图2(c))，生物量分布不足30 t/hm2的森林区域约占10%，30～90 t/hm2的森林区域约占80%，其余区域约占10%。

3.4 3种数据组合进行森林地上生物量反演结果

3种数据组合提取的参数中，重要性排序前20的特征为：LS_G9_mea、LS_S1_7_mea(LS_G9_mea表示Landsat 8绿波段9×9窗口的均值纹理信息；LS_S1_7_mea表示Landsat 8短波红外1波段7×7窗口的均值纹理信息，以此类推)、LS_S2_9_var、LS_S1_9_mea、LS_S2_9_mea、LS_S2_7_mea、LS_G3_mea、LS_S2_3_mea、LS_S1_3_mea、gf_B9_mea3、s_R9con、gfR7_mea3、gf_RVI、gf_SR、gf_NDVI、s_R9mea、gfR9_mea3、gfG9_mea3、gfB1_7mea3、gfR3_mea3，其“%IncMSE”和“IncNodePurity”的值范围在4.47～20.65和854.28～10 095.17，可以看出Landsat 8纹理特征占比较大，GF-1次之。

基于3种光学数据组合、马氏距离的KNN-FIFS生物量反演模型获得的最优特征组合为MEAB_7*7_GF，CON_B_9*9_SEN，MEA_R_9*9_LS，NDVI_GF(MEAB_7*7_GF表示GF-1蓝波段7×7窗口的均值纹理信息；MEA_R_9*9_LS表示Landsat 8红波段9×9窗口的均值纹理信息，以此类推)，当K=9，窗口大小为9×9时，研究区森林AGB估测结果最优，其估测精度为R2=0.42，RMSE=23.86 t/hm2，制图结果见图2(d)，图中生物量分布不足30 t/hm2的森林区域约占10%，30～90 t/hm2的森林区域约占85%，其余区域约占5%。

综合分析基于Landsat 8、Sentinel-2A、GF-13种光学数据森林AGB的反演结果，可知Sentinel-2A和GF-1数据的反演结果精度均较好，RMSE值较低，而将3种数据组合进行森林AGB的反演结果反而低于仅采用单一数据源的反演结果。由图2和图3可知，Landsat 8数据在森林AGB水平较低时(70 t/hm2)反演结果精度较高，高于70 t/hm2时出现了饱和现象；基于Sentinel-2A数据则在研究区森林AGB变化范围内未见明显的饱和现象，其在整个森林AGB变化范围，均有较好的估测精度；基于GF-1的数据在低生物量和高生物量水平其估测结果精度均较高，但是当森林AGB在70～100 t/hm2变化时，其反演结果精度较低。

图3 KNN_FIFS估测森林AGB交叉验证结果

4 结论与讨论

在林业工作上，森林生物量这一重要评价指标一直存在基础资料不全、传统地面实测人力物力消耗较大的问题，是在进行区域尺度上森林生物量遥感估测面临的最大难题。从1978年开始，我国在全国范围内建立了以省、自治区、直辖市为总体的森林资源连续清查体系，各县慢慢积累了全面的固定样地调查数据[24]。2008年后，随着航空航天技术的发展，包括Landsat 8、国产高分系列、Sentinel-2A/B等在内的各种免费光学数据源日趋丰富。本研究中，基于几种常用的光学遥感影像研究了采用KNN-FIFS方法进行森林AGB反演的结果。

研究结果中，Landsat 8 OLIs数据估测结果与实测值之间的R2为0.47，RMSE为23.29 t/hm2，研究与朱妍[25]的研究结果类似，该研究同样采用Landsat 8 OLIs，但区分了森林类型，研究将森林类型划分为针叶林、阔叶林和混交林，各类森林类型的森林AGB反演结果与地面测量结果的R2在0.44～0.48；此外，曹林等[26]基于机载小光斑LiDAR数据提取亚热带森林参数信息，运用逐步回归方法拟合的林分蓄积模型R2为0.46～0.55，略高于本文研究结果。许振宇等[24]基于Landsat 8数据的研究中，采用随机森林模型估测区域森林AGB，R2达到了0.65，将该参数与Sentinel-1数据联合反演，结果显示R2提高到0.72，其R2明显提高可能是由于SAR信息的输入以及非参数估测方法的使用。李云等[19]利用光学和SAR数据进行森林AGB的反演也证实了引入SAR数据可以有效提高森林AGB的反演结果的精度。

Sentinel-2A数据估测结果与实测值之间的R2为0.6，RMSE为21.40 t/hm2，该结果略高于曹霖等[27]采用最优估测模型(随机森林)和Sentinel-2A数据估测的结果，该研究结果的RMSE为35.36 t/hm2。在陈松等[28]的研究中，选用Sentinel-2A数据和机载LiDAR作为研究数据，采用四种估测模型做对比，分别是普通回归模型、误差变量联合方程模型、随机森林模型和KNN估测模型，得出最优模型为误差变量联合方程组R2为0.60，RMSE为48.64 t/hm2；单一Sentinel-2A数据源的研究中，郭正齐等[17]的森林AGB反演精度最高，其R2为0.77，RMSE为39.49 t/hm2，其特别加入了生物物理参数作为建模因子之一，均采用多元回归模型进行估测反演，得到精度较高的反演结果；在Sentinel-1和Sentinel-2A联合反演的研究中，单一数据反演精度R2为0.50，联合SAR数据可将R2提高到0.57，这些研究的研究结果表明Sentinel-1的SAR影像纹理特征值含有一定的辅助信息，对于模型精度的提高具有一定的作用[29]。

GF-1数据估测结果和实测值之间的R2为0.59，RMSE为22.11 t/hm2，在曾晶等[30]基于GF-1数据，采用多元回归模型估测森林AGB的结果相比，该研究结果的R2为0.65；张少伟等[31]结合SAR数据协同反演，采用KNN-FIFS模型估测结果与实测值之间的R2为0.56，RMSE为25.95 t/hm2，该估测结果与本文研究结果接近；周俊宏等[32]研究表明GF-1数据的纹理及波段信息能较好表达森林生物量，在估算生物量中具有明显优势，其R2为0.49～0.77，优于线性模型(R2为0.39)；构建的随机森林经验模型在遥感估算森林地上生物量中具有重要潜力(R2为0.77)。

在本研究中，以Landsat 8、Sentinel-2A、GF-13种光学数据，采用KNN-FIFS进行森林AGB反演，结果表明：Sentinel-2A和GF-1数据的反演结果精度均较好，RMSE值较低，Landsat 8反演森林AGB结果最差。在参数提取方法相同的情况下，数据源的空间分辨率大小对森林AGB估测结果影响明显，Sentinel-2A数据和GF-1数据空间分辨率相近，反演精度差异不大，而Landsat 8数据空间分辨率为30 m，图像存在大量混合像元，在森林AGB较低的区域，Landsat 8像元值很容易受到其他地物影响，不能正确反映实际的生物量信息；而在森林AGB较高的区域，由于生物量反演普遍存在饱和问题，导致像元实际生物量被低估。从特征组合上来看，3种数据单独的特征组合中，均值和非相似性都被筛选参与模型算法中，从原理上分析因为均值反映了像元纹理规律性，值越大表示纹理规律明显，易于描述，说明均值与像元内反映森林AGB规律相匹配；非相似性表示像元内灰度差别越大，图像的视觉效果越清晰，说明像元内植被覆盖部分与非植被部分分异性大，能较好的提高模型反演精度。3种数据结合的特征组合中，3种数据的均值、对比度指标都被筛选加入模型反演，特别是GF-1的NDVI指数被筛选，而其他两种数据没有NDVI加入，在一定程度上更加突出我国高分数据在定量反演上的优势。从特征重要性排序上来看，7×7和9×9的均值纹理特征在3种数据中重要性排序中相对靠前，说明其对森林AGB的解释能力强及在森林AGB反演中具有一定优势。研究表明：高分辨率光学数据有助于提高森林AGB反演精度，国产高分系列数据可以替代现有免费国外数据源，其反演森林AGB结果精度优于或等同于目前常用的两种免费国外数据。由于文中仅采用3种光学数据且没有SAR数据进行比较，而且在不同地区、不同温度带、不同森林类型中，其反演精度差异有待进一步研究。