青藏高原退化草地的恢复潜势研究

王瑞泾，冯琦胜，金哲人，刘洁，赵玉婷，葛静，梁天刚

（兰州大学草地农业生态系统国家重点实验室，兰州大学农业农村部草牧业创新重点实验室，兰州大学草地农业教育部工程研究中心，兰州大学草地农业科技学院，甘肃 兰州 730020）

草地生态系统是全球范围内分布最为广泛的生态系统之一，约占地球陆地面积的25%［1］。同时，草地也是一种十分重要的自然资源，在气候调节、防风固沙、生物多样性保育、水土保持和维护生态平衡等方面具有不可替代的作用。如今世界各地的草地都面临着退化，尤其是以高寒草地（约占比49%）为主要草地类型的青藏高原［2］，虽然近年来青藏高原高寒草地生态系统整体上呈现改善的状态，但仍有部分草地存在不同程度的退化［3］。草地生态系统恢复的关键是草地植被的恢复［4］，目前青藏高原草地恢复措施主要有施肥、补播、围栏封育、人工草地建植等。近年来，近自然恢复（close-to-nature restoration）理念受到了广泛关注并在植被恢复中得到了有效利用，近自然恢复并不是完全摒弃传统人工恢复，而是借助传统人工恢复措施，主要依靠生态系统的自我调节进行可持续生态恢复，这也是一种较为切合高寒草地恢复的措施［5］。但青藏高原高寒草地近自然恢复潜势的大小及其空间分布格局尚不明确。因此，开展青藏高原草地恢复潜势研究对青藏高原草地生态系统保护与恢复具有重要意义。

近年来，关于恢复潜势的研究已有一些进展，Venter 等［6］用Sen 斜率和Mann-Kendall 检验计算1986−2019年的增强植被指数（enhanced vegetation index，EVI）变化趋势，分析南非的植被退化与恢复潜力；高海东等［7］与赵广举等［8］通过地形、土壤、植被、气候、坡度等因素叠加分区，基于“生境越相似的区域，植被恢复潜力越接近”的相似生境法原则，以某区内植被覆盖度的最大值评估黄土高原的恢复潜力；李海东等［9］采用层次分析法和模糊综合评价方法，建立了高寒河谷沙地植被恢复潜力综合评价模型来研究不同类型沙地植被的恢复潜力。潘竟虎等［10］基于潜在归一化植被指数（normalized difference vegetation index，NDVI）数据，采用改进的CASA 模型模拟得到中国潜在植被净初级生产力。

植被净初级生产力（net primary productivity，NPP）可以较为准确地反映草地植被的生长状况，且对气候变化与人类活动的影响较为敏感［11］，适合作为评价草地恢复潜势的指标。本研究拟在估计青藏高原草地现实与潜在NPP 的基础上，将恢复潜势定义为草地恢复现状与恢复顶级之间的“距离”，依此分析青藏高原高寒草地恢复现状和恢复潜势。

基于此，本研究拟通过CASA 模型估计现实NPP（actual net primary productivity，NPP），通过Thornthwaite Memorial 模型估计潜在NPP（potential net primary productivity，PNPP），采用最大值合成法将19年来的年PNPP进行合成得到PNPP 最大值（maximum potential net primary productivity，PNPPm），用PNPPm与NPP 之间的差值评估草地的恢复潜势，定量研究2001−2019年青藏高原草地恢复潜势及其演变趋势。本研究结果可以判定不同草地恢复现状、恢复潜力与恢复价值，为青藏高原退化草地恢复政策的制定提供科学与理论支撑，对指导青藏高原草地保护与生态系统恢复具有重要意义。

1 材料与方法

1.1 研究区概况

青藏高原是我国的重点牧区之一，约占我国草地总面积的1/3。介于26°00′12″−39°46′50″N，73°18′52″−104°46′59″E，涉及青海、西藏、新疆、四川、甘肃和云南6 个省区共201 个县级行政区，面积约为2.5724×106km2，约占中国陆地总面积的26.8%。青藏高原草地资源丰富，草地类型众多，其中包括高寒草原、高寒草甸、高寒草甸草原、高寒荒漠等17 种草地类型［12］（图1）。

图1 青藏高原草地类型分布Fig.1 Distribution of grassland types on the Qinghai-Tibet Plateau

