全球气候变化与陆地生态系统的关系是全球变化科学研究中的核心内容之一^[1]。IPCCAR5显示，1951—2012年全球平均气温上升0.72 ℃(0.49—0.89 ℃)，其中1983—2012年成为北半球过去800年来最暖的30年^[2]，气候的剧烈变化无疑会对陆地生态系统造成严重影响。植被作为陆地生态系统的主体，既是气候变化的敏感指示器^{[3, 4]}，同时也通过与大气之间的能量、水分和物质交换反作用于气候^{[5, 6]}。目前，植被与气候变化的关系研究受到了国内外学者的广泛关注与研究，其中Myneni等^[7]和Tucker等^[8]基于NDVI数据研究了北半球高纬度地区植被变化情况，结果表明：在气候变暖背景下，北半球植被活动显著增强，而同样的植被变化趋势在中国也被检测出来^[9]。相反，植被覆盖的变化引起了地表下垫面性质及其他生物化学过程的改变，从而对区域气候产生影响^[10]，已有研究表明，我国退耕还林区植被覆盖的增加导致昼间气温下降，夜间气温升高，从而造成昼夜温差减小^[11]。

在全球变暖背景下，近50年(1960—2009)我国地表平均气温升高了1.38 ℃，远高于全球或北半球同期平均增温速率^[12]。显著的气候变暖势必会对陆地生态系统造成严重影响，尤其表现在植被覆盖的变化^{[13, 14]}。我国学者针对中国东部^[15]、西北地区^{[16, 17]}、青藏高原^[18]、黄土高原^[19]、内蒙古^[20]以及三江源^[21]等地区植被覆盖变化做了大量研究工作。以上研究成果从区域尺度很好地揭示了植被覆盖变化及其与气候因子的关系，并得出了有意义的结论。然而，关于全国尺度植被变化的宏观格局认识仍较为缺乏，尤其是长时间尺度的时空变化特征更是鲜有报道。此外，另有研究表明，欧亚大陆的植被覆盖增加趋势在20世纪80—90年代期间逐渐变缓并趋于停滞、甚至部分地区呈下降趋势^{[22, 23]}，如1996年前后东亚地区植被覆盖变化趋势发生翻转，由增加趋势变为下降趋势^[24]。植被对气候变化的这种响应引起了学界的广泛关注，其中朴世龙^[10]和张镱锂^[25]分别对欧亚大陆和青藏高原部分地区的植被变化趋势进行了探讨，二者均检测出了植被变化的转折现象。综上可以看出，上述研究成果加深了我们对植被覆盖变化的认识水平。植被覆盖状况作为重要的生态指示因子，有必要进一步加强其变化特征及归因的研究，以全面认识植被覆盖的时空变化特征及其影响因素，保障生态安全。

我国位于东亚季风区，对全球气候变化的响应极为敏感，全球气候变化已对我国植被覆盖造成显著影响^[26]。与此同时，近年来我国植被覆盖经历了最强烈的人类活动干扰，如退耕还林还草工程、三北防护林工程的实施。然而，关于我国植被覆盖变化的宏观特征及其长期趋势尚不明晰。基于上述认识，本文利用GIMMS-NDVI、MODIS-NDVI数据和气象数据，辅以分段线性回归模型、趋势分析以及相关分析等手段，分析了1982—2012年我国植被生长季NDVI的时空变化特征及驱动因素。综合评价我国近30 a植被NDVI的时空变化特征不仅能够正确认识气候变化对我国植被生长的影响，并且能够正确评价我国退耕还林还草等生态恢复工程的实施效果，为进一步生态恢复工程的实施及生态文明的建设提供基础数据。

本文所用GIMMS NDVI数据为美国马里兰大学GLCF(Global Land Cover Facility)研究组生产的15d最大合成数据^[27]，其空间分辨率为8 km × 8 km，该数据是目前时间序列最长的NDVI数据集，时间跨度为1982—2006年。在数据制备过程中已经过辐射校正和几何校正、除云、除火山气溶胶等预处理，并且利用沙漠控制点进行验证，保证了数据质量，与其他NDVI数据相比，其误差小，精度较高，已被广泛应用于全球及区域植被变化研究^[12]。在数据处理过程中，本文采用最大合成法(MVC)进一步消除了云、大气等异常值的影响^[28]。

