黄土高原植被NEP时空变化及对气温降水和退耕还林（草）的响应

Spatiotemporal Variation of Vegetation NEP on Loess Plateau and Its Responses to Temperature， Precipitation， and Grain-for-Green Program

摘要

目的 净生态系统生产力（net ecosystem productivity， NEP）是衡量陆地生态系统碳收支的重要指标，定量评估净生态系统生产力及其影响因素，对于深入理解区域碳循环及其驱动机制具有重要意义。 方法 以黄土高原为研究区，基于MOD13A1 NDVI和MOD17A3 NPP数据产品、气温、降水量、DEM和土地覆盖等数据，利用Carnegie-Ames-Stanford-Approach（CASA）模型、Theil-Sen趋势分析、Mann-Kendall显著性检验和相关性分析等方法，分析2001—2022年黄土高原植被NEP时空分布变化特征，并探讨其对气温、降水量和退耕还林（草）工程实施的响应。 结果 1）2001—2022年黄土高原植被NEP多年平均值为327.06 g/m²（以C计），总体呈显著增加趋势，增加速率为7.51 g/（m²·a） （p<0.05）。其中，草地的增速最快。2）空间上，黄土高原植被NEP以增加趋势为主，其中显著增加区域占77.38%，主要分布在黄土高塬沟壑区的中部和黄土丘陵沟壑区的中部；显著减少区域占比仅为0.98%，主要分布在土石山区及河谷平原区的西部、沙地和农灌区的西南部及黄土高塬沟壑区的西部。不同分区中，黄土丘陵沟壑区植被NEP呈显著增加趋势的区域面积占比最高，占该分区面积的93.00%。3）黄土高原植被NEP与年平均气温、年降水量均呈正相关为主，其中呈显著正相关的区域分别占区域总面积的6.28%和26.38%，前者集中分布在黄土高塬沟壑区南部，后者集中分布在黄土高塬沟壑区的北部和黄土丘陵沟壑区东部、沙地和农灌区中部、东北部、土石山区及河谷平原区东北部；呈显著负相关的区域分别占1.38%和0.13%，分别集中在黄土高塬沟壑区西部和西北部。不同分区中农业植被、森林和草地植被NEP与年降水量呈显著正相关的面积比例高于年平均气温，且沙地和农灌区植被NEP与年降水量关系最为明显。4）黄土高原坡度>25°的区域中，耕地转为林地和草地区域的植被NEP增速，高于耕地保持不变的区域。 结论 退耕还林（草）工程的实施及气候暖湿化对黄土高原植被NEP增加起到重要作用。

Abstract

Objective Net ecosystem productivity （NEP） is a key indicator for measuring the carbon budget of terrestrial ecosystems. Quantitatively assessing NEP and its influencing factors is of great importance for deeply understanding the regional carbon cycle and its driving mechanisms. Methods Taking the Loess Plateau as the study area， the research used data products such as MOD13A1 NDVI， MOD17A3 NPP， climate data （temperature and precipitation）， digital elevation model （DEM）， and land cover data. The Carnegie-Ames-Stanford-Approach （CASA） model， Theil-Sen trend analysis， Mann-Kendall significance test， and correlation analysis were employed to analyze the spatiotemporal distribution changes of vegetation NEP on the Loess Plateau from 2001 to 2022 and to explore its responses to temperature， precipitation， and the implementation of the Grain-for-Green Program. Results 1） From 2001 to 2022， the multi-year mean vegetation NEP on the Loess Plateau was 327.06 g/m² （in terms of C）， showing a significant increasing trend at a rate of 7.51 g/（m²·a） （p<0.05）. Among land cover types， grassland had the fastest growth rate. 2） Spatially， vegetation NEP on the Loess Plateau mainly exhibited an increasing trend. Significantly increasing areas accounted for 77.38% of the total， primarily distributed in the central part of the loess tableland gully region and the central part of the loess hilly gully region. Significantly decreasing areas accounted for only 0.98%， primarily located in the western part of the rocky mountainous and river valley plain region， the southwestern part of the sandy land and agricultural irrigation region， and the western part of the loess tableland gully region. Among the subregions， the loess hilly gully region had the highest proportion of areas where vegetation NEP showed a significant increasing trend， accounting for 93.00% of its total area. 3） Vegetation NEP on the Loess Plateau was overall positively correlated with both mean annual temperature and annual precipitation. Areas showing significant positive correlation accounted for 6.28% and 26.38% of the total area， respectively. The former concentrated in the southern part of the loess tableland gully region， while the latter was mainly distributed in the northern part of the loess tableland gully region and the eastern part of the loess hilly gully region， the central and northeastern parts of the sandy land and agricultural irrigation region， and the northeastern part of the rocky mountainous and river valley plain region. Areas showing significant negative correlation accounted for 1.38% and 0.13%， respectively， and were primarily distributed in the western and northwestern parts of the loess plateau gully region. Among different subregions， the proportion of areas where the vegetation NEP of cropland， forest， and grassland showed a significant positive correlation with annual precipitation was higher than that with mean annual temperature. The relationship between vegetation NEP and annual precipitation was most pronounced in the sandy land and agricultural irrigation region. 4） In areas of the Loess Plateau with slopes greater than 25°， the increase rate in vegetation NEP was higher in areas where cropland had been converted to forest or grassland than in areas where cropland remained unchanged. Conclusion The implementation of the Grain-for-Green Program and the warming-wetting climate play critical roles in increasing vegetation NEP on the Loess Plateau.

Graphical abstract

关键词

净生态系统生产力 / CASA模型 / 气候变化 / 退耕还林（草） / 黄土高原

Key words

net ecosystem productivity / CASA model / climate change / Grain-for-Green Program / Loess Plateau

