2000-2019年青海地区植被总初级生产力遥感估算及时空变化分析

林小丁，常乐，冯丹

（伊犁师范大学化学与环境科学学院，新疆伊犁835000）

植被总初级生产力（gross primary productivity，GPP）是指单位时间内生态系统中绿色植物通过光合作用，吸收太阳能同化二氧化碳固定的有机碳总量［1］。GPP决定了进入陆地生态系统的初始能量与物质总量，反映了陆地生态系统植被生产力，对陆地生态系统碳循环和全球气候变化研究具有深远意义。大约40%的陆地生态系统植被生产力被人类直接或间接利用［2］，因此对GPP的估算及其变化研究是人类赖以生存的重要基础，也是评价生态系统可持续发展的重要指标。对其进行准确的模拟、时空分布和动态变化分析及其与气候因子相关性分析对陆地生态系统和地球可持续发展有很好的评估意义。目前，经验模型是估算区域和全球GPP的有效途径［3］，即利用地面实测数据与遥感植被指数构建经验关系得出的遥感植被指数GPP模型。成熟的GPP模型有增强植被指数（enhanced vegetation index，EVI）［4-11］、归一化植被指数（normalized difference vegetation index，NDVI）［7］、光化学植被指数（photochemical reflectance index，PRI）［12］、日光诱导植被荧光（solar-induced fluorescence，SIF）［13］以及植被近红外反射率指数NIRv（near-infrared reflectance vegetation，NIRv）［14］。其中，NIRv值在一定程度上考虑到像元内植被比例和土壤反射率等因素的影响，与GPP具有高度一致性［15］。

青海是青藏高原重要的组成部分，也是我国降水量跨度最大的省份，其生态环境对全球气候变化较为敏感。同时，青海作为三江源所在地，监测其生态系统变化对我国的生态文明建设具有战略意义。近年来，许多研究对青海地区的净初级生产力（net primary productivity，NPP）进行了遥感估算。刘鹏［16］利用GLOPEM-CEVSA模型对2000-2012年三江源地区NPP进行了估算，并分析了其对气候因素的响应；张继平等［17］分析了三江源地区草地生态系统碳汇量及碳储量时空变化；刘凤等［18］利用CASA改进模型分析了近16年青海地区植被NPP空间变化及气候与人为因素响应；贾俊鹤等［19］利用CASA模型分析了中国西北部NPP时空变化及气候响应。然而，NPP是植物光合作用固定的碳（GPP）与自养呼吸作用消耗的碳之差，其中自养呼吸作用又包含生长呼吸和维持呼吸，影响因素众多，且与植被指数关系较弱，因此NPP估算具有较大不确定性。相比之下，使用植被指数进行GPP估算则更为精确，因此也更能反映出植被生长状况的时空变化。然而，目前只有少数研究估算GPP，比如，丁少文等［20］通过基于植被NIRv模型与NDVI指数对比进行了GPP估算，分析了2001-2016蒙山地区植被动态变化；吴力博等［21］对青海三江源地区的栽培草地GPP进行了简要分析。因此，基于新型植被指数NIRv的青海地区总初级生产力遥感估算和进一步长时间序列时空变化分析及其对气象因素响应分析深化，都对总体把握青海地区生产力潜力评估具有重要研究意义。

为此，本研究采用2000-2019年青海地区MODIS遥感数据，结合土地利用资料与气象资料，利用NIRv模型，估算了20年青海地区植被GPP，并分析了其时空分布规律，及其对气候因素变化的响应。以期为青海地区生态系统碳循环和评价当地生态系统可持续发展能力做贡献，并为大尺度范围内应用新型植被指数NIRv可行性提供重要方向。

1 材料与方法

1.1 研究区概况