本文采用另一遥感数据为月最大合成MODIS NDVI数据，空间分辨率为1 km×1 km，时间跨度为2000—2012年，取自美国国家航空航天局(NASA)的EOS\\MODIS数据产品(http://e4ftl01.cr.usgs.gov)。与GIMMS NDVI相比，该数据具有更高的分辨率和更稳定的影像质量，在监测和评价陆地植被变化的能力要优于前者^[27]。在数据处理过程中，采用MODIS Reprojection Tools(MRT)软件对下载的原始数据进行格式转换、投影转换、数据拼接、裁切等预处理。为与GIMMS NDVI数据格式保持一致，本文将MODIS NDVI数据重采样为8 km × 8 km分辨率的影像。另外，为消除低植被覆盖区的影响，本文将多年平均NDVI值大于0.1的区域定义为植被区域，并进一步定义每年的4—10月为该年份的生长季。

由于两种NDVI数据采用了不同的传感器，因此在数据融合之前需要对两种数据进行一致性检验^[19]。本文采用的GIMMS NDVI时间范围为1982—2006年，MODIS NDVI数据的时间范围为2000—2012年。首先将两种数据源的重合年份进行时序分析，2000—2002年之间具有较好的一致性，因此文中采用以上3个年份的生长季逐月数据进行相关性分析，并求得相关系数分别为0.984、0.989和0.998，均通过了0.001的置信度检验。事实上，以上相关分析存在年内生长周期的影响，为消除此影响，本文进一步对生长季的逐月数据进行7步长滑动平均处理，将处理后的数据重新进行相关性分析，发现2001和2002年相关性系数分别为0.45和0.99，仅2002年通过了0.001置信度水平的检验。因此，本文采用2002年的两种NDVI数据进行所有像元的拟合，即MODIS NDVI(MNDVI)作为自变量，GIMMS NDVI(GNDVI)作为因变量，获取二者回归方程，从而实现两种数据的融合。根据计算，拟合方程为：GNDVI = 0.7153×MNDVI+0.1171(样本量为143150，R²=0.7281，P<0.001)。需说明的是，由于2003—2006年的GIMMS NDVI数据与MODIS NDVI和SPOT NDVI均表现出较大的差异，因此，本文采用上述拟合方程对2003—2012年的GIMMS NDVI数据进行插补。

为保证资料的完整性和连续性，建立均一、稳定的气候序列，本研究共选取研究区578个气象台站的月平均气温和月降水量数据，资料时间跨度为1982—2012年，来源于中国气象科学数据共享服务网(http://cdc.cma.gov.cn)。所选站点均经过了严格的质量检查和控制，保证了数据的质量。

一元线性回归分析能够模拟每个栅格的变化趋势，该方法是指在一定时间内，采用最小二乘法逐像元拟合年均NDVI的斜率，用以综合反映植被覆盖的时空演变特征^[29]：

点击浏览公式

为了检测1982—2012年我国植被覆盖变化趋势的持续性，本文采用分段线性回归模型进行分析^[30]。由于该模型能够较好地提取长时间序列数据变化的转折点，从而有效弥补普通线性回归模型无法真实反映长时间序列变化趋势的缺陷，已被广泛应用于植被与气候变化趋势分析^[31]：

点击浏览公式

地理系统是一个多要素构成的复杂巨系统，系统中任何一个要素的变化必然影响到其他要素的变化，而偏相关分析可以有效解决这个问题。它是指研究某一要素与另一要素的相关程度时，将其他要素的影响视为常数的方法^[32]：

点击浏览公式

1982—2012年中国植被NDVI呈缓慢增加趋势，增速为0.2%/10a，未通过显著性检验(图 1)。通过分段线性分析，我们发现近30年中国植被NDVI变化主要分为两个阶段：(1)1982—1997年呈显著增加趋势，增速达1.2%/10a(P<0.05)；(2)1997—2012年增加趋势小于前一时期，为0.6%/10a(P<0.05)。1997—1998年植被NDVI迅速下降，将两年NDVI进行相减，统计得出有61%的像元呈下降趋势，其原因可能是该时期强厄尔尼诺事件导致大气环境的改变，进而对植被覆盖造成影响。由于1997年NDVI值对后一阶段植被变化趋势存在显著影响，为此，进一步计算了1998—2012年的植被NDVI变化趋势，结果显示该时段植被覆盖呈极显著增加趋势，增速为1.0%/10a(P<0.01)。以上分析表明，近30年我国植被NDVI呈增加趋势，但是不同阶段增速有所差异。

图 1 1982—2012年中国生长季NDVI年际变化 Fig.1 Inter-annual variations of growing season NDVI over China during 1982—2012

图选项

空间平均值可以表征植被NDVI整体变化趋势，但由于存在不同区域变化趋势相反，进而相互抵消的情况，所以不能很好地描述不同区域的变化特征。因此，基于一元回归模型，在像元尺度上分析了近30年植被NDVI变化趋势，并对计算结果进行显著性检验。1982—2012年植被NDVI呈增加和减少趋势的面积分别占53.79%和46.21%，其中极显著增加区域(P<0.01)占29.32%，主要分布在黄土高原、青藏高原中西部、祁连山西段、新疆准葛尔盆地以及云南的部分地区；而极显著减少地区(P<0.01)占23.89%，主要分布在东北地区的大、小兴安岭和长白山、新疆的阿尔泰山和天山、青藏高原的东部以及东南沿海等地区(图 2a—b)。

依据分段结果，进一步分析了1982—1997年(时期1)和1997—2012年(时期2)的变化趋势，发现，两个时间段内植被NDVI变化的空间格局具有显著的差异。时期1内植被NDVI呈增加和减少趋势的面积分别占61.09%和38.91%，其中不显著增加的面积占45.90%，而极显著增加的面积仅占4.97%，且主要分布在40°N以北的东北平原和黄土高原北部地区(图 2c—d)。时期2内植被NDVI呈增加和减少趋势的面积分别占54.69%和45.31%，其中极显著增加(P<0.01)区域占27.89%，主要分布在110°E以西地区；极显著减少(P<0.01)区域占19.96%，主要分布在东北地区的大兴安岭和小兴安岭、新疆的天山和阿尔泰山、祁连山东段、黄河源以及东部沿海等地区(图 2e—f)。上述分析表明，近30年我国中高纬度山区植被活动呈下降趋势。

图 2 1982—2012年中国植被覆盖变化趋势（a,c,e）及显著性（b,d,f） Fig.2 Trend (a,c,e) and siginificant (b,d,f) of vegetation coverage change in China during 1982—2012

图选项

在空间上，我国植被NDVI整体呈东南、东北高、西北低的空间格局，高值区主要分布在东北地区、华北地区、华中和华南地区以及东南沿海等地，其原因是该地区主要植被类型为落叶针叶林、针叶阔叶混交林、温带落叶阔叶林、常绿阔叶林区域以及热带季雨林和雨林；低值区主要分布在内蒙古中西部、新疆大部分地区以及青藏高原中西部，这些地区主要为草原区、荒漠区以及青藏高原高寒植被区(图 3a)。由频度分布图可以看出(图 3b)，我国植被NDVI呈“双峰”结构，且NDVI值介于0.6—0.8的区域所占比例高达48.80%，远高于NDVI值介于0.1—0.3的区域(20.17%)，而值域在0.3—0.5的区域共占25.33%。这种空间分布形式说明我国高植被覆盖与低植被覆盖区域差异显著。

图 3 中国植被NDVI空间分布(a)及频度分布(b) Fig.3 Spatial distribution of NDVI (a) and its frequency distribution (b) in China

图选项

1982—2012年中国气温呈增加趋势，增速为0.04 ℃/a(-0.05—0.15 ℃/a)。在空间上，97.06%的站点呈增加趋势，其中增温趋势较大的地区主要分布在我国西北地区和西藏等地；而气温呈减少趋势的站点主要分布在我国东北地区的大兴安岭、新疆北部的天山、华北平原等地(图 4)。同期全国平均降水呈减少趋势(-0.54 mm/a)，其中呈减少和增加趋势的站点分别占53.11%和46.89%，且减少站点集中分布在100° E以东地区，包括东北地区、华北北部、长江中下游、西南地区以及东南沿海等地，而新疆的天山和阿尔泰山部分地区也呈减少趋势(图 4)。

图 4 1982—2012年中国年均气温和年降水量变化趋势空间分布图 Fig.4 Spatial distribution of change trend of annual mean temperature and annual cumulative precipitation in China during 1982—2012

图选项

为保证数据精度及结果的可靠性，文中植被NDVI对气候因子的响应采用逐站点偏相关的方法，具体步骤为：1)以气象站点为中心，逐月提取站点周边3 km×3 km范围内的NDVI平均值，作为该站点的NDVI值；2)基于NDVI提取结果，逐站点计算NDVI与气温和降水的偏相关系数。

