川西高原植被NDVI动态变化特征及对气候因子的响应

谢慧君，张廷斌,2，易桂花，秦岩宾，李景吉，别小娟，范 馨

〔1.成都理工大学 地球科学学院,四川 成都 610059;2.国家环境保护水土污染 协同控制与联合修复重点实验室(成都理工大学),四川 成都 610059;3.成都理工大学 旅游与城乡规划学院,四川 成都 610059;4.成都理工大学 生态环境学院,四川 成都 610059〕

随着社会经济发展，全球气候发生显著变化，生态环境问题日益突出，生态系统的稳定在气候变化和人类活动干扰下受到影响[1-2]。全球气候变化导致高原冰川消退、雪线上升以及湿地空间分布发生改变、草地退化、土地沙化等一系列生态环境问题[2-5]，给草原生态带来严重的破坏。川西高原处于我国地势第一阶梯向第二阶梯过渡地带，是我国西部重要的生态屏障，也是四川省重要的生态保护红线区[6]。区内自然环境独特，生物多样性丰富，具有水源涵养、水土保持和气候调节等重要生态功能[6-7]。同时，由于其特殊的地理环境，区域内植被类型多样，生态环境脆弱，对气候变化响应敏感[7-9]。植被是陆地生态系统的重要组成部分，影响着全球的物质循环和能量交换，对维持气候系统、碳循环、水资源和生态环境的稳定发挥着重要作用[10-12]。植被的改善有助于提高水源涵养和水土保持功能[13]。因此，准确掌握川西高原植被动态变化特征，探索植被对气候变化的响应关系[14]，对指导区域生态环境建设和水土保持具有重要意义。

归一化差异植被指数(normalized difference vegetation index，NDVI)对气候变化极为敏感，是评价植被生长状况的重要指标[15-18]。利用NDVI在不同时间和空间尺度上评估植被变化特征及其对气候变化的响应已成为全球研究热点之一[19-21]。宋怡等[22]基于NDVI数据对中国寒旱区植被变化特征进行了分析，并研究了其对气候变化的响应，发现总体上植被呈改善趋势，植被变化与降水和温度均呈正相关关系；韩炳宏等[23]利用NDVI数据对青藏高原地区植被覆盖变化及其驱动因素进行了探讨，认为青藏高原植被整体呈稳定恢复状态，与气温和降水均显著相关，与气温的关系更密切。多年以来川西高原积极推进生态文明建设，开展退耕还林、生物多样性保护、水土流失治理等生态保护工程[24-25]，在这一过程中植被的生长状况怎样？如何掌握植被变化情况？植被变化趋势以及对气候变化的响应具有怎样的特征？这是川西高原进一步加强生态建设和水土流失治理所面临的关键问题。但目前针对川西高原植被动态变化特征及其对气候响应的研究并不多见。本文基于MODIS (moderate-resolution imaging spectroradiometer)MOD09A1数据反演川西高原植被NDVI，综合气象数据，DEM(digital elevation model)和植被类型等资料，开展2001—2017年川西高原植被的NDVI时空动态变化特征及其对气候因子的响应研究，为指导川西高原生态环境保护和水土保持提供科学依据。

1 材料与方法

1.1 研究区概况

川西高原位于四川省西部，处于27°57′—34°11′N和97°21′—104°26′E之间，总面积2.33×105km2，分布有31个县。区内地形复杂，起伏较大，海拔在790～6 500 m之间，地势从西北向东南逐渐降低(图1a)。区内水资源丰富，雅砻江、大渡河、岷江等水系发育；植被类型多样(图1b)，主要分布有高寒草甸、灌丛、针叶林、阔叶林等，其中针叶林类型之多为全国之冠[7,26]；川西高原气候类型多样，主要以寒温带气候类型为主，光照充足，气温较低，干湿季分明，降水较少；受地形影响，气候具有明显的立体变化特征[27-28]。

图1 川西高原DEM，气象站点分布和植被类型

1.2 数据来源及处理

