基于随机森林模型的黄土高原草地净初级生产力时空格局及未来演变趋势模拟

刘欢欢， 陈 印， 刘 悦， 刚成诚

（1.西北农林科技大学草业与草原学院，陕西 杨凌 712100；2.西北农林科技大学水土保持研究所，陕西 杨凌 712100；3.中国科学院水利部水土保持研究所，陕西 杨凌 712100）

草地是陆地生态系统主要组成部分之一，约占陆地面积30%［1］。草地生态系统在维持生态系统功能以及服务社会经济中均发挥重要功能，包括保护生物多样性、水土保持以及减缓温室效应等［2］。植被净初级生产力是衡量生态系统固碳能力的重要指标之一［3］。草地净初级生产力（Net Primary Pro⁃duction，NPP）不仅能够直接反映草地生态系统生产力的高低，而且是评价草地生态系统碳收支及调节草地生态过程的主要因子［4］。因此，科学评估草地NPP及其时空动态对草地生态系统管理及草地资源保护至关重要。

传统野外采样估算草地NPP 的方法虽然结果较为准确，但由于其需要耗费大量的人力、物力和财力，在大空间尺度草地NPP 的估算中较少使用。遥感技术和生态模型的发展为此提供了重要的方法和技术手段。遥感技术通过利用遥感指标数据与草地生物量的相关关系来估算草地NPP，如马文红等［5］基于NDVI 以及地上生物量数据模拟了1982—2006 年中国北方草地生物量空间分布及动态；王莺等［6］利用EVI 与草地干物质产量回归模拟了2006—2008 年甘南草地NPP 变化。模型模拟同样是大空间尺度草地NPP 估算的常用方法，如Guo等［7］利用CASA模型模拟了2001—2018年内蒙古草地NPP；张赟鑫等［8］利用CASA模型和Miami 模型模拟了中亚草地NPP；Zhang 等［9］利用BIOME-BGC 模型模拟了呼伦贝尔草原NPP。近年来机器学习方法在植被生产力估算中亦有广泛的应用。Jung 等［10］利用机器学习中模型树集成算法模拟了全球总初级生产力；Yao 等［11］利用模型树集成算法模拟了1982—2015 年中国0.1°空间分辨率的总初级生产力。随机森林（Random Forest，RF）是基于决策树的集成机器学习算法，目前已被应用于全球［12］、我国西北地区［13］及青藏高原［14］不同类型草地NPP的估算。

黄土高原是世界上最大的黄土覆盖区，气候干旱，水土流失严重，生态环境极其脆弱［15-16］。随着退耕还林还草工程等生态工程的实施，黄土高原植被变化及其环境效应备受关注［17］。目前，在黄土高原草地NPP 时空动态及驱动力分析方面已开展了一定的研究［18-19］。但由于黄土高原地区地形复杂，不同研究结果之间存在较大差异。在全球变暖背景下，黄土高原草地NPP 未来的演变趋势仍不清楚。因此，本研究基于气候、植被、土壤和地形因子，利用随机森林模型模拟了2002—2020 年黄土高原草地NPP 的时空动态，并基于共享社会经济路径（Shared Socioeconomic Pathway，SSPs）4 个未来气候情景数据，估算了黄土高原草地NPP的未来演变趋势，以期为黄土高原草地生态系统管理及可持续发展提供科学依据。

1 研究区概况与方法

1.1 研究区概况

黄土高原地处我国黄河流域中部，地理位置位于100°54′～114°33′E，33°43′～41°16′N，总面积达63×104km2［20］（图1），横跨陕西、甘肃、宁夏、山西、青海、河南以及内蒙古7 个省份，是我国四大高原之一，同时也是世界上最大的黄土覆盖区［15］。黄土高原包括了干旱、半干旱以及半湿润气候区，年降雨量在150～800 mm，海拔在200～3000 m，不同区域在气候上有着明显差异［21］。黄土高原的植被类型以森林、草地以及灌木地为主，草地面积约占总面积60%［22］。

图1 研究区概况及采样点分布Fig.1 Study area and sampling points

1.2 数据获取及处理

