基于MODIS数据的珠三角地区NDVI时空变化特征及对气象因素的响应

2019-11-04何全军

生态环境学报 2019年9期

何全军

1. 广州气象卫星地面站，广东广州 510640；2. 广东省生态气象中心，广东广州 510640

植被是地球表面植物群落的总称，作为地表状况的重要表征，在地球系统中扮演着重要的角色，能影响地气系统的能量平衡。植被生长过程具有明显的年纪变化和季节性变化，是气候和人文因素对环境影响的敏感指标（Bounoua et al.，2000；郭铌，2003；朱建华等，2007），卫星遥感可以完成大范围的陆表植被生态环境的持续监测，从而能够揭示全球或区域环境状况的演变（Zhao et al.，2018；Guan et al.，2018；Chen et al.，2019）。通过遥感手段获得的归一化植被指数 NDVI（Normalized Difference Vegetation Index）（Tucker，1979；Hilker et al.，2015）已被作为一种直观量化指标，广泛用于定性和定量评价植被覆盖及生长活力，以及进行不同尺度上的植被物候、生态环境和气候变化的研究（刘绿柳等，2006；刘世梁等，2016；Liu et al.，2019）。然而，目前针对 NDVI与气候变化的相互关系研究中，不同区域的研究结论并不完全一致，甚至有的观点完全相反。白淑英等（2012）发现长江流域NDVI对气温变化的响应大于降水，孔冬冬等（2017）的研究表明温度对青藏高原植被物候的影响占主导地位，而曲学斌等（2018）则认为呼伦贝尔草原NDVI年变化主要受降水驱动，与气温和总辐射的相关性较小。此外，刘艳等（2018）发现不同草地类型的水热影响差异显著，天山草地约41%的区域与降水量呈显著正相关，31%区域与温度呈显著负相关。以上结果表明NDVI与气象因素的相关关系存在重大的区域差异。

珠江三角洲（Pearl River Delta，PRD）地区是中国城市化水平最高、经济最发达、经济增长速度最快的地区之一，但是，随着人口的迅速膨胀，以及城市化和工业化的迅猛发展，珠三角地区的土地利用发生巨大变化（黎夏，2004；李慧等，2011），大量耕地、林地和草地转变为城镇工交用地，植被覆盖迅速减少（闫小培等，2006）。该区域的地表植被时空变异性极明显，分布格局更复杂，出现了一系列生态环境和安全问题（王柳等，2010；徐庆勇等，2011；胡志仁等，2018），阻碍了经济的可持续发展以及影响了人类居住环境的适宜性。植被是陆地生态系统的主要组分，是生态系统变化的指示器，生态环境保护首先是地表植被的保护。植被的种类、数量和分布是衡量区域生态环境安全及宜居性的重要指标。因此，开展绿化行动、恢复植被、修复生态是珠三角地区的重点工作。

本研究利用 2001－2017年覆盖珠三角地区的卫星遥感NDVI时间序列及气象站地面观测数据，进行研究珠三角区域NDVI的时空变化趋势及空间自相关特征，分析NDVI对区域气象因素变化的响应规律，为利用遥感手段进行珠三角地区植被生态环境监测提供有力方法，对开展生态环境保护以及促进生态文明建设提供科学建议。

图1 珠三角行政区及气象站分布Fig. 1 Distribution of PRD administrative region and meteorological stations

1 研究区概况

珠江三角洲，简称珠三角，是由西江、北江和东江入海时冲击沉淀而成的一个三角洲，位于中国广东省中南部、珠江入海口处，毗邻港澳，与东南亚地区隔海相望。珠三角既是地理区域，也是经济区域，它是拥有全球影响力的先进制造业基地和现代服务业基地，是中国人口集聚最多、创新能力最强、综合实力最强的三大城市群之一。珠三角地区范围内的主要城市的行政规划包括广州、深圳、佛山、东莞、中山、珠海、江门等7个市以及惠州市和肇庆市的部分地区（如图 1所示），总面积约41000 km2。珠三角大部分地区位于北回归线以南，属亚热带海洋季风气候，雨量充沛、热量充足，区域地带性植被为南亚热带季风常绿阔叶林。珠三角地区南面临海，东、北、西三面被山地、丘陵包围，中部平原上水系丰富。因此，珠三角地区的水热条件优越，日照时间长，植被覆盖度较高。

