河北坝上地区2000-2019年植被绿度动态及其土地利用/覆被变化归因分析

王彦芳, 裴宏伟, 赵 超

(1.河北地质大学 河北省农业干旱遥感监测国际联合研究中心/河北省高校生态环境地质应用技术研发中心, 河北 石家庄050031; 2.河北建筑工程学院 市政与环境工程系, 河北 张家口 075000; 3.河北省气象科学研究所, 河北 石家庄 050021)

植被是陆地生态系统的重要组成部分，是连接大气、水圈、生物圈和土壤圈的天然纽带，通过光合作用、蒸散和地表反照率调节全球碳循环、水循环和能量交换等[1-4]。在全球范围内，植被活动的变化主要受生物地球化学效应和社会经济因素的影响，前者主要包括区域气候变化、二氧化碳的“施肥效应”和氮沉降等，人类活动主要包括土地利用变化和社会经济因素[5]。研究区域植被的连续变化及探讨其影响因素，能够反映环境演变的趋势和原因，对区域生态评价、生态安全和经济可持续发展具有重要意义[6]。

自然界的绿度与地表植被覆盖率直接相关，在中分辨率的遥感真彩色合成影像上，植被覆盖度越高，图像颜色越绿[4]。随着不同时间、空间及光谱分辨率的遥感卫星资源的持续丰富，众多研究者借助遥感手段来实现对地球表面的植物种类的识别以及地表植被的覆盖度等的研究，如基于遥感反演能灵敏地反映植被生长特征(覆盖度、叶面积指数、生物量等)[7-8]、生物物理生化特性(叶含水量、叶绿素含量等)[9-10]和生态系统参数(高度、长势等)等[11]。其中，植被指数(vegetation index)选用对绿色植物强吸收的可见光红光波段和对绿色植物高反射的近红外波段构建，成为衡量地表植被绿度的重要指标[4]。在众多植被指数中，归一化植被指数(normalized difference vegetation index, NDVI)是植被生长状态和植被绿度的最佳指示因子之一，在全球范围内广泛应用[5,12-13]。

已有的研究在区域及大陆尺度植被绿度变化规律、植被对气候因子的响应规律，以及气候变化与人类活动对植被变化的贡献展开广泛的研究[13-17]。研究表明，20世纪80年代以来，全球植被绿度不断改善，尤其是北半球中高纬地区[13]，中国西北、华北、青藏高原和黄土高原等地区植被绿度也在不断增长，不同区域植被变化受气候变化和人类活动影响程度不同[2,15-17]。土地是人类活动的空间载体，改变土地利用方式是人类活动影响植被绿度的重要手段，将土地利用数据与遥感NDVI相结合可以评估不同土地利用方式下区域植被生态状况。例如，Chen等[14]研究发现中国和印度的陆地植被面积对全球绿化增量的贡献约1/3，其中，中国占1/4，中国的绿化主要来自森林(42%)和农田(32%)，而印度的绿化主要来自农田(82%)，森林的贡献较小(4.4%)。景振华等[18]分析了山西六大煤田NDVI对土地覆盖类型的响应特征；道日娜[19]研究认为1995—2010年，三北地区NDVI减小的主要原因是大量的草地转为耕地，NDVI增加与大量未利用地转为耕地有关。但目前的研究，缺少结合土地利用方式变化的植被绿度变化的来源定量分析，明确区域植被绿度变化所发生的主要土地类型。因此，本研究以京津冀生态环境支撑区中重要的区域之一——河北坝上地区为例，采用2000—2019年的MODIS MOD13Q1 NDVI数据，结合2000，2010和2018年土地利用数据，分析坝上近20 a植被绿度变化及土地利用/覆被变化归因，旨在为地区生态建设和京津冀生态环境支撑区的建设提供一定参考，为土地合理利用、生态环境保护和地区可持续发展提供科学依据和决策支持。

1 研究区概况