青海地区位于我国西北内陆，河流湖泊分布众多，其东部地区为青藏高原向黄土高原过渡地带，西部为高原和盆地，属于高原大陆性气候。地理位置介于北纬31°36′-39°19′，东经89°35′-103°04′，全省面积72.23万km2，地形复杂，地貌多样（图1a）。年均温分布北高南低，全省年降水量15～750 mm，总分布趋势由东南向西北逐渐减少［22］。其中，中覆盖草地、低覆盖草地为青海地区主要植被类型，广泛分布在南部地区，占全省面积的48.9%；高覆盖草地在东北部地区、南部黄河流域和中部长江流域广泛分布；灌木林地主要集中于东部地区、三江源地区和格尔木市周围，占全省面积的2.3%；疏林地在东部地区、南部澜沧江流域和格尔木市周围稀疏分布；林地仅在西宁东部地区少量分布；沙地、裸岩石砾地和盐碱地则主要覆盖在西北地区（图1b）。

图1 青海地区位置图（a）及土地利用图（b）Fig.1 Location map（a）and land use map（b）of Qinghai Province

1.2 数据与处理

采用2000-2019年的h25v05和h26v05两景MOD13A 1 16 d合成500 m分辨率植被指数产 品（https：//e4ftl01.cr.usgs.gov/MOLT/MOD13A 1.006/）和BESS全球逐日0.05度分辨率光合有效辐射（photosynthetically active radiation，PAR）产品（http：//environment.snu.ac.kr/bess_rad/）。首先从MOD13A 1产品计算月平均NIRv。从原始文件中提取出NDVI和NIR数据集并计算16 d合成NIRv（式1），再提取出合成日期（composite day of the year）数据集，该日期指示了每个像元NIRv所对应的具体日期，以此计算月平均NIRv。

其次，采用Jiang等［23］的方法对NIRv进行土壤背景校正。假设土壤背景NIRv值只在0～0.1变化，将0～0.1分为10个等间距（0.01）的区间，对每个像元的时间序列NIRv数据统计落在各区间内的频次，频次最高的区间（众数）作为土壤背景NIRv的初始估计。鉴于稀疏植被的NIRv值也可能落在0～0.1，对每个像元只统计其小于平均NIRv的值，以此去除生长季对土壤背景NIRv估算的干扰。考虑众数区间及其相邻的两个区间，计算落在这3个区间的NIRv的平均值得到土壤背景NIRv（NIRv，soil）。同时，从每个像元的时间序列NIRv数据中获取最大NIRv（NIRv，max），则土壤背景校正NIRv（SANIRv）可计算为：

另一方面，从BESS产品计算月平均PAR［24］，并使用地理空间栅格数据处理库GDAL进行裁减、投影和重采样，使其空间参考与MOD13A 1一致。

采用青海地区3个站点的地面观测GPP数据作为基准。3个站点分别是海北灌丛生态系统观测站（101°18.77′E，37°36.77′N，海拔3216 m），海北湿地生态系统观测站（101°19.64′E，37°36.51′N，海拔3357 m）和阿柔高寒草原观测站（100°27′51.6″E，38°02′50.3″N，海拔3033 m）。海北灌丛2003-2010年日均GPP数据来源于中国通量观测研究联盟（http：//www.chinaflux.org）。海北灌丛2002年和海北湿地2003-2005年月均GPP来源于FLUXNET 2015数据集（https：//fluxnet.fluxdata.org/data/fluxnet2015-dataset）。阿柔站2015-2017年日均净生态系统碳交换量（net ecosystem exchange，NEE）和气象数据来源于国家青藏高原科学数据中心（https：//data.tpdc.ac.cn），然后用REddy ProcWeb网上工具对日均NEE进行插补和拆分得到日均GPP。最后，所有站点GPP都换算到月均GPP。

采用青海地区1∶10万土地利用数据集［25］、中国1 km分辨率逐月平均气温数据集（1901-2017）［26］和中国1 km分辨率逐月降水量数据集（1901-2017）［27］，采用偏相关分析法进行GPP时空变化分析。数据均来源于青藏高原科学数据中心（https：//data.tpdc.ac.cn）。矢量格式的土地利用数据在ArcGIS中进行区域合并并栅格化，然后使用GDAL统计每个MODIS像元内出现频次最大的地类作为该像元的地类。气温和降水数据使用GDAL进行裁减、投影和重采样，使其空间参考与MOD13A 1一致。