2 资料与方法

2.1 数据来源与预处理

2.1.1 NDVI数据

NDVI数据来自美国 LAADS（Level-1 and Atmosphere Archieve & Distribution System）DAAC（Distributed Active Archive Center）发布的2001－2017年 MODIS（Moderate Resolution Imaging Spectroradiometer）数据 MOD13A3 产品（https：//ladsweb.modaps.eosdis.nasa.gov/），该产品是采用限制视角的最大值合成（Constrained View angle Maximum Value Composite，CV-MVC）月产品。MOD13A3是正弦曲线投影，需要将其转换为等经纬度投影，并利用珠三角地区的矢量边界数据进行裁剪。

MOD13A3数据中存在NDVI的低异常值及空白无效值，需要对其进行重构处理。通过对比NDVI时间序列谐波分析方法HANTs（Harmonic Analysis of NDVI Time-Series）（Roerink et al.，2000）、修订的最佳指数斜率提取方法 M-BISE（Modified Best Index Slope Extraction）（Lovell et al.，2001）、S-G滤波方法（Savitzky-Golay filter）（Chen et al.，2004）以及基于数据质量标志数据集的样条曲线插值方法（SPLINE）的重构结果，最终采用 M-BISE和S-G滤波的组合方法实现月NDVI时间序列重构，并在重构的月时间序列数据基础上实现年时间序列产品以及多年合成产品。

2.1.2 气象数据

气象数据为广东省气象局提供的 2001－2017年珠三角及周边地区的国家地面观测站的逐月的平均气温、累积降水以及日照时数资料，气象站点分布如图1所示。克里金（Kriging）方法是气象要素空间插值的常用方法（李军龙等，2006；于海敬等，2019），本文采用该方法对站点数据进行插值，对气温插值时考虑了高程影响（潘耀忠等，2004），气温垂直递减率设置为0.0065 ℃·m-1。气象数据插值时设置的空间分辨率与NDVI一致，并按照珠江三角洲的范围对插值结果进行裁剪，最后进行多年的月和年数据的合成处理。

2.2 研究方法

2.2.1 变化趋势分析

线性回归方程被广泛用于植被指数在时间上的变化速率的评估（Stow et al.，2003；武永利等，2008；Baniya et al.，2018），采用常规的数理统计和线性回归拟合方法分析珠三角地区植被的年际变化趋势。此外，逐像元构建NDVI随年际变化的线性回归方程，分析珠三角地区植被变化趋势的空间分布特征。

由于年合成数据消除了NDVI的季节内变化影响，体现了NDVI的年际总生长趋势，而植被生长具有显著的季节变化特征（何全军等，2008），通过移动平均方法即月NDVI与年NDVI的比值获得季节指数 SI（Seasonal Index）（Ittig，2004；贾俊平等，2009）来分析珠三角地区植被的季节变化趋势。SI是进行时间序列数据变化趋势预测的重要手段，它能刻画时间序列在一个年度内各月份的典型季节特征。

变异系数CV（Coefficient of Variation）是标准差与平均值的比值，是衡量资料中各观测值变异程度的统计量（毛德华等，2012），通过逐像元计算珠三角地区年均值合成NDVI的CV，分析珠三角地区植被时间序列的稳定性。

2.2.2 空间自相关分析

