1999-2018年黄河流域NDVI时空变化及驱动力分析

付含培, 王让虎, 王晓军

(1.山西大学 环境与资源学院, 太原 030006; 2.山西大学 黄土高原研究所, 太原 030006)

植被是联结大气圈、土壤圈、水圈和生物圈的重要纽带，对于监测生态环境的变化具有重要指示作用[1]，因此植被覆盖的时空变化特征及驱动因子分析在区域生态环境变化研究中具有重要意义[2-3]。归一化植被指数(NDVI)[4]作为表征植被生长状况的重要指示因子，与叶面积指数、生物量存在明显的相关性，可以对地表植被生长状况进行快速评估[5-8]。基于MODIS，SPOT/VEGETATION等卫星遥感影像得到的NDVI数据具有易获取、分辨率高的特点，是长时序植被覆盖研究的重要数据来源。目前国内外众多学者已利用不同时段的NDVI数据采用不同方法对区域尺度植被覆盖的时空变化规律及对气候变化的响应进行了深入研究[9-10]。

黄河流域是我国第二大流域，生态环境较为脆弱，近年来在气候变化和人类活动等因素作用下，流域植被覆盖状况发生了转变[11]。目前，黄河流域植被覆盖变化监测及驱动机制研究取得了一定的进展。张亚玲等[12]采用1998—2012年SPOT-NDVI数据发现黄河流域NDVI平均值与纬度变化存在一定的负相关性，并随着时间增长呈现出增加的趋势。袁丽华等[13]利用2000—2010年的MOD13Q1数据发现，黄河流域NDVI呈现出北部较低，而东南部、西部较高的特征。赵安周等[14]基于MODIS-NDVI数据发现在2000—2014年期间，黄河流域NDVI保持了增长的趋势，并且达到了6.93%/10 a的增速。刘启兴等[15]利用2000—2016年的MODIS-NDVI数据，综合采用多种方法(肯德尔法检验法等)探讨了气候因子对于黄河流域植被覆盖特征的变化，研究结果显示该区域NDVI在2000—2016年保持了较高的稳定性；其中大部分区域的植被NDVI保持增长趋势，其比例达到了70.4%，而年增长率基本保持在0～0.004范围内。郭帅等[16]分析1982—2015年GIMMS NDVI数据及相关气象数据，得出黄河流域NDVI与年平均气温、年降水量以正相关为主，其中显著正相关分别为22.39%，21.99%。叶培龙等[17]通过分析1980—2018年黄河上游各气象站点的数据，发现黄河上游呈现暖湿化趋势，其有利于黄河上游植被生长。但鲜有学者对影响黄河流域植被覆盖的自然和人文因子展开系统性分析和研究，且缺少相关驱动力分析。地理探测器是探测空间分异性，以及揭示其背后驱动力的一组统计学方法[18]，已有学者将该统计方法应用于植被NDVI驱动机制的研究[19-21]。因此，本研究基于SPOT NDVI遥感数据，采用Sen趋势分析、Manna-Kendall检验、Hurst指数、空间转移矩阵分析1999—2018年黄河流域NDVI时空变化特征，并采用地理探测器模型对黄河流域NDVI空间分异特征与驱动因子进行探测分析，为黄河流域植被恢复与生态建设提供科学依据。

1 研究区概况

黄河流域位于95°50′29″—119°6′26″E，32°6′53″—41°48′18″N，流域面积为79.46万km2，地跨9省。海拔0～6 255 m，由西向东依次为西部青藏高原区，多为山地，海拔较高；中部地区以黄土高原区为主，黄土地貌，地表破碎，水土流失严重；东部黄淮海平原区，主要由黄河冲击而成，海拔较低，整体地势西高东低。黄河流域西部为大陆性干旱气候，中部为季风性半干旱气候，东部属海洋性半湿润气候。黄河流域年均温9℃，流域大部分地区年降水量在200～650 mm，西南部和东南部多于650 mm，深居内陆的西北部部分地区则少于150 mm。除气候及地形等方面的因素外，该区域存在多种类型的地貌，同样促进了不同植被的生长。

