王亚晖，唐文家，李森*，赵鸿雁，谢家丽，马超，颜长珍

（1.中国科学院西北生态环境资源研究院，甘肃兰州 730000；2.中国科学院大学，北京 100049；3.青海省生态环境监测中心，青海 西宁 810007）

草地是地球上重要的生态系统，约占全球陆地面积的30%［1］。有研究表明，世界草地生态系统每年产生9.06×1011美元的生态价值，占所有生态系统服务总价值的2.72%［2］。我国草地资源丰富，总面积约4×106km2，占土地面积的41.7%［3］，因此草地生态系统质量变化对我国生产、生活和生态等影响巨大。青海省位于青藏高原东北部，是我国重要的天然草地分布区，同时也是我国重要的放牧区［4］。自1980年以来，由于过度放牧、草地开垦以及采挖虫草等因素，青海省草地发生不同程度退化［5］。

草地净初级生产力（net primary productivity，NPP）作为表征草地生产能力的重要指标［6］，其时空变化和驱动机制已经成为草地相关研究的重要内容［7-8］。有研究［9］表明，草地生产力变化受到气候变化和人类活动的影响。为定量区分气候因素和人类活动对草地生产力的影响，相继提出了残差趋势分析法（RESTRENT）［10-11］、分段残差分析法［12］、人类NPP法［13-15］等方法，用以分离人类活动对草地的影响，但是结果对驱动因素的构成和空间差异表达并不直观。为解决这一问题，Leroux等［16］和王巨［17］先后构建多层判断的驱动因素分类框架，分别定量分析非洲萨赫勒地区和河西地区植被变化驱动因素的空间结构差异，结果表明该方法在降水主导植被生长的干旱半干旱区［18］具有良好的适用性。在植被与温度相关性更强的高寒地区［19-20］，该方法能否定量区分植被变化的主导因素，需要深入研究。

为阐明近年来青海省草地生产力变化趋势，分析其驱动机制的空间异质性，本研究基于2001-2017年MODIS NPP时间序列数据产品，分析青海省草地NPP时空变化趋势，并用优化后的驱动因素分类决策树，量化草地NPP驱动因素的空间差异。

1 材料与方法

1.1 研究区概况

青海省地处89°35′-103°04′E和31°30′-39°19′N（图1），东邻四川，南接西藏，西靠新疆，北连甘肃，是青藏高原重要组成部分，总面积为72.23×104km2；地势呈西高东低、南北高中间低，平均海拔大于3000 m；气候为典型的高原大陆性气候，全省年均气温在-5.1～9.0℃之间，绝大部地区的年均降水量低于400 mm［21］。

图1 研究区地理位置与草地分布Fig.1 The location and grassland distribution of the study area

青海省天然草地面积41.92×104km2，其中74.70%的草地存在不同程度的退化［22］。为扭转草地退化趋势，国家和地方政府先后实施了退耕还草、退牧还草、三江源生态保护和建设工程等生态工程［23］，建立了三江源国家公园、祁连山国家公园等11处生态保护地。近年来，随着高原气候的暖湿化，加之多项生态工程的实施和自然保护区的建立，草地质量呈现好转趋势［24］。

1.2 主要数据和预处理

1.2.1 草地时空数据 青海省草地数据来源于中国土地覆被数据集（ChinaCover）［25］，包括2000、2005、2010和2015年等4期30 m空间分辨率的数据，其草地类型有温带草甸、温带典型草原、温带荒漠草原、高寒草甸、高寒草原、高寒荒漠草原等。为减少土地覆被类型变化带来的NPP突变影响，将4期数据均为草地的区域作为研究区，并重采样为500 m×500 m分辨率。经处理后，研究区草地面积为38.80×104km2。