此外，还采用2000-2019年的h25v05和h26v05两景MOD17A 2H 8 d合成500 m分辨率GPP产品作为对比。该产品与MOD13A 1空间参考一致。研究计算其月平均值并用地面观测GPP对其进行验证分析，并与基于SANIRv计算的GPP进行精度比较。

1.3 模型与方法

1.3.1 GPP估算模型 根据Wu等［28］和Jiang等［23］的研究，表明NIRv与GPP高度线性相关，且经过土壤背景校正后的SANIRv则进一步与GPP呈比例关系：

式中：c为比例系数。使用遥感观测SANIRv数据和地面观测GPP数据对式（3）进行标定，得到研究区最优的c值。进而使用式（3）对地面观测GPP数据进行升尺度，得到覆盖研究区的2000-2019年月平均500 m分辨率GPP数据。

1.3.2 GPP模拟精度评价 采用100次重复5折交叉验证对式（3）进行精度评价。在每次交叉验证中，将成对的遥感观测SANIRv和地面观测GPP数据随机分为5份，取其中的4份作为训练数据计算c值，剩下1分独立于模型标定的数据作为验证数据。采用确定系数（R 2）和均方根误差（RMSE）作为精度评价指标。100次试验的平均R2和RMSE作为基于SANIRv的GPP估算精度。也通过比较MOD17A 2H GPP和地面观测GPP，计算相应的R2和RMSE作为MODIS GPP产品的估算精度作为参照［28］。

1.3.3 时空分析方法GPP年际趋势分析使用一元线性回归法研究青海地区植被GPP的年际变化规律，通过计算每个像元的趋势线斜率（slope）来反映GPP年际变化规律。当slope＞0时，表示在研究时段中GPP增加，反之减少。研究采用Mann-Kendall方法对趋势分析进行显著性检验，因显著水平（α）具有任意性且考虑到本研究站点精确数据数量较少，将显著水平（α）小于0.1视为变化显著［29］。

在GPP与气候因素相关性分析中，主要通过分别建立气温、降水与GPP的偏相关系数来分析气候因素与GPP变化的关系。采用t检验法对各相关系数、纯相关系数进行显著性检验［30］（将其显著性水平α取值为10%）。对研究区域根据其相关系数不同进行分类，把r＞0，P＜0.1、r＞0，P≥0.1、r＜0，P≤0.1、r≤0，P≥0.1分为显著正相关、不显著正相关、显著负相关和不显著负相关的不同区域。

2 结果与分析

2.1 GPP结果估算精度评价

通过对青海地区2000-2019年海北湿地、海北灌丛、阿柔3个地面站点实测数据与本研究采用土壤背景校正NIRv模型估算的GPP值进行相关性分析及交叉验证。本研究NIRv模型估算GPP与地面实测GPP数据呈线性关系，确定系数R2=0.91（P＜0.001），均方根误差（RMSE）为3.0 g C·m-2·d-1（图2a）。PAR×SANIRv估算GPP交叉验证确定系数R2=0.91（图2b），两者表明NIRv模型估算结果稳健，具有可靠性，能代表青海地区总初级生产力并向多地区推广。此外，本研究还将NIRv模型估算GPP值与MODIS GPP进行对比分析（图2c），MODIS GPP在本研究区内有明显低估情况，且在GPP较低值区域表现不好。综上所述，本研究NIRv模型估算GPP在青海地区的应用结果从站点至区域均具有较高精度，且具有适用性。

图2 GPP估算的精度评价Fig.2 Accuracy evaluation of GPP estimation

2.2 青海地区植被GPP空间分布格局

