河北省张承地区2001—2020年植被动态变化及驱动因素分析

王 倩，金晓媚，张绪财，殷秀兰，金爱芳，罗绪富

(1.中国地质大学(北京)水资源与环境学院，北京 100083；2.中国地质环境监测院，北京 100081)

0 引 言

植被是感知生态环境的敏感因子[1]，是生态系统赖以生存的基础，在气候、水文和生化循环中起非常重要的作用[2]，对风沙防治、水质改善和水资源调节发挥重要作用[3]。植被具有年际变化和季节变化的特点，其动态变化在一定程度上反映土地覆盖变化和气候变化趋势。而气候因子、地下水位埋深、土壤湿度和人类活动等也影响植被生长状况。因此，全面系统分析植被变化及其潜在驱动因素，成为近些年来人们研究陆地表层系统的主要内容之一[4]。

归一化植被指数(NDVI)被人们用来表示植被生长状况的好坏，也是目前应用最广泛的植被指数。国外一些学者基于NDVI对植被变化特征与气候因子等因素的关系进行了许多探讨[5-9]，国内也有相关的研究[10-12]。目前也有学者对研究区进行植被变化研究，郭伟伟等利用张家口2006年和2010年两景TM数据提取植被指数分析该市植被变化[13-14]。张晓凤利用2000—2013年的MODIS-NDVI数据分析坝上地区植被覆盖变化及气候因子对其影响[15]。邵雅琪等研究了张承地区2000—2014年植被指数的变化特征及其与气候因子的相关关系[1]。影响植被生长变化的往往是多方面的综合因素。但以上研究都只是分析植被与气候因子的关系，针对张承地区植被变化更全面的影响因素的研究很少，特别是植被与地下水位埋深及土壤湿度的关系研究。

2022年北京冬奥会成功举办后，为合理利用张家口赛区体育场馆，目前政府正在推进“京张体育文化旅游带”建设，并且张承(张家口—承德)地区又是北京市及京津冀地区重要的生态屏障，其生态环境的变化关系到首都地区的生态安全，因此对其生态环境和水源涵养都有更高水准的要求。本文基于MODIS数据，对张承地区2001—2020年的植被时空变化特征进行分析，并探讨NDVI与气候因子、用地类型、地下水位埋深和土壤湿度因素之间的响应关系，为区内的地下水涵养的进一步研究提供参考依据，也为推动研究区后奥运经济发展提供生态建设参考。

1 研究区概况

张家口市和承德市(以下简称“张承地区”)位于河北省西北部，处于内蒙古高原与华北平原过渡带[16]。全区地理位置为东经113°50′—119°30′、北纬 39°30′—42°40′，总面积约76276.8 km2(图1)。张承地区主要分坝上、坝下两个不同的自然地貌单元。

图1 张承地区地理位置及地形地貌Fig.1 Map showing the geographic location and topography of Zhang-Cheng area

张承地区属于温带大陆性季风气候，一年四季特征明显。由于受地形和气流影响，张承地区气候类型存在差异，坝上高原地区气候受蒙古高压控制，坝下气候主要受地形和纬度影响，雨季多集中在夏季，整体上降水量西北少，东南多，张家口市降水量少于承德市，坝上降水量少于坝下。张承地区多年平均气温在7 ℃左右，多年平均降水量在430 mm左右。

2 研究数据及方法

2.1 数据来源

2.1.1 MODIS数据