1.2.2 MODIS产品数据集 MOD17A 3HGF.006数据集由美国国家航空航天局（National Aeronautics and Space Administration，NASA）生产分发，目前提供有2001-2020年全球500 m×500 m分辨率的年尺度NPP数据。MODIS NPP产品的反演算法相对成熟［26］，其精度经过验证，已用于区域性植被生产力的相关研究［27-28］。本研究使用的青海省MOD17A 3HGF.006 500 m NPP数据经由Google Earth Engine（https://earthengine.google.com）平台导出下载，其时间尺度与气象数据年限保持一致，为2001-2017年。

1.2.3 气象数据 气象数据集为青藏高原及周边地区气温和降水格点数据（1998-2017年）［29］，来自国家青藏高原科学数据中心（https://data.tpdc.ac.cn/zh-hans/）。该数据集包括1998-2017年青藏高原逐年的平均气温（0.1℃）和降水量（0.1 mm）格网数据，空间分辨率为1 km。为与NPP数据集分辨率一致，将该数据集重采样为500 m×500 m，用以表示气温和降水因素。

1.3 研究方法

1.3.1 趋势分析 Theil-Sen斜率估计（又称Sen斜率估计）和Mann-Kendall检验法（MK检验）是无须数据正态分布的非参数估计与检验方法，对测量误差和小的离群点不敏感［30-31］。因此，Sen-MK趋势分析在长时间序列数据趋势分析中被广泛使用［32］。本研究NPP时间序列的Sen斜率计算公式如下：

式中：NPPi和NPPj为NPP时间序列数据中第i年和第j年的NPP值；Medain（）为取中值函数；β为所有数据对斜率的中值，即时间序列的Sen斜率。当β>0时，时间序列呈增加趋势；β<0时，时间序列呈减少趋势。时间序列趋势的显著性按不考虑值重复（NPP数据在年际间很难出现重复）和时间序列长度n≥10（本研究n=17）时的MK方法检验［33］，其统计量S的计算公式为：

式中：x i和x j分别是时间序列的第i和第j项，sgn（θ）为符号函数，定义为：

当n≥10时，对S标准化后的检验统计量Z可按以下公式计算：

根据Sen-MK趋势分析的显著性水平|Z|>1.96（P<0.05）和|Z|>2.58（P<0.01），将NPP变化趋势分为5类（表1）。

表1 草地NPP变化趋势分级标准Table 1 Classification of gr assland NPP tr ends

1.3.2 相关性分析 相关性分析常用来衡量要素间关系的密切程度。Pearson相关系数是对两个要素x与y相关程度的度量，定义为：

式中：Rxy为要素x与y的Pearson相关系数，x i和yi分别表示第i组样本中x和y的值，和分别为要素x与y的平均值，n为样本量。

在多元相关分析中，偏相关系数可以控制其他变量的影响，突出所研究的两个变量间的真实关系。在有3个变量时，控制第3变量后其他两个变量间的偏相关系数为：

式中：r12·3为变量1和变量2的偏相关系数，r12、r13、r23分别为3个变量中两两间的Pearson相关系数。

基于本研究数据的自由度15（n-2），在95%置信度水平上，即|Rxy|>0.482时，两个要素之间显著相关。而偏相关的显著性检验采用标准的t检验（双边），置信度水平同样为95%。

1.3.3 残差趋势分析 在未经过相关性检验确定NPP与年均气温、年降水量同时相关之前，使用一元回归模型拟合NPP与单气候因素间的关系；而在确定NPP与两个气候因素都相关之后，使用多元回归模型拟合其关系。具体的回归模型如下：

式中：a和b为一元回归模型的系数，c、d和e为多元回归模型的系数，W为某一种气候因素，P为年降水量，T为年均气温。用OLS方法计算两种回归模型的系数，并用标准的F检验衡量模型的显著性，置信度水平取95%。