在年尺度上，中国植被NDVI与同期气温并未呈现出较好的相关性，仅有38.06%的站点呈正相关关系，其中6.06%和9.69%的站点分别通过0.01和0.05的显著性检验，且主要分布在青藏高原和内蒙古中西部地区；而NDVI与降水量有53.98%的站点呈正相关关系，且主要分布在中国北方和青藏高原中西部地区，通过0.01和0.05的显著性检验的站点分别占5.88%和11.25%(图略)。

由于植被生长对水热条件的改变存在一定的滞后性，因此年尺度的相关性会掩盖部分信息，并不能很好地解释植被与气温、降水的相关性，因此本文进一步从月尺度分析了中国植被NDVI与气温、降水的相关性。植被NDVI与当月气温具有显著相关性的区域主要分布在东北地区、新疆北部、甘肃中东部以及华中、华南等地(图 5a)；NDVI与前1个月相关性空间格局与当月相似，明显不同的地区主要为河北、山西以及陕西北部等农牧交错带地区(图 5b)；而植被NDVI与前2个月气温的相关性有所减弱，存在较高相关性的区域主要分布在东北地区、北方农牧交错带、西藏南部、青海东部、四川西部以及东南沿海部分地区(图 5c)；NDVI与前3个月气温的相关性与前两个月相似，仅东南沿海地区相关性进一步减弱(图 5d)。综上可以看出，中国北方地区植被对气温具有较长的响应时间，尤其表现在东北地区和北方农牧交错带地区；而南方地区植被对气温则存在0—1个月的响应时间。