根据地理学第一定律（Tobler，1970），地物之间的相关性与距离有关，距离越近，地物间相关性越大；距离越远，地物间相异性越大。莫兰指数Moran's I（Moran，1950）是应用最广泛的度量空间自相关的指标之一，可以反映空间邻接或临近的区域单元植被参数特征的相似程度，又分为全局Moran's I和局部 Moran's I（Anselin，1995）。全局Moran's I是对研究区域内属性数据空间自相关的综合反映和度量，但就区域内部而言，各局部区域的空间自相关并非完全一致，而是常表现出不同性质与程度的空间异质性。因此，Anselin（1995）提出空间局部自相关分析方法，并定义了局部Moran's I，用 Moran 散点图和 LISA（Local Indicators of Spatial Association）集聚图来描述局部空间相关特性。利用空间数据分析软件GeoDa（Anselin et al.，2005）对珠三角地区的年NDVI及多年合成NDVI进行全局和局部空间自相关分析，分别通过Moran's I、Moran散点图以及LISA集聚图来分析珠三角地区植被的空间集聚程度及变化特征。

图2 珠三角地区年平均NDVI与SI的时间序列趋势Fig. 2 Time series trends of annual average NDVI and SI in PRD

2.2.3 相关性分析

皮尔森相关系数（Pearson Correlation Coefficient）是被广泛用来反映两个变量线性相关程度的统计量，其定义为两个变量之间的协方差和标准差的商（仲晓春等，2016）。分别针对月和年时间序列数据计算NDVI与气温、降水及日照时数之间的皮尔森相关系数，并进行显著性检验，分析珠三角地区植被生长对气候变量的响应。

3 结果与分析

3.1 珠三角地区NDVI的时间序列变化特征

图 2反映了 2001－2017年珠三角地区 NDVI在年和月时间尺度的变化趋势。其中图 2a是利用年合成数据统计得到的年平均NDVI的变化曲线，可以看出珠三角地区的 NDVI呈现“升-降-升”的持续波动的年际上升趋势。根据一元线性拟合得到的趋势方程可知，2001－2017年中珠三角地区的NDVI整体上的平均年增长速率为0.0051，说明珠三角地区的植被一直在缓慢增长。图 2b是珠三角地区 SI的统计变化曲线，SI的最低值和最高值分别出现在2月和9月。SI通过振幅变化的相对程度来描述植被随季节变化的强度的变化趋势，因此SI在如实反映NDVI季节变化特征的同时，还能反映NDVI季节变化的增长幅度随着时间变化的趋势。从总趋势来看，2001－2017年SI的振幅呈现越来越小的趋势，说明随着时间的发展，植被生长在遵循自然周期规律的同时季节内变化差异逐渐减小。此外，2004年NDVI出现最小值，SI的振幅在2005年出现最大值，其原因与 2004年广东省发生特大干旱有关，2004年9月下旬开始广东发生严重秋冬连旱，抑制了植被的生长，直到 2005年汛期之后旱情缓解（罗晓玲，2005），持续降水及气温回升使植被恢复后快速生长。

由于下垫面的土地利用类型不同以及区域发展程度不同，珠三角范围内不同地区的植被分布和变化存在空间差异性（图3）。图3a是珠三角地区的多年NDVI合成结果，可以看出珠三角地区呈现显著的中部低NDVI被周边高NDVI围绕的空间分布特征，这与图1所示的珠三角城市和地形的分布一致。中心区域属于城市群集区，NDVI偏低，随着向外围扩展，NDVI逐渐升高。图3b是珠三角地区的NDVI增长趋势的空间分布。珠三角范围内大部分地区的NDVI增长趋势表现为正值，只有小部分地区为负值。增长趋势为正值说明植被越来越好，负值说明植被发生退化。图 3b中，中山市北部、肇庆四会市区、广州花都区、东莞市珠江沿岸和中部地区、深圳西北部以及佛山靠近广州的地区出现了不同程度的植被退化。对增长趋势进行统计，珠三角地区植被增长的区域占总范围面积的90.84%，其中趋势值大于 0.01的区域占总面积的7.12%；植被退化的区域占总范围面积的9.16%，其中趋势值小于-0.01的占 0.19%；整个珠三角地区NDVI的区域趋势均值为0.0051，与图2结果一致，说明植被状况好转是整个珠三角地区的总趋势。图3c是珠三角地区年NDVI时间序列数据的CV空间分布图。CV值越大，表明数据分布越分散，说明该地区植被的波动较大；反之则表明数据分布越紧凑，说明该地区植被时间序列较为稳定（张雪艳等，2009）。对CV进行分级统计的数据显示，珠三角地区NDVI的CV小于0.1的区域占83.80%，大于0.1的占16.20%（其中大于0.2的占0.64%），说明大部分地区的植被变化幅度不大，植被时间序列较为稳定。但从CV的空间分布来看，围绕珠江口周边分布的城市集聚区植被变化幅度最大，说明该地区植被时间序列稳定性较差。此外，通过对比图3b、图3c可知，NDVI的增长趋势和CV的高值区域并不能完全对应。增长趋势是对长期植被变化结果的描述和未来预测，反映的是一种累积的结果，而 CV则可以将过程中发生过的显著的变化记录下来，反映了长期过程中曾出现的阶段性变化。

