基于MGWR的长江流域植被演化及其影响因素

李泳君,陈青长*,方 贺,李 建(.上海应用技术大学城市建设与安全工程学院,上海 048；.浙江省气候中心,浙江 杭州 30057)

植被作为陆地生态系统的重要组成部分之一[1],在维持和优化生态系统服务等过程中发挥着关键作用,并被视为全球环境变化的敏感指标[2].因此,进行区域和全球尺度的植被变化研究对于深入理解生态系统的物质循环和信息传递具有重要意义[3].

在受到气候变化、人类活动和地形差异等多种因素的干扰下,植被往往会呈现出复杂的空间异质性演变[4-5].从全球尺度来看,气候变暖加速了北方高纬度地区植被的恢复速度[6],但对热带地区的植被却产生了负面影响[7].在流域尺度上,雅鲁藏布江流域上游植被与降水、干旱呈显着正相关,中游植被与气温、下游植被与降水呈显着负相关[8].总之,不同影响因素对植被的影响关系会随着空间位置的变化而发生变化(即空间非平稳性).因此,揭示植被活动与影响因子之间的空间非平稳关系,能够为应对气候变化、城市扩张以及提升生态系统适应能力提供理论依据[9-10].

长江流域作为亚洲第一大流域[11],其植被覆盖状况及动态变化对中国、周边国家乃至全球生态平衡具有重要影响[11].因此,研究长江流域植被的时空变化,以及驱动这些变化的主要影响因素,对于区域生态环境保护和植被资源管理具有重要现实意义.近年来,已有许多学者对长江流域的植被变化及其影响因素进行了研究,他们利用不同分析方法探讨了植被的驱动机制[11-14].这些研究从不同的角度对长江流域的植被变化及其驱动机制进行了深入探讨,但也存在以下局限性:1)影响因素考虑不充分,以往研究侧重于研究气候地形等自然因素对植被的影响.2)影响因素作用尺度差异,大多数学者假定各影响因子对于植被的驱动过程均是在相同尺度下.3)驱动分析方法的不足,相关性分析、偏相关分析或者地理探测器没有顾及植被与相关影响因子的空间非平稳性.

自1996年Brunsdon 等[15]提出了空间非平稳性分析工具--地理加权回归模型(GWR)以来,该模型已成为应用广泛的局部空间统计分析方法之一[16-17].尽管许多学者已经应用GWR 模型对植被与外部环境因素的空间非平稳性进行了深入研究[1,18-19],然而经典的GWR 假定各因子对于植被的驱动过程均是在相同尺度下.这个假设明显违背现实地理空间关系的异质性或非平稳性规律,不可避免地忽略了不同关系可能在不同尺度上发生的事实[20].为解决多尺度的空间非平稳性和异质性问题,Fotheringham 等[21]提出了顾及多尺度的地理加权回归模型(MGWR).该模型能够很好的模拟自变量与因变量的地理空间响应过程,从而精确识别不同影响因子对于因变量作用的不同尺度,程度及方向[22].

综上所述,针对以往研究对于长江流域植被变化与影响因子尺度效应考虑不全面,以及各因素与植被的空间关系的非平稳性研究不足等问题.本文借助GEE 云平台探讨了2000～2022年长江流域植被覆盖度的时空演化特征,应用MGWR 等方法,研究流域植被与不同类型要素(地形、气候、社会等方面)的空间非平稳关系,以定量分析不同影响要素对植被活动的控制作用.

1 数据与方法

1.1 研究区概况

长江是中国第一大河,也是世界第三大河,全长约6300km,流经中国青海、西藏、云南、四川、重庆、湖北、湖南、江西、安徽、江苏等十多个省市(图1).长江流域地处中国中部地区,90°～122°E,24°～35°N 之间,总面积约180 万km2,占中国陆地总面积的1/5.长江流域的气候多样,从温带到亚热带,不同地区的气候、水文特征各异.由于山地、丘陵、平原交错,地势逐渐平缓,长江流域的地貌变化也相当复杂,形成了长江三角洲、洞庭湖平原、四川盆地等地貌类型.长江流域是中国重要的农业和工业基地,拥有众多的水力发电站和水运枢纽,是中国最重要的内河航运干线.同时,长江流域也是中国珍稀濒危野生动植物的重要栖息地,保护生态环境对于长江流域的可持续发展至关重要.总的来说,长江流域是中国最为重要、多样化和广泛的流域之一,对中国的经济、文化和自然资源都有着重要的影响.

1.2 数据源

