相关研究^[1-2]表明，1951—2012年全球平均气温上升0.72 ℃，其中近30 a平均温度是过去1400 a来的最高水平。植被作为地球陆地生态系统的主体，植被与气候之间的动态耦合对陆地生态系统影响深远^[3-6]。同时，植被对气候变化的反应极其敏感，所以植被动态变化常被作为地球气候变化的生物指示器^[7-9]。研究发现，在北半球中高纬度地区，近几十年植被生长季延长和覆盖度的增加主要是由于气候变暖引起的，而近年来部分地区植被覆盖出现减缓甚至下降，很大可能是与气温上升的减缓和降水量的减少有关^[9-13]。随着对全球气候研究的不断深入，植被变化及其与气候变化的关系已成为全球变化研究的核心内容之一^[14-15]。

在全球变化背景下，中大尺度的植被覆盖变化受气候变化和人类活动的影响日趋显著。Sidi等^[16]对近18 a科拉马盆地植被变化进行研究，得出降水和温度等气候因素是影响植被生长的主要限制因子的结论，并发现在雨季和旱季，NDVI与气温和降水的相关性不同。陈晨等^[17]分析了1982—2015年黄河流域不同气候区植被的影响因素，表明降水对半干旱区NDVI影响最大，而气温对半湿润区NDVI影响最大。易浪等^[18]利用SPOT VGT数据分析黄土高原植被覆盖变化与气候和人类活动的关系，认为降水变化是影响植被覆盖变化的主要因素及人类活动对植被NDVI的正影响大于负影响。卢乔倩等^[19]利用1982—2012年植被生长季GIMMS NDVI 3g数据集对中国8种典型植被覆盖类型NDVI与气温和降水的时空响应特征进行研究，结果表明多数地区多年生长季平均NDVI与平均气温呈显著正相关。谢宝妮^[20]基于LTDR NDVI和同期气象插值数据对黄土高原1982—2014年不同气候时空变化以及气候和人类活动对植被覆盖变化的相对贡献进行研究，发现气候变化和人类活动共同影响黄土高原近30 a植被覆盖变化。综上分析，不同时空尺度下植被NDVI演变格局存在显著差异性，同时同一研究区不同历史时期植被NDVI变化的主导驱动因子也存在较大不同。基于此，本研究利用GIMMS NDVI和MODIS NDVI数据和气象站点数据等，引入地理探测器和重心模型等，从地理分区的视角定量分析和探讨了不同子区植被NDVI时空演变格局，并明确了各子区不同历史时期的植被NDVI演变主导驱动因子。

1 数据源及研究方法 1.1 数据源及预处理

