太行山区NDVI时空变化及其与气候因子的关系

刘 琪，景海涛，刘盼盼，汤金颖，王 磊

(河南理工大学 测绘与国土信息工程学院，河南 焦作 454000)

植被是陆地生态系统的核心组成部分，是重要的自然资源，起到连接生态系统中大气、土壤、水分等自然因子的纽带作用[1,2]。植被的改善有助于提高水源涵养和水土保持功能，对植被的时空变化趋势监测在全球气候变化背景下的陆地生态系统研究中具有重要意义[3]。植被遥感数据凭借优异的空间和时间特征，现已成为监测全球及区域尺度植被的生长状态、地表覆盖变化等植被信息的关键数据源。当前用于度量植被变化的植被指数有很多，其中归一化植被指数(NDVI)既能够表征植被覆盖的整体状况，且与植被覆盖度之间存在显著的线性关系，是最常用的表征植被生长状态和覆盖度的指标，被广泛用于植被变化及其驱动因子的研究当中[4]。

目前，国内外学者针对植被覆盖度的变化及其影响因子已经做了大量的研究[ 5-7]。研究结果表明，NDVI与气候因子有着显著的相关关系[8,9]，气候条件是影响植物生长的关键性条件[10]。伴随着经济快速发展，各地生态系统频繁遭到损伤，全球气候变暖趋势日益严重[11]。基于此，理解和揭示植被与气候因子之间的内在联系，对于区域生态保护治理具有重要的意义。太行山地处我国第2阶梯的东沿，是平原到山地高原的过渡带。同时，又是东部经济发达地区到中西部经济欠发达地区的过渡带[12]。目前关于太行山区的植被变化研究已经取得相当一部分的成果，李新元等[4]对太行山区NDVI变化及驱动因素分析发现气象因素对NDVI变化起促进作用；
李薇等[13]对太行山区不同坡度NDVI变化趋势分析发现，太行山区不同坡度NDVI变化趋势存在显著差异；
范松克[14]对于太行山南麓进行归一化植被指数变化及成因分析发现，NDVI与气温和降水呈正相关等，已有的研究结果虽能一定程度反映植被与气象因子的相关关系，但大多为整体上分析太行山区植被覆盖变化驱动因子的研究，很少单独考虑气象因子对太行山区植被覆盖变化趋势的影响，缺乏植被对气象因子响应机制的深层分析，并且未考虑植被对气象因子的时滞效应。

本研究基于1998-2018年21 a时间序列的NDVI数据，采用一元线性回归及趋势分析等方法，分析太行山区植被的时空变化特征；
同时，结合气候数据，利用相关分析、滞后性分析和统计分析等方法，在时间尺度和空间尺度探讨植被与气候因子的响应关系，分析NDVI与气温和降水的时间滞后关系，为太行山区的植被变化监测、生态环境保护和治理提供有力的依据。

1.1 研究区概况

太行山区即太行山地及其毗邻区域[15](34°34′-41°06′N，110°13′-114°33′E)，总面积约12.7万km2，纵跨北京、河北、山西、河南4省市，山脉北起北京市西山，向南延伸至河南与山西交界地区的王屋山，西接山西高原，东临华北平原，呈东北-西南走向。太行山区是中国地形第2阶梯的东缘，也是黄土高原的东部界线，海拔39～2 994 m，最高峰是五台山(图1)。太行山区地势北高南低，地貌类型复杂，从气候上看，太行山属于温带半湿润大陆性季风气候，年均气温为10.3 ℃，年降水量为501.8 mm。全年冬无严寒，夏无酷暑，雨热同季，虽四季分明，但冬长夏短。

1.2 数据来源与处理

