青海省植被NEP时空变化及驱动因素分析

冶晓娟， 王永辉， 潘红忠， 白 钰， 董得福， 姚华明,3

（1.长江大学油气地球化学与环境湖北省重点实验室，湖北 武汉 430000；2.黄河上游水电有限责任开发公司，青海 西宁 810000；3.中国长江电力股份有限公司，湖北 宜昌 443000）

陆地生态系统是全球气候变化研究的核心内容之一，在全球碳平衡中发挥着重要作用［1］。定量估算陆地碳汇时空变化是预测气候变化和执行温室气体控制协议的基础。通过比较区域间的碳平衡才能深刻理解区域间循环过程。通过人工生态工程的实施来提高植物固碳量，是区域内增汇的主要手段。净生态系统生产力（NEP）是单位时间单位面积上植物的碳收支状况，可作为定性和定量化分析陆地植被碳循环的关键参量［2］，是近几年的研究热点，如基于气候变量的统计模型，基于光能利用率相关的参数模型和生态系统机制的过程模型等，尤其是3S 技术推动遥感反演获取相关地表参数和植被信息发展，使得3S技术支撑下快速发展的生态系统机制模型成为了当前研究的主流之一［3］。

青海省位于“第三极”东北部，被誉为“三江之源，中华水塔”，是对全球气候变化响应最敏感的地区之一，故青海省植被固碳能力变化对青藏高原乃至全球控制温室气体排放和碳循环均有巨大影响［4］。在青海省城市化和工业化进程中，人口分布、地表覆盖也随之发生改变，这必将在一定程度上影响陆地上的植被碳循环［5-6］。虽然，已确定陆地生态系统将在未来持续长时间发挥碳汇的作用，但植被NEP 驱动机制和植被固碳量变化仍存在较大争议，如何在复杂的变化环境中阐明气候敏感区植被NEP近几十年的变化特征以及驱动机制有待探究［7-8］。刘凤等［9］和周夏飞等［10］对青海高原植被、青藏高原草地的NEP进行气候驱动力分析，指出气候变化对该区域NEP的演变产生影响；张新中等［11］通过分析石羊河流域的NEP与NDVI、降水、气温和太阳辐射的关系，指出流域内NEP 变化主要受降水和气温影响，受太阳辐射的影响不显著；刘春雨［12］研究了西南地区植被NEP 的时空演变，并结合人类强度指数，采用地理加权回归法探析不同视角下的人为因素对植被NEP 影响。已有对植被NEP 的驱动因素研究主要是基于气候因素，对人为因素的驱动研究还限于单因子的影响。为探析双因素交互作用对植被NEP的影响，本文将地理探测器用于NEP驱动因素的定量分析，该方法既可进行单因子影响分析，亦可探测多种人为因子和气候因子交互作用对NEP的影响［13-15］。

近年来为应对西北地区“气候转型”［16］和人为活动的影响，研究环境变化下的植被固碳能力成为实现“碳中和”的迫切需求［17］。本文利用经验公式和MODIS数据估算青海省NEP，并探析青海省NEP时空演化规律及驱动因子，以期科学评估青海省植被碳储量变化和对青海省植被增汇政策提供理论参考。

1 研究区概况与数据

1.1 研究区概况

青海省位于89°35′～103°04′E，31°9′～39°19′N之间，属高原温带半干旱地区，多风，日照时间长，降雨少，年降水量在50～550 mm，年均温在-5.1～9.0 ℃；地势总体西北向东南逐渐倾斜；区域降水稀少，主要集中于春、夏季，径流主要来源于冰雪融水和西南暖湿气流。

青海省林地资源主要分布于江河源区附近，全省林地面积2.65×106km2［18］，天然草地面积4.2×105km2。草地类型主要为高寒草甸和草原（图1），主要分布在祁连山东北部和柴达木盆地东部地区。由于特殊的气候、高海拔和风化侵蚀地貌，导致青海省植被易受人类和气候变化影响［6］。

图1 研究区植被概况Fig.1 Vegetation map of study area

1.2 数据来源与处理