图 5 月NDVI与当月(a)、前1个月(b)、前2个月(c)、前3个月(d)气温偏相关系数空间分布 Fig.5 Spatial distribution of partial correlation coefficients between monthly NDVI and the current temperature (a),the preceding month temperature (b),temperature of two preceding months (c),and temperature of three preceding months (d)

图选项

植被NDVI与当月降水具有显著正相关的区域主要分布东北地区、内蒙古中部、河北、山西、陕西、新疆北部、甘肃和青海东部、西藏南部以及云南等地区(图 6a)；全国大部分地区NDVI与前1个月降水呈显著正相关关系(图 6b)。值得注意的是，除内蒙古中西部、新疆、青海北部等地区之外，中国北方大部分地区与前2个月和前3个月的降水均呈负相关关系(图 6c—d)，而华中、华南以及东南地区则与同期降水呈显著正相关关系，尤其是NDVI与前3个月降水关系更为显著。由此说明，我国北方地区植被对降水的响应主要有0—1个月的响应时间；而除云南外，我国南方地区植被对降水的响应时间存在1—3个月的响应时间，并且随着滞后时间的延长，相关性逐渐加强。

图 6 月NDVI与当月(a)、前1个月(b)、前2个月(c)、前3个月(d)降水量偏相关系数空间分布 Fig.6 Spatial distribution of partial correlation coefficients between monthly NDVI and the current temperature (a),the preceding month temperature (b),temperature of two preceding months (c),and temperature of three preceding months (d)

图选项

通过分析植被对气候因子的滞后时间，可以发现，植被对气温的响应不存在滞后时间的站点占40.83%，且主要分布在东北地区的大、小兴安岭和长白山、新疆的天山和阿尔泰山、甘肃中东部以及长江中下游地区；存在1个月滞后时间的站点占38.24%，主要分布在东北平原、山西、四川西部以及东南沿海等地区；存在2个月滞后时间的站点仅占4.67%，零星分布在内蒙古西部、广东和广西的南部以及云南部分地区；而存在3个月滞后时间的站点共占16.26%，主要分布在内蒙古中部、华北平原、西藏中西部以及云南等地(图 7a)。植被对降水的响应不存在滞后时间的站点占27.34%，主要分布在大兴安岭、河北和河南北部、山西、陕南地区、甘肃中东部以及四川东部等地区；存在1个月滞后时间的站点占53.81%，并广泛分布在东北、华北、西北、西南等地区；存在2个月和3个月滞后时间的站点分别占9.00%和9.86%，主要分布在东南沿海地区(图 7b)。

图 7 NDVI对气温(a)、降水(b)滞后时间分布 Fig.7 Spatial pattern of time lag of NDVI to temperature (a) and precipitation (b)

图选项

