北方农牧交错带植被NPP的时空变化及其驱动因子分析

薛晓玉, 王晓云, 段含明,2, 颉耀文

(1.兰州大学 资源环境学院, 兰州 730000; 2.西华师范大学 国土资源学院, 四川 南充 637002)

近年来，全球范围的气候异常变化、人类活动不断增强，这些问题正影响着陆地生态系统的平衡，导致区域生态环境的破坏，直接或间接地影响到植被NPP变化，使其变化驱动的研究显得尤为重要[1-2]。NPP不仅可以直接表征陆地生态系统的质量状况，而且在保持生态系统平衡和调节区域生态环境中扮演着重要的角色[3]。净初级生产力(Net Primary Productivity，NPP)是指植物在单位时间和面积上由光合作用产生的有机物总量减去自身呼吸作用消耗有机物的剩余量，这部分能量是其自身用于生长和繁殖的能量，也是生态系统中其他成员生存和繁衍的物质基础[4]。

植被NPP的估算方法很多，根据模型对输入参数的选择，可将现有模型大体分为3类：气候相关模型、机理模型和光能利用率模型[5-7]。气候相关模型是根据野外获取的采样点实测数据来建立NPP与气候因子的简单统计回归模型，但该模型由于缺乏生态学理论作为依据或只能对潜在NPP进行模拟，估算的准确性不高[5]。机理模型是以气候、土壤和植被等因子作为输入参数，模拟生态系统中的光合作用、呼吸作用和土壤微生物的分解过程，从而模拟植被的能量流动与物质转换，涉及到的输入参数较多，且其中部分参数不易获取，单纯的机理模型发展具有局限性[7]。光能利用率模型是以光合作用原理为依据，以NDVI数据、太阳辐射数据、气温和降水数据及植被类型数据为输入因子，与其他模型相比，对地面实测数据的依赖较小，且模型的输入参数相对较少，易于收集，在区域尺度和大尺度测算NPP的研究以及全球碳循环研究中被广泛采用[8-13]。

NPP的产量不仅仅受到植被本身的影响，而且还会受到外部条件的限制[14]。目前，对于影响NPP的自然因素研究主要集中在气温、降水和地形等方面的研究[15-17]。随着社会生产力的不断发展，人类活动对生态环境产生的影响逐渐表现出全球性，其对植被生长的影响也日益显著，因此，定量分析气候变化和人类活动如何影响植被变化成为当前的研究热点[18-20]。北方农牧交错带作为中国北方典型的农、牧业系统过渡带，地理位置突出，对外界的抵抗力较弱，且对气候变化和人类活动的干扰敏感。20世纪70年代末期以来，国家先后实施了“三北防护林”、“退耕还林”和“京津沙尘暴源治理”等生态工程，出台了草地生态保护补偿奖励政策，努力促使北方农牧交错带的植被恢复和好转[21]。深入研究这一区域NPP时空变化及其驱动因子，不仅可以揭示气候变化、人类活动与植被变化的相互作用机制，还可为改善区域的生态环境提供参考。

1 研究区概况与数据

1.1 研究区概况

本研究参考赵松乔[22]、赵哈林等[23]的研究成果，将北方农牧交错带界定为多年平均年降水量300～450 mm、降水年变率15%～30%和干燥度指数0.2～0.5范围内的区域[21]。为了保持空间的连贯性，并考虑不同时期的北方农牧交错带边界变化的特点，构建了宽50 km的缓冲区，由此形成了本研究所采用的研究区[21]，其介于34°48′—47°19′N，101°43′—126°00′E，面积为83.5万km2，海拔高度为100～4 700 m。该区位于东部季风区与西北大陆性气候区的过渡带，年平均气温为0～9℃，其西北地区干旱少雨，植被稀少；东南地区降水丰沛，气候湿润。由于特殊的地理位置，北方农牧交错带在地貌和气候的双重作用下，生态环境的承载力较小。

1.2 数据来源及预处理