通过文献、国家地球系统科学数据中心数据以及GEDI生物量数据（https：//daac.ornl.gov/）获取黄土高原草地地上生物量数据1788个，根据根冠比估算地下生物量，使用草地生物量乘以转换系数（0.475）估算草地NPP（g C·m-2），作为实测NPP数据［23］。

为训练模型，选取19个特征参数包括土壤理化因子（S）、地形因子（T）、气象因子（A）和植被因子（B）4 类。气象因子中温度和降水数据来源于国家地球系统科学数据中心-黄土高原分中心（http://loess.geodata.cn），根据月平均温度数据合成年平均温度（TEM）数据以及年生长季平均温度（TEM4-10）数据，根据月最低温度数据合成为年最低温度（TMN）数据，根据逐月降水量数据合成年降水量（PRE）数据以及年生长季降水量（PRE4-10）数据［24］。数字高程模型（DEM）数据来自NASA 的SRTM 数据（http://srtm.csi.cgiar.org/），坡度（Slope）数据由DEM 数据通过ArcGIS 计算得来。土壤理化因子数据均来自世界土壤数据库（https://data.isric.org/）。日光诱导叶绿素荧光（SIF）数据来自GOSIFv2 数据［25］（https://globalecology.unh.edu），植被有效光合辐射吸收比例（FPAR）数据、归一化植被指数（NDVI）以及蒸散量（ET）来自MODIS 数据产品（https://modis.gsfc.nasa.gov/）。未来降雨量和温度数据来源于WorldClim（https://www.worldclim.org/）。黄土高原空间范围数据以及黄土高原草地分布数据来自中国科学院资源环境与数据中心（https://www.resdc.cn/）空间分辨率为1 km。以上数据分辨率均通过ArcGIS 重采样为1 km，采用D_WGS_1984 投影。本文所有环境因子见表1。

表1 随机森林（RF）中使用的环境因子Tab.1 Environmental factors in RF

SSPs 情景在典型浓度路径（RCPs）情景基础上发展而来，用于定量描述气候变化与社会经济发展路径之间的关系，反映未来社会面临的气候变化适应和减缓挑战［26］。本研究分别选择SSP126、SSP245、SSP370 和SSP585 4 个 情 景 下，ACCESSESM1-5、BCC-CSM2-MR、CanESM5、CMCC-ESM2和CNRM-CM6-1，共计5 个气候系统模式2021—2040年（2030s）、2041—2060 年（2050s）、2061—2080 年（2070s）、2081—2100 年（2090s）4 个时间段的降水量和温度作为输入数据，利用RF 模型预测未来不同气候情景下黄土高原草地NPP 的演变趋势。RF模型中的其他环境因子仍使用2020年数据。

1.3 研究方法

1.3.1 建模方法 本研究采用RF 进行建模，RF 是基于分类和回归树（Classification and Regression Tree，CART）而开发的集成算法，可用于分类问题以及回归问题［27-28］。RF 的主要原理为通过结合多个决策树，并平均其结果使得决策树泛化误差收敛从而产生更好的预测结果［29］。RF 可广泛用于预测问题，并且应用简单，无需复杂地调整参数，通过调整RF 中回归树的数量以及决策树每处节点预测变量的数量来优化模型，提高预测精度［29］。此外，RF 具有处理高维特征空间数据的能力，在遥感应用中已有较多研究［30］。RF 模型性能通过决策系数（R2）以及均方根误差（RMSE）评价其拟合能力以及预测性能。本研究使用R 4.0.5进行随机森林建模、参数优化、预测以及空间可视化。

1.3.2 Theil-Sen Median趋势分析和Mann-Kendall检验 Theil-Sen Median 趋势分析是一种稳健的非参数统计趋势计算方法［31］，计算公式如下：

式中：xj和xi为样本数据值。当Senslope＆gt;0 时，表明该数据在此时间序列呈增长趋势，反之则表示呈下降趋势。

Mann-Kendall检验方法是一种非参数统计检验方法，用于判断趋势的显著性［30］。计算公式如下：

式中：S为Mann-Kendall 检验需要构建的统计量，ZMK是标准正态变量。本文在置信度水平α=0.05上来检验NPP变化趋势的显著性。当P小于0.05时表示NPP变化显著，反之则表示不显著。