模型通过检验后，基于模型自变量和模型参数获得NPP预测值，用以表示单一气候因素或降水+气温对NPP的影响。而NPP残差是NPP实际值与NPP预测值之间的差值，用来表示自变量以外的其他因素的影响。当NPP残差序列在时间上有显著的趋势时，认为NPP受回归模型中自变量以外的因素影响；否则认为NPP的变化仅受回归模型内部自变量的影响。本研究用扣除气候因素影响后的NPP残差表示人类活动对NPP的影响部分。为提高非线性趋势的检测，NPP残差时间序列的趋势显著性同样使用Sen-MK方法检验。

1.3.4 基于决策树的驱动因素分析 在自然条件下，除气候因素外，地形、土壤等自然因素，在短期内难以发生较大变化，因此假定引起草地NPP变化的自然因素主要为气温和降水两个气候因素，而其他变化是由人类活动影响造成的。基于上述假定，综合利用Sen-MK趋势分析、相关性分析、残差趋势分析，构建影响草地NPP变化的主导因素的分类决策树框架（图2）。通过该框架，可结合草地NPP趋势分析结果，将人类活动、气温和降水3个驱动因素组合为人类活动、气温、降水、人类活动+气温、人类活动+降水、气温+降水、人类活动+气温+降水等7种主导因素类型，并检验各区域草地NPP显著变化是何种主导因素类型造成的。

图2 草地NPP显著变化的主导因素类型分类决策树Fig.2 The decision tree to distinguish combination types of dominant factors driving significant grassland NPP changes

同时，基于草地NPP变化趋势及其主导因素类型空间分类结果，可以从决策树分类过程中获取人类活动、气温和降水3个因素分别对草地NPP变化的作用类型：在气候因素主导或参与主导草地NPP显著变化的区域，叠加相关性分析结果和草地NPP变化趋势，获得气候因素对草地NPP变化的4种作用类型，包括正相关促进、负相关促进、正相关抑制、负相关抑制。其中，正相关促进表示气候因素的升高（增加）导致草地NPP增加；负相关促进表示气候因素的降低（减少）导致草地NPP增加；正相关抑制表示气候因素的升高（增加）导致草地NPP减少；负相关抑制表示气候因素的降低（减少）导致草地NPP减少。在人类活动主导或参与主导草地NPP显著变化的区域，根据Sen-MK趋势分析和残差趋势分析结果，得到人类活动对草地NPP变化的2种作用类型，包括促进和抑制。其中，促进表示无气候因素参与主导时草地NPP的Sen-MK趋势为增加或扣除气候因素后的NPP残差趋势为增加；抑制则反之。

以上研究方法及其结果中各指标作用区域面积和面积占比等的统计分析，基于Python 3.6，使用gdal和numpy等库实现。

2 结果与分析

2.1 2001-2017年青海省草地NPP时空变化

2001-2017年，青海省草地NPP呈自东向西逐渐减低的趋势，具有明显的地域分异性（图3），除柴达木盆地外，其空间分布与年均气温相对一致。其中，NPP高值区主要分布在青海湖周边和龙羊峡上游黄河干流以东地区；NPP低值区主要分布在可可西里、昆仑山以及柴达木盆地边缘。研究时段内，全省平均草地NPP缓慢波动上升（图4），总体上与平均气温波动趋势近于一致，比平均降水滞后一个阶段。同时期内，青海省草地NPP与气温的关系更密切。

图3 2001-2017年青海省草地多年平均NPP空间分布Fig.3 Spatial distribution of multi-year average grassland NPP in Qinghai Province during 2001-2017

图4 2001-2017年青海省平均草地NPP、平均气温和平均降水的年际变化Fig.4 Inter-annual var iations of aver age grassland NPP，aver age temperatur e and average pr ecipitation in Qinghai Pr ovince during 2001-2017

