基于GEE和PLSR的2001-2020年长三角地区植被覆盖时空变化及驱动力分析

李星媚,严 军

(南京林业大学 风景园林学院,江苏 南京 210037)

植被覆盖是陆地生态系统的基础和最重要的组成部分[1-2],受气候变化和人类活动的强烈影响,植被覆盖可以迅速响应周围环境的变化。近年来,随着全球气候变暖和人类活动的不断加剧,地面植被生态系统经历了前所未有的破坏[3-4],因此植被覆盖的动态变化及其对环境的响应成为了生态学等领域的研究热点[5]。

遥感技术能够实现长时序、大面积的植被监测。相比于易饱和的归一化植被指数NDVI,植被增强型指数EVI能够进一步降低噪声影响,更好地表征高植被覆盖区的植被生长状况,并且与部分气候因子的相关性更显著[6-7]。越来越多的学者开始利用EVI数据对区域植被覆盖变化进行研究。王虎威等[8]基于EVI数据分析了山西省2000-2014年植被动态变化及与气候因子的响应。肖骁等[9]通过一元线性回归方法,定量估算2000-2015年东北老工业区的植被覆盖度时空变化特征,并定性分析其与城市化水平的相关性。马昊翔等[10]利用2005-2014年的EVI数据,结合相关分析、残差趋势分析等方法,对青海省生长季草地植被覆盖动态变化及其驱动因素进行分析。然而,一方面,当前遥感影像数据获取成本高、更新慢、分辨率低[11],处理相对费时。google earth engine(GEE)云平台能够提供一系列免费、高精度、长时序的遥感数据集,并且其强大的云计算能力极大地缩短了影像处理时间,提高了工作效率[12]。另一方面,目前对植被覆盖变化的驱动研究多侧重于定量分析其与气候因子的响应关系,对人类活动响应的研究较少,且方法上多为定性或以简单回归为基础的定量研究。偏最小二乘回归(PLSR)分析能够在不排除多重共线变量的前提下有效反映出多个自变量或因变量之间的相关性[13],是定量分析植被覆盖变化对气候、人类活动响应的可行方法。

长江三角洲(简称长三角)地区是我国经济发展最活跃的地区之一。近几十年,尤其是2000年以来,长三角区域的气候变化和人类活动均对植被覆盖和生态环境产生了较大的影响[14-15]。因此,本研究借助GEE平台研究2001-2020年长三角地区EVI的时空演变特征,并结合PLSR模型定量分析气候和人类活动对植被覆盖的影响,以期为长三角地区的植被覆盖监测、生态环境保护与建设提供参考。

1 材料与方法

1.1 研究区概况

长三角(27°02′-35°08′N,114°54′-123°10′E)位于长江下游地区,以冲积平原为主,地势低平,平均海拔137.29 m。行政区包含上海市、江苏省、安徽省和浙江省,共41个城市,总面积为35.8万 km2(图1)。区域属于亚热带季风气候,四季分明,雨热同期。截至2020年底,长三角地区总人口为2.35亿人,生产总值为24.45万亿元[16],常住人口城镇化率超过60%[17]。

图1 研究区地形Fig.1 Topographic map of the study area

1.2 数据来源及处理

