武陵山区植被生长季NDVI时空变化及其对气候变化的响应

刘 恒, 汤弟伟, 孙 毅, 宋鄂平, 常 胜, 刘小芳

(1.湖北民族大学 林学园艺学院, 湖北 恩施 445000; 2.恩施土家族苗族自治州自然资源和规划局, 湖北 恩施 445000)

植被作为陆地生态系统的重要组成部分，是连接大气、水体和土壤等生态要素的自然纽带[1]，其动态变化能够表征地球表面的能量流动、碳循环和人类活动，在全球变化中起到“指示器”的作用[2-3]。同时，植被是生态环境变化的承受者，很大程度上代表了生态环境总体状况[4-5]。归一化植被指数(Normalized Difference Vegetation Index，NDVI)与叶面积指数、生物量及净初级生产力等有较好的相关性[1,6]，能够很好的反映植被生长状况和植被覆盖变化，是大尺度植被研究的重要指标。

国内外学者利用NDVI数据分析了不同时空尺度下植被时空特征及其对气候变化的响应。全球变暖使得北半球中高纬度地区植被覆盖有增加趋势[7]，也有学者认为北半球中高纬度地区植被活动增强与温度有关，而降水减少导致了南半球干旱地区NDVI呈减少趋势[8]。有研究表明欧亚大陆植被生长季NDVI显著增加[9]，而在不同数据源、不同时间尺度下我国植被NDVI整体也表现出增加趋势，并且地区差异明显，其中国家重点生态治理的青藏高原[10]、黄土高原[11]和三江源等[3]地区显著增加。植被生长主要受到气候变化和人类活动的影响，其中气候变化是影响植被年际活动变化的重要影响因素[12-13]。植被对气候变化响应研究聚焦于与水热因子关系，其关系具有高度的空间异质性，如青藏高原[10]和北方地区[14]植被生长对气温变化更敏感，而在宁夏[15]降水和湿润指数的影响最为显著。因此，植被NDVI的变化特征及水热因子的决定力存在区域差异。

武陵山区位于我国中南部，区内森林覆盖率高、自然资源丰富，《全国主体功能区规划》将其规划为限制开发区，基本功能为生物多样性及水土保持。由于多山、地形起伏大，恶劣的自然条件和人类不合理的开发活动使得该区面临着生物多样性受到威胁、土壤侵蚀严重、地质灾害多等问题，同时作为亚热带森林生态系统核心区、长江流域重要的水源涵养区和生态屏障，武陵山区具有重要的生态保护和科学研究价值。然而，目前鲜有武陵山区植被覆盖的相关研究，随着该区区域性的扶贫开发，以及退耕还林(草)、石漠化治理等工程的实施，其植被生长和分布状况必然发生变化，该区植被覆盖的变化特征及对气候变化的响应特征有待厘清。基于此，本研究旨在利用MODIS NDVI数据揭示武陵山区近20 a植被生长季NDVI的时空变化特征，分析影响植被生长的关键气候因子及其滞后性，以期为气候变化背景下该区的植被保护、生态治理和合理开发提供科学参考。

1 研究区域及数据

1.1 研究区概况

武陵山区地跨湖南省、湖北省、重庆市和贵州省，包含71个县(市、区)，总面积约17.18万 km2，地理范围为25°52′—31°24′ N，107°04′—112°02′ E，是云贵高原的东部延伸地带，在我国地形阶梯分级中，属二、三级阶梯过渡带，地势西北高、东南低(图1A)，平均海拔1 000 m左右。其位于季风区，气候属亚热带向暖温带过渡类型的山区类型气候，夏凉冬冷。区内土地利用多为林地和耕地(图1B)，森林覆盖率平均达56%，植被类型以常绿阔叶林和落叶阔叶林为主。

图1 武陵山区地形及土地利用

1.2 数据来源及处理