2001-2017年，青海省70.59%的草地区域NPP变化不显著，主要分布在玉树州东部以及果洛州南部（图5）。NPP呈极显著增加的草地面积为4.61×104km2，占全省草地面积的11.88%，在空间上分布较为集中，主要在可可西里三江源区、扎陵湖和鄂陵湖的北岸及东岸、龙羊峡上下游河段、青海湖北部汇水区以及祁连山山区；这些区域多位于国家公园或者生态保护工程区内，草地生态得到较好的保护和恢复。NPP呈显著增加的草地面积为6.69×104km2，占全省草地面积的17.25%，其空间分布与极显著增加的区域一致，但主要在极显著增加草地的周边，而且更为破碎。NPP极显著减少和显著减少的草地面积分别为0.03×104和0.08×104km2，占全省草地面积的0.08%和0.20%；二者在空间上分布少而集中，主要在治多县与曲麻莱县交界的通天河两侧、黄河源的西北和东南部分区域、玛多县北部边界、久治县境内以及青海湖东南部日月山口附近，而且NPP极显著减少的草地多位于NPP减少草地的内部，草地生产力退化区域的中心比边缘严重。

图5 2001-2017年青海省草地NPP变化趋势空间分布Fig.5 Spatial distr ibution of NPP trends in Qinghai Pr ovince gr asslands during 2001-2017

2.2 草地NPP变化主导因素类型空间分布

青海省11.41×104km2草地NPP显著变化区域中，不同因素类型主导的草地NPP显著变化的面积和空间分布存在明显差异（图6）：草地NPP显著变化的主导因素类型为气温的区域面积最大，其占比为60.66%，主要分布在唐古拉山、祖尔肯乌拉山和可可西里等高寒地区，以及鄂陵湖南岸和青海湖西北汇水区；其中高寒地区，人类活动稀少，温度较低，气温是草地生长的限制性因子，而在鄂陵湖南岸和青海湖西北汇水区，由于水分充足，相对降水来说，草地生长对气温更加敏感。其次是人类活动主导的区域，面积占比为23.45%，主要分布在扎陵湖与鄂陵湖北岸、柴达木盆地西侧、共和盆地南部、湟水谷地和祁连山山区；这些区域主要为人口相对密集，并且人类活动相对强烈的地区。第三是人类活动+气温主导的区域，面积占比为9.49%，主要分布在沱沱河以北到可可西里、哈拉湖东南部、黑河上中游等地区。降水、人类活动+降水和人类活动+气温+降水3种结构类型分别主导的区域面积较小，其占比均不足3.00%，主要在共和盆地、青海湖周边沿岸和柴达木东侧边缘一带。气温+降水主导草地NPP显著变化的面积最小，仅有3.50 km2。

图6 2001-2017年青海省草地NPP显著变化的不同主导因素类型空间分布Fig.6 Spatial distribution of the combination types of dominant factors driving significant grassland NPP changes in Qinghai Province during 2001-2017

2.3 不同驱动因素对草地NPP变化的影响

2001-2017年，气温、降水和人类活动3个因素对青海省草地NPP显著变化的主导区域和作用类型有着明显差异（图7）。

图7 2001-2017年青海省气温（a）、降水（b）和人类活动（c）分别对草地NPP变化影响的空间分布Fig.7 Spatial distribution of effect of temperature（a），precipitation（b）and human activities（c）on grassland NPP in Qinghai Province dur ing 2001-2017

气温主导或参与主导草地NPP显著变化的面积为8.11×104km2。其中，气温升高导致草地NPP增加的面积占99.81%，主要分布在高寒地区；气温对青海省草地NPP的正相关促进作用十分显著。温度降低导致草地NPP增加的面积占0.16%，主要分布在柴达木盆地的东南边缘；这些地区气候干旱、人类活动相对频繁，相对于温度，降水和人类活动对草地生长的限制性更强。而气温变化导致草地NPP减少的区域极少，仅有0.03%。