青海地区2000-2019年平均GPP整体由西向东呈显著增长趋势，空间内存在明显差异，主要与植被分布和地理气候特征差异相关。如图3所示，西北部柴达木盆地区域年均GPP在0～450 g C·m-2·yr-1，该区域植被类型主要为GPP水平低下的沙地、裸岩石砾地及盐碱地，但区域边缘处分布有GPP水平较高的疏林地、高覆盖草地和少量灌木林，因此年均GPP差异显著；南部三江源区域年均GPP主要在0～660.1 g C·m-2·yr-1，该区是我国最大的自然保护区，植被类型以低覆盖草地为主，但生态环境脆弱，出现了严重的草场退化、土地沙化和土壤侵蚀现象，因而草地覆盖度低，多年GPP低下。同时，随该区海拔由西向东逐渐降低，二氧化碳浓度增高且温度适宜，植被类型逐渐多样，主要为疏林地、灌木林和中覆盖草地，且该区实施了大面积退耕还林生态工程，因而GPP得到显著提升；东部地区植被类型丰富，灌木林地和草地类型分布面积较大，沼泽湿地和森林资源丰富，水热条件优渥，总初级生产力水平高，大致在150.0～980.7 g C·m-2·yr-1，是青海地区高植被总初级生产力的主要贡献区域。

图3 青海地区多年平均GPPFig.3 Aver age GPP in Qinghai Pr ovince for many year s

2.3 青海地区植被GPP时间变化特征

青海地区2000-2019年GPP值在116.3～165.9 Tg C·yr-1，平均值为140.5 Tg C·yr-1，单位面积GPP为195.1 g C·m-2·yr-1。区域内GPP总量的年际变化如图4所示，近20年青海地区GPP总体呈波动增长趋势，年增长率为1.25 Tg C·yr-1（P＜0.05），说明近20年中青海地区植被情况得到改善。其中，2001年GPP均值最低，仅为116.3 Tg C·yr-1，2010年GPP均值最高，为165.9 Tg C·yr-1。并且该趋势具有明显阶段性变化，与中覆盖草地和低覆盖草地GPP变化趋势大致契合。2003-2006年和2008-2010年GPP呈快速上升趋势，主要因为中覆盖草地和低覆盖草地GPP的显著提高与该年间降水量、年均温显著增加相关；2006-2008年和2013-2014年显著降低，主要与中覆盖草地和低覆盖草地GPP下降相关（图5）。

图4 青海地区GPP总量的年际变化Fig.4 The interannual change of GPP in Qinghai Province

如图5所示，不同植被类型的GPP值年际变化表现出很大差异。其中，中覆盖草地是青海地区总初级生产力最高的植被类型，其总初级生产力在2010-2011年达到峰值，为52.4 Tg C·yr-1。其次是低覆盖草地类型，总初级生产力在2010年达到峰值，为51.5 Tg C·yr-1。高覆盖草地和灌木林GPP趋势波动较小，在11.9～16.0 Tg C·yr-1和8.5～10.7 Tg C·yr-1。

图5 青海地区不同植被类型GPP的年际变化Fig.5 Inter annual var iation of GPP of different vegetation types in Qinghai Pr ovince

2.4 青海地区GPP的空间变化特征

本研究通过一元线性回归对青海地区近20年年均GPP进行变化分析（图6），GPP由西向东呈增加趋势，全省有54%的区域呈显著上升趋势。其中，东部地区因灌木林地、疏林地和高覆盖草地分布广泛，使得年均GPP增长率较高，约为0.007273～0.063460 g C·m-2·yr-1；西北地区多以沙地、戈壁、盐碱地为主，边缘分布有灌木林地和疏林地，使得年均GPP增加，增长率约为0.004932～0.048760 g C·m-2·yr-1；西南地区，植被类型主要以低覆盖草地和中覆盖草地为主，海拔较高，常年有积雪覆盖，因受全球气候变暖影响，年均GPP略有增加，增长率约0.0002636～0.0232900 g C·m-2·yr-1；东南地区，植被类型丰富，主要以中覆盖草地和高覆盖草地为主，由于该植被类型受水汽因素影响较大且该区水土流失较严重，因此年均GPP变化较大，约为-0.049530～0.000663 g C·m-2·yr-1，呈现出显著下降趋势。