2 数据与方法

2.1 数据来源与预处理

在本次研究中采用的植被NDVI数据来自于1 km植被指数(NDVI)空间分布数据集，并通过最大值合成法形成；气象气候数据来自于气象背景数据集，基于全国1 915个站点的气象数据插值生成500 m×500 m的空间分布数据集；地貌、土壤数据分别来自于中国1∶100万地貌类型数据、土壤类型空间分布数据；90 m DEM来自于SRTM V4.1(Shuttle Radar Topography Mission，SRTM)数据，通过拼接重采样生成1 km栅格数据，并利用ArcGIS进一步获取坡度和坡向数据；土地利用类型数据主要来自于中国土地利用现状数据集，分辨率1 km；人口密度、GDP密度数据分别来源于中国人口空间分布公里网格数据集、GDP空间分布公里网格数据集，均为1 km网格。所有的数据均通过中国科学院资源环境科学数据中心(http:∥www.resdc.cn/)获取，并通过ArcGIS软件利用黄河流域范围矢量数据进行裁剪、投影变换等处理。

2.2 研究方法

2.2.1 趋势分析 采用Sen氏趋势分析与M-K检验结合，分析像元尺度上NDVI变化趋势。Sen氏趋势分析对数据误差有较强的抵抗能力，且分析结果较为科学可信[22]。其计算公式如下：

当斜率S大于0时，表示呈现增长趋势，小于0则呈现下降趋势。

Manna-Kendall检验(M-K检验)可判断趋势的显著性，且能够排除少数异常值的干扰。因此，本文将二者相结合，用来判断黄河流域内NDVI在像元尺度的变化趋势及其显著性。

2.2.2 Hurst指数 Hurst指数多用于描述自然界中长时间序列的自然现象[23]。具体公式见相关文献[13]。若通过上述公式证明时间序列中存在Hurst现象，可利用最小二乘法拟合得到H值。

NDVI序列的完全随机性和持续性由H值决定。当满足条件0

2.2.3 地理探测器 目前在空间分异性研究中广泛采用了地理探测器模型[18]。在本研究中即利用该工具分析黄河流域NDVI变化的驱动因子单个及其交互关系。已有研究表明，气候因子对NDVI有显著影响[15]。黄河流域面积广阔，地形复杂多样。西部高原山地，植被覆盖较为原始，受人类影响较小，西北部内蒙古高原及北部河套平原，历史上为农业灌溉区，受人类活动影响较大。因此，在考虑系统性、典型性等因素下，本研究选取气候、地表及人类活动三大类因素，14个因子，探测其对黄河流域NDVI的影响。在ArcGIS 10.2中，按照1 km×1 km网格，随机生成7 520个采样点，结合空间位置关系对采样点NDVI数据进行关联。在本次研究中共选取14个影响因子(表1)。本文将平均降水量、湿润指数、干燥度、积温(≥10℃)、积温(≥0℃)、年均温分为6个等级，将高程分为7个等级，将坡度、地貌、土壤、GDP、人口、土地利用类型分为9个等级，将坡向分为10级，划分过程中采用自然断点法(图1)。

表1 黄河流域NDVI影响因子

地理探测器具体包括如下4种类型：

(1) 分异及因子探测。其可以对NDVI的空间分异性进行探测，在此基础上可以进一步分析不同因子产生的影响。具体公式如下：

(2) 交互作用探测。基于交互作用探测可以对人类活动、气象以及地表等因素之间的交互作用进行分析，可以评价多因子的交互作用是否能够影响到对于NDVI的解释能力。

(3) 风险区探测。据此可以评价不同子区域属性均值的差异性，用于鉴别各区域植被覆盖情况，用t统计量来检验[20]。具体公式如下：

(4) 生态探测。生态探测需要借助F统计量，衡量两因子X1和X2对NDVI空间分布的影响是否有显著差异。