1.2.1 数据来源 2001—2017年的地表反射率数据(MOD09A1)来自NASA(National Aeronautics and Space Administration)(https:∥ladsweb.modaps.eosdis.nasa.gov/)搭载于Terra卫星上的EOS/MODIS传感器，包含了7个波段，其中Band1和Band2的空间分辨率为250 m，其余波段为500 m，时间分辨率8 d，该产品经过了辐射校正、大气校正等预处理。气象数据是由中国气象科学数据共享服务网(http:∥data.cma.gov.cn)提供的2001—2017年川西高原及周边37个气象站点逐日数据。植被类型数据采用中国科学院资源环境科学数据中心(http:∥www.resdc.cn)发布的1∶100万全国植被类型数据集。DEM数据来自地理空间数据云平台(http:∥www.gscloud.cn/)的SRTM(Shuttle Radar Topography Mission)数字高程数据产品，空间分辨率为90 m。

1.2.2 数据预处理 利用MRT(MODIS reprojection tools)工具对MOD09A1数据进行格式转换和重投影等预处理，根据质量控制文件剔除数据中的异常值，利用MOD09A1的红外(Band1)和近红外(Band2)波段计算植被NDVI值：

(1)

式中：NIR为近红外波段反射率；R为红波段反射率。对获得的植被NDVI值乘以0.000 1缩放因子并按川西高原边界范围进行裁剪。由于川西高原海拔5 000 m以上的区域大多覆盖着永久冰川[20]，本研究不予讨论，利用DEM数据剔除5 000 m以上的区域。相对于已有的16 d NDVI数据产品，MOD09A1数据计算的NDVI数据时间分辨率更高，采用MVC(maximun value composite)方法计算2001—2017年逐月植被NDVI值，进而计算逐年年均植被NDVI值可以更加有效的削弱云的影响。

利用同质性检验对不可用气象站数据进行剔除[29-30]。现有研究中常采用克里金法(Kriging)和反距离权重法对降水量插值[31-32]，采用薄板样条函数法(thin plate spline，TPS)和Kriging对气温插值[7,33]。本文利用气象站数据验证法对空间插值方法进行验证[34]，并考虑到气温和降水在空间分布上存在差异，且气温受海拔影响[34-35]，最终选择Kriging方法对降水进行插值，对气温的栅格化选择以DEM作为协变量的TPS方法进行空间插值[34]。研究中所使用的植被类型、DEM和气象数据等均经过投影变化、裁剪、重采样等数据处理，将空间分辨率统一为500 m。

1.3 研究方法

1.3.1 Theil-Sen media趋势分析与Mann-Kendall统计检验 采用Theil-Sen media趋势分析与Mann-Kendall统计检验[36]相结合，逐像元分析2001—2017年川西高原植被NDVI变化趋势。其中，Theil-Sen media趋势分析是一种非参数统计趋势分析方法[37]，有助于减少缺失或异常值的干扰[38]，计算公式为：

(2)

式中：xi，xj分别为第i和j年的NDVI值；n为数据集合的时间长度；β用于量化单调趋势，当β>0，表明植被NDVI值随时间呈上升趋势，β<0，表明植被NDVI值随时间呈下降趋势。

Mann-Kendall统计检验是一种非参数检验方法，该方法所用样本不需服从一定分布[39-40]，不受少数异常值干扰[38]，用于评估趋势的重要性[41]。在进行趋势检验时，将逐年每个像元的NDVI值看作一组时间序列，参数Z作为统计检验量[40]，具体计算方法为：

(3)

其中：

(4)

(5)

(6)

式中：sgn为函数符号。本文在α=0.01置信水平上判断植被NDVI变化趋势的显著性[37]，当|Z|>Z1-α/2，表示数据变化趋势显著；当|Z|≤Z1-α/2，表示数据变化趋势不显著。

1.3.2 Hurst指数分析 基于Hurst指数方法分析川西高原植被NDVI未来趋势变化，该方法可判断时间序列未来变化的趋势，用于定量描述长时间序列信息相关性[39,42]。Hurst指数预测未来趋势分为3种情况：当0.5

1.3.3 偏相关分析 应用偏相关分析研究川西高原植被NDVI对气温和降水的响应，计算公式为[44]：

(7)

式中：Rxy,z为自变量z固定后因变量x与自变量y的偏相关系数[52]；Rxy，Rxz，Ryz为两因子间的相关系数；x为NDVI值；y为气温值；z为降水值。

2 结果与分析

2.1 川西高原植被NDVI时空分布特征