研究区域矢量图来自中国科学院资源环境科学与数据中心(https://www.resdc.cn/)。EVI遥感影像数据来自Google Earth Engine(GEE)云计算平台(https://earthengine.google.com)的MOD13Q1数据集,时间分辨率为16 d,空间分辨率为250 m,已经过几何、大气校正等预处理工作。在GEE平台应用最大值合成法[18]和移动窗口平滑算法得到逐月EVI数据,能够减少云雨等低值噪声产生的影响。采用生长季(4-9月)EVI的平均值代替年度EVI,得到EVI年际变化。土地利用数据来源于武汉大学遥感与信息工程学院黄教授团队制作的1985-2020年中国30 m年度土地覆盖产品。该数据集精度高,质量好,改善了时空一致性问题[19]。本研究选取2001、2020年2期土地利用数据,并通过重采样统一分辨率。高程数据采用NASA SRTM Digital Elevation 30 m分辨率产品[20]。

本研究所涉及的气候以及人类活动因子数据均来源于国家统计局(http://data.stats.gov.cn/)、中国气象科学数据共享服务网(http://data.cma.cn/),以及各地国民经济和社会发展统计公报。

1.3 研究方法

1.3.1 Theil-Sen Median趋势分析和Mann-Kendall显著性检验 Theil-Sen Median趋势度是经过计算序列的中值,抗噪性强,但本身不能实现序列趋势显著性判断,而Mann-Kendall显著性检验法[21]可用于检验时间序列变化趋势的显著性,对样本分布无要求且对异常值不敏感[22]。因此两者结合已成为判断植被、气候等长时间序列数据趋势的重要方法[23]。

Theil-Sen Median趋势计算n(n-1)/2个数据组合的斜率中位数,计算公式如下

(1)

式中:β是植被变化趋势,当β>0时,表示EVI呈上升趋势,反之则为下降趋势。

Mann-Kendall显著性检验计算公式如下

设定,j=2001,2002,…,2020,定义Z统计量为

(2)

(3)

(4)

(5)

根据β、Z,并结合前人研究[24],本研究将显著性结果分为5个等级,分别为严重退化(β<-0.000 5,Z<-1.96)、轻微退化(β<-0.000 5,-1.960.000 5,-1.960.000 5,Z>1.96)。

1.3.2 偏最小二乘回归 偏最小二乘回归(PLSR)是一种集合了多元线性回归、典型相关分析和主成分分析3种研究方法的基本功能和优点、能够有效解决共线性问题,且适用于小容量样本的多元统计方法[25]。目前已在土地利用、景观格局、绿色空间变化等领域的研究中得到应用[26]。鉴于本研究数据容量较小,不适用于传统线性回归,因此选用PLSR定量分析气候与人类活动两大类因素(表1)与植被覆盖度变化之间的关系。

表1 2001-2020年长三角植被覆盖度变化驱动因子Table 1 The driving factors of vegetation cover in the Yangtze River Delta during 2001-2020

2 结果与分析

2.1 EVI的时空演变特征

2.1.1 EVI的时间变化特征 由图2可知,20 a来研究区植被覆盖整体上得到显著改善,EVI整体增长速率约为0.16%·a-1,最低值为2001年的0.433,最高值为2014年的0.474。

图2 2001-2020年长三角EVI年均变化Fig.2 Annual variation of enhanced vegetation index (EVI) in the Yangtze River Delta during 2001-2020

2.1.2 EVI的空间分布特征 根据长三角地区2001-2020年EVI多年平均空间分布格局(图3),研究区植被覆盖格局整体呈现“西高东低、南北高中部低”的分布特征。西部和南部多为安徽省和浙江省的山地区域,植被生长良好,EVI值较高,介于0～0.729,均值为0.486。北部地区以安徽省和江苏省的平原为主,作物种植广泛,EVI值也相对较高,介于0～0.694,均值为0.478。中部为3省的核心城市群,和东部的上海市,均为经济发展较为迅速、人工扰动较大的区域。因此EVI较低,介于0～0.679,均值为0.374。

图3 2001-2020年长三角EVI空间分布Fig.3 Spatial distribution of EVI in the Yangtze River Delta during 2001-2020

2.1.3 EVI的动态变化趋势 根据图4、表2来看,植被覆盖以稳定、改善为主,其中明显改善区域占比最大(43.26%),主要分布在江苏省北部、安徽省西部及南部,和浙江省绝大部分区域;轻微改善和稳定不变区域约占34%,主要包括安徽省西部、江苏省东部,以及浙江省各市外围;轻微退化和严重退化区域占比相当,均为11%左右,零星分布在各大城市的中心城区,并向外辐射。

表2 2001-2020年长三角EVI动态变化趋势分类表Table 2 Classification table of dynamic trend of EVI in the Yangtze River Delta during 2001-2020

图4 2001-2020年长三角EVI动态变化趋势Fig.4 Dynamic trend of EVI in the Yangtze River Delta during 2001-2020

2.2 EVI时空变化的驱动分析

2.2.1 驱动力主成分分析 借助SPSSAU平台进行最小二乘回归分析,根据交叉有效性(Qh2)、解释能力R2(表3)以及投影重要性(VIP)(表4),确定最佳主成分个数为1,因此选取1个主成分再次进行回归分析。

表3 偏最小二乘回归结果Table 3 Results of partial least squares regression

表4 不同主成分时投影重要性(VIP)指标汇总Table 4 Summary of variable projection importance (VIP) with different principal component number

提取的1个主成分(用A表示)与原始驱动因子的表达式为

A=0.165X1-0.113X2+0.059X3-0.065X4+0.318X5-0.314X6+0.359X7+0.319X8+0.354X9+0.374X10+0.383X11+0.329X12+0.044X13

(6)

由式(6)可知,主成分A主要受非农人口、生产总值、一二产增加值的正向影响,以及农业人口的负向影响,是一个充分反映城镇化发展水平的指标。

2.2.2 EVI变化驱动分析 经过计算,最终的偏最小二乘回归方程如下

EVI=0.041X1-0.028X2+0.015X3-0.016X4+0.080X5-0.079X6+0.090X7+0.080X8+0.089X9+0.094X10+0.096X11+0.093X12+0.011X13

(7)

投影重要性(VIP)代表了自变量对因变量的解释程度。由式(7)和表4可知,人口和社会经济因素VIP基本均>1,而自然因素VIP基本均<0.5,因此,人类活动因素对长三角地区植被覆盖度变化的影响程度大于气候因素。

2.2.2.1 气候因素对植被覆盖度的影响 气候因素中日照时数和相对湿度对植被覆盖度起负向影响(回归系数RC<0),年降水量和平均气温则起正向影响(RC>0),年降水量的作用相对较大(VIP>0.5),说明降水对植被的影响超过气温。

2.2.2.2 人类活动因素对植被覆盖度的影响 人类活动因素中二产增加值>一产增加值>非农人口>生产总值>三产增加值>人口密度>总人口,均起显著的正向影响(VIP>1,RC>0)。造林总面积也起着正向影响,但作用较小(VIP<0.5,RC>0)。而农业人口则起显著的负向影响(VIP>1,RC<0)。即20 a间长三角地区的快速城镇化进程成为了引发植被覆盖度变化的关键因素。

土地利用变化被认为是植被动态变化最重要的驱动力[27],而快速的城镇化进程极大地改变了土地利用类型,并导致相应的植被覆盖变化[28]。2001-2020年,长三角地区土地利用变化较大,变化面积达35 374 km2(表5),最显著的是耕地的转出与建设用地的转入。耕地转出面积为23 751.75 km2,主要转为建设用地(贡献率73.6%)、林地(贡献率13.5%)及水体(贡献率12.9%)。建设用地转入面积为19 271 km2,主要来自于耕地(贡献率90.7%)和水体(贡献率6.3%)。林地和水体也有较大变化,主要是大部分转出为耕地,但同时也由部分耕地转回。

由表6可知,耕地、林地、灌木、草地以及水体的EVI变化趋势均以改善为主(占比>50%),建设用地则两者兼有(退化46%、改善44%)。结合表7和图4、图5、图6可以发现,占用城镇周边原有耕地、林地等的建设用地,植被退化严重,而原有的城镇建设用地则多数呈现植被改善的结果。同时,植被在耕地、林地、水体之间的动态转化包括:1)退化,主要集中在退耕还湖的区域,如石臼湖、太湖周边等;2)改善,主要集中在退耕还林的区域,尤其是浙江省。

