植被作为地理环境的重要组成部分，在陆地表层系统中所占的比例很高，是生物圈及其生态系统的核心和功能部分，其变化不仅对全球及区域能量循环、物质的生物化学循环具有重要的影响，同时也对区域及全球气候变化有重要的贡献^[1-2]。相较于其他数据，遥感数据具有时间和空间上的连续性，故常被用来监测全球和区域植被时空分布及变化趋势^[3-6]，为研究植被的空间分布格局和预测未来发展趋势奠定了坚实的基础。

归一化植被指数(normalized difference vegetation index，NDVI)对植被生物物理特征十分敏感，作为植被生长状况的重要度量参数，被有效使用在监测植被动态变化当中。国内外学者利用NDVI数据在区域和全球尺度研究植被覆盖时空演变特征和植被覆盖对气候变化的响应机制，并取得了丰硕的成果。de Jong等^[7]利用GIMMS NDVI3 g数据反演了全球植被覆盖变化趋势，并在像元尺度上探测了植被覆盖变化的转折点。结果表明，全球超过26%地区的植被呈单调显著变化。方精云等^[8]发现生长季延长和植被生长速度加快是中国植被覆盖状况改善的重要原因，且温度上升和夏季降水增加是导致植被活动增强的主要气候因子。孟丹等^[9]发现京津冀地区植被覆盖呈改善趋势，且整体上植被生长与降水相关关系为正，与气温相关关系为负。晏利斌等^[10]利用GIMMS NDVI3 g分析了京津冀地区生长季植被覆盖时空变化特征及其与气候因子的相关关系。结果表明，京津冀地区植被覆盖总体呈上升趋势，且降水对植被生长的作用强于气温。

在全球气候变暖和人类活动加强的背景下，研究植被覆盖动态变化及其与气候变化之间的相关关系和响应机制具有重要的现实意义。已有研究大多侧重于探究植被覆盖时空变化特征^[11-14]，而对预测其未来变化趋势鲜有研究，且在分析植被NDVI与降水和气温相关关系时，并未顾及地理位置信息给降水和气温插值精度带来的影响，降低了研究结果的精度。因此，本文基于2001—2019年MODIS NDVI时间序列，利用Theil-Sen Median趋势分析法、Mann-Kendall显著性检验法、R/S分析法和Person相关分析等数学方法，结合ANUSPLIN气象插值模型，反演2001—2019年京津冀地区植被覆盖时空演变特征，并预测其未来变化趋势，进一步分析植被NDVI对降水和气温变化的响应机制和时滞效应，研究结果可为区域乃至全球植被覆盖时空演变特征及其与降水和气温响应关系研究提供理论支撑，对京津冀地区林业生态工程的实施效果评估和生态环境建设结果评估具有重要的现实意义。

1 研究区概况、数据来源及研究方法 1.1 研究区概况

京津冀地区(113°04′—119°53′E，36°01′—42°37′N)位于华北北端，地表形态复杂、地形起伏较大，主要表现为东南低、西北高。京津冀地区由北京和天津两个直辖市，以及河北省组成，总面积约为2.17×10⁵ km²。受半湿润半干旱大陆性气候的影响，四季分明，冬季低温干燥，夏季高温湿润。年平均气温10.4~11.9 ℃，年累计平均降水量375.5~684.7 mm。西北部土地利用类型主要以林地和草地为主，东南部主要以农用地为主。京津冀地区位于环渤海经济圈核心区，经济发展水平、城市化速度和人口密度均高于其他地区。人口总数占全国总人口的8.07%(中国统计年鉴2019)，而面积仅占全国总面积的2.25%。受气候变化和人类活动的影响，该地区生态环境较为脆弱，植被生长对气候变化的响应敏感。

1.2 数据来源 1.2.1 MODIS NDVI遥感数据

MODIS NDVI数据来源于美国国家航空航天局发布的MOD₁₃Q₁C₆数据集，时间跨度为2001年1月至2019年12月，时间分辨率为16 d，空间分辨率为250 m，每年23个时相。首先，利用MRT(MODIS Re-Projection Tools)，对原始格式为HDF的MOD₁₃Q₁C₆数据集进行NDVI波段提取、批量镶嵌、重采样和投影转换，输出数据格式为Geo-tiff的影像文件，投影为Alberts Equal Area。然后，基于最大化合成法(maximum value composition, MVC)得到月最大NDVI时间序列，该处理可以减少大气中云、颗粒以及太阳高度角的影响。最后，利用均值法计算得到年平均NDVI时间序列，并裁剪出2001—2019年覆盖研究区的NDVI时间序列。

1.2.2 DEM数据和气象数据