1.2 数据来源与预处理

1.2.1 气象数据 2001−2017年的逐月降水和温度插值数据采用的是Ge 等［13］研究中使用的数据。2018−2019年的逐月降水和温度插值数据来自国家科技基础条件平台−国家地球系统科学数据中心（http：//www.geodata.cn）。

1.2.2 太阳总辐射量数据 在ArcGIS 软件中，通过Area Solar Radiation 工具计算得到2001−2019年青藏高原每月太阳总辐射量（solar radiation，SOL）。该过程使用到的数字高程模型（digital elevation model，DEM）数据来源于NASA 的航天飞机雷达地形测绘任务数据（shuttle radar topography mission，SRTM）（http：//srtm. csi.cgiar.org/），空间分辨率为90 m［14］。

1.2.3 MODIS 数据 本研究使用了来自美国航空航天局对地观测系统数据与信息系统（NASA’S earth observing system data and information system，EOSDIS）开发的MOD13A3 和MCD43A4 产品。MOD13A3 产品提供了2001−2019年逐月NDVI 数据，空间分辨率为1 km；MCD43A4 产品提供了2001−2019年逐日经过BRDF 校准的地表反射率数据，空间分辨率为500 m，它包括MODIS1～7 波段，本研究仅使用第2 和6 波段计算所得地表水体指数（land surface water index，LSWI）作为CASA 模型的输入变量计算NPP。

对于以上数据本研究使用MODIS 数据重投影工具（MODIS reprojection tool，MRT）进行了格式转换、拼接与转投影，处理后的数据为Geo-Tiff 格式，采用WGS 1984 Albers 投影。并将所有数据重采样至1 km，裁剪出青藏高原地区的影像数据。与此同时，本研究还采用最大值合成法（maximum value composite，MVC），合成了逐月地表水指数（land surface water index，LSWI）数据，用于后续研究［15］。

1.3 现实净初级生产力（NPP）

本研究采用CASA 模型对2001−2019年青藏高原的NPP 进行估算，该模型是当前应用最广泛的净初级生产力估计模型［16］。该模型主要基于光合有效辐射（absorbed photosynthesis active radiation，APAR，MJ·m−2）和光能利用率（ε，g C·MJ−1）估计NPP，其公式为：

式中：SOL 为太阳总辐射量（MJ·m−2），是用数字高程模型（DEM）通过ArcGIS 软件计算得出［14］；FPAR 为植被冠层对光合有效辐射的吸收比例，由各植被类型的NDVI 及其最大值、最小值，比值植被指数（SR）及其最大值、最小值以及FPAR 最大值、最小值确定，计算过程参照朱文泉等［16］的研究；ε 为光能转化率，由月平均气温、植物生长最适月平均气温、LSWI 以及理想条件下最大光能转化率（ε*，本研究选取全球通用值0.389 g C·MJ−1）得出，计算过程参考Xu 等［17］的研究。

1.4 潜在净初级生产力（PNPP）

本研究采用经验证的Thornthwaite Memorial 模型对2001−2019年青藏高原的PNPP 进行估算，该模型是在Miami 模型的基础上对Thornthwaite 潜在蒸散发模型加以修正而建立的。是利用蒸散与碳吸收的关系，以降水和温度为输入，来估计PNPP，这一方法比Miami 模型更加精确［18］。其表达式为：

式中：PNPP 为年潜在净初级生产力（g C·m−2）；v为年平均实际蒸散量（mm）；L为年平均蒸散量（mm）；t和r分别表示年平均温度（℃）与年降水量（mm）。年平均温度（t）和年降水量（r）使用2001−2019年研究区降水和温度空间插值数据［13］，空间分辨率为1 km，数据为Geo-Tiff 格式，采用WGS 1984 Albers 投影。

1.5 恢复潜势

本研究认为天然草地的恢复潜势（R）可以定义为潜在净初级生产力的最大值（PNPPm）与现实净初级生产力之间的差值，PNPPm是基于2001−2019年PNPP 数据，采用最大值合成法得到的，即：

1.6 趋势分析