表6 2001-2020年长三角不同土地利用类型植被变化趋势占比Table 6 Percentage of vegetation trends in different unchanging land use types in the Yangtze River Delta during 2001-2020 %

表7 2001-2020年长三角不同土地利用变化类型植被变化趋势占比Table 7 Percentage of vegetation trends in different land use change types in the Yangtze River Delta during 2001-2020 %

1.耕地;2.林地;3.灌木;4.草地;5.水体;7.未利用地;8.建设用地,1→2即为耕地转向林地。

总的来说,尽管城镇化快速发展而导致生态用地被占,以及部分区域大力发展渔业等经济社会活动,对植被覆盖产生负向影响,但通过实行城市绿化建设以及生态保护工程,20 a来区域的快速城镇化进程总体上对植被覆盖起到了显著的正向影响。

与此同时,土地利用与植被覆盖的变化也导致了区域景观格局以及生态质量的变化[29],借助景观格局指数能够较好的量化这一过程。本研究选择的景观指数定义及计算公式见文献[30-31]。景观水平上(表8),2001-2020年长三角地区的斑块密度、周长面积分维数和香农多样性指数小幅上升,最大斑块指数、蔓延度指数明显下降,表明区域的景观破碎度加深,边界形状和景观结构更加复杂,同时景观异质性增强,景观类型呈多样化趋势。

