2000-2020 年石羊河流域草地净初级生产力时空动态及其对气候的响应

李怀海，李纯斌，吴 静，吉珍霞，秦格霞，马娟娟

（甘肃农业大学资源与环境学院, 甘肃 兰州 730070）

草地作为生态系统的重要组成成分，涵盖了地球陆地表面约20%的地域，是全世界广泛分布的主要植被类型之一。因其覆盖面积大、种类繁多，在区域乃至世界碳循环研究中具有特殊的地位和非凡的生态意义[1]。植被净初级生产力(net primary productivity, NPP)是绿色植物光合作用产生的有机物总量扣除自养呼吸后的剩余部分[2]。NPP 作为考量气候变化的核心内容之一，是决定草地生态系统健康稳定和生产能力高低的重要指标[3]。与此同时，草地NPP 的研究对于合理利用草地资源，开发草地生产潜力，实现草地产量最大化具有重要意义[4]。草地NPP 受人为和自然等因素的影响，造成草地NPP 年际波动的主要驱动力之一是气候变化[5]。近100 年来，全球范围内气温呈现升高趋势，干旱半干旱地区气温升高尤为明显。研究表明，全球变暖贡献率的44%来自于半干旱地区的气温上升[6]。因此，探讨河西典型干旱内陆河流域的气候变化对西北草地生态系统的科学管理具有重要指导作用。

CASA (carnegie-ames-stanford approach)模型是随着遥感技术的发展而产生的草地NPP 估计模型[7]，具有所需参数少，与植被生理特征高度相关等优点，是目前世界上最常用的草地NPP 估算模型之一[8]。近年来，有很多学者将CASA 模型用于陆地生态系统碳循环及其对气候变化响应等研究领域，如朱莹莹等[9]对我国西北地区草地NPP 进行因子分析，发现草地NPP 对温度和降水的响应存在一定的空间差异；李金珂等[10]对秦巴山区的NPP 时空变化特征及相关因子进行解析，得出自然因子对NPP 的贡献率最高的是实际蒸散量，其次为降水、气温，坡度对其贡献率最小；人为因子对NPP 具有积极和消极两方面作用。也有学者对CASA 模型进行了不同程度的改进，如朱文泉等[11]引入分类精度和植被覆盖度类型图改进CASA 模型，估计了中国陆地上不同植被覆盖度下的时序草地NPP；杨勇等[12]改进CASA 模型最大光能利用率及简化模型参数实现了内蒙古小尺度植被NPP 估测。但这些改进都没能实现草地NPP 与草地类型的结合。

草原综合顺序分类法(comprehensive and sequential classification system, CSCS)是任继周[13]提出的一种针对草地的植被分类方法，其依据草地水热条件变化对草地分类体系进行了改进，具有特有的先进性与创新性。以≥0 ℃年积温、湿润度(k)等生物气候指标为基础，对积温和润湿度进行耦合分类，并将具有相同热量和湿润度水平的集成为一类，避免主观因素对草地划分的影响[14-15]，其划分出的潜在植被(potential nature vegetation, PNV)类型能够反映出区域内现有植被的发展趋势[16-17]。与另一种传统分类系统(中国草地分类系统)相比，具有定量化优势，在科学性上处于世界草地分类科学的前沿[18]。

基于CSCS 法这种生物气候指标的草地分类系统对CASA 模型进行改进，能够很好地弥补草地类型与草地NPP 脱节问题。

鉴于此，为进一步深入探讨石羊河流域长期序列草地NPP 时空变化特点及其对气候变化的响应，利用张美玲等[19]改进的CASA 模型，结合CSCS 严谨的植被发生学理论与严谨的分类体系，不仅可以体现草地发生的水分平衡因素(湿润度)和热量因素(积温)，而且能够揭示草地NPP 与草地类型之间所存在的内在关联。此外，CSCS 方法能够避免因主观因素对草地类型划分结果的不同而影响到草地生物量大小，其为地带性草地潜力研究提供了新思路，更重要的是其对石羊河流域草地生态系统碳循环等相关研究具有重要意义，同时也能为草地保护和草地生态系统稳定发展提供科学支撑。

