生态工程背景下黄土高原植被变化时空特征及其驱动力

修丽娜， 颜长珍, 钱大文， 幸赞品

(1.西北师范大学 地理与环境科学学院， 甘肃 兰州 730070; 2.中国科学院 西北生态环境资源研究院，甘肃 兰州 730000； 3.中国科学院 西北高原生物研究所， 青海 西宁 810008)

黄土高原作为我国历史上开发时间较早、人工改造力度较大的区域之一，深受人类活动干扰的影响，尤其是过去几十年来，毁林开荒、过度放牧等行为严重的破坏当地的生态环境[1-2]，再加上黄土高原地处半干旱半湿润气候带，风力和水力侵蚀严重，导致了水土流失、沙漠化加剧等一系列的生态环境问题[3-4]。针对这些生态环境问题，1990年以来，我国开始实施退耕还林还草的一系列生态恢复工程，而黄土高原则是重点实施区域之一[5-6]。植被是陆地生态系统的重要组成部分，在陆地表面能量交换、水分循环和生物地球化学循环过程中起着至关重要的作用[7]。植被的时空变化会影响全球碳循环和气候系统的稳定性，监测植被的时空变化趋势对研究全球气候变化背景下的陆地生态系统具有重要的意义[8]。一方面，植被变化是区域生态变化的有效指标，其时空变化特征可以作为生态恢复工程的评价依据，另一方面，植被改善能够通过改变地表反射率和下垫面粗糙度以达到减缓径流冲蚀、风力侵蚀和保持水土的作用，对区域生态环境起促进作用[9]。因此黄土高原地表植被的变化研究对于该地区环境保护和生态工程建设具有重要的指示意义。

归一化植被指数(normal difference vegetation index, NDVI)是目前最常用表示植被生长状况的指标，它与植被覆盖度、生物量和光合作用强度等都密切相关[10-11]。很多学者利用NDVI对于黄土高原地区植被状况的研究已经做了大量的工作，并取得了丰富的成果。高健健等[12]利用1981—2012年的NDVI数据遥感反演黄土高原的植被覆盖度，并发现2001年以后该区生长季平均植被覆盖度大于之前，呈现增加趋势。李双双等[11]发现2000—2009年陕甘宁黄土高原地区NDVI增速远快于三北防护林工程区，这一结果主要是有人类活动和气候变化共同驱动，且人类活动作用明显。肖强等[7]通过NDVI数据估算了黄土高原植被覆盖度，其研究结果发现2000—2010年植被覆盖度总体呈现东南高西北低，由东南向西北递减的特征。然而，这些研究通常使用单一数据集来评估相对较短时期内或较粗分辨率的植被动态，导致缺乏对植被变化的长期监测，使这些研究的代表性有限。此外，NDVI作为表征植被生长状况的指标，以往研究仅从NDVI单方面对黄土高原地区植被进行评价，无法获取全面的植被变化信息。植被生长受到环境驱动因素(例如气候变化)和人类活动的影响，在已有研究中区分这些因素的方法主要分为3种： ①利用气候与NDVI之间的相关性，并辅助相应的统计数据间接的区分气候因素和人类活动对NDVI变化的影响[4,11]; ②使用水分利用效率变化模型提取人类影响，该方法只适合于植被变化主要受降水因素影响的干旱和半干旱区域[13]； ③利用回归残差法能较好的区分气候因素和人类活动对植被的影响，目前应用较为广泛[1,14]。但此方法只是单纯的得出人类活动对植被NDVI的正或负面影响，对于人类和气候共同影响的区域没有明确得出。为此，本研究利用GIMMS NDVI和MODIS NDVI相关关系将MODIS NDVI时序延伸，获得长时间序列较高分辨率的NDVI数据集，从植被质量角度，以生长季NDVI作为指标对1990—2015年黄土高原植被的时空演变规律进行分析，辅助土地利用数据从植被面积变化角度全面分析黄土高原植被时空变化特征，并从气候变化和人类活动两个方面分析本地区植被变化的驱动力，以期为生态环境保护相关政策的制定、区域可持续发展提供理论依据。

1 材料与方法

1.1 研究区概况

