赵翊含，侯蒙京，冯琦胜，高宏元，梁天刚*，贺金生，2，钱大文

（1. 兰州大学草地农业科技学院，草地农业生态系统国家重点实验室，兰州大学农业农村部牧草创新重点实验室，兰州大学草地农业教育工程研究中心，甘肃 兰州 730020；2. 北京大学城市与环境学院，北京 100871；3. 中国科学院西北高原生物研究所，青海 西宁 810008）

草地生物量是反映草地生长状况和生态环境评估的一项重要指标［1］，准确地估测草地产量一直是畜牧业科研工作者持续探索的目标［2］。在高寒地区，由于过度放牧等因素引起的草地退化、水土流失等环境恶化问题日益严重；并且生态系统受干扰后植物生长和恢复速度缓慢［3］。为了确保高寒地区畜牧业的可持续发展，对于准确估算高寒草地地上生物量十分必要［4］。

传统地面草地调查的方法耗费成本太高［5］，刈割法对草地的破坏性大，不利于草地的再生与可持续。基于遥感技术对草地资源监测省时、省力，能快速且客观地对草地生长现状做出评估，适合大尺度估测草地地上生物量（above-ground biomass，AGB）。自20 世纪80 年代，有很多国内外学者利用多源遥感数据开展AGB 估算研究。如姚兴成等［6］利用中分辨率成像光谱仪（moderate-resolution imaging spectroradiometer，MODIS）和地面实测资料，以云南省为例建立草地AGB 遥感估算模型，用草地植被群落特征优化了模型精度；于惠等［7］构建古浪县生物量估算模型，结果表明基于MODIS 短波红外波段的归一化耕作指数（normalized difference tillage index，NDTI）与草地AGB 的相关性最好；冯琦胜等［8］基于青藏高原地面实测样本数据与MODIS-NDVI 构建生物量风干重估算模型，研究表明MODIS-NDVI 的指数函数模型是最优模型。尽管利用MODIS 遥感数据可以估测草地生物量，但是中分辨率影像存在空间分辨率较粗（大于250 m）和混合像元等问题，致使草地AGB 模型的精度不高，不适合县域尺度高精度精细化的遥感监测。

Landsat-8 OLI（operational land imager，OLI）于2013 年发射，空间分辨率达30 m，与MODIS 相比其空间分辨率有实质性的提升，并且包含了9 个波段，其中包括一个15 m 的全色波段，适合对草地AGB 进行更精细的遥感估算。Landsat-8 OLI 改进了波段数量、光谱范围及辐射分辨率等方面，在性能上相比之前几代Landsat 卫星拥有明显改善，对草地AGB 估测的准确性有提升作用［9］。李斌等［10］利用Landsat 8 遥感影像对不同处理阶段［数字量化值（digital number，DN）、辐射定标、大气纠正］的归一化植被指数（normalized difference vegetation index，NDVI）进行对比研究，表明应该采用经过大气校正后的值作为植被覆盖的定量指标；Li 等［11］利用Landsat-8 增强型植被指数（enhanced vegetation index，EVI）建立了与植物干物质含量的回归模型，评价分析了整个青藏高原草甸和草原群落地上生物量的变化特点，并比较了植物性状。乌如汗［12］以内蒙古正蓝旗为试验区，利用Landsat-8 OLI 遥感影像和地面实测数据，确定了NDVI 植被指数的一元二次回归模型作为草地地上生物量反演模型。以上研究表明Landsat-8 可以实现对草地AGB 的有效估测。然而，现有研究多采用传统参数化回归方法构建模型，反演误差相对较大，因此在提高草地AGB 监测精度上仍然有必要进一步开展深入研究。近年来，机器学习方法在植被遥感监测方面取得了良好的效果［13］。它通常基于完整的光谱集进行建模，能够充分利用光谱信息，具有非线性、预测准确率高、泛化能力强的特点［14］，因此在遥感领域基于非参数化的机器学习算法估测AGB 逐渐受到重视［13］。

基于以上因素考虑，本研究利用Landsat-8 OLI 遥感数据结合青海省门源县的地面实测数据，将计算出的植被指数作为自变量，构建传统回归模型和随机森林模型，对比分析模型精度，确定遥感反演生物量的最佳模型，以期为当地草地资源可持续利用和草畜科学管理提供理论依据。

1 材料与方法

1.1 研究区概况

