植被净初级生产力(net primary production, NPP)指绿色植物在单位面积单位时间内所累积的有机物数量，不仅是陆地生态系统中物质与能量转换的基础，而且能反映植物群落的生产力。此外，植被NPP还可用于评估生态系统的功能协调性及其与其他环境因素的相互作用^[1]。因此，在全球气候变化的大背景下，植被NPP对气候变化的响应机制方面发挥着重要作用。随着“退耕还草”等一系列生态保护政策和生态恢复工程的推行，长时序、定量化的植被NPP研究分析对于评价生态系统质量状况，监测区域植被变化具有重要意义。

随着国内外学者对NPP研究的逐渐深入，NPP获取方法也在不断革新。实地测量法^[2-3] 最先被应用，但由于工作量大，获取周期长等诸多因素不能被推广。随后利用模型模拟间接估计的多种方法被提出，其中属于遥感反演法之一的Carnegie-Ames-Stanford Approach (CASA)模型具有模型简单、数据获取方便、精度较高等优势，被广泛应用^[4-5]。朴世龙等^[6]采用CASA模型估算了中国1997年植被NPP分布情况；张美玲等^[7]将顺序分类系统和温湿度指标引入CASA模型，得到2004－2008年中国草地NPP时空分布，并分析得到影响NPP的决定因素。此后，对于NPP的研究^[8-13]更倾向于基于典型区域、多源数据的NPP时空分布特征和驱动力响应机制研究，旨在揭示NPP和独特区域特征间(高程、坡向、气候因子等)的内在联系。近年来我国针对类似甘南藏族自治州的脆弱生态区推行了“退耕还草”等一系列生态保护政策和生态恢复工程，这一系列举措对区域植被变化造成了一定的改变，引起了学者们的兴趣。学者们开始对甘南藏族自治后植被NPP进行研究，采用趋势分析法和变异系数法探究了NPP时空变化规律与气温和降水的响应关系，并分析了区域内NPP在不同植被类型下的变化趋势分布^[14-16]。但针对甘南藏族自治州植被NPP的已有研究中，长时序下该地区植被NPP重心有无迁移缺乏探讨，此外，不同植被类型下探究甘南藏族自治州NPP与归一化植被指数(normalized difference vegetation index, NDVI)、太阳辐射、气温、降水各因子之间的响应机制及趋势分布未见报道。

本研究基于MODIS遥感数据、气象站点数据、植被类型数据，利用CASA模型计算了甘南藏族自治州2000－2019年的植被NPP；利用相关系数法探究NPP与NDVI、气温、降水、太阳辐射间的响应关系及不同区间下的变化和分布。通过一元线性回归趋势分析NPP的变化趋势，基于重心模型^[17-18]揭示长时序下植被NPP在空间上的差异和重心迁移规律，并进一步探究NPP在不同植被类型下与太阳辐射、气温、降水的变化趋势，明确甘南藏族自治州NPP与各因子间的相应机制。本研究旨在为认识区域生态环境、生态措施制定与实施实效提供理论支撑。

1.   研究区概况

甘南藏族自治州(100°46′～104°44′ E，33°06′～36°10′ N)地处甘肃省西南部，青藏高原和黄土高原过渡地带，全州下辖夏河、玛曲、碌曲、卓尼、迭部、临潭、舟曲和合作共7县1市；东部山地交错，农牧业兼营；西北部为广阔草原，是重要的农牧业生产基地和国家级生态主体功能区(图1)。甘南藏族自治州年平均气温在1～13 ℃，年内高于10 ℃的月份不足3个月，昼夜温差极其明显，年均降水在440～800 mm间且集中在7月－9月，其气候为典型的高原大陆性季风气候^[19-20]。甘南藏族自治州独特的气候特征和境内丰富的植被资源为植被NPP的研究提供了自然基础。

图  1  研究区行政区划图及植被类型分类图

Figure  1.  Administrative division map and vegetation cover classification map of the study area

下载: 全尺寸图片 幻灯片

2.   材料与方法

2.1   数据来源及处理

