基于延安市不同植被类型NDVI对气象因子影响及空间异质性

张家政，闵志强，王得军，李清顺，孙景梅，李宏韬，李崇贵*

(1.西安科技大学 测绘科学与技术学院，陕西 西安 710054；2.国家林业和草原局西北调查规划设计院 旱区生态水文与灾害防治国家林业局重点实验室，陕西 西安 710048)

植被与气候变化响应研究是目前生态与环境研究中的重要内容之一[1-3]，对土壤系统与大气系统的能量传输具有至关重要的作用[4-5]。目前，研究方法更多是基于植被指数(normalized difference vegetation index，NDVI)进行不同尺度分析。如白旭阳等[6]通过采用对不同驱动因子进行贡献率量化的方法分析1982-2015年气候和人类活动对于新疆玛纳斯流域植被变化的影响，气候变化对于NDVI变化的贡献率高于人类活动对于NDVI变化的贡献率，并且气温和降水对NDVI存在显著性的影响；M.Lamchin等[7]通过分析1982-2014年亚洲地区气象因子(温度、降水和蒸散发等)与植被指数NDVI相关关系，结果表明温度是导致亚洲地区植被指数NDVI变化最主要的因素；许玉凤等[8]从年、季、月等不同尺度研究NDVI动态变化趋势以及对气候因子(日照时数、气温、降水、地表温度等数据)响应，结果显示贵州植被生长是受温度调控，且气温升高、气候变暖是促使植被生长季延长的气候效应之一。很多研究将全区NDVI均值与气象因子做相关性分析[9-10]，有些研究将NDVI在时间尺度作为一个连续变量，与气象因子做空间相关性分析[7,11]。前者只考虑到了相关性关系，但是忽略了NDVI空间变化对于气象因子响应；后者也没有分析不同植被类型与气象因子的响应。因此，为了解不同植被类型NDVI对于气象因子响应，本研究基于长时序Landsat数据集与植被覆盖数据，探索延安市不同植被类型时空变化以及与气象因子(平均温度、总降水、最高温度和最低温度)的响应，并采用最小二乘法+时空地理加权回归模型[12-13]揭示NDVI与气象因子拟合的空间异质性。

Google Earth Engine(GEE)平台不仅提供存储海量的遥感数据，而且广泛用于变化监测[14-15]、土地利用动态监测[16]、农业应用[17]和灾害监测与评估[18]等。GEE提供了JavaScript 和Python的API和基于Web的交互式开发环境，可以轻松地利用Google云的强大空间地理分析功能，使得长时序植被动态分析成为可能。

1 材料与方法

1.1 研究区概况

延安市(107°41′-110°31′E，35°21′-37°31′N)位于我国黄河流域中部，属于黄土丘陵沟壑区，是黄土高原乃至全国水土流失最严重的地区之一。该市总面积为3.67万km2，平均海拔1 200 m。延安地区夏天炎热少雨，冬季干燥寒冷，年均气温为9～12℃，年均降水量在100～800 mm，气温年、日较差大。为遏制严重的水土流失，1999年国家开始实施生态修复工程。

1.2 研究方法

