气候变化下珠穆朗玛峰自然保护区南北坡植被变化差异

卞梨交, 李景吉,3, 徐彬妮, 向 莹

(1.成都理工大学 生态环境学院, 成都 610059; 2.成都理工大学 环境与土木工程学院, 成都 610059;3.地质灾害防治与地质环境保护国家重点实验室, 成都理工大学, 成都 610059)

气候变化对生态系统的影响和反馈已成为国内外学者关注的重要生态学问题之一[1-2]，全球气候变暖背景下的极端气候事件将更加剧烈和频繁，这将严重威胁植被分布[3-6]。开展长期观测和模型模拟，阐明植被对气候变化的响应，是当前全球气候变化研究的热点问题之一[7-10]。

植被是指示陆地生态系统的重要指标[11-12]，植被生长与气候变化密切相关[13]，研究植被变化对全球气候变化响应的时空动态特征，对预测未来生态系统动态至关重要[14]。归一化植被指数(NDVI)是表征植被覆盖度和植被变化的重要指标之一[15-16]，与植被覆盖度呈正相关，能够直观指示植被覆盖度和植被生长状态[17]，已成熟应用在全区不同地区的植被对气候变化响应关系研究中[15-16]。Kalisa等[15]分析了东非地区NDVI在不同的空间格局上的断裂点以及NDVI趋势的系统性变化；Xu等[14]基于AVHRR NDVI数据讨论了1982—2011年中国的植被生长动态及其与气候因子的关系；Zhou等[13]基于NDVI探讨了中国植被的历史动态变化，并基于多元回归模型构建了植被的预测模型。

青藏高原是全球气候变化的“敏感区”和“响应器”[18-21]，是研究陆地生态系统对气候变化的响应的最理想地区之一[22]。已有研究发现，青藏高原气候变化程度超过北半球大部分地区甚至整个世界的气候变化[23-24]。近年来，青藏高原已成为国内外学者研究全球气候变化的热点地区，不同学者分别选择雅鲁藏布江流域[25]、西藏地区[26]、藏北高寒草原[27]、珠峰保护区等[28]青藏高原典型区，开展植被NDVI动态变化研究，以尝试揭示全球高原植被变化对气候变化的响应机制，但对受地形条件限制下空间异质性强烈的复杂地理单元植被动态对气候变化的差异性响应机制的解释尚缺乏实际证据。珠穆朗玛峰国家级自然保护区(以下简称“保护区”)位于青藏高原南缘，平均海拔约4 200 m，南北坡地形、气候和植被差异非常大，南坡以深切峡谷地貌为主，气候类型以亚热带、温带湿润气候为主，植被分布具有明显的“三向地带性”[29]；北坡多为高原地貌类型，气候类型为干旱、半干旱气候，植被分布以高原灌丛、草甸和草原为主。珠峰保护区是研究复杂地理单元植被动态对气候变化的差异性响应典型区，目前对该区南北坡植被NDVI与气候响应差异的研究鲜有报道。鉴于此，本研究基于趋势分析和相关性分析，研究近20年来(2000—2018年)保护区南北坡植被动态变化及其对气候变化的差异性响应特征，对解释复杂地貌单元植被动态对气候变化的差异性响应机制具有重要意义。

1 研究区概况

保护区地处我国西藏与尼泊尔交界处，位于84°27′—88°E，27°48′—29°19′N，分布海拔1 448～8 844 m，全区面积约3.4×104km2，辖定日、吉隆、聂拉木和定结4县(图1)。保护区南、北坡气候、植被差异显著，南坡受印度洋湿暖气流影响，降水充沛，形成山地森林生态系统，年均气温7～10℃，年降水量大于1 000 mm；北坡受“雨影区”影响，气候干燥，具有大陆性高原气候特点，形成半干旱灌丛、草原生态系统，年均气温2～5℃，年均降水小于300 mm。

图1 研究区地理位置

2 数据与方法

2.1 数据源

MODIS传感器能够在宏观格局上提供植被信息[30]，本文获取的NDVI数据为MODIS 13A3产品，该产品是由美国国家航空航天局(https:∥ladsweb.modaps.eosdis.nasa.gov)提供，空间分辨率为1 km×1 km，时间分辨率为30 d，本次下载了2000年2月—2018年12月的所有影像。