川西高原近17 a植被NDVI均值介于0～0.85之间，平均值为0.486，植被具有较好的长势，空间上呈现由西北向东南逐渐升高的分布格局(见封2，附图2)。由于川西高原地形差异大，海拔由西北向东南逐渐降低，植被呈现明显空间异质性。NDVI低值区多出现在4 000 m以上的高海拔区域，主要分布在石渠县北部、德格县中部、新龙县与白玉县交界边缘、理塘县与稻城县交界处、康定南部、黑水县与红原县交界处、松潘县东部、九寨沟南部以及小金县、理县和汶川县的交界地带，植被类型多为高山植被、高寒草甸和灌丛。NDVI高值区多出现在研究区东部和南部海拔3 000 m以下地区，主要分布在若尔盖西南部、九寨沟东南部、雅江县中部、丹巴县东北部、泸定县东部以及大渡河、雅砻江和岷江流域，植被类型以针叶林、阔叶林和灌丛为主。结合川西高原地形特征，将川西高原高程划分为5个高程区间(<1 000 m，1 000～2 000 m，2 000～3 000 m，3 000～4 000 m，4 000～5 000 m)，分区统计植被NDVI均值随高程变化情况(见表1)。结果表明，近17 a来川西高原植被NDVI均值随高程增加呈现先增加后降低的特点(表1)。3 000 m高程以下区域面积占比7.07%，该区域植被NDVI随高程增加呈逐渐增加趋势；3 000 m高程以上区域面积占比较大为92.93%，该区域植被NDVI随高程增加呈降低趋势。植被NDVI均值在2 000～3 000 m高程区间出现最高值(0.63)，主要植被类型为针叶林和阔叶林；最低值出现在4 000～5 000 m高程区间(0.41)，区内植被类型以高寒草甸和高山植被为主。在1 000～2 000 m和3 000～4 000 m高程区间，植被NDVI均值分别为0.61和0.58，该区植被类型分别以灌丛和高寒草甸为主；在<1 000 m高程区间，由于人类活动干扰较大，植被NDVI相对较低。川西高原近17 a植被NDVI均值处于波动上升趋势(见图2)，增长速率为0.01/(10 a)(p=0.001)，最小值出现在2012年(0.475)，2015年达到最高(0.508)。其中2006年、2009—2011年、2013—2017年植被NDVI均值高于多年平均值，在0.489～0.508之间波动。

图2 川西高原2001-2017年植被NDVI均值年际变化

表1 川西高原植被NDVI垂直分布特征

2.2 川西高原植被NDVI动态变化特征

2.2.1 川西高原2001—2017年植被NDVI变化趋势分析 为分析川西高原2001—2017年植被NDVI变化趋势特征，根据Theil-Sen media趋势分析方法得到川西高原近17 a植被NDVI趋势图(见图3)，基于Theil-Sen media趋势分析和Mann-Kendaill统计检验结果将研究区植被NDVI趋势划分为显著退化(β<0，α<0.01)、不显著退化(β<0，α>0.01)、不显著改善(β>0，α>0.01)和显著改善(β>0，α<0.01)4种类型。

由图3，表2可知，2001—2017年川西高原植被NDVI整体以改善趋势为主，显著改善面积占比12.09%，集中分布在若尔盖县和九寨沟县、松潘县中部、茂县和黑水县交界处等研究区东北部边缘地带；不显著改善面积占比达71.41%，分布在整个川西高原。不显著退化面积占比15.88%，呈零星分布特征，主要分布在石渠县、汶川县、雅江县和稻城县等区域；显著退化区主要分布于雅江县北部、康定县与雅江县交界处，面积占比0.61%。其中，2008年汶川地震导致灾区原生生态环境遭到大面积破坏，植被生态环境退化，雅江县中部、康定县与雅江县交界处、石渠县中部以及汶川县东部植被NDVI表现出明显退化趋势。另外，色达县中部和稻城县中部、南部的植被NDVI表现出轻微退化趋势。

表2 川西高原2001-2017年植被NDVI变化趋势统计

图3 川西高原2001-2017年植被NDVI趋势和植被NDVI变化趋势类型