归一化植被指数(normalized difference vegetation index，NDVI)数据集包括GIMMS NDVI(1981—2006年，时间分辨率为15 d，空间分辨率为8 km)和MODIS NDVI数据(2000—2019年，时间分辨率为16 d，空间分辨率为1 km)，上述两种数据分别来自资源环境科学与数据中心(https://www.resdc.cn/)和地理空间数据云(http://www.gscloud.cn/)，并利用最大值取值法获取了逐年最大NDVI数据，在较大程度上消除了云霾的影响以及作物的物候期差异导致的不确定性。考虑到两种指数数据集的来源不同，本研究基于2000—2006年重叠年份数据，利用回归分析法构建了数据校正模型：Y_MODISNDVI=1.077 7X_GIMMSNDVI+0.004 6(拟合样本2 899，R²=0.925，p＜0.001)。使用上述模型对2007—2019年的MODIS NDVI数据进行了校正，并利用ArcGIS 10.2 Resample工具将该数据集重采样为8 km。1980年和2020年1∶10万土地遥感监测数据来自中国科学院资源环境科学数据中心，空间分辨率为1 km，投影类型为Krasovsky_1940_Albers，包括6个1级类型，包括了25个2级类型，总体数据精度为97.15%。DEM数据来自SRTM3产品，由地理空间数据云免费下载，空间分辨率为90 m。中国土壤类型空间分布数据由《1∶100万中华人民共和国土壤图》数字化生成，该数据集包括12土纲，61个土类，227个亚类，由资源环境科学与数据中心(https://www.resdc.cn/)免费获取。气象站点数据来自中国气象数据网(http://data.cma.cn/)，包括中国824个地面基准气象站的20—20时降水量(0.1 mm)、日最低气温(0.1 ℃)、日最高气温(0.1 ℃)、日均温(0.1 ℃)、日照时数(0.1 h)等要素的逐日数据。利用C^#编程和ArcGIS 10.2的Cokriging插值方法得到1981—2019年的年均温、年降水量、年日照时数、年积温的栅格数据，其空间分辨率为8 km。

1.2 研究方法 1.2.1 相关性系数法

相关性分析主要用来反映要素之间的相关方向和相关程度，本研究采用基于像元的皮尔森相关系数，考虑p＜0.05和p＜0.01两种置信水平来探讨植被年NDVI对气温、降水的响应关系，其计算公式为^[21]：

$ {r_{x, y}} = \frac{{\sum\limits_{k = 1}^n {\left( {{x_k} - \bar x} \right)} \left( {{y_k} - \bar y} \right)}}{{\sqrt {\sum\limits_{k = 1}^n {{{\left( {{x_k} - \bar x} \right)}^2}} \sum\limits_{k = 1}^n {{{\left( {{y_k} - \bar y} \right)}^2}} } }} $

(1)

式中：r_{x, y}为x, y两变量的相关系数; x_k, y_k为x, y两变量的第k年取值; x, y分别为两变量的平均值; r的取值范围为[-1, 1]；r>0表示正相关; r < 0表示负相关, |r|值越大, 说明要素之间的相关性越高; r=0表示零相关。

1.2.2 重心模型

重心概念源于力学，是指区域空间上的某一点在前后左右各个方向上的力量对比保持相对平衡。在地理学研究中，“重心”可以表征地理要素的时空分布特征，广泛应用于经济、人口、生态、区域发展等领域。区域重心定义为：设某一区域由n个子区域构成，其中，第i个子区域的重心坐标为(x_i, y_i)，M_i为该子区域的属性值，则该区域重心的地理坐标表示为^[22]：

$ \bar x = \frac{{\sum\limits_{i = 1}^n {{M_i}} {x_i}}}{{\sum\limits_{i = 1}^n {{M_i}} }} $

(2)

$ \bar y = \frac{{\sum\limits_{i = 1}^n {{M_i}} {y_i}}}{{\sum\limits_{i = 1}^n {{M_i}} }} $

(3)

1.2.3 地理探测器

地理探测器是探测空间分异性，以及揭示其背后驱动因子的一种新的统计学方法^[23]。地理探测器q统计量，可用以度量空间分异性、探测解释因子、分析变量之间交互关系，已经在自然和社会科学多领域应用。

探测植被NDVI的空间分异性，以及探测某因子在多大程度上解释植被NDVI的空间分异。用q值度量，q的取值范围为[0, 1]，值越大说明植被的空间分异性越明显，反之则越弱。表达式为：

$ q = 1 - \frac{{\sum\limits_{h = 1}^L {{N_h}} \sigma _h^2}}{{N{\sigma ^2}}} = 1 - \frac{{{\rm{SSW}}}}{{{\rm{SST}}}} $

(4)

$ {\rm{SSW}} = \sum\limits_{h = 1}^L {{N_h}} \sigma _h^2 $

(5)

$ {\rm{SST}} = N{\sigma ^2} $

(6)

式中：h=1, …；L为变量Y或因子X的分层(strata)或分类; 为层h的单元数, N为全区的单元数; σ_h²是层h的Y值的方差, σ²是全区的Y值得方差。SSW为层内方差之和(within sum of squares)；SST是全区总方差(total sum of squares)。

2 中国自然地理分区

依据地形地貌、气候及自然地理特征，参照前人分区成果将研究区划分为六大区域(该分区方案来自资源环境科学与数据中心(https://www.resdc.cn/data.aspx? DATAID=276)，图 1)。

分别是东北、华北、西北、华东、西南以及中南地区。东北地区包括黑、吉、辽3个省区，为是我国重要的农业、工业基地。华北地区包括京、津、冀、晋、内蒙古5个省区，主要位于我国地势的第二级阶梯，是我国小麦的主产区。西北地区深居内陆，包括陕、甘、青、宁、新5个省区，主要为温带大陆性和高原气候，降水相对稀少。华东地区位于中国东部，包括沪、苏、浙、皖、闽、赣、鲁、台8个省区，地形以平原、丘陵为主，是我国经济发展的重心。西南地区包括川、黔、滇、渝、藏5个省区，而中南地区则包括豫、湘、鄂、粤、桂、琼、港、澳8个省区。

3 结果与分析 3.1 植被NDVI重心时空变化分析

为了分析和探讨研究区植被NDVI重心的空间变化规律，本文利用ArcGIS 10.2度量地理分布工具集的MeanCenter工具计算和分析了1981—2019年5 a尺度(5 a的植被NDVI均值)下全国和分区尺度上的重心迁移情况(图 2)。由图 2可知：①1981—2000年期间，东北地区植被NDVI重心呈现西南—东北—西南的迁移趋势，2006—2010年重心相比2001—2005年向西南方向迁移，说明该时间段内西南部植被NDVI的增速和增量大于东北部，而2011—2015年和2016—2019年相比2006—2010年出现向北迁移的回旋趋势。总体上近40 a东北地区植被NDVI重心呈向南移动的趋势，其原因在随着人类活动干扰强度的加剧，如三北防护林的建设及退耕还林还草的实施，该地区南部植被的恢复及增长速度和体量高于北部地区^[24-26]。②华北地区近40 a植被NDVI重心表现为向南迁移，其中2006—2010年→2011—2015年重心向北迁移，其他时间段重心均向南迁移，说明近40 a南部植被NDVI的增量要高于北部。总体上，华北地区植被NDVI重心向南迁移，其主要原因在气候变暖导致了该地区北部发生一定程度的干旱，不利于植被的恢复和生长，同时，随着南部地区植树造林工程的实施，其植被增量要高于北部地区^{[3, 17, 20]}。③西北地区植被NDVI重心在前3个时段表现为向西北方向迁移，分别是1981—1985年→1986—1990年、1986—1990年→1991—1995年和1991—1995年→1996—2000年。而在后4个时段均表现为向东南方向迁移，分别是1996—2000年→2001—2005年、2001—2005年→2006—2010年、2006—2010年→2011—2015年和2011—2015年→2016—2019年。综上分析，西北地区植被NDVI重心总体呈现向东南方向移动的趋势，其原因在受降水增加和人类活动(退牧还草及植树造林)的综合影响，东南部地区植被的恢复速率及增加总量要高于西北部地区(该地区干旱、少雨，植被很难恢复)，特别是十八大以来(2012年)，各项生态保护措施及工程的实施，一定程度上对于该地区的植被恢复产生了促进作用^{[6, 9]}。④1986—1990年华东地区植被NDVI重心相比1981—1985年向西北方向迁移。1986—2015年的5个时间段里，植被NDVI重心迁移轨迹呈现“西南—东北—东北—西南—南”的环形迁移趋势，而2016—2019年植被NDVI重心相比2011—2015年继续向南迁移。近40 a，华东地区植被NDVI重心向南迁移，其主要原因在于该地区南部地区多灌林地，随着各项生态环境保护措施的有效实施(2012年以来)，该地区植被的恢复和增加速率高于北部地区^{[1, 8]}。⑤中南地区植被NDVI重心迁移轨迹呈现“之”字形，1981—1985年→1986—1990年重心往东北方向迁移，1986—1990年→1991—1995年、2006—2010年→2011—2015年和2011—2015年→2016—2019年重心均是往西南方向发生迁移，而1991—1995年→1996—2000年、1996—2000年→2001—2005年和2001—2005年→2006—2010年3个时段重心均向西北方向迁移。综上分析，近40 a中南地区的植被NDVI重心向西移动，其主要原因在于随着退耕还林还草工程的有效实施(特别是2012年以来)，西部地区大量陡坡开垦的耕地转为林草地，极大了促进了该地区植被的恢复^{[3, 17]}。⑥1986—1990年和1991—1995年西南地区植被NDVI重心相比1981—1985年向西北方向迁移且1981—1985年→1986—1990年迁移幅度最大，1991—1995年→1996—2000年和2001—2005年→2006—2010年植被NDVI重心都是向东南方向迁移，而1996—2000年→2001—2005年重心向西南方向迁移，2011—2015年和2016—2019年植被重心相比2006—2010年向西北方向移动。近40 a西南地区植被NDVI重心向西南方向移动，原因在于随着全球气候变暖趋势的日益加剧，青藏高原地区的增温和增湿效应尤为明显，进而极大地促进了该地区植被的生长和恢复，从而导致了该地区西部和南部的植被增量和增速高于其他地区^{[8, 19]}。⑦1981—2019年中国植被NDVI重心主要集中于陕西省延安市及铜川市。1981—2019年，全国植被NDVI重心共发生了6次重心方向的改变，在南北方向上，除了1981—1985年→1986—1990年和2006—2010年→2011—2015年重心向北迁移，其他时间段重心均向南迁移。综上分析，近40 a中国植被NDVI重心向南迁移，表明南部地区植被NDVI的增量和增速要高于北部地区，这与区域气候条件的时空分布及变化相关^[10]。

3.2 气候单因子对植被NDVI变化的影响

全球气候变暖背景下，植被NDVI时空变化格局受气候变化的影响尤为显著。本研究在p＜0.01的置信水平下计算了近40 a植被NDVI与气温、降水和日照的相关系数(R)，进而分析了植被NDVI与气候单因子的相关性(图 3)。

由图 3可知，①植被NDVI与气温呈负相关(-1＜R＜0, p＜0.01)的区域占总面积的46.77%，主要分布于中国的北部地区，如西北地区的新疆中部，甘肃和青海省的西北部等，华北地区的内蒙古自治区以及东北平原的东部地区，其原因在于该地区气温的上升，加剧了植被的蒸腾和地表土壤水分的蒸发，进而加剧了区域干旱状况，不利于植被的生长^{[5, 8]}。植被NDVI与气温呈正相关(0＜R＜1, p＜0.01)的区域占总面积的53.23%，主要分布于黄土高原地区、四川盆地的中部、云贵高原以及长江中下游地区，该地区降水充沛，高温有助于植被的光合作用，加速植被的生长和恢复^[15]。②植被NDVI与降水呈正相关(0＜R＜1, p＜0.01)的区域占总面积的59.54%，呈正相关的地区分布范围较广，其中相关系数较高的区域包括华北地区的中部，西北地区的东北部，其原因在于该地区水热不均衡，地表水资源短缺，降水的增加有助于改善土壤水分状况，利于植被生长^[10]。植被NDVI与降水呈负相关(-1＜R＜0, p＜0.01)的区域占总面积的40.46%，主要分布于云贵高原和长江中下游地区等，其原因在于该地区降水密集，多暴雨，水土流失及滑坡、泥石流严重，极大破坏地表植被状况^[19]。③植被NDVI与日照呈正相关(0＜R＜1, p＜0.01)的区域占总面积的42%，其中显著正相关区域主要集中于新疆中部、甘肃和青海省的西北部以及云贵高原西部，其原因在于该地区海拔较高，充足的光照可以改善地表热环境，进而促进植被的生长^[14]。植被NDVI与日照呈负相关(-1＜R＜0, p＜0.01)的区域占总面积的58%，主要集中在中部地区(105°—115°E之间)，如贵州、山西和河南等地区，其原因在于该地区水热时空分布不均衡，过多的光照可增加区域温度，加大旱灾的风险^[20]。

3.3 不同历史时期—不同分区因子交互作用对植被NDVI变化的影响

相关研究表明，不同地理分区的植被生态系统演变的主导影响因子随着气候变化和人类活动作用强度的变化而变化。本研究利用地理探测器定量分析和探讨了1981年(1981—1983年平均植被状态)和2019年(2017—2019年平均植被状态)的驱动因子差异性(图 4)。由图 4可知：①1981年东北地区不同因子之间的交互作用对植被NDVI变化解释力依次为(前6组)：土地利用∩日照(0.657)＞高程∩日照(0.585)＞土壤类型∩日照(0.568)＞日照∩气温(0.566)＞坡向∩日照(0.528)＞积温∩日照(0.525)，日照与其他因子交互作用的解释力占据主导地位(图 4a₁)。而在2019年(图 4a₂)，解释力较高主导交互因子(前6组)分别为：土地利用∩坡向(0.473)＞土地利用∩气温(0.446)＞土地利用∩降水(0.356)＞土壤类型∩土地利用(0.347)＞土壤类型∩坡向(0.312)＞土地利用∩日照(0.305)，土地利用与其他因子交互作用的解释力却占据主导地位，而气候因子之间的交互作用则对植被NDVI的变化影响相对较弱。随着城市化的不断推进，人类活动成为影响地表生态系统演变的主导因素之一，其中土地利用变化是反映人类活动强度的重要指标之一，因此，土地利用变化在地表植被变化过程中发挥重要作用^{[4, 8]}。②1981年华北地区交互因子解释力最大的为日照∩气温(0.752)，最小的是坡向∩坡度(0.111)，同时，该时期土壤类型与其他因子之间的解释力都较大(图 4b₁)。2019年交互因子解释力最大的为土地利用∩降水(0.770)，最小的是坡向∩坡度(0.149)(图 4b₂)。华北地区多平原，地势平坦，地形因子对植被的时空分布和变化影响较弱，而气候因子和人类活动则是影响该地区植被生态系统演变过程的主导因素^{[15, 18]}。③1981年西北地区自然因子之间的交互作用对植被NDVI变化解释力依次为(前6组)：土壤类型∩降水(0.669)＞土地利用∩降水(0.642)＞降水∩气温(0.625)＞土壤类型∩日照(0.624)＞高程∩降水(0.621)＞日照∩气温(0.616)，降水与其他因子交互作用的解释力占据主导地位(图 4c₁)，其原因在于该地区干旱少雨，地表水资源是影响区域植被生长的主导因素，而降水是该地区水资源的重要补给^{[8, 16]}。2019年的主导交互因子解释力依次(前6组)为：土地利用∩降水(0.699)＞土地利用∩日照(0.681)＞土壤类型∩土地利用(0.648)＞土壤类型∩降水(0.635)＞土壤类型∩日照(0.589)＞高程∩降水(0.578)，其中土地利用∩降水之间的交互作用解释力最大(图 4c₂)，这与区域生态环境特征与人类活动相关^{[19, 21]}。④1981年华东地区交互因子解释力最大的为坡向∩降水(0.374)，最小的是积温∩气温(0.074)(图 5a₁)。2019年交互因子解释力最大的为土地利用∩坡向(0.412)，而最小的为积温∩气温(0.073)(图 5a₂)。由于该地区多为经济发达地区，城镇化水平较高，人类活动强度剧烈，已成为影响区域生态系统演变的主导因子^{[7, 14]}。⑤1981年中南地区各因子间的交互作用对植被NDVI变化解释力依次为(前6组)：土壤类型∩坡向(0.303)＞坡向∩降水(0.276)＞土壤类型∩日照(0.268)＞积温∩日照(0.255)＞土地利用∩降水(0.251)＞日照∩降水(0.250)，坡向与土壤类型、降水之间的交互作用解释力较大(图 5b₁)。2019年解释力高的6组为：土地利用∩坡向(0.476)＞土地利用∩降水(0.442)＞土壤类型∩土地利用(0.407)＞土地利用∩高程(0.392 0)＞土地利用∩积温(0.391 9)＞土地利用∩气温(0.391)(图 5b₂)，该时期土地利用与其他因子交互作用的解释力占据主导地位。该地区植被NDVI变化的主导驱动因子由自然因子转变为人类活动因子，这与区域经济发展和城镇建设相关^{[17, 19]}。⑥1981年西南地区交互因子解释力最大的为高程∩日照(0.810)，最小的是坡向∩坡度(0.087)(图 5c₁)，2019年交互因子解释力最大的为高程∩日照(0.807)，最小的是坡向∩坡度(0.084)(图 5c₂)。该地区位于西南山区，海拔较高，多高山、高原，其地形因素决定了区域生态环境特征及变化，因此高程是影响地表植被生态系统演变的主导因素之一。此外，西南地区多阴雨天气，光照不充足，极大限制了植被的生长，因此，日照成为影响植被变化的另一重要因子^{[23, 27]}。同时，西南山区生态环境极其脆弱，生态环境的保护一直是区域经济发展中亟待解决的前瞻性问题，特别是党的十八大以来，生态建设以保护为主，人为对生态的扰动逐渐变小，自然因素在植被生态系统演变中占据主导地位。

综上分析，近40 a不同地区的主导交互因子类型发生了一定变化并且存在显著差异：东北地区，日照∩其他因子(1981年)→土地利用变化(2019年); 华北地区，日照∩气温(1981年)→土地利用∩降水(2019年)；西北地区，降水∩其他因子(1981)→土地利用∩其他因子(2019年)；华东地区，坡向∩降水(1981年)→土地利用∩坡向(2019年)；中南地区，坡向∩土壤类型(1981年)→土地利用∩其他因子(2019年)；西南地区，地形因子∩其他因子→土地利用∩其他因子(2019年)。随着人类活动干扰强度的增加，土地利用及其相关交互因子逐渐成为各子区植被生态系统演变的主导因子。

4 结论

(1) 近40 a中国植被NDVI重心向南迁移，而不同分区的植被NDVI重心迁移轨迹存在一定差异性：东北地区植被NDVI重心表现为先向南移动继而向北移动的趋势；华北地区植被NDVI重心表现为向南迁移；西北地区植被NDVI重心表现为先向西北移动继而东南移动的趋势；华东地区植被NDVI重心呈向南移动趋势；中南地区植被NDVI重心迁移轨迹呈现“之”字形，总体表现为向西南方向移动；西南地区植被NDVI重心表现为向西南方向移动趋势。

(2) 植被NDVI与气温呈负相关区域主要分布于西北地区和华北地区，呈正相关区域则主要分布于黄土高原地区、四川盆地的中部、云贵高原以及长江中下游地区。植被NDVI与降水呈正相关区域主要分布于华北地区的中部，西北地区的东北部，呈负相关区域则主要分布于云贵高原和长江中下游等地区。植被NDVI与日照呈显著正相关区域主要分布于新疆中部、甘肃和青海省的西北部以及云贵高原西部，而显著负相关区域则主要分布于中部，包括贵州、山西和河南等地区。

(3) 近40 a不同分区、不同历史时期的植被NDVI演变主导驱动因子存在显著差异：东北地区，日照∩其他因子(1981年)→土地利用变化(2019年); 华北地区，日照∩气温(1981年)→土地利用∩降水(2019年)；西北地区，降水∩其他因子(1981年)→土地利用∩其他因子(2019年)；华东地区，坡向∩降水(1981年)→土地利用∩坡向(2019年)；中南地区，坡向∩土壤类型(1981年)→土地利用∩其他因子(2019年)；西南地区，地形因子∩其他因子→土地利用∩其他因子(2019年)。