羌塘国家自然保护区植被生长季NDVI时空变化及其对气候变化的响应

袁 雷，周刊社，张东东

（西藏自治区气候中心，拉萨 850000）

0 引言

IPCC第六次评估报告第一组工作报告指出，2011—2020年全球地表温度比工业革命时期上升了1.09℃。从未来20年的平均温度变化预估来看，全球温升预计将达到或超过1.5℃。在全球变暖的背景下，陆地生态系统正在受到气候变化的影响，作为陆地生态系统的主要组分，植被对气候变化十分敏感[1]。归一化植被指数(normalized differential vergetation index，NDVI)是反映植被生长状况的最佳指标之一[2]，广泛用于植被生长状况演变分析[3-4]，植被NDVI在不同时空尺度下的变化特征及与气候变化的关系也有很多研究[5-7]。青藏高原属于气候变化的敏感区和生态脆弱带，被认为是研究陆地生态系统对气候变化响应机制的理想场所[8]。羌塘国家级自然保护区位于青藏高原腹地，由于其独特的高寒荒漠生态系统，生态更为脆弱和敏感[9]，已经有很多学者对羌塘的植被生态[9-11]、荒漠化[12-13]、土壤温湿度[14-15]、降水[16]、持续干旱日数[17]、冰雹日数[18]、积雪[19]、土壤冻结天数[20]、地表湿润度[21]等进行了研究。其中，吴晓萍等[9]利用2001—2010年的NDVI数据及气象资料对羌塘国家级自然保护区东、中、西部各县NDVI、气温、降水变化进行了分析，认为保护区东部NDVI呈显著上升趋势，中部和西部呈下降趋势。东部年均NDVI与气温呈显著正相关、与降水呈显著负相关，中部反之；西部年均NDVI与气温、降水均呈显著负相关。杜军等[11]分析了羌塘气温、降水、水体以及NDVI变化趋势，认为气温、降水和NDVI都呈增加趋势。这些研究以区域平均NDVI分析植被变化，很难体现小区域的植被NDVI较详细的变化特征，对于不同功能区植被NDVI时空变化特征及其差异没有进行定量化研究，也少有羌塘植被NDVI对气候变化响应的研究。鉴于此，笔者以MODIS传感器的NDVI数据作为基础，以羌塘国家级自然保护区各功能区为研究对象，逐像元定量计算其NDVI时空变化趋势及其稳定性，期望能够为羌塘国家级自然保护区生态环境保护提供科学参考。

1 数据与方法

1.1 研究区概况

羌塘国家级自然保护区（以下简称保护区）是仅次于格陵兰国家公园的世界第二大陆地自然保护区，位于西藏自治区西北部，平均海拔5000 m以上，被称为“世界屋脊的屋脊”，保护区内植被种类较少、群落结构简单，从东南向西北草地植被大体呈高寒草甸草原、高寒草原、高寒荒漠草原、高寒荒漠的分布[11]。高寒植被以紫花针茅(Stipa purpurea)、藏羊茅(Festuca wallichanica)和高山蒿草(Kobresia pygmaea)等为主[10]，有实验区、缓冲区和核心区3个大的功能分区（图1）。

图1 羌塘国家级自然保护区地理位置

保护区年平均气温大都在0℃以下，年平均地温为-5.0～-1.5℃，降水量为95.6～294.9 mm，光照条件充足，年日照时数为2800～3600 h。杜军[11]研究表明1971—2017年自然保护区附近气象站点年平均气温以0.46℃/10 a的速率显著升高，明显高于同期全球和亚洲地表温度的升温率。

1.2 数据来源与处理

NDVI数据为2000—2020年的MODIS NDVI月数据，空间分辨率为0.01°，来自于国家气象中心，已经进行了大气校正、辐射校正和几何校正。在研究植被变化时，使用最大值合成法(maximum value composites，MVC)计算植被生长季（6—9月）的NDVI数据。气象数据使用中国第一代大气和陆面再分析产品(CRA/LAND)的2000—2020年月气温和月总降水量数据，该数据来自于国家气象信息中心，使用该数据计算出植被生长季平均气温和植被生长季总降水量。由于使用的CRA/LAND数据空间分辨率为0.25°，因此在分析NDVI与气象因子相关性时，将植被生长季NDVI数据重采样为0.25°。植被类型数据来源于中科院1:100万中国植被数据集。

1.3 分析方法