将Theil-Sen Median 趋势分析、Mann-Kendall 检验与Hurst 指数相结合，可以评估2001−2019年青藏高原NPP 变化趋势分布特征。在进行Theil-Sen Median 趋势分析时，由于几乎不存在SNPP完全等于0 的区域，所以本研究根据实际情况，认为SNPP＜−2 的区域为退化区域；−2＜SNPP＜2 的区域为稳定不变的区域；SNPP＞2 的区域为改善区域。本研究将Mann-Kendall 检验在0.05 置信水平上的显著性检验结果划分为趋势显著变化（Z＞1.96 或Z＜−1.96）和不显著变化（−1.96≤Z≤1.96）。将Hurst 指数小于0.5 的区域定义为不可持续性变化，Hurst 指数大于0.5 的区域为可持续性变化。具体计算方法参考袁丽华等［19］的研究。

2 结果与分析

2.1 NPP 的时空变化特征

利用2001−2019年的年NPP 数据计算得到多年平均NPP 空间分布图（图2）。从图2 可以看出，青藏高原NPP 表现为东南部的甘孜藏族自治州、阿坝藏族羌族自治州以及甘南藏族自治州等地明显较高，西北部的阿里、那曲以及海西西部等地区偏低，由此可见青藏高原草地NPP 呈现出由西北向东南逐渐增加的趋势，即青藏高原东南部草地状况较西北部更好。

图2 2001-2019年青藏高原年平均NPP 空间分布Fig.2 Spatial distribution of annual mean NPP over the Qinghai-Tibet Plateau from 2001 to 2019

基于2001−2019年NPP 数据，用Theil-Sen Median 趋势分析、Mann-Kendall 检验与Hurst 指数进行分析，得到了2001−2019年NPP 的变化趋势与可持续性耦合信息（图3 和表1）。

图3 2001-2019年青藏高原NPP 变化特征空间分布Fig.3 Spatial distribution of NPP variation characteristics over the Qinghai-Tibet Plateau from 2001 to 2019

表1 NPP 变化趋势统计Table 1 Statistics of NPP trend（%）

从图3 和表1 中可以看出，耦合共分为6 种情况，即未来趋势不确定、持续性轻微退化、持续性显著退化、持续性稳定不变、持续性轻微改善以及持续性显著改善。除去无法预测未来变化趋势的区域，青藏高原草地NPP 状况改善的区域明显多于退化区域。其中，无法预测区域占36.21%，持续恢复区域占40.98%，持续稳定不变区域占12.72%，而持续退化区域仅占3.47%。相对而言，退化较为严重的草地类型为热性灌草丛和温性草甸草原。而青藏高原中面积较大的草地类型−高寒荒漠、高寒草甸草原、高寒草甸、温性草原化荒漠以及温性荒漠均呈现出较为明显的恢复趋势。由此可见，整个青藏高原天然草地还是以可持续的恢复状态为主。

为了分析不同草地类型的NPP年际变化，以便更好地了解青藏高原NPP 的变化情况，本研究制作了各草地类型NPP年际变化图（图4）。从图4 可以看出，在2001−2019年间大部分草地类型年NPP 都呈现波动上升趋势，且在2018年有大幅增加，多类草地增幅达到100 g C·m−2，这可能与当地气候变化或政策制定有关。

图4 2001-2019年青藏高原各草地类型NPP年际变化Fig.4 Interannual variation of NPP of different grassland types in the Qinghai-Tibet Plateau from 2001 to 2019A：热性灌草丛类Thermal shrub tussock；B：热性草丛类Thermal tussock；C：暖性灌草丛类Warm-temperate shrub tussock；D：暖性草丛类Warmtemperate tussock；E：低地草甸类Lowland meadow；F：山地草甸类Mountain meadow；G：温性草甸草原类Temperate meadow steppe；H：温性草原类Temperate steppe；I：温性荒漠草原类Temperate desert grassland；J：温性草原化荒漠类Temperate steppe desert；K：温性荒漠类Temperate desert；L：高寒草甸类Alpine meadow；M：高寒草甸草原类Alpine meadow steppe；N：高寒草原类Alpine grassland；O：高寒荒漠草原类Alpine desert grassland；P：高寒荒漠类Alpine desert；Q：沼泽类Marsh.下同The same below.

2.2 PNPPm的空间分布特征