政府间气候变化专门委员会（Intergovernmental Panel on Climate Change， IPCC）气候变化评估报告^［1］指出，气候变暖导致极端气候事件发生的频率和强度增加，碳排放是全球气候变化及其系列生态环境问题的重要驱动因素。陆地生态系统作为碳排放的重要载体，其碳源、碳汇功能研究已成为国内外学者关注的热点问题^［2］。净生态系统生产力（net ecosystem productivity， NEP）是衡量生态系统碳收支状况的关键指标，是生态系统通过光合作用固定的碳量与呼吸作用释放的碳量差值，其数值的正负反映生态系统在一定时期内作为碳汇或碳源的功能状态。近20 a来中国陆地生态系统整体呈碳汇形势^［3-5］，植被NEP呈上升趋势^［5-7］，如京津冀^［8］和黄土高原延河流域^［9］NEP均呈上升趋势。

目前，植被NEP估算的传统方法主要包括样地清查法和涡度相关法等，这些方法精度较高但周期较长，难以在大尺度区域进行估算。相较于传统方法，基于遥感技术的模型方法已被广泛应用于植被NEP估算^［10］。植被NEP估算涉及植被净初级生产力（net primary productivity， NPP）和土壤异氧呼吸2个内容。其中与NPP关系较为密切的模型为光能利用效率模型，该模型是基于植被对太阳辐射的有效转换和利用，使用植被指数、光合效率和光合有效辐射来估算NPP^［5］。CASA模型（carnegie-ames-stanford approach）是应用最广泛和最具代表性的光能利用效率模型之一^［11］。目前，土壤异氧呼吸模型方面，主要是PEI等^［12］为估算青藏高原高寒草原碳平衡所建立的关于气温、降水量等环境因子之间的关系模型，该方法在区域碳汇分布估算方面得到广泛应用^［13-15］。

黄土高原是我国“两屏三带”生态安全建设的重要组成部分^［16］，其植被NEP研究得到较多关注。周怡婷等^［3］研究发现，气候因素是影响黄土高原植被NEP空间变化的主要驱动因子，近年来实施的山水林田湖草沙等保护和修复工程增强了黄土高原生态系统的固碳能力；裴宏泽等^［4］研究表明，降水和土地类型对植被NEP的空间分布影响最为显著；张怡等^［5］研究表明，植被NDVI是驱动黄土高原植被NEP时空变化的主要因素，农田与裸地转为林草地对区域生态系统增汇贡献率最大；管亚兵等^［9］研究指出，其他土地类型向林地的转变对碳汇量提升效果更显著，耕地、草地转为林地后的植被NEP平均值均大于耕地转为草地和草地转为耕地后的植被NEP平均值。

已有研究多侧重于宏观层面黄土高原植被NEP时空变化研究，在退耕还林（草）工程重点实施的25°以上坡耕地退耕与植被NEP变化关系研究方面尚显不足。基于此，本研究以MOD13A1 NDVI和MOD17A3 NPP数据产品、气温、降水量、DEM和土地覆盖数据为基础，采用改进的CASA模型、土壤异养呼吸模型、Theil-Sen趋势分析、Mann-Kendall显著性检验和相关性分析等方法，开展2001—2022年黄土高原植被NEP时空变化及对气温降水和退耕还林（草）的响应研究，定量评估气温、降水量及退耕还林（草）工程实施对黄土高原植被固碳能力的影响，以期为黄土高原植被生态保护与恢复政策科学调整提供支撑。

1 材料与方法

1.1 研究区概况

黄土高原（107°28′~111°15′E，35°21′~39°34′N）位于中国西北部，涉及青海、甘肃、宁夏、内蒙古、山西、陕西、河南7个省区，总面积为62.46×10⁴ km²。黄土高原由东南向西北气候类型依次为温带大陆性半湿润季风气候、温带大陆性半干旱气候和高原气候等，夏季炎热多雨，冬季寒冷干燥。多年平均气温5.8~15.2 ℃，多年平均降水量为108~919 mm，年降水量由东南向西北递减。根据已有研究^［17］成果，将黄土高原划分为4个生态分区：黄土高塬沟壑区、黄土丘陵沟壑区、沙地和农灌区、土石山区及河谷平原区，分别占黄土高原总面积的35.03%、20.44%、20.77%、23.76%（图1）。其中，黄土高塬沟壑区是黄土高原面积最大、地貌破碎程度较高的区域，沟壑密布，水土流失显著，是区域水土保持与生态修复的重点；黄土丘陵沟壑区地形起伏较大，耕地分布广泛，土地利用敏感性强；沙地和农灌区以风沙地貌与农业灌溉区为主，生态系统脆弱；土石山区及河谷平原区地势相对低平，河流冲积平原发育，是黄土高原中耕作条件最好的区域之一。

1.2 数据来源

1）遥感数据。为便于与已有研究结果进行对比，本研究使用MOD13A1 NDVI（normalized difference vegetation index）数据产品，为2001—2022年500 m分辨率16 d数据，并利用最大值合成方法获取月值和年值。MOD17A3 NPP数据产品为2001—2022年500 m分辨率逐年数据，用于验证利用CASA模型估算的NPP的结果。以上数据均下载自（https：//search.earthdata.nasa.gov/）。

2）气象数据。包括气温、降水数据和太阳辐射数据，其中气温、降水数据为2001—2022年1 km分辨率逐月栅格数据（下载自国家青藏高原科学数据中心，https：//data.tpdc.ac.cn/home）；太阳辐射数据来源于ECMWF（欧洲中期天气预报中心）的ERA5再分析数据集，时间为2001—2022年，分辨率为0.1°×0.1°。

3）土地覆盖数据。土地覆盖数据采用CLCD（China Land Cover Dataset）数据集，空间分辨率为30 m，选取2001年和2022年的土地覆盖数据，下载自（https：//zenodo.org）。