1.3.1 Sen趋势分析 使用Sen趋势分析来逐像元描述植被生长季NDVI及气候因子的长期变化趋势[22-23]。Sen趋势分析又被称为Sen斜率估计，计算时间序列斜率对的中值，抗噪声能力强，是一种稳健的非参数统计的趋势计算方法。该方法计算效率高，对于测量误差和离群数据不敏感，常被用于长时间序列数据的趋势分析中。Sen斜率计算如式(1)所示。

式中，xj和xi为时间序列数据，使用趋势度β来判断时间序列趋势的升降，当β＞0时时间序列呈上升的趋势，当β＜0时时间序列呈下降的趋势。

1.3.2 Mann-Kendall趋势检验 使用Mann-Kendall(MK)检验法对趋势分析结果进行检验[22-23]，趋势检验法过程如下：对于时间序列Xt=( )x1,x2,...,xn，做如下假设：H0—假设序列中的数据为独立同分布样本，即无显著趋势，H1—假设序列存在上升或下降单调趋势。检验统计量S由式(2)～(3)计算。

本文中时间序列长度为21（2000—2020年），采用Sen趋势逐像元计算保护区植被生长季NDVI、植被生长季气温和植被生长季降水量变化趋势；使用MK检验法逐像元计算保护区植被生长季NDVI、植被生长季气温和植被生长季降水量变化趋势显著性。最终结果分为5级[4]，即极显著下降(slope＜0，P≤0.01)、显著下降(slope＜0，0.01＜P≤0.05)、变化不显著(P＞0.05)、显著上升(slope＞0，0.01＜P≤0.05)、极显著上升(slope＜P≤0，P≤0.01)。

1.3.3 变化趋势率 使用最小二乘法分析研究区及各功能区平均NDVI的变化趋势率，其计算如式(5)[24]。

式中，n为研究序列长度，在本研究中为21；i表示第i年，xi表示第i年的NDVI，θ表示变化趋势率。

1.3.4 稳定性 变异系数CV即标准差与平均值的比值[24]，反映变异程度，广泛运用于波动水平的分析中，计算如式(6)。

式中，SDNDVI为逐年植被生长季NDVI的标准差，为逐年植被生长季NDVI的平均值。CV值消除了单位和平均值不同对2个或多个变量变异程度比较的影响。本研究采用CV分析保护区植被生长季NDVI逐个像元21年间的变异情况，以揭示植被生长季NDVI的稳定性。并按几何间隔法将稳定性分为5类[4]，即高稳 定 (0.017＜CV≤0.055)、较 高 稳 定 (0.055＜CV≤0.058)、中 等 稳 定 (0.058＜CV≤0.096)、较 低 稳 定(0.096＜CV≤0.692)、低稳定(0.692＜CV≤9.954)。

1.3.5 相关性检验 采用皮尔逊(Pearson)相关分析方法研究气象要素与植被生长季NDVI之间的关系，使用T检验来检验相关系数显著性。相关系数如式(7)所示[24]。

式中，rxy相关系数，n为年份，x为自变量（气温、降水等气象要素），y为因变量（植被生长季NDVI）。

我小学五年级因抗战爆发、家乡沦陷、学校停办而失学，一直在老家农村种地。实际上即使我在上学的时候，也是一直跟着大人下地劳动的，农村的孩子，一般十来岁早就下地劳动了。

2 结果与分析

2.1 保护区植被生长季NDVI空间分布特征

从整个保护区来看，植被生长季NDVI空间分布差异较大（图2），表现为由东南向西北逐渐减小的空间分布特点。2000—2020年植被生长季NDVI的平均值为0.138，根据植被生长季NDVI值大小将保护区格点像元划分为5个等级[25]，分别为≤0.1、0.1～0.2、0.2～0.3、0.3～0.4、＞0.4。植被生长季NDVI≤0.1为无植被覆盖区，植被生长季NDVI处于0.1～0.2区间为低植被覆盖区，植被生长季NDVI处于0.2～0.4为中等植被覆盖区，植被生长季NDVI大于0.4为高等植被覆盖区[25]，无植被覆盖区占保护区面积的22.97%，主要分布于北缓冲区昆仑山南麓，植被群系为垫状驼绒藜荒漠；低植被覆盖区占保护区面积的64.39%，植被群系为青藏薹草草原、紫花针茅草原和部分沙生针茅草原；中等植被覆盖区占保护区面积的12.56%，主要分布于实验区东部和南缓冲区中东部，植被群系为紫花针茅草原、三指雪莲花和西藏扁芒菊稀疏植被；高等植被覆盖区只占保护区面积的0.08%，零星分布于保护区东部。