图6 青海地区GPP的变化趋势（a）与显著性检验（b）Fig.6 The change trend（a）and significance test（b）of GPP in Qinghai Province

2.5 青海地区GPP与气候因素的相关性

植被生长以及积累GPP能力受到气象因子影响，本研究通过分析气温和降水时空变化及分布，与近20年青海地区植被GPP建立关系，采用逐像元偏相关分析法进一步探讨。首先，研究区近20年年均温度和年降水量呈上升趋势（图7）。其中，青海地区研究期内温度年平均值为-2.5℃，最大值为9.3℃，最小值为-19.6℃；总降水量年平均值为283.6 mm，最大值为650.2 mm，最小值为14.0 mm；年均温年际变化在0～0.05℃·yr-1，东、南部地区变化较明显；总降水量年际变化在2.214～6.686 mm·yr-1，除南部少数地区外整体呈上升趋势（图8a～d）。其次，温度显著性上升区域约占总面积的48%，主要为东、南部地区，不显著上升区域约占总面积的51%，多为西北、东部边缘地区（图8e），年均降水量变化除极少部分地区不显著下降外，大部分区域降水量为不显著上升（图8f）。

图7 青海地区年均温度和年降水量的变化趋势Fig.7 The interannual change of annual mean temper ature and annual total pr ecipitation in Qinghai Pr ovince

图8 青海地区温度与降水的时空分异Fig.8 Spatiotemporal variation of temperature and precipitation in Qinghai Province

逐像元偏相关分析结果如图9所示，表明在研究区尺度上温度和降水与GPP的相关系数（去趋势后）分别为0.45（P＜0.1）和0.09（P＞0.1），植被GPP与温度呈正相关性，与降水量相关性较微弱。且在不同区域内，不同植被对气温和降水因子响应具有显著空间异质性。其中，植被GPP在主要植被类型为中覆盖草地地区与温度呈现出显著正相关，相关系数为0.6，原因为该区域本身水汽条件充足，该植被类型对气温变化较为敏感，因此气温的增加使得中覆盖草地在适宜温度下生长加速，GPP不断积累；在主要植被类型为高覆盖草地和灌木林地地区表现为不显著正相关，该植被类型稳定性较强，不易受到气温因子的干扰。植被GPP与降水的相关性由东北向西南递减，且分化现象严重。南部低覆盖草地区域大部分呈不显著负相关，相关系数约为-0.29，一是因为青海地区南部海拔较高，植被蒸腾耗水量较少，低覆盖草地积累GPP主要受温度限制。二是由于长江、澜沧江、黄河均位于青海地区南部，河流众多，本身水汽条件充沛，且该区平均年总降水量为青海全区域最高，降水量超过低覆盖草地生长所需水分最大值，反而会抑制低覆盖草地生长；只有东北部高覆盖草地和灌木林地地区呈现出部分显著正相关，相关系数约为0.55，原因是东北部地区海拔较低，水分主导植被GPP年际变化，随着降水量的不断增加，高覆盖草地和灌木林地生长发育得到促进，有利于GPP的不断累积。

图9 GPP同温度与降水的相关性Fig.9 The correlation of GPP with temperature and precipitation

3 讨论