研究区位于青海省海北藏族自治州门源回族自治县，地处青藏高原腹地，地理位置为北纬37.5°-38.0°，东经102.0°-102.5°（图1）。门源回族自治县地势西北高，东南低，平均海拔2800 m 以上，属于高原大陆性气候，夏季凉爽短暂，冬季寒冷漫长，年平均气温-1.7 ℃，年平均降水500 mm。该地区天然草地类型主要为高寒草甸类，主要以莎草科牧草嵩草（Kobresia myosuroides）、黑褐苔草（Carex atrofusca）和禾本科牧草为优势种［15］。门源县是河湟地区和河西走廊重要的水源涵养区和补给地，基于农牧资源丰富的区位优势，积极贯彻“立草为业、畜牧强县”的发展方针，大力发展高寒地区现代高效畜牧业及很多特色种养业。

1.2 数据来源

1.2.1实测数据获取 本研究共获得202 个地面实测调查样本数据，包括青海海北高寒草地生态系统国家野外科学观测研究站（简称海北站）多年连续的观测数据和全县范围内开展的草地AGB 外业调查数据（图1）。观测站数据获取时间为2019 年5-8 月、2020 年7 月和2021 年7-8 月，共得到178 个样本数据，全县范围内AGB 外业调查共获得24 个样点数据，调查时间为2021 年7 月1 日-7 月8 日。在对草地生长期内的AGB 进行采集时，每个样点选取3 个随机样方，样方大小为50 cm×50 cm，齐地剪取地上生物量，并记录各样方测量时间、经纬度、海拔、草地盖度、草层高度、鲜重，将样本装袋编号带回实验室进行65 ℃烘干48 h 后称重。最终，取3 个样方的平均值作为样本生物量，总计有202 个观测样本。

图1 研究区海拔及采样点空间分布Fig.1 Spatial distribution of altitude and sampling points in the study area

1.2.2遥感数据获取 遥感影像数据选用Landsat-8 OLI 数据，是美国陆地卫星计划的第8 颗卫星，传感器类型为陆地成像仪（OLI），时间分辨率为16 d，空间分辨率为30 m。下载的陆地产品级别为二级产品（Collection 2 Level-2），该产品已经经过辐射定标、大气校正和几何精校正。为了使影像数据的成像时间与地面实测数据的调查时间差异较小，下载影像时选择影像完全覆盖样地，成像时间与采样时间最接近的影像。选取2019-2021 年的7-8 月覆盖门源县采样点及全境、云量覆盖<30%的Landsat-8 OLI 影像，总计7 景（表1）。研究区实测草地生物量样本数据有202 个，剔除受天气状况影响的样本，实际参与模型构建的样本数为199 个。

表1 所用Landsat-8 OLI 影像信息Table 1 Landsat-8 OLI image information used

1.3 模型构建及精度验证

1.3.1单变量回归模型建立 利用ArcGIS 10.2 软件结合样地经纬度坐标，导出Landsat-8 OLI 波段2~7，并利用Extraction 工具提取与199 个地面采样点对应的像元值。根据公式计算29 种植被指数（表2）。分别构建29种植被指数与AGB 的单因素回归模型，将植被指数作为自变量，AGB 作为因变量，二者之间建立线性、对数、指数、乘幂4 种拟合模型；采用10-fold 交叉验证的方式对模型进行精度评价。

表2 植被指数计算公式Table 2 Calculation formula of vegetation indexes

1.3.2变量筛选 为了筛选出合适的建模因子，本研究在RStudio 中利用“glmnet”程序包实现基于最小绝对收缩筛选方法（the least absolute shrinkage and selection operator，LASSO）的变量选择。这种方法可以将变量的数量显著减少，从而实现降维，被广泛应用于高维数据中选择敏感变量［16］。LASSO 变量筛选方法以模型系数的绝对值函数作为惩罚项施加在普通最小二乘法的损失函数上，以此来压缩模型系数，使一些解释因变量效果不好的系数变小甚至压缩至0 从而移出模型，因此LASSO 方法可以提供一种稀疏解，能够同时进行变量筛选和参数估计［17-18］。

1.3.3随机森林模型构建 随机森林算法（random forest，RF）是一种新型机器学习算法，是运用多棵决策树训练样本并集成预测的一种非参数机器学习算法，是从原始样本中，利用bootstrap 重抽样技术随机抽取数据后构造多个样本，然后对每个重抽样样本采用节点的随机分裂技术构造N 棵决策树［19］。在每棵树生长过程中，从全部特征变量中随机抽选mtry 个进行内部节点划分；最后，将N 棵决策树的预测结果集合起来，采用投票的方式决定新样本的类别［20］。ntree 参数值上限一般设置为1000，大量研究已证明该值对许多RF 程序有效［21］。利用“random forest”程序包实现随机森林算法，需要定义2 个参数：决策分类树的数目（ntree）和节点分割的特征数目（mtry），参数影响着模型的精度，本研究采用的参数决策树数量ntree 是500、特征个数mtry 是4。