NDVI数据来源于中国科学院资源环境科学与数据中心(www.resdc.cn/)，NDVI数据集是基于SPOT/VEGETATION NDVI卫星遥感数据，在旬数据基础上采用最大值合成法生成的月度植被指数数据集，时间跨度为1998-2018年，空间分辨率为1 km。所获取的NDVI数据经过水、云、气溶胶等处理，有效减弱了其他因素对数据的影响，提高了数据的准确性。对于获取的NDVI数据进行简单预处理，在ArcGIS中通过太行山区矢量边界进行掩膜裁剪，使用栅格计算器计算出月、年植被指数NDVI数据集。

图1 研究区地理位置

气象数据来源于中国气象网(http://www.cma.gov.cn/)，时间跨度1998-2018年。通过Python编程从全国所有的站点中选取太行山区内的18个基准气象监测站点的日均气温和日累计降水量，进一步处理获取月均气温/降水和年均气温/降水数据。根据各气象站的经纬度将气温、降水信息导入ArcGIS空间分析模块运用反距离权重法进行空间插值。通过太行山区矢量边界进行掩膜裁剪、重采样、重投影得到与NDVI数据具有相同像元大小、相同坐标系的栅格影像数据集。

DEM数据由美国奋进号航天飞机的雷达地形测绘SRTM(shuttle radar topography mission，SRTM)数据经重采样生成，空间分辨率为1 km。将DEM数据导入ArcGIS中并计算出高程，使用太行山区矢量边界进行掩膜裁剪、重采样、重投影得到与NDVI数据具有相同像元大小、相同坐标系的栅格影像数据集。

1.3 研究方法

1.3.1 NDVI变化趋势分析 采用一元线性回归方程的斜率来模拟栅格点的变化趋势[16,17]，斜率的正负表示变化趋势的增减。进而模拟分析NDVI的年际变化趋势，计算公式为：

(1)

式中:slope为1998-2018年太行山区NDVI的变化速率；
n为研究的时间序列长度，取n=21；
i为1～21的年序号；
NDVIi为第i年NDVI值。在Arc GIS中进行栅格数据的计算与统计分析得出太行山区21 a NDVI年际变化趋势。

1.3.2 相关分析性 采用皮尔逊(Pearson)相关性分析方法，对NDVI数据与气温、降水空间插值数据进行逐像元相关性分析[18-20]。由此来反映出气温与降水等气候要素与NDVI序列的时间相关程度和空间分布规律。相关系数计算公式为:

(2)

通过查阅相关系数表，结合研究区具体情况，进行显著性等级划分。当n=21时，rxy的相关性等级划分见表1。

表1 显著性等级划分

1.3.3 偏相关性分析 相较于简单相关分析，偏相关分析是探索植被生长与单一气候因子之间的联系，可剔除其他气候因子的影响，有助于更好地分析影响植被变化的气候因子主控因素，被广泛用于检测影响植被生长的主要气候驱动因子[21]。计算方法如下:

(3)

式中:x为NDVI；
y为气温(℃)；
z为降水(mm)；
rxy,z为将降水固定后气温与NDVI的偏相关系数，降水与NDVI的相关系数同理。

1.3.4 滞后性分析 利用滞后相关系数进行太行山区NDVI对气候因素的滞后性分析[17,18]。将月均NDVI与同期、前1个月以及前2个月的月均气温/降水数据进行相关性分析，根据相关系数和显著性检验结果进行统计分析，相关系数越大，显著性越强则代表2个因素间相关性越强。通过比较相关性的大小来确定时间滞后关系。计算公式为：

Rmax=max{|R0,R1,R2,R3……|}

(4)

式中:Rmax表示最大的相关系数；
R0、R1、R2、R3……分别表示与同期、前1个月、前2个月、前3个月……的相关系数，以此类推。若Rmax=RA，即代表滞后A个月。

1.4 研究流程

在遥感、GIS技术的理论支撑下，采用NDVI卫星遥感数据和气象插值数据，运用统计学思想并结合Python语言进行数据处理分析。在时间、空间2个尺度展开对太行山区植被的变化特征及其对气象数据变化的响应关系进行研究分析，以明确太行山区21 a来植被变化与气候因子的响应关系，为太行山区生态保护治理提供科学支撑。技术路线见图2。

图2 本研究技术路线

2.1 NDVI的时间变化趋势

对逐年NDVI数据集(1998-2018年)进行统计分析年际NDVI变化总体趋势(图3)。太行山区的NDVI值介于0.62～0.78，2001年植被NDVI为最小值0.62，2018年植被NDVI出现最大值0.78，21 a间太行山区NDVI涨幅约25%。结合图像和拟合的线性函数可知，太行山区的NDVI呈波浪式的增长趋势，平均增长速率为0.03/10 a(线性增长率为0.067/10 a)，其中R2=0.849，通过0.01显著性，说明NDVI的总体变化趋势与时间相关性显著。

图3 NDVI年际趋势

利用Mann-Kendall方法对NDVI变化趋势进行显著性检验和突变点检测[22](图4)。由M-K检验可知年均NDVI的Z值为4.74，通过0.01的显著性检验，说明NDVI的年际增加趋势明显，与R2的结果呼应。由突变点检测可知，1998-2018年太行山NDVI的突变点是2010年，且超出临界值1.96，说明突变并不显著。结合温度与降水的数据分析可知，2010年NDVI的突变可能是由于气温降低和降水异常增加等气候因素所导致的。综合分析可知，在时间尺度上太行山区的整体植被生长状态良好(0.67～0.78)，且增长趋势显著性明显，结合NDVI的时间变化趋势方程(y=0.006 7x+0.635 5)，预测太行山区植被生长状态未来会进一步提升。

图4 M-K趋势检验

2.2 NDVI的空间变化趋势

对1998-2018年太行山区年NDVI采用均值法得到太行山区21 a的NDVI 的空间分布图(图5A)。由图5可知，研究区的植被覆盖状况在空间分布上有明显的区域性差异，呈现出西北低、东南高、中部区域高低值交叉分布的特点；
太行山区NDVI的空间分布值介于0.16～0.87。通过阅读文献[4,6,20]和查找资料，结合太行山区的实际情况，将太行山区的NDVI进行等级划分，划分结果为 <0.2、0.2～0.4、0.4～0.6、0.6～0.8、>0.8这5个等级(图5B)。对各个等级的NDVI进行逐像元统计分析，各部分像元所占的比例为3∶83∶3 200∶24 666∶3 865。其中太行山区89.67%区域的NDVI>0.6，77.5%区域的NDVI介于0.6～0.8，整体看，太行山区植被覆盖度良好。为更好地描述太行山区NDVI在空间尺度的分布，结合县级区划图在ArcGIS中进行讨论分析，可知太行山区西北沿线区域(怀来、延庆、古县以及浮山县等)的NDVI为0.4～0.6；
在太行山东南区域(夏县、阳城、泽州等)的NDVI>0.8；
而小于<0.2及0.2～0.4的NDVI区域主要呈零星状分散在太行山区。这也正好符合太行山区NDVI西北低、东南高的区域差异性分布规律。

A.空间分布；
B.等级划分。

对太行山区1998-2018年NDVI逐像元做一元线性回归分析，根据逐像元NDVI的变化率得出NDVI变化趋势的空间分布格局(图6A)。根据回归显著性检验结果(P<0.01)，将太行山区植被NDVI变化情况划分为显著减少(slope<0)、无显著变化(slope=0)和显著增加(slope>0)3类(图6)。在ArcGIS中统计分析可知，显著增加(slope>0)的像元个数占总像元个数的93.46%，占太行山区的大部分面积；
显著减少(slope<0)的像元个数占总像元个数的6.52%，且主要集中在太行山区的东部边缘沿线(昌平—房山—易县—顺平—曲阳—平山—鹿泉，武安—邯郸—林州—鹤壁—辉县—焦作—济源等)、西部边缘沿线(永济—运城—侯马—临汾,延庆—怀来—广灵—浑源—应县—代县)以及中部的阳泉市、昔阳市、中南部的长治市、高平市和晋城市等地，其中东部沿线显著减少(slope<0)的像元个数和聚集程度明显大于西部沿线。

为更好探究太行山区NDVI的空间变化趋势，结合DEM数据进行叠加统计分析可知，太行山区NDVI变化呈显著减少区域(slope<0)主要集中在海拔高度500 m左右的区域，这些区域地势平坦，随着经济发展和城市的扩张，大量的耕地和林地转为建设用地，导致这些区域的植被覆盖度减少。而太行山区NDVI变化呈显著增加区域(slope>0)主要集中在海拔1 500 m左右的区域，这些区域高程较高，地势较陡，人类活动相对较弱，所以植被可以自然生长，植被覆盖度逐年增高。太行山区NDVI变化趋势的分析进一步验证了前面Mann-Kendall方法的结果，太行山区植被覆盖逐年变好并且显著性明显。

A.年际变化趋势；
B.变化趋势分类。

2.3 NDVI与气温/降水的相关性分析

2.3.1 NDVI对气温的响应关系 为研究NDVI与气温在时间尺度上的相关性，对太行山区1998-2018年的年均NDVI与年均气温进行统计分析，绘制年均NDVI和年均气温变化趋势图，并添加年均气温线性趋势线和趋势方程(图7)。基于变化趋势图分析太行山区NDVI与气温的年际关系，并在SPSS中进行显著性检验，结果见表2。

结合图7和表2分析可知，太行山区年际平均气温在10.0～11.5 ℃，且太行山区年均气温变化趋势平缓。年均气温在2017年最大(11.46 ℃)，2012年最小(10.08 ℃)，年均气温的变化趋势为0.039/10 a。在显著性检验中，相关系数R2=0.026，显著性P=0.91，未通过0.05的显著性检验。年均NDVI与年均气温在年际时间尺度上不存在显著的相关性。

太行山区不同地方不同的自然条件、植被的生长状态和种类不同，对于气温和降水的响应程度各不相同，故采用偏相关性分析来研究太行山区NDVI对气温和降水在空间尺度的响应关系，并根据表1对偏相关性进行显著性等级划分。

图7 年均气温变化趋势

太行山区年均NDVI与年均气温的偏相关系数和显著性等级见图8。太行山区植被NDVI与气温在空间上呈现负相关趋势，但趋势并不显著。大体上可以概括为，太行山区中部地区NDVI与气温呈负相关，西南、西北地区呈正相关。NDVI与气温在空间上呈负相关的面积占总体面积的54.37%，主要是太行山区的中部地区即山西省的东南部分区域(武乡县、襄垣县、沁县、长治市等)和南部边缘区域(辉县、焦作、博爱、济源等)以及东北区域沿线(延庆、昌平、房山、易县等)。其中通过0.05显著性检验的区域仅有0.96%，未通过显著性检验的面积占53.41%。NDVI与气温在空间上呈正相关的面积占总体面积的45.63%，主要是太行山区的西南区域(永济、运城、芮城、平陆、夏县等)和西北区域(涿鹿、蔚县、广灵、灵丘、繁峙等)，其中通过0.05显著性检验的区域有1.34%，未通过0.05显著性检验的面积是44.29%。

A.偏相关系数；
B.显著性检验。

2.3.2 NDVI对降水的响应关系 同理绘制年均NDVI和年累计降水量变化趋势图，并添加年累计降水量的线性趋势线和趋势方程(图9)，并在SPSS中将结果进行显著性检验(表2)。

结合图8和表2可知，太行山区年际年累计降水量在385～685 mm，且太行山区年累计降水变化波动明显。年累计降水量在2001年最少(385.59 mm)，2003年最多(664.45 mm)，年均累计降水量趋势为21.279 mm/10 a。在显著性检验中，相关系数R2=0.286，显著性P=0.209 ，未通过0.05的显著性检验，即年均NDVI与年累计降水量在年际时间尺度上不存在显著的相关性。但相较于年均气温与年均NDVI的相关系数和显著性检验结果，年均NDVI 与年累计降水量的相关系数较大、显著性更强，故年均NDVI与年累计降水量在年际时间尺度上相关性更强。

图9 年累计降水量变化趋势

表2 月均NDVI与不同时期气温和降水的相关性

太行山区年均NDVI与年累计降水量的偏相关系数和显著性等级见图10。由图10A可知，太行山区植被NDVI与降水相关性在空间上呈现以正相关为主的趋势。整体上可以概括为从西北到东南相关系数逐渐变小。NDVI与降水在空间上呈正相关面积占总体面积的81.89%，呈负相关面积占总体面积的18.11%。除太行山区东南区域沿线(左权—涉县—平顺县—林州市—陵川县—焦作市)和北部边缘区域(延庆县、北京市部分区域)与降水是负相关外，其余区域与降水都呈现正相关。由图10B可知，NDVI与降水呈正相关区域中通过0.01显著性检验的区域有8.42%，通过0.05显著性检验的区域有12.53%，未通过显著性检验的有区域有60.94%，呈不显著正相关。呈负相关区域中通过0.05显著性检验的区域仅有0.1%，其他区域都呈不显著负相关。

A.偏相关系数；
B.显著性检验。

2.3.3 NDVI变化驱动因子分区 为了更好地分析NDVI对气温和降水的响应关系，探究造成NDVI在不同地理空间上的差异，参考陈云浩等[23]的研究，对太行山区NDVI变化进行驱动因子分析研究，空间分布结果见图11，统计分析结果见图12。

将驱动因子分区空间分布结果(图11)与高程图(图1)进行叠加分析可知，太行山区有约57%的区域受气温和降水的共同影响，其中受气温、降水强影响的区域([T+P]+)面积约56%(T表示气温驱动型区域，P表示降水驱动型区域)，在空间上呈从西南到东北逐渐减少趋势，并主要集中在海拔0～1 000 m；
受气温、降水弱影响的区域([T+P]-)面积约1%，主要呈零星点状分布在太行山区的南部，并主要分布在海拔1 000～1 500 m；
受降水驱动型的区域(P)面积约27%，主要分布在太行山区的北部，在不同海拔上并没有明显的分布差异(500～1 500 m)；
受气温驱动型的区域(T)面积约12%，主要分布在太行山区的南部以及太行山区北部的边缘区域(延庆、昌平等)，海拔1 000～1 500 m；
只有南部的少部分地区(4%)是非气候因子驱动型(NC)，即人类活动对太行山植被的影响，海拔在500 m左右。

2.4 滞后性分析

为了更深入地研究太行山区NDVI对气温/降水的响应关系，考虑到气候因子对植被覆盖度存在滞后性，故对太行山区的NDVI与气温和降水进行年内关系分析。将月均NDVI与同期的月均气温/降水、前1个月的月均气温/降水、前2个月的月均气温/降水进行相关性分析(表2)，并对相关系数进行显著性检验。以此来验证太行山区NDVI与气温和降水间的滞后性关系。

图11 驱动因子分区

图12 驱动因子分区饼图

对结果进行联合统计分析可知，月均NDVI与同月气温不存在明显的相关性(图13A)。月均NDVI与同期月均气温的相关系数是0.355，其中显著性P>0.05的面积占总面积的99.8%；
而月均NDVI与前1个月的平均气温存在极显著的相关性(图13B)，月均NDVI与前1个月的平均气温的相关系数是0.916，其显著性达到0.01的面积占总面积的98.44%，显著性达到0.05的面积占总面积的0.54%，只有太行山区的南部边缘区域(济源、孟州)和西南部分区域(临汾)显著性P>0.05，占总面积的1.02%，未通过0.05显著性检验；
月均NDVI与前2个月的气温的相关性较显著(图13C)，月均NDVI与前2个月的平均气温的相关系数是-0.560，其中通过0.01显著性的区域面积占总面积的10.05%，通过0.05显著性的区域面积占总面积的47.76%，未通过0.05显著性检验的区域面积占总面积的42.19%。

综上可知，月均NDVI与前1个月的气温相关性最强，即太行山区月均NDVI对月均气温存在1个月的滞后性。

对结果进行联合统计分析可知，月均NDVI与同月累计降水存在明显的相关性(图14A)。月均NDVI与同月累计降水的相关系数为0.580，其中通过0.01显著性检验的面积占总体面积的32.13%，通过0.05显著性检验的面积占总面积的28.81%，未通过0.05显著性检验的面积占总面积的39.06%；
月均NDVI与前1个月累计降水不存在明显的相关性(图14B)。月均NDVI与前1个月累计降水的相关系数为0.214，其中通过0.01显著性检验的面积占总体面积的0.02%，通过0.05显著性检验的面积占总面积的0.71%，未通过0.05显著性检验的面积占总面积的99.27%；
月均NDVI与前2个月累计降水存在极明显的相关性(图12C)。月均NDVI与前2个月累计降水的相关系数为0.853，其中通过0.01显著性检验的面积占总体面积的94.89%，通过0.05显著性检验的面积占总面积的3.7%，未通过0.05显著性检验的面积占总面积的1.41%。

图13 月均NDVI与不同时期气温的相关系数和显著性检验

图14 月均NDVI与不同时期降水的相关系数和显著性检验

综上可知，月均NDVI与前2个月累计降水的相关性最强，即太行山区月均NDVI对月均累计降水存在2个月的滞后性。

1)在时间尺度上，太行山区NDVI变化与时间关系显著，呈波浪式的增长趋势，从1998年的0.67增加到2018年的0.78，平均增长速率为0.03/10 a(线性增长率为0.067/10 a)，其中R2=0.849，通过0.01显著性，这说明太行山区NDVI的增长与时间相关性显著。这一结果与李新元等[4]、李晓荣等[12]的研究结果大体上一致，但由于研究时间跨度更大，增幅也相对较大，表明随着太行山区生态修复工作以来，植被覆盖状况累年变好。