本研究数据主要分为两大类,分别为因变量数据(植被覆盖度 FVC),是利用 GEE 平台(https://code.earthengine.google.com/)进行计算得到,数据源来自美国国家航空航天局(NASA)官网(https://www.nasa.gov/)的 MODIS 植被指数产品数据MOD13Q1.数据空间分辨率为250m,时间分辨率为16d,时间范围:2000年2月～2022年12月.自变量(影响因子)包括地形因子、气候因子、社会相关因子(表1),各因子的空间分布如图2所示.

表1 影响因子数据Table 1 List of impact factor data

图2 影响因子的空间分布Fig.2 Spatial distribution of Impact factor

1.3 植被覆盖度反演方法

植被覆盖度(FVC)采用像元二分模型进行计算,该模型假设像元的地表由有植被覆盖部分地表与无植被覆盖部分构成,而传感器观测到的光谱信息也由这2 个组分因子线性加权合成,各因子的权重是各自的面积在像元中所占的比率,如其中FVC 可以看作是植被的权重[23-24].

计算FVC 的公式如下:

式中:NDVIsoil为纯裸地像元的NDVI 值;NDVIveg为纯植被像元的 NDVI 值[25].其中,裸地像元值NDVIsoil和纯植被覆盖像元值NDVIveg的理论值应分别接近0 和1.本研究中采用5%置信度截取NDVI上下阈值,将NDVI 数值最大、最小的5%区域分别做平均值[24],得到NDVIsoil和NDVIveg.

1.4 时空演化分析方法

1.4.1 变化趋势分析 本文使用 Sen 氏趋势分析[8]与MK 检验[26]相结合,分析FVC 在像元尺度上的变化趋势及显著性;在区域尺度上则采用线性回归分析FVC 变化的总体趋势.计算公式如下:

式中:FVCi和FVCj分别表示像元i年和j年的FVC值;n 表示时间序列的长度;sgn 是符号函数;统计量Z的取值范围为(-∞,+∞).在给定显著性水平α 下,当|Z|>u1-α/2时,表示研究序列在α 水平上存在显著的变化[27].

1.4.2 变异系数法 变异系数是反映观测数据的相对波动程度,具体衡量数据的稳定性.当变异系数值较大时,说明数据波动较明显,表示植被变化较大;反之,变异系数值较小表示数据波动较平缓,植被变化较为稳定[28-29].变异系数的计算公式如下:

1.4.3 Hurst 指数 本研究使用Hurst 指数来分析未来FVC 的演变趋势.该方法最初由英国水文学家Hurst于1965年提出[30-31],现已被广泛应用于各个领域的时间序列分析,该方法的计算公式详见文献[4-5].当0.5

1.5 驱动分析方法

地理加权回归模型可有效地解释空间非平稳性,其通过借用附近样本点的数据来进行参数估计,然而GWR 模型中的带宽是所有自变量不同带宽的平均值,进而导致估计偏差[21,33].作为GWR 模型的最新改进版本,多尺度地理加权回归(MGWR)模型同时考虑多种带宽以解决该缺陷,进而识别不同变量影响的尺度效应[33-34].因此,本研究拟使用MGWR模型以研究长江流域植被空间分布与其他因子之间的空间非平稳性与尺度效应,其公式如下:

式中:yi、xij、εi分别代表空间上i 点的因变量、自变量和随机误差;(ui、vi)为i 点的空间位置;j 代表自变量的个数;β是i 点上的回归系数[35].bwj 代表了第j 个变量回归系数使用的带宽.MGWR 的每个回归系数βbwj都是基于局部回归得到的,且带宽具备特异性.MGWR 的核函数和带宽选择准则依然延续使用经典GWR 的几种经典的核函数和带宽选择准则.本文使用最为常用的二次核函数和AICc 准则[36].

2 结果与分析

2.1 长江流域植被时空演化特征

以植被覆盖度(FVC)为植被指标,对长江流域植被的时空演化特征进行综合评估和定量研究.首先是对长江流域植被的空间演化特征进行定量分析,在此基础上,进一步对长江流域植被在时间维度上的变化特征进行深度分析.

2.1.1 长江流域植被空间分布 本研究基于GEE 云平台,计算得到2000～2022年长江流域平均年度FVC 的空间分布(图3),并参照已有研究将FVC 划分为5 个等级[6,29].由图3 可知,长江流域FVC 总体呈现出西北低东南高的空间分布趋势.其中,流域上游的青海东部以及四川的西北部FVC以<0.2 为主(图3a),面积占比约为7.04%.长江流域中部地区以0.4～0.8 的植被覆盖度为主,由图3b 可知,陕西南部,湖北西部、四川和重庆东北部植被覆盖相对较高(0.6～0.8),面积占比约为32.35%.长江中游城市群植被覆盖度相对较低(图3c),尤其以武汉、长沙、南昌等城市的植被覆盖较低(<0.2).此外,流域下游的上海、江苏南部、安徽东南部植被覆盖也以<0.2 为主.然而,整体来看,长江流域植被相对较好,FVC>0.4 的区域面积占比约为79.56%.

图3 长江流域植被覆盖度的空间分布Fig.3 Spatial distribution of FVC in the Yangtze River Basin

2.1.2 长江流域植被时间变化 为研究FVC 在2000～2022年间的变化趋势,本文将逐年 FVC 数据作为数据源,利用Matlab 编程实现Theil Sen-MK 趋势分析,得到植被变化趋势结果.将其划分为5 个等级,详见表2.将其划分为五个等级.

由图4 可知,长江流域植被覆盖度得到了显著改善,显著改善面积占流域总面积的56.56%.然而,严重退化区域仅占长江流域总面积的2.35%,主要分布在流域上游的青海东部(图4a)、中游长株潭城市群(图4c)以及下游长三角城市群(图4d).由图4b可知,流域中部是显著改善区域聚集区,包括甘肃南部、四川东部、陕西南部和湖北西部等地区.从植被覆盖度的年际变化情况来看(图5),2000～2022年长江流域植被覆盖呈现波动上升趋势,增速为0.245%/a.其中,FVC 最大值为2021年的53%,最小值为2005年的45%.其中,2012年前长江流域植被覆盖度增速较慢,2012年后植被覆盖度增速显著,增速为0.636%/a.但总体上看,近23a 长江流域植被覆盖度得到显著改善.

图4 长江流域植被覆盖度的趋势变化空间分布Fig.4 Spatial distribution of trend changes of FVC in the Yangtze River Basin

图5 长江流域年际植被覆盖度变化曲线Fig.5 Annual changes of FVC in the Yangtze River Basin

2.1.3 长江流域植被变化稳定性 根据长江流域2000～2022年植被覆盖度结果,逐像元统计变异系数值,依照计算结果将稳定性分为五个程度,低等波动变化(<0.05),中低等波动变化(0.05～0.10),中等波动变化(0.10～0.15),中高等波动变化(0.15～0.20),高等波动变化(>0.20).如图6所示,2000～2022年间长江流域植被变化在大部分地区的波动较平稳.其中,中低波动变化区域占比最高,为64.10%,广泛分布于流域的各个区域.中等波动变化区域面积占比为18.43%,主要分布在流域上游的四川东部,重庆西部以及贵州北部等地区(图6b).低等波动变化区域面积总占比为12.93%,主要位于四川西部(图6a)、湖北中西部(图6c)以及安徽南部(图6d)等地区.植被变化不稳定地区(>0.20)主要位于大中型城市,如上海,苏州,武汉、南京等.

2.1.4 长江流域植被未来趋势变化 为进一步了解长江流域FVC 未来变化趋势,将FVC 变化趋势与Hurst 指数进行耦合分析(表3),得到FVC 变化趋势与持续性组合的空间分布(图7).总体上看,长江流域植被覆盖度的未来发展趋势以由改善向退化为主,其面积占比为69.24%,广泛分布于流域各个区域.然而,FVC 未来由退化向改善发展的区域面积约占流域总面积的6.98%,主要分布于上游的青海东部(图7a),云贵川渝地区也有零星分布(图7b).持续退化区域面积占比为3.65%,主要位于流域中游的湖北中东部(图7c),以及流域下游的长三角城市群(图7d).持续改善区域占比相对较多,为17.10%,主要位于湖北东北部、河南南部,此外流域内其他省份也都有较为稀疏分布.

表3 植被未来变化趋势分级Table 3 Classification of future changes trend of vegetation

图7 长江流域植被未来变化趋势的空间分布Fig.7 Spatial distribution of future changes trend of vegetation in the Yangtze River Basin

2.2 植被变化影响因素的空间分异特征

研究植被动态及其决定因素对于进一步深入理解区域和全球尺度的气候-植被相互作用/反馈、指导可持续管理和自然资源保护至关重要.因此,在针对长江流域植被时空演化特征定量分析的基础上,本节对于长江流域植被演化的背后驱动机制进行探讨分析.首先,分析地形、气候及社会三种类型因素对于长江流域植被变化的驱动力.其次,讨论各类型因子对于长江流域植被变化驱动作用的尺度效应.

2.2.1 影响因素的空间异质性 本研究借助MGWR 模型在不同空间位置(县级)实现了驱动因子与FVC 关系空间异质性的特点,得到由模型输出的各驱动因子回归系数的空间分布(图8).其中回归系数的绝对值表明驱动因子对于植被的影响程度,回归系数的正负表明影响的方向.

图8 长江流域植被影响因素的回归系数Fig.8 Regression coefficients of factors influencing vegetation in the Yangtze River Basin

地形方面,高程对于长江流域中下游的江西、安徽、江苏、浙江及上海等地区植被具有较强的负向作用,回归系数绝对值聚集于-23.33～-4.84(图8a);而流域上游的青海、四川、重庆等地区主要受到了高程的较小的正向作用.坡度对于长江流域植被的作用呈现出自西向东逐渐降低的趋势(图8b),与高程相似,流域的中下游同样受到了坡度较大的负向作用;但是流域植被受到正向作用的地区面积占比明显高于负作用.

气象方面,降水对于长江流域植被以负向作用为主(图8c),主要位于流域中上游的青海,四川、重庆及湖南等地区;然而成都及其周边各县却受到较小的正向作用.与降水不同,气温对于植被的影响程度表现出西正东负的趋势(图8d),其中长三角地区植被受到气温的负向作用相对较强(-15.52～-8.96),然而上游地区受到的正向作用相对较低(0.00～1.67).相对湿度对于长江流域植被的影响以负向作用为主(图8e),只有湖北北部、陕西和河南南部等地区受到正向作用.蒸散发对于长江流域植被的影响出现显著的南北空间分异性(图8f),其中南部的云贵川渝以及湘鄂赣地区受到负向作用,而北部的青海、四川、陕西、河南以及长三角地区主要受到正向作用.如图8g所示,长江流域植被主要受日照时间的负向作用,其中流域上游大部分地区主要受到较小的负向作用(-1.19～0.00),只有成都及附近等区县受到较小的正向作用(0.00～1.13).

社会方面,GDP 对于长江流域植被的影响为单一负向作用,且由东北向西南呈下降趋势(图8h),影响程度也相对较小(-0.128～-0.116).与GDP 不同,人口密度对于长江流域植被的影响为单一正向作用,呈现出自西向东逐渐升高的趋势(图8i),影响力也相对较小(0.08～0.15).由图8j 可知,长江流域植被受到夜间灯光数据的单一负向作用,且影响程度相对较低.

2.2.2 影响因素的尺度差异 图9 总结了GWR 和MGWR 之间所有变量带宽的差异.MGWR 能够直接反映不同变量对于长江流域植被的差异化作用尺度,而传统GWR 则只能反映出各个变量作用尺度的平均值.由图9 可知,经典GWR 的带宽为66,占总样本数量的7.4%.然而,MGWR 模型的计算表明不同影响因素对于长江流域植被的作用尺度差异较大.其中,地形(高程和坡向)、气候(降水、气温、相对湿度、蒸散发及日照时间)因子对于长江流域植被的作用尺度相对较低,为43,仅占样本数的4.8.然而,社会因子对于长江流域植被的作用尺度相对较大,其中人口密度的作用尺度为886,为全局尺度.GDP的作用尺度仅次于人口密度,为872,而夜间灯光数据对于植被的作用尺度为237.

图9 经典地理加权回归与多尺度地理加权回归模型带宽Fig.9 Bandwidth of GWR and MGWR

3 讨论

3.1 长江流域植被时空变化特征分析

本研究表明,近22a 来长江流域的植被时空演化呈现出显著的空间异质性.从空间分布来看,长江流域植被呈现出东西低中部高的趋势,这与以往研究结果一致[12,37].从时间变化的角度来看,2000～2022年间,长江流域的FVC 总体呈现出波动上升的趋势.图5 也从空间分布上证实了这一趋势,这也与徐勇等学者的研究结论一致[11].长江流域植被改善区域明显高于退化区域,占流域总面积的82.36%,主要位于流域中部,这里是生态修复工程实施的区域[38],主要包括防护林工程、岩溶地区石漠化综合治理工程、天然林保护工程等[11].以往研究也得出相同结论[11,38-39].

3.2 长江流域植被变化驱动分析

3.2.1 地形因素对于植被的驱动分析 研究发现,长江流域的植被受地形因素的影响较大(图8 和表4),其中坡度对流域植被的影响明显大于高程.坡度作为植被变化的关键因素,通过调节植物所需的水热条件和营养有机质,间接影响着植被的生长变化[24].此外,海拔也会通过影响人类活动,进而间接影响植被变化[16].图8 显示,长三角城市群等低海拔地区的植被受到高程和坡度的负面影响较大,主要原因是海拔较低、坡度适中,适合人类活动和社会发展[40],同时也促进了农业开垦[41],进而导致植被覆盖相对较低[29].邓元杰等[13]的研究结果也支持了这一观点.

表4 多尺度地理加权回归系数统计描述Table 4 Statistical description of MGWR coefficient

3.2.2 气候因素对于植被的驱动分析 研究结果显示,气温对长江流域植被的影响最大,回归系数绝对值均值为-2.674(表4),Qu 等[42]学者的研究也得出类似的结论.图8d 表明,流域上游地区的植被覆盖度随着温度升高逐渐增加,这可能是由于气温升高增强了高海拔地区森林植被的光合作用[1,16].此外,青藏高原地区的温度上升缩短了冻融时间,改善了土壤水分,增加了植被的生长期,这也有利于植被繁盛[43-44].相对湿度对植被的影响次于气温,主要表现为负向作用,这也与之前的研究一致[45].降水对植被的影响也主要是负向的,但在中下游地区出现了部分正向作用,可能是由于三峡水电站开发和水库建设[11]导致下游水资源相对较低,增加的降水量促进了下游植被的生长.

3.2.3 人为因素对于植被的驱动分析 研究发现,人为因素和自然因素对植被NDVI 的影响方式存在显着差异(图8).与以往研究认为的人为因子对植被的较强影响相比,本研究发现人为因素对植被的影响程度明显较低,自然因素的作用更为显著.在影响方向上,社会经济因素对长江流域的植被呈现出正负空间非平稳性.例如,GDP 和夜间灯光数据对植被的影响呈现负向关系(图8h 和图8j),而人口密度与植被的关系则呈正向关系.这可能由于GDP 和夜间灯光数据的空间趋势与城市化发展的空间格局相关,以往研究认为GDP 和夜间灯光数据也具有相当的空间分布模式[46].随着经济水平和城市化程度的提高,城市建设不断扩张,持续占用植被所需的土地资源,从而导致植被的减少.

3.2.4 不同影响因素的尺度效应 从MGWR 模型导出的尺度系数可用于识别不同因子对于植被的不同作用区域.在使用GWR 模型时,忽略了部分变量的全局鲁棒性[20],所有因子对植被的作用尺度仅仅是单一的均值尺度(图9).不同因子的作用尺度差异性往往被忽视[16],这导致了空间非平稳性不显着的误导性结论[21].与GWR 不同,MGWR 的研究结果表明,流域不同因子对于植被的作用尺度显著不同.植被受到自然因子(地形因子和气候因子)作用的尺度相对较低,尺度为43.然而人为因子对于植被的作用尺度相对较大,如GDP 的尺度为872.分析结果表明,与人为因素相比,自然因子对植被的影响非常局部(带宽较小),这种对空间尺度关系的精确诊断支持具有成本效益的土地覆被规划[16,47],从这个角度制定植被恢复政策应侧重于区域间特征.

4 结论

4.1 长江流域植被覆盖度时空演化差异显著.空间分布上,植被覆盖度总体呈现出东西低,中部高的空间分布格局.植被覆盖度在0.4～1.0 之间的面积占比约80%,植被呈较高覆盖状态.时间变化上,近22a 长江流域植被整体呈现出波动上升的变化趋势,增速为0.245%/a.植被改善区域(82.36%)明显高于退化区域(7.12%),在空间上表现为中部改善强于东部,东部改善强于西部.未来发展来看,长江流域有退化风险,其中由改善向退化面积占比最大,为69.24%.

4.2 植被和各驱动因素之间的响应过程表现出显著的空间非平稳性.值得注意的是,与人为因子相比,地形和气候因子对于长江流域植被的影响更为显著,其中坡度、高程、气温及相对湿度回归系数均值都超过1.100.这些结果进一步证实了MGWR 模型可以应用于处理人类活动和环境变化等非平稳特征的类似研究,以同时捕捉时间和空间上的非平稳性.

4.3 各因子对于植被的作用尺度差异明显.其中植被对自然因子非常敏感,且存在很强的空间非平稳性,地形和气候的影响尺度是变量中最小的,为43.影响因子按其空间尺度由小到大分别为高程、坡度、降水、气温、蒸散发、日照时间、相对湿度、夜间灯光数据、GDP、人口密度.

致谢：部分数据和技术支持由中国气象局“风云卫星遥感产品校验”青年创新团队提供,在此表示感谢.