本研究所用气温和降水等气象数据均来源于国家气象科学数据中心中国地面气候资料月值数据集(http://data.cma.cn/data/)。由于甘南藏族自治州气象数据有限，获取时空数据插值误差较大，为提高精度选取西北地区110个气象站点的月平均气温和降水量数据插值得到西北地区气候时空数据，最后裁剪出研究区范围数据，站点分布如图2所示。数据预处理过程中，首先对缺测值和异常值进行预处理，将缺测值超过6个月的数据予以剔除，未超过6个月的通过MATLAB R2018b进行线性插值后得到，将降水量低于0 mm的数据用前后月份的均值作为替换。之后利用 ANUSPLIN 软件对研究区域的气温、降水数据在西北地区尺度上插值处理后裁剪出研究区2000－2019年逐月气温和降水分布图，进行波段组合后得到年平均气温、降水分布图，通过ANUSPLIN方法插值时引入了数字高程模型(digital elevation model, DEM)数据，以提升局部和全局的插值精度。

图  2  气象站点分布图

Figure  2.  Distribution of meteorological stations

下载: 全尺寸图片 幻灯片

太阳辐射数据来自于国家气象科学数据中心(http://data.cma.cn/)的中国辐射气候标准值月值数据集。经筛选和预处理后共选取118个气象站点数据，处理方法与生成气温、降水图一致。

NDVI数据来自于美国地质调查局土地过程分布活动档案中心(https://lpdaac.usgs.gov//) 发布的MOD13Q1全球产品，数据周期为16 d，空间分辨率为500 m。选取覆盖研究区甘南州的影像，利用MODIS Reprojection Tool (MRT)软件进行影像的批量拼接、转投影、最大值合成处理。将NDVI 在0.10以上的区域定义为植被覆盖区，去除NDVI小于0.10的无植被覆盖区。借助ArcGIS 10.3软件提取甘南2000－2019年各月份NDVI数据，并统一定义为 Albers 投影。最后与气象数据处理方法相同，将全年12个月的数据进行波段组合后得到年均 NDVI分布图。

数字高程模型(DEM)数据来源于NASA (National Aeronautics and Space Administration)的SRTM数据(http://srtm.csi.cgiar.org//)，空间分辨率为500 m。利用嵌套在ANUSPLIN包中的薄板平滑样条进行空间插值，其中薄板平滑样条包含了温度和降水数据的高度依赖性，协变量的引入能够使得插值精度更高，插值结果更为平滑。

用来与CASA模型反演NPP结果进行一致性分析的MOD17A3 NPP产品数据来源于MODIS MOD17A3H版本产品中年度净初级生产力数据集，时间分辨率为1年，空间分辨率为500 m，且已经通过大气校正、辐射校正、几何校正和去云等处理。数据下载完成后使用ArcGIS 10.3软件进行拼接、裁剪、重投影后剔除其异常值，得到2000－2020年的NPP精度验证数据集。

植被类型数据来自中国科学院资源环境科学与数据中心(https://www.resdc.cn/)提供的共分为22种植被类型的栅格数据，掩膜提取出甘南藏族自治州植被类型数据，经整合后重分类为研究区主要的5种植被类型(林地、灌丛、草地、湿地、耕地)，并定义为Albers投影，用于CASA模型中最大光能利用率的确定。

2.2   空间插值方法

本研究中站点数据采用国际上认可并被广泛应用的ANUSPLIN软件进行插值，由于引入数字高程模型数据，其插值结果精度更高，结果的细节展示也更为突出。本研究中主要使用薄盘光滑样条函数和LAPGRD (局部薄盘光滑样条表面估值并计算贝叶斯标准误差)等主要程序模块。ANUSPLIN软件的核心是局部薄盘光滑样条法，其理论模型为：

This page contains the following errors:

Below is a rendering of the page up to the first error.

2.3   CASA模型估算NPP方法

选取改进后的CASA模型^[5]来估算NPP。该模型以NDVI为驱动，借助气温、降水、太阳辐射为主要驱动因子，在对相关参数(ε、SR、SR_min、SR_max)进行了适宜性调整后对NPP进行反演。模型表达如下：

This page contains the following errors:

Below is a rendering of the page up to the first error.

式中：S_(x，t)为像元x在t月份的太阳辐射量(MJ·m⁻²)；FPAR_(x，t)为植被层对光合有效辐射的吸收比例；常数0.5表示植被所能利用的有效辐射占总辐射的比例。其中FPAR、NDVI、比值植被指数(SR)之间的关系公式如下：

式中：SR_min和SR_max分别对应各植被NDVI的5%和95%下侧百分位数，SR_{(x, t)}可由NDVI_{(x, t)}求得，公式如下：

式中：NDVI_(x，t)为像元x在t月份的NDVI。

2.4   相关性分析

相关性分析主要用来反映各个因子之间相关程度和相关方向，本研究采用皮尔森相关系数法探讨Z-score标准化后的NPP与各因子的响应关系，公式如下：

式中：x_i、y_i为其他因子(x)和NPP (y)第i年的值。R_xy为变量x和y年际值的相关性系数，取值范围为−1.0到1.0。采用F显著性检验结果，将相关性分为为8个等级：显著负相关(r < −0.8)、中度负相关(−0.8 < r < −0.5)、轻度负相关(−0.5 < r < −0.3)、弱负相关(−0.3 < r < 0)、弱正相关(0 < r < 0.3)、轻度正相关(0.3 < r < 0.5)、中度正相关(0.5 < r < 0.8)、显著正相关(0.8 < r)。

2.5   趋势分析

一元线性回归分析通过拟合一条数据趋势线来模拟每一个像元的年际变化趋势，其优点是可消除特定年份由于异常因素对NPP变化趋势的影响，真实反映NPP的演化趋势。计算公式：

式中：θ_slope为拟合得到的年际NPP数值变化速率，负值表示NPP呈下降趋势，正值则表示NPP呈增长趋势。n为所研究时间序列长度，NPP_i为第i年的NPP值。根据Slope和显著性检验结果(P值)将NPP空间变化趋势分为5个等级：显著减少(Slope < 0, P < 0.05)，不显著减少(Slope < 0, P ≥ 0.05)，不显著增加(Slope > 0, P ≥ 0.05)，显著增加(Slope > 0, 0.01 ≤P < 0.05)和极显著增加(Slope > 0, P < 0.01)。

2.6   Z-score标准化

在研究不同量纲的因子间变化趋势分布时，将它们的趋势值分别做Z-score标准化处理用于后续的研究分析，表达式如下：

式中：X为要进行Z-score标准化的一元线性回归后的趋势像元值，其中X(i)为第i年的趋势像元值，Mean (X)和Std(X)是研究期间内X的平均值和标准偏差。Z-score标准化在多个目标值上界和下界已知的情况下，将有量纲表征的数据规约变换成相同尺度内无量纲的数据，消除样本奇异值。

2.7   NPP重心分析

“重心”的迁移方向和距离可以反映某一地理现象的程度和趋势的变化。重心模型被广泛应用于研究经济社会问题，包括人口分布及其迁移、区域经济增长趋势、粮食分布及其变化、土地利用分布及其变化、生态环境变化等。该模型对事物内部变化趋势有着较强的展示，本研究利用重心模型分析植被NPP增速在空间上的差异和植被NPP重心迁移规律。

重心定义为一对坐标(x, y)，模型内笛卡尔坐标设为(x_i , y_i) 。然后将n个空间单元(像元)组成的区域的空间均值定义为笛卡尔坐标。结合植被NPP像元特点，应用于研究区NPP重心模型迁移时的公式为：

式中：x_i和y_i是要素i的坐标，n等于要素总数。NPP_i为第i个像元处NPP的值。

3.   结果与分析

3.1   甘南藏族自治州植被NPP时空变化特征

甘南藏族自治州2000－2019年植被NPP总体呈稳定增长趋势且年均植被NPP具有明显的时间分异性(图3)。植被NPP年均增长率约为7.5 g·(m²·a)⁻¹，变化范围为556～691 g·(m²·a)⁻¹，其中表现为增长趋势的时间段为2000－2002年、2003－2008年、2011－2015年、2018－2019年，且增长率大致相同，约为15 g·(m²·a)⁻¹。表现为减少趋势的时间段为2002－2003年、2008－2011年、2015－2018年，其中2002－2003年和2008－2011年下降趋势大致相同，下降率为12.3 g·(m²·a)⁻¹，2015－2018年下降较为缓慢，下降率约为6.7 g·(m²·a)⁻¹ (此处增长率和下降率是由CASA模型反演得出的2019年NPP与2000年NPP的差值与时间跨度得出，并非是拟合后的趋势)。总体上，甘南藏族自治州植被NPP在2000－2019年间呈现出稳定增长趋势，虽然时有下降趋势但并不明显。

图  3  2000－2019年甘南藏族自治州CASA模型模拟年际NPP值及MODIS NPP产品年际值

Figure  3.  Interannual net primary production (NPP) value simulated by CASA model and interannual value of MODIS NPP products in Gannan Tibetan Autonomous Prefecture from 2000 to 2019

下载: 全尺寸图片 幻灯片

This page contains the following errors:

Below is a rendering of the page up to the first error.

图  4  甘南藏族自治州2000－2019年NPP分布图

Figure  4.  Net primary production (NPP) distribution map of Gannan Tibetan Autonomous Prefecture from 2000 to 2019

下载: 全尺寸图片 幻灯片

3.2   甘南藏族自治州NPP重心迁移

年际NPP空间分布差异和增速的空间异质性可通过植被NPP重心的迁移有效呈现。在CASA模型反演植被NPP结果的基础上，通过计算2000－2019年年均NPP和逐年NPP的重心，并考虑到极端气候等因素的影响以10年为时间尺度分为2000－2009年、2010－2019年2个时间段来研究甘南藏族自治州植被NPP的重心迁移规律(图5)。重心的迁移不仅反映出植被NPP聚集的转移，亦是植被NPP增速空间异质性的一种表现。如图5a所示，整体重心迁移呈“自西向东，从北到南”跳动偏移的趋势。2000－2002年、2005－2006年间NPP重心整体出现向南迁移的趋势，说明该时间段内南部植被NPP增速大于北部；2003－2009年NPP重心整体呈现从西南向东北迁移的趋势，且西向迁移的趋势大于北向迁移的趋势，说明在此期间东部地区植被NPP增速远超西部地区。2010－2019年甘南藏族自治州NPP重心迁移幅度较2000－2009年相比更大，跳跃也更明显，NPP重心更为明显的跳跃说明研究区内植被NPP的变化较2000－2009年更为剧烈。其中2010－2011年、2012－2013年、2014－2015年、2016－2017年和2018－2019年这些时段内植被NPP出现由西北向东南方向迁移的趋势，说明在此时间段内研究区东南地区植被NPP增速大于东北地区(图5b)。2011－2012年、2013－2018年植被NPP重心出现明显的由东南向西北方向迁移的趋势。

图  5  甘南藏族自治州2000－2019年植被NPP重心迁移图

Figure  5.  Center of gravity migration of vegetation net primary production (NPP) in Gannan Tibetan Autonomous Prefecture from 2000 to 2019

下载: 全尺寸图片 幻灯片

3.3   甘南藏族自治州植被NPP对各因子的响应

This page contains the following errors:

Below is a rendering of the page up to the first error.

图  6  植被NPP与各因子相关性系数空间分布

Figure  6.  Spatial distribution of correlation coefficient between vegetation net primary production (NPP) and various factors

下载: 全尺寸图片 幻灯片

This page contains the following errors:

Below is a rendering of the page up to the first error.

3.4   不同植被类型下NPP对气候因子的响应

本研究中将研究区划分为林地、灌丛、草地、湿地、耕地共5种植被类型(图1)，分析2000－2019年间NPP与各驱动因子间的变化特征，并与其他区域NPP对气候因子响应研究对比。NPP和主要气候因子在不同植被类型下的分布存在明显差异(图7)。其中降水和辐射由于研究区相同的原因，在5种植被类型下的大体趋势一致且不同类型下差异不大。与降水和辐射相比，气温在5种类型下的差异尤为明显，结合研究区植被分布特点和相关文献结论，气温的明显差异主要是由于不同植被的蒸腾与呼吸作用的差异导致^[7]。相比气温、降水、辐射这3种气候因子在不同类型下的不同趋势，NPP在不同类型下的差异相对较小并稳定维持在575～650 g·m⁻²，符合NPP与气候因子内在相关性，即多种气候因子共同影响NPP的变化，NPP是响应各个因子变化后的结果^[10]。

图  7  不同植被类型下NPP与各因子分布

Figure  7.  Distribution of net primary production (NPP) and factors under different vegetation types

下载: 全尺寸图片 幻灯片

在对NPP和各个因子做长时序趋势分析后，对其趋势分析结果(各因子年际增量)分别进行Z-score标准化处理后发现：研究区林地NPP高度集中在0.4～0.6趋势区间内(84.62%)，NPP、太阳辐射、气温、降水趋势分布集中在0.2～0.6趋势区间内；灌丛NPP与降水分布高度集中在0.4～0.6趋势区间内(分别为86.62%、61.59%)，NPP、太阳辐射、气温、降水多分布在0.2～0.6趋势区间内；草地NPP的分布仍高度集中在0.4～0.6趋势区间内，其余因子的分布集中于0.4～0.8趋势区间内，降水亦偏向0.4～0.6趋势区间(61.92%)；湿地的各因子分布与林地、灌丛、草地分布差异明显，NPP高度集中在0.4～0.6趋势区间内(87.56%)，太阳辐射集中分布于0～0.2、0.2～0.4区间内(分别占比66.98%、20.09%)，气温、降水分别高度集中在0.6～0.8、0～0.2趋势区间内(分别占比91.83%、98.87%) (图8)。在2000－2019年间，林地、灌丛、草地的NPP的年际变化趋势大致相同呈稳定增长趋势，集中在0.4～0.6趋势区间，其他因子的趋势分布差异并不明显且在各个趋势区间上分布相对均匀。湿地的各因子趋势变化差异较大，其中气温在研究时段增长相对明显，降水和辐射主要集中在低增长率区间内，NPP呈现平稳增长趋势。耕地的降水和气温变化不明显，NPP呈现平稳增长趋势，辐射呈现相对明显的增长趋势。

图  8  不同植被类型下NPP和各因子的趋势分布

Figure  8.  Trend distribution of net primary production (NPP) and factors under different land use types

下载: 全尺寸图片 幻灯片

4.   讨论

4.1   CASA模型反演NPP与MOD17A3H系列的NPP数据的一致性分析

本研究将CASA模型反演的NPP数据集与MOD17A3H系列的NPP数据集对比分析以验证CASA模型反演的NPP数据集的精度。结果表明，两者的年际NPP值差异不大，在1.2～72.6 g·m⁻²间波动，整体变化趋势都表现为稳定上升趋势，线性拟合后的起始最值也差异较小保持在50.0 g·m⁻²内(图2)。出现这些差异主要是由于获取机理的不同导致。MOD17A3H系列NPP数据集侧重于地表NPP的反馈，受地物类型影响大，而作为NPP光能利用率反演模型的CASA模型是以遥感数据NDVI为驱动因子，结合主要驱动因子(气温、降水、太阳辐射等)，在遥感数据的基础上对NPP值进行反演。结合以上两种NPP数据源获取机理的不同^[10]，本研究认为出现的数值上的和整体变化趋势的较小差异是正常的，也从侧面印证所用CASA模型反演的NPP数据集的可靠性。

4.2   甘南藏族自治州植被NPP时空变化及重心迁移

基于CASA模型估算甘南藏族自治州植被NPP，发现在2000－2019年间甘南藏族自治州植被NPP整体呈稳定增长趋势，表现出明显的时间分布上的差异。其中2000－2002年、2003－2008年、2011－2015年、2018－2019年呈增长趋势，2002－2003年、2008－2011年、2015－2018年呈现出下降趋势。相关政策资料显示：2003年以来我国开始实行退耕还林还草政策，同年甘肃、四川等8个省(区)开始推行退牧还草计划^[21]；甘肃省在2005年继续推进退牧还草工程实施；2011－2015年间，中央每年在内蒙、宁夏、甘肃等8个省(区)拨款以推进草原生态保护补助奖励政策。对比本研究分析结果，这些政策实施的时间节点2003年、2005年、2011－2015年内甘南藏族自治州NPP都处于明显的增长趋势。研究认为表现为增长和下降趋势的年份分布的差异性主要是由于降水、气温和太阳辐射的时滞效应造成^[22]。这与叶许春等^[13]对气候因子时滞性效应的研究结论是一致的，不同降水量、气温、太阳辐射对年际NPP变化均会产生滞后性的表现。适宜的降水量、气温和太阳辐射对随后时间内NPP的累积有决定性作用。

This page contains the following errors:

Below is a rendering of the page up to the first error.

重心理论在众多领域被广泛应用^[17-18]，本研究通过构建NPP重心，将研究区内2000－2019年年均NPP重心和逐年NPP重心以10年尺度分段，探究甘南藏族自治州近20年植被NPP重心迁移规律和甘南藏族自治州NPP的空间分异性。甘南藏族自治州NPP重心迁移方向及幅度的不同说明了在此期间地区间植被长势的差异，可以得出西北部植被长势不如东南部，或因近20年内气候变化、人类活动、政策实效不同所导致。

4.3   甘南藏族自治州植被NPP与各因子间的响应关系

近年来众多专家学者对NPP与各影响因子(NDVI、辐射、气温、降水)的相关性和潜在的响应机制研究分析^[23-24]。本研究发现，2000－2019年间甘南藏族自治州植被NPP与NDVI、辐射、气温、降水间均存在明显的相关关系，且相关关系在空间上有明显差异。结果显示绝大部分地区上植被NPP与四者之间显著正相关，而Zhang ^[23]研究显示，NPP与降水负相关，这种差异可能是由于地域和植被类型不同导致。仅有少数地区呈负相关关系且并不显著，呈负相关关系的区域有着植被稀疏、高海拔的共同特征。光热条件和水分状况对NPP的影响占据主导地位，高海拔对植物产生的胁迫作用或是造成负相关的主要原因^[25]。在本研究中，NPP与降水和太阳辐射显示出较高的相关性，对研究区特征分析后得出区域高海拔特点和独特植被覆盖特质是导致差异的主要因素。区域特征往往是影响研究结果的决定因素^[26]，目前以生态指标为基础的气候模型很难具有区域普适性，在之后的研究中会聚焦此问题。

研究期间人口增长导致农田和建设用地增加，林地减少。随着人口和粮食消费的增加，土地资源变得有限，通过化肥、机械和农用化学品等生产要素的投资，土地利用资源利用效率得到提高，单位面积内植被NPP 也得到了提高，在此背景下分别研究不同植被类型下NPP和各因子的变化和响应十分必要^[27]。本研究将研究区划分为林地、灌丛、草原、湿地、耕地5种类型后，分别探究NPP和各驱动因子在不同利用类型下的变化和响应机制后发现：辐射在不同植被类型下变化趋势大致相同，气温在不同植被类型下差异明显，或与不同植被的蒸腾作用和呼吸作用的不同有关；NPP在不同植被类型下的差异相对较小，稳定维持在575～650 g·m⁻²。以上结论与贾俊鹤^[28]等得出的结论基本一致，并符合NPP与各气候因子内在相关性，且与多种气候因子共同影响着NPP的变化，NPP是响应各个因子变化后的结果。在不同植被类型下对NPP及各个因子做归一化趋势分析后发现：在林地、灌丛、草原植被类型上，NPP高度集中在中度增长趋势的区间(0.4～0.6)，其他因子趋势分布差异并不明显，在各个趋势区间上分布相对均匀而湿地区域内各因子趋势分布差异较大。与相关研究结论^[18]不同的是，本研究中降雨在不同土地类型下变化分布差异相对明显，可能是植被类型不同导致^[21]。对于趋势分布差异较大的区域因地制宜地丰富植被种类，实行适宜的生态恢复计划是植被活动增加的主要驱动力^[17]。

5.   结论

本研究基于MODIS遥感数据、气象站点数据、植被类型数据等，采用CASA模型反演2000－2019年甘南藏族自治州植被NPP，反演结果得到MOD13A系列NPP数据集的验证。获得如下主要结果：

This page contains the following errors:

Below is a rendering of the page up to the first error.

2) NPP重心迁移呈由东向西、从北到南跳动偏移的趋势，研究区东部植被NPP增速远超西部；2010－2019年NPP重心迁移幅度较2000－2009年相比更大，跳跃更为明显。

3)研究区内植被NPP与各因子大体呈现明显的正相关关系，仅有极少区域呈负相关关系；辐射和气温对NPP的影响占据主导地位。林地、灌丛、草地NPP与各因子均呈稳定增长趋势且趋势分布差异不明显，湿地类型下气温处于明显增长趋势区间，降水和辐射处于低增长率区间，NPP呈平稳增长趋势。整体表明甘南藏族自治州所实施的生态保护政策效果良好。