2000-2015年太行山地区NDVI变化及驱动因素分析
2020-08-19李新元王莉景海涛范晨雨余汛
李新元 王莉 景海涛 范晨雨 余汛
摘要：基于2000-2015年MODIS NDVI数据，结合同期气象、土地利用等相关数据进行研究。根据相关计算获得NDVI变化数据和土地利用类型变化数据，对NDVI变化数据进行反向研究，寻找影响NDVI变化的驱动因素。结果表明，2000-2015年研究区植被覆盖情况整体得到改善，平均增速為0.0075/年;NDVI增加区域的驱动因素是气象因素对植被生长的促进作用;人类活动使耕地的农作物长势更加旺盛，直接导致了耕地范围内的NDVI提高。太行山绿化工程、退耕还林工程、农村人口向城市转移等生态保护措施对草地、林地、水域的植被生长也产生了积极影响，人类活动对离城市较远的自然区域植被起促进作用;NDVI减少区域的驱动因素是人类活动使太行山城市周边地区的耕地、草地面积减少，城市面积增加，对城市及其周边地区的植被生长起减弱作用。
关键词：NDVI;ArcGIS;植被覆盖度;气候因素;人类活动;太行山地区
中图分类号：Q948.1 文献标识码：A
文章编号：0439-8114（2020）09-0070-07
DOI：10.14088/j.cnki.issn0439-8114.2020.09.014
植被作为陆地生态系统的主体，它与地质、地貌、气候、水文、土壤、动物和微生物共同构成了自然地理环境，是最能反映其他要素性质的指示者[1]。同时，植被还是连结土壤、水分与环境的自然纽带，在陆地碳平衡和气候系统调节中发挥了重要作用。变化环境下的植被生态系统演变规律、驱动机制和预测调控研究是全球变化研究中的核心内容之一[2-4]。监测植被动态变化，评估人类活动和气候变化对植被变化的影响程度，对确定合理的生态工程布局和适应性管理对策具有重要的实用价值，已成为变化环境下植被生态系统可持续管理研究领域的前沿和热点[5-7]。
归一化植被指数（NDVI）是反映植被生长状态的最佳指示因子，能够表征植被覆盖的整体状况，且与植被覆盖度之间存在着极为显著的线性关系，因此利用NDVI计算植被覆盖度的方法被广泛使用。
太行山脉作为京津和华北大平原的天然屏障，既可抵御西北寒潮的袭击，又可接纳东南暖湿气流，是海河水系的发源地及黄河中下游部分支流的上游。太行山区山高坡陡，岩石裸露，植被稀少，由于森林资源曾遭到严重破坏，水土流失严重，易暴发山洪。政府相关部门为改善太行山地区的植被覆盖度，实施了太行山绿化工程等一系列植被恢复计划。目前有关太行山地区的植被变化已经有一部分研究成果，例如李薇等[8]关于太行山区不同坡度NDVI变化趋势差异的分析研究，李晓荣等[9]对太行山区植被NPP时空变化特征及其驱动力的分析，胡实等[10]对中国典型山地植被垂直地带性特征及其影响要素的研究，以及其他有关太行山地区的零散的、具有小区域特征的研究。通过研究分析发现已有的研究结果没有深入分析太行山地区的植被覆盖变化信息和植被变化的驱动因素。因此，本研究利用MODIS NDVI数据，结合气温、降水、DEM、居民地分布、土地利用类型等数据，对太行山地区2000-2015年16年间植被覆盖度的时空变化情况进行分析，并对影响植被变化的气候因素和人类活动进行相关分析。
1 研究区概况、数据来源与方法
1.1 研究区概况
太行山地区位于北纬34°34′-40°43′、东经110°14′-114°33′，是中国东部地区的重要山脉。太行山脉从东北向西南延伸，北起拒马河谷，南至晋豫边境黄河沿岸，西临山西省汾河，东濒华北平原，横跨山西、河南、河北、北京4省市。太行山地理环境复杂，地势西缓东陡，受河流切割，多横谷，为东西交通要道，气候属于暖温带半湿润大陆性季风气候。
本研究中研究区的划分与传统意义上以行政区边界的划分方法不同，研究区的选择更多考虑自然环境因素。研究区面积为12.7484万km2，以山区为主，不包括太行山脉周边的平原地区，在行政区域上无明显界限（图1）。
1.2 数据来源与预处理
采用的NDVI数据来源于中国科学院计算机网络信息中心国际科学数据镜像网站（http：//www.gs-cloud.cn）MODNDIM中国500m NDVI月合成产品，数据是由MODNDID计算得到，计算方法为取月内每天最大值。坐标系EPSG：4326（WGS84），空间分辨率500m，时间分辨率每月。时间跨度为2000-2015年。
根据研究区矢量边界数据对NDVI数据进行裁剪。本研究所使用的MODIS NDVI数据已经经过水、云、气溶胶等处理，减弱了部分因素对数据的干扰。为进一步消除其他异常值对研究的影响，采用最大值合成法（Maximum value composites，MVC），将月最大NDVI合成年最大NDVI，以年最大NDVI来表示该地区的植被覆盖情况，在一定程度上提高了数据的可靠性[11]。由此获得研究区2000-2015年16幅影像图。
气象数据来源于中国气象数据共享服务网，包括太行山区域内18个气象站2000-2015年日降水量数据。对气象数据进行剔除异常值等预处理，利用反距离权重插值法将每年的气温、降水数据插值为与MODIS NDVI数据集相同空间分辨率的栅格影像。研究区边界矢量数据来源于中国科学院山地表生过程与生态调控重点实验室[12]。太行山DEM数据下载于地理空间数据云平台。太行山绿化工程相关数据来源于中国林业网。土地利用类型数据来源于ESACCI-LC的2000年和2015年全球土地利用类型覆盖图，空间分辨率为300m。根据研究区范围进行裁剪后按照联合国土地覆盖分类系统（LCCS）进行分类，分为耕地、草地、林地、城市、水域、裸地6类土地利用类型。
1.3 研究方法
1.3.1 像元二分模型 归一化植被指数是一种由遥感传感器所接收的地物光谱信息推算而得的反映地表植被状况的定量值。根据像元二分模型，一个像元的ND VI可以表达为由绿色植被部分所贡献的信息NDVIveg与由无植被覆盖（裸土）部分所贡献的信息NDVIsoil这两部分组成，像元植被覆盖度（NDVIFVC）表达式如下：
式中，NDVImax为该像元上校正过的归一化植被指数;NDVIsoil为全部是裸土或无植被覆盖区域的归一化植被指数;NDVIveg为全部被植被所覆盖的归一化植被指数[13]。NDVIveg和NDVIsoil取給定置信度的置信区间内的最大值与最小值，分别取累积频率为99.5%与0.5%时对应的值作为NDVIveg和NDVIsoil，小于NDVIsoil的像元赋值为0，大于NDVIveg的像元赋值为1。
1.3.2 Mann-Kendall检验Mann-Kendall（简称为M-K）非参数检测法所用样本不需要遵从一定的分布规则，也不受少数异常值的干扰，因此其被广泛应用于长时间序列的趋势和突变检验分析[14，15]。在M-K检验中，将时间序列植被覆盖率X作为独立分布的样本数据，Xi表示年际植被覆盖率由区域内所有像元对应年份的NDVI求平均值计算得到，S表示秩序列。M-K检验计算公式为：
趋势变化由OF检验值判断：
式中，E（S）和var（S）分别为S的均值和方差。
若OF>0，表明序列呈上升趋势UF<0，则表明序列呈下降趋势。将植被覆盖度的逆序列按上述步骤进行计算，使|UB|=-UF，给定显著性水平α=0.05，当|UB|曲线超出临界线时变化趋势显著，而OF和UB在临界线之间出现交点，则交点对应的时刻或区域就是突变的开始。
1.3.3 趋势分析 一元线性回归分析可以模拟每个栅格在研究时间内的变化趋势，对每个像元进行基于空间位置的时间变化特征能够反映整个研究对象的时空变化特征，综合反映出区域时空格局的演变过程[16，17]。Slope为像元NDVI线性回归方程的斜率（s），表示随时间变化的植被指数变化趋势[18]。其计算公式如下：
式中，s为像元NDVI回归方程的斜率;n为年跨度，i为年序号，NDVIi为第i年的NDVI。s>0表明NDVI的变化趋势是增加的，反之则是减少的。采用F检验进行显著性分析，公式为：实观测值;yi为未能拟合部分的yi;y为覆盖度的16年平均值，并将结果划分为极显著相关（P<0.01）、显著相关（0.01
1.3.4 相关分析 利用研究区处理过的年NDVI数据分别与插值后的气象数据进行相关分析，计算相关系数，分析植被覆盖变化与气象因子变化之间的相关性[9]。
式中，i为年序号;NDVIi为第i年的NDVI;NDVI是对应像元16年的年平均NDVI;n为统计年数（16年）;xi为气象要素;x为气象要素多年平均值;r为Pearson相关系数。r>0表示气象要素对植被覆盖度的影响呈正相关;r<0表示气象要素对植被覆盖度的影响呈负相关。采用T检验进行显著性分析，并将结果划分为极显著相关（P<0.01）、显著相关（0.01
1.3.5 残差分析 利用复直线相关回归模型拟合出受气候因子影响的植被覆盖度，进而计算其与真实植被覆盖度之间的差值，忽略其他不确定因素，把差值作为人为因素来考虑，该方法即为残差分析[19]。
对NDVI与温度、降水之间进行复直线回归分析，回归方程为：
Z=a+bx+cy（8）
其中的参数公式为：
式中，x、y、z分别为气温数据x、降水数据y和NDVI数据z的均值;σx、σy、σz分别为各自数据系列的均方差;rxz、rxy、ryz分别为x与z，x与y、y与z的相关系数。
NDVI模拟残差计算公式：
NDVIε=NDVI-NDVIP
对图3进行统计分析，2000-2015年太行山地区Slope值大于0的区域面积为104728.1km2，占总面积的82.15%，植被变化表现为增加;Slope值小于0的区域面积为22755.9km2，占总面积的17.85%），植被变化表现为减少。采用F检验进行显著性分析，结果为植被变化极显著（P<0.01）区域占总面积的18.00%，植被变化显著（0.01
2.2 太行山地区2015年气象因素对NDVI贡献值分析
利用2000-2015年NDVI的时间序列数据以及气温、降水数据求出參数a、b、c的值，通过确定NDVI与气温、降水之间的复直线相关回归方程公式带入相应的气温、降水数据得到2015年的NDVI预测值（图4）。这个预测值就是气象因素对NDVI的贡献值。
将2015年气温和降水对植被NDVI的贡献值空间分布图（图4）结合太行山地区的DEM图像（图1）进行分析发现，高海拔地区气候因素对NDVI的贡献值比低海拔地区更大，人类居住地及其周围的区域气候因子对植被NDVI的贡献值较小。NDVI贡献值在0.5以上的区域较多分布在地形较为复杂的山地，这源于地形复杂的区域人类活动较少，有比较完整的生态系统，所以在气候条件适宜的情况下植被自然生长更加迅速。北部地区由于长期的矿物开采导致生态环境破坏严重，植被稀少，植被恢复比较困难，所以该地区自然状态下的植被生长比其他地区恢复的更慢。
统计分析气候因子对植被NDVI的贡献值，结果（表2）显示，2015年99%以上的区域气候因子对NDVI的贡献为正，气候因子对NDVI贡献值在0.5～1.0的区域占总面积的85.22%;而2015年NDVI实际值在0.5～1.0的区域占总面积的86.84%，二者在相同区间（0.5～1.0）面积相近。这说明气象因素是太行山地区植被覆盖变化的主要驱动因素之一，表现出对植被生长的积极促进作用。
2.3 太行山地区2015年人类活动对NDVI的影响分析
为定量化分析人类活动对太行山地区植被覆盖的影响，将实际观测值减去利用复相关回归分析计算得到的预测NDVI值，即得到残差值，认为该残差值就是人类活动对NDVI贡献值（图5）。
统计分析人类活动对NDVI贡献值，结果（表3）显示，对植被NDVI起降低作用的区域面积为35237.68km2，占太行山地区总面积的27.64%：对太行山地区植被NDVI起促进作用的区域面积为92246.32km2，占太行山地区总面积的72.36%。
根据从欧空局获取的预处理后获得的2000年和2015年的不同土地利用类型数据进行叠加分析得到不同土地利用类型之间的流转变化情况，进行统计分析得到2000-2015年土地利用类型变化的转移矩阵（表4）。结果显示，在6种土地利用类型中耕地的总面积减少了844.33km2，变化率为-2.31%;草地的总面积减少了1877.47km2，变化率为-3.76%;林地的总面积增加了206.99km2，变化率为0.57%;城市的总面积增加了2437.06km2，变化率为244.14%：水域的总面积增加了77.74km2，变化率为18.86%;裸地的总面积无变化。由表4还可以看出，耕地、草地、林地、水域都有转入和转出变化，但城市无转出变化。耕地、草地面积减少显著，林地面积略微增加，城市面积增加显著，水域、裸地变化不明显。耕地、草地、林地、水域都有部分区域转入为城市用地，其中草地和耕地转入面积最多。
对2000年和2015年各土地利用类型的矢量区域进行提取，结合研究区在研究时段内的NDVI和变化趋势Slope值进行统计分析获得不同土地类型的平均NDVI变化结果（表5）和Slope值小于。的区域在2015年不同土地类型分布图（图6）。其中裸地面积太小，平均NDVI值变化不具有代表性，不进行统计。由表5可知，2000-2015年耕地的平均NDVI增加了0.15;草地的平均1VDV1增加了0.10;城市的平均NDVI增加了0.10;水域的平均NDVI增加了0.06。由此可知，16年来为使农作物更加高产，人们通过品种改良、科学施肥、科学管理等手段使耕地种植的农作物长势更加旺盛，直接导致了耕地范围内的NDVI最大值提高，此结论与李晓荣等[9]对太行山地区植被NPP的研究结果相同。
由图6可知Slope值小于。的区域在南部都是围绕城市相对集中分布，在北部地区的分布比较零散。将城市的分布信息与各种植被变化信息结合分析可知，2000-2015年NDVI降低的区域大部分集中在城市及周围区域，除北部自然环境较恶劣的区域外，气象因素对NDVI贡献值较低的区域都是城市及其周围地区，人类活动的NDVI贡献值为负的区域也大都分布在城市及其周围地区。
3 小结与讨论
3.1 小结
本研究基于2000-2015年MODIS NDVI数据，结合气温、降水、DEM、行政区域分布和土地类型数据，对太行山地区2000-2015年16年间的植被覆盖度的时空变化特征进行研究，对影响其变化的气候因素和人类活动对植被覆盖度深层变化情况进行相关分析，再结合分析结果对太行山地区的植被变化情况及其原因做初步总结，结果如下。