本文所使用的DEM数据来源于美国国家航空航天局发布的SRTM DEM数据，空间分辨率30 m。为了与NDVI时间序列保持相同空间分辨率，将DEM重采样为250 m空间分辨率。气象数据来源于国家气象科学技术中心提供的《中国地面气候资料月值数据集》(http://data.cma.cn/data/cdcdetail/dataCode/SURF_CLI_CHN_MUL_MON.html)，主要包括2000—2019年降水和气温数据，时间分辨率为月，气温的精度为0.1 ℃，降水量精度为0.1 mm，数据经过严格的精度控制，质量良好。为了提高研究结果的精度，本文选取覆盖研究区及其周边地区的65个气象站点数据，以DEM作为协变量，基于ANUSPLIN插值模型生成降水和气温时间序列。

1.3 研究方法 1.3.1 Theil-Sen Median趋势分析和Mann-Kendall显著性检验

趋势分析是一种非参数统计的斜率估计方法。由于该方法不需要数据服从特定的分布特征，且受异常值的影响较小，故常用于估算水文、气象、植被指数等长时间序列数据的变化趋势^[15-18]。Mann-Kendall显著性检验是一种科学的数学统计方法，可用于检验时间序列数据变化趋势的显著程度^[19]。本文将两种方法结合起来，用于估算植被NDVI变化趋势，并根据变化斜率判断其变化趋势的显著程度。

Theil-Sen Median趋势分析用来计算植被NDVI的变化趋势，用slope表示NDVI变化趋势，当slope>0时，植被NDVI呈上升趋势；当slope=0时，植被NDVI基本保持不变；当slope < 0时，植被NDVI呈下降趋势。Mann-Kendall显著性检验法被用来检测NDVI时间序列变化趋势的显著性，在给定显著性水平α下，当|Z|>Z_1-α/2，表示时间序列在α水平上变化显著，反之，则变化不显著。本文定义变化趋势在α=0.05下显著时，为显著变化；在α=0.01下显著时，为极显著变化。

1.3.2 R/S分析及Hurst指数

R/S分析法(rescaled rang analysis method)，又称之为重新标度极差分析法，最早是由Hurst在总结尼罗河的多年水文观测资料时提出的一种分析方法^[20]，后来经过Mandelbrot和Wallis进一步补充和完善，将其发展成一种研究时间序列的分型理论^[21]，目前在水文学、经济学、气候学、地质学等领域有着广泛的应用。Hurst指数的取值范围为(0＜H＜1)。不同的H，意味着改时间序列是完全随机的或者是有不同趋势的变化。若0＜H＜0.5，表示NDVI时间序列数据具有反持续性；若H=0.5，表示NDVI时间序列为相互独立的随机序列，若0.5＜H＜1，表明NDVI时间序列具有长期相关性，过程具有持续性。

1.3.3 Pearson相关分析

Pearson相关分析可用来分析两列数据之间是否存在相关关系及相关程度^[22]，相关系数越大，代表两个因子间的相关性越好，反之，则代表相关性越差。以往研究表明，植被NDVI对降水和气温变化具有时滞效应，且滞后期往往为0~3个月^{[6, 23-24]}，故本文基于MATLAB在像元尺度上计算植被NDVI与同期及前期1~3个月降水和气温的相关系数，并根据其最大相关系数，得到植被NDVI对降水和气候变化响应的滞后期。

2 结果与分析 2.1 植被覆盖时空演变特征 2.1.1 植被NDVI时间变化特征

如图 1所示，2001—2019年京津冀地区植被NDVI整体呈波动上升趋势，上升速率为0.002 2/a。京津冀地区植被NDVI最大值为0.432，最小值为0.381，分别出现在2018年和2001年。北京、天津和河北植被NDVI变化趋势表现出较大的差异，北京和河北植被NDVI呈上升趋势，且北京植被NDVI上升趋势明显高于河北，而天津植被NDVI呈波动下降趋势。以上结果表明，2001—2019年京津冀地区植被覆盖整体呈现改善态势，尤以北京改善程度最为显著，但天津植被覆盖呈现退化趋势。

2.1.2 植被NDVI空间变化特征

如图 2a和表 1所示，2001—2019年京津冀地区植被NDVI变化趋势呈现出明显的空间异质性。植被NDVI上升区域占总面积的76.19%，其中极显著上升(p＜0.01)约占总面积的49.31%，从京津冀地区东南部呈带状延伸至西北部，主要包括北京、保定和石家庄西北地区、张家口和承德南部地区、秦皇岛北部以及沧州、邢台和衡水东南部分地区。以上地区土地利用类型以林地和草地为主，地势较高，受人类负向扰动较小。自1998年以来，得益于退耕还林还草工程、三北防护林工程和环北京和天津地区防沙治沙工程等一系列生态林业工程的实施，以上地区在植被群落恢复和生态系统建设等方面取得了显著的成果，植被覆盖呈明显上升趋势，荒漠化程度得到了有效的控制，荒漠化面积大大缩减^[25-28]。

植被NDVI下降区域占总面积的23.81%，其中，极显著下降(p＜0.01)和显著下降(p＜0.05)区域约占总面积的10.52%，主要分布在京津冀地区东南部各城市中心及其周边地区，尤以唐山、天津、廊坊、石家庄和邯郸为著。以上地区主要位于地势平坦的华北平原，地势较低，受社会经济发展、人类活动增强、人口密度增加等的负向影响，各城市中心及其周边地区植被活动减弱、植被覆盖呈下降趋势。

2.1.3 植被NDVI未来变化趋势

R/S分析常用来探究时间序列数据变化趋势的持续性与反持续性。京津冀地区植被NDVI Hurst指数均值为0.507，整体呈弱持续性变化。植被NDVI呈持续性变化的区域约占京津冀地区总面积的50.18%，略大于呈反持续性变化的区域(49.82%)。

为揭示京津冀地区植被NDVI未来变化的持续性特征，本文将趋势分析和R/S分析结果进行重分类并进行叠置分析，得到各像元未来植被NDVI变化趋势，并将结果分成4个类别(详见图 2b和表 2)：持续下降、下降、上升、持续上升。统计分析表明，未来植被NDVI下降和持续下降区域占总面积的一半以上(53.11%)，其中下降的比例占总面积的39.83%，均匀分布于京津冀地区；呈持续下降的区域占总面积的13.28%，主要分布在环渤海湾地区的天津、廊坊和沧州及其周边地区。未来植被NDVI呈持续上升的区域分布较为集中，主要分布张家口和承德大部分地区，其他城市零星分布，约占京津冀地区总面积的36.90%。

由以上结论可以推断，由于多项林业生态工程的实施，构建起了一张乔木、灌木和草地相结合的防护网，不仅起到涵养水源、保持水土的作用^[29-30]，也为张家口和承德及其周边地区植被生长和繁衍提供了良好的环境，使得该地区未来植被NDVI继续以上升为主；而天津、廊坊和沧州及其周边地区，由于受城市人口密度增加和城区面积扩张的影响，加速了耕地、林地和草地向城市建设用地的转换^[31]，使得未来植被覆盖继续保持退化态势，影响区域生态环境保护和生态文明建设。

2.2 植被NDVI对降水和气温的响应 2.2.1 植被NDVI与降水相关关系

为探究京津冀地区植被生长对气候变化的响应机制，计算2001—2019年京津冀地区植被NDVI与同期及前1~3个月降水和气温之间的相关系数，根据最大相关系数揭示京津冀地区植被NDVI对降水和气温变化最大响应的空间分布特征。在以上研究基础上，利用T检验法对植被NDVI与降水和气温最大相关关系的显著性进行检验，结果如图 3所示。

如图 3a所示，京津冀地区植被NDVI与降水的最大相关系数为-0.832至0.945，相关系数均值为0.319。植被NDVI对降水变化的响应强度呈现明显的地域差异，主要表现为西北高，东南低。由表 3可知，京津冀地区植被NDVI主要与降水呈正相关，相关系数为正的区域约占京津冀地区总面积的87.35%。

如图 3b所示，保定西北部、北京西部和北部、张家口东南部、沧州东部、以及承德西南部地区，植被NDVI与降水呈极显著和显著正相关，分别占总面积的17.80%和16.44%。以上地区由于受生态林业工程的影响，植被覆盖程度呈逐年递增的趋势，且该地区属半湿润半干旱地区，水分是限制植被生长的主要气象因子^[10]。一方面，区域植被覆盖状况的改善能够起到防止水土流失、涵养水源的作用。另一方面，水分条件的改善能够促进区域植被的生长和繁衍进程。植被NDVI与降水呈极显著和显著负相关的区域约占总面积的0.63%，主要集中在廊坊和邯郸各城市中心及其周边地区。以上地区由于受高速经济发展和快速人口增长的影响，城区面积由中心向边缘扩展，导致植被覆盖呈下降趋势，而该地区研究时段内降水呈上升趋势，故植被NDVI与降水呈负相关。

2.2.2 植被NDVI与气温相关关系

如图 4a所示，京津冀地区植被NDVI与气温的最大相关系数为-0.804至0.906。植被NDVI对气温变化的最大响应强度由东南到西北依次呈现“弱—强—弱”的空间分布格局。由表 3可知，京津冀地区气温对植被生长的正向促进作用远小于降水，植被NDVI与气温的最大相关关系均值仅为0.132，由此可以推断降水是京津冀地区限制植被生长的主要气候因子，这与晏利斌等^[10]研究结果一致，即京津冀地区植被NDVI与降水和气温变化均呈正相关，但与降水的相关性更加显著。

如图 4b和表 3所示，植被NDVI与气温最大相关系数为正的区域约占研究区总面积的69.14%。其中，植被NDVI与气温呈极显著正相关和显著正相关的区域约占总面积的9.90%，呈带状零星分布于京津冀地区腹地。植被NDVI与气温呈负相关的面积约占京津冀地区总面积的30.86%，主要分布在京津冀地区两翼。其中，呈显著负相关和极显著负相关的面积约占总面积的3.23%，主要分布在京津冀地区东南沿线各城市中心及其周边地区。由于受城市扩张和人类活动的影响，以上地区植被覆盖状况较差且抵御外界扰动的能力较弱，且城市周边地区土地利用类型以耕地为主，受人类农业管理活动的影响较大，京津冀地区部分区域植被NDVI对气温变化呈不显著相关甚至负相关^{[19, 32]}，表明气温对植被生长具有较强的抑制作用。

2.2.3 植被NDVI对降水和气温变化最大响应的滞后效应

考虑到植被生长对降水和气温变化具有时滞效应，本文通过植被NDVI与降水和气温的最大相关系数，得到京津冀地区植被NDVI对降水和气候变化最大响应滞后期，揭示植被生长对降水和气温变化时滞效应的空间分布特征。

如图 5a和表 4所示，京津冀地区植被生长对降水变化存在明显滞后效应，且滞后期表现出明显的地域分异格局。植被NDVI对前2月和前3月降水变化最大响应的区域约占总面积的58.74%；其中，植被NDVI对前3月降水最大响应的面积最大，约为44.22%，主要分布在张家口西北部、承德东北部、唐山东部、沧州、邯郸、石家庄西部和保定中部；植被NDVI对前1月降水最大响应的区域约占京津冀地区总面积的27.33%，仅次于前3月，呈片状分布在张家口东南部、承德西南部和北京。

如图 5b和表 4所示，相较于降水，京津冀地区植被生长对气温变化的滞后期较短。植被NDVI对当月和前1月气温变化最大响应的区域约占总面积的68.41%；其中，植被NDVI对当月气温变化最大响应的区域约占总面积的32.81%，呈片状分布在张家口及其周边地区；植被NDVI对前1月气温变化最大响应的面积最大，约为35.60%，呈带状从京津冀地区东北部延伸至西南部。植被NDVI对前3月气温变化最大响应区域仅占总面积的15.23%，主要集中分布在保定、石家庄和秦皇岛西北部。

由以上可知，京津冀地区植被NDVI对降水和气温变化均存在明显的时滞效应，植被NDVI对降水变化的滞后期略长于气温，这一结论与前人的研究结果一致^[23-24]。降水过程中的部分降水经地表径流流出，不能通过渗透作用进入土壤；而另一部分可被土壤吸收，补偿土壤水分，改善土壤墒情。只有保留在土壤部分的降水，可以被认为是植被生长的“有效降水”。降水从地表经过渗透作用被土壤吸收，而后通过根部吸收转移到植物的各个组织，最后成为植物生长可用的水资源需要花费一定的时间^[33-34]。因此，植物生长往往对降水变化具有一定的滞后效应。相较于降水，温度的变化可以直接影响植物的呼吸作用强度和光合作用效率，对植被生长具有更加直接的作用，故植被生长对气温变化的滞后期更短。

3 结论

本文利用多源数据和多数学分析模型，在不同时空尺度上定量反演了近19 a京津冀地区植被覆盖时空演变特征，并预测其未来变化趋势，进一步探究了植被NDVI对降水和气温变化的时空响应特征及时滞效应。主要研究结论如下：

(1) 得益于林业生态工程的实施，京津冀地区植被覆盖整体呈现改善态势。近19 a来京津冀地区植被NDVI上升区域约占总面积的76.19%。其中，极显著上升约占总面积的49.31%，主要集中在京津冀地区西北部。

(2) 京津冀地区未来植被呈下降趋势的面积略大于上升趋势的面积，其中呈持续下降的区域约占总面积的13.28%，主要分布在天津、廊坊和沧州及其周边地区，以上地区受社会经济发展、城市扩张和人口密度增加的负向影响较大。

(3) 植被NDVI对降水和气温变化的响应强度呈现明显的地域差异，植被NDVI与降水的相关性主要表现为西北高，东南低，而植被NDVI与气温的相关性由东南到西北依次呈现“弱—强—弱”的空间分布格局。植被NDVI与降水的相关关系略强于气温。因此，降水被认为是京津冀地区控制植被生长的主要气候因子。

(4) 京津冀地区植被生长对降水的气温的变化存在明显的时滞效应。植被NDVI对前3月降水最大响应的面积最大，约占总面积的44.22%；植被NDVI对前1月气温变化最大响应的面积最大，约占总面积的35.60%。