坝上地区指京津冀地区北部，系内蒙古高原的南缘，其中本研究涉及的河北坝上地区包括张家口市的4县两区(张北县、康保县、尚义县、沽源县、察北管理区、塞北管理区)和承德市的丰宁满族自治县和围场满族蒙古族自治县一部分。该区域总面积近1.90×104km2，占河北省总面积约10%，海拔多在1 200～1 800 m，地形为丘陵和平原，呈东南高、西北低格局(图1)；河北坝上人口约112万人，占河北省总人口的1.7%[20]。气候条件方面，坝上系大陆性季风气候，自然条件恶劣，表现在寒冷、多风、干旱等，坝上年均气温1～2 ℃，无霜期90～120 d，年降水量400 mm左右[20]。坝上以闪电河为界，东部属于坝上外流区，降雨较多，地表覆盖以森林草原和草甸草原为主；西部属于封闭的内流区，气候干旱，以干草原为主，历史上曾是一个原生态的自然生态系统，曾经河网密布，水淖丰富[21]。由于受气候变化和人类开滩种地和过度放牧等活动的影响，草场逐步退化，湖淖不断萎缩，地下水位下降明显。

图1 研究区范围及地形

2 数据与方法

2.1 数据来源

本研究所用的主要数据包括遥感NDVI数据和土地利用数据。NDVI数据为MODIS (moderate-resolution Imaging Spectroradiometer， MODIS)MOD13Q1数据，时间序列为2000—2019年，空间分辨率为250 m。为了提高数据精度，该数据采用最大值合成方法(MVC)减少了云和气溶胶等的影响，合成后数据的时间分辨率为16 d。土地利用数据来源于中国科学院资源环境科学与数据中心(http:∥www.resdc.cn/)，包括2000，2010和2018年的3期坝上地区土地利用数据，该数据来源于国家尺度1∶10万比例尺多时期土地利用/土地覆盖遥感监测数据库(CNLUCC)，以美国陆地卫星Landsat遥感影像数据作为主信息源，通过人工目视解译获得，数据空间分辨率为30 m[22]。土地利用类型根据土地资源及利用属性包括耕地、林地、草地、水域、居民地和未利用土地6个1级类型，分类结果的综合评价精度在94.3%以上[22]。

2.2 分析方法

2.2.1 植被绿度 NDVI可以很好地反映区域植被绿度。根据植被绿度季节性变化差异，本研究同时采生长季平均NDVI和年最大NDVI分别评价区域年平均和最大绿度。由于受气候条件的影响，该区域的植被生长主要集中在4—10月。本研究以该时期内平均NDVI表示年生长季平均植被绿度。首先将MOD13Q1 NDVI每年16 d间隔的23个波段合成一个文件，生长季平均NDVI即波段范围为7～20的平均值。2000—2019年逐像元生长季平均NDVI和年最大NDVI计算公式为：

(1)

NDVImax,k=max(NDVIk,i)

(2)

2.2.2 线性倾向估计 本研究利用线性倾向定量的估算各栅格植被绿度变化趋势。对于样本量为n的变量xi来说，线性倾向估计通过建立变量序列xi和时间ti之间的一元线性回归关系[4]，即

(3)

式中：回归系数b表示变量x的趋势倾向，其绝对值大小上升或者下降的速率，符号表示变化方向，b可以通过最小二乘法拟合得到：

(4)

(5)

(6)

因此，每种土地利用类型对于区域植被绿度变化的贡献量为Rl

(7)

年最大绿度变化及其定量分析与生长季平均值相同。

3 结果和分析

3.1 河北坝上地区植被绿度时空变化特征

