山西省植被NDVI时空变化及驱动力研究

窦永静，王让虎，付含培，张越，刘瑞杰

（1.太原师范学院 地理科学学院，山西 晋中 030619；2.山西大学 黄土高原研究所，山西 太原 030006；3.黄土高原生态恢复山西省重点实验室，山西 太原 030006；4.山西大学 环境与资源学院，山西 太原 030006）

0 引言

植被在联系大气圈、土壤圈、水圈和生物圈的过程中起着枢纽的作用，其与各种自然环境要素密切相关，对生态环境的变化高度敏感［1］。因此地表植被覆盖的时空变化特征及影响机制是全球生态变化研究的热点领域［2-4］。卫星遥感影像具有覆盖面积广、时空分辨率高、获取难度低等优势，使得人类能在更大范围的宏观尺度上对长时间序列的植被覆盖变化进行监测［4］。其中，归一化植被指数（NDVI）是常见植被指数中的一种，可作为检测植被生长状态的重要指示因子。目前国内外学者基于MODIS、SPOT/VEGETATION、Landsat等卫星遥感影像利用NDVI指数开展了大量的植被覆盖变化监测及驱动机制研究［5-9］。

山西省位于黄土高原东部，拥有海河流域和黄河流域两大水系，被称为“华北水塔”，也是三北防护林的重要组成部分，是京津冀地区的水源涵养地和风沙源生态屏障，但目前也面临着生态脆弱、水土流失等问题。近年来，已有学者针对山西省的植被覆盖时空变化特征及影响因素开展了研究。刘庚等［10］以MODIS遥感影像为数据源，对山西省2005−2006两期植被指数进行计算发现，全省植被覆盖度处于中等水平的区域有所增加，但植被覆盖度较高的区域增幅并不明显；武永利等［11］利用山西省1982−2006年的 NASA/GIMMSNDVI数据，得出25年来植被呈现总体变好的趋势，NDVI受降水的年际变化影响最大；张茹等［12］利用1981−2011年的遥感数据，研究了人类活动和气候因子对山西省植被覆盖的影响，结果表明温度、降水等气候因子对植被覆盖影响较大；卫宇婷等［13］发现2000年以后山西省的植被变化对气候变化的响应更为敏感。但目前已有研究多以气候因素为主要驱动因子，未对地貌、地形、土壤等其他非气候因子的影响进行深入讨论，山西省植被覆盖空间分异及变化的驱动机制尚不明确［14］。地理探测器是探测空间分异性，以及揭示其背后驱动力的一组统计学方法［15］，已有众多学者将其应用在NDVI时空变化的驱动力研究中［16-17］。因此本文基于SPOT NDVI 遥感数据，采用趋势分析、空间转移矩阵分析了2000−2015年山西省NDVI 时空变化特征，并采用地理探测器模型对影响植被变化的气候因素与非气候因素进行探测分析，明确了山西省NDVI空间分异及变化的主导驱动因素，全面地分析了山西植被动态变化规律及驱动机制，探讨了山西省自2000年以来实施退耕还林等生态恢复工程的成效，为山西植被恢复与生态建设提供科学依据。

1 数据来源及研究方法

1.1 研究地区概况

山西省国土面积为156 748.01 km2，位于我国 黄 土 高 原 东 部（北 纬 34°36′～40°44′，东 经110°15′～114°32′），属于水土 流失严重地区。全省以温带季风气候为主，年降水量大约为400 mm～620 mm，夏季炎热多雨，暴雨频发，冬季寒冷干燥。地势地貌复杂多样，总体特征为两山夹一川，山地丘陵面积占全省国土面积的67.52%。山西省第一次全国地理国情普查结果显示，全省植被覆盖（包括种植土地、林草覆盖）面积为144 801.59 km2，占国土面积的92.38%。

1.2 数据来源及处理