黄土高原位于我国中部偏北，黄河流域中部，地理位置为33°41′—41°07′N, 100°54′—114°07′E。黄土高原地跨山西、宁夏、陕西、甘肃、内蒙古、河南和青海省，总面积为6.24×105km2。区域内地势西北高、东南低，千沟万壑、地形破碎。气候为大陆性季风气候，冬季寒冷干燥多风沙，夏季炎热多暴雨。植被由东南向西北可划分为森林带、森林草原带、典型草原带、荒漠草原和草原化荒漠带。土壤类型主要有黄绵土、黑钙土等[12,15-16]。

1.2 研究数据

1.2.1 NDVI数据来源及处理 本研究利用的遥感数据为NOAA/AVHRR和TERRA/AQUA MODIS的NDVI数据集，NOAA/AVHRR NDVI数据集来自戈达德航天中心(Goddard Space Flight Center， GSFC)全球监测与模型研究(Global Inventor Modeling and Mapping Studies， GIMMS)数据中心。影像以32bit geotiff格式存放，采用Albers Conic投影方式，空间分辨率为8 km×8 km，时间分辨率为15 d，时间序列从1990年1月至2006年12月(可在ftp.glcf.umiacs.umd.edu免费下载)。MOD13A3是MODIS陆地产品组开发的月植被指数产品，包括月最大合成的NDVI和EVI植被指数，空间分辨率为1km，可以通过NASA对地观测系统数据共享平台(https:∥earthdata.nasa.gov/)下载。

已有研究[17]发现GIMMS和GIMMS 3g数据集在北半球存在较大差异，因此建议可以结合GIMMS和MODIS数据对该区域的植被变化进行分析。通过GIMMS NDVI和MODIS NDVI相关关系建立将MODIS NDVI时序延伸。首先将15 d GIMMS NDVI数据采用最大值合成法生成月NDVI数据集，然后将MODIS 1km月NDVI(MOD13A3)数据集按照GIMMS NDVI的8km格网分别计算每个格网的月NDVI均值和最小值。计算最小值主要是为排除非植被的区域做数据准备，求得的均值是与GIMMS NDVI数据做匹配。选取研究区内，2000—2006年7年间逐月最小值均>0的格网，分别计算7年每个月GIMMS NDVI同MODIS NDVI均值之间的线性相关系数[18-20]，结果详见表1。选取最小值大于0的栅格主要是筛选出植被分布区域，由于1 km分辨率相较8 km的空间分辨率较精细，便于排除非植被区域NDVI值参与计算。然后依据此关系，使用1990—1999年GIMMS NDVI扩展MODIS NDVI的时间序列。为了反映植被NDVI的年际变化特征，本研究中采用生长季NDVI值(4—10月NDVI的平均值)来表征植被生长[8]。

表1 GIMMS NDVI与MODIS NDVI之间的相关关系

1.2.2 气象数据来源及处理 气象数据来源于中国气象局提供的中国气象科学数据共享服务网(http:∥cdc.cma.gov.cn/home.do)，包括1990—2015年黄土高原及其周边地区共84个气象台站观测的日平均气温和日降水量数据。本文利用ANUSPLIN 4.3软件对整理好的月降水和温度进行空间插值，空间分辨率为1 km。该方法已经广泛应用于气象因素的插值[21]，其精度要高于反距离权重(inverse distance weighted, IDW)和克里金插值(Kriging)[22]。

1.2.3 土地利用数据的解译 本研究中用于土地利用解译的影像包括没有云覆盖的1990，2000，2005，2010的TM影像和2015年的OLI影像，这些数据下载于美国地质调查局(United States Geological Survey，USGS)网站(http：∥glovis.usgs.gov/)。大多数选定的图像是从6月到9月获得的，丢失的图像或质量差的图像被替换为同期最近的月份或年份的图像。在进行遥感解译之前，影像已经通过几何校正，辐射校正和大气校正等预处理。本研究中使用的6种土地利用类型分别是：林地，草地，耕地，人工表面和其他类型，其中其他类型包括裸岩、戈壁、裸土、沙漠和盐碱地[23]。通过eCognition 8.6软件从影像中提取土地覆盖信息，该软件使用面向对象的分类方法[24-25]。

1.3 研究方法