图3 珠三角地区多年平均NDVI（a），增长趋势（b）和CV（c）的空间分布特征Fig. 3 Spatial distribution of multi-year average NDVI (a), growth trend (b) and CV (c) in PRD

图4 珠三角地区年NDVI空间自相关特征变化趋势Fig. 4 Trends of spatial autocorrelation characteristics of annual NDVI in PRD

3.2 珠三角地区NDVI的空间自相关特征

在空间数据分析软件GeoDa 1.12中选择Queen邻接类型建立空间权重矩阵，分别对 2001－2017年的年NDVI时间序列数据进行单变量全局和局部Moran's I分析，并统计珠三角地区年NDVI空间自相关特征的变化趋势（图4）。Moran's I的值在-1和1之间，正值表示空间正相关，值越大说明空间相关性越明显，而负值表示空间负相关，值越小则空间差异越大。图4a是2001－2017年全局Moran's I变化趋势，可以发现珠三角地区历年的全局Moran's I都在0.91以上，说明长期以来珠三角地区植被都表现为较强的空间集聚正相关。虽然年际变化上存在一定的波动性，但从其变化趋势来看，珠三角地区植被的空间格局保持向高集聚发展的趋势。图4b是LISA空间相关类型的分类数据的时间变化趋势。其中“High-High”表示高NDVI集聚区，“Low-Low”表示低NDVI集聚区，“Low-High”表示低NDVI被高NDVI包围，“High-Low”表示高NDVI被低 NDVI包围，“Not Significant”表示未通过P=0.05显著性检验的区域。由曲线变化图可以看出，珠三角地区NDVI的空间自相关性主要集中在“High-High”、“Low-Low”和“Not Significant”3个类型，仅有很少量的“Low-High”和“High-Low”类型（图中对“Low-High”和“High-Low”两种类型的百分比放大了1000倍）。对比图中5种空间类型的变化曲线，“Low-Low”类型的时间变化不明显，“High-High”类型稍有升高，“Not Significant”类型持续降低，说明珠三角地区的低植被集聚区的变化较为稳定，而高植被集聚区的范围增加了，这种增加在很大程度上可以归因于“Not Significant”类型的转变。“High-Low”和“Low-High”两种类型所占的百分比虽然小但变化却比较活跃，表明这两种类型所代表的区域的植被变化具有较高的不稳定性。

图5是对珠三角地区多年平均NDVI进行局部空间自相关分析得到的Moran散点图和LISA集聚图，可以直观展示出珠三角地区NDVI的局部空间集聚的规律和空间分布格局。图5a中Moran散点图的X坐标和Y坐标分别表示标准化后的NDVI和空间滞后NDVI（Lagged NDVI），4个象限分别表示4种空间关系：右上和左下的I和Ⅲ象限为相似值的空间聚类，其中右上为“High-High”，左下为“Low-Low”；左上和右下的Ⅱ和Ⅳ象限为不同值的空间关系，左上为“Low-High”，右下为“High-Low”。珠三角地区的 NDVI主要位于“High-High”和“Low-Low”象限，“Low-Low”象限的散点数目明显多于“High-High”象限，说明珠三角地区的植被在空间分布上集聚度高，并且低值区的空间集聚度高于高值区的集聚度。此外，珠三角地区多年平均NDVI的全局Moran's I是0.932765，说明珠三角地区的植被在空间上呈现显著的正相关，整体集聚度高。