2000−2015年NDVI数据来自于中国植被指数（NDVI）空间分布数据集（空间分辨率为1 km），其基于SPOT/VEGETATION 卫星遥感数据采用最大值合成法生成；气象气候数据来自中国气象背景数据集（空间分辨率为500 m），其利用反向距离加权平均的方法对全国1915个站点的气象数据插值生成；地貌数据来自于中国1∶100万地貌类型空间分布数据集（空间分辨率为1 km），其来源于《中华人民共和国地貌图集（1∶100万）》；土壤数据来自于中国土壤类型空间分布数据集（空间分辨率为1 km），其来源于《1∶100万中华人民共和国土壤图》；植被类型数据来自于中国100万植被类型空间分布数据集（空间分辨率为1 km），其来源于《1∶1 000 000中国植被图集》；数字高程模型（DEM）来自于SRTM V4.1数据（空间分辨率为90 m）。以上栅格数据均通过中国科学院资源环境科学数据中心（http：//www.resdc.cn/）获取，并利用ArcGIS软件完成裁剪、投影变换等处理。

1.3 研究方法

1.3.1 趋势分析

Sen氏斜率趋势分析能够较好的解决数据误差问题，且分析结果较为科学可信［18］。其计算公式如下：

式中，NDVIi、NDVIj分别为 i、j时间序列的 ND⁃VI值。当斜率S大于0时，表示NDVI呈现增长趋势，反之则呈现下降趋势。Manna-Kendall检验（MK检验）可判断趋势的显著性，且能够排除少数异常值对数据的干扰［19］。因此，本文将二者相结合，用来判断山西省NDVI在像元尺度的变化趋势及其显著性。

Hurst指数多用于描述自然界中长时间序列的自然现象［20］。对于NDVI的时间序列｛NDVI（t）｝，t=1，2，…，n，基本原理和表达式如下：

对比值 R(T)/S(T)≅R/S，如果有 R/S∝TH，则可以证明序列中体现出Hurst现象。可根据log(R/S)n=a+H×log(n)利用最小二乘法拟合得到 Hurst指数H 值［21］。当 0

1.3.2 空间转移矩阵

空间转移矩阵能够定量表征某土地利用或植被斑块在某一时间段内从T时刻向T+1时刻的转化过程，可以直观反映不同类型或不同等级斑块的面积变化和转入转出情况，已广泛应用于土地利用和植被覆盖度变化等方面［22］。其表达式为：

式中，S为研究地区面积，i、j分别为2000年与2015年的NDVI等级序号，n为NDVI分级数。本研究中将NDVI按照≤0.2、0.2～0.4、0.4～0.6、0.6～0.8、≥0.8分为5个等级。

1.3.3 地理探测器

地理探测器模型能够用来探测地理要素空间分异性和揭示其背后驱动因子，它既可以度量自变量对因变量的解释度，又可以分析两因子的交互作用对因变量的影响［15］。本文选取了年平均气温、年平均降水量、≥0 ℃积温、≥10 ℃积温、干燥指数、湿润指数、高程、地貌、土壤、植被类型等10个自然因子指标（表1）。本文利用自然断点法将年平均气温、年平均降水量、≥0 ℃积温和≥10 ℃积温划分8级；将干燥度指数、湿润指数和高程划分6级；将地貌类型重分类划分为6类：平原、台地、丘陵、小起伏山、中起伏山、大起伏山；土壤类型为10类：淋溶土、半淋溶土、钙层土、初育土、半水成土、水成土、盐碱土、人为土、湖泊河流及养殖场；植被类型9类：针叶林、针阔混交林、阔叶林、灌丛、草原、草丛、草甸、沼泽、栽培植被。在ArcGIS中利用Create Fishnet工具，生成研究区范围7 km格网，取每个网格的中心点为采样点，共3023个，然后再利用ArcGIS的Extract Values to Points工具将NDVI及各影响因子数据提取到对应采样点上，最后导入地理探测器模型中计算。

本研究采用了3个探测器工具：

（1） 因子探测

主要用于探测自然因子能够多大程度地反映植被NDVI的变化，也就是解释力度的大小，用q来衡量，公式如下：

式中，L为因变量NDVI或自变量影响因子X的分类或分区，Nh和σh2分别为层h的单元数和方差，N和σ2分别为区域内整体的单元数和方差。q值度量影响因子对NDVI的解释力，其值域为［0，1］，q值越大，表示该自然因子对 NDVI的解释力就越强。

（2） 交互因子探测