1 材料与方法

1.1 研究区概况

石羊河流域属于河西走廊地区三大主要内陆河流域之一，地处河西走廊东段、祁连山北麓(101°41′～104°16′ E, 36°29′～39°27′ N)，位于黄土高原、蒙新高原和青藏高原重要交汇过渡地带，流域范围为4.16 × 104km2(图1a)。行政区划包括3 市8 县(区)，分别为金昌市1 区(金川区) 1 县(永昌县)、武威市1 区(凉州区) 3 县(古浪县、民勤县、天祝县)、张掖市2 县(肃南裕固族自治县和山丹县部分地区)。地势西南高、东北低，地貌类型自南向北由山地向平原、丘陵与荒漠过渡(图1b)。气候以大陆性干燥气候为主，蒸发较强，全年内降水少且集中，水资源供求矛盾突出。

图1 石羊河流域位置(a)及研究区信息(b)Figure 1 Shiyang River Basin location (a) and study area information (b)

1.2 数据来源及预处理

1.2.1 高程数据及气象数据

高程数据采用SRTM 90 m DEM 产品，来源于地理空间数据云(http://www.gscloud.cn/)。河流系统和流域边界数据来自中国西部环境与生态科学数据中心提供的石羊河流域信息系统数据集(http://westdc.westgis.ac.cn)。气象资料来源于中国气象科学数据共享服务网(http://data.cma.cn)，包含气温、降水、日照时数等，时间跨度为1991-2020 年，利用MATLAB 9.5 软件进行处理，逐日累加降水合成年降水量，而年积温是连续5 d 或超过5 d 大于0 ℃的日平均气温的累积。利用专业气象插值软件ANUSPLINE4.3 将积温、降水、太阳辐射等插值为空间分辨率为90 m × 90 m 的栅格数据集[20]。所有数据均采用GCS-WGS-1984 地理坐标系统且投影为WGS-1984-Albers，并在ArcGIS10.2 中按照研究区范围进行镶嵌、裁剪。

1.2.2 遥感数据

MOD13Q1 NDVI 产品源于NASA 对地观测系统数据共享平台(EOS data gateway，http://delenn.gsfc.nasa.gov/-imswww/pub/imswel-come/)，空间分辨率为250 m，时间分辨率为16 d，时间范围是2000-2020 年。借助MRT (modis re-projection tool, MRT)完成数据格式转换和投影，在ArcGIS10.2 中重采样成与气象数据等同的栅格大小，采用栅格计算器将其乘以比例因子0.000 1 得到真实的归一化植被指数(normalised vegetation index，NDVI) 值，并运用最大值合成法(maximum value composites, MVC)对NDVI 数据进行月值合成，该方法能有效去除大气、云、太阳高度等因素的影响[21]。

1.3 研究方法

1.3.1 CASA 模型

CASA 模型是Potter 等[7]建立的光能利用率模型的典型代表，CASA 模型是由温度、降水、太阳辐射、植被类型等因素驱动的光能利用模型[22]。CASA模型之中的NPP 是植被吸收的光合有效辐射(APAR)、最大光能利用率( εmax)、温度胁迫系数(Tε1、Tε2)和水分胁迫系数(Wε)的函数[7]，其计算公式如下：

式中：x为对应栅格像元；t表示对应月份。

本研究采用张美玲等[19,23]改进的CASA 模型，该模型主要是从水分胁迫系数(Wε)和最大光能利用率( εmax) 两个方面进行改进，将计算Wε的复杂土壤参数仅用 ＞ 0 ℃的年累积温、湿度k 值代替，实现了CASA 模型和CSCS 分类系统的耦合。最大光能利用率 εmax参考张美玲等[23]对草原综合顺序分类最大光能利用率的模拟值，其中草甸类草地的最大光能利用效率利用包刚等[24]对内蒙古草原植被最大光能利用效率结果校正。水分胁迫系数Wε计算公式为：

1.3.2 Sen 的斜率估计分析和Mann-Kendall (M-K)突变检验

基于栅格尺度，采用Sen 的斜率估计法对流域2000-2020 年NPP 的斜率进行了测算，斜率可以代表NPP 的平均变化率和此时的变化趋势。当β＞0 时，序列呈递增趋势；当β= 0 时，序列的趋势不明显；当β＜ 0 时，序列呈下降趋势[25]。

Mann-Kendall (M-K)突变检验是一种无分布式检验方法，选取的样本具有易于计算、受其他因素影响较小且不必满足一定的分布原则等优点，能够反映样本的长期变化趋势和变异情况，是一种常用的突变检测方法[26]。将得到的NPP 趋势栅格基于M-K 统计检验法进行NPP 变化趋势显著性检验[27]。

1.3.3 NPP 稳定性分析

基于像元尺度，利用变异系数(coefficient of variation, CV)来表征草地NPP 的受扰动幅度，值越小表示草地生长状态越稳定；数值越大，草地的生长越不稳定[28]。计算公式如下：

式 中：CV为NPP 变 异 系 数；i表 示 第i年 对 应的NPP 值；NPPmean为 基 于2000-2020 年 得 到 的NPP数据多年均值。

1.3.4 未来趋势分析

Hurst 指数被广泛应用于气候、水文、地质等研究领域，它可以定量描述时序草地NPP 的自相似性和长期依赖性，并可以预测未来NPP 的变化趋势[29]。Hurst 指数(H)与时序草地NPP 之间的趋势表现特征如下：当0 ＜ H ＜ 0.5 时，说明时序草地NPP 具有反持续性，即未来的趋势与过去的趋势相反，值越接近0，其反持续性越强；当H = 0.5 时，时序草地NPP 具有随机序列、有限方差和无长期相关性的特点；当0.5 ＜ H ＜ 1 时，说明时序草地NPP 是具有长期相关特征的连续序列，未来变化趋势与过去趋势一致，其值越接近1，持续性越强[29]。

1.3.5 偏相关分析

偏相关分析可以在排除其他变量影响前提下，分析两个变量之间的线性相关关系[28]。草地NPP受积温和降水影响程度大小可通过相应的偏相关系数来表达，公式如下：

式中：x表示NPP；y代表积温；z为降水量；rxy表示年均草地NPP 与积温的相关系数；rxz表示NPP 与年降水量的相关系数；ryz表示年积温与年降水量的相关系数。rxy(z)为变量降水固定后年NPP 与积温的偏相关系数，以此类推。

2 结果与分析

2.1 石羊河CSCS 草地分类结果及NPP 验证

表1 石羊河流域CSCS 草地类划分Table 1 Classification of CSCS grassland in Shiyang River Basin

利用2006 年甘肃草原15 个调查采样点(包括经纬度、草地类型、总产草量等)和2020 年的16 个实测采样点对NPP 数据进行了验证(图2b)。实测采样时，在相距大于3 km 的不同草原类型上随机建立样地(荒漠类7 个、半荒漠类7 个、草原类9 个、草甸类8 个)，在10 m × 10 m 的样地中设置5 个样方，第1 个位于样地中心，其余在四角，样方按1 m × 1 m布置。采用对角线采样法[30]，将在地面上获取的鲜草放入烘箱，烘干至恒重，取平均值作为地上生物量。由于总生物量包括地上和地下两部分，地上生物量乘以比例系数可得到地下生物量，其中，不同草地类型的地上地下部分生物量相关比例系数参考文献[31-32] 得出。最后将两者相加后乘以系数0.475[33]转 换 成 以 碳 为 单 位[g·(m2·a)-1] 的NPP 总量。相关性分析显示，石羊河草地NPP 模拟值与实测值的Pearson 相关性系数为0.89，存在显著正相关关系(R2= 0.787 2，P＜ 0.01)，基本满足研究精度要求。

图2 CSCS 草地分类结果(a)及NPP 模拟值和实测值相关性分析(b)Figure 2 CSCS grassland classification results (a) and correlation analysis between NPP simulated and measured values (b)

2.2 石羊河草地NPP 的年际变化特征

2000-2020 年间，石羊河流域草地年均NPP 值为170.24 g·(m2·a)-1(图3a)，0～70 g·(m2·a)-1所占面积比例最大，＞ 350 g·(m2·a)-1的面积占比最小。年均NPP 的波动区间为137.63～211.18 g·(m2·a)-1(图3b)，其中，年均NPP 最高值和最低值的年份分别出现在2019 年和2001 年。

图3 草地NPP 时空变化与检验Figure 3 Temporal and spatial variation and test of grassland NPP

草地类型上(图3c)，山地草甸类(ⅡE30)年均NPP 最高，为548.74 g·(m2·a)-1；山地草甸草原类(ⅡD23)次之，为454.50 g·(m2·a)-1；温带荒漠类(ⅣA4)年 均NPP 最 低，为91.65 g·(m2·a)-1。多 类 草 地在2009 年出现NPP 低值，高值区出现年份2012 年、2019 年、2020 年，总体变化趋势整体呈现NPP 波动上升趋势。

采用M-K 突变检验进一步探究21 年间石羊河流域草地NPP 变化(图3d)：UF 统计量呈现波动上升，但在2016 年之前均未通过0.05 的显著性水平线，说明在此期间的草地NPP 相对稳定。在2013 年后UF 曲线呈现直线上升，UF 和UB 曲线在2015和2016 年间相交，出现突变节点，打破了原有的波动格局，且在2016 年后UF 曲线超过0.05 显著水平检验，出现草地NPP 快速增加的态势。结合NPP 离差可看出，在2014 年之前，NPP 离差值多为负，2014年后都为正值，印证了M-K 突变检验的结果。

2.3 石羊河草地NPP 空间分布及动态

基于Sen’s 和M-K 检验法，得出21 年间草地类NPP 变化趋势(图4a)，其NPP 平均增量为28.96 g·[m2·(10 a)]-1，整体呈现波动上升趋势。NPP 增加的面积(不显著增加、显著增加、极显著增加)占整个流域面积的72.82%，其中，极显著增加面积占27.63%，主要分布在祁连山及流域上游外围地区。NPP 减少区域(不显著减少、显著减少、极显减少)主要分布在流域中部及北部民勤盆地，特别是在武威市的中西部地区NPP 呈现显著性减少趋势，主要受巴丹吉林沙漠和腾格里沙漠的侵袭，导致植被生长受限，出现一定程度的退化特征(图4b)。

不同类型草地NPP 变化趋势所占比例上(图4c)，NPP 极显著增加表现出：ⅢC17 最大，占比89.64%；ⅢB10 次之，占比77.25%；ⅣA4 最小，占比6.10%。NPP 显著增加：ⅠE29 最大，占比44.60%；ⅡE30 次之，占比38.49%；ⅣA4 最小，占比4.37%。NPP 极显著减少和显著减少的草地类都表现出：ⅣA4 最大，总占比4.72%，退化趋势明显。

图4 NPP 变化趋势显著性检验与各草地类变化统计Figure 4 Significance test of NPP trends and statistics of various grassland types

2.4 石羊河草地NPP 稳定性分析

石羊河流域草地NPP 变异系数CV 介于0.01～4.12，平均变异系数为0.24，在空间上存在显著差异(P ＜0.05) (图5a)。高稳定区(CV＜0.1)的面积占比0.31%，仅分布在下游荒漠地区，该地植被覆盖低，整体变化小；低波动区(0.1 ≤ CV ≤ 0.15)面积占比为4.74%，主要分布在流域下游区域，常年受两大沙漠影响，植被稀疏，整体变化较小；中波动区(0.15＜CV≤0.3)面积占比为85.30%，广泛分布在流域的中下游地区；高波动区(CV ＞ 0.3)面积占比为9.64%，集中分布在武威盆地与民勤盆地交汇地带，受自然和人为干扰强烈，NPP 变化显著(P ＜0.05)。

图5 草地NPP 变异系数与各草地类稳定性占比Figure 5 Variation coefficient of grassland NPP and proportion of stability of each grassland type

草地类型上(图5b)，仅ⅢC17 存在高波动，中波动占据主体地位，说明草地变化较稳定，对高波动区域应持续关注，防止其恶化。

2.5 草地NPP 变化趋势预测

2000-2020 年流域草地NPP Hurst 指数(H)均值为0.63 (图6a)，表明未来一段时间草地NPP 总体呈现与此21 年相同变化趋势(P＜ 0.05)，即未来一段时间的总体草地NPP 将会增加。

图6 2000-2020 年草地NPP 可持续分析Figure 6 Sustainable analysis of grassland NPP from 2000 to 2020

在ArcGIS 中叠加NPP 变化趋势栅格和Hurst指数栅格得到未来草地变化预测图(图6b)。结合表2 得出：H 指数小于0.5 的区域中，反持续极显著增加(P＜ 0.01)面积占比为10.86%，主要分布在流域西南部四县交汇地带，该地未来草地NPP 增加放缓。反持续显著减少(P＜ 0.05)和反持续极显著减少(P＜ 0.01)占比都很小，在武威市中部零星分布；H 指数大于0.5 的地区，持续未发生变化的区域主要分布在流域下游，占流域一半以上区域(51.8%)，其次为持续极显著增加(P＜ 0.01)的区域占比最大(16.72%)，集中分布于古浪县，未来该地草地会有一定程度上的改善；NPP 持续显著增加(P＜ 0.05)的区域主要分布在肃南县和天祝县西南部。NPP 预测类型结果编码如表2 所列。

表2 石羊河流域草地NPP 变化趋势预测Table 2 Prediction of grassland NPP change trend in Shiyang River Basin

2.6 气候变化对草地NPP 的影响

石羊河流域草地NPP 与积温的关系图中可看出(图7a)，除2009、2018 和2019 年NPP 与积温呈相同变化趋势外，其他年间均呈现相反变化趋势，说明石羊河流域积温升高在一定程度上对草地NPP增加产生不利影响；草地NPP 与降水的变化整体呈现一致性，协同性较好(图7b)，即降水量的增加有利于石羊河流域草地NPP 的增加，具有明显的促进作用。

图7 2000-2020 年草地NPP 与气候年际变化关系Figure 7 Relationship between grassland NPP and interannual climate change from 2000 to 2020

通过偏相关分析探讨了草地NPP 与气候变化的响应关系。分析得出NPP 与降水、积温的偏相关系数分别分布在-0.84～0.94 和-0.79～0.88。在控制积温的情况下(图8a、b)，NPP 与降水呈正相关区域面积占比为70.89%。其中，显著相关(P＜ 0.05)区域面积占24.93%，主要分布在流域上游的祁连山地区。NPP 与降水呈负相关的区域占29.11%，通过显著性检验区域仅占流域总面积的2.78% (P＜ 0.05)，这可能是因为降水增加幅度相对较小，未能抵消积温升高，蒸散发较大造成的水分损耗，而呈现负相关；将通过显著性检验的栅格与降水栅格做叠加分析发现，显著性正相关的区域基本和380 mm 降水线持平(P＜ 0.05)，说明在大于该降水量的地区草地NPP 主要受降水控制；在控制降水的情况下，草地NPP 与积温的关系也存在正负相关并存的情况(图8c、d)，NPP 与积温呈现正相关的区域占流域面积的54.62%，其中通过显著性检验的面积占4.43% (P＜0.05)，主要分布在流域上游；呈现负相关的面积占总面积的45.38%，通过显著性检验的面积占0.95%(P＜ 0.05)，主要分布在流域的中部地区。

图8 NPP 与降水、积温偏相关分析及显著性检验Figure 8 Partial correlation analysis and significance test of NPP with precipitation and accumulated temperature

草地类型上，NPP 对积温响应最为敏感的是草甸类草地，敏感性大小依次为ⅠE29 ＞ ⅡE30 ＞ ⅡD23，通过显著性检验的面积占比依次为51.69%、7.95%和2.64%。荒漠类草地与积温呈现一定的负相关关系，如ⅣA4 草地类。除两类荒漠草地(ⅢA3 和ⅣA4)对降水响应关系不明显外，其他草地对降水都存在较为明显的相关关系，其中ⅠE29 显著性面积占比为85.88%，其次为ⅡD23 (84.44%)和ⅡE30 (81.37%)。

3 讨论

草地NPP 变化受气候和人类活动的影响，特别是温度和降水是影响植被光合作用的主要因素。NPP 大小与水热条件的时空分布差异密切相关[34]。流域草地NPP 主要受降水影响，在降水量大于380 mm的地区，NPP 与降水有显著相关性(P＜ 0.05)。当降水减少时，流域会有高达70.89% 的区域呈现NPP减少，24.93%的区域NPP 出现显著减少(P＜ 0.05)，这与在干旱区降水对NPP 响应更敏感有关，印证了朱莹莹等[9]和Xu 等[35]得出的在西北干旱区降水与NPP 关系相比气温更加密切的结果。温度胁迫在干旱区具有两面性[36]，当积温较高时，在武威盆地、金昌盆地等中下游缺水地区，会增加失水，草地因缺水而生长受限，NPP 降低；在湿润度较高的高寒草甸区，水分较为充足，积温成为草地NPP 增加的限制因素，故该地区积温与NPP 之间呈现正相关关系。另外，流域上游的祁连山区出现NPP 极显著增加，与该地区暖湿气候具有很大关联，这与张雪蕾等[37]研究相一致。此外，祁连山草地NPP 的增加，也可能与近年来造林、退耕还林还草、祁连山国家公园等重大工程的实施有关[38]。

4 结论

本研究基于草原综合顺序分类系统(CSCS)改进的CASA 模型，模拟2000-2020 年石羊河草地NPP，分析了草地NPP 与气候因子的关系。得出以下结论：1) 从时间上看，2000-2020 年流域草地NPP 波动上升，10 年NPP 增量为28.96 g·m-2，尤其2016 年后，NPP 出现明显上涨趋势。在空间上，NPP 呈现南高北低的分布格局，未来一段时间内NPP 还会有所增加，其中，持续极显著增加面积占比16.72%，集中分布在古浪县；反持续极显著增加面积次之，占流域总面积的10.86%，主要分布在流域西南部四县交汇地带。2)不同草地类型中，山地草甸类(ⅡE30) 年均NPP 最高，为548.74 g·(m2·a)-1；山地草甸草原类(ⅡD23)次之，为454.50 g·(m2·a)-1；温带荒漠类(ⅣA4)年均NPP 最低，为91.65 g·(m2·a)-1。流域草地整体稳定，中波动草地占据主体地位，仅ⅢC17 存在较高波动。3)多数草地NPP 增加的主导因素是降水，对降水响应敏感区域面积占比为24.93%，仅两类荒漠草地(ⅢA3 和ⅣA4) 对降水响应关系不明显；草甸类草地NPP 对积温响应最为敏感，敏感性大小依次为ⅠE29 ＞ ⅡE30 ＞ ⅡD23，荒漠类草地与积温呈现一定的负相关关系。

致谢：由衷地感谢甘肃农业大学花立民教授为本研究提供草地调研数据。