保护区2000—2018年气温、降水栅格数据[31]来源于国家科技基础条件平台—国家地球系统科学数据中心(http:∥www.geodata.cn)，月平均气温和月累计降水数据是来源于中国气象数据网(http:∥data.cma.cn/data/cdcdetail/dataCode/SURF_CLI_CHN_MUL_DAY_V3.0.html)的“中国地面气候资料日值数据集(V3.0)”，并通过计算统计得到。高程数据采用ASTER GDEM 30 M分辨率数字高程数据，来自地理空间数据云(http:∥www.gscloud.cn/home)。讨论人类活动影响中的人工生态系统数据是基于中科院资源环境科学与数据中心(https:∥www.resdc.cn)，并根据Google Earth高分历史影像和野外调查数据调整，并将农田和聚落生态系统合并为人工生态系统以表示人类活动区域。

野外调查数据为2009年、2010年、2019年对保护区植被调查数据，共计278个样方(图1)，用于讨论不同NDVI值区间对应分布的主要植被种类。调查方法为：在已设的调查样线基础上，以样方调查为主，同一植被类型做3个重复样方，调查记录：群落生境、乔木层建群种、郁闭度等；灌木层建群种、丛数、盖度等；草本层建群种、盖度、高度；样方布设根据群落的特点，乔木林群落设置10 m×10 m的调查样方；灌木林群落2 m×2 m的调查样方；草本群落1 m×1 m的调查样方。

2.2 方 法

通过最大值合成法(MVC)处理NDVI数据，在30 d的NDVI的基础上得到年最大NDVI(NDVImax)，用来表征当年植被生长的最佳状况，并将NDVI<0的值设置为nodata并处理为0～1[32]。把2000—2018年按等时间间隔分为2000—2009年和2009—2018年两个时间段进行分析，采用Theil-Sen(TS)趋势分析估算进行趋势分析，Mann-Kendall检验法进行检验，对两个时间段温度、降水和NDVI栅格数据进行趋势分析，并采用Hurst指数法结合TS趋势分析来预测未来植被变化趋势。将研究区两个气象站点(聂拉木、定日)分别做1 km缓冲区，统计出2000年2月—2018年12月缓冲区内各月平均NDVI，将2000年02月至2018年12月平均气温、降水量分别与月均NDVI做Person偏相关分析，结果见表1，并基于像元对2000—2018年NDVImax与温度、降水栅格数据进行偏相关分析，student t检验法进行检验。

2.2.1 趋势分析 Theil-Sen(TS)中位数趋势分析和Mann-Kendall(MK)检验可以结合在一起用于分析长期植被序列的变化趋势，该方法不受样本分布方式的影响，也不受少数异常值的干扰，是检验时间序列趋势的有效方法[33-35]。TS公式为：

(1)

当βNDVI>0时，NDVI呈现上升趋势，否则NDVI呈现下降趋势。

Mann-Kendall是一种非参数统计检验方法，标准统计量Z值用于指示趋势的强度。

对于NDVIj时间序列，

(2)

(3)

(4)

(5)

式中:NDVIi、NDVIj分别代表第i年和第j年的NDVI值,而n代表时间序列的长度。若|Z|>1.96,则表明NDVI时间序列在0.05水平上有显著变化;|Z|<1.96,则表明NDVI时间序列不显著变化。因此,本文将植被覆盖趋势可分为4种类型:在0.05水平上显著增加(βNDVI>0,|Z|>1.96);不显著增加(βNDVI>0,|Z|<1.96);不显著减少(βNDVI<0,|Z|<1.96);在0.05水平上减少(βNDVI<0,|Z|>1.96)。

2.2.2 相关性分析 为了探讨NDVImax与气候因子的关系，在像元尺度上采用了偏相关系数法[36]。偏相关系数[37]是研究多元变量的相关性的一种方法，该变量用来度量在其他变量控制下两个变量直接线性相关性[38-39]。根据定义，在的影响保持不变的情况下，变量和变量的偏相关系数计算为：

(6)

式中:rxy,rxz,ryz分别表示3个变量和z的Pearson相关系数。Pearson相关系数的计算公式为:

(7)

式中：N是总年份数；xt是第t年的x变量的值；yt是第t年的y变量的值；x是所有年份的平均值；y是所有年份的平均值。