主要是识别不同自变量Xi两两交互作用对因变量NDVI的空间分异解释力是增强还是减弱。确定因子间交互影响类型，是通过对比q（X1）和 q（X2）与交互时的 q（X1∩X2）这三个量来实现的。如果 q（X1∩X2）小于 q（X1）和 q（X2）的最小值则为非线性减弱，大于两者最大值则为双因子增强，等于两者之和则为独立作用，大于两者之和则为非线性增强，如果介于q（X1）和q（X2）之间则为单因子非线性增强。

（3） 风险区探测

用于评价不同子区域属性均值的差异性，可鉴别各区域植被覆盖情况，用t统计量来检验［15］。具体公式如下所示：

2 结果分析

2.1 NDVI变化趋势分析

2.1.1 区域尺度

本研究采用ArcGIS的Cell Statistics工具，对2000−2015年山西省区域尺度的NDVI月、季、年平均值进行统计（图1）。结果表明，ND⁃VI月变化曲线呈单峰特征，2月的NDVI值最低，随后逐渐升高，并在8月达到峰值，之后迅速降低。NDVI季节平均值统计结果同样表明，夏季NDVI最高，多年平均植被指数高达0.64，其次是秋季 0.59，春季 0.43，冬季 NDVI最低只有0.22。这种特征主要受气候的影响，山西省处于温带大陆性季风气候，夏季高温多雨，冬季寒冷干燥，故6−9月植被指数值较高，且山西省处于温带落叶阔叶林带，冬春季大部分树木无树叶覆盖，导致冬春季节植被指数低。

从NDVI季节变化特征上看，春季植被指数增长的趋势最突出，平均变化率0.006/a；夏季，秋季次之，平均变化率分别为0.004/a、0.003/a；冬季植被指数增长幅度最小，变化率为0.001/a。从NDVI年际变化特征上看，2000−2015年山西省年平均NDVI为0.65，总体上处于波动上升的趋势，变化范围在0.546 8（2001年）～ 0.730 7（2013年）之间，年NDVI平均增长率为0.006/a。整体上呈现出四个周期性特征，四个时段分别为2000−2004，2004−2008，2008−2013和 2013−2015，2000−2001有 所 下降，2001−2004为增加阶段；2004−2005有所下降，2005−2008为增加阶段；2008−2009有所下降，2009−2013为增加阶段；2013−2015有所下降。

2.1.2 像元尺度

根据Sen趋势分析结果及MK检验结果（0.05 置信水平），将NDVI变化趋势与显著性分为5级（图2）。结果表明，15年间山西省NDVI变化趋势整体向好，以缓慢增长为主。极显著缓慢增长区域面积最大，占山西省面积的51.92%；极显著快速增长区域其次，占山西省面积的36.08%，主要分布在沿黄河地区、吕梁山以及系舟山和太岳山区；不显著变化区域占10.08%；极显著缓慢退化和极显著快速退化占比最小，共占1.83%，多散布在盆地地区，如大同盆地、太原盆地、临汾盆地等。极显著快速退化区大多位于各地市的核心城区，以太原市小店区最为显著；极显著缓慢退化区位于核心城区外围，是城市化主要推进区，植被覆盖度缓慢下降，以晋城市区周边最为明显；不显著变化区主要位于城郊区域，连接着盆地与山脉丘陵，植被覆盖以农作物为主，其受季节影响大，如清徐县、太谷区。

2.2 NDVI变化可持续性分析

山西省Hurst指数分析结果如图3所示，H值在0～0.3的区域零星分布在全省各处，仅占山西省面积的1.18%；0.4～0.5占比最大，占全省的面积的39.07%，特别是山西省西部呈现出明显的反持续特征；H值0.5～0.6占33.66%，0.6～0.7占10.04%，大于0.7的区域占1.67%，表现为可持续性特征的区域呈破碎状分布于吕梁山区、太岳山区和运城盆地。反持续性特征区域共占全省面积的53.91%，这表明山西省ND⁃VI受自然或人为因素影响较大，特别是西部植被较好的区域尤为明显，需要进一步加强生态保护工作，增强可持续性。