3.1.1 坝上地区植被绿度的空间分布特征 不同土地利用状况下，地区的植被绿度存在明显差异。2000—2019年坝上生长季多年平均植被绿度和年最大NDVI的多年平均态如图2a,2b所示。由图2可知，坝上生长季平均NDVI为0.39，年最大值均值为0.63，空间上呈现坝东高于坝西。植被最差的区域主要分布在干涸的湖淖区，其次是康保县北部，尚义县南部造林初期区域植被。NDVI变异系数可以体现植被稳定性，2000—2019年NDVI变异系数结果显示，坝东植被较坝西稳定，坝西地区变异系数浮动在15%之内，植被绿度变异最大为坝西湖淖区，其次是尚义南部和北部地区(图2c)。结合土地利用状况，除水域区，植被变异较大的区域主要是人类活动影响较大的区域，如耕地植被和建设用地的植被，其次是草地和未利用地，林地最为稳定。统计2000，2010和2018年不同土地类型的平均NDVI植被生长曲线(图2d)，结果表明，林地植被绿度最高，其次草地、耕地、未利用地，除水域区域外，建设用地区域植被绿度最低。从5月中旬起，人工植被即耕地和建设用地NDVI增长速率明显高于林地、草地等自然植被。

图2 坝上地区植被绿度空间分布特征

3.1.2 近20 a坝上地区植被绿度变化趋势 2000—2019年间，河北坝上地区植被绿度(生长季均值或年最大值)呈现整体增加的趋势，其中，年生长季均值增速为0.044/10 a，年最大值平均增速更大，达到0.063/10 a，植被绿度增加的区域占总面积的90%以上(图3)。增幅较大的区域主要分布在尚义县中南部、康保县和张北县西部、沽源县大部地区以及丰宁和围场的中北部地区，这些区域增加趋势的显著性水平也较高。康保县北部零星区域为下降趋势，但并不显著(图4)。

图3 河北坝上地区2000-2019年生长季平均NDVI和年最大NDVI的变化趋势

图4 河北坝上地区2000-2019年生长季平均NDVI和年最大NDVI变化趋势的显著性分析

生长季平均绿度呈显著(p<0.05)增加趋势的比例为83%，不显著上升趋势的为15%，不显著下降的比例为1.6%，显著下降的区域面积不足0.5%，主要为建设用地扩张。年最大植被绿度变化趋势和生长季平均状况大致相同，但程度不同，其中，绿度显著上升的面积占比为60%，不显著上升趋势的面积比例次之为34.8%，不显著下降区域的占比为4.8%，显著下降趋势不足0.4%。因此，无论是生长季平均绿度还是最大绿度，近20 a均呈现上升趋势，显著上升区占比在60%～83%。

3.2 河北坝上地区土地利用/覆被变化特征

2000—2018年的河北坝上地区土地利用分布和转移情况如图5和表1所示，自“游牧区”转变为“半农半牧区”以来，坝上地区耕地面积不断增加，目前耕地面积已达到区域土地面积的1/2以上，成为河北坝上地区主要的土地利用类型，其次是草地和林地，3种类型占到区域总面积的91%。但耕地主要分布在中西部地区，林草地主要分布在坝东地区。2000—2018年坝上耕地变幅不大(图5)，耕地面积转移主要表现在转为林地、草地和建设用地，变化的原因主要是由于农田开垦、退耕还林退耕还草、耕地撂荒和城镇扩张(表1)等。与草地转换密切的土地类型主要有耕地、林地和未利用地，但总体来说，草地面积近20 a减少了10%。同样减少的还有水域面积，但坝上河流和湖泊多属于季节性水域，由于数据来源的季节性差异导致土地利用数据的水域面积提取具有较大不确定性，但从实地调查和多源遥感数据看，近20 a河北坝上水域面积大幅萎缩。坝上林地以增加为主，主要来源于耕地和草地的退耕还林、植树造林活动。

图5 河北坝上地区2000，2010和2018年土地利用状况

表1 河北坝上地区2000-2018年土地利用转移矩阵 hm2