1.3.4精度验证 为了比较分析传统回归模型和RF 模型对AGB 的反演效果，这两种模型均采用十折交叉验证的方法［16］，将所用的样本数据随机分成10 份，其中9 份作为训练集，1 份作为验证集，分别计算均方根误差（root mean square error，RMSE）、决定系数（coefficient of determination，R2）评价AGB 估测模型的精度［22］。计算公式如下：

式中：Yi表示样本实测值，Yi表示模型反演值，yi表示样本生物量实测平均值，n表示样本数量；RMSE 越小，表示拟合精度越好；R2越接近1，拟合程度越好，参考价值越高［1］。

2 结果与分析

2.1 生物量样本数据分析

2019-2021 年门源县用于建模的199 个实测样本的草地生物量统计结果见表3。从该表可以看出，海北观测站 生物量最大值为4142 kg·hm-2，最小值为187.2 kg·hm-2，平均值 为2061 kg·hm-2，变异系 数（coefficient of variation，CV）为0.47。外业调查数据覆盖门源县主要草地类型，生物量最大值为4666.67 kg·hm-2，最小值为637.07 kg·hm-2，平均值为1824.01 kg·hm-2，变异系数（CV）为0.56。

表3 观测样本生物量汇总Table 3 Biomass statistics of observed samples

2.2 单因素植被指数回归模型

续图2 单因素植被指数回归模型预测值与实测值关系Continued Fig.2 Relationship between predicted value and measured value of single-factor vegetation indexes regression model（n=199）

从结果来看，NDVI 与AGB 的对数模型精度达到最优，验证集R2达0.50，RMSE 为702.89 kg·hm-2。单因素植被指数构建的模型的验证集R2介于0.37~0.50，RMSE 介于702.89~792.95 kg·hm-2。所有植被指数最优模型 统 计 分 析 的 结 果 表 明，NDVI、EVI、RVI、GNDVI、DVI、RDVI、SAVI、MSAVI、OSAVI、PVI、TVI、TVI2、MSR、GSAVI、GOSAVI、TSAVI、RGBVI、VDVI、GRVI、GCI、RBNDVI、BNDVI 最优拟合模型为对数，其他植被指数的最优拟合模型为线性。图2 为所有植被指数最优模型预测值和实测值的关系。验证集R2较高的植被指数前10 个分别是NDVI、RBNDVI、TVI、GNDVI、MSR、TSAVI、GCI、GRVI、RVI、panNDVI。

图2 单因素植被指数回归模型预测值与实测值关系Fig.2 Relationship between predicted value and measured value of single-factor vegetation indexes regression model（n=199）

2.3 随机森林模型分析

本研究对比了RF 模型变量筛选前后的模型精度。以29 种植被指数建立的模型，其验证集R2为0.61，RMSE为621.14 kg·hm-2（图3a）。利用LASSO 变量筛选方法，采用十折交叉验证的方法确定模型均方误差最小时λ0=0.25，并取距离λ0一个标准差位置的λ1（1.22）作为最终的惩罚系数（图4）。从29 种植被指数中选出11 种植被指数，其重要性排序分别是MNLI、MSAVI、GRNDVI、MVI、SAVI、GOSAVI、GRVI、RVI、RGBVI、GSAVI、GCI（图5）。通过变量筛选所构建的AGB 估测模型为最优，验证集R2为0.62，RMSE 为621.95 kg·hm-2（图3b）。由此可见，在保证精度的基础上LASSO 变量筛选对草地AGB 反演模型可以有效进行降维和简化。以上2 种RF 模型的精度均高于单因素回归模型，其中最优RF 模型的R2比最优单因素模型提高0.12，RMSE 降低了80.95 kg·hm-2。

图3 随机森林模型预测值与实测值关系Fig.3 Relationship between predicted value and measured value of RF model

图4 LASSO 筛选结果Fig.4 LASSO results

图5 筛选出的植被指数的重要性Fig.5 The importance of screened vegetation indexes

2.4 门源县AGB 空间分布