2 结果与分析

2.1 NPP与环境因子相关性分析

在植被因子中，NPP 与SIF 相关性最强（r=0.82），其次为NDVI 和FPAR（图2a）。在气象因子中，NPP与ET相关性最高（r=0.91），与PRE4-10、PRE、TMN 相关性也较强，而与TEM 以及TEM4-10相关性较弱。NPP 与土壤因子整体相关性较弱，其中NPP与SAND_S 相关性最强（r=-0.27），与SOC_T 相关性最弱（r=0.11）。整体而言，NPP 与气象因子、植被因子相关性较强，而与气象因子中的TEM 和TEM4-10、土壤因子以及地形因子相关性弱。而在随机森林模型因子重要性中（图2b），对模型估算NPP影响最大的是ET。PRE、PRE4-10以及植被因子中的SIF 对NPP估算影响也较大。结合相关性分析和重要性排序，影响黄土高原草地NPP的环境因子主要为植被因子以及PRE、PRE4-10和ET。

图2 NPP与环境因子相关性（a）及随机森林因子重要性（b）Fig.2 The correlations between NPP and environmental factors(a),and the importance of random forest factors(b)

2.2 RF模型精度评价

2.3 草地NPP空间动态

黄土高原草地NPP 的空间分布呈现较强的异质性（图4）。2002—2020 年黄土高原草地年平均NPP 为276.55 g C·m-2。黄土高原草地NPP 的分布呈现从南向北逐渐递减的趋势，其中黄土高原西北地区草地NPP 最低，主要位于陕西北部、宁夏以及甘肃兰州地区附近，最低值为76.52 g C·m-2；草地NPP 最高的区域主要分布在黄土高原南部地区，如陕西关中地区，最高值为634.86 g C·m-2；而在黄土高原东部地区草地分布较少，草地NPP 约为450 g C·m-2。

图4 2002—2020年黄土高原草地NPP空间分布Fig.4 Spatial pattern of grassland NPP in the Loess Plateau during 2002-2020

2002—2020年，黄土高原草地NPP呈现总体增加的趋势（图5），其中16.59%的地区呈显著增加趋势，38.52%的地区呈轻微增加趋势，主要分布在黄土高原东部以及南部；1.89%的地区草地NPP 变化不显著；35.57%的地区呈轻微减少趋势，7.43%的地区呈显著减少趋势，主要分布在黄土高原西北地区。

图5 2002—2020年黄土高原草地NPP空间变化Fig.5 Spatial dynamic of grassland NPP in the Loess Plateau during 2002-2020

图3 随机森林模型精度验证Fig.3 Accuracy test of Random Forest model

2.4 未来不同气候情景下黄土高原地区草地NPP演变趋势

到本世纪末，4 个未来气候情景下黄土高原草地NPP 均呈现增加的趋势，其中SSP126 情景草地NPP增加幅度最少，为1.20 g C·m-2·（10a）-1；SSP585情景增幅最大，为8.10 g C·m-2·（10a）-1。SSP245 和SSP370 情景，草地NPP 的增速分别为3.90 g C·m-2·（10a）-1和4.00 g C·m-2·（10a）-1（图6）。草地NPP 增加的区域主要集中在黄土高原南部地区，SSP585中草地NPP 在黄土高原南部地区增加最多。而草地NPP 下降的区域主要集中在黄土高原西部地区，其中在SSP126中草地NPP将下降最多。

图6 不同气候情景下2090s黄土高原草地NPP空间分布Fig.6 Spatial pattern of grassland NPP in the Loess Plateau in 2090s under different climate scenarios

3 讨论

3.1 随机森林模型估算草地NPP适用性