1.2.1 数据来源 Landsat影像 本研究Landsat遥感影像来自于GEE平台(https://earthengine.google.com/)上landsat surface reflectance data(陆地卫星地表反射率数据)。该数据集已经经过大气辐射校正处理，其空间分辨率为30 m，时间分辨率为16 d。首先在GEE平台上利用空间和时间获取1993-2017年延安市Landsat数据集[1993-2011年Landsat 5 TM(‘LANDSAT/LT05/C01/T1_SR’)、2012年Landsat 7 ETM+(‘LANDSAT/LE07/C01/T1_SR’)和2013-2017年Landsat 8 OLI(‘LANDSAT/LT08/C01/T1_SR’)]；然后批量对每一景影像进行去云掩膜处理和裁剪处理，并计算归一化植被指数(NDVI)；最后采用最大值合成法(maximum value composition,MVC)对每年NDVI影像进行最大值合成，以此消除大气和云等干扰误差，得到的年最大NDVImax代表相应各年的NDVI。

气象数据 降水(PRE)、平均气温(Tmean)、最高温度(Tmax)和最低温度(Tmin)数据来源于国家科技基础条件平台——国家地球系统科学数据中心(http://www.geodata.cn)，其空间分辨率为1 km，选取时间序列从1993年1月至2017年12月逐月数据。因数据格式为.nc，所以先利用Python进行数据预处理，使之变成可以使用的栅格图像，再将其合并成年尺度数据[19]。

土地覆盖数据 土地覆盖数据来源于欧洲航天局的1993-2017年共25期土地利用/覆盖(ESA CCI-LC)数据集(http://maps.elie.ucl.ac.be/CCI/viewer.)，空间分辨率300 m。延安市植被覆被类型分为针叶林(coniferous forest)、阔叶林(broad-leaved forest)、灌木/草地(shrub/grassland)、镶嵌林地(mosaic tree)(乔木分布>50%)、镶嵌草地(mosaic grassland)(草本分布>50%)、人类占用地(human occupied land)和水体(water bodies)共7类(图1)。

图1 研究区植被覆盖情况

1.2.2 数据处理方法 线性分析研究表明[20-21]，一元线性回归模型(Slope)能真实有效地反映植被的年际变化。因此，本研究采用一元线性回归对延安市每个像元1993-2017年植被NDVImax进行趋势分析，计算公式为：

(1)

Mann-Kendall(M-K)检验 与以往研究者所采用的基于最小二乘法的时间序列线性回归模型的检验方法相比较，M-K属于非参数检验方法，该方法不需要样本遵循一定的分布，也不受少数异常值的干扰[22]。M-K检验已经在气象、水文等领域的时间序列分析中得到广泛的应用。

重心迁移模型 为深入研究各植被类型时空集聚和迁移特征，本研究采用重心迁移模型进行植被时空演化的定量表达，以此揭示延安市不同区域植被类型活动的演变轨迹。其计算公式如下：

(2)

最小二乘回归模型+时空地理加权回归模型(OLS+GTWR)最小二乘法(ordinary least squares，OLS)常常在遥感变化监测中应用，也是常用的求解无约束最优化问题的回归分析方法之一。具体表达式为：

(3)

式中,y代表因变量，本研究指的是NDVImax；β0为常数项，βj为第j个自变量的回归系数；xij代表第i个样本的第j个自变量；ζ为符合正态分布的随机误差项。

地理加权回归模型(GWR)是空间异质性研究的经典模型，而时空地理加权回归模型(GTWR)则是将时间属性联系到GWR的空间属性中，以更好地反映研究区时空变化信息，使得模型估算结果更加有效[23-24]。具体表达式为：

(4)

式中,yi代表第i个样本点的被解释变量；ui，vi分别为第i个样本点的经、纬度坐标；ti为第i个样本点的时间坐标；β0(ui,vi,ti)为第i个样本点的常数项；βk(ui,vi,ti)为第k个解释变量在第i个样本点的回归系数。

1.2.2.3 独立生活能力训练:每周1次，每次40 min。向患者讲解良好的生活习惯的重要性，对患者日常生活能力进行基础评估，与患者共同制定个体化日常行为规范，督促其按时睡觉、按时起床、整理病床、洗漱、自动排队进食服药及定时参加适当的体育锻炼，组织患者看电视、报刊和杂志。帮助患者建立良好的生活习惯，培养患者有规律的生活和自理能力。

本研究通过借助ArcGIS10.3平台Zonal Statistics工具提取每种植被类型的年NDVImax、PRE、Tmean、Tmax和Tmin等数据，再采用OLS+GTWR模型的方法建立不同土地覆盖类型年NDVImax与各种气象因子的回归模型。

2 结果与分析

2.1 NDVImax和PRE、Tmean、Tmax和Tmin年际变化特征

利用一元线性回归方法分析延安市NDVImax和PRE、Tmean、Tmax和Tmin的年际变化趋势(图2)。NDVImax整体呈现波动增长趋势，NDVImax均值由1993年的0.239增加到2017年的0.386，并在2013年达到最大值0.511，之后有所降低。NDVImax变化率为0.008/a(P<0.01)，植被生长整体呈现较好的生长趋势；PRE、Tmax和Tmin总体呈现波动增长趋势，其增长速率分别为4.619 mm/a(P<0.05)，0.035℃/a(P<0.01)和0.013℃/a，其中波动最大的是Tmax，由1993年的14.01℃增长到2017年15.49℃；Tmean总体则呈现平稳增长趋势，其增长率为0.023℃/a。

图2 1993-2017年延安市NDVImax和PRE、Tmean、Tmax和Tmin年际变化趋势

2.2 土地覆盖面积及土地覆盖变化对植被指数NDVImax的影响

除水体在25 a基本保持不变(20.146 km2)以外，其他土地覆盖类型面积均发生了不同程度的变化(图3a)。针叶林和阔叶林面积均呈现增加趋势，其变化率分别为0.185 km2/a(P<0.01)和0.007 km2/a(P<0.01)。特别的，阔叶林面积变化最大(增加3.85×103km2)；灌木/草地的面积呈现先减少后增加趋势，1993-2007年草地/灌木呈现减少趋势，变化率为-5.279/a(P<0.01)，而在2008年急剧增加为12.345×103km2；镶嵌林地和镶嵌草地面积均呈现减少趋势，其变化率分别为-0.014 km2/a(P<0.01)和-7.223 km2/a(P<0.01)；人类占用地呈现总体增加趋势[0.001 km2/a(P<0.01)]。特别的，在2006年人类占用地面积达到最大值(0.989×103km2)，而在2008年面积由2007的0.989×103km2急剧减少为0.986×103km2，其后保持稳定增加状态。不同土地覆盖类型NDVImax变化介于0.1～0.7，其NDVImax大小分别为针叶林>阔叶林>镶嵌林地>镶嵌草地>灌木/草地(图3b)，且逐年呈现增加趋势。其中针叶林增加幅度最大，分别为0.009 1/a；阔叶林增加幅度最小，为0.007 1/a；灌木/草地、镶嵌林地和镶嵌草地增加幅度一致，为0.008/a。

图3 不同植被覆盖的面积以及NDVImax年变化

2.3 不同植被覆盖类型重心演变轨迹

重心分布变化能够从空间上反映不同植被覆盖类型的时空演变过程。为进一步探讨延安市植被土地覆盖类型重心的空间变化规律，本研究分析了1993-2017年植被覆盖类型重心演变情况(图4)。针叶林重心主要集中在洛川县北部地区，其重心总体向西南部偏移2 516.03 m，其中，2004-2005年重心转移距离最大(776.55 m)；阔叶林重心主要集中在富县西部地区，其重心变化幅度最小，重心总体向东偏移675.96 m，其中，1994-1995年重心转移距离最大(186.40 m)；灌木/草地重心主要集中在安塞县中南部地区，其重心总体向东南偏移1 186.70 m，其中，1994-1995年重心转移距离最大(331.29 m)；镶嵌林地重心主要集中在富县西部地区，其重心变化幅度最大、最复杂，重心总体向南偏移4 157.74 m，其中，2002-2003年重心转移距离最大(1 240.62 m)；镶嵌草地重心主要集中在甘泉县东部地区，其重心总体向西南偏移1 896.36 m，其中2007-2008年重心转移距离最大(1 235.66 m)。

图4 重心转移

2.4 植被指数NDVImax空间变化特征

基于Slope+M-K趋势检验的分析方法得到NDVImax变化趋势的空间分布(图5a)。NDVImax以增加为主，面积占研究区总面积的80.40%[(显著增加(Slope≥0.02)和增加(0.002≤Slope<0.02)的总和)]，远大于NDVImax减少趋势面积[(显著减少(Slope≤-0.02)和减少(-0.02≤Slope<-0.002)的总和)]，表明植被覆盖状况呈现不断改善趋势。NDVImax呈现减少趋势地区主要集中在延安市中部地区(宝塔区南部、甘泉县北部和富县中北部等)。该地区主要以灌木/草地为主，随着城镇化建设的加快推进，部分地区灌木/草地转化为人类占用地加剧了植被退化。显著增加区域主要分布在黄龙县东部、宜川县中北部和延长县南部等区域。

通过对NDVImax进行显著性检验(图5b)，1993-2017年，除1993年外，NDVImax在其余时间段的UF曲线均>0，呈上升趋势，且2004-2017年UF曲线超过95%显著水平线，说明此时间段延安市NDVImax增长趋势显著，植被改善明显。

图5 年NDVImax变化趋势及M-K检验

2.5 全区植被指数NDVI对气象因子的响应

本研究基于像元对植被年NDVImax与各气象要素进行相关性分析(图6)。图6a的均值为0.262，从像元尺度来看，研究区NDVImax与PRE呈正相关的面积89.53%，主要分布在安塞县、延川县、延长县以及宜昌县北部地区，尤其在延川县延川镇和关庄镇、安塞县张坪乡等地NDVImax与PRE之间相关性较高，呈极显著正相关(P<0.01)，PRE是这些地区植被生长的主要限制因子，其他地区植被与降水相关性不明显；图6b的均值为0.085，从像元尺度来看，研究区NDVImax与Tmean呈正相关的面积71.38%，主要分布在延川县中部、宜川县南部、黄龙县中部以及洛川县中部等地区。研究区NDVImax与Tmean呈显著负相关(P<0.01)的地区占总面积的6.02%，主要分布在志丹县南部桥镇和永宁镇等地区，Tmean是影响该地区植被生长的主要限制因子；图6c的均值为0.119，研究区NDVImax与Tmin呈正相关的面积75.57%，主要分布在子长县中南部、延川县西部和延长县中北部等地区，尤其在延川县贾家坪镇、子长县王家湾镇等地NDVImax与Tmin之间相关性较高，呈极显著正相关(P<0.01)，占研究区总面积的1.61%；图6d的均值为0.033，研究区NDVImax与Tmax呈正相关的面积57.02%，主要分布在洛川县中部、富县中南部以及黄龙县。研究区NDVImax与Tmax呈极显著负相关(P<0.01)的地区占总面积的3.02%，主要分布在志丹县南部等地区。

图6 1993-2017年NDVImax与PRE(a)、Tmean(b)、Tmin(c)和Tmax(d)相关系数空间分布

有研究表明，水热条件的增加会促进植被生长，且气候的变化对植被生长影响存在一定的滞后性[23]，这些因素都有可能造成气候变化对植被无显著性影响，也有可能是不同植被覆盖对气候变化存在差异影响，从而掩饰了气候变化对于植被的作用。

2.6 不同植被类型NDVImax对气象因子的响应机制

在不同植被覆盖类型区气象因子对各种植被覆盖类型的作用有明显差异。由表1知，针叶林与Tmean和Tmin显著性相关，年NDVImax与Tmean和Tmin呈正相关，相关系数分别为0.431(P<0.05)和0.505(P<0.05)；阔叶林仅与Tmin显著性相关，年NDVImax与Tmin呈正相关，相关系数分别为0.396(P<0.05)；灌木/草地和人类占用地与PRE和Tmin显著性相关，年NDVImax与PRE和Tmin呈正相关，相关系数分别为0.402(P<0.05)、0.409(P<0.05)和0.397(P<0.05)、0.408(P<0.05)；镶嵌林地和镶嵌草地仅与Tmin显著性相关，年NDVImax与Tmin呈正相关，相关系数分别为0.433(P<0.05)和0.414(P<0.05)；水体与各气象因子均无显著性关系，Tmax与各种覆盖类型年NDVImax均无显著相关性，说明Tmax在年尺度上对不同类型植被作用均不显著。

表1 不同植被类型NDVImax与气象因子相关系数

对不同植被覆盖类型年NDVImax与气象因子建立OLS模型，剔除方差膨胀因子(Variance Inflation Factor，VIF)>7.5的气象因子，然后采用不同植被类型NDVImax的变化量作为因变量，剔除完VIF>7.5的因子作为自变量，对NDVImax的变化进行GTWR分析(表2)。相比OLS模型，GTWR模型对于各植被类型与气象因子的拟合效果显著提高，不同植被类型的GTWR模型拟合平均R2(0.419)远大于OLS模型的拟合平均R2(0.100)，与OLS模型的AICc之差(ΔAICc=5254.78)远>3，二者拟合差异明显，GTWR模型拟合度高于OLS模型。

表2 GTWR与OLS模型结果比较

对GTWR模型各因子拟合系数采用箱线图描述(图7)。针叶林箱线图中，PRE、Tmean、Tmax和Tmin对NDVImax平均回归系数分别为0.040、-0.017、0.030和0.075，表明PRE、Tmax和Tmin对针叶林NDVImax呈正向影响，Tmean对针叶林NDVImax呈负向影响；上下四分位数范围分别为-0.053～0.126、-0.152～0.081、-0.076～0.135和0.024～0.136，表明各年份PRE和Tmin密度系数变化幅度相对较大，而Tmean和Tmax密度系数变化差异不显著；阔叶林箱线图中，PRE、Tmean、Tmax和Tmin对NDVImax平均回归系数分别为-0.037，-0.008，-0.127和0.112，表明PRE、Tmean和Tmax对阔叶林NDVImax呈正向影响，Tmin对阔叶林NDVImax呈负向影响；上下四分位数范围分别为-0.093～0.077、-0.067～0.038、-0.185～-0.076和0.065～0.158，表明各年份PRE和Tmin密度系数变化幅度相对较大，而Tmean和Tmax密度系数变化差异不显著；灌木/草地箱线图中，PRE、Tmax和Tmin对NDVImax平均回归系数分别为0.087、0.011和0.030，中位数分别为0.075、0.013和0.033，表明PRE、Tmax和Tmin对灌木/草地NDVImax呈正向影响；上下四分位数范围分别为0.06～0.168、-0.074～0.094和-0.038～0.092，表明各年份PRE密度系数变化幅度相对较大，而Tmax和Tmin密度系数变化差异不显著；镶嵌林地箱线图中，PRE、Tmean、Tmax和Tmin对NDVImax平均回归系数分别为-0.017、-0.034、0.050和0.112，表明PRE和Tmean对镶嵌林地NDVImax呈正向影响，Tmax和Tmin对镶嵌林地NDVImax呈负向影响；上下四分位数范围分别为-0.035～0.063、-0.038～-0.030、0.038～0.06和0.111～0.113，表明各年份PRE、Tmean、Tmax和Tmin密度系数变化差异不显著；镶嵌草地箱线图中，PRE、Tmax和Tmin对NDVImax平均回归系数分别为0.061、-0.094和0.097，表明PRE和Tmin对镶嵌草地NDVImax呈正向影响，Tmax对镶嵌草地NDVImax呈负向影响；上下四分位数范围分别为-0.008～0.117、-0.216～-0.004和0.025～0.167，表明各年份Tmax密度系数变化幅度相对较大，而PRE和Tmin密度系数变化差异不显著。

图7 GTWR拟合回归系数箱型

3 结论与讨论

随着国家对延安市实行大规模退耕还林等生态修复工作以来，延安市植被动态变化监测已成为目前研究的热点[25-26]，但目前很少研究不同植被覆盖类型变化及其与气象因子关系，由于单一植被生态系统在水土保持、调节气候等方面与多植被类型存在很大差距，所以研究不同植被覆盖类型变化具有较深的意义。对此，本研究定量探究了不同植被覆盖类型变化及其与气象因子关系，揭示了不同植被类型对气象因子的响应规律。与MODIS数据相比，Landsat卫星具有时间跨度长、空间分辨率高(30 m)的优势，能够为长时序植被动态监测提供更好的数据支撑。GEE平台可以提供大量的数据，这样为今后长时序、大幅度研究提供了便捷。

本研究基于ESA CCI-LC植被覆盖分类数据以及Landsat数据集，使用Slope+Manna-Kendall非参数检验方法，分析全区和不同植被覆盖类型区内NDVImax与气象因子时间变化特征与变化趋势，并应用重心迁移模型分析其空间变化特征，最后应用最小二乘法和时空地理加权回归模型进一步揭示不同植被覆盖类型区内NDVImax分异特征及其驱动力，得到以下结论。

1)1993-2017年，延安市全区植被年NDVImax整体呈现波动增长趋势，NDVImax均值由1993年的0.239增加到2017年的0.386，针叶林、阔叶林、灌木/草地面积呈现增加趋势，而镶嵌林地和镶嵌草地面积则呈现减少趋势。不同植被覆盖类型年NDVImax也呈现增长趋势，其中针叶林增加幅度最大，为0.009 1/a。这与朱会利等[25]、何立恒等[26]的研究结果相一致，这表明自延安市生态修复工程以来，植被覆盖状况极大改善。

2)延安市全区NDVI与气象因子(Tmean、PRE、Tmax和Tmin)之间没有显著相关性，但在不同植被覆盖类型情况下，气象因子对NDVImax存在显著性作用。这说明对植被类型进行分类更加有利于揭示气象因子对于植被的作用机制。其中，针叶林年NDVImax与Tmean和Tmin呈现显著性正相关，相关系数分别为0.431(P<0.05)和0.505(P<0.05)；阔叶林年NDVImax仅与Tmin呈现显著性正相关，相关系数分别为0.396(P<0.05)；灌木/草地和人类占用地年NDVImax与PRE和Tmin呈现显著性正相关，相关系数分别为0.402(P<0.05)、0.409(P<0.05)和0.397(P<0.05)、0.408(P<0.05)；镶嵌林地和镶嵌草地年NDVImax仅与Tmin呈现显著性相关，相关系数分别为0.433(P<0.05)和0.414(P<0.05)。

3)针叶林主要分布在洛川县，迁移速率分别为100.64 m/a；灌木/草地重心主要分布在安塞县，迁移速率为47.46 m/a；镶嵌草地主要分布在甘泉县，迁移速率为75.85 m/a；镶嵌林地和阔叶林主要分布在富县，迁移速率为分别为166.31 m/a和27.04 m/a。

4)Tmean对阔叶林和镶嵌林地NDVImax呈现正向影响，对针叶林NDVImax呈现负向影响；PRE对各种植被类型NDVImax均呈现正向影响；Tmax对针叶林、阔叶林和灌木/草地NDVImax呈现正向影响，对镶嵌草地和镶嵌林地NDVImax呈现负向影响；Tmin对针叶林、灌木/草地和镶嵌草地NDVImax呈现正向影响，对阔叶林和镶嵌林地NDVImax呈现负向影响。