图5b的LISA集聚图从空间分布上反映NDVI的空间集聚特征，在通过 P=0.05显著性检验条件下，珠三角的中心城市区表现为显著的低值自相关特征，即低植被集聚，最外围表现为高值自相关。对代表低植被集聚的“Low-Low”类型以及代表高植被集聚的“High-High”类型的像素进行统计，其中“Low-Low”类型占珠三角地区总面积的28.77%，“High-High”类型占33.69%，二者之和达到62.46%，充分体现了珠三角地区植被覆盖在空间关系上呈现较高的集聚度。虽然“High-High”类型的空间集聚度比“Low-Low”类型低，分布格局相对破碎，但其占据区域的面积大于“Low-Low”类型，说明珠三角地区植被覆盖总体良好。此外，在低值集聚区和高值集聚区之间是未通过显著性检验区域（“Not Significant”）以及少量的“Low-High”以及“High-Low”，这部分区域类型如何转变是影响珠三角地区植被生态环境变化的关键。

图5 珠三角地区多年平均NDVI的Moran散点图（a）和LISA集聚图（b）Fig. 5 Moran scatterplots (a) and LISA cluster map (b) of multi-year average NDVI in PRD

3.3 珠三角地区NDVI与气象因素的相关性

对2001－2017年珠三角地区的月NDVI与气温、降水及日照时数计算皮尔逊相关系数，并进行显著性检验。考虑到植被对气候变化的响应存在滞后性，从当前月开始共计算了滞后时差6个月的月NDVI与各气象因素的相关系数（表1）。相关系数的值介于-1和1之间，正值表示正相关，负值表示负相关。相关系数的绝对值越大，相关度越强；相关系数越接近 0，相关度越弱。从计算结果可以发现，珠三角地区的NDVI和气温、降水以及日照时数的月变化显著相关。通过对相关系数的大小进行比较可以发现，NDVI与3种气象因子的相关性从大到小依次为气温、降水和日照时数，说明在月时间尺度上，气温是影响珠三角地区植被生长的最主要因素，其次是降水，最后是日照时数。

从滞后时间来看，NDVI与气温的相关系数的最大值出现在滞后1个月，随后相关系数降低，滞后4个月时呈现负相关；NDVI与降水的最大相关系数出现在滞后1个月，但与滞后2个月的相关系数相差不大，第3个月开始有明显的降低，并且降低的速率低于气温，直到滞后5个月时开始出现负相关；NDVI与日照时数的最大相关系数出现在当前月，第3个月开始为负相关。此外，NDVI与气温、降水量从当前月到滞后3个月的相关系数都通过了P=0.05的显著性检验，而NDVI与日照时数仅在当前月通过显著性检验。以上分析说明珠三角地区NDVI对气温和降水的响应分别存在1个月和1－2个月的滞后，气温和降水对NDVI的影响持续时间较长，其中气温的影响变化幅度大，而降水的影响变化相对平缓。NDVI对日照时数的响应不存在滞后，日照时数对NDVI的影响持续时间短。

珠三角地区的NDVI对气温、降水和日照时数的响应除在时间上存在差异之外，在空间分布上也存在区域差异，图6是NDVI与3种气象因素之间的最大相关系数所对应的滞后月数的空间分布。在珠三角范围内，NDVI与气温的最大相关系数主要出现在当前月及滞后1个月，而NDVI对降水的响应较气温迟钝一些，最大相关系数主要出现在滞后1个月和2个月，并且在空间上都表现出沿“西南-东北”走向滞后月数逐渐变大的现象。在整个珠三角地区范围内NDVI与日照时长的最大相关系数的滞后时间都为0个月，在空间分布上不存在区域差异。