降水主导或参与主导草地NPP显著变化的面积相对较少，为0.73×104km2。其中，由降水增加导致青海省草地NPP增加的面积占95.96%，主要分布在青海湖北岸、共和盆地至柴达木盆地东部一线以及龙羊峡上下游黄河两侧等海拔相对较低的地区。由降水增加导致青海省草地NPP减少的面积占3.74%，主要分布在久治县境内。由降水减少导致青海省草地NPP增加的面积占0.30%，大多零星分布在唐古拉山附近，并且伴有人类活动因素的主导。而降水对NPP变化的正相关抑制作用并没有被检测到。

人类活动主导或参与主导草地NPP显著变化的草地面积为4.15×104km2。其中，人类活动对草地NPP变化起促进作用的面积占95.40%，在空间上占了人类活动和人类活动+气温两种因素类型主导区域的绝大部分。人类活动对草地NPP变化起抑制作用的面积占4.60%，大于草地NPP减少趋势区域的面积，而在空间上除了与草地NPP减少趋势区域的重叠部分，主要分布在共和盆地南侧与西侧、青海湖北岸、祁连县东部的黑河北侧等地。

3 讨论

本研究基于MODIS NPP数据产品，使用Sen-MK趋势分析，发现2001—2017年青海省29.41%的草地NPP存在显著变化趋势，其空间分布与刘洋洋等［34］的研究结果相一致，但范围更大，这是本研究时间范围更长导致的；与刘旻霞等［35］对青海省未来植被NPP变化的预测相符。然而，2000年以前，青海省草地已经发生严重的退化，草地NPP的显著增加并不表示草地退化的消失，而且70.59%的草地NPP变化不显著表明青海省的草地质量没有完全恢复，仍存在不同程度的草地退化。

将驱动因素分类决策树框架应用到青海省草地NPP变化研究中，为确定气温优先作为第二层考虑因素的可行性，本研究同时以降水优先作为第二层考虑因素进行了计算。结果发现两种方式获得的草地NPP主导因素类型图中有95.21%的相同区域，具有高度的相似性：其中主导因素类型为气温的区域面积分别以60.66%和57.11%的占比均位列第一；其次都是人类活动单一主导的区域，且其面积大小与空间分布完全一致，表明人类活动的主导作用突出时，主导因素类型为人类活动的区域结果对气温和降水的顺序并不敏感；而人类活动+气温主导的区域面积占比仅相差0.32%；人类活动和气温在两种方式的结果中，均表现为引起草地NPP显著变化的主要因素。而两种方式结果有差异的区域，主要分布在青海湖周边和共和盆地部分地区。使用Pearson相关分析和偏相关分析方法分析有差异地区的草地NPP与气温或降水间的关系，结果表明65.00%区域的NPP与气温的相关系数更高；但如果将降水作为决策树框架的优先考虑因素，其中有42.69%区域的气温主导作用被完全忽略。另外，已有学者研究表明气温对青藏高原植被生长的驱动作用强于降水［36］。因此，本研究首先分离气温因素的结果是可靠的，而且优于降水优先的结果。而偏相关分析与多元回归分析的引入，则可以充分利用上一层的判断结果，使各层的分析方法更有针对性。同时，对比金凯等［37］和李传华等［38］使用贡献率区分气候因素和人类活动的方法，两者有所不同：用贡献率可以量化同一空间上3个驱动因素的贡献并进行比较，但不易突出驱动植被变化的主导因素；而驱动因素分类决策树框架建立在整体Sen-MK趋势显著之上，主要探究主导草地生产力发生显著变化区域的驱动因素及其空间上的差异，并在空间上得到一幅较为可靠和更加直观的主导因素类型图，对理解人类活动、气温和降水在草地长期趋势中的驱动力构成有很大帮助。而且，相对于金凯等［37］研究得到的青海省人类活动对植被恢复基本无影响的结果，本研究方法更能体现出气候不可控的情况下，国家和地方政府一系列的生态保护与修复工程对草地恢复的人类活动作用，肯定了生态工程对草地恢复的正面意义。但是本研究未对草地变化不显著区域和短期内的波动变化做深入研究，也未分解多种驱动因素共同主导草地变化区域的单因素分量。对此，可以基于本研究的结果，在不同类型的驱动因素结构中使用适应的贡献率计算方法，以进一步细化各因素间的差异。