本研究的偏相关系数的显著性是通过使用Student t检验在0.05显著性水平上来检验各像元NDVI与气温和降水的偏相关系数，并将相关程度分为以下4类：在0.05水平上显著正相关、不显著正相关、不显著负相关和在0.05水平上负相关。

2.2.3 Hurst指数趋势预测 本研究采用Hurst指数法[40-41]结合TS趋势分析来预测珠峰保护区植被未来变化趋势，并将其分为5种变化类型：持续增加、由增到减、持续减少、由减到增、独立变量(H=0.5)。用H表示Hurst指数值，当0

3 结果与分析

3.1 NDVI随时间变化特征

3.1.1 NDVI年际趋势 6分析近20 a来(2000—2018年)保护区NDVI年际变化趋势发现(图2)，保护区年均NDVI呈波动上升趋势，平均增长率为0.000 8/年(R2=0.15，p=0.100 7)。2000—2009年的年平均NDVI呈波动下降趋势，平均为-0.001 9/年(R2=0.15，p<0.001)，2009—2018年NDVI以0.003 0/年的平均增率呈显著上升趋势(R2=0.16，p<0.001)。

图2 2000-2018年年均NDVI变化

3.1.2 年内NDVI与气候变化的线性相关性 保护区月际NDVI与温度的相关性大于与降水的相关性(表1)。在聂拉木1 km缓冲区内，以降水为控制变量时，NDVI与气温的偏相关系数为0.853(p<0.01)；以气温为控制变量时，NDVI与降水未表现出显著相关性。在定日站1 km缓冲区内，以降水为控制变量时，NDVI与气温的偏相关系数为0.401(p<0.01)；以气温为控制变量时，NDVI与降水的偏相关性系数为0.134(p<0.05)。

表1 2000-2018年月际温度和降水与NDVI的偏相关系数

3.2 NDVI空间变化特征

3.2.1 NDVI分布格局 NDVI值大于0.6的地区主要分布在保护区南坡，主要位于陈塘镇、曲当乡、绒辖乡、樟木镇、聂拉木镇和吉隆镇，少量分布在贡当乡、宗噶镇和日屋镇的南部(图3B中R2，R3，R4，R5区域)，植被类型以乔木为主，优势物种包括：糙皮桦(Betulautilis)、乔松(Pinusgriffithii)、尼泊尔桤木(Alnusnepalensis)、高山栎(QuercussemicarpifoliaSmith)、喜马拉雅冷杉(Abiesspectabilis)等。NDVI值范围为0.1～0.4的地区主要分布在保护区北坡(图3B中R1和R7区域)，植被类型以灌草丛为主，优势物种包括：杜鹃(Rhododendronspp.)、忍冬(LonicerajaponicaThunb.)、蔷薇(Rosaspp.)、金露梅(P.fruticosa)、苔草(Carexsp.)等；NDVI值范围为0.4～0.6的地区主要分布在北坡的加措乡中部的盆吉乡、琐作乡、差那乡和折巴等地的北部海拔较低区域，NDVI值在0.1以下的区域主要分布在海拔高的区域(表2)。

2000—2009年保护区内整体植被覆盖率低，大部地区NDVI值均在0.1～0.4，对气候变化比较敏感；其中NDVI值0.1～0.2的区域面积占比最大(达37.45%)，NDVI值为0.2～0.4的区域面积占比为33.23%；NDVI值大于0.6的地区面积占比仅为6.02%，；NDVI值为0.4～0.6的区域面积占比为8.11%，主要分布在北坡的加措乡中部，盆吉乡、琐作乡、差那乡和折巴乡的北部。2009—2018年保护区NDVI分布情况与2000—2009年基本一致，但NDVI总体上有所上升。表现为：NDVI值范围在0～0.1，0.1～0.2的区域均有减少，NDVI值范围在0.2～0.4，0.4～0.6，0.6～1的区域均有增加(图3)。

3.2.2 NDVI变化趋势及其预测分析 2000—2009年，保护区NDVI总体呈下降趋势，变化速率为0.001 7/a速率；其中70.88%的区域NDVI下降，显著下降区面积占比5.05%；29.11%的区域NDVI增加，显著增加区面积占比仅有1.47%；退化区主要分布在北坡地区，改善区主要分布在南坡地区。2009—2018年，保护区NDVI总体呈上升趋势，变化速率为0.003 0/a，有77.92%的区域NDVI增加，显著增加区域面积占比为7.06%；退化区主要分布在东北部，改善区全域零散分布，其中显著增加区主要零星分布在折巴乡、差那乡及保护区中部(图4)。