本研究采用最新的土壤背景校正NIRv模型，估算了青海地区植被GPP，结果显示，研究期内年均GPP处于增长趋势，温度升高为主要因素，这与刘刚等［31］、许洁等［32］、韩炳宏等［33］、刘青瑞［34］等众多研究结果一致。即全球陆地生态系统总体呈“变绿”趋势［35］，我国陆地生态系统也以变绿为主［36］，主要原因为CO2浓度增加及温度升高［37］。目前，全球GPP遥感产品主要采用过程模型［38］，机器学习模型［39］和光能利用率模型［40］。过程模型较为严谨，但结构复杂，参数众多，运算速度慢。机器学习模型最为精确，但在缺乏训练数据的区域表现较差。光能利用率模型应用最为广泛（如MODIS GPP产品），但其预设了GPP对温度和湿度的函数，精度较差。本研究采用的土壤背景校正NIRv模型虽然是经验模型，却能够精确模拟美国几十个站点的GPP［23］，且仅有一个参数，具有可靠性。本研究证明该模型在青海地区同样具有很好的表现（R2=0.91），是一种有前景的方法。此外，本研究估算的青海地区年均GPP为140.5 Tg C·yr-1，考虑到GPP与NPP的比值通常在0.4～0.6，则青海地区年均NPP约为56.2～84.3 Tg C·yr-1，与刘凤等［18］通过对青海地区植被NPP时空变化的研究结果（79.4 Tg C·yr-1）相近。本研究与刘凤等［18］的研究采用了不同的数据源与方法，却得到一致的结果，反映了结果的可信度。

本研究分析青海大部分地区年均植被GPP与年均气温存在显著性不高的正相关，及大部分地区植被GPP与年降水量相关关系微弱，这可能存在两方面原因，一是时间序列不够长，二是气象数据的插值精度影响了与植被GPP的相关性。因此，通过优化气象插值精度的方法将是提高与植被GPP相关性的有效手段。其次，虽然温度和降水是影响GPP的主要气象因子，但研究区属高原大陆性半干旱气候，日照充足，光合有效辐射对GPP也起主导作用，本研究在与气候因素关系分析中，仅分析了气温和降水，未考虑光合有效辐射因子对GPP的影响。原因有：1）由于研究区样本量的限制，较精确的气象数据仅有气温和降水。其他一些因素，如冰冻期、大气湿度、土壤水分、二氧化碳、云覆盖度等也会对GPP的年际变化产生影响，但缺乏数据，从而表现出一定的局限性。2）由于是模型输入，故没有直接纳入分析。因此，本研究虽采用最新优化的土壤背景校正NIRv模型提高了区域植被GPP的估算精度，但植被GPP变化由自然和人为等多种因素共同作用产生，因此在今后研究中提高植被指数分辨率，充分考虑太阳辐射及多种相关因子的响应则是未来研究方向。

4 结论

本研究采用土壤背景校正NIRv模型估算青海植被2000-2019年GPP，并对结果进行估算精度评价。在此基础上分析青海地区GPP时空变化特征，不同植被类型年际变化GPP和与同期气象因素的响应关系。初步揭示不同植被对地区GPP的影响和对气候的响应机制，并取得5点结论：

1）土壤背景校正NIRv模型估算的青海地区（海北灌丛生态系统观测站、海北湿地生态系统观测站、阿柔高寒草原观测站）GPP与地面实测GPP数据呈线性关系，相关系数R2=0.91（P＜0.001），模型估算均方根误差（RMSE）为3.0 g C·m-2·d-1。相较于MODIS GPP，土壤背景校正NIRv模型估算GPP在青海地区的应用结果从站点至区域均具有较高精度，更具有适用性。

2）从空间分布角度看，青海地区植被GPP多年平均值为140.5 Tg C·yr-1，整体趋势由西向东呈显著增加。空间内存在明显差异性，源于不同植被类型分布、地理气候特征差异较大。

3）从时间跨度角度看，2000-2019年青海地区年均GPP整体处于上升趋势，年增长率为1.25 Tg C·yr-1（P＜0.05）。

4）青海地区GPP时空变化具有明显阶段性变化特征，不同植被类型GPP的年际变化也表现出较大差异，其中，中覆盖草地和低覆盖草地总初级生产力与青海地区GPP变化趋势大致契合。

5）气象因素对不同植被GPP的影响存在明显空间异质性，整体上，年均气温与GPP变化的相关性高于平均降水量。其中，气温因素是以中覆盖草地为主的南部地区的GPP变化的重要因素，而降水因素是柴达木盆地东部和龙羊峡水库周围地区GPP变化的重要因素。