式中：L1，L2分别为X1，X2分层的数量；Nx1，Nx2分别为二者的样本数目；SSWx1，SSWx2对应着二者的分层、层内方差和[18]。

3 结果与分析

3.1 植被NDVI空间分布格局

基于1999—2018年的NDVI数据，结合黄河流域实际情况，本文利用等间距法将植被状况划分为5种类型(表2)。由植被覆盖空间格局图(图2)可见，黄河流域NDVI分布具有明显的空间分异特征，总体表现为南高北低，与纬度表现出一定的负相关性，即纬度越低，则NDVI值越高。其中，较高及以上植被覆盖度占黄河流域面积的47.13%，主要分布于流域的东部和南部地区，包括山西西部吕梁山脉、陕南及关中地区；中植被覆盖度占26.80%，多位于中部地区，呈东北—西南走向带状分布；较低植被覆盖度占21.41%，集中在流域的北部；低植被覆盖度面积仅占4.65%，分布在河套平原两侧。黄河流域西南部主要为山地，植被以森林、草地以及草本湿地为主，有较高的NDVI值；东南部属于海洋性半湿润气候，水热较充足，且林地和农作物分布较为广泛，所以NDVI值较高。北、中部为内蒙古高原、河套平原，生态环境脆弱，气温高、降水少，植被破坏严重，现有植被多为低矮的农作物，因此NDVI值较低。

图1 黄河流域因子空间分布

表2 植被覆盖类型

图2 1999-2018植被覆盖空间格局

3.2 植被NDVI时间变化趋势

3.2.1 区域尺度NDVI的时间变化趋势 由区域尺度植被NDVI的时间变化趋势图(图3)可见，1999—2018年黄河流域整体呈增加态势，平均年际变化率为0.007。20 a间，年均NDVI最高值出现在2012年，为0.44，最低值出现在2000年，为0.34。就季节来看，夏季NDVI值最高，秋季、春季其次，冬季最低。黄河流域春、夏、秋季NDVI变化趋势与年际变化趋势基本一致，呈波动上升趋势，平均年际变化率分别为0.005，0.007，0.007，而冬季NDVI的年际变化相对平稳，年际变化率只有0.000 07，这可能因为当地植被主要以落叶植被为主，冬季植被大都凋落，增长不够明显。从月份NDVI变化趋势来看，最低值出现在冬季的2月份，为0.16；NDVI最高值出现在夏季的8月份，为0.53；2—8月NVDI值上升明显，8月—次年2月下降明显，呈波峰和波谷状。

图3 1999-2018年区域尺度NDVI时间变化趋势

3.2.2 像元尺度NDVI的时间变化趋势 根据Sen趋势分析结果将NDVI变化趋势分为7级(图4)，根据M-K检验结果(0.05置信水平)将NDVI变化显著性分为6级(图5)。由变化趋势结果可见，黄河流域20 a间NDVI时间变化趋势以极显著缓慢增长为主，整体变化趋势向好。其中，极显著缓慢增长区域所占面积最大，占整个流域面积的53.12%；不显著变化的区域其次，占29.98%；极显著较快增长的区域占15.47%；快速增长的占1%；而极显著较快退化和极显著缓慢退化合计共占黄河流域面积的1.43%。极显著较快增长和极显著快速增长的区域主要位于黄河流域的东部，黄河两岸，地理位置上属山西西部吕梁山脉和陕西北部黄土高原，这与该区域持续退耕还林，人工在造林有密不可分的关系，一定程度上减少了黄土高原的水土流失，改善了黄河生态环境。基本不变的区域则位于黄河流域北部，河套平原周边地区以及黄河下游地区，该区域城镇、耕地面积较大，因而基本没有变化。

图4 1999-2018年NDVI趋势变化

图5 1999-2018年NDVI变化显著性