1.3.1 变化趋势分析 为了研究过去26 a来黄土高原植被的生长状况及在空间上的差异，使用一元线性回归对每个像元进行斜率分析[7]。当斜率为正时，表示随时间变化植被指数升高；反之，斜率为负时，表示随时间变化植被指数呈现下降趋势。对所得的斜率进行显著性水平检验，通过置信度检验的区域则认为植被变化的趋势显著。其计算公式为：

(1)

式中：n——模拟的时间序列长度(n=26)；i——时间序列的序号； NDVIi——第i年的NDVI值；Θslope——变化曲线的斜率。

1.3.2 土地利用类型净变化率 年均土地利用类型净变化率定量反映该土地利用类型面积变化速率，其对比较各土地利用类型变化的区域差异具有重要作用[26]。其计算公式为：

(2)

式中：K——某土地利用类型年均净变化率；Ua，Ub——研究初期与末期某一种土地利用类型的面积；T——研究期时长，时段设定为年。

1.3.3 土地利用类型转移矩阵 利用土地利用类型面积作为土地利用类型状态转移矩阵中的向量，分别统计1990—2000年和2000—2015年2期的转移矩阵[27]，计算过程为：

(3)

式中：A——面积；i，j(i，j=1,2，…，n)——转移前与转移后的土地利用类型；Aij——土地利用类型从类型i变为j的面积；n——转移前后土地利用类型数。

1.3.4 多元回归分析 气候因素(温度和降水)之间通常是相互作用的，并且对植被的影响不是对立的。本研究利用多元回归分析方法来研究温度、降水对NDVI的共同影响[28]。

(4)

利用F检验对上述多元线性回归方程进行显著性检验。

1.3.5 改进的残差趋势分析 残差趋势法由Evans和Geerken提出[13]，该方法利用每个像元全时间序列的年NDVI与降水数据建立响应关系，从而得到每个像元在全时间序列上实际NDVI与预测NDVI的残差，从而判断人类因素对植被生长状况的影响[29]。曹鑫等[30]考虑到温度对植被生长的影响，在已有模型基础上加入了温度因子。

ΔNDVI=NDVI真实值-NDVI模拟值=

(5)

为了研究更深入地研究气候和人类活动对于植被变化的影响程度，田海静[31]等在以上模型的基础上，将年度NDVI的变化趋势与气候变化之间相关性不显著的区域相对应。 ①由气候变化引起的植被显著恢复。 NDVI随时间的变化趋势显著增加，NDVI与气候因素显著相关且NDVI残差趋势变化不显著; ②由气候变化引起的植被显著退化。NDVI随时间的变化趋势显著下降，NDVI与气候因素显著相关且NDVI残差趋势变化不显著； ③由人类引起的植被显著恢复。NDVI随时间的变化趋势显著增加，NDVI与气候因素不显著相关且NDVI残差趋势变化显著增加; ④由人类引起的植被显著退化。NDVI随时间的变化趋势显著下降，NDVI与气候因素不显著相关且NDVI残差趋势变化显著下降; ⑤气候因素和人类活动共同引起的植被显著恢复。NDVI随时间的变化趋势显著增加，NDVI与气候因素显著相关且NDVI残差趋势变化显著增加； ⑥气候因素和人类活动共同引起的植被显著退化。NDVI随时间的变化趋势显著下降，NDVI与气候因素显著相关且NDVI残差趋势变化显著下降。

2 结果与分析

2.1 黄土高原地区植被质量时空变化

2.1.1 黄土高原地区植被NDVI年际变化特征 由图1可以看出，1990—2015黄土高原地区植被的NDVI总体表现为上升趋势，且变化趋势较明显的分为两个时期：1990—1999年期间，黄土高原地区的NDVI值波动较大，在0.3304～0.367 0之间，有轻微的减少趋势，特别是1998和1999年；2000年以后，该区域的NDVI表现为增加趋势，范围为0.286 4～0.372 3，最大值出现在2012年。

图1 黄土高原1990-2015年生长季NDVI年际变化

2.1.2 黄土高原地区植被NDVI空间变化特征 为了便于分析1990—2015年黄土高原NDVI变化趋势的空间分布规律，对这26 a以来生长季的NDVI逐像元进行一元线性回归分析(图2)，将所得相关系数进行0.05显著性水平检验，所得变化趋势分布图(图3)。黄土高原地区植被NDVI呈现增加趋势的区域面积较大，主要分布在黄土高原的东部、中部和南部地区，包括陕西省的北部，甘肃和宁夏地区南部，内蒙古自治区的东部，以及山西省的西部，其他区域的植被也有不同程度的增加；NDVI降低的区域主要分布在宁夏和内蒙古的北部，陕西的南部，以及山西的部分区域。年均增长率最大的像元值达到0.020 0，年均负增长率最大的像元值达到-0.027 4(图2)。