本文将Hurst指数分析图与趋势分析图在ArcGIS中叠加，得到了山西省NDVI未来趋势变化图（图3）。结果表明，持续快速退化和持续缓慢退化面积极小，仅占1.14%，大部分区域位于城市中心或城镇周边，受城市建设和城市扩张的影响较为明显；持续基本不变的区域占5.33%；波动变化的区域占全省面积的54.54%，是占比最大的区域。持续缓慢增长区域占24.18%；持续快速增长占14.71%，主要分布于吕梁山中南部地区、太岳山区、太原市东部山区（系舟山脉），这些山区受地形地貌影响较大，人类活动较少，同时受退耕还林、封山育林、植树造林等政策的影响，植被快速增长，持续性强，生态功能逐步恢复。

2.3 NDVI空间变化分析

2000−2015年山西省植被状况总体呈现出东南部植被覆盖较高，西北部植被覆盖较低的分布，区域内NDVI的空间分布差异比较显著，中部的盆地地区和吕梁山以西的黄土丘陵区植被覆盖相对稀疏（图4）。2000−2015年NDVI变化以西部、南部的山地和中部盆地为主。其中植被NDVI增加的区域主要分布在吕梁山和太行山为代表的山地地区；植被NDVI值减少的区域主要分布于太原盆地、忻州盆地、大同盆地和长治盆地为主的人类聚集区。

在ArcGIS中通过叠加分析得到了2000−2015年植被NDVI变化转移矩阵表（表3），结果表明，16年来NDVI≥0.8的区域面积大幅增加，从占比0.35%增加到25.78%，主要是由低一级植被覆盖转化而来。而0.4

2.4 NDVI空间分异及变化的地理探测

2.4.1 NDVI空间分异的驱动分析

利用地理探测器模型研究不同因子对NDVI空间分异特征的驱动作用。因子探测结果显示（表3），不同驱动因子对NDVI空间分异特征的解释力度由高到低排序为：地貌类型>年平均降水量>湿润指数>土壤类型>植被类型>年平均气温>高程>≥0 ℃积温>干燥指数>≥10 ℃积温，这说明地貌因素主导着山西省植被的空间分布状况，其次是年平均降水量和湿润指数。交互探测结果表明（表4），在两两交互的组合当中，交互作用最强的因子组合为年平均降水量和≥0 ℃积温，他们双因子交互作用的q值为0.48，体现出较强的非线性增强交互作用。

通过风险探测可以得出（图5），总体上不同分区内的NDVI均值存在差异，但并不显著。其中，NDVI随着高程的增加而增加，随着年平均降水量的增加也呈逐渐上升的趋势。

2.4.2 NDVI变化趋势的驱动分析

利用地理探测器对2000−2015年NDVI变化趋势归因，从因子探测器结果（表3）可以看出，各驱动因子对山西NDVI变化趋势的影响程度排序为：地貌类型>湿润指数>植被类型>年平均降水量>土壤类型>≥10 ℃积温>高程>≥0 ℃积温>干燥指数>年平均气温。从表5可以看出，在NDVI变化趋势影响因素的交互作用中，大部分因子两两组合后表现出非线性增强作用，也有因子组合后表现出双因子增强作用。其中地貌类型和湿润指数组合的交互作用解释力最强，q值为0.154 3；其次是地貌类型和年平均降水量（q值为0.151 1）以及年平均降水量和高程组合，q值为0.150 6；此外，年平均气温和其他影响因子的交互作用均表现出非线性增强作用。

如图6所示，NDVI的变化率随着年平均降水量的增加呈上升趋势，在610.6 mm～672.1 mm之间（分区7）变化率最高；在不同高程上，NDVI变化率非线性波动，增长最快的区域高程出现在1435 m ～ 1753 m之间（分区5）；不同的地貌类型对NDVI变化趋势也存在着空间差异，中起伏山地（分区5）的NDVI变化最快，而平原（分区1）的NDVI增长缓慢；在土壤类型方面，人为土（分区8）上生长的植被出现了减少的趋势，高山土（分区9）上生长的植被增长最快；植被类型分区中，针叶林（分区2）中的NDVI增长较快。

3 讨论

山西省2000−2015年间植被覆盖状况整体呈改善趋势，空间分布存在差异，其中山地地区植被覆盖度大于盆地地区。NDVI≥0.8的区域面积显著增加，说明退耕还林还草政策的实施发挥了重大成效，这与张珺等［23］的研究结果基本一致。