2.2.2 川西高原植被NDVI未来变化趋势 川西高原近17 a的Hurst指数空间分布如图4所示。川西高原植被NDVI的Hurst指数介于0.16～1之间，平均值为0.69，研究区植被NDVI变化趋势总体具有较强的持续性。通过叠加川西高原植被NDVI变化趋势和Hurst指数，可将川西高原植被NDVI未来变化趋势划分为未来退化趋势、持续性改善趋势、未来改善趋势和持续性退化趋势4种类型(图4)。由图4和表3可知，川西高原植被持续性变化趋势主要以持续性改善为主，面积占比78.33%，分布在川西高原全区；约1.95%的植被NDVI区域面积由退化向改善方向变化，主要分布在石渠县西部。约5.18%的植被NDVI由改善向退化方向变化，零星分布在全区；14.54%的植被NDVI呈持续性退化趋势，主要分布在汶川县东部、雅江县北部、石渠县中部、色达县中部和稻城县全境，若尔盖县和红原县等地分布较少。

表3 川西高原植被NDVI持续性变化趋势统计

图4 川西高原植被NDVI Hurst指数空间分布和持续性变化趋势分布特征

2.3 川西高原植被NDVI变化对气候因子的响应

2.3.1 气候因子年际变化特征 图5表示川西高原2001—2017年平均气温和年累积降水量变化特征。由图5可以看出，川西高原近17 a平均气温总体呈波动上升趋势，平均气温在2009年出现最高值2.48 ℃，在2004年出现最低值1.26 ℃，多年平均值为1.94 ℃，整体以0.42 ℃/(10 a)(p=0.023)的速率增温。研究区年累计降水量在2001—2017年整体上升趋势不明显(图5)，波动较大。其中年累计降水量最小值出现在2002年为545 mm，在2012年达到最大值806.64 mm，多年平均累计降水量为713.49 mm，整体以29.7 mm/10 a(p=0.286)的速率增加。

图5 川西高原2001-2017年平均气温和年累计降水量

2.3.2 植被NDVI对气候因子的响应分析 川西高原2001—2017年植被NDVI与气温的偏相关系数值介于-0.92～0.9之间，正、负相关区域面积占比分别为75.7%，24.3%。研究区约有12.8%的区域通过了p<0.05显著性水平检验，显著正相关和显著负相关区域分别占研究区总面积的11.6%和1.2%(见图6)，这表明整体上川西高原植被NDVI与气温呈正相关关系。植被NDVI与气温呈显著正相关关系的区域主要分布在石渠县、若尔盖县、九寨沟县、红原县、阿坝县以及汶川县，主要植被类型以高寒草甸为主，该区域属高寒区，生态环境极为敏感脆弱，气温变化对植被生长的影响较大。植被NDVI与气温呈显著负相关关系的区域主要集中在雅江县和新龙等中部海拔相对较低处(见图6)，主要植被类型为高寒草甸、灌丛和阔叶林。川西高原植被NDVI与降水的偏相关系数分布在-0.93～0.93之间，正、负相关区域分别占研究区总面积的44.6%和55.4%。有11.6%的区域通过了p<0.05显著性水平检验，显著正相关和显著负相关的区域面积分别占比6.3%和5.3%(见图6)，整体上川西高原植被NDVI与降水呈负相关关系，表明降水对有些地区的植被生长起到一定抑制作用。植被NDVI与降水呈显著正相关关系的区域主要分布在若尔盖县、红原县、阿坝县、色达县、甘孜县、以及德格县等北部地区，显著负相关关系主要分布在雅江县、丹巴县、乡城县和巴塘县等中部和西南部区域(见图6)。

图6 川西高原植被NDVI与气温和降水偏相关系数的t显著性检验

3 讨 论

本研究基于NDVI对川西高原植被动态变化特征进行分析，发现川西高原NDVI整体呈改善趋势，并且大部分区域呈持续改善趋势。这一结论说明该地区多年来积极推进生态文明建设，实施退耕还林、封山育林等生态保护工程[24-25]效果显著，对植被生长起到了积极作用。另外川西高原建立了湿地公园、水源涵养、水土保持功能区和森林公园等各类保护地，这可能是相应区域植被得到改善的主要原因[20,45-47]。但在川西高原局部地区植被NDVI呈现退化趋势，尤其是在冰川覆盖以及人类聚集区。这可能由三方面影响造成的，首先，高寒区的植被对气候变化响应敏感，近来气温升高、冰川冻土消融使得高海拔区域植被生长受到影响[48-50]。川西高原海拔5 000 m以上区域分布有格聂山和贡嘎山等大面积冰川[20]，地形起伏大，植被覆盖较低，易发生土壤侵蚀，造成水土流失[50]，植被呈现出退化趋势。其次，由于川西高原生态环境脆弱，地质构造复杂，地质灾害频发，原生生态环境极易遭到破坏[20,25,53]。2008年汶川地震导致大面积地表裸露，植被根系生长环境遭到破坏，造成植被生态系统稳定性变差，水土流失严重。最后，川西高原部分城镇及周边等人口密集地区人类活动干扰加剧，植被生长环境受到一定破坏[52-54]。川西北高原石渠县境内超载放牧，土地沙化严重[54]，造成植被覆盖较低。雅江县等城镇周边由于城镇化扩张、基础设施建设以及旅游开发等，挤占了植被生长空间[20]，植被呈现出持续退化趋势。