人类活动对草地NPP变化起抑制作用的面积为0.19×104km2。但由于气温和降水的影响，这些区域内草地NPP变化并非完全为减少趋势，其中，有41.56%区域的主导因素类型为单一的人类活动，并贡献了草地NPP减少趋势区域的72.36%，主要分布在居民点周边，过度放牧等人类活动对草地的破坏作用掩盖了气候因素的影响；2.94%的区域与气候因素起抑制作用的区域重合，共同贡献了草地NPP减少趋势区域的5.13%；55.42%的区域与气候因素起促进作用的区域重合，在这些区域内人类活动对草地NPP变化的抑制作用小于气候因素对草地NPP变化的促进作用，最终使草地NPP变化仍表现为增加趋势。人类活动对草地NPP变化的促进作用，主要源于对人类破坏草地生态行为的限制和人类对草地生态的主动修复。其中，可可西里地区在建立自然保护区和国家公园后，已退化草地一方面在人类活动减少的影响下缓慢恢复，另一方面受气温升高的影响，草地生产能力得到提高，在主导因素类型分类结果中表现为人类活动和气温的共同主导。扎陵湖和鄂陵湖北部、共和盆地、龙羊峡上游、柴达木盆地东部边缘等地区，是青海省沙漠化十分严重的地区［39-41］，自20世纪70年代以来，前人对这些沙漠化地区做出了众多研究，并实施了一系列有效的沙漠化治理措施［42］，人类活动对草地的恢复作用得到突出体现；这些地区在本研究结果中主要表现为以单-人类活动为主导的草地NPP显著增加。河湟谷地农业发达，人类活动密集，生态政策对人类活动的约束效果明显，该地区退化草地在良好的水热条件下快速恢复，表现为人类活动主导的草地NPP显著增加。

4 结论

1）2001-2017年，青海省有11.41×104km2的草地NPP存在显著变化趋势，占全省草地面积的29.41%。其中NPP极显著增加和显著增加的草地面积占比分别为11.88%和17.25%，主要分布在可可西里、扎陵湖-鄂陵湖流域、祁连山山区、柴达木盆地东侧以及青海湖周边各流域；NPP极显著减少和显著减少的草地分布较少，面积占比分别为0.08%和0.20%，主要在曲麻莱县和久治县境内。17年来，全省草地生产力明显提高。

2）青海省草地NPP显著变化地区的7种主导因素类型中，气温主导的区域最大，占比60.66%；其次是人类活动和人类活动+气温，分别占比23.45%和9.49%。气温和人类活动是推动青海省草地NPP显著变化的最主要的两个驱动因素，而降水的作用只在柴达木盆地东侧、共和盆地、龙羊峡下游黄河两侧等局部区域较为明显；并且三者对青海省草地NPP变化的作用均以促进为主。

3）人类活动对青海省草地NPP变化起促进作用的区域占草地NPP增加趋势区域面积的35.06%。在三江源国家公园与祁连山国家公园等自然保护区和沙漠化治理生态工程区等地区的人类活动正向效果显著，人类活动对生态保护的作用不容忽视。但同时，草地NPP减少趋势区域的77.49%与人类活动密切相关，人类活动仍是青海省草地退化的主要因素。

本研究能够直观可靠地表现出草地发生显著变化区域的驱动因素空间异质性，对植被变化因素分析研究、生态工程效果评价、生态环境精准保护等具有重要意义。但在植被变化趋势不显著和驱动因素构成复杂的地区有待进一步的因素分解研究，以便更好地理解各因素间以及与植被变化间的作用机制。