本文选取2001—2020年间MOD13Q1遥感数据，其时间尺度为16 d，空间分辨率为250 m，共460景数据。数据来源于美国NASA网站(http：//reverb.echo.nasa.gov/reverb)。

2.1.2 土地利用数据

为分析土地利用类型变化，本文选用2000年、2005年、2010年、2015年和2020年共5期土地利用数据，空间分辨率为30 m，数据来源于中国科学院资源环境科学数据中心(http：//www.resdc.cn)。

2.1.3 气候数据

为分析植被和气候因子之间的相关性，本文选择2001—2020年间ERA5数据中气温和降水的栅格数据，分辨率为0.1 rad，数据来源于GEE平台。

2.1.4 其他数据

为研究地下水位与植被的关系，收集了2020年张承地区59个潜水地下水位监测孔数据，数据来源于中国地质环境监测院。

为研究土壤湿度与植被的关系，2021年7月份在研究区内进行实地测量，得到99个测量点位的土壤湿度数据。

2.2 研究方法

2.2.1 植被覆盖率计算

本文植被覆盖率是由植被指数NDVI计算得到，计算公式如下：

(1)

式中：fc为植被覆盖率，NDVI是所求像元的植被指数，NDVImax是纯植被覆盖的植被指数，NDVImin是裸土表面的植被指数。结合实际情况本文设置NDVImax为0.8，NDVImin为0.1。

2.2.2 Sen+Mann-Kendall 趋势分析

Sen+Mann-Kendall趋势分析法是将Sen趋势分析法和Mann-Kendall方法相结合，这样可以减少噪声干扰和实现序列趋势的显著性检验[17-18]。首先根据趋势度β来判断时间序列的趋势变化，若β>0时，NDVI呈上升趋势，β<0时，NDVI呈下降趋势。然后进行Mann-Kendall趋势分析，因为本文时间序列数为20，所以要计算相应统计量Z，当|Z|>Z1-α/2时，表示序列趋势通过0.05显著性水平检验，序列有一定的变化趋势，而且趋势明显，反之趋势不明显[19-21]。

2.2.3 偏相关性分析

简单相关性分析表示两个独立变量之间的相关程度，用于分析NDVI与气温、NDVI与降水的相关关系[22]，偏相关性分析是假设在其中一个变量不变时，另外两个变量之间的相关关系，用偏相关性系数表示，其计算公式如下：

(2)

Rxyz的取值区间为[-1，1]，绝对值越接近于1，相关性越强。rxy、rxz、ryz分别代表NDVI与降水、NDVI与气温、降水与气温的相关系数。然后进行检验，偏相关系数的检验一般采用t检验方法。

2.2.4 插值分析

空间插值是将实际测量的数据点转换成连续数据面。反距离权重(IDW)插值法是一种全局插值法，其结果较为合理，不会出现无意义的结果。本文用到59个潜水埋深监测孔数据和99个土壤湿度实测点数据，将点数据导入到软件ArcGIS中，通过反距离权重(IDW)插值，得到与NDVI分辨率(250 m)一致的地下水位埋深空间分布和土壤湿度空间分布。

3 研究结果

3.1 植被覆盖率的时空分布

研究区年均NDVI值在0.1～0.8之间，其中6—9月NDVI值较高，为植被的生长季，近20年NDVI年均值为波动上升趋势，增长速率约为0.4 ha-1。根据年均NDVI计算得到研究区2020年的植被覆盖率分布图(图2)。本文参考其他学者研究，结合张承地区实际情况，将该地区植被覆盖率的分级标准确定为：低植被覆盖(fc=0～20%)、较低植被覆盖(fc=20%～40%)、 中等植被覆盖(fc=40%～60%)、较高植被覆盖(fc=60%～80%)、高植被覆盖(fc=80%～100%)。

图2 张承地区2020年植被覆盖率空间分布图Fig.2 Spatial distribution map of vegetation coverage in Zhang-Cheng area in 2020

从图2可以看出，研究区植被覆盖率在空间分布上存在西北低、东南高的特点，由于张承地区地形存在差异，海拔由西北向东南逐渐降低，植被呈现出较明显的空间异质性。总体上看，坝下地区的植被覆盖率明显高于坝上地区，且承德地区明显高于张家口地区，在张家口地区东部的植被覆盖好于西部。低植被覆盖和较低植被覆盖主要分布在张家口的康保县、张北县、沽源县、尚义县、阳原县、蔚县和怀安县。表1为2020年各级覆盖率的面积统计结果，植被覆盖率小于0的为水体，面积占比为0.14%；较低植被覆盖、中等植被覆盖和较高植被覆盖的面积占比分别是32.54%、35.40%和28.90%；低和高植被覆盖率面积分别是1380.54 km2和923.47 km2，共占研究区总面积的3.02%。

表1 2020年植被覆盖率分级面积占比统计结果Table 1 Statistics of classified area of vegetation coverage in 2020

在此基础上，本文进一步统计了2001—2020年不同植被覆盖率面积的变化趋势，如图3所示：低植被覆盖率的面积在2001—2004年逐年下降，在2009年突然增大，之后呈波动下降趋势；较低植被覆盖面积为下降趋势；中等植被覆盖的面积在2001—2003年逐年增加，之后波动下降；较高植被覆盖面积则呈现明显增加的趋势，高植被覆盖面积在2013年后呈现大幅增加。

3.2 植被变化趋势分析

利用Sen+Mann-Kendall法对研究区2001—2020年NDVI平均值进行计算，并对计算结果划分等级，分析研究区植被指数的空间变化趋势。由表2和图4可知，研究区内近20年植被整体上明显改善，占据总研究区85.79%的面积，特别是坝下地区；轻微改善和基本不变的地区主要分布在坝上，分别占总面积的7.37%和5.49%；在城镇地区由于扩建，存在植被退化的情况，植被退化面积占总面积的1.35%。

表2 2001—2020年Sen+Mann-Kendall法趋势分析结果Table 2 Results of Sen+Mann-Kendall trend analysis from 2001 to 2020

图4 张承地区20年间植被覆盖变化趋势分布Fig.4 Spatial distribution of vegetation cover variation trend in Zhang-Cheng area in 20 years

为分析研究区不同时间段内植被指数的变化特征，本文以5年为一周期将研究时间段分为2001—2005年、2006—2010年、2011—2015年和2016—2020年4个时间段，由图5和表3可知，2001—2005年间坝上及城镇地区存在植被退化状况，面积占比为4.72%，整个研究区植被改善情况面积占比达92.14%；2006—2010年间，这五年间研究区内存在大范围的植被明显退化情况，面积占比达到29.4%，其中康保、张北县、怀安县、蔚县、丰宁、围场和平泉市等地植被退化严重；2011—2015年，康保县、张北县北部、沽源县北部和丰宁西北部的植被存在严重退化情况，明显退化面积占比为7.53%，相比2006—2010年植被退化面积有所减少；2016—2020年间植被整体改善，明显改善情况主要分布在张家口中部和承德中北部，面积占比为66.07%，只有在研究区南部零星分布严重退化情况，面积占比为3.50%。

表3 张承地区每五年植被覆盖空间变化趋势面积占比(%)Table 3 Proportion of vegetation cover spatial variation trend in Zhang-Cheng area for every five years (%)

图5 张承地区每五年植被空间变化趋势图Fig.5 Spatial variation trend of vegetation cover in Zhang-Cheng area for every five years(a)2001—2005年；(b)2006—2010年；(c)2011—2015年；(d)2016—2020年

4 植被影响因素分析

4.1 气候因子

本文采用2001—2020年ERA5数据中气温和降水量数据，分别与对应的NDVI建立偏相关性分析。近20年的张承地区植被与降水偏相关系数分布在-0.72～0.80之间，植被与气温相关系数分布在-0.73～0.88之间。由表4和图6可知，张家口地区的植被分布受降水的影响大于承德地区，其中显著正相关的区域主要分布在尚义、怀安、逐鹿、蔚县以及沽源西部地区的，面积占比为5.92%；不显著正相关的区域面积占总面积的78.41%；相关性较差的区域主要分布在承德的隆化县、滦平县、平泉市及承德县北部地区的，面积占比为15.67%。研究区大部分地区的植被分布与气温呈不显著正相关关系，面积占比达72.16%；呈现显著正相关关系的区域面积占比为25.06%，主要分布在张家口地区的尚义县、张北县的东部、崇礼区的东北部、怀安县南部、阳原和蔚县的大部分地区以及承德的围场西北部地区、兴隆县的南部；相关性较差的区域面积占比仅有2.77%。综上，张承地区的植被分布受气候因子影响较大，植被指数随着气温和降水的增加而增大。

表4 植被与降水、气温相关性显著分级统计Table 4 Classification statistics of significant correlation between NDVI and precipitation and temperature

图6 张承地区植被指数与降水(a)和气温(b)的偏相关性分析检验结果Fig.6 Spatial distribution of partial correlation analysis between NDVI and precipitation (a) and temperature (b) in Zhang-Cheng area

4.2 用地类型

用地类型对植被的影响间接反映了人类活动对植被的影响，本文统计研究区2000年、2005年、2010年、2015年和2020年的不同用地类型的面积，由表5可知，研究区用地类型主要是耕地、林地和草地。2000—2005年用地类型变化较小，草地和水域面积略有减少，其它用地类型面积皆略有增加，此阶段各用地类型处于自然转化中，受人类活动影响较小，此时期的植被也在改善(图5(a))。2006—2010年，用地类型变化相对较大，主要受国家政策影响，由于退耕还林等政策和人工植树造林[23-24]，耕地面积减少了638.7 km2，林地面积增加了378.3 km2；受城市扩张的影响，建筑用地面积急剧增长了1019.7 km2；草地、水域和未利用土地面积均有所减少，这五年人类活动对土地利用类型变化影响较大，植被也在此阶段有明显退化(图5(b))。2011—2015年，这一阶段人类活动的影响相对减弱，除建筑用地面积增加了212.5 km2、草地减小了118.8 km2外，其它用地类型面积均变化不大，相比较前面5年此阶段植被退化减少，植被改善面积增加(图5(c))。2016—2020年，建筑用地面积增加了248.7 km2，草地和水域面积有少许增加，耕地面积减少了363.9 km2，林地和未利用土地面积变化较小，这一阶段随城市扩张人类活动对土地利用变化的影响又有所加大。结合图5可以发现，研究区植被变化与土地利用变化具有响应联系，所以人类活动对研究区植被变化作用明显。

4.3 地下水位埋深

植被与地下水的相互关系较复杂，地下水位埋深只有在一定范围内对植被生长有影响，且不同植被对地下水的敏感度不同[25-27]。张承地区属于典型大陆性气候，自然植被与地下水有一定关系，因此为研究张承地区植被与地下水之间的相互作用关系，将研究区2020年59个潜水地下水位监测孔的数据进行反距离权重(IDW)插值，得到与NDVI分辨率一致地下水位埋深栅格图像，然后统计2020年对应不同地下水位埋深的NDVI像元点的个数，得到NDVI像元值与地下水位埋深的统计图。由图7(a)可知，NDVI像素点主要集中在地下水位埋深5～12 m范围内，特别在6～9 m内NDVI像素点较多，随着地下水位埋深增加，NDVI像素点个数逐渐减小，在大于12 m后，植被生长受其影响较小。因此将NDVI分布图与地下水位埋深小于12 m范围的栅格图进行叠加，得到研究区NDVI平均值与地下水位埋深的关系图(图7(b))。由图7(b)可以看出，NDVI值随着地下水位埋深的增加呈现出先增加后减小的趋势，在地下水位埋深到6.2 m时，NDVI值出现峰值0.53，在地下水位埋深大于8 m后，NDVI值随着地下水位埋深增加而减小，因此研究区影响植被生长的地下水位埋深范围为5～8 m，在此范围内植被生长受地下水影响较大，当地下水位埋深超过8 m后，植被与地下水位埋深的关系变弱。

4.4 土壤湿度

土壤湿度作为地表特征参数之一，是植被蒸腾作用和光合作用的限制因素，是植被生长、植被恢复的重要影响因素。特别是人工植被直接受土壤水分的影响作用，同时地表植被本身的保水耗水及覆盖状况也会对土壤湿度进行反馈作用[28]。根据2021年7月野外实测平均土壤含水量数据，进行反距离权重(IDW)插值，并将分辨率重采样为和NDVI一样，得到0～10 cm (图8(a))和0～20 cm (图8(b))的表层土壤湿度空间分布图。由图可知土壤湿度范围为3%～39%，空间分布也存在“西北低，东南高”的特点，坝上土壤湿度相比坝下较低，张家口中北部地区除张北县的安固里淖附近土壤湿度较高外其它地方均较低，承德市的土壤湿度较高于张家口市，整体上的分布与植被覆盖率空间分布较契合。然后将NDVI分布图分别与0～10 cm和0～20 cm的表层土壤湿度栅格图进行叠加，在ArcGIS里以1%为间隔重分类土壤湿度，提取每个区间内的NDVI算数平均值，绘制NDVI算数平均值与土壤湿度的变化曲线。由图9可知，在10 cm和20 cm土层深度的情况下，NDVI值整体趋势上都是随着土壤湿度的增加而增大，两者拟合方程的判定系数(R2)分别达到0.7277和0.7606，说明NDVI与各层土壤湿度的数据拟合较好，表层土壤湿度对植被生长具有一定影响作用。

图8 实测张承地区土壤湿度空间分布Fig.8 Spatial distribution map of measured soil moisture in Zhang-Cheng area(a)0～10 cm表层土；(b)0～20 cm表层土

5 讨 论

20世纪60年代，风沙紧逼北京城，国家为治理风沙，实施一系列的生态保护政策，如建设国有林场等，其中成果最显著的是在承德塞罕坝上建成百万亩林海，因此在国家政策的大背景下，张承地区近20年年均NDVI为波动上升，植被变化趋势也是整体在改善。由于张承地区坝上坝下地形差异较大，在空间分布上植被覆盖度承德地区高于张家口地区、坝下高于坝上，呈现为西北低、东南高的特点，在2001—2004年低植被覆盖面积相对明显降低，在2013年后高植被覆盖面积急剧增加，整体植被覆盖度在增加，这与邵雅琪等人在张承地区研究结果一致[1]。虽然植被总体在改善，但在2006—2010年间植被退化严重，这与城乡建筑扩张以及人们为追求经济效益过度放牧与开垦耕地有关[24]，之后植被覆盖又逐步改善，也与当地经济生态化的政策[29]及成为2022年北京冬奥会比赛场地之一有关。

影响张承植被变化的因素中，气候因子与用地类型、地下水位埋深、土壤湿度共同作用，其中气候因子降水与气温对张承地区的植被均起到促进作用。张家口地区相较于承德干旱，植被对降水的正相关性高于承德地区，这与前人研究结果一致[15]。人类活动作用使土地利用类型发生变化，建筑用地面积近20年逐年增加，尤其在2006—2010年急剧增加，而此时期植被退化严重。耕地面积在2005年后持续减少；近20年林地面积总的在增加，草地面积总体在减少，这与国家退耕还林、以及防护林建设政策有关。用地类型的变化与前文植被变化趋势相一致，说明人类活动也是驱动植被变化的重要因素。崔浩楠等研究表明长江经济带由于经济建设和生态修复工程使人类活动成为植被覆盖变化的重要驱动力[4]。不同地区的地下水位埋深对植被影响不同，适宜植被生长的埋深也不同。张承地区由NDVI与地下水位埋深曲线关系可以得出影响植被生长的地下水位埋深范围为5～8 m。土壤湿度的空间分布与地形有响应关系，具有明显的坝上土壤湿度低于坝下地区，植被覆盖与土壤湿度较强的响应关系，NDVI随土壤湿度增大而增大，两者之间有相互促进作用。与张京等人研究格尔木河流域植被对土壤湿度有积极的响应作用结果基本一致[30]。

上述结果表明，张承地区植被变化由降水、气温、用地类型、地下水位埋深和土壤湿度共同作用，降水和气温对植被影响较为明显和直观，从年度降水和气温对植被影响上分析，气温影响程度大于降水；在长时间序列中人类活动转变用地类型对植被影响也较为明显，根据野外观测可知研究区多是人工种植林，植被覆盖明显受人类活动影响；研究区特别是张家口地区地下水位埋深较深，因此地下水位埋深对植被生长影响程度较小。土壤湿度与植被有较明显的联系，承德的土壤湿度明显大于张家口，而承德的植被覆盖也明显高于张家口。综合以上各因素对植被的影响程度，按其大小做如下排序：气温>降水>用地类型>土壤湿度>地下水位埋深。从结果上看，人类活动对植被动态变化的影响较明显，张承的生态经济发展需要侧重对土地利用的合理规划。本研究所用的地下水位埋深数据和土壤湿度均为实测数据，具有真实可靠性，也使本文更具有真实性。

6 结 论

本文基于MODIS遥感数据，对张家口—承德地区的植被指数和植被覆盖进行了计算，并对其时空分布特征及影响因素进行分析和探讨，得到了以下几点主要结论：

(1)研究区2001—2020年NDVI呈增加趋势，增长速率约为0.4 ha-1。植被覆盖率在空间分布上呈现“西北低，东南高”的特点，坝下植被覆盖高于坝上植被覆盖；受国家生态环境保护政策和当地发展规划影响，近20年低植被覆盖面积呈明显下降趋势，较高植被覆盖率的面积显著增加。

(2)2001—2020年研究区内绝大部分区域植被在改善，特别是在坝下地区，而在城镇区域存在植被退化。2001—2005年植被改善最为显著；2006—2010年植被退化明显，占29.4%的总面积；2011—2015年植被退化面积减少，改善面积增加；2016—2020年植被改善面积占比有85.14%。

(3)在年际尺度上植被指数与降水量和气温之间相关性较好，正相关面积占比分别为84.33%和97.22%，研究区植被随降水和气温的增加而增加。用地类型变化与5年为一周期的植被变化有较强的响应联系，说明人类活动使用地类型发生变化，进而构成研究区植被变化的驱动因素。

(4)研究区内影响植被生长的地下水位埋深为5～8 m，地下水位埋深超过8 m时，植被与地下水位埋深关系较弱，在地下水位埋深6.2 m的地方植被长势最好。土壤湿度在空间上呈现坝上低、坝下高的特点，NDVI随着土壤湿度的增大而增大，且拟合方程判定系数在0.7以上。表明研究区地下水位埋深与土壤湿度皆为植被变化影响因素。