植被净初级生产力(Net Primary Productivity, NPP)是反映自然条件下植被生产力的重要因子^[1], 是指单位面积和单位时间内绿色植物积累有机物的数量^[2]。NPP不但反映了自然条件下植被固定CO₂的能力, 并且表征了生态系统的质量状况和生产力, 也是评估陆地生态系统可持续发展的重要要素^[3]。国际地圈生物圈计划(IGBP)、全球变化与陆地生态系统(GCTE)以及京都协议书等都将植被NPP的研究作为核心内容之一^[4], 分析植被NPP时空格局及其对气候和人类活动的响应, 对陆地生态系统保护与可持续发展具有重要意义。

光能利用率模型(CASA模型)是由遥感数据驱动的一种NPP估算模型, 由于其能较高精度估算不同尺度不同植被类型的NPP, 因此被广泛应用^[5]。1993年Potter等^[6]估算了全球植被NPP, 并提出了CASA模型。Yuchao Yan等^[7]定量分析了2000—2015年影响中国北方草地NPP驱动因素的贡献度。GUO Bing等^[8]分析了2000—2015年中国青藏高原NPP的时空特征及驱动因素。相较于国际, 国内对NPP的研究虽然起步较晚, 但发展迅速。朴世龙等^[9]利用CASA模型对我国1997年植被NPP进行估算并对其进行时空特征分析, 得出了我国1997年NPP的分布趋势。张锋等^[10]基于CASA模型评估了20年间内蒙古草原NPP的时空特征及其对气候因子的响应特征。张美玲等^[11]结合草原综合顺序分类系统(CSCS)对CASA模型进行改进并对中国草地NPP的时间变化及其影响因子展开研究。

黄河流域是我国重要的生态屏障, 在我国经济社会发展和生态安全方面具有十分重要的地位^[12]。黄河流域生态系统复杂多变, 构成主要以草地为主。作为生态系统的敏感区, 近些年的气候变化和人类活动, 尤其是“退耕还林还草”等生态工程的实施, 使得整个流域的生态系统发生了显著变化。研究黄河流域草地NPP时空动态, 定量分析气候因素、人类活动对NPP影响, 对于黄河流域生态环境保护与高质量发展至关重要。

本文采用MODIS NDVI、太阳辐射、土地覆盖等多源遥感数据及气象数据, 基于改进的CASA模型估算黄河流域草地NPP, 探讨分析其对气候和人类活动等驱动因素的响应。

1 研究区概况

黄河流域介于95°53′E—119°05′E, 32°10′N—41°50′N之间, 地处我国中北部, 西起巴颜喀拉山, 东临渤海, 北抵阴山, 南达秦岭。地势西高东低, 由西向东逐级降低, 大致可分为三级阶梯。黄河流域土地覆盖类型主要以草地、耕地及林地为主, 如2018年流域草地面积占流域总面积的69.44%, 耕地占21.17%, 林地占5.47%, 其它土地覆盖类型占3.92%(图 1)。气候属大陆性气候, 区域差异较明显, 中部属半干旱气候, 西北部为干旱气候, 流域内水资源匮乏, 径流量年际变化大。

2 数据与方法 2.1 数据来源及处理 2.1.1 多源遥感数据

基于改进的CASA模型反演NPP, 需要归一化植被指数(NDVI)、太阳辐射及植被覆盖等数据资料, 本文采用的NDVI、土地覆盖数据资料均来自MODIS影像, 空间分辨率均为0.05°(5600m×5600m), 由NASA/EOS LPDAAC数据分发中心(search.earthdata.nasa.gov)提供。

MODIS NDVI数据, 采用MOD13C2v061数据集, 该数据集时间分辨率为月, 数据格式为HDF-EOS, 提供以每个像素为基础的植被指数, 包含归一化植被指数(NDVI)与增强型植被指数(EVI), 本文提取归一化植被指数。该套数据覆盖黄河流域区域包括216景影像。数据获取后基于Python完成数据转格式、重投影、重采样、裁剪等一系列预处理; 最后采用最大值合成法得到黄河流域2001—2018年间的0.1°×0.1°空间分辨率的逐年NDVI数据资料。