建设用地增幅最大，耕地和草地为主要的转入类型，转出类型主要为耕地，主要是对坝上部分空心化严重村庄，实施整体拆迁和易地扶贫搬迁，对“空心村”拆迁后的土地进行复垦，从转出类型看大部分转移为耕地，少部分转为草地。与未利用地转移较多的土地类型为耕地和草地，主要表现为农业开垦和草地沙化，另外部分湖淖干涸形成盐碱地。总的来说，2000年以来河北坝上地区土地利用变化以政策性引导为驱动，随着国家“三北”防护林建设工程、退耕还林退耕还草等生态项目的落地实施和不断推进，包括草地、林地等生态用地小幅增加，生态环境有所改善。例如尚义县南部，自2002年实施退耕还林以来，甲石河乡和下马圈乡共退耕地3 000 hm2以上[23]，2017年林草覆盖率达70%以上。

但同时，近十几年蔬菜的规模化、集约化种植加大对区域水资源的消耗，耕地面积虽然变化不大，但是耕作强度增加，从雨养农业转为灌溉农业，加大区域生态安全威胁。

3.3 区域植被绿度动态的土地利用/覆被变化归因

考虑土地利用变化，采用公式(5)—(7)分析不同时期各土地利用/覆被的面积变化和强度变化对区域植被绿度的贡献(表2)。从表2可以看出，除了2000—2010年未利用地的贡献为负之外，其他时段和类型的贡献均为正值。2000—2018年，坝上地区生长季平均NDVI的变化量中，1/2以上发生在耕地，包括耕地开垦、撂荒、退耕还林还草等措施导致的耕地面积的变化，以及种植结构变化、旱地水浇地变化和气候变化引起的耕地植被绿度的变化；其次，林地和草地的贡献率分别为28.62%和15.07%，水域、建设用地和未利用地的总贡献率在7%左右。在不同发展时期，各土地利用方式的贡献有一定差异，2000—2010年耕地的影响明显大于其他土地类型，耕地总面积和耕地绿度都有明显增加，对区域生长季平均NDVI的贡献量为0.034，贡献率达到65.74%，耕地和林地面积和绿度变化的影响占区域植被变化的90%以上。2010年后，虽然耕地植被绿度仍呈增加趋势，但由于面积出现下降，因此，耕地的影响逐渐下降。相反，在退耕还林退耕还草的政策下，林草面积不断上升，耕地、林地和草地影响趋于均等(表2)。

表2 河北坝上地区不同时期各土地利用/覆被变化对区域植被绿度变化的贡献率 %

耕地对年最大NDVI的影响较生长季平均NDVI更大，结果表明，2000—2018年的年最大NDVI的变化中57.22%来自于耕地，高于生长季平均值的50.52%，其次，仍然是林地和草地，贡献率分别为21.73%和14.41%，较对生长季平均NDVI的影响略小。不同阶段，耕地的贡献率差异明显，2000—2010年耕地的绿度贡献率达67.29%，2010—2018年降至43.74%，林地的影响较为稳定，草地的贡献率上升明显。水域、建设用地和未利用地由于面积占比较少，对区域植被绿度总趋势影响不大，但呈现明显的增加趋势，从前10 a的2.18%～3.18%增加到2010—2018年的9.62%～10.44%。

4 讨 论

作为首都和京津冀重要生态环境支撑区和水源涵养区，河北坝上地区植被状况对于其生态功能的发挥起着重要作用，因此对坝上地区的植被绿度的分析具有重要意义。本研究以MODIS MOD13Q1 NDVI遥感数据为反映植被绿度的数据源，结合Landsat土地利用数据，分析2000—2019年坝上植被绿度的年际变化趋势，并定量分析其变化的土地利用/覆被变化归因。植被绿度空间特征与时间序列变化方面，本研究的结果与前人研究结果类似[24-28]，坝上植被绿度(覆盖)和稳定性均表现为从坝西到坝东依次递增，2000—2019年坝上地区的植被绿度呈现明显的增加趋势，60%以上的区域呈显著增加趋势[27]，土地沙化程度减轻，生态状况好转，相应地，区域生态系统服务价值呈增加态势，2000—2015年生态系统服务总价值共增加334.83亿元[24]。