多数已有研究表明降水是植被活动的主要驱动力［24-26］。在温带季风气候显著、区域内地形复杂的山西省，地貌类型对NDVI变化的解释度最高，这可能是因为地貌类型不仅影响水热条件，并在一定程度上影响人类的活动范围，从而影响植被覆盖程度的空间分异以及变化趋势，这与同英杰等人在陕西省的研究结果类似［27］，中起伏山地NDVI变化最快，大起伏山地次之，平原地区增长最慢。地形对植被覆盖度影响较大，山脉、丘陵区植被覆盖度高且变化快，如太岳山、中条山等海拔较高的山区；盆地、河谷地区植被覆盖度低且变化慢，如大同盆地、太原盆地和运城盆地。高程在1435 m～1753 m间NDVI增速最快，表明该高程区间有利于植被生长，这也对应了山区、丘陵区的地形；另外，湿润指数对NDVI空间分异和变化趋势的影响程度都比较大，这可能是因为植被的生长不仅跟降水有关，还与蒸发和径流有关。本次研究采用交互探测器对双因子交互作用进行探讨，结果表明年平均降水量和≥0 ℃积温双因子作用下对NDVI空间分异的解释力近50%；地貌类型和湿润指数双因子作用下对NDVI变化趋势的解释力达15%。

从人类活动情况来看，山西省在第一轮退耕还林中（2000−2015）累计完成了156.70万hm²的退耕还林任务，主要包括在山西北部的风沙区以及西部的黄土丘陵区实施退耕还林，在汾河两岸和吕梁山脉一带实施荒山造林，在太行山脉一带实施封山育林［28］。本研究结果也表明NDVI极显著快速增长区域主要分布在沿黄河地区、吕梁山以及系舟山、太岳山、中条山区。另外太原市实施了环城绿化带的建设，新建多处城郊森林公园，使得太原市东西山植被覆盖度迅速提高。而太原市、大同市、朔州市、长治市、晋城市等城市的城郊区域植被指数有所下降，主要是由于城市扩张占用土地，且聚落较多、人员密集，受人类活动影响较大，植被减少明显。山西城郊地区的农作物以玉米、小麦为主，NDVI较低且受季节明显。这说明山西省植被覆盖度的增加和减少都与人类活动密不可分，由于积极的生态政策使得山区丘陵等受人类活动影响小的区域植被迅速恢复，但城市化和工业化的快速发展导致了城郊地区的植被减少。

以往对山西省NDVI的研究大多仅以气候变化为响应，本研究采用地理探测器更为全面深入地分析了5类共10种自然因子对山西省NDVI变化的驱动作用。近年来人类对自然的干预有所增加，下一步NDVI变化驱动机制的深入研究还应该结合人类活动因素如人口密度、GDP、土地利用类型、粮食播种面积等进行综合分析［29］，同时进一步厘定退耕还林等重大生态工程对NDVI变化的贡献率［30］。

4 结论

本文基于SPOT NDVI 遥感数据，采用趋势分析、空间转移矩阵分析了2000−2015年山西省NDVI 时空变化特征，NDVI年平均增长率为0.006/a，变化趋势以极显著缓慢增长为主，极显著快速增长其次，整体趋势向好，但未来趋势以波动变化为主，而NDVI≥0.8的区域面积大幅增加，从占比0.35%增加到25.78%，主要是由低一级植被覆盖转化而来。同时，本文采用地理探测器模型对影响植被变化的气候因素与非气候因素进行了探测分析，明确了山西省NDVI空间分异及变化的主导驱动因素，地貌类型、年平均降水量和湿润指数对NDVI空间分异的单因子解释力均超过了15%，而地貌类型对NDVI变化趋势的单因子解释力也是最高的。因子交互作用均能增强其对植被空间分异和变化趋势的解释力度，对NDVI空间分异的影响最显著的是年平均降水量和≥0 ℃积温双因子组合，而地貌类型和湿润指数的非线性增强交互作用对NDVI的变化趋势影响最为显著。未来山西省NDVI变化驱动机制的深入研究还应充分考虑人类活动因素，并厘定退耕还林等重大生态工程对NDVI变化的贡献率。