注：(A)为2000—2009年NDVImax多年平均值，将NDVI分为0～0.1，0.1～0.2，0.2～0.4，0.4～0.6和大于0.6这5个区间；(B)为2009—2018年NDVImax多年平均值。根据样方调查点的分布情况和NDVI值的分布情况，将该区域划分出7个区，R1，R2，R3，R4，R5，R6，其余地区划分为R7。

表2 主要植被类型分布及其对应的NDVI值

利用Hurst指数和趋势分析预测保护区未来NDVI变化趋势发现(图5)，保护区未来NDVI变化不稳定，趋势不一致的像元占总像元的70.29%，且58.68%的区域可能会出现NDVI下降趋势。

图4 NDVI趋势变化空间异质图

图5 趋势预测类别空间分布

3.3 NDVI与气候因子的关系

对气候因子变化的趋势分析发现，保护区气温在2000—2009年呈增加趋势，2009—2018年呈下降趋势。2000—2009年北坡聂拉木县和定日县降水呈上升趋势，其他地区呈下降趋势；2009—2018年降水在大部分地区呈上升趋势，而定日县和定结县的南部呈下降趋势(图6)。

对保护区NDVI与气候因子相关性分析发现，NDVI与降水的相关性大于与温度的相关性(图7)。保护区周缘区NDVI与温度表现出正相关性，聂拉木县中部NDVI与温度表现出高的负相关；59.02%的区域NDVI与温度呈现出正相关，其中显著正相关区域面积占比仅为1.79%；40.98%的区域NDVI与温度呈现出负相关，显著负相关区域面积占比为14.57%，主要分布在北坡地区。保护区大部地区NDVI与降水呈现出高的正相关性，尤其是在东部的定结县；62.13%的区域NDVI与降水呈现出正相关，显著正相关区面积占比为8.21%；37.87%的区域NDVI与温度呈现出负相关，显著负相关区面积占比为2.94%。

4 讨 论

4.1 NDVI分布特征

研究发现保护区南坡NDVI值相对较高，主要原因有两个方面：(1) 受印度洋湿暖气流的强烈影响，保护区南坡具有海洋性季风气候特征，区内主要分布山地森林生态系统[28]，多分布半常绿阔叶林、常绿阔叶林、常绿针叶林、硬叶常绿阔叶混交林、落叶阔叶桦木林、高寒杜鹃灌丛和高寒常绿针叶灌丛[43]；(2) 受喜马拉雅地质隆升影响，南坡地区地形复杂、海拔落差大、坡度陡，受人类活动的干扰小[28]。

北坡NDVI值相对较低，主要原因有：(1) 受喜马拉雅山脉的屏障作用，北坡形成“雨影区”，降水较少，受全球气候变化影响，北坡呈现干旱化趋势，土地沙化、盐碱化现象明显[44]，植被生长受到限制；(2) 北坡多为高原区，植被类型以高山草甸、草原为主，区内畜牧业相对发达，人类放牧和城镇化建设增加了对植被NDVI不利干扰[28]。

4.2 NDVI趋势及其预测

保护区南坡NDVI在2000—2018年总体上呈一定的上升趋势，该区多为保护区核心区，受人类活动影响小，大面积分布森林生态系统对气候变化具有较强的缓冲能力[43]；北坡地区NDVI在2009—2018年呈现上升趋势，这与2002年以来西藏自治区大力开展退耕还林、还草工程有关；在2009年之前，受生态工程时间较短影响，退耕还林、还草工程尚未发挥效果[44]，2009年之后，生态工程的生态效益逐渐体现。Hurst指数预测结果表明，未来珠峰保护区NDVI变化不稳定，NDVI由增加到减少的面积占比最大，这意味着极端的气候条件可能将持续影响珠峰植被。

4.3 NDVI对气候因子的响应

4.3.1 年内NDVI在时间上对气候因子的响应 本研究发现，2000—2018年聂拉木站、定日站1 km缓冲区的月平均NDVI与气温均表现出了显著的相关性，高于月均NDVI与降水的相关性，说明区域年内植被生长受气温的影响显著于降水，结果与中国北方相关研究结果表现一致[45]。