经过变量筛选的随机森林草地AGB 估算模型最优，利用该模型对研究区2019-2021 年的AGB 进行了反演（图6）。可以看出，草地AGB 高值主要集中在西北部，东南部相对较低，总体呈中部高，四周低的趋势。2019 年草 地AGB 值 主 要 在0~800 kg·hm-2和1500~2000 kg·hm-2；2020 年 主 要 在0~800 kg·hm-2和1000~1500 kg·hm-2；2021 年主要分布在1000~1500 kg·hm-2和1500~2000 kg·hm-2。2019-2021 年全县天然草地总产草量介于4.2827 万~8.9776 万t，平均单产介于1063.49~1484.82 kg·hm-2。门源县草地类型以山地草甸、温性草原和高寒草甸三类为主，其中分布最广泛的是高寒草甸，2019-2021 年产草量介于4.0825 万~5.6653 万t，平均AGB为1060.38~1471.94 kg·hm-2。 山 地 草 甸2019-2021 年 产 草 量 介 于973.7687~1571.6790 t，平 均AGB 介 于1036.81~1637.43 kg·hm-2；温 性 草 原2019-2021 年产草量介于746.6281~1112.7140 t，平均AGB 介于1198.72~1786.63 kg·hm-2（图7）。

图6 门源县2019（a）、2020（b）和2021（c）年天然AGB 空间分布Fig.6 Spatial distribution of AGB of natural grassland in 2019（a），2020（b）and 2021（c）in Menyuan County

图7 门源县3 种主要草地类型的平均地上生物量统计Fig. 7 Average aboveground biomass of three main grassland types in Menyuan County

3 讨论

现阶段国内外研究基于植被指数构建生物量估算模型的方法较为常见，其数值可以代表植被活力，比单波段更具有灵敏性［23］。NDVI 是目前在遥感领域应用最广泛的一种植被指数，可以捕获植被整个生长季的动态，可用于植物生物量的预测，在各类植被指数中，NDVI 可以较好反映AGB，对植被生长的相关信息敏感，如水分含量和植被覆盖等［24］。从本研究对29 种植被指数单因素建模的研究结果来看，NDVI 用于草地AGB 的估测模型在传统统计回归模型中精度最优，这与许多学者［25-27］的研究结果一致。也有很多学者研究认为SAVI 比NDVI 更适合估算草地地上生物量，因为SAVI 适应植被密度变化及消除土壤影响的能力较强［28-29］。但是，本研究区分布范围最大的草地类型是高寒草甸，其植被覆盖度较高且均一，所以土壤背景对植被估测的影响不如其他学者的研究区域明显［4］，并且这是单因素估算模型，其反演结果不如NDVI 与AGB 估测效果，这与杨鹏万等［4］的研究结果类似。

除NDVI 以外，单因 素 回 归 模 型R2较高的前10 个植被指数依次分 别 是NDVI、RBNDVI、TVI、GNDVI、MSR、TSAVI、GCI、GRVI、RVI、panNDVI。经过LASSO 变量筛选后构建研究区草地AGB 反演模型的植被指数有RVI、SAVI、MSAVI、MNLI、GSAVI、GOSAVI、RGBVI、GRVI、GCI、MVI、GRNDVI。由此可见，传统回归结果与RF 建模筛选出的变量有一定差异。值得一提的是，在29 种单因素植被指数回归模型中，NDVI 与AGB 的相关性最高，但是在随机森林模型构建筛选变量的过程中NDVI 没有被筛选出来，筛选出来的相近植被指数是GRNDVI。这是因为GRNDVI 是基于NDVI 基础上改进的植被指数，GRNDVI 与NDVI 具有较好的线性相关关系，利用这两种植被指数估算叶绿素可以达到相似的精度［30］。GRNDVI、RVI、GRVI、GCI 在变量筛选中被筛选出来，同时也是单因素回归模型中与生物量拟合精度较好的几种植被指数。但在RF 模型中其他被筛选出的植被指数（SAVI、MSAVI、MNLI、GSAVI、GOSAVI、RGBVI、MVI），其单因素模型的精度并不高。多种植被指数性能具有交互影响和互补作用可能是导致这种现象的重要原因。两种模型采用的是完全不同的变量筛选方法，单因素模型是植被指数各自分别与AGB 的相关性，不能涵盖多种因素的综合影响，基于单因素构建的回归模型的误差偏大；而在RF 模型中是所有筛选出来的植被指数共同影响并决定了模型的精度，明显优于传统回归模型。