3.2.3 NDVI变化的可持续性分析 从黄河流域NDVI的Hurst指数(图6)可以看出，H值低于0.3的区域零星分布在流域的中部及南部地区，仅占流域面积的1.46%；H值大于0.7的区域位于中部和北部地区，占流域面积的0.67%；0.5～0.6的区域占28.78%，0.6～0.7占6.32%，主要分布于山西西部和陕西中部；0.3～0.4的区域占流域的18.13%，0.4～0.5的区域占黄河流域的45.15%，其反退化占比较高，主要分布于甘肃陇中黄土高原、宁夏地区和陕西、山西黄土高原区，这些区域本身为黄土高原或干旱沙地，植被覆盖差，多为季节性降水，加之人类活动频繁，其植被生长持续性不强。0.5～0.6的区域占流域面积的28.78%，0.6～0.7占6.32%。

本文将Hurst指数图与趋势分析图在ArcGIS中叠加，得到黄河流域NDVI未来趋势变化图(图7)。结果表明，持续快速退化和持续较快退化面积极小，仅占0.01%，几乎可以忽略不计；持续较缓慢退化的面积占0.14%，大多分布于城镇及周边，是城区扩张造成的；持续基本不变的面积占30.19%；波动变化占整个流域面积的44.86%，是占比最大的区域；缓慢退化区和缓慢增长区域是波动变化的主要发生地；持续缓慢增长占17.67%，与Hurst指数0.5～0.6分布基本一致；持续较快增长占6.66%，持续快速增长面积仅占0.47%，二者占比较小，多分布于持续缓慢增长区内。

图6 1999-2018年NDVI Hurst指数

图7 未来变化趋势

3.3 植被NDVI空间分布变化

在ArcGIS中通过叠加分析得到了1999—2018年植被覆盖空间转移矩阵表(表3)，结果表明，1999—2018年期间高植被覆盖度和较高植被覆盖度大幅增加，中植被覆盖度增加明显，主要是由低一级植被覆盖度转化而来。演化趋势为：较低植被覆盖度→中植被覆盖度，中植被覆盖度→较高植被覆盖度，较高植被覆盖度→高植被覆盖度。20 a来，植被覆盖度保持稳定的区域占31.6%，持续增加的占65.99%，减少的占2.41%。总体来看，黄河流域植被增长明显，生态恢复成效显著，生态环境整体向好。

表3 黄河流域植被覆盖空间转移矩阵 %

3.4 NDVI空间分布的地理探测

3.4.1 不同驱动因子对NDVI空间分异的驱动分析 用地理探测器对NDVI的空间分异及变化趋势进行归因，因子探测结果见表4，不同因子对黄河流域NDVI空间分布的解释能力如下：年平均降水量>湿润指数>干燥度>土壤类型>年均温>土地利用>积温(≥0℃)>积温(≥10℃)>高程>地貌>坡度>人口>GDP>坡向，这也反映了地处内陆半湿润、半干旱、干旱交汇区的黄河流域，水分为主的气候因素是影响NDVI的主导因素。其中，年平均降水量、湿润指数、干燥度、土壤类型的解释力均超过30%，分别为59.6%，46.1%，34.9%，33.9%，是主要影响因子；年均温、土地利用、积温(≥0℃)、积温(≥10℃)、高程、地貌解释力大于20%，分别为27.9%，26.9%，26.2%，25.7%，24.2%，20.7%，是次要影响因素；坡度、人口、GDP、坡向解释力均低于20%，影响力较低。

表4 因子探测结果

3.4.2 交互作用探测分析 由表5可知，交互作用解释力最强的是平均降水量和土地利用，它们双因子交互q值达到了0.704，这说明了土地利用的变化影响降水量,并以此作为因变量解释黄河流域NDVI较强。解释力最差的是坡向与GDP的交互，这可能是由于本身坡向及GDP对NDVI解释力较差，在这样两个因子交互下，解释力不足。同时，降水量与大多数因子相结合都能产生较高的q值，最低值是平均降水量与坡向的组合，却也达到了0.638。坡向、GDP、人口密度与多数因子交互解释力较低。

表5 交互作用探测结果