为了直观地反映青藏高原2001−2019年PNPP 的空间分布以及比较各草地类型的具体情况，基于2001−2019年PNPP 数据，用最大值合成法制作了PNPPm的空间分布图（图5），从图中可以看出，PNPPm在青藏高原的分布呈现出明显的南高北低的态势，其中日喀则地区、阿坝藏族羌族自治州、甘南藏族自治州和甘孜藏族自治州等地是青藏高原PNPPm较高的几个地区，这与NPP 的空间分布相近。统计不同草地类型的PNPPm（图6），可以看出热性灌草丛、热性草丛、暖性灌草丛和暖性草丛是青藏高原中PNPPm较大的4 种草地类型，其PNPPm约为1300～1500 g C·m−2，而低地草甸、温性荒漠草原、温性草原化荒漠、温性荒漠、高寒荒漠草原以及高寒荒漠PNPPm较小，约为250～750 g C·m−2，其他草地类型PNPPm在1000 g C·m−2左右。这表明在仅考虑气象因素的情况下，青藏高原东南与西南部草地会具有更好的生长状况，热性灌草丛、热性草丛、暖性灌草丛和暖性草丛应为长势较好的4 种草地类型。

图5 2001-2019年青藏高原PNPPm空间分布Fig.5 Spatial distribution of PNPPm over the Qinghai-Tibet Plateau from 2001 to 2019

图6 2001-2019年青藏高原各草地类型PNPPm统计Fig.6 PNPPm statistics of grassland types in the Qinghai-Tibet Plateau from 2001 to 2019

2.3 恢复潜势的空间分布特征

基于2001−2019年的PNPPm与NPP 数据得到恢复潜势分布图，并分析了19年平均恢复潜势的空间分布（图7）。从图7 可以看出，恢复潜势的分布与PNPPm的分布情况相似，东南与西南两部分较高，北部普遍偏低。其中日喀则地区与甘南地区恢复潜势最大，而阿里北部以及那曲北部等地区恢复潜势较低。因此在日喀则地区与甘南地区开展草地恢复与保护的价值最大，前景最优。

图7 2001-2019年青藏高原恢复潜势空间分布Fig.7 Spatial distribution of recovery potential over the Qinghai-Tibet Plateau from 2001 to 2019

基于2001−2019年青藏高原恢复潜势数据，分草地类型制作了19年平均恢复潜势的柱状图（图8）。从图中可以看出，热性灌草丛、热性草丛、暖性灌草丛、暖性草丛、山地草甸以及沼泽是青藏高原恢复潜势较大的6 种草地类型，而低地草甸、温性草原化荒漠、温性荒漠、高寒荒漠草原、高寒荒漠、温性荒漠草原则恢复潜势较小，这一结果与PNPPm的情况相近，在开展青藏高原草地恢复时应优先考虑热性灌草丛、热性草丛、暖性灌草丛、暖性草丛、山地草甸以及沼泽这6 种草地类型。

图8 2001-2019年青藏高原各草地类型恢复潜势统计Fig.8 Recovery potential statistics of grassland types in Qinghai-Tibet Plateau from 2001 to 2019

2.4 气温、降水与恢复潜势的关系

本研究基于2001−2019年的年均降水与年均温数据得到青藏高原降水与气温空间分布图，并分析了它们与恢复潜势的关系（图9 和图10）。从图中可以看出，青藏高原气温较高的区域为东部、南部以及西南部，降水则具有明显的由东南向西北逐渐降低的态势，这与青藏高原西南与东南两部分恢复潜势较高的情况不谋而合，这说明青藏高原恢复潜势或许受到降水和温度的影响，且温度的影响可能较降水更大一些，这与德吉央宗等［20］研究得出的结论一致。

图9 2001-2019年青藏高原气温空间分布Fig.9 Spatial distribution of temperature over the Qinghai-Tibet Plateau from 2001 to 2019

图10 2001-2019年青藏高原降水空间分布Fig.10 Spatial distribution of precipitation over the Qinghai-Tibet Plateau from 2001 to 2019

2.5 恢复潜势的变化趋势