2)在空间尺度上，太行山区植被覆盖在空间分布上有明显的区域性差异，呈现出西北低、东南高、中部区域高低值交叉分布的特点。太行山区93.46%的区域面积呈现增加趋势，且集中在中高海拔地区，而低海拔地区NDVI呈现减少趋势，其主要原因是经济发展和人类活动的增加导致大量林地和耕地转向建设用地。

3)太行山区植被NDVI与气温在空间上呈现不显著负相关。其中，太行山区中部区域NDVI与气温呈负相关，西南、西北区域呈正相关；
NDVI与降水呈现显著正相关，整体上可以概括为从西北到东南相关系数逐渐变小。总体上看气温和降水等气候因子会促进植被覆盖率的增长，保护太行山区的植被覆盖必须考虑气候因子的影响。

4)对NDVI变化驱动因子分区研究发现，气温和降水共同影响的区域占太行山总面积的57%，且主要是低、中海拔区域；
气温驱动型区域占太行山总面积的12%，主要分布在太行山区的南部中、高海拔区域；
降水驱动型区域占太行山总面积的27%，主要分布在太行山区的北部，在海拔上没有显著的特征。非气候因子驱动型(NC)区域占太行山总面积的4%，主要分布在太行山区南部低海拔区域。

在对太行山区进行生态环境保护和修复时，低海拔区域要重点关注人类活动的影响，控制太行山区植被受人类活动影响的范围区域，进行合理规范的利用；
在中、高海拔区域，随着海拔的升高，降水量的增加以及气温的降低，影响植被生长的主要因素变为气温，所以要加大对气温的监控，有效预防严重自然灾害的发生。

5)对NDVI与气温、降水进行滞后性分析发现，太行山区NDVI 与气温和降水有明显的时间滞后性关系，其中NDVI与气温的滞后性是1个月，与降水的滞后性是2个月。在进行太行山脉生态系统稳定性研究和太行山区生态环境治理与保护时，应充分考虑NDVI与气温和降水的滞后性关系。