净初级生产力(NPP)数据来源于MODIS(https://ladsweb.modaps.eosdis.nasa.gov/)所提供的MOD17A3-HGF产品，时间分辨率为年，空间分辨率500 m，时间跨度2000—2020 年。归一化植被指数（NDVI）数据来源于MOD13Q1 产品，时间分辨率16 d，空间分辨率250 m，时间跨度2001—2020年，采用最大值合成月数据。土地利用数据来源于MCD12Q1 产品，时间分辨率为年，空间分辨率500 m，时间跨度2001—2019 年。以上数据均利用MRT 及ArcGIS 进行数据拼接，并投影裁剪研究区范围。

青海省气象数据来自2000—2020 年国家气象科学数据中心(http://data.cma.cn/)月值数据集,空间分辨率为0.5°×0.5°，对该数据用ArcGIS进行数据还原并采用样条函数法插值。

人均GDP 和人口密度数据均来自中国科学院资源环境科学与数据中心(https://www.resdc.cn)，分辨率为1 km。以上所有数据均采用Albers 等面积圆锥投影进行统一分析。

用地理探测器探测空间分异性时，选取40 km的空间网格来分析变化环境对植被净生态系统生产力（NEP）的影响，降水、气温采用自然间断法划分为6类；人均GDP、人口密度采用分位数划分为6类［19］。

2 模型与方法

2.1 NEP的估算

NEP 常被用来描述区域碳平衡动态变化，可指示单位时间单位面积上植物的碳收支状况，不仅可作为衡量植被活动的关键因子，还能用来表征植被的健康程度和固碳能力［20］。在不考虑其他自然及人为因素影响下，可用NEP值表征陆地和大气生态系统的碳交换率。NEP＞0，说明植被固碳量多于土壤呼吸（RH）排放的碳量，表现为植被碳汇［21］。计算公式如下：

本文土壤呼吸RH采用裴志永等［22］建立的经验公式，由气温和降水间接估算。经已有研究实地采样验证，该公式可用于监测青海地区、西北地区土壤呼吸。计算公式如下［12,23］：

式中：T指年均温(℃)；P指年降水量(mm)；RH为年土壤呼吸（g C·m-2）。

2.2 驱动因素分析

2.2.1 趋势分析 采用Mann-Kendall 法［24］及Theil-Sen法分析2000—2020年植被NEP的年际变化规律与动态变化［25］，计算公式如下：

式中：xk，xi为连续的NEP数据序列；n为年份。标准化检验统计量ZC，其中ZC＞0 时表示上升趋势，反之，则表示下降趋势。由于Mann-Kendall 法强调单调趋势，结果不一定是线性，为体现变化幅度引入Theil-Sen斜率法表征变化幅度Q的大小［26］：

式中：1＜i＜k＜n。当Q＞0时，反映上升趋势；反之，则反映下降趋势。

2.2.2 相关性分析

（1）偏相关系数。运用偏相关分析可定量刻画植被NEP 对其他影响要素的响应。偏相关系数计算公式如下：

假设有i、j、k3 个变量，Pij,k是自变量k与自变量i和j的偏相关系数；Pij，Pjk，Pik为两因子间的相关系数［27］。

（2）地理探测器。地理探测器是一种探测空间分异性，并揭示其背后驱动力的统计学方法［28］。根据研究需要，本文应用因子探测器，交互探测器等4种探测器探测NEP的驱动力及多因子交互作用［15］。因子X对属性Y的空间分异解释程度用q值度量，q值的计算公式如下：

式中：q为某因子对NEP 空间分异的解释力，q的值域［0，1］，q值越大表征该因子的空间分异性愈明显。在极端情况下，q=1 表明该因子完全控制了NEP的空间分布。h取值范围为［1，L］，表征因子X的分区，σh2和σ2分别为第h层和全区NEP 的方差，Nh和N分别为层h和全区的单元数；SSW和SST分别为层内方差之和和全区总方差。

交互探测器可用来评估不同因子间交互作用对NEP分异解释力的强弱。评估方法是先分别计算2种影响因子X1和X2对NEP的解释力q，再计算其交互作用时的q值，最后进行比较。判断依据如表1。

表1 交互作用判断依据Tab.1 Basis of interaction judgment

3 结果与分析

3.1 NEP精度评价

由于NEP的实测数据难以获取，本文采用其他学者在青海省实测NEP结果、青藏高原数据中心三江源气象站点实测NPP 结果和站点气象数据模拟RH结果计算NEP［29-33］。将实测数据与估算结果进行相关性分析、显著性检验，表明两者的相关系数为0.864（图2），实测与估算数据显著相关（R2=0.7821，P＜0.01），显示结果具有可靠性。

图2 青海省NEP模拟结果验证Fig.2 Verification of NEP simulation results in Qinghai Province

3.2 青海省NEP时空分布特征

3.2.1 青海省NEP空间分布特征 近20 a青海省年均NEP的空间分布差异显著，呈现出东南高西北低的特征（图3）。青海省年均植被NEP为134.05 g C·m-2·a-1，其中NEP＞0 区域高达4.89×105km2（占比74.15%），NEP＜0 区域可达890 km2。NEP 高值区（312～735.90 g C·m-2·a-1）面积总占比5.07%，主要位于祁连山附近、青海湖湿地附近、江河源区和黄河上游谷地等地。植被碳源区主要有格尔木市、治多县和杂多县等地。植被在空间分布上过渡明显的原因，主要与当地的降水状况以及植被分布有关。NPP 多年均值分布与NEP 分布相似。RH呈西北向东南递增之势，具体分布为森林＞草原＞荒漠，空间差异较明显。

图3 青海省NEP、RH和NPP空间分布Fig.3 Spatial distribution of NEP,RH and NPP in Qinghai Province

3.2.2 青海省NEP的年际变化特征 近20 a青海省NEP 年际变化较明显（图4），2000—2020 年青海省生态系统的固碳能力呈波动增加趋势。年均植被NEP 在114.55～155.49 g C·m-2·a-1，NEP 年均增幅为1.54 g C·m-2·a-1。植被NEP 的极大值、极小值分别出现在2018年、2001年。按照植被NEP近20 a的变化态势，划分3个时段讨论其年际变化特征：2000—2008年间植被NEP波动增加，2006年出现第1个峰值144.93 g C·m-2·a-1；2008—2015年间，NEP围绕多年均值上下波动；2015—2020 年间，植被NEP 值逐年上升，该阶段也为增幅最大的时段。

图4 青海省NEP的年际变化Fig.4 Interannual variation of NEP in Qinghai Province

3.2.3 青海省NEP的动态变化 为研究青海省NEP的空间动态变化规律，采用逐像元Theil-Sen斜率估计法，得出青海省21 a NEP 变化趋势（图5）。植被NEP 空间分布增长格局与其空间分布格局相似，从变化趋势可知，青海省植被固碳能力整体上有增长趋势，少部分呈负增长趋势的区域主要分布在海西市和玉树市内，约占青海省面积的2.17%，主要为植被稀疏地区；增长幅度大的高值区（占比3.17%）主要位于青海湖湿地、贵南县和共和县附近等地，这些区域的降水较多，植被覆盖丰富。将Theil-Sen 趋势分析和M-K显著性检验结果叠加，通过P＜0.05显著性水平的区域面积为3.55×105km2，其中显著增加区域占全省面积的52.95%，主要分布在柴达木盆地荒漠、祁连山森林与高寒草原、江河源区-甘南高寒草甸草原等生态区；显著减少的区域仅占0.23%，零星分布于天峻县、玛多县、曲麻莱县、湟中区、治多县等地。

图5 青海省Mann-Kendall和Theil-Sen植被NEP的空间分布Fig.5 Spatial distribution of vegetation NEP of Mann-Kendall and Theil-Sen in Qinghai Province

3.3 青海省NEP驱动机制分析

3.3.1 NEP对气候因子的响应 气候因子通过影响植被的生长态势、净初级生产力以及土壤呼吸，进而影响植被净生态系统生产力的变化；尤其是气温、降水这2个主要气象因素，在一定程度会影响植被NEP的变化。

由图6 可知，近20 a 青海省植被NEP 的极大值出现在2010 年，降水量极大值出现在2009 年（463.54 mm），在2010年左右年均温也有1个显著的峰值。说明青海省在2010 年前后水热组合相对较好，有利于提高生态系统固碳能力。NEP 的极小值出现在2001 年，该年降水量显著偏低，气温虽相比2000年略有升高，但长历时干旱易造成土壤呼吸增强和植被退化，从而导致固碳能力下降。对上述3个因素间的响应程度进行分析，2000—2014年间植被NEP变化相对降水存在一定的滞后性，滞后期约为1 a；2014 年后两者变化相对一致。近20 a，气温与植被NEP的变化保持一致。上述分析表明，青海省植被NEP对年均温和年降水的响应程度较高。

图6 2000—2020年青海省NEP与年均温、年降水量的年际变化Fig.6 Interannual variations of NEP,annual temperature and annual precipitation in Qinghai Province from 2000 to 2020

空间上分析青海省植被NEP与年均温、年降水量的偏相关关系，并对其进行P＜0.05的显著性检验（图7、图8）。区域内植被NEP与降水的偏相关系数均值为0.134，与气温的偏相关系数均值为0.429，这表明气温、降水与NEP 均呈正相关，且降水的偏相关小于气温的偏相关。初步分析干旱和半干旱区气温对植被的分布及长势影响更大。

2000—2020 年，青海省降水介于325～460.85 mm，年平均降水量为395.26 mm。降水与植被NEP的偏相关系数介于-0.24～0.42 之间（图7a）。植被NEP与降水呈负相关的区域占全省面积的23.32%，主要分布在玉树、囊谦、曲麻莱、甘德、久治、班玛、杂多、河南县等地。植被NEP与年降水量呈正相关的区域占青海省面积的65.75%，主要分布在海南市、海北市、海西、海东等地；其中与降水呈现显著正相关（P＜0.05）的区域面积达7.81×104km2，约占全省面积的11.7%（图8a）。

图7 NEP与降水和年均温偏相关关系Fig.7 Correlation of NEP with precipitation and annual mean temperature deviation

2000—2020 年，青海省年均温-2.35～-0.70 ℃，年均温均值-1.44 ℃。植被NEP与气温呈负相关的区域面积为1.142×104km2，仅占总面积的1.71%，零星分布于各市；植被NEP与年均温呈正相关的区域占青海省面积的71.5%，这些区域位于山麓处、江河源区以及河流谷地，主要为森林、草原。其中，通过显著性检验（P＜0.05）的植被NEP与气温正相关的区域面积达2.572×105km2，占全省面积的38.53%（图8b）。

图8 NEP与降水和年均温偏相关显著性分析Fig.8 Correlation analysis of NEP with precipitation and annual mean temperature deviation

3.3.2 NEP与NDVI的相关性 分析表明，青海省年均植被NEP 与NDVI 的空间分布具有高度相似性，在青海湖湿地、贵南县和共和县等地是植被NEP（312～735.90 g C·m-2·a-1）与NDVI（0.68～0.91）的高值分布区。在空间分布上（图9），植被NEP与NDVI呈正相关的区域占全省面积的65.83%，其中通过显著性检验（P＜0.05）的区域占比达31.43%，这些区域主要分布在海南、海东和海西东南部。NEP与NDVI呈负相关的区域占全省面积的7.37%，主要分布在曲麻莱、玉树、杂多、久治、同仁和门源等地；通过P＜0.05 显著性检验的区域仅占0.12%，这些区域零星分布于各地。

图9 植被NEP与植被NDVI相关性系数分布Fig.9 Distribution of correlation coefficient between vegetation NEP and NDVI

3.3.3 NEP的人为驱动力 土地利用类型转变能揭示人为因素对植被NEP的影响。参考《土地利用现状分类》和原数据IGBP标准，将MCD12Q1重新划分为6大类（耕地、林地、草地、水域、城乡居民用地、未利用地）。为研究2000—2020年土地利用类型对青海省NEP 的影响，结合数据的时间跨度（2005 年、2010 年、2015 年、2020 年），以5 a 为1 个研究期，探究人为因素对植被NEP的影响。分析表明，2015—2020年的土地利用类型转换最为明显，转换面积达研究区总面积的2.7%，该阶段主要是未利用地和草地之间的类型转换，主体为荒漠向草地转换，其面积占本阶段转换总面积的49.48%，主要位于柴达木盆地荒漠区、昆仑山东段高寒荒漠草原、北羌塘高寒荒漠草原区，主要为荒漠与草原过渡区域。耕地面积在2005—2015年呈下降趋势并转换为草地，主要是政府在该时期内采取了退耕还草等生态保护措施。2001—2020 年荒漠与草地间的转换较为频繁，2010—2015年荒漠化面积略大于绿洲化面积,该阶段存在微弱的荒漠化趋势，植被NEP呈下降趋势（图10）。

图10 青海省土地利用转移桑基图Fig.10 Sangji map of land use transfer in Qinghai Province

3.3.4 因子影响力探测 为探究不同因子对青海省植被NEP 的影响，选取6 种典型因子并对其进行地理探测器分析，探测其空间分异的驱动力（其中2020年人均GDP数据缺失）。单因子对植被NEP的影响分析结果表明（图11），在该研究期内，NDVI是影响植被NEP 的主导因子，其q值均在0.7 以上，并且年际波动较小，说明植被覆盖度对植被NEP贡献较大。气温与降水的解释力较强，均超过了0.3，降水在2015年的影响力超过气温，主要是由于年降水偏少，限制了植被生长，从而抑制了植被的固碳能力。土地利用和人均GDP 影响力相对稳定且解释力较低，人口密度解释力较高，在2010 年的影响力高于气温；降水在2020年解释力略高于土地利用。

图11 代表年份植被NEP单因子影响力和交互作用q值Fig.11 Single factor influence and interaction q value of vegetation NEP in representative years

地理探测器交互作用结果表明（图11）：区域内两因子交互作用的q值显著大于单因子的q值，表明交互作用呈非线性增强，两因子交互作用后对植被的NEP解释力会显著增强。对4个代表年份的主导交互因子进行排列，最大主导作用的交互因子是该年最强主导因子与其他因子的组合。2000—2020年最强主导交互因子均为NDVI 与人均GDP、气温与NDVI，解释力均达到了0.8以上。虽然GDP的单因子影响力小，但当人均GDP与NDVI交互作用时，在特定年份会超越其他因素。

4 讨论

MOD17A3HGF 模型已被用于估算青海省NPP,通过作物产量验证，证实估算结果可靠［18］。估算RH的检验模型在青海高原和青藏高原通过了实地验证［9,18,22］。

本文估算青海省近20 a NEP 年均值为134.69 g C·m-2·a-1，与周夏飞等［10］（120.8 g C·m-2·a-1）、刘凤等［9］（128.40 g C·m-2·a-1）的估算结果较为接近。虽然本文采用的NEP估算方法不同且更新了时间序列，但估算结果与前人的较为接近，说明估算方法可行。该方法存在的不确定性，主要源自土壤呼吸的不确定性。将RH定义为与气象因子相关的函数，容易忽视呼吸底物的分解对异养呼吸的贡献，容易造成高植被覆盖度区域的碳汇高估，低植被覆盖度区域对碳汇的低估，因此需慎重采用［34］。此外，RH究竟是否适用于青海省全省植被NEP 的估算需要大量的地表观测实测值作参考，尽管以多年平均值做比较，但时间跨度的区别难免产生一定误差［23］。本文NEP与气温的关系，与已有研究存在一定差异［9］，主要由于两者对NPP 的估算方式及采用的气象原数据插值方法存在差异所致，在中亚地区和石羊河流域的研究结论与本文相似［35-36］。

诸因子对植被NEP的影响分析表明，影响较大的是以气温、降水为代表的气候因子和以土地利用为代表的人为因子。分析结果显示，气温对NEP的影响程度要大于降水，原因是适度的升温对NPP的促进作用要大于RH的促进作用［17］，而且当海拔超过2000 m 时，水分对于植物NPP 的限制作用越来越弱［24］。与气温呈偏负相关的主要是荒漠区域（如柴达木盆地荒漠区），该区域地表缺水、蒸发旺盛，气温降低会导致区域内蒸发减弱，有利于缺水状况缓解并促进植被NEP的积累，进而提升该地区的植被固碳能力。受三江源头太阳辐射强烈、冰川资源丰富、海拔高等因素影响，在青海南部NEP 与降水呈负相关、与气温呈正相关。NDVI 与植被NEP 在空间分布上具有极大相似性，植被覆盖率高的区域光合效率高，植物群落固定在地面的有机碳含量及植被净碳交换量同样也高。

运用地理探测器探测青海省NEP 的人为及气候驱动力，整体上气候对植被NEP地理分异性的解释力要超过人为因素，但在某些特殊年份，情况迥异。在气候趋于“暖湿”背景下［16］，由于气温逐年上升，陆地生态系统的植被储碳能力对气候的响应将呈现新的格局，人为因素对NEP空间分异解释力较弱，但与其他因子交互时均超过0.3。双因子的协同交互作用远超单因子的作用，未来在提升生态系统固碳能力时，需在关注当前气候变化的同时，采取适当的人为干预。降水对NEP 的解释力在2020 年最低的原因尚待进一步探究。

由于生长季内农作物固定的碳将迅速进入大气碳循环，鉴于研究区内农田比例较大，若能剔除作物NEP 的影响将对区域增汇更有指导意义［37］。青海省作为五大牧区之一，在后续研究中还应将牲畜饲养量纳入人为因素，分析其对青海植被NEP的影响。

5 结论

本文详细分析了青海省近20 a NEP 的时空变化规律及驱动因素，并利用地理探测器定量分析不同因子对NEP的影响，具体结论如下：

（1）时间尺度上2000—2020 年青海省NEP 呈增加趋势，增幅为1.07 g C·m-2·a-1，多年平均NEP为134.69 g C·m-2·a-1。空间上，NEP 呈自东向西递减的特征，NEP分布与植被覆盖基本一致。

（2）在空间与年际尺度上，植被NEP与气温、降水均呈正相关，气温与NEP的相关性大于降水。

（3）地理探测器表明，近20 a 植被覆盖度对NEP 起主导作用，气温的影响作用相对较强，在干旱年份降水的影响力要超过除NDVI以外的其他因子，双因子的交互作用超过单因子影响力。青海省植被NEP 空间分布是气候和人为因素共同作用的结果，且不同影响因素对其作用的强弱不同，在特定年份人类活动与其他因素交互作用是植被NEP变化的主导因素。