黄河流域土地覆盖数据来自MCD12C1v006数据集。该套数据时间分辨率为年, 数据格式为HDF-EOS。该套土地覆盖数据资料基于Terra/Aqua MODIS反射率数据进行监督分类得出, 采用国际地圈生物圈计划(IGBP)分类系统定义了17种主要土地覆盖类型。本文根据需要将研究区域土地覆盖类型重新分类为9个主要类别(表 1)。

太阳辐射数据采用CERES_EBAF_Edition 4.1数据集, 来源于Atmospheric Science Data Center(https://eosweb.larc.nasa.gov/project/ceres/ceres_table), 数据格式为netCDF, 时间分辨率为月。该数据集是云和地球辐射能量系统能量平衡和填充(EBAF)大气层顶部(TOA)和地表月平均值数据。数据获取后采用Python进行月尺度数据提取及格式转换, 并通过ArcGIS平台采用最邻近法重采样得到0.1°空间分辨率数据。

为验证CASA模型反演NPP精度, 本文采用MODIS NPP (MOD17A3HGFv061)数据集对模型反演结果进行精度验证分析。MOD17A3数据集时间分辨率为年, 数据格式为HDF-EOS, 共18幅图像。该数据集是基于BIOME-BGC生态系统过程模型估算得到的全球每年植被净初级生产力。数据获取后进行影像拼接及重采样, 得到逐年0.1°辐射影像数据。

2.1.2 气象数据及基础地理数据

本文所使用的降水、气温等气象数据来自于国家青藏高原科学数据中心^[13], 该数据由Yang等^[14—15]提供。中国区域地面气象要素驱动数据集是以国际上现有的Princeton再分析资料、GEWEX-SRB辐射资料、GLDAS同化资料, 以及TRMM降水资料为背景场, 融合了中国气象局常规气象观测数据制作而成。该数据集已去除非物理领域的值, 使用ANU-Spline统计插值得到, 格式为NETCDF, 时间分辨率为月, 空间分辨率为0.1°。数据精度优于国际上已有再分析数据。

研究区域的省份边界数据来源于中国气象局国家气象信息中心(http://data.cma.cn)。

2.2 研究方法 2.2.1 实际NPP(ANPP)的估算

本文选用朱文泉等^[16]改进的Carnegie-Ames-Stanford approach (CASA)模型估算黄河流域草地NPP。CASA模型是由NDVI、气象数据、植被类型等共同驱动的光能利用率模型, 计算参数少, 计算较简单, 已经广泛应用于NPP估算^[17], 其计算公式如下：

(1)

式中, APAR(x, t)代表植被的有效光和辐射, 其计算公式如下：

(2)

式中, SQL(x, t)、FPAR(x, t)分别代表像元x在t月份的太阳总辐射、吸收的光合有效辐射; 常数0.5表示植被所能利用的太阳有效辐射占太阳总辐射的比例。其中FPAR(x, t)的计算公式如下：

(3)

式中, NDVI(x, t)表示像元x, 在t月份的NDVI值; NDVI_{i, min}及NDVI_{i, max}代表第i种植被类型的NDVI最小值与最大值; FPAR_max, FPAR_min取值为常数, 分别为0.01、0.95。

已有研究表明^[18]FPAR与比值植被指数(SR)存在一定的线性关系, 且基于NDVI与SR方法计算的FPAR的平均值与实际值的误差最小。计算公式如下：

(4)

(5)

(6)

光能利用率(ε)表示植被利用光合作用把所吸收的有效光合辐射转换为有机碳的效率。主要由气温和水分确定, 其计算公式如下：

(7)

式中, T_ε1x, t、T_ε2x, t、W_ε(x, t)分别代表气温和水分对光能利用率的影响; ε_max表示植被在理想状态下的最大光能利用率, 其取值和植被类型相关。本文ε_max采用朱文泉等^[19]利用中国NPP实测数据模拟出的草地最大光能利用率(0.54)进行计算。

2.2.2 潜在NPP(PNPP)及人为NPP(HNPP)的计算

潜在NPP采用Miami模型^[20—21]计算, 其计算结果表示气候因子对实际NPP (ANPP)的影响。Miami模型由年降水量和年均气温驱动, 计算公式如下：

(8)

(9)

(10)

式中, NPP_T、NPP_R分别表示通过年均气温和年降水量得到的潜在植被净初级生产力; T表示年均气温; R表示年降水量。

人为NPP (HNPP)通过植被实际生产力与植被潜在生产力之差进行计算, 基于该指数衡量人类活动对实际NPP (ANPP)的影响, 其计算公式如下：

(11)

2.2.3 一元线性回归分析及其检验方法

本文采用基于最小二乘理论的回归分析方法, 计算逐像元草地NPP年变化趋势。该方法已被广泛应用于时间序列数据的变化趋势计算^[22], 计算公式为：

(12)

采用F检验法对一元线性回归分析结果进行显著性分析：

(13)

2.2.4 突变检验

非参数统计检验Mann-Kendall(M-K)方法^[23]被用来检测黄河流域2001—2018草地NPP突变点。M-K方法计算过程如下：

首先, 对具有n个样本的时间序列x₁, x₂, .., x_n构造一个秩序列, 秩序列S_k是第i个时刻数据大于j时刻时, 数值个数的累加, 其计算公式为：

(14)

(15)

在时间序列随机独立的假设下, 定义统计量：

(16)

式中, 当k=1时, UF₁=0;E(s_k)、Var(s_k)分别是s_k的均值和方差, 计算公式分别为：

(17)

(18)

重复上述过程, 计算UF_k的反序列得到UB_k, 同样UB₁=0:

(19)

若UB_k和UF_k两条曲线出现交点, 且交点在临界线之间, 可认为交点对应的时刻便是突变开始的时间。本文取显著性水平0.05, 则当|UF_k|>1.96时, 说明有显著趋势。

2.2.5 滑动T检验

滑动T检验是通过检验两个样本的均值是否有显著性差异来检验突变^[24], 基本思想是：将两个子集均值是否有显著性差异看作两个总体均值是否有显著性差异的问题进行检验。设滑动点前后, 两个子序列的分布函数为F₁(x)和F₂(x), 从中分别抽取大小为n₁和n₂的两个样本, 要求检验假设：F₁(x)=F₂(x), 则其计算公式为：

(20)

其中：

(21)

(22)

(23)

(24)

自由度计算公式为：

(25)

T服从t(v)分布, 确定显著性水平后, 查询t分布表得到临界值t_α/2。当T>t_α/2时, 拒绝原假设, 存在显著性差异; 当T < t_α/2时, 则接受原假设, 不存在显著性差异。

2.2.6 偏差分析

偏差分析可以表示各年NPP偏离平均值的程度^[25]。定义为逐像元计算NPP值与多年NPP均值的差值, 其计算公式为：

(26)

式中, G_ij表示第i行, 第j列的NPP偏差值, NPP_ij^t指的是t年第i行, 第j列的NPP值, NPP_ij表示第i行, 第j列多年NPP均值, 本文数据时间跨度为18a。

2.2.7 相关性分析

本文采用Pearson方法来分析黄河流域2001—2018年草地NPP与气候因子的相关性, 其计算公式为：

(27)

式中, r_xy为要素x与y之间的相关系数, n为研究时间跨度(18a), x_i为第i年的NPP值, y_i为第i年的年均气温或降水量。

3 CASA ANPP模拟结果检验

由于研究区域NPP实测值较难获取, 本文采用MODIS NPP影像对改进CASA模型模拟实际净初级生产力(ANPP)结果进行验证。对比分析改进CASA模型模拟黄河流域草地区域ANPP和对应时段的黄河流域草地区域的MODIS NPP(表 2)可知改进CASA模型模拟结果精度较高, 模拟的2001—2018年ANPP结果与MODIS NPP的差值集中分布在1.67—29.91gC m^-2 a^-1间, 最小差值出现在2015年, 最大差值出现在2003年。整体上, 本文基于改进CASA模型模拟实际净初级生产力与MODIS定量反演NPP值的最大差值(29.91 g C m^-2 a^-1)显著低于其他学者的研究结果, 如黄悦悦等^[26]基于CASA模型估算华北平原NPP与MOD17A3 NPP对比分析的最大差值为147.95 gC m^-2 a^-1, 王建文^[27]模拟图们江跨国界地区NPP与MOD17A3 NPP对比分析的最大年均差值为114.19 g C m^-2 a^-1。此外, 改进CASA模拟ANPP结果与MOD17A3 NPP呈显著相关性(图 2), R²为0.70(P < 0.01)。整体上, 本文基于改进CASA模型估算的黄河流域草地实际净初级生产力精度较高。

4 结果与分析 4.1 黄河流域草地NPP突变分析

对CASA模拟实际净初级生产力ANPP和潜在净初级生产力PNPP进行MK和滑动T突变分析, 精确确定草地NPP突变点。根据图 3UF及UB曲线交点可知ANPP在2012、2013、2016年可能发生突变, 因此为进一步确定突变点对ANPP进行滑动T检验(图 3)得出：ANPP显著突变点为2011及2013年, 但2011年不为可能突变点, 结合MK和滑动T突变分析结果确定黄河流域草地实际净初级生产力ANPP在2013年发生突变。对PNPP进行MK突变分析发现黄河流域潜在净初级生产力PNPP在2013年发生突变(图 3), 即ANPP与PNPP均于2013年产生突变。综上, 结合MK突变和滑动T突变分析最终确定黄河流域草地净初级生产力NPP(包括ANPP和PNPP)在2013年发生突变。

4.2 黄河流域草地NPP时空动态变化特征 4.2.1 ANPP, PNPP, HNPP时空趋势分析

逐像元进行NPP年变化趋势及趋势显著性分析。由实际净初级生产力年变化趋势(图 4)可知2001—2013年间, 约92.66%区域的ANPP呈增加趋势, 主要分布在陕西及山西省北部、内蒙古自治区东南部、青海省东南部以及四川省北部, 其中约15.08%区域通过显著性检验(P < 0.05);ANPP减少区域零星分布在河南省北部及内蒙古自治区南部, 约占流域的7.34%, 其中约0.50%区域通过显著性检验。2014—2018年间, ANPP增加区域与减少区域分别为57.77%、42.23%, 其中通过显著性检验的区域分别为2.36%, 4.31%。相对2001—2013时段, 2014—2018时段ANPP减少区域增加, 主要分布在内蒙古自治区东南部、山西省及河南省北部、甘肃省西南部。

2001—2013年潜在净初级生产力PNPP以增加趋势为主(图 4), 约69.00%区域呈增加趋势, 主要分布在宁夏回族自治区南部、甘肃省北部、陕西省中部以及内蒙古自治区东南部, 其中约7.51%区域通过显著性检验; 减少区域占流域总面积的31.00%, 其中显著减少区域约占0.42%。2014—2018年PNPP呈增加趋势, 增加区域占比75.33%, 减少区域占比24.67%, 二者较前一阶段前者区域占比增加6.33%, 后者区域占比减少6.33%。

2001—2013年人为净初级生产力HNPP整体呈增加趋势(图 4), 但趋势较弱, 增加区域约占研究区域的37.21%, 主要分布在山西省及陕西省北部, 其中约1.44%区域通过显著性检验; 减少区域占62.79%, 主要分布在青海省及甘肃省北部、宁夏回族自治区及内蒙古自治区南部, 其中约1.35%区域通过显著性检验。2014—2018年, HNPP呈增加趋势, 约占研究区域的77.73%, 较前一阶段增加40.52%;减少区域约占22.27%, 较前一阶段减少40.52%。

4.3 气候和人类活动对黄河流域NPP动态的影响 4.3.1 气候和人类因素对NPP影响的空间分析

根据表 3定义的方法, 逐像元比较3种不同类型NPP的年际变化斜率, 确定黄河流域2001—2018年ANPP变化的驱动因素并将分析结果定义为五类情况：气候因素导致ANPP增加(CAI)、人类活动因素导致ANPP增加(HAI)、气候因素导致ANPP减少(CAR)、人类活动因素导致ANPP减少(HAR)以及ANPP无变化(NV)。

由图 5可以得出：2001—2013年ANPP呈增加趋势, 其中由气候因素导致ANPP增加的区域占51.63%, 由人类因素导致ANPP增加的区域占41.03%;ANPP减少区域占7.34%, 气候因素占比4.37%, 人类活动因素占比2.97%。由图 5可知：2014—2018年由人类活动导致ANPP减少的面积较前一阶段增加了34.89%, 主要分布在四川省东北部、青海省北部及南部、甘肃省西南部及陕西省中部。整体上, 2014—2018年人类活动对ANPP增长的抑制作用逐渐增强。

4.3.2 气候因子对NPP的影响

已有研究表明降水量是影响植被生长的重要因素。气候因素对中国植被NPP的影响存在区域差异, 植被NPP与降水、气温、光照等气候要素相关, 如降水量是影响内蒙古植被净初级生产力的重要因素^[28]; 长江流域NPP、NDVI与气温有较强相关性^[29]。基于前人的研究结论, 本文重点分析黄河流域降水、气温等气候要素对其草地NPP的影响。

2001—2018年间, 黄河流域降水量呈增加趋势(图 6), 速率为2.94 mm/a, 且在2003年达到最大值(593.76 mm), 2001年达到最小值(417.01 mm); 气温亦呈增加趋势, 年变化速率为0.17 ℃/10a, 于2013年达到最大值(4.83℃)。2001—2013年间ANPP与降水呈正相关关系的区域占研究区域总面积的97.79%, 显著正相关区域约占45.01%, 负相关关系区域占比2.21%(图 7); ANPP与气温呈正相关关系的区域约占研究区域总面积的29.87%, 负相关关系的区域约占70.13%, 其中显著负相关区域约占3.39%(图 7); 2014—2018年间ANPP与降水呈正相关关系的区域占79.08%, 显著正相关区域约占14.09%, 负相关区域占比20.92%, 显著负相关区域占比20.07%(图 7); ANPP与气温呈正相关关系的区域占17.11%, 负相关关系的区域占82.96%(图 7)。

整体上黄河流域草地实际净初级生产力ANPP与降水量的相关关系更为显著, 这一研究结果与陈强等^[30]研究发现2001—2010年黄河流域NPP变化与气候因素关系相一致。

4.3.3 人类活动对NPP的影响

气候因子是影响NPP的重要因素, 人类活动对其影响也不可忽略, 如已有研究表明耕地面积变化、草地恢复和退化面积等因素也是影响NPP的重要人类活动要素^[31—32]。本文基于黄河流域逐年土地利用/覆盖数据资料, 分析2001—2018年黄河流域草地、耕地面积动态, 讨论其对NPP的影响。

黄河中上游地区是退耕还林还草工程的重点实施区域。由图 8黄河流域2001—2018年不同土地利用类型的面积年际分布变化及表 4黄河流域2001—2018年土地利用转移矩阵可知, 2001—2013年间黄河流域草地面积快速减少, 2014—2018年间草地面积呈微弱增多趋势; 2001—2018年土地利用类型的转换主要发生在城市建设用地、草地、耕地及其它植被类型之间, 草地主要向耕地转出, 耕地面积的增加主要来自于草地, 林地的转入, 城市和建设用地于2013年之后增加趋势加强。2012—2018年间城市和建设用地及耕地呈现增加趋势, 且耕地在2005—2013年间增加更为迅速, 城市和建设用地于2013—2016年迅速增加, 另外黄河流域草地面积逐年减少, 2014—2018年ANPP增加速率较2001—2013年降低, 人为活动驱动ANPP减少的占比明显增加。这一现象说明城市和建设用地及耕地影响了黄河流域草地的生长。

5 结论

本文基于CASA模型估算黄河流域草地实际净初级生产力, 讨论气候变化和人类活动对黄河流域草地生态系统的影响, 主要结论：

2001—2018年间, 黄河流域净初级生产力(ANPP)、潜在净初级生产力(PNPP)及人为净初级生产力(HNPP)均呈现增加趋势, 其中HNPP增加趋势较弱。时序上在2013年发生显著突变。2001—2013年间, ANPP及PNPP明显增加, 而HNPP呈现减少趋势, 减少区域主要分布在黄河流域中部; 2014—2018年间ANPP增加减慢, HNPP及PNPP明显增加, 增加区域主要分布在黄河流域北部。

黄河流域草地净初级生产力时空动态受气候变化和人类活动共同作用。2001—2018时段, 气候变化是影响实际净初级生产力的主要因素, 而不同时段影响要素有些许差异, 2001—2013年间气候因素及人类活动均促进ANPP的增加; 2014—2018年间由于城市和建设用地面积的显著增加以及草地面积的逐渐减少, 人类活动因素对实际净初级生产力的抑制作用逐渐增强, 进而导致了ANPP增加速率的降低。

本文结合气象数据及遥感数据估算黄河流域草地NPP时空变化特征, 量化了气候因素和人类活动对草地NPP的影响, 可为草原管理政策及今后的研究提供有价值的参考。相对已有研究, 本文考虑了研究区域逐年草地面积变化对其净初级生产力的影响, 但土壤、地形等自然要素, 及城市扩张、国家政策等人类活动因素都是影响NPP的主要因子。因此NPP时空动态的驱动因素分析仍需进一步研究。