表8 2001-2020年研究区景观水平景观指数Table 8 Landscape indices at landscape level from 2001 to 2020

从类型水平来看(表9),20 a来灌木、草地、未利用地的斑块密度较低,且变化较小,受人为干扰较少,因此重点讨论耕地、林地、水体以及建设用地的景观格局指数变化。4种地类的斑块密度、边缘密度、景观形状指数和周长面积分维数均呈上升趋势(除林地斑块密度小幅下降),表明4种地类的景观破碎度增加,景观形状更加复杂。最大斑块指数方面,水体和建设用地有所增加,表明两者的景观结构更加复杂,人为干扰程度更加严重;耕地和林地则有所下降,说明两者的景观优势度和完整性受到一定破坏。斑块聚集度指数方面,耕地、林地和水体小幅降低,而建设用地则显著提升,意味着前三者的连通度下降,建设用地的连通度明显提高,即其在扩张过程中表现为集聚化、连片化的发展趋势。

表9 2001-2020年不同土地利用类型的景观指数Table 9 Landscape indices of different land use types from 2001 to 2020

3 结论与讨论

3.1 结论

时序上,研究区EVI总体呈上升趋势,最低值0.433,最高值0.474,增长速率为0.16%·a-1,植被状况较好。空间尺度上,EVI总体呈现“西高东低、南北高中部低”的分布特征,高值区均值0.486,低值区均值0.374。

植被覆盖变化趋势以改善(76.29%)和稳定(8.0%)为主,主要分布在以林地为主的西部和南部。植被退化区域(22.71%)集中在城市扩张地区,并向外辐射。

人类活动因素(VIP>1)对研究区20 a来的植被覆盖变化影响程度大于气候因素(VIP<0.5)。气候因素中,年降水量和平均气温起正向影响,降水成为最主要的气候驱动因子。主导人类活动因素的快速城镇化进程引起土地利用的较大变化,但由于生态保护工程等措施,人类活动总体上对植被覆盖起到了显著的正向影响。

3.2 讨论

以往在运用EVI等植被指数对区域植被覆盖动态变化进行监测和研究时,经常会受到遥感图像质量和时序连续性的制约,而GEE云平台强大的云计算能力和海量的数据库能够较好的解决这些问题[32]。本研究通过PLS模型对长三角地区植被覆盖变化的驱动力进行研究,定量分析了气候和人类活动的影响力,但驱动因子是基于前人研究成果进行选取的,往往存在信息遗漏或冗余的问题。同时,因子之间往往也存在交互作用,普通的回归模型难以分析多个因变量之间复杂的相互作用,今后可以通过机器学习方法筛选驱动因子、构建结构方程模型等,从而更好地理解植被覆盖变化的复杂机制。

尽管城镇化快速发展总体上对植被覆盖起到了显著的正向影响,但通过景观格局指数的计算发现,区域的景观异质性有所增强,而景观破碎度增加,景观形状和结构更加复杂。耕地和林地的景观优势度和完整性受到一定破坏,建设用地则在快速扩张的同时成片聚集。因此改善区域植被覆盖度不仅需要恢复植被退化区域,更要注意保护耕地和林地,控制用地类型的转移,并通过构建生态安全格局、土地集约利用等方式减少景观破碎,提升连通性,从而提升生态质量。