MODIS NDVI数据采用NASA发布的2000—2019年MOD13Q1植被指数数据集，时间和空间分辨率分别为16 d和250 m。利用MRT软件对数据进行格式转换、重投影等操作，采用最大值合成法消除异常值影响并剔除NDVI值小于0.1的无植被区域，生成NDVI月值数据集，进一步选取植被生长季(4—10月)NDVI平均值作为年NDVI平均值，最后利用研究区边界裁剪获得NDVI逐年、逐月数据。

气象数据(气温、降水)为国家气象科学数据中心提供的2000—2019年中国地面气候资料月值数据集。选取研究区内及周围共49个气象站点，其中44个用于插值，5个用于验证。气象数据插值采用Anusplin软件，以经纬度为自变量，并引入高程作为协变量来提高插值精度。交叉验证法用于评估插值结果，选用平均绝对误差(MAE)、平均相对误差(MRE)和均方根误差(RMSE)作为检验标准[16-17]。验证结果见表1，气温插值精度高于降水，两者的插值精度均能够满足研究需求。地表覆盖数据GlobaLand30(2010年)来源于全国地理信息资源目录服务系统，在该分类体系中，研究区主要包括耕地、林地、草地、灌木地、湿地、水体和人造地表等7种类型(图1B)。此外， DEM数据(分辨率90 m)来源于地理空间数据云(表2)。所有数据定义相同的空间参考(UTM WGS 84 49N)和分辨率(250 m)，数据处理与分析通过ArcGIS和Python完成。

表1 气象插值交叉验证结果

表2 研究数据来源

2 研究方法

2.1 稳定性分析

变异系数可用于表示NDVI的相对波动程度，其值越大，说明植被受到的干扰多、不稳定，反之，植被处于相对稳定状态。其计算公式为：

(1)

2.2 趋势分析与检验

Sen趋势分析的优点是不需要数据服从特定的分布，并且不受异常值的干扰，对测量误差或离散数据具有较强的规避能力[18-19]。其计算公式为：

2000≤i

(2)

式中:β为NDVI变化趋势;β>0说明植被覆盖呈上升趋势,β< 0则为下降趋势。变化趋势的显著性通过Mann-Kendall统计检验法判断,相关公式如下:

检验统计量

(3)

符号函数

(4)

标准化统计量

(5)

(6)

趋势检验的方法[20]：零假设H0，β=0，当|Z|>Z1-α/2，拒绝零假设。Z1-α/2为标准正态方差，α为置信度水平。

2.3 Hurst指数

基于重标极差(R/S)分析方法的Hurst指数是定量描述时间序列信息长期依赖性的有效方法[3, 19]，它能够依据长时间序列数据的过去趋势预测未来发展趋势。其基本原理是：

考虑一个NDVI时间序列{NDVI(τ)}，对于任意正整数，计算公式如下所示。

均值序列：

(7)

累计离差：

1≤t≤τ

(8)

极差：

(9)

标准差：

(τ=1,2,…，n)

τ=1,2,…，n

(10)

若存在R/S∝τH，说明时间序列存在Hurst现象，H称为Hurst指数，H值在双对数坐标系(lnτ，ln(R/S))中用最小二乘法拟合获得。

3 结果与分析

3.1 植被生长季NDVI变化特征分析

3.1.1 植被生长季NDVI的时间变化特征 武陵山区近20 a植被生长季NDVI的变化范围为0.64～0.76(图2)，多年NDVI均值是0.71，整体呈现波动增加态势，增速为0.52%/a(p<0.01)，植被条件不断改善。NDVI的年际变化存在2001年(0.64)，2006年(0.66)和2014年(0.71)3个明显谷值，其中2001年为研究时段最低值，而2019年(0.76)为最高值。

图2 武陵山区植被生长季NDVI年际变化趋势