影响因素方面，通常研究主要分析气候变化和人类活动的影响，尤其是气候条件的波动。例如陈辉等[28]研究表明1987—2000年冀北地区土地利用类型变化、错季蔬菜生产等人为因素对NDVI的变化作用不明显，降水、特别是生长季降水分布特征是影响NDVI变化的重要因素；类似地，2000—2013年坝上地区年均NDVI也与年降水的相关性较大[29]。坝上地区作为典型的水分限制区，植被对降水的变化较为敏感。

本研究聚焦区域植被绿度发生的空间，即通过不同土地利用类型的绿度和面积变化的分解，定量得到每一种土地利用/覆被类型的变化(植被绿度和面积)对于整个区域植被绿度的影响。其中，某种土地利用方式的植被绿度变化，既有气候因素的影响也有人类活动的影响。研究发现，2000—2019年坝上区域植被绿度变化的1/2以上的贡献发生在耕地，其次是林地和草地，贡献率分别为21.73%～28.62%，14.41%～15.07%。结合土地利用强度与土地利用类型转移特征，可以得出耕地的贡献量中，除了气候原因之外人类活动影响最大的是农业种植结构变化，如传统雨养耕地到水浇地的转变导致利用强度的转变，另一方面在于耕地开垦、退耕还林还草工程、易地搬迁旧宅基地腾退复垦项目等引起的耕地面积和使用强度的变化。林地和草地植被绿度的变化类似，除了受降雨、气温等自然要素因素的影响之外，植树造林工程、退耕还林还草、人工林抚育等项目都是影响林地和草地面积和绿度变化的关键因素。2010年后，耕地的影响逐渐下降，耕地、林地和草地影响趋于均等。为建设好首都水源涵养功能区和生态环境支撑区，坝上地区仍需加快生态环境恢复，加大恢复力度，进一步增加林草生态用地面积，对“空心村”治理后的原宅基地及耕地复垦后实施还林还草，避免流转作为水浇地。

本研究结果定量的指出区域植被绿度发生的主要土地利用/覆被类型，掌握了区域绿度变化来源，为地区生态建设和京津冀生态环境支撑区的建设提供一定参考，为京津冀地区可持续发展、土地合理利用及生态环境保护提供科学依据和决策支持。但该研究仍具有一定的局限性和不确定性。一方面未能对归因进行更细致的分析，如2000—2018年耕地引起的区域植被绿度50.51%～57.22%的变化中，来自于耕地面积的变动和耕地绿度变化所占比例各为多少，耕地与其他土地利用类型之间的转移各对区域植被绿度产生多大的影响，还缺少深层次的分析。以后的研究可以结合具体的研究内容和应用领域，可以针对某地类，展开更加深入的分析。

5 结 论

(1) 2000—2019年河北坝上地区植被生长季平均NDVI为0.39，年最大值均值为0.63，空间上呈现坝东高于坝西，坝东植被稳定性较坝西高，林地植被绿度最高，其次草地、耕地、未利用地，除水域区域外，建设用地区域植被绿度最低。

(2) 2000—2019年河北坝上地区植被绿度呈现增加趋势，其中，生长季平均NDVI增速为0.044/10 a，年最大值增速更大，达到0.063/10 a，植被绿度增加的区域占总面积的95%以上，其中显著增加趋势占60%～83%。

(3) 土地利用/覆被变化归因分析表明，2000—2018年，坝上地区植被绿度变化量中，耕地贡献了50.51%～57.22%，林地和草地分别贡献21.73%～28.62%和14.41%～15.07%，水域、建设用地和未利用地的总贡献率在6%左右。其中，2000—2010年耕地植被变化对区域植被绿度变化的贡献率在65%以上，林地和草地植被贡献不足1/3；2010—2018年，耕地的贡献在降低，其他土地利用类型的贡献呈现逐步上升趋势，耕地、林地和草地的贡献趋于均等。