图2 黄土高原1990-2015年生长季NDVI变化趋势

从黄土高原NDVI显著性变化趋势分布图(图3)可以看出，1990—2015年黄土高原植被NDVI整体呈现的上升趋势，显著上升的区域主要集中在中部和东部地区。研究结果表明(表2)，1990—2015年期间，黄土高原地区53.45%的地区，植被NDVI呈上升趋势，其中33.05%的地区植被为显著上升。植被NDVI下降的区域主要分布在甘肃和宁夏地区的北部，陕西省的南部和山西省的东部等地，有46.55%的地区植被NDVI呈下降趋势，其中24.61%的地区植被NDVI呈显著下降。1990—1999年，黄土高原的植被NDVI总体呈现降低趋势，上升区域面积占总面积的44.92%，其中显著上升仅占6.45%；2000年以后，相较之前，黄土高原植被NDVI迅速上升，上升区域面积达到总面积的91.90%，显著上升为65.78%。

图3 黄土高原1990-2015年生长季NDVI显著性变化分布

表2 各时间段黄土高原生长季NDVI变化类型像元个数及百分比

2.2 黄土高原地区植被面积变化特征

黄土高原地区植被面积总体表现为增加趋势(表3)，林地主要分布在黄土高原的南部和东部地区，陕西和山西省境内，草地分布于黄土高原的大部分区域(附图2)。

林地自1990年开始一直表现为持续增加，且2000年之后增加速度较快，2000—2005年林地的净变化率达到0.31%，之后减慢，但是仍快于2000年以前，直到2010，林地面积开始有轻微减少。草地面积2000年以前总体表现为减少，净变化率为-0.10%，之后表现为与林地相同的变化趋势，2000—2005年面积的增加速度较快。

表3 黄土高原植被面积变化

1990—2000年增加的林地主要来自于耕地的转入，面积为470.73 km2,其次为草地和其他，面积分别为68.46和68.86 km2；减少的草地主要转化为耕地，面积达到4 104.05 km2, 其次为人工表面，面积为270.87 km2(表4)。随着人口数量的增多，人工表面的迅速扩张，黄土高原地区的开荒行为持续增加，导致大量的草地被开垦为耕地，或退化为裸土等其他类型。

2000—2015年，增加的林地和草地主要来自耕地的转入，面积分别为3 723.29和14 278.59 km2；此外还有大量的裸土等其他类型转入草地，面积达到1 197.67 km2(表4)。尽管城镇化进程持续加剧的行为会破坏植被，但是也导致农村农民对土地、山林的依赖性降低，大量农田撂荒，砍伐减小。加之随着退耕还林(草)生态工程的实施，大面积的退耕还林(草)、植树造林等一系列措施的实行都会导致黄土高原的植被面积持续增加。

表4 黄土高原不同时期土地利用类型转移矩阵 km2

2.3 植被对气候变化和人类活动的响应

通过改进的残差趋势分析得到，人类和气候共同引起植被恢复的区域主要分布在黄土高原的中部，主要集中在陕西的北部，占总面积的9.41%；人类引起植被恢复的区域主要分布在内蒙古的东部和北部、甘肃和宁夏的南部，以及陕西和山西省的中部地区，占总面积的21.74%；气候引起植被恢复的区域仅占总面积的1.91%，主要分布于青海省。人类和气候共同引起的植被退化主要分布在内蒙古、山西和宁夏的北部，占总面积的4.71%；人类引起的植被退化面积占11.99%，主要分布于内蒙古的西北部、山西和陕西省的部分地区；气候引起的植被退化主要分布于青海和甘肃省(图4，表5)。

表5 黄土高原1990-2015年生长季NDVI变化驱动因素解析

图4 黄土高原1990-2015年植被 显著恢复和显著退化区域分布