3.1.2 植被生长季NDVI的空间变化特征 2000—2019年，武陵山区各地多年植被生长季NDVI均值的范围为0.1～0.86，其中NDVI>0.6的地区比例达到94.65%，表明该区植被覆盖整体较好。但从图3A可以看出，区内NDVI空间差异明显，空间格局总体表现为“东北高、东南和西部低”。高值区主要是湖北境内的鹤峰、五峰和长阳等地，该区气候湿润，植被长势好、覆盖度高，同时建有国家级自然保护区，植被不易受到侵扰。而区域的东南和西部多以耕地为主(图1B)，特别是西部喀斯特地貌分布广泛，存在严重的石漠化问题，致使植被稀疏、NDVI低值集中。

武陵山区植被生长季NDVI变异系数的变化范围为0.01～0.62(图3B)，均值是0.07。变异系数低于0.15的区域比例达到99.54%，表明研究时段武陵山区植被生长季NDVI空间变化整体处于较稳定状态。区域东北部海拔高、以林地为主，受到的干扰少，变异系数低值聚集。区域西部，尤其是贵州境内，可能是进行退耕还林(草)、封山育林等生态工程治理石漠化，致使NDVI变化相对剧烈、高值集中，但植被生态系统趋向良性循环。东南部整体变化较稳定，由于城镇化建设对植被的破坏，少数不稳定区域零散分布于城镇周围。

图3 武陵山区多年生长季NDVI及变异系数空间分布

NDVI的垂直变化特征可以分为3个阶段(图4)：第一段(≤500 m)，低海拔地区人类活动剧烈，用地多为耕地，植被覆盖低，随着海拔的升高，植被类型变得多样，NDVI急剧上升；第二段(500～2 100 m)，人类活动的影响减弱，耕地比例减少而林草地增多，NDVI平稳上升；第三段(>2 100 m)，高海拔地区虽然林地占比多，但植被类型单一，同时生境变恶劣、植被生长受到限制，NDVI下降。总体而言，NDVI随着海拔上升表现出先增后减的特征。

图4 不同海拔NDVI的分布

3.1.3 植被生长季NDVI空间变化的趋势 图5将Sen趋势与Mann-Kendall检验结合，分析武陵山区植被生长季NDVI变化趋势的空间分布：趋势分析β值介于0.033～0.022/a，均值是0.004/a，其中94.79%的区域植被覆盖呈上升趋势(β>0)，说明武陵山区植被生长季NDVI总体趋势为微弱上升，与3.1.1线性变化趋势结论一致。

图5 NDVI变化趋势及其显著性的空间分布

表3呈现了NDVI变化趋势的分级标准及结果，可以看到植被覆盖虽然以上升趋势为主，但多为不显著上升(80.58%)，而下降趋势也主要表现为不显著(4.96%)。从空间分布上看，不显著上升广泛分布于整个区域，显著上升区域集中于研究区中、西部，极显著上升和不显著下降区域分布均较为零散，其中不显著下降区域主要在城镇周围及低植被覆盖地区，而显著下降和极显著下降区域较少。结合高程来看，下降趋势主要发生在海拔1 000 m以下的地区(62.39%)，其次是1 000～2 000 m(37.29%)。

表3 NDVI变化趋势统计

为分析不同用地类型NDVI的变化特征，统计其变异系数及变化趋势(表4)。各用地类型变异系数均较小，从小到大排序为林地、草地、耕地、灌木地和湿地：林地较稳定、抗干扰力强，NDVI变化最小；而湿地易受干扰、NDVI变化相对大，但整体上NDVI处于稳定状态。退耕还林(草)等生态工程实施以后，不同用地类型NDVI都有上升趋势，最多的是林地(95.82%)，然后依次是耕地(94.16%)、草地(93.86%)、灌木地(90.22%)和湿地(76.47%)，但多为不显著上升。

表4 不同用地类型的变异系数及变化趋势统计 %