比较聂拉木站、定日站1 km缓冲区的月平均NDVI与气温相关性发现，聂拉木站大于定日站，这可能与区域降水差异性有关。已有研究表明，土壤水分含量相对高，降水充沛，温度的升高有利于植被生长季的延长和干物质量的积累，从而增强温度变化对植被的影响[46-47]；干旱区、半干旱区生长季升高的温度通常会增加地表水的蒸发，这对植被的生长有所限制，尤其是灌木和稀疏植被的生长[48-50]。聂拉木站1 km缓冲区月平均NDVI与降水并未表现出相关性，这可能与南坡足够充沛的降雨有关，且有关研究显示植被生长对降水的响应具有时滞性[15-17]；定日站1 km缓冲区NDVI与降水表现出一定相关性，说明干旱区、半干旱区的植被生长主要受生长季的降水量和蒸散发的影响，区域降水增加能够促进植被生长[48]。

注：A，B分别表示温度在2000—2009年和2009—2018年的趋势空间分布；C，D分别表示降水在2000—2009年和2009—2018年的趋势空间分布。

图7 NDVI与气候因子偏相关显著性空间异质图

4.3.2 年际NDVI在空间上对气候因子的响应 保护区南坡在2000—2009年温度上升、降水减少，2009—2018年温度下降、降水增加，在2000—2018年温度、降水与NDVI呈现出的相关性均不显著，保护区NDVI在大部分地区呈现一定程度的增加趋势，这一现象与自然植被分布和人类扰动较少有关。南坡主要分布山地森林生态系统[28]，植被自然恢复能力和对气候变化缓冲能力较强[51]，分布区域多为核心区，受人类活动的干扰小[28]。

保护区北坡地区在2000—2009年气温上升，多地降水量减少，NDVI大部分地区有下降趋势；在2009—2018年气温总体下降，降水增加，NDVI在多地有上升趋势；NDVI与气温呈现负的相关性，与降水呈现正的相关性。这一现象与北坡植被对气候的敏感性较高有关，北坡气候干燥、降水量较少，温度增加促进降水蒸发，植被对气候变化较为敏感[44]。

4.4 NDVI对人类活动的响应

根据趋势分析结果，2000—2009年保护区NDVI大部分地区呈现下降趋势。为明确保护区人类活动对NDVI影响程度，筛选该时间段内NDVI显著下降但与气温和降水相关性不显著的区域(图8筛选区)，并叠加2000年、2010年人工生态系统类型数据，发现橘黄色区域多零星分布在县域居民地附近(尤以北坡定日和定结县居民地附近最为明显)，且与人工生态系统分布较为一致；此区域植被NDVI表现出显著退化趋势，但NDVI变化与气候因子的相关性较弱，且人工生态系统面积在此处明显增加，这说明人类活动是导致2000—2009年保护区NDVI退化的一个重要因素，区域植被NDVI对人类活动十分敏感。

图8 人类活动敏感区

5 结 论

(1) 2000—2009年，保护区植被NDVI总体呈下降趋势，2009—2018年，呈上升趋势；年内NDVI与温度相关性强于降水。

(2) 保护区南北坡植被差异性明显，整体植被覆盖率低。南坡地区NDVI值基本大于0.6且变化稳定，主要分布着喜马拉雅冷杉(A.spectabilis)等高大乔木，对外界干扰具有较强的抵御力；北坡地区NDVI值基本保持在0.1～0.4且容易波动，主要分布着苔草(Carexsp.)等灌草丛，对气候变化比较敏感。

(3) 南坡地区NDVI在2009年前后两个时段总体上呈现上升趋势；北坡地区NDVI在2000—2009年呈下降趋势，在2009—2018年呈现上升趋势。Hurst指数预测保护区未来NDVI变化不稳定，58.68%的区域可能会出现下降趋势。

(4) 总体上，保护区年际NDVI与降水的相关性大于温度。南坡地区NDVI与气候因子的相关性弱，北坡地区NDVI与温度呈现负相关，与降水呈现正相关。NDVI显著下降但与气温和降水相关性不显著的区域在北坡定日和定结县居民地附近分布最为明显，人类活动是导致此区域NDVI在2000—2009年退化的重要因素。