NDVI数据来源于中国科学院资源环境科学数据中心(http:∥www.resdc.cn/)。该产品经过了辐射校正、大气校正、几何校正和多波段合成技术等预处理，采用最大值合成法生成1998年以来的中国月度(1—12月)、季度(春、夏、秋、冬)、年度的植被指数数据集，数据的空间分辨率都为1 km。气象数据为1998—2016年北方农牧交错带及其周围166个站点的日平均气温、日降水量和30个站点的逐月辐射数据，来自于中国气象数据网(http:∥data.cma.cn/site/index.html)。对该数据进行累加统计得到相应站点的月平均气温和月降水量数据，与辐射月值数据均采用ANUSPLIN方法[24]进行空间插值获得与NDVI数据有相同投影和空间分辨率的栅格数据。植被类型数据来源于国家青藏高原科学数据中心(http:∥data.tpdc.ac.cn/zh-hans/)。该数据是基于IGBP分类系统的植被类型数据，生产年份为2000年，空间分辨率为1 km，中国地区的总体分类精度为71%，高于其他数据分类产品，可以作为光能利用率模型的输入数据。本文对以上数据进行了裁剪、投影变换、镶嵌等预处理，作为CASA模型的主要输入数据，最终的分析区域是将研究区内非植被区域进行掩膜处理的植被覆盖区。

2 研究方法

2.1 基于CASA模型的NPP估算及验证

CASA模型[25]是以NDVI数据、气象数据(气温、降水、太阳辐射)以及植被类型数据为输入信息来驱动模型。在模型中，NPP的值由植物吸收的光合有效辐射(APAR)和光利用率(ε)2个因子来表示，其估算公式如下：

NPP(x,t)=APAR(x,t)×ε(x,t)

(1)

式中:t为时间；x为空间位置；APAR(x,t)为像元x在t月份吸收的光合有效辐射[MJ/(m2·月)]；ε(x,t)为像元x在t月份的实际光利用率(g C/MJ)。具体解算过程参见文献[25]。

在NPP估算过程中，会由于一些不可控因素引入误差，因此需要对估计结果进行精度评价。本文采取与其他学者的成果进行对比的方法，从而验证CASA模型在本研究区的适用性。发现模拟结果与其他学者的总体相似(表1)，与李萌萌[26]基于改进的CASA模型估算的陕西省植被NPP结果相比较，各类植被的NPP值基本一致，其中差别较大的是落叶阔叶林，可能的原因是研究区内阔叶林种类不同所致；其与CEVSA模型[27]和GLO-PEM[28]模型估算的结果存在一定的差异，但其模拟结果均表现为林地NPP偏高，草地、农田偏低，误差的存在是必然的，但是其值的大小均在可以接受的合理变动范围之内，说明本研究利用CASA模型估算的NPP是具有较高准确性的，并且具有一定价值，可作为区域NPP研究的基本数据。

表1 本文植被NPP模拟结果与其他模型研究对比

2.2 Theil-Sen Median趋势分析

Theil-Sen Median趋势分析[29]是计算样本数据的中值，计算公式为：

(2)

式中:SNPP为Theil-Sen Median趋势；NPPi，NPPj分别为第i年和第j年像元的NPP值。该计算方法计算的是n(n-1)/2个数据组合斜率的中位数。当SNPP>0时；NPP呈上升趋势,反之则呈下降趋势。

2.3 偏相关方法

本研究的NPP是基于气象因子作为输入数据利用CASA模型模拟得到的，其与降水、气温具有一定的关系，因此通过相关分析方法[30]来分析NPP与降水量或气温的相关性，相关性系数公式如下：

(3)

在分析多个因素对NPP的影响时，暂不考虑其他因素的影响而单独研究某一个因素与NPP的相关程度即为NPP的偏相关分析，偏相关系数计算公式如下：

(4)

式中:rxyz为变量z固定后变量x与变量y的偏相关系数;rxy，rxz，ryz分别为变量x与变量y、变量x与变量z、变量y与变量z的相关系数。