植被覆盖变化是气候变化和人类活动共同作用的结果，尤其是退耕还林还草等大面积生态恢复工程的实施，我国植被变化已深深记录了人类活动的印记。本文分析了我国退耕还林工程的历年累积造林面积与植被NDVI变化情况(图 8)，可以发现，自1999年退耕还林还草工程实施以来，我国造林面积持续上升，同时植被NDVI也表现出升高趋势，说明近年来生态恢复工程的实施对我国植被覆盖的增加起到了极大的促进作用，并取得了较好的效果。NDVI变化趋势并未与造林面积严格一致，而是存在年际波动，尤其是2008年之后NDVI波动较大，表明近年来气候波动对植被覆盖具有重要影响。同时，城市化、过度放牧等人类活动势必会对植被覆盖造成负面影响。以上分析表明，我国植被覆盖在气候变化和人类活动等多重因素共同驱动下，呈现出波动上升趋势。

图 8 退耕还林以来NDVI与历年累积造林面积变化 Fig.8 Change of NDVI and accumulated afforestation area in China

图选项

(1)1982—2012年中国植被NDVI呈缓慢增加趋势，增速为0.2%/10a，未通过显著性检验。近30年植被变化的阶段性特征显著，其中1982—1997年和1997—2012年均呈显著增加趋势，增速分别为1.2%/10a和0.6%/10a，均通过了显著水平0.05的检验。

(2)空间上，我国陕北黄土高原、西藏中西部以及新疆准格尔盆地等地植被NDVI呈显著增加趋势；而东北地区的大、小兴安岭、新疆北部的天山和阿尔泰山以及黄河源和秦巴山区等关键生态区和过渡带植被NDVI呈显著下降趋势，其中东北地区和新疆北部山区下降尤为显著，说明近年来我国中高纬度山区植被活动呈下降趋势。

(3)不同区域植被对气温、降水的响应存在差异，我国北方地区植被对气温具有较长的响应时间；而除云南外，我国南方地区植被对降水的响应时间存在1—3个月的响应时间，并且随着滞后时间的延长，相关性逐渐加强。

(4)我国植被覆盖增加是气候变化和人类活动共同驱动的结果，尤其是1999年之后人类活动影响逐渐加强，其中我国东北地区和新疆北部等植被覆盖的下降可能是由降水的减少导致，东南沿海地区则受城市化的影响更为显著。

已有研究指出，欧亚大陆的植被覆盖增加趋势在20世纪80—90年代期间逐渐变缓并趋于停滞、甚至部分地区呈下降趋势^{[22, 23]}，与本文研究结果存在差异，其原因可能是前人研究采用的是GIMMS NDVI，而本文对两种数据进行了融合处理，由此表明不同数据源监测结果存在较大差异，甚至是相反的结论。因此，应加强不同数据源监测结果的对比研究，以获得更加可靠的结论。近年来中高纬度山区、关键生态过渡带、关键生态区的植被活动在下降，而这些地区对我国生态安全具有举足轻重的作用。然而，目前对这些区域的本底知识仍然相对缺乏，有待进一步加强研究，以揭示其时空变化机理。同时，我国植被对气温和降水的滞后响应存在差异，主要表现为区域之间、不同地形之间的差异显著。本文从宏观尺度上揭示了植被对气候因子的响应格局，结合已有的区域尺度研究成果，从多尺度视角，全面认识植被对气候变化的响应规律。

本文存在的不确定性主要表现在：(1)我国山地众多，地形起伏较大，受空间分辨率的限制(8 km×8 km)，本文研究结果对区域细节的解释能力有限，但对于植被宏观变化格局及其对气候因子响应的空间格局仍具有重要的指示意义；(2)数据融合方法上采用了通用的一元回归分析，然而不同数据之间可能存在非线性的关系，因此在融合处理上会存在一定的不确定性。虽然气候变化是植被生长的关键控制要素，但是人类活动也对植被快速恢复起到重要作用，尤其是近年我国一系列生态恢复工程的实施。然而，在快速气候变化和强烈人类活动干预的背景下，如何定量区分气候变化与人类活动对植被变化的贡献率仍将是植被变化对气候响应研究的重要内容。本文仅仅分析了植被变化的长期趋势，并对其驱动机制进行归因，在全球变化背景下，极端气候事件的频度和强度在不断增多和增强^[33]，如高温热浪、暴雨、干旱等事件，而极端事件对植被覆盖变化及其空间格局势必会造成严重的影响。