比较分析最优RF 模型涉及的多种植被指数，可以看出：1）基于近红外和红光波段组合计算的植被指数，是反演草地AGB 的重要变量。本研究由LASSO 筛选出的11 个植被指数作为建模因子构建了最优RF 模型，其中，有4 个植被指数是基于Landsat-8 近红外和红光波段计算的，分别为RVI、SAVI、MSAVI、MNLI。已有研究表明，利用红边和近红外光谱计算的植被指数可以最大限度地减少大气和水分吸收等因素的影响［31］，研究区具有高寒阴湿的气候特点，这可能是这类植被指数对高寒地区草地AGB 有显著影响的关键因素之一。2）由绿波段参与构建的植被指数对高寒草地AGB 也具有重要影响。大部分植被指数多基于近红外和红光波段反射率组合计算，然而近年来研究发现以近红外光和绿光波段构建的植被指数有时更加敏感［32-33］。本研究筛选出的GSAVI、GOSAVI、RGBVI、GRVI、GCI、GRNDVI 这几种植被指数的计算均有绿波段，这表明基于绿波段构建的植被指数对高寒地区的草地生物量的精准估算也具有显著影响。3）MVI 和MNLI 这两个植被指数对干旱区的植被更为敏感，这可能与本研究区地处青藏高原腹地，属半干旱地区密切相关。首先短波红外波段参与构建的植被指数对特定区域的草地植被的AGB 更为敏感。MVI 植被指数是基于红光、近红外和短波红外波段计算的植被指数，于惠等［7］研究表明短波红外波段对干旱区草地植被更为敏感；此外，MNLI 是非线性植被指数的改进指数，它与叶绿素含量具有较高的相关性，在半干旱环境下，MNLI 可以将植被指数与地表生物物理参数的非线性化关系线性化［34］。故MVI 和MNLI 被筛选出来，对处于半干旱区的植被生物量估算也具有一定的作用。综上所述，本研究通过变量筛选后，用选出的多种类型的植被指数建模，这不仅简化了模型，也从叶绿素含量、消除水分影响、植被密度等方面综合反映了植被特征，多种植被指数与草地AGB 之间具有显著的相关性，并且不同植被指数之间也具有一定的互补性，从而可以综合地反映草地生物量的状况。

相较于传统回归分析方法，机器学习算法更适用于较复杂的运算，可以更好地进行变量筛选和组合，很大程度地提升草地AGB 估测模型的精度。随机森林模型可以组合不同含义的变量特征，且有效解决“过饱和”和共线性的问题［35］。有很多学者利用随机森林模型进行了农作物、森林、竹林的生物量估算和茶园提取［36-37］，均取得了不错的研究成果。但是，RF 模型也有其局限性，特别是在构建回归分类树的方式上，会低估超出训练集范围内的高生物量值［32］。此外，机器学习基于数据驱动，通常需要大量的数据，若研究的数据量较小，可能对模型准确性造成一定的影响［38］。

Landsat-8 OLI 遥感数据因其相较于低分辨率卫星具有更高的分辨率，相较于之前的Landsat 卫星具有更多的波段和覆盖范围，有避免大气吸收特征的干扰，可用于海岸带观测和云检测，近红外波段、短波红外波段与MODIS 对应波段接近等优势，在我国草地资源的监测方面有较广泛的应用［22］。本研究通过机器学习算法和Landsat-8 OLI 数据提升了传统草地生物量估算模型的精度，但是还存在一定的局限性。Landsat-8 OLI 的空间分辨率是30 m，比起MODIS 的影像空间分辨率（250~1000 m）高很多，但其时间分辨率是16 d。另外，研究区地处青藏高原，受云量和当地天气影响较大，卫星过境重访周期长，不能很好地获得高质量的影像，可能会对模型精度和AGB 动态监测造成一定的影响。在未来研究中，可以获取更高时间分辨率的Landsat 数据，进一步提高模型的精度。另外，草地生物量受地理位置、气候、土壤、地形等多种因素的影响［39］。其中气候因素包括光照、气温和降水等；土壤因素包括土壤营养元素、土壤结构和肥力等。并且还受到草地类型与分布、物种多样性等生物因素，以及放牧、围栏封育、轮牧等管理因素的共同影响［40］。未来可以考虑土壤、地形、气象因素等变量参与模型的构建，优化现有RF 模型。

4 结论

本研究基于Landsat-8 OLI 遥感数据，分别构建并比较了29 个植被指数与草地AGB 的单因素回归模型和RF模型，主要得出以下结论：1）综合并确立了29 个国内外广泛应用的植被指数与高寒草地AGB 的线性最优关系，表明多种植被指数与草地AGB 之间均具有显著的相关性。2）多种植被指数之间具有一定的互补性，机器学习算法可以更好地进行变量筛选和组合，因此基于机器学习的多种植被指数的组合应用在很大程度上可以提升草地AGB 模型的反演精度，从而更加精准地反映草地生物量的时空变化状况。3）得到了2019-2021 年门源县30 m分辨率的草地AGB 分布。空间分布特征为西北部较高，东南部相对较低；大体呈中部高，四周低的状况，其中高寒草甸分布最多。