3.1.4 植被生长季NDVI空间变化的可持续性 由Hurst指数空间分布图可知(图6)，武陵山区生长季NDVI的Hurst指数介于0.09～0.98，均值为0.41。Hurst指数大于0.5表明NDVI变化趋势是持续的上升或下降，小于0.5则是反持续，未来趋势与过去趋势相反。参考相关研究[21]，进一步将变化趋势的持续性划分为强反持续(0

图6 Hurst指数空间分布

图7 Hurst指数分级及统计

将空间变化趋势与可持续性图层叠加，分析植被覆盖的未来变化特征。由图8可以看出，武陵山区生长季植被NDVI未来趋势主要表现为下降(81.55%)，其中上升—下降占比80.27%，意味着80.27%的区域上升趋势将转变为下降，全区广泛分布。未来表现为上升趋势的区域占比18.45%，持续上升的区域多在湖南境内，而由下降转变为上升的主要为重庆的丰都、彭水、武隆及贵州的江口、余庆、湄潭、石阡等地。值得注意是，1.28%的区域过去和未来都将处于下降状态，主要分布于湖南的怀化、洞口、武冈、邵阳和湖北的恩施等地，多集中在城镇用地周边，城镇化建设大量侵占耕地和林草地导致NDVI持续下降。

图8 NDVI的未来变化特征

总体而言，武陵山区生长季植被覆盖虽然以上升为主，生态保护工作取得了一定成效，但未来变化特征却主要表现为下降。这些区域多为海拔1 000 m以下(76.78%)，坡度相对平缓(<15°占比78.58%)，容易受到人类活动的影响，可能是未来下降趋势的原因。需在巩固现有保护成果基础上，持续关注重点区域，防止植被覆盖下降趋势的发生。

3.2 植被生长季NDVI对气候因子的响应

3.2.1 气温和降水时空变化格局 2000—2019年，武陵山区植被生长季平均气温变化范围为-0.03～0.07℃/a(图9A)，平均值是0.017℃/a，呈不显著增温趋势。除中部和东北部少数地区气温呈下降趋势外，大部分地区气温增加，北部和西南部相对明显，空间差异显著。植被生长季累积降水量变化介于-21.56～17.31 mm/a(图9B)，平均值是1.70 mm/a，呈不显著增加趋势。降水变化趋势也存在明显的空间差异，东北、西南部降水减少，湖南中部降水增加显著。总的来说，不同用地类型水、热均有所增加，武陵山区植被生长季气候趋向暖湿。

图9 2000-2019年武陵山区生长季气温、降水变化趋势的空间分布

3.2.2 植被生长季NDVI与气温、降水的相关性 采用偏相关性分析，定量识别植被生长季NDVI对气候因子的响应。NDVI与气温的偏相关系数介于-0.84～0.89(图10A)，均值是0.07，大部分区域(62.75%)表现为正相关。NDVI与降水的偏相关系数在-0.90～0.87(图10B)，均值为0.1，正相关区域占比66.64%。NDVI与气温和降水均以正相关为主，且与降水的关系更密切，区域差异显著：NDVI与气温呈负相关主要是重庆的秀山、酉阳、黔江与贵州的沿河等地，在湖南的会同、芷江、靖州、通道、洪江和辰溪等地明显正相关；NDVI与降水在贵州的务川、德江、沿河、印江和重庆的酉阳、秀山等地以负相关为主，明显正相关主要集中于湖南的怀化、溆浦、洪江、辰溪、麻阳和重庆的丰都等地。

图10 植被生长季NDVI与气温、降水偏相关性

根据NDVI与气候因子偏相关系数正负及大小，将偏相关系数分为8个类别(图10C，D)，分析区域NDVI主要气候因子，结果表明：

RT+>RP+，RT+>|RP-|，|RT-|>RP+和RT-RP+主要是湖北北部、湖南东北和西南部及贵州西南部分地区，NDVI与气温、降水均为正相关，占整个区域的21.44%，这些区域气温呈上升趋势(图9A)，少数地区(湖北、贵州境内)降水呈下降趋势(图9B)，适当增温能够促进植被光合作用，气温是主要气候因子。湖北北部和贵州西部等地，RT+>|RP-|，NDVI与气温和降水分别为正、负相关，比例为5.72%，多为高海拔地区，气温低而降水充足，气温对植被生长的影响高于降水。重庆中部、贵州东部及湖北东部，|RT-|>RP+，NDVI与气温和降水分别为负、正相关，比例为6.02%，该区气温高，进一步升温会增强植被呼吸作用从而加速营养物质消耗，同时加快地面水分蒸发不利于干物质积累，为主要限制因子。重庆南部和湖北东北部，RT-

降水主导包括RT+≤RP+，RT+≤|RP-|，|RT-|≤RP+和RT-≥RP-。30.85%的区域RT+≤RP+，包括湖南大部分地区(尤其是中部)及重庆西北部，该区降水量较少，气温与降水均呈增加趋势，水热条件改善，特别是充足的水分利于植被生长。RT+≤|RP-|零散分布于贵州中部和湖北中部，占比4.99%，水分增加对植被生长起抑制作用，降水是主要气候因子。重庆北部和湖南东南部，|RT-|≤RP+，比例为8.30%，水分增加促进植被生长，降水是主要影响因素。贵州中部和湖北东部，RT-≥RP-，NDVI与气温和降水均为负相关，比例为9.85%，水热因子对NDVI具有负效应且主要受水分影响。

综合来看，气温和降水作为武陵山区生长季NDVI主导气候因子的比例分别为46.01%，53.99%，整体上降水的影响范围更大。此外，气温或降水主导的区域与NDVI的关系以正相关为主(63%)，两者的比例分别为27.16%和35.84%。

3.2.3 植被生长季NDVI对气候响应的时滞效应 植被生长对气候响应具有一定时滞效应，滞后程度对于不同气候因子、不同地区存在差异。计算植被生长季NDVI与气候因子当月(4—10月)、前1月(3—9月)、前2月(2—8月)的偏相关系数，分析植被对气候响应的时滞效应。由图11可知，月尺度上，武陵山区植被生长季NDVI整体上对气温的响应未见明显滞后性，占比达到98.50%，对降水则以滞后1个月为主，比例为63.79%。但仍有1.13%的区域NDVI对气温响应存在1个月滞后期，主要位于湖南的城步、隆回等地，响应滞后2个月的区域则零星分布。同时湖北东北部及重庆的丰都等地NDVI对降水响应没有滞后性，比例为7.23%；高于1个月响应滞后期的区域主要集中在湖北利川、恩施、来凤、咸丰及重庆的石柱、黔江等地，占整个区域的28.98%。

图11 NDVI对气候响应滞后性的空间分布

将NDVI与当月、前1月、前2月气候因子的偏相关系数分类，得到0～2个月滞后期主控气候因子的空间分布(图12)。当月，降水为NDVI主控气候因子，类型以RT+≤RP+为主(79.01%)，整个区域广泛分布，其次是|RT-|≤RP+(20.27%)，主要包括贵州的务川、湄潭、道真和重庆的彭水、酉阳及湖南的城步、通道等地。前1月，NDVI与气温、降水均以正相关为主，58.28%的区域气温主控(RT+>RP+占56.63%)，主要位于研究区中部和南部，而降水主控比例为41.72%(RT+≤RP+占39.91%)，集中在区域北部和东南部。前2月，NDVI主控气候因子为气温，类型主要是RT+≤|RP-|(52.90%)，包括北部的大部分区域，其次是RT+>RP+(37.72%)，主要分布在湖南，尤其是东南部。总体而言，0～2个月滞后期主控气候因子依次是降水、气温和气温，降水的影响随着时间推移逐渐减弱，而气温逐渐增强。

图12 NDVI对气候响应滞后性主控因子的空间分布及统计

4 讨论与结论

4.1 讨 论

植被在区域生态保护和恢复中扮演着重要角色，武陵山区生态环境脆弱，复杂的地形地貌和多变的自然环境，对植被覆盖产生了巨大影响，因此，对该区植被进行长时间监测和分析具有重要的现实意义。随着退耕还林、天然林保护、石漠化治理和水土保持等重点生态修复工程的实施，武陵山区生长季植被NDVI整体呈现波动增加态势，植被条件不断改善，与相关研究结果一致[22-23]。空间上，绝大多数区域(94.79%)植被覆盖呈上升趋势，下降区域主要是城镇周围及低植被覆盖地区。2001年、2006年和2014年3个明显NDVI谷值可能与该年份水热不匹配有关：2001年气温低、降水少，干冷不利于植被生长；2006年气温多年最高、蒸发旺盛，但降水少，缺水抑制了植被生长；2014年气温低、降水过多，达到1 319.42 mm，易引起植被根系缺氧。

植被受地形因素影响在垂向结构上呈现明显的空间分异特征[24]，武陵山区NDVI随海拔升高表现为先增后减，与植被NPP随海拔变化节点大致相同[25]。低海拔地区易受城镇化建设和农业生产等人类活动的影响，是植被退化易发生区，研究区的退化区域主要集中于2 000 m以下，特别是1 000 m以下。

气温和降水通过影响有效积温和可利用水分来调控植物光合作用、呼吸作用及土壤有机碳分解等进而影响植物的生长和分布[26]，是主要的非生物因素。受季风气候和地形分异的影响，武陵山区气温和降水时空差异明显，研究时段内气候总体趋向暖湿，有利于植被生长。生长季尺度，近20 a武陵山区植被生长季NDVI与气温、降水的关系均整体表现为正相关。NDVI与气温呈较显著正相关主要是西南部(湖南西部)和高海拔的湖北北部，这些区域气温相对低，气温呈升高趋势促进了植被生长。研究区中部和西北部多为低海拔地区，气温高而降水不足，水分成为制约植被生长的主要因素，而降水在研究时段呈明显增加趋势，因此NDVI与降水具有较显著的正相关关系。武陵山区NDVI最大值多为7月，而7月是最高气温月份，此外降水多集中在5—9月，特别是6月，月尺度而言，武陵山区对气温未见明显滞后性，对降水以滞后一个月为主。

4.2 结 论

(1) 植被生长季NDVI的变化范围为0.64～0.76，多年NDVI均值是0.71，整体呈现波动增加态势，增速为0.52%/a(p<0.01)，植被条件不断改善。

(2) 武陵山区各地多年植被生长季NDVI均值的介于0.1～0.86，94.65%的区域NDVI>0.6，植被覆盖整体较好。NDVI空间差异明显，空间格局呈现“东北高、东南和西部低”的特点。空间变化整体处于较稳定状态，东北部低值聚集，高值主要位于西部。NDVI的垂直变化特征分为急剧上升、平稳上升和下降3个阶段，随海拔升高总体表现出先增后减。

(3) 绝大多数区域(94.79%)植被生长季NDVI呈上升趋势，但多不显著，而呈下降趋势的区域也以不显著为主。下降主要发生在海拔1 000 m以下，其次是1 000～2 000 m。NDVI空间变化整体趋向反持续性，未来大部分区域NDVI变化趋势将会发生改变，但改变幅度比较小。未来变化特征主要表现为下降。

(4) 生长季尺度，武陵山区气候趋向暖湿，NDVI与气温、降水均以正相关为主。气温和降水对区域NDVI的影响有明显的空间差异性，整体上看降水的影响范围更大。

(5) 月尺度，NDVI整体上对气温的响应未见明显滞后性，而对降水的响应大都滞后一个月。0～2个月滞后期主控气候因子依次是降水、气温和气温，其中降水的影响随着时间推移逐渐减弱，而气温的影响逐渐增强。