图2 生长季NDVI空间分布图

在保护区的3个功能区中，低植被覆盖区面积占比均达到58.50%以上，高植被覆盖区面积占比均在0.12%以下，表明3个功能区整体植被条件均较差；实验区、缓冲区和核心区的中等植被覆盖区域面积占比分别为21.98%、11.21%和8.30%，表明实验区的植被条件在保护区中最好。3个功能区植被生长季NDVI均是东部好于西部。

2.2 趋势分析

2.2.1 趋势变化空间特征 从近21年保护区植被生长季NDVI变化趋势空间分布来看（图3a），保护区86.39%区域植被生长季NDVI呈缓慢的上升趋势，变化趋势率为(0～0.01)/10 a的区域占保护区面积的57.51%，在保护区广泛分布；变化趋势率为(0.01～0.03)/10 a的区域占保护区面积的27.06%，主要分布于保护区西北部。从变化趋势显著性检验来看（图3b），植被生长季NDVI上升达到显著性水平(P＜0.05)的区域占保护区面积的45.15%，主要分布在保护区中部和北部，增加趋势多在(0～0.03)/10 a。植被生长季NDVI下降达到显著性水平(P＜0.05)的区域占保护区面积的1.11%，零星分布在缓冲区南部。总体来看，保护区植被状况明显改善。结合图2～3可以看出，植被生长季NDVI变化趋势率基本上表现为：植被生长季NDVI值小的区域，变化上升趋势率小；植被生长季NDVI值大的区域，变化上升趋势率大；植被生长季NDVI值处于0.1～0.2范围时保护区西北部变化趋势率高于东南部。实验区、缓冲区和核心区植被生长季NDVI上升达到显著性水平（显著上升和极显著上升）的面积占比分别是24.47%、47.43%、57.27%，表明核心区植被生长季NDVI改善程度优于实验区和缓冲区。

图3 生长季NDVI变化趋势(a)及显著性检验(b)

2.2.2 NDVI年际趋势变化特征 从保护区植被生长季各等级NDVI面积比例逐年变化分析（图4），植被生长季NDVI值处在0.1～0.2范围的面积比例在每年均最高，变化区间在57.7%～67.8%，总体表现出明显的增加趋势，平均每年增加0.33个百分点(P＜0.001)；植被生长季NDVI≤0.1的面积比例呈现明显的减少趋势，平均每年减少0.56个百分点(P＜0.001)，到2020年只占保护区面积的21.0%，在2018年甚至只占到保护区面积比例的17.7%；对于植被生长季NDVI值在0.2～0.3范围的面积比例呈现震荡上升趋势，在2017年和2018年达到保护区植被面积的16.7%和16.3%，分析表明保护区无植被覆盖区减少，低植被覆盖区和中等植被覆盖区增加，植被总体呈现变好趋势。

图4 保护区生长季NDVI各等级面积比例

2.3 时空稳定性分析

使用CV来分析2000—2020年植被生长季NDVI的年际间稳定性程度，保护区植被生长季NDVI时空稳定性见图5。较低稳定度分布面积最大，占整个保护区植被面积的47.86%；中等稳定度分布面积也较大，占整个保护区植被面积的38.8%；较低稳定度和中等稳定度总共占整个保护区面积的86.66%，表明2000年以来保护区的NDVI稳定性较差，这是可能由于保护区植被以草地植被和灌丛为主，受气象条件影响明显，气象条件的较小波动就能引起植被状态的较大变动。较高稳定度和高稳定度主要分布于保护区中东部及西部边沿。保护区中东部和南部气象条件较好，利于植被生长，因此植被生长季NDVI稳定性较高；保护区西部植被生长季NDVI稳定性较高的原因可能是该区域雪山融水较丰富，利于植被生长。

结合图2和图5可以看出，植被生长季NDVI≤0.1的大部分区域处于较低稳定度，表明该区域易受环境影响，对外界干扰反应敏感，这与王青霞[25]、徐维新[26]的研究结果较为一致。

图5 生长季NDVI稳定性空间分布