另外，植被动态变化与气候变化密切相关[55]。本文分析了川西高原植被NDVI对气温和降水的响应关系，整体上川西高原植被NDVI受气温和降水共同影响，在空间上存在一定差异。由于川西高原特殊的地形地貌和复杂的气候条件，在北部的高寒地带，气候类型为高原高寒气候，植被生态环境脆弱，植被生长受水热条件限制，但随着近年来气温升高，降水增加，植被与气温和降水主要呈正相关关系；而西南部分地区全年气温均值相对偏高，降水较少，多为干热河谷气候，在植被生长较为适应气温状况的前提下，降水成为影响植被生长的主要因素。整体上川西高原气温升高，降水增加趋势不显著，地表蒸发有所增加，土壤水分含量降低，并且受川西高原地形影响，降水虽然有保湿作用，但易造成水土流失[49]，破坏植被生长环境，导致降水对有些地区的植被生长起到一定抑制作用。另一方面，极端气候的发生影响植被生长，2012年川西高原气温较低，降水量最大，该年份发生涝灾[56]而导致当年植被NDVI值较低。此外，文中植被NDVI与气候因子偏相关分析的结论与郑杰等[51]、曾波等[57]、荣欣等[20]的研究结果基本一致，但本文与荣欣等[20]川西高原植被EVI(enhanced vegetation index)与气候因子偏相关分析在空间上存在一些差异。大量研究表明NDVI和EVI描述高植被覆盖和低植被覆盖植被生长情况的能力存在差异[58]。川西高植被类型多样，包括针叶林、灌丛、草甸和高山植被等植被；西北部地区海拔相对较高，植被覆盖相对较低，类型单一；东南部海拔相对较低，植被覆盖较高，类型多样。由此可见，两种植被指数的不同可能是导致对气候因子的响应存在差异[58-59]的主要原因。另外，参与计算的气象站数量不同[11,59]也可能是导致差异的原因之一。

4 结 论

本文基于MOD09 A1数据，反演了川西高原2001—2017年植被NDVI，结合气象数据、植被类型数据和DEM数据，利用Theil-Sen media趋势分析、Mann-Kendaill统计检验、Hurst指数分析以及相关分析等多种方法，分析了川西高原年际植被NDVI的时空分布特征、变化趋势、未来趋势变化特征以及植被NDVI变化对气候变化的响应。

(1) 川西高原2001—2017年植被NDVI均值为0.486，受地形影响明显，在空间上呈现由西北向东南逐渐升高的分布格局，在垂直分布上随海拔升高呈现先增加后降低的特点，低值区多分布在4 000 m以上的高海拔区域，高值区主要分布在研究区海拔3 000 m以下的东部和南部，近17 a植被NDVI均值以0.01/10 a的速率波动上升。

(2) 川西高原近17 a植被NDVI整体以不显著改善趋势为主，面积占比达71.41%，有12.09%的区域呈现显著改善特征，不显著退化占15.88%，显著退化面积仅占0.61%。研究时段植被NDVI的Hurst指数均值为0.69，总体变化趋势持续性较强，持续性改善面积占比78.33%，未来改善占1.95%，未来退化占比5.18%，持续退化占比14.54%。

(3) 川西高原近17 a气温增温速率为0.42 ℃/10 a，多年平均累计降水量增幅为29.7 mm/10 a。研究时段内植被NDVI与气温偏相关关系的正、负相关区域面积分别占75.7%和24.3%，与降水偏相关关系的正、负相关区域面积分别占44.6%和55.4%。整体上植被NDVI受气温的影响大于降水，与气温呈正相关关系，与降水呈负相关关系；空间上存在差异，研究区北部植被NDVI主要受气温和降水影响，西南部植被NDVI主要受降水影响。