对黄土高原植被主要受气候影响的区域进行进一步分析。气候引起植被变化区域的温度和降水有着相同的趋势，都呈上升趋势(图5)。气温在过去26 a间升高幅度较大，且植被退化区域比恢复区域温度上升较快。降水变化趋势表现为轻微上升，且植被恢复区域较退化区域明显。

图5 1990-2015年气候引起植被恢复和退化区域的年均气温和年降水变化趋势

3 讨论与结论

1990—2015年，黄土高原地区植被的NDVI总体表现为上升趋势，这一趋势与已有的研究结果一致[8,32]。具体来说，1990—1998年黄土高原的植被处于相对稳定的时期，之后植被NDVI迅速下降，2002年植被NDVI进入迅速上升时期，这主要是由于随着气候干旱化趋势发展，植被NDVI表现为减少且幅度小，1999—2001年降水明显偏少(图6)，导致植被NDVI减少，自2002年以来，随着降水量的恢复和退耕还林还草政策的大规模实施，植被NDVI呈现明显上升趋势[4]。自2000年起，黄土高原植被面积表现为增加趋势，且主要来自于耕地的转入，这说明退耕还林(草)工程的实施，在黄土高原地区效果明显，不仅植被面积持续增加，而且生长状况也得到了改善。但值得注意的是，城镇化进程的加剧，虽然会破坏植被，但是也导致农村农民对土地、山林的依赖性降低，大量农田撂荒，砍伐减小，进而对黄土高原植被的恢复起到积极作用。

图6 黄土高原地区1990-2015年生长季 年均温度和年降水变化图

26 a间黄土高原植被NDVI整体呈现的上升趋势，显著上升的区域主要集中在中部和东部地区。结合人类残差分析的结果得知，黄土高原植被恢复的原因主要是气候和人类共同影响，且人类影响的程度较大，区域占总面积的21.74%，本研究中大量耕地转入草地印证这一结果。结合李双双等[11]的研究可知，陕北地区2000—2009年退耕还林(草)面积为9.81×105hm2，榆林、延安退耕还林还草面积分别为3.09×105hm2，2.91×105hm2，也说明中东部地区的植被恢复主要得益于生态工程的实施。人类引起植被退化的区域主要是由于城镇化扩张。罗娟等[33]对毛乌素沙地近10 a的土地利用动态变化进行了分析，发现居民地和工矿用地面积增加，且主要是由草地转入；马瞳宇等[34]选取黄土高原风蚀水蚀交错区典型代表流域进行分析，1990—2010年期间，该地区煤矿开采活动活跃，影响流域内林草生长，并对侵蚀环境的人为恶化起着巨大的促进作用；卞鸿雁等[15]发现1980—2005年黄土高原南部地区建设用地净增了1 238.29 km2。除了以上原因，对比黄土高原土地利用分布图可以看出，部分人类引起植被退化区域是耕地分布的区域，由于研究中所使用的遥感数据分辨率有限，不可避免的会将一些区域耕地的NDVI值的变化带入模型中，随着退耕政策的实施，导致耕地的减少，近而影响到最后的模拟结果。气候因素主导植被退化区域主要是由于气温升高较快，但降水增加不明显，加速了相应区域气候的暖干化，抑制了植被的生长。以上研究结果也侧面的验证了人类残差分析在黄土高原的有效性。

尽管本研究中使用的NDVI相较其他长时间序列的NDVI产品具有更高的分辨率，但1 km的空间分辨率对于黄土高原地区而言仍然相对较低。在一些复杂的地形中，坡度和其他地形因素对植被的反射率有一定的影响，进而影响NDVI的观测值。此外，黄土高原所处的地理位置，使得当某一地区的植被稀疏时，植被变化趋于退化，这些退化或暴露的斑点可能会影响遥感监测的反射率，进而影响NDVI。

残差趋势分析法被广泛用于分离气候和人为因素对植被变化的影响[35]。残差法以像素为分析单元，提出了一种充分考虑坡度，土壤和植被空间变异性对温度—降水—植被NDVI的影响的分析方法。然而，这种方法也有一些缺点[36]，潜在的降水—温度—NDVI关系的趋势分析是基于相同的时间序列。由于遥感数据的限制，残差趋势法只能检测到遥感数据时间序列内的退化。考虑到影响植被NDVI变化的因素还有很多，比如蒸散量、土壤水分等，这些因素都需要在未来的研究中仔细考虑，进而改进模型。