由表6可知，绝大多数因子交互呈现相互增强或非线性增强，表现为两因子q值相加大于单因子q值。其中以年平均降水量为代表的主导因子表现更为明显，X4∩X12(0.704)>X4∩X11(0.698)>X4∩X10(0.688)>X4∩X7(0.687)>X4∩X6(0.669)>X4∩X1(0.665)>X4∩X2(0.654)>X4∩X8(0.640)>X4∩X14(0.632)>X4∩X5(0.627)>X4∩X13(0.624)>X4∩X3(0.615)>X4∩X9(0.609)，结果表明，年平均降水量与土地利用类型、土壤类型、地貌、高程、年均温等因子的相互作用呈现相互增强的关系。X9∩X4(0.609)>X9∩X5(0.479)>X9∩X3(0.368)>X9∩X11(0.359)>X9∩X6(0.295)>X9∩X12(0.290)>X9∩X1(0.279)>X9∩X2(0.276)>X9∩X7(0.254)>X9∩X10(0.224)>X9∩X8(0.128)>X9∩X14(0.114)>X9∩X13(0.047)，且X9∩Xn均大于X9+Xn，可见坡向与年平均降水量、湿润指数、干燥度、土壤类型等因子的相互作用效应呈现非线性增强的关系。

研究结果表明，仅有两对因子交互对NDVI的影响是相互独立的，分别是年平均降水量与土地利用类型交互及高程与人口密度的交互(表6)，X4∩X13=0.624=X4+X13，X7∩X14=0.339=X7+X14。

3.4.3 生态探测分析 由表7可知，气候要素中各因子大多与其他气候要素因子存在显著性关系，年平均降水量尤为突出，而与地表因素和人类活动因素极少存在显著性关系。比如积温(≥0℃)与干燥度、年平均降水量、湿润指数和土地利用类型存在显著性关系，年平均降水量与积温(≥0℃)、积温(≥10℃)和干燥度存在显著性关系。在地表因素中，大多数因子均与土壤和土地利用类型之间存在显著性差异，比如高程与土壤类型和土地利用类型存在显著性关系；坡向与地貌、土壤类型、土地利用类型和人口密度有显著性差异；地貌与坡度、坡向、土壤类型和土地利用类型存在显著性差异。由此可见，NDVI受三大类要素共同影响，气候因素仍为主导。其中，年平均降水量、土壤和土地利用类型对NDVI影响最大，GDP对NDVI影响最小。

表6 影响植被NDVI主导因子交互作用分析

表7 生态探测结果

3.4.4 影响因子的适宜范围或类型 基于风险探测，我们能够得出各类影响因子中，植被NDVI最适宜的范围(表8)。气象因子方面，年平均降水量在658～1 039 mm，年均温-2.6～1.8℃时，NDVI均值最高，可见水热条件较好的地区植被覆盖度更好。高程在3 216～4 005 m，坡度在15.29°～27.56°时植被覆盖度最高，与地貌的大起伏中山、大起伏高山，土地利用的林地、草地相对应，说明在山高坡陡的地区，植被覆盖度更高。草甸土和砂姜黑土广泛分布于山地和山前洼地中，从侧面反映出海拔对植被覆盖的影响。

表8 不同影响因子的适宜类型或范围

总体来看，由风险探测我们可知，气象因素、地表因素是影响黄河流域NDVI的主要因素，而人类活动对NDVI的影响主要体现在土地利用方式的改变上，退耕还林与植树造林工程发挥了重要作用。

4 讨 论