在各功能区中，实验区、缓冲区和核心区高稳定度区域面积占比分别是11.47%、9.48%和9.49%，中等稳定度区域面积占比分别是40.44%、37.65%和39.47%，较低稳定度区域面积占比分别是44.82%、49.39%和47.61%，实验区高稳定度和中等稳定度区域面积占比均略高于缓冲区和核心区，而较低稳定度面积占比略低于缓冲区和核心区，说明实验区植被稳定度好于缓冲区和核心区。原因可能是实验区在保护区南部和东南部，气象条件较好，植被稳定度高。

2.4 相关性分析

2000—2020年，保护区植被生长季平均气温变化率为 0.01～0.87℃/10 a（图 6a），平均值是 0.43℃/10 a，增温趋势明显，东南角至西北角的条带增温最明显，每10年增加0.5～0.6℃，该结论与杜军[11]的研究结果基本一致；植被生长季累积降水量变化率为-3.55～27.52 mm/10 a（图6b），平均值为9.92 mm/10 a，保护区中部降水量呈减少趋势、东部与西部降水呈增加趋势，其中，东南部达到每10年增加17～27 mm，西部大部分区域也达到了每10年增加10～17 mm。总体来看，保护区植被生长季气候趋向于“暖湿化”。

图6 保护区2000-2020年生长季气温和降水量变化趋势

逐像元分别计算生长季NDVI与生长季平均气温、生长季总降水量的相关系数。结果显示空间分布差异较为明显，保护区大部分区域生长季NDVI与气温呈正相关性（图7a），占保护区总面积的80.9%。相关系数集中分布在0～0.5范围，占保护区总面积的72.9%。与气温呈负相关性的区域主要分布在保护区西南部。保护区所有像元生长季NDVI均与气温达到显著相关(P＜0.05)，99.4%的区域达到极显著相关(P＜0.001)。保护区大部分区域生长季NDVI与降水量呈正相关性（图7b），占保护区面积的80.7%。相关系数集中分布在0～0.5范围，占保护区总面积的77.4%。与降水量呈负相关性的区域空间分布特征不明显。保护区所有像元生长季NDVI均与降水量均达到极显著相关(P＜0.001)。表明气温与降水量对保护区生长季NDVI都有重要的正向影响，在保护区西南部的部分区域，可能是气温升高导致蒸发加大使生长季NDVI与气温呈现出负相关性。

图7 保护区2000—2020年逐年生长季NDVI与气温(a)和降水量(b)的相关性

3 结论与讨论

（1）保护区多年植被生长季NDVI均值为0.02～0.55，平均值为0.138，呈现由东南向西北逐渐减小的空间分布特征。

（2）保护区86.39%区域植被生长季NDVI呈缓慢的上升趋势，并且大部分达到显著水平，虽然有13.61%的植被生长季NDVI处于下降趋势，但只有很少的区域达到了显著性水平；植被总体呈现上升变好趋势。

（3）2000年以来保护区的生长季NDVI稳定性较差，较低稳定度和中等稳定度占整个研究面积的86.66%，植被生长季NDVI≤0.1的大部分区域处于较低稳定度。

（4）在植被生长季，保护区气候趋向于“暖湿化”；气温与降水量对保护区生长季NDVI都有重要的正向影响。

青藏高原植被生长季NDVI在0.1～0.9，多年均值为0.49[1]，青藏高原草地生长季NDVI变化范围在0.431～0.471[27]，羌塘保护区植被生长季多年NDVI的平均值为0.138，明显低于青藏高原整体植被NDVI和青藏高原草地NDVI，说明羌塘保护区植被状况较差。

从全球来看，有56.3%的陆地区域植被NDVI呈增加趋势，且具有显著季节变化趋势[28]。在中国，NDVI呈增加趋势的面积大约占53.8%，达到显著水平的面积占29.3%[29]。在青藏高原，NDVI升高区域面积约占57.1%[1]。羌塘保护区植被生长季NDVI呈增加趋势的面积占植被总面积的86.4%，其中达到显著水平的面积占保护区植被面积的45.15%，明显高于青藏高原整体水平。

综合NDVI变化趋势和气象要素变化趋势来看，NDVI下降区域与气温升高、降水量减少区域重叠较大，而NDVI上升区域主要与气温升高、降水量增加区域重叠较大，徐增让[30]的研究表明，保护区居民点和人口在实验区，缓冲区和核心区基本上无分布，说明保护区影响植被生长的人为干扰小，水热条件是保护区植被生长最重要影响因素，气候“暖湿化”有利于保护区植被改善。