表 2是珠三角地区的年平均 NDVI和年最大NDVI分别与年平均气温、年总降水以及年总日照时数之间的相关系数。可以发现，珠三角地区年NDVI和气象因素之间相关性很弱，并且都未能通过P=0.05的显著性检验，说明在年尺度上，珠三角地区NDVI的变化与气候因素的相关性不显著。但是，仅从相关系数来看，降水是3个气候因素中唯一表现有滞后效应的因素，年总降水会对当年和次年的植被具有一定影响，且次年的影响较大。

表1 珠三角地区月NDVI与气温、降水量及日照时数的相关系数Table 1 Correlation coefficients between monthly NDVI and air temperature, precipitation and sunshine duration in PRD

表2 珠三角地区年NDVI与气温、降水及日照时长的相关系数Table 2 Correlation coefficients between annual NDVI and air temperature, precipitation and sunshine duration in PRD

图6 珠三角地区的月NDVI与气温（a）、降水（b）和日照时长（c）的滞后效应空间分布Fig. 6 Spatial distribution of lag effect between NDVI and air temperature (a), precipitation (b) and sunshine duration (c) in PRD

4 讨论与结论

4.1 讨论

植被变化是一个自然与人类活动交互作用的过程，基于遥感数据的NDVI是综合了所有自然条件和人为因素共同作用于植被生长的客观真实反映，对植被覆盖在时间和空间上的变化监测具有极大优势。珠三角地区NDVI的整体生长趋势在年际上向好发展，年平均增长速率为0.0051，但在增长过程中具有显著的“升-降”波动变化。在空间上，珠三角地区植被增长的区域占总面积的90.84%，有9.16%的地区发生植被退化。植被退化区主要集中在珠三角的核心城市区域，这与城市化和工业化过程中土地利用类型的转变有关（闫小培等，2006），王柳等（2010）认为珠三角核心地带植被地表由于强烈的人为干扰退化较明显。李慧等（2013）指出，珠三角地区经济发展过程中大量耕地转为经济价值更高的建设用地，同时耕地和建设用地的扩展和收缩具有阶段性，这也是导致NDVI时间序列产品不稳定的原因。

植被是具有一定空间分布的区域变化量，空间自相关分析方法能从空间关联角度对相邻地物的相似程度进行判断，有助于揭示植被覆盖的空间集聚模式和规律。李慧等（2011）认为空间自相关分析方法可直观地显示分析结果，是一种深入理解区域土地利用格局和变化的重要和有益的方法。珠三角多年平均NDVI的全局Moran's I达到0.932765，表现为高空间集聚状态，呈现显著的中部低 NDVI被周边高NDVI围绕的空间分布特征，与珠三角的城市群和地形的分布一致（王柳等，2010）。这种植被覆盖的集聚模式与珠三角地区的生态脆弱区分布一致（徐庆勇等，2011；胡志仁等，2018），同时体现了地形对NDVI分布具有控制作用（潘霞等，2019）。

NDVI和气象因素的关系与研究区所处的地理位置、时间尺度以及下垫面类型有很大关系（白淑英等，2012；刘艳等，2018；曲学斌等，2018），气候变化已对植被物候产生显著影响（孔冬冬等，2017；闫敏等，2019）。对珠三角地区而言，植被生长受多种气象因素的综合影响，NDVI和气候因素的年际变化相关性很弱，而年内变化月尺度上NDVI与气温、降水和日照都具有显著相关，同时NDVI对不同气候因素的响应存在时间和空间上的差异。虽然 MODIS数据对判断珠三角地区 NDVI的发展趋势以及识别地区间的变化差异不存在影响，但在对更长期的植被的生长产生影响的气候变化分析方面存在不足。通过使用不同数据进行时间序列数据的扩展，有助于对长期NDVI与气候因子之间的关系进行更深入的探讨（毛德华等，2012）。此外，植被覆盖度变化不仅受气候因素的影响，还受到人为活动的影响（张凯选等，2019），研究珠三角区域气候变化和人为活动对植被变化的贡献也是以后的研究方向。