黄河流域地跨我国三级阶梯，地貌类型复杂多样，自然与人为因素共同作用，对流域内植被覆盖产生重要影响。本研究表明黄河流域自实施退耕还林、生态修复工程实施以来，流域NDVI变化趋势以极显著缓慢增长为主，植被覆盖空间格局整体向好，而气象因素、地表因素与人类活动因素共同影响了植被分异，其中气象因素依旧为主导因素，平均降水量与湿润指数解释力均大于45%，占主导地位。由此可见降水是植被活动的主要驱动力，这与前人的研究结果一致[24-25]。土壤类型是植被NDVI变化的重要因子[20]，在以降水为主要限制因素的地区，土壤类型对植被生长和雨水再利用效率影响显著，如淋溶土、半淋溶土地区的NDVI值较高；而漠土与干旱土NDVI值较低，不适宜植被生长；初育土广泛分布于陕北黄土高原区，该区是黄河流域生态修复与水土保持的重点区域，因此NDVI增速较快。人类活动中，土地利用方式对NDVI影响较大，自1999年退耕还林还草工程实施以来，黄河流域NDVI持续增长，但增速较慢，这可能与城市化和过度放牧有关。

气候因素与地表因素、人类活动因素相交互，会显著提高对NDVI的解释力。本文研究发现，降水量与土地利用解释力最强，q值为0.704，降水量与其余各因子交互q值均大于0.6，湿润指数与其余各因子交互q值均大于0.47，也反映了水分的主导作用。这与相关研究[26]认为的社会经济因子与自然因子交互，会显著提高其对植被NDVI空间分布的影响力一致。

像元尺度上的区域NDVI时空变化及其驱动因素分析一直是地理学和生态学研究的重点，本文利用SPOT NDVI数据可观揭示了黄河流域20 a间植被覆盖度变化的时空特征，本研究表明黄河流域NDVI变化的可持续性不强，主要以波动变化为主，与其他已有研究的相同时期对比，基本趋势是一致的[13]。已有研究表明，SPOT系列卫星NDVI数据在光谱波段设计、空间精度和几何纠正等方面都比其他数据更有优势，因此其应用也较为广泛[27]。因此，本研究是可信的。但是NDVI存在饱和现象，特别在夏季或植被覆盖度高的林区尤为明显。EVI能很好地反映植被生长，即使在生物含量较高的月份也能反映出植被生长的季节性特征[28]，因而EVI更适用于长序列季节性研究。本文大多基于年际NDVI数据展开分析，因此受饱和现象影响较小。

地理探测器是探测空间分异性，以及揭示其背后驱动力的一组统计学方法[18]。地理探测器在各领域应用广泛，特别是在对NDVI的探测上。本文采用地理探测器对黄河流域NDVI空间分异特征展开探测，结论与同流域研究基本一致[29]。但是，在选取因子指标和对各因子数据进行离散化处理时，存在一定的主观性。

5 结 论

(1) 1999—2018年黄河流域植被覆盖整体较好，多年平均NDVI值为0.39，呈现南高北低的分布格局，其中较高及以上植被覆盖度区域占流域面积的47.13%。从NDVI空间变化特征来看，1999—2018年期间，高植被覆盖度和较高植被覆盖度大幅增加，中植被覆盖度增加明显，主要是由低一级植被覆盖度转化而来。

(2) 20 a来黄河流域NDVI整体呈增加态势，年增率为0.007；季节变化特征明显，最低值出现在冬季的2月份，NDVI值为0.16，最高值出现在夏季的8月份，NDVI值为0.53。从NDVI像元尺度的时间变化趋势来看，1999—2018年黄河流域NDVI变化趋势以极显著缓慢增长为主，占整个流域面积的53.12%，整体变化趋势向好，但可持续性不强；极显著较快增长和快速增长的区域占15.47%，主要位于山西西部吕梁山脉和陕西北部黄土高原，这与当地生态恢复工程紧密相关。需要注意的是，未来持续波动区域所占面积比例达到44.86%，其中以缓慢退化和缓慢增长为主。

(3) 黄河流域的NDVI空间差异主要由年平均降水量、湿润指数、干燥度、土壤类型决定，其因子解释力均超过30%，这表明气候因素仍然是影响黄河流域NDVI的主导因素。绝大多数驱动因子间交互呈现相互增强或非线性增强，交互解释力最强的是年平均降水量和土地利用，q值为0.704。只有年平均降水量与土地利用类型交互及高程与人口密度的交互是相对独立的。