本文采用T检验方法对以上的偏相关系数进行显著性检验，计算公式如下：

(5)

2.4 基于Thornthwaite模型估算植被潜在生产力

本研究选用Thornthwaite纪念模型[1]来估算潜在NPP，模型计算公式如下：

NPP=3000×[1-e-0.0009695(E-20)]

(6)

式中:NPP为潜在净初级生产力[g C/(m2·a)]；E为年实际蒸散量(mm)。

另外，E的具体计算式如下：

(7)

L=2000+25T+0.05T3

(8)

式中:r为年总降水量(mm)；L为年最大蒸散量(mm)；T为年平均气温(℃)。

2.5 气候变化和人类活动对植被NPP变化相对影响的评估方法

本文参考郭继凯[31]、周妍妍[20]等有关NPP变化的情景设定方案，通过每年不同类型NPP的变化量来评估气候变化和人类活动对植被生产力的相对贡献。其中，第1种是实际NPP的变化量(NA)；第2种是潜在NPP的变化量(Np)；第3种为人类活动对NPP的影响(NH)，通过计算潜在NPP变化量跟实际NPP变化量的差值来获得(NH=NP-NA)。

通过计算NA,NP,NH的斜率KA，KP，KH来区分气候变化和人类活动对NPP的改善和退化作用。KA>0代表实际NPP趋于增加状态，KA<0表示实际NPP逐渐减少。KP>0代表气候变化有利于植被恢复，KP<0表明气候变化导致植被退化。KH>0代表人类活动导致植被退化；KH<0表示人类活动促进植被生长。

计算NA,NP和NH的变化趋势，采用的方法是Theil-Sen Median方法。具体情景设定方案见表2。

表2 各种可能情景下气候变化与人类活动在植被净初级生产力变化中相对作用的评价方法

3 结果与分析

3.1 北方农牧交错带植被NPP空间分布及变化趋势

3.1.1 NPP的空间分布 本文以CASA模型为基础，制作了1998—2016年北方农牧交错带的植被NPP数据，空间分辨率为1 km。植被NPP多年平均值空间分布见图1，空间上呈西南、东北高，中间低的分布特征，其值为13.87～1 251.06 g C/(m2·a)，主要集中分布在200～500 g C/(m2·a)，总体平均值为364.93 g C/(m2·a)。植被多年平均NPP的高值区主要位于西南部的祁连山山麓地带、中部的燕山山脉、吕梁山、恒山等地、东部的大兴安岭地区，而低值区主要分布在中西部植被稀少的干旱地区，植被NPP多年平均值为13.87～200，200～500，500～1 251.06 g C/(m2·a)时，分别占总面积的11.58%，74.43%，13.99%。

3.1.2 NPP的时空变化特征 1998—2016年北方农牧交错带植被NPP年均值变化见图2，近19 a间，植被NPP年均值变化范围在302.70～459.76 g C/(m2·a)，呈显著的上升趋势，2013年的植被NPP年均值为最高，达到了459.76 g C/(m2·a)，2001年的植被NPP年均值最低，每年的植被NPP的平均增加量约为3.95 g C/(m2·a)。

图1 1998-2016年北方农牧交错带植被NPP多年平均值空间分布

图2 1998-2016年北方农牧交错带NPP年均值变化

采用Theil-Median趋势分析方法模拟了每个像元的净初级生产力变化趋势，反映了19 a间研究区植被NPP变化的空间特征(图3)。研究区内植被NPP增加的区域明显大于植被NPP减少的区域，说明植被生长状况较好。显著增加的区域主要分布在研究区中、西部的偏南地区及东北平原的部分零散区域；显著减少的区域较少，主要分散分布在内蒙古自治区东北部。其中，植被NPP增加速率大于3 g C/(m2·a)的区域面积占总面积的56.95%；植被NPP减小速率大于3 g C/(m2·a)的区域面积占总面积的5.42%；植被NPP增加速率和减小速率小于3 g C/(m2·a)的区域面积分别占总面积的13.18%，24.55%。