相比其他NPP 估算模型，随机森林的优势在于可以更准确地描述NPP与其他环境因子之间的非线性复杂关系，且随机森林有较好抗过拟合能力［32-33］。本研究利用随机森林模拟得到黄土高原草地的平均NPP（276.55 g C·m-2）。刘洋洋等基于CASA模型模拟的2000—2015 年黄土高原草地年均NPP 为202.93 g C·m-2，原因可能与使用草地分布数据不同有关［19］。魏建洲基于MOD17A3 NPP 数据模拟的黄土高原草地NPP 在208.04～262.66 g C·m-2之间，其高值与本研究结果较为接近［34］。Li等利用CASA模型分别模拟了石羊河流域及我国西北干旱半干旱区草地年均NPP，其结果分别为236.99 g C·m-2和277.97 g C·m-2［35-37］。相比于其他研究结果，本研究模拟草地NPP结果偏大，一方面在于不同研究中所使用的草地分布数据不同，如刘洋洋使用的草地分布数据来自Global Land Cover 2000 数据产品［19］，本研究所使用草地分布数据来源于中科院土地利用数据产品；另一方面基于随机森林模型模拟草地NPP，受训练数据影响很大，本研究用于训练模型的实测草地NPP数据主要集中在黄土高原西部，东部地区的样本点数量相对较少。但总体而言，本研究利用随机森林模拟的黄土高原草地NPP 结果较为可信。

3.2 黄土高原草地NPP的动态变化

本研究发现2002—2020年黄土高原草地NPP呈增加趋势，这与前人在该区域研究结果一致［19,38-40］。由于黄土高原气候干旱，降水是主要的植被生长限制因子［18］。2002—2020 年黄土高原降水量呈增加趋势（图7），促进了黄土高原草地NPP 的增加。刘洋洋等［19］也认为黄土高原降水增加导致的湿润环境为草地NPP 增加提供了有利条件。Gang 等［40］的研究结果同样表明气候变化对黄土高原草地NPP增加有着促进作用。从空间变化上来看，黄土高原西部和南部草地NPP增加，与刘铮等［38］的研究结果相似。这些地区位于半干旱以及半湿润地区，降水量的增加导致了NPP的快速积累［41］。此外，自1999年以来，退耕还林还草工程通过围栏封育、封山禁牧等措施使草地生态环境得到了较大改善，是导致黄土高原草地NPP增加的重要因素［40-42］。到本世纪末，4个情景中黄土高原草地NPP仍呈增加趋势，原因在于未来黄土高原地区降水量和最低温度将继续增加，暖湿化的气候条件有利于草地NPP 的积累，尤其是在SSP585情景中。

图7 2002—2020年黄土高原年降水量变化趋势Fig.7 Temporal trend of mean annual precipitation in the Loess Plateau from 2002 to 2020

3.3 RF估算结果的不确定性分析

本研究使用RF 模型模拟黄土高原草地NPP 精度较高，但仍存在一定的不确定性。首先，RF 模型性能受数据集数量及质量影响较大［43］，黄土高原区域地形复杂，本文使用1 km的空间分辨率数据模拟草地NPP 导致部分区域草地NPP 模拟结果精度较低。其次，RF 模型作为经验模型，受训练数据质量以及数量影响较大，本研究草地NPP实测数据空间分布不均匀，在黄土高原西北地区以及西南地区较多，而其余地方较少，未来应考虑使用更多的实测数据进行模拟。此外，RF 模型属于黑箱模型，通过多次调试来确定合理模拟结果，但无法控制模型内部运行［44］。

4 结论

本研究基于野外实测数据，以植被、气候、土壤和地形因子作为环境因子，利用随机森林模型模拟了黄土高原草地NPP的时空动态，并估算了未来不同气候情景下，草地NPP 的演变趋势，得出了以下主要结论：

（1）2002—2020 年，黄土高原草地平均NPP 为276.55 g C·m-2·a，空间上呈“东南高西北低”的分布特征。近19 a 黄土高原草地NPP 整体呈现上升趋势，主要得益于20世纪以来黄土高原地区降雨量的增加及“退耕还林还草”等一系列生态工程。

（2）到本世纪末，在暖湿化的气候背景下，黄土高原草地NPP 均将继续保持增加的趋势，其中在SSP585情景中将增加最多。

（3）基于气候、环境、土壤和地形因子构建了黄土高原草地NPP随机森林模型，尚未包含人为活动指标因子，构建综合气候、环境、土壤和人为因素等多因子的要素体系，将会进一步提高草地植被生产力的估算精度。