4）DEM数据。来源于日本宇宙航空研究开发机构（JAXA）的高级陆地观测卫星，数据下载自（https：//search.asf.alaska.edu），空间分辨率为12.5 m。

以上数据统一投影为WGS1984坐标系，空间分辨率重采样至250 m，以便于后续数据的处理与分析。

1.3 研究方法

1.3.1 植被NEP估算模型

NEP为植被NPP与土壤异养呼吸（Rh）的差值，其计算公式^［18］为：

\mathrm{N}\mathrm{E}\mathrm{P}(x,t)=\mathrm{N}\mathrm{P}\mathrm{P}(x,t)-\mathrm{R}\mathrm{h}(x,t)（1）

式中：NEP（x，t）、NPP（x，t）和Rh（x，t）分别为在像元x处t年的植被净生态系统生产力、植被净初级生产力和土壤异养呼吸量，g/m²（以C计）^［18］。

1）植被净初级生产力（NPP）采用CASA模型进行估算，计算公式为：

\mathrm{N}\mathrm{P}\mathrm{P}(x,t)=\mathrm{A}\mathrm{P}\mathrm{A}\mathrm{R}(x,t)\times \epsilon (x,t)（2）

式中：APAR（x，t）为像元x在时间t内获取的光合有效辐射，MJ/m²；\epsilon（x，t）为像元x在时间t内实际光能利用率，g/MJ。

\mathrm{A}\mathrm{P}\mathrm{A}\mathrm{R}(x,t)=\mathrm{S}\mathrm{O}\mathrm{L}(x,t)\times \mathrm{F}\mathrm{P}\mathrm{A}\mathrm{R}(x,t)\times \mathrm{0.5}（3）

式中：SOL（x，t）为在像元x处时间t内的太阳辐射总量，MJ/m²；FPAR（x，t）为植被层对入射光合有效辐射的吸收比例；常数0.5为植被所能利用的太阳有效辐射占太阳总辐射的比例。

由于FPAR、NDVI与比值植被指数（simple ratio index， SR）存在线性关系，故可通过NDVI和SR计算FPAR。计算公式为：