3.2 气候因子的空间分布及变化趋势

3.2.1 气候因子的时空变化特征 根据降水和气温的线性拟合结果可以看出(图4)，表明1998—2016年北方农牧交错带的降水量呈上升趋势，速率为3.945 4 mm/a(图4A)。年平均气温呈上升趋势，速率较小，线性回归增加速率是0.007 3℃/a(图4B)。近19 a间，2000年的总降水量最低，为316.29 mm，2012年降水量最高，年平均气温最低，分别为509.31 mm和4.74℃，2015年的年平均气温最高，为6.94℃。

图3 1998-2016年北方农牧交错带NPP的Theil-Sen Median趋势分析结果

根据北方农牧交错带及其周边166个气象站降水资料的插值数据进行趋势分析，1998—2016年的年降水量空间分布如图5A所示，可以看出研究区内的降水量具有明显的空间差异，总体上东南高，西北低。研究区的最高年均降水量为715.01 mm，最低年均降水量为194.74 mm。受海陆位置和地形等地理环境的影响，研究区内年均降水量具有明显的区域性差异，自东南向西北逐渐减小，大兴安岭、阴山山脉、黄土高原北缘一线的西北地区年均降水量多在200～400 mm，东南季风区距离海洋较近，雨水充足，降水量在400～600 mm。由研究区的年平均气温空间分布图(图5B)可以看出，该地区的气温也有明显的空间差异，呈现南高北低，东高西低的趋势。全年年均气温小于5℃的区域主要位于阴山山脉以北和大兴安岭以西以及研究区西南部青海省和甘肃省的祁连山地区，其中部分地区小于0℃。年平均气温在10℃以上主要分布在无定河上游的小部分区域，其余地区年平均气温为5～10℃。

图4 1998-2016年北方农牧交错带降水气温变化

图5 1998-2016年北方农牧交错带年降水量和年平均气温空间分布

3.2.2 气候因子的空间变化特征 近19 a间，研究区大部分区域年降水量呈增加趋势(图6A)，呈增加、减少趋势的面积分别占总面积的96.72%，3.28%，增加较为明显的区域分布于阴山山脉以南的陕西省东北部和山西省西北部，西辽河平原吉林省铁岭市北部、松嫩平原等地。年降水量呈下降趋势的区域主要位于阴山山脉以东的内蒙古自治区中东部和六盘山以西的甘肃省东南部等地。在全球变暖的大趋势下，1998—2016年研究区内70.60%的区域年平均气温呈上升趋势(图6B)。增温明显的区域主要分布在大兴安岭的东北部、燕山北部、吕梁山北部、黄土高原西南部、祁连山东北部山麓地带。年平均气温呈下降趋势的区域占研究区总面积的29.40%，主要位于松嫩平原、浑善达克沙地东南部、鄂尔多斯高原等地。

3.3 NPP变化与影响因子的关系

3.3.1 NPP变化与气候变化的关系 从植被NPP与降水的相关分布及相关显著性分布(图7A、图8A)可以看出，研究区植被NPP与降水量之间主要呈正相关，植被NPP与年降水量的总体平均偏相关系数为0.369 2，研究区内97.70%的像元NPP值和降水量正相关。NPP和降水量的正相关的区域主要分布在黄土高原、鄂尔多斯高原、河套平原、大青山、浑善达克沙地、大兴安岭、松嫩平原、西辽河平原等大部分地区，降水的增加可促进植被生长状况的改善，但这些区域的正相关关系并不显著；而NPP和降水呈负相关关系极显著的区域主要位于研究区西南部祁连山东北部山麓、吕梁山和太行山北部、燕山山脉周边等地，这些区域海拔较高或者纬度较高，气温较低，说明降水的增多并不能改善该地的植被生长。图7B和图8B为研究区植被NPP与相应时期平均气温的相关关系及相关显著性分布，NPP与平均气温的整体平均偏相关系数为0.004 8，呈负相关的像元占比为50.49%，且负相关关系是极显著的；NPP与平均气温呈正相关关系不显著的区域较多，集中分布在青海省东北部、甘肃省东南部及松嫩平原南部等地区，研究区内中部和东北部的大部分区域海拔较高、纬度较高，冬季长，热量累积较少，从而表明热量不足会抑制植被生长。