基于2001−2019年的青藏高原恢复潜势数据得到青藏高原恢复潜势变化趋势图，并分析了青藏高原恢复潜势的变化趋势以及成因（图11）。从图中可以看出，青藏高原草地恢复潜势提升的区域主要集中在西南和东北部，尤其是西南部，降低的区域主要为青藏高原中部，而青藏高原恢复潜势整体为提升、降低与不变三者均有的局面。青藏高原恢复潜势变化趋势的这种空间分布格局也可以作为研究青藏高原草地恢复价值的一个标准，即西南部与东北部草地恢复价值更大。

图11 2001-2019年青藏高原恢复潜势变化趋势Fig.11 Change trend of recovery potential of Qinghai-Tibet Plateau from 2001 to 2019

3 讨论

2001−2019年青藏高原NPP 呈现东南高、西北低的空间分布特征，且各草地类型NPP 均呈现波动上升的趋势，这与陈舒婷等［21］的结论一致。本研究基于CASA 模型计算NPP 时，光能转化率（ε）选用了全球通用值0.389 g C·MJ−1，然而，目前对于光能转化率的看法不一，如王保林等［22］认为我国草原最大光能转化率的取值区间为0.608～1.000 g C·MJ−1；Raymood 等［23］认为其最大值为3.5 g C·MJ−1。本研究的取值可能会造成对NPP 的低估，因此还需在后续研究中进一步探讨更适合估算青藏高原NPP 的最大光能转化率。

本研究使用Thornthwaite Memorial 模型估算青藏高原2001−2019年PNPPm为822.70 g C·m−2。潘竟虎等［10］利用改进的CASA 模型对全国植被PNPP 进行了模拟，其均值为468.94 g C·m−2；吕振刚等［24］采用改进的CASA 模型估算2001−2010年中国华北落叶松的PNPP，其均值为342.7 g C·m−2；任正超等［25］采用改进的综合顺序分类系统以及改进的CASA 模型模拟1982−2012年中国自然植被PNPP，其均值为278.7 g C·m−2。本试验的研究结果较上述结果偏大，最主要的原因是本研究采用最大值合成法合成了19年年最大PNPP 得到了本研究所使用的PNPPm；最后本研究采用Thornthwaite Memorial 模型，而前人的研究多采用改进的CASA 模型，模型的选择差异也是导致研究结果有差距的重要原因。因此，在之后的研究中还可以探讨其他计算PNPP 的模型对于青藏高原草地恢复潜势的计算精度，如改进的CASA 模型、综合自然植被净第一生产力模型、NPP 分类指数模型等。

本研究认为恢复潜势应为某地有可能达到的最佳恢复程度，即可能达到的PNPP 与NPP 之间的差值，因此计算恢复潜势时使用了2001−2019年最大值合成的PNPPm，而前人研究大多认为恢复潜势为年PNPP 与年NPP之间的差值，这会导致本试验的PNPPm较其他相关试验偏大，恢复潜势也会因此偏大，但这并不影响研究区恢复潜势的空间分布以及得出的青藏高原各地恢复价值的大小以及恢复措施。

本研究发现青藏高原恢复潜势呈现明显的西南部最高、东南部次之而北部偏低的分布情况，分析认为形成这种分布格局的原因可能与青藏高原的气温、降水、海拔以及人类活动有关，在今后的研究中应着重探讨青藏高原恢复潜势分布的成因并采用多种方法估算青藏高原恢复潜势，减少系统误差。

4 结论

本研究采用PNPPm与NPP 之间的差值作为恢复潜势，着重探讨青藏高原天然草地NPP、PNPP 与恢复潜势的空间分布与时空变化，以期为指导青藏高原草地恢复提供理论支持，主要结论如下：

1）2001−2019年青藏高原草地NPP 由西北向东南逐渐增加，青藏高原草地整体以可持续的恢复状态为主，其中热性草丛、低地草甸、温性荒漠以及高寒荒漠草原等草地类型恢复态势明显。

2）2001−2016年间青藏高原草地NPP 波动上升，上升趋势并不明显，而在2016−2018年期间青藏高原各草地类型NPP 大幅增加。

3）2001−2019年青藏高原草地恢复潜势呈现西南与东南部较高，北部偏低的分布情况。由此可见日喀则地区、阿里南部地区、阿坝藏族羌族自治州以及甘南藏族自治州等地草地具有较高的恢复价值，在这些地区开展草地恢复的前景更好。