\mathrm{F}\mathrm{P}\mathrm{A}\mathrm{R}(x,t)=\left(\frac{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}(x,t)-\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{i,\mathrm{m}\mathrm{i}\mathrm{n}}}{\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{i,\mathrm{m}\mathrm{a}\mathrm{x}}-\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{i,\mathrm{m}\mathrm{i}\mathrm{n}}}\right)\times (\mathrm{F}\mathrm{P}\mathrm{A}{\mathrm{R}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-\mathrm{F}\mathrm{P}\mathrm{A}{\mathrm{R}}_{\mathrm{m}\mathrm{i}\mathrm{n}})+\mathrm{F}\mathrm{P}\mathrm{A}{\mathrm{R}}_{\mathrm{m}\mathrm{i}\mathrm{n}}（4）

\mathrm{F}\mathrm{P}\mathrm{A}\mathrm{R}(x,t)=\left(\frac{\mathrm{S}\mathrm{R}(x,t)-\mathrm{S}{{{\mathrm{R}}_{\mathrm{i}}}_{,}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}{\mathrm{S}{{{\mathrm{R}}_{\mathrm{i}}}_{,}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-\mathrm{S}{{{\mathrm{R}}_{\mathrm{i}}}_{,}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}\right)\times (\mathrm{F}\mathrm{P}\mathrm{A}{\mathrm{R}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-\mathrm{F}\mathrm{P}\mathrm{A}{\mathrm{R}}_{\mathrm{m}\mathrm{i}\mathrm{n}})+\mathrm{F}\mathrm{P}\mathrm{A}{\mathrm{R}}_{\mathrm{m}\mathrm{i}\mathrm{n}}（5）

\mathrm{S}\mathrm{R}(x,t)=\left[\frac{\mathrm{1}+\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}(x,t)}{\mathrm{1}-\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}(x,t)}\right]（6）

\mathrm{F}\mathrm{P}\mathrm{A}\mathrm{R}(x,t)=a\times \mathrm{F}\mathrm{P}\mathrm{A}{\mathrm{R}}_{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}+(\mathrm{1}-a)\times \mathrm{F}\mathrm{P}\mathrm{A}{\mathrm{R}}_{\mathrm{S}\mathrm{R}}（7）

式中：NDVI _i，max、NDVI _i，min分别为第i种植被类型的年NDVI最大值和最小值；FPAR_min、FPAR_max的取值与植被类型无关，分别为0.001和0.95；SRi,\mathrm{m}\mathrm{i}\mathrm{n}、SRi,\mathrm{m}\mathrm{a}\mathrm{x}分别对应第i种植被类型NDVI的5%和95%百分位数。

光能利用率主要受气温和水分的影响，其计算公式为：

\epsilon (x,t)=T_{\mathrm{1}}(x,t)\times T_{\mathrm{2}}(x,t)\times W(x,t)\times {\epsilon }_{\mathrm{m}\mathrm{a}\mathrm{x}}（8）

式中：T₁（x，t）、T₂（x，t）分别为低温和高温对光能利用率的胁迫作用；W（x，t）为水分胁迫影响系数；ε_max为理想条件下的最大光能利用率。

2）土壤异养呼吸量采用PEI等^［12］提出的经验公式计算：

\mathrm{R}\mathrm{h}=\mathrm{0.22}\times [(\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{0.091}\mathrm{3}T(x,t)+\mathrm{l}\mathrm{n}(\mathrm{0.314}\mathrm{5}P(x,t)+\mathrm{1})]\times \mathrm{30}\times \mathrm{46.5}\mathrm{\%}（9）

式中：T（x，t）为像元x在第t月的月平均气温，℃；P（x，t）为像元x在第t月的月降水量，mm。

1.3.2 Sen趋势分析与MK检验

基于Theil-Sen趋势分析方法分析逐像元植被NEP变化趋势，并通过Mann-Kendall检验对变化趋势进行显著性检验。Theil-Sen计算公式为^［19-20］：

\beta =\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}\left(\frac{x_j-x_i}{j-i}\right)\hspace{1em}\forall j>i（10）

式中：x_i 为第i年的植被NEP值；\beta为趋势度，当\beta<0时，呈下降趋势，当\beta>0时，呈上升趋势，当\beta=0时，无显著变化。

Mann-Kendall检验公式为^［19-20］：

Q={\displaystyle \sum _{i=\mathrm{1}}^{n-\mathrm{1}}}{\displaystyle \sum _{j=i+\mathrm{1}}^n}\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(x_j-x_i)（11）

（12）

（13）

式中：Q为检验的统计量；Z为进行标准化过后的检验统计量；\mathrm{V}\mathrm{a}\mathrm{r}\left(Q\right)为方差，在给定的a显著性水平下，若Z<a/2则表示通过相应置信度的显著性检验。

通过Mann-Kendall统计检验法对变化趋势进行显著性检验，当MK检验统计量的绝对值>1.96时，表明变化趋势在0.05水平上显著（p<0.05）；当MK检验统计量的绝对值<1.96时，变化趋势不显著（p≥0.05）。将变化斜率\beta和MK检验结果进行叠加分析后，共得到4种变化趋势类型：显著增加、不显著增加、不显著减少和显著减少^［21］。

1.3.3 相关性分析

采用基于像元的空间分析法，分析植被NEP与气温、降水的相关性，确定气候因素与植被NEP的相关系数^［22］。

R_{xy}=\frac{{\displaystyle \sum _{i=\mathrm{1}}^n}\left((x_i-\overline{x})(y_i-\overline{y})\right)}{\sqrt[]{{\displaystyle \sum _{i=\mathrm{1}}^n}(x_i{-\overline{x})}^{\mathrm{2}}{\displaystyle \sum _{i=\mathrm{1}}^n}(y_i{-\overline{y})}^{\mathrm{2}}}}（14）

式中：R_xy 为x（2001—2022年平均气温和年降水量）和y（2001—2022年植被NEP）之间的相关系数，用于衡量二者之间相关关系的强度与方向，其取值为［-1，1］。当R>0时，表示二者呈正相关；当R_xy <0表示负相关；而R_xy 绝对值越大，表明变量间的相关性越紧密。x、y为2个变量的样本值；\overline{x}\mathrm{、}\overline{y}分别为变量x和y的样本平均值。根据查表（95%置信度，a=0.05时）得到的阈值为0.423，当|R_xy |>0.423，表示x与y之间相关性显著，其中R_xy >0为显著正相关，R_xy <0为显著负相关；当|R_xy |<0.423，表示x与y之间相关性不显著，其中R_xy >0为不显著正相关，R_xy <0为不显著负相关^［23］。

1.3.4 平均相对误差

平均相对误差（mean relative error， MRE）用于衡量估算值与参考值之间的相对偏差^［24］。在研究区内随机选取有效值点（同时具有本文基于CASA模型估算的NPP和MOD17A3 NPP数值的点），计算本文基于CASA模型估算的NPP与MOD17A3 NPP平均相对误差，计算公式为：

\mathrm{M}\mathrm{R}\mathrm{E}=\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}\frac{|{\widehat{y}}_i-y_i|}{y_i}\times \mathrm{100}\mathrm{\%}（15）

式中：y_i为第i个样本的参考值（MOD17A3 NPP数值）；{\widehat{y}}_i为第i个样本的估算值（本文基于CASA模型估算的NPP数值）；n为样本数量。MRE数值越小，表明基于CASA模型估算的NPP精度越高。

2 结果与分析

2.1 NEP估算结果验证

由于植被NEP模拟结果无法通过实测数据进行直接验证，可通过对植被NPP和土壤呼吸量的验证，间接验证植被NEP模拟结果的准确性，故采用MOD17A3 NPP数据产品验证CASA模型估算NPP^［25-27］的准确性，已有研究^{［13，27］}证明该方法具有一定的可靠性。在研究区内随机生成10 000个点，并最终筛选出6 000个有效点用于验证。结果表明，基于CASA模型估算的NPP与MOD17A3 NPP具有较高的相关性（R²=0.91，p<0.05）（图2），平均相对误差为22.33%，说明基于CASA模型估算的NPP精度较高。土壤异养呼吸量估算具有不确定性，本文采用PEI等^［12］提出的经验公式估算土壤异养呼吸量，此方法具有普遍性及可靠性。在已有黄土高原植被NEP计算结果中，周怡婷等^［3］计算的2000—2021年黄土高原植被NEP的多年平均值为104.62 g/m²，本文计算的2001—2022年多年平均值为327.06 g/m²，二者采用的计算方法基本一致，差异表现在前者使用的为500 m空间分辨率的MOD17A3 NPP，并最终采样为1 km后计算的多年平均值，而本文使用基于CASA模型计算的250 m分辨率的NPP，结合本文基于CASA模型计算的NPP结果高于MOD17A3 NPP数值，故本文计算的黄土高原植被NEP结果是基本可信的，可用于后续分析。

2.2 黄土高原植被NEP时空变化特征

2.2.1 时间变化过程

2001—2022年黄土高原植被NEP多年平均值为327.06 g/m²，2019年植被NEP平均值最大，为399.98 g/m²，2001年植被NEP平均值最小，为184.39 g/m²。黄土高塬沟壑区、黄土丘陵沟壑区、土石山区及河谷平原区植被NEP多年平均值高于黄土高原总体，分别为357.12、336.54、468.46 g/m²，沙地和农灌区植被NEP多年平均值低于黄土高原总体，为74.68 g/m²。以2011年为界，可将黄土高原及各分区植被NEP变化过程划分为2001—2011年的快速增加阶段（年平均增长率为4.94%）和2012—2022年的缓慢增加阶段（年平均增长率为0.24%）。变化趋势方面，2001—2022年黄土高原植被NEP呈显著线性增加，变化速率为7.51 g/（m²·a）（p<0.05）。不同分区中，黄土高塬沟壑区、黄土丘陵沟壑区、沙地和农灌区、土石山区及河谷平原区植被NEP均呈显著线性增加，变化速率分别为8.12、11.34、4.39、6.05 g/（m²·a） （p<0.05），均通过0.05显著性水平检验（图3a）。选取2022年黄土高原土地覆盖类型中占百分比较高的草地（占48.48%）、农业植被（占29.65%）和森林（占15.1%），分析其2001—2022年植被NEP变化过程表明，黄土高原不同土地覆盖类型植被NEP变化均呈线性增长趋势，变化速率为草地>森林>农业植被（图3b）。

不同分区中不同植被类型植被NEP对气候变化的响应规律在时间变化上皆呈波动上升趋势，并通过显著性检验，黄土高塬沟壑区和黄土丘陵沟壑区植被类型中草地植被NEP变化速率较森林大，分别为8.072、12.153 g/（m²·a），沙地和农灌区、土石山区及河谷平原区植被类型中森林植被NEP变化速率最大，分别为6.967、7.075 g/（m²·a）（图4）。

2.2.2 空间变化过程

2001—2022年黄土高原植被NEP以增加趋势为主，占区域总面积的95.54%，其中呈显著增加趋势的区域占77.38%，集中分布在黄土高原中部的黄土高塬沟壑区和黄土丘陵沟壑区；植被NEP呈显著减少趋势的地区仅占0.98%，主要分布在土石山区及河谷平原区的西部、沙地和农灌区的西南部及黄土高塬沟壑区西部（图5、图6）。黄土高塬沟壑区、黄土丘陵沟壑区、沙地和农灌区、土石山区及河谷平原区内植被NEP均以增加趋势为主，分别占区域总面积的96.67%、98.84%、94.26%、92.16%，其中呈显著增加的区域占各自分区总面积的81.18%、93.00%、71.91%、63.17%（图6）。

2.3 黄土高原植被NEP对气候变化的响应

2001—2022年黄土高原植被NEP与年平均气温以正相关为主，占区域总面积的67.33%，其中呈显著正相关的区域占6.28%，主要分布在黄土高原西南部及东北部，呈显著负相关的区域仅占1.38%，主要分布在黄土高原西部。按分区来看，黄土丘陵沟壑区、沙地和农灌区、土石山区及河谷平原区植被NEP与年平均气温均以正相关为主，显著正相关分别占各分区总面积的5.59%、2.26%、8.06%，主要分布在黄土丘陵沟壑区的中部和东部、沙地和农灌区的东部、土石山区及河谷平原区的东部区域；黄土高塬沟壑区植被NEP与年平均气温以负相关为主，显著负相关占分区总面积的2.72%，主要分布在黄土高塬沟壑区的西部及中部部分区域。各分区中植被NEP与年平均气温呈显著正相关或显著负相关的比例均较低。

2001—2022年黄土高原植被NEP与年降水量以正相关为主，占区域总面积的92.35%，其中呈显著正相关的区域占26.38%，主要分布在黄土高原西北部和东部，呈显著负相关的区域仅占0.13%，主要分布在黄土高原西部。按分区来看，黄土高塬沟壑区、黄土丘陵沟壑区、沙地和农灌区、土石山区及河谷平原区植被NEP与年降水量均以正相关为主，显著正相关分别占各分区总面积的21.00%、26.27%、43.24%、19.69%，主要分布在黄土高塬沟壑区的北部、黄土丘陵沟壑区东部、沙地和农灌区中部及东北部区域、土石山区及河谷平原区东北部区域。各分区中植被NEP与年降水量呈显著正相关且比例较高（图7）。

在黄土高原不同分区中，不同土地覆盖类型植被NEP与年平均气温和年降水量的相关性总体以正相关为主（黄土高塬沟壑区草地与年平均气温相关关系除外），且不同土地覆盖类型在不同分区中与年降水量的正相关（显著正相关）关系面积占比均高于与年平均气温（黄土高塬沟壑区森林与年降水量相关关系除外）（图8）。其中，黄土高塬沟壑区、黄土丘陵沟壑区、土石山区及河谷平原区不同土地覆盖类型植被NEP与年平均气温的显著正相关比例较高，而沙地和农灌区不同土地覆盖类型植被NEP与年降水量的显著正相关比例较高。

2.4 黄土高原植被NEP对退耕还林（草）的响应

2001—2022年黄土高原坡度>25°区域植被NEP多年平均值为527.64 g/m²，其中2019年平均值最大，为604.37 g/m²，2001年平均值最小，为367.20 g/m²。不同分区中，黄土高塬沟壑区、黄土丘陵沟壑区、沙地和农灌区坡度>25°区域植被NEP多年平均值低于黄土高原总体，分别为485.83、451.29、191.44 g/m²，土石山区及河谷平原区坡度>25°区域植被NEP多年平均值高于黄土高原总体，为689.37 g/m²。变化趋势方面，2001—2022年黄土高原坡度>25°区域植被NEP呈显著线性增加趋势，变化速率为7.93 g/（m²·a）（p<0.05）。不同分区中，黄土高塬沟壑区、黄土丘陵沟壑区、沙地和农灌区、土石山区及河谷平原区坡度>25°区域植被NEP均呈显著增加过程，变化速率分别为7.76、12.45、3.72、5.64 g/（m²·a）（p<0.05）（图9）。

从分区来看，土石山区及河谷平原区的植被NEP最高，平均值达689.37 g/m²；沙地和农灌区的植被NEP平均值最低，仅为191.44 g/m²。在变化速率方面，黄土丘陵沟壑区的植被NEP增速最大，为12.45 g/（m²·a），而沙地和农灌区的植被NEP增速最小，仅为3.72 g/（m²·a）。

为探讨退耕还林还草工程实施对植被NEP变化的影响，选取2001—2022年黄土高原坡度>25°区域中耕地转为林地、耕地转为草地和耕地保持不变区域，提取其植被NEP年平均值并绘图（图10）。2001—2022年黄土高原坡度>25°区域中，耕地转林地、耕地转草地和耕地保持不变区域的植被NEP均呈增加趋势，变化速率分别为13.659、12.419、10.100 g/（m²·a） （p<0.05）。

综上，黄土高原分区中，作为退耕还林还草工程实施的核心区域，黄土丘陵沟壑区植被NEP变化速率最大，且黄土高原坡度>25°区域中，耕地转林地和耕地转草地的区域植被NEP增速高于耕地保持不变区域，反映出退耕还林还草工程在促进黄土高原植被恢复和提高植被生态系统生产力方面成效显著。

3 讨 论

本文在利用CASA模型和土壤异养呼吸模型对2001—2022年黄土高原植被NEP进行估算过程中，采用空间分辨率为0.1°的太阳辐射数据，可能对估算结果产生部分影响。但通过与广泛使用的MOD17A3 NPP结果及已有黄土高原植被NEP估算结果^［3-5］对比可知，采用空间分辨率为0.1°的太阳辐射对本文估算的植被NEP结果影响较小，不影响估算结果的使用。

2001—2022年黄土高原植被NEP呈显著线性增加趋势，与全国植被NEP变化趋势一致^［28］。黄土高原植被NEP呈显著增加趋势的区域集中在黄土高塬沟壑区和黄土丘陵沟壑区中部，与张怡等^［5］研究结果一致，反映出退耕还林（草）等生态保护和建设工程实施有效提升黄土高原林草覆盖水平，改善生态系统结构和功能^［29］，促进植被固碳能力上升。不同分区中，植被NEP均呈增加趋势，其中黄土丘陵沟壑区增加速率最大，沙地和农灌区最小，与裴宏泽等^［4］和周怡婷等^［3］研究结果相近。2001—2022年黄土高原年平均气温和年降水量呈波动上升趋势（图11），为植被生长发育提供了有利的水热条件。黄土高原植被NEP与年平均气温和年降水量均以正相关为主，其中对降水的响应较为敏感。HE等^［30］、侯金龙等^［8］对我国华北等地区的研究也表明，降水是影响植被NEP空间格局变化的关键气候因子。但刘凤等^［18］研究表明，青海高原中气温对植被NEP以抑制作用为主，其影响高于降水；WANG等^［10］对我国西北干旱地区的研究发现，气温上升显著延长生长季，促进植物生产力的提高。可见，植被NEP与气温、降水的相互作用关系具有空间异质性。

已有研究^［21］表明，退耕还林（草）工程实施对陆地植被碳循环具有重要影响，在增加陆地植被碳汇和减缓气候变暖方面可发挥重要作用。本研究表明，2001—2022年黄土高原坡度>25°区域及其中耕地转为林地、耕地转为草地区域植被NEP均呈显著线性增加趋势，反映出黄土高原退耕还林（草）工程实施成效显著，改善植被覆盖水平，提升植被NEP。

4 结 论

1）2001—2022年黄土高原植被NEP多年平均值为327.06 g/m²，增加速率为7.51 g/（m²·a）（p<0.05）。不同分区中，植被NEP多年平均值为土石山区及河谷平原区>黄土高塬沟壑区>黄土丘陵沟壑区>沙地和农灌区。变化速率上，黄土丘陵沟壑区增加速率最快，为11.34 g/（m²·a） （p<0.05），沙地和农灌区增加速率最慢，为4.39 g/（m²·a） （p<0.05）。

2）空间上，2001—2022年黄土高原植被NEP以增加趋势为主，占区域总面积的95.54%，其中呈显著增加的区域占77.38%，主要分布在黄土高塬沟壑区中部和黄土丘陵沟壑区中部；呈显著减少的区域仅占0.98%，主要分布在土石山区及河谷平原区的西部、沙地和农灌区的西南部及黄土高塬沟壑区西部。不同分区中，黄土丘陵沟壑区植被NEP呈显著增加趋势的区域面积占比最高，占该分区面积的93.00%。

3）2001—2022年黄土高原植被NEP与年平均气温、年降水量均以正相关为主，其中呈显著正相关的区域分别占区域总面积的6.28%和26.38%，前者集中分布在黄土高塬沟壑区南部，后者集中分布在黄土高塬沟壑区的北部和东部、沙地和农灌区中部及东北部、土石山区及河谷平原区东北部；呈显著负相关的区域分别占1.38%和0.13%，分别集中在黄土高塬沟壑区西部和西北部。在黄土高原不同分区中，不同土地覆盖类型植被NEP与年平均气温和年降水量的相关关系总体以正相关为主（黄土高塬沟壑区草地与年平均气温相关关系除外）。且不同分区中农业植被、森林和草地植被NEP与年降水量呈正相关（显著正相关）的面积比例高于年平均气温，且沙地和农灌区植被NEP与年降水量关系最为明显。

4）2001—2022年黄土高原坡度>25°区域植被NEP多年平均值为527.64 g/m²，呈显著线性增加过程，变化速率为7.93 g/（m²·a）。不同分区中黄土丘陵沟壑区的植被NEP增速最大，为12.45 g/（m²·a），而沙地和农灌区的增速最小，仅为3.72 g/（m²·a）。坡度>25°区域中耕地转为林地和耕地转为草地的区域，植被NEP增加速率均高于耕地保持不变区域，表明退耕还林（草）工程实施对生态系统碳汇提升具有显著作用。本研究认为退耕还林（草）工程的实施及气候暖湿化对黄土高原植被NEP增加起到重要作用。受黄土高原气候暖湿化影响，后续退耕还林（草）等生态建设与恢复工程实施，可以考虑以自然恢复为主，人工恢复为辅的实施策略。

参考文献

[1] SHI L J， FENG P Y， WANG B， et al. Quantifying future drought change and associated uncertainty in southeastern Australia with multiple potential evapotranspiration models［J］.Journal of Hydrology，2020，590：e125394.

[2] 翟涌光，王晓妮，郝蕾，等.2001—2020年内蒙古净生态系统生产力格局多时间尺度分析［J］.生态环境学报，2024，33（2）：167-179.

[3] ZHAI Y G， WANG X N， HAO L， et al. Multi-time scale analysis of net ecosystem productivity pattern in Inner Mongolia from 2001 to 2020［J］.Ecology and Environmental Sciences，2024，33（2）：167-179.

[4] 周怡婷，严俊霞，刘菊，等.2000—2021年黄土高原生态分区NEP时空变化及其驱动因子［J］.环境科学，2024，45（5）：2806-2816.

[5] ZHOU Y T， YAN J X， LIU J， et al. Spatio-temporal variation in NEP in ecological zoning on the Loess Plateau and its driving factors from 2000 to 2021［J］.Environmental Science，2024，45（5）：2806-2816.

[6] 裴宏泽，赵亚超，张廷龙.2000—2020年黄土高原NEP时空格局与驱动力［J］.干旱区研究，2023，40（11）：1833-1844.

[7] PEI H Z， ZHAO Y C， ZHANG T L. Analysis of spatial and temporal patterns and drivers of local regional NEP in the Loess Plateau from 2000 to 2020［J］.Arid Zone Research，2023，40（11）：1833-1844.

[8] 张怡，王志慧，卢小平，等.黄土高原生态系统碳汇时空变化及其影响因素［J］.水土保持研究，2025，32（1）：266-274.

[9] ZHANG Y， WANG Z H， LU X P，et al. Evolution of ecosystem carbon sink and its driving factors in the Loess Plateau［J］.Research of Soil and Water Conservation，2025，32（1）：266-274.

[10] 刘文利，姜亮亮，刘冰，等.中国植被碳源/汇时空演变特征及其驱动因素［J］.生态学报，2024，44（4）：1456-1467.

[11] LIU W L， JIANG L L， LIU B， et al. Spatio-temporal evolution characteristics and driving factors analysis of vegetation carbon sources/sinks in China［J］.Acta Ecologica Sinica，2024，44（4）：1456-1467.

[12] ZHANG K， ZHU C M， MA X D， et al. Spatiotemporal variation characteristics and dynamic persistence analysis of carbon sources/sinks in the Yellow River basin［J］.Remote Sensing，2023，15（2）：e323.

[13] 侯金龙，马志强，杨澄，等.京津冀地区植被碳源/汇的时空变化特征及影响因素分析［J］.生态环境学报，2024，33（9）：1329-1338.

[14] HOU J L， MA Z Q， YANG C， et al. Analysis of spatio-temporal variation of vegetation carbon sources and sinks in the Beijing-Tianjin-Hebei region and influencing factors［J］.Ecology and Environmental Sciences，2024，33（9）：1329-1338.

[15] 管亚兵，王军，覃莉，等.黄土高原植被碳源/汇估算及其对土地利用变化的响应：以延河流域为例［J］.环境科学，2025，46（3）：1657-1665.

[16] GUAN Y B， WANG J， TAN Li， et al. Estimation of vegetation carbon source/sink and its response to land use change in the Loess Plateau： A case study of Yanhe River basin［J］.Environmental Science，2025，46（3）：1657-1665.

[17] WANG C， ZHAO W Z， ZHANG Y Y. The change in net ecosystem productivity and its driving mechanism in a mountain ecosystem of arid regions， northwest China［J］.Remote Sensing，2022，14（16）：e4046.

[18] GUO D， SONG X N， HU R H， et al. Grassland type-dependent spatiotemporal characteristics of productivity in Inner Mongolia and its response to climate factors［J］.Science of the Total Environment，2021，775：e145644.

[19] PEI Z Y， OUYANG H， ZHOU C P， et al. Carbon balance in an alpine steppe in the Qinghai-Tibet Plateau［J］.Journal of Integrative Plant Biology，2009，51（5）：521-526.

[20] 宋双双，秦世姣，孙彭成，等.基于多源遥感数据的汾河流域植被碳汇量估算及预测［J］.生态环境学报，2025，34（3）：345-357.

[21] SONG S S， QIN S J， SUN P C， et al. Estimation and prediction of vegetation carbon sinks in Fenhe River basin based on multi-source remote sensing data［J］.Ecology and Environmental Sciences，2025，34（3）：345-357.

[22] 干靓，朱佩露，杨颖.陆地生态系统碳汇能力估算及空间分布特征：以嘉兴西南三县为例［J］.生态学杂志，2025，44（4）：1393-1408.

[23] GAN J， ZHU P L， YANG Y. Estimation and spatial distribution of carbon sink capacity in terrestrial ecosystems： A case study of three counties in southwest Jiaxing［J］.Chinese Journal of Ecology，2025，44（4）：1393-1408.

[24] 周姝含，曹永强，么嘉棋，等.东北三省碳源/汇和碳盈亏时空分布与影响因素［J］.生态学报，2023，43（22）：9266-9280.

[25] ZHOU S H， CAO Y Q， YAO J Q， et al. Spatio-temporal distribution and influencing factors of carbon source/sink， carbon surplus and deficit in three northeast provinces［J］.Acta Ecologica Sinica，2023，43（22）：9266-9280.

[26] 周日平.黄土高原典型区土壤保持服务效应研究［J］.国土资源遥感，2019，31（2）：131-139.

[27] ZHOU R P. Assessing the soil erosion control service in the typical area of Loess Plateau［J］.Remote Sensing for Land and Resources，2019，31（2）：131-139.

[28] 杨梅焕，王添晴，李扬，等.黄土高原植被水分利用效率时空变化及其对不同影响因子的响应强度［J］.水土保持研究，2025，32（3）：159-169.

[29] YANG M H， WANG T Q， LI Y， et al. Spatiotemporal variation of vegetation water use efficiency on the Loess Plateau and its response intensity to different influencing factors［J］.Research of Soil and Water Conservation，2025，32（3）：159-169.

[30] 刘凤，曾永年.2000—2015年青海高原植被碳源/汇时空格局及变化［J］.生态学报，2021，41（14）：5792-5803.

[31] LIU F， ZENG Y N. Analysis of the spatio-temporal variation of vegetation carbon source/sink in Qinghai Plateau from 2000—2015［J］.Acta Ecologica Sinica，2021，41（14）：5792-5803.

[32] SHAO Y K， ZHU Q， FENG Z K， et al. Temporal and spatial assessment of carbon flux dynamics： Evaluating emissions and sequestration in the three northern protection forest project areas supported by google earth engine［J］.Remote Sensing，2024，16（5）：e777.

[33] 黄琪，彭立，李赛男，等.基于GAM的喀斯特植被覆盖与驱动因素非线性关系分析［J］.中国环境科学，2023，43（5）：2489-2496.

[34] HUANG Q， PENG L， LI S N， et al. Analysis of the nonlinear relationship between karst vegetation cover and driving factors based on GAM［J］.China Environmental Science，2023，43（5）：2489-2496.

[35] 赖金林，齐实，廖瑞恩，等.2000—2019年西南高山峡谷区植被变化对气候变化和人类活动的响应［J］.农业工程学报，2023，39（14）：155-163.

[36] LAI J L， QI S， LIAO R E， et al. Vegetation change responses to climate change and human activities in southwest alpine canyon areas of China from 2000 to 2019［J］.Transactions of the Chinese Society of Agricultural Engineering，2023，39（14）：155-163.

[37] 高萌萌，李小磊，杨楠，等.黄河流域植被时空变化及其与土壤湿度的相关性分析［J］.水文地质工程地质，2023，50（3）：172-181.

[38] GAO M M， LI X L， YANG N， et al. Spatio-temporal variation of vegetation and its correlation with soil moisture in the Yellow River basin［J］.Hydrogeology and Engineering Geology，2023，50（3）：172-181.

[39] 谢艳玲，夏正清，王涛，等.黄河流域植被NPP时空变化及其对水热条件和退耕还林还草工程实施的响应［J］.测绘通报，2023（2）：15-20.

[40] XIE Y L， XIA Z Q， WANG T， et al. Temporal and spatial variation of vegetation net primary product and its response to hydrothermal conditions and grain for green project in the Yellow River basin［J］.Bulletin of Surveying and Mapping，2023（2）：15-20.

[41] 张雪纯，马毅，张靖宇，等.基于分段自适应算法的浅海水深遥感反演融合模型研究［J］.海洋科学，2020，44（6）：1-11.

[42] ZHANG X C， MA Y， ZHANG J Y， et al. Research on the remote sensing inversion fusion model of shallow water depth based on the piecewise adaptive algorithm［J］.Marine Sciences，2020，44（6）：1-11.

[43] 徐勇，黄雯婷，郭振东，等.2000—2020年我国西南地区植被NEP时空变化及其驱动因素的相对贡献［J］.环境科学研究，2023，36（3）：557-570.

[44] XU Y， HUANG W T， GUO Z D， et al. Spatio-temporal variation of vegetation net ecosystem productivity and relative contribution of driving forces in southwest China from 2000 to 2020［J］.Research of Environmental Sciences，2023，36（3）：557-570.

[45] LI Z， CHEN Y N， ZHANG Q F， et al. Spatial patterns of vegetation carbon sinks and sources under water constraint in Central Asia［J］.Journal of Hydrology，2020，590：e125355.

[46] 王江涛，杨永崇，杨梅焕.基于地理探测器的黄土高原NPP时空变化及驱动力研究［J］.西安理工大学学报，2023，39（1）：12-20.

[47] WANG J T， YANG Y C， YANG M H. Spatial and temporal variation and driving forces of NPP on the Loess Plateau based on Geodetector［J］.Journal of Xi´an University of Technology，2023，39（1）：12-20.

[48] HUANG Y T， WANG F， ZHENG L J， et al. Corrigendum： Changes and net ecosystem productivity of terrestrial ecosystems and their influencing factors in China from 2000 to 2019［J］.Frontiers in Plant Science，2023，14：e259137.

[49] 白萌，莫淑红，莫兴国，等.退耕还林背景下黄土高原蒸散量时空演变特征及归因［J］.生态学报，2023，43（20）：8344-8358.

[50] BAI M， MO S H， MO X G， et al. Spatio-temporal variation of evapotranspiration and its attribution over the Loess Plateau since the implementation of the Grain for Green Project［J］.Acta Ecologica Sinica，2023，43（20）：8344-8358.

[51] HE H L， WANG S Q， ZHANG L， et al. Altered trends in carbon uptake in China´s terrestrial ecosystems under the enhanced summer monsoon and warming hiatus［J］.National Science Review，2019，6（3）：505-514.

基金资助

国家自然科学基金项目(42271309)

2025年度陕西省哲学社会科学研究专项(2025YB0239)