图6 1998-2016年北方农牧交错带年降水量和年平均气温变化趋势空间分布

图7 1998-2016年北方农牧交错带植被NPP与气象因子的空间相关性分布

图8 1998-2016年北方农牧交错带植被NPP与气象因子的相关显著性分布

3.3.2 气候变化和人类活动对NPP的相对影响 从气候变化和人类活动影响植被NPP变化的空间分布可以看出，气候变化和人类活动对植被影响的积极作用明显大于其对植被的消极贡献(表3，图9)。气候变化主导植被恢复的面积占总面积的47.64%，主要分布在研究区东部的西辽河平原和松嫩平原，主要土地利用类型为农田，气候变化主导植被恢复的情况在研究区中西部的农牧交错带呈零星分布。人类活动使植被NPP增加的面积占总面积的2.16%，集中分布在甘肃省兰州市周边、内蒙古自治区乌兰察布市和呼和浩特市边缘交界处。二者共同作用促进植被改善的面积占总面积的30.56%，主要分布研究区中西部的农牧交错地区。

表3 各影响因素主导植被NPP变化的比例分布

图9 气候变化和人类活动主导植被NPP变化的空间分布

气候变化导致植被NPP退化的面积占总面积的0.56%，主要分散分布在甘肃省白银市、内蒙古自治区乌兰察布市等地。人类活动主导植被NPP减少的面积占总面积的17.67%，集中分布在内蒙古自治区中、东部等地，主要土地利用类型为草地和林地。二者共同作用导致植被退化的面积占总面积的1.41%，主要离散分布在甘肃省白银市西部、内蒙古自治区呼和浩特市东部、乌兰察布市西部等地。由上述分析可知，在研究区内的甘肃省兰州市和白银市交界处、内蒙古自治区呼和浩特市和乌兰察布市交界处，引起植被变化的影响因素交叉分布。

4 讨论与结论

4.1 讨 论

(1) 在分析NPP与影响因素的关系时，仅考虑了降水、气温和人类活动的影响；由于研究区范围较大，在分析自然因素影响时，没有考虑到地形地貌、土壤质地等方面，对人类活动主控因素的确定没有具体化。需要进一步对研究区的NPP变化成因进行分段分析，并结合多方面的影响因子来说明变化成因。

(2) 限于研究区的范围较广、研究尺度较大等原因，没有进行实地植被NPP的野外调查，缺少实地验证，在以后的工作中，应加强设计合理的野外观测试验，从而进一步提高模型的模拟精度。

4.2 结 论

(1) 1998—2016年，研究区的NPP年均值空间上呈现西南、东北高，中间低的特点，研究区植被的NPP年均值主要集中分布在200～500 g C/(m2·a)。植被NPP增加的区域明显大于植被NPP减少的区域，显著增加的区域主要分布在研究区中、西部的偏南地区及东北平原的零星地区；显著减少的区域分散分布于研究区的东北边缘。

(2) 研究区内气象要素(降水和气温)的分布具有明显的时空差异，年降水量呈东南高、西北低分布，而年平均气温分布呈现南高北低、东高西低，在研究时段内，年降水量整体呈增加趋势，大部分地区的年平均气温也呈上升趋势。

(3) 北方农牧交错带的植被NPP与降水量之间存在正相关关系，而与年平均气温具有明显的负相关关系。研究时段内，气候变化和人类活动对植被影响的积极贡献明显大于对植被变化的消极作用，气候变化、人类活动以及气候变化和人类活动共同作用造成植被恢复的面积分别占恢复总面积的59.28%，2.68%，38.04%；其三者造成植被退化的面积，分别占退化总面积的2.87%，89.96%，7.17%。