成渝双城经济圈NPP 时空变化及气候因子驱动分析

梁帅, 辛宇, 辜寄蓉

四川师范大学地理与资源科学学院, 四川 成都 610068

政府间气候专门委员会(IPCC)第六次评估报告提出：自第一次工业革命以来，人类活动所排放的温室气体是导致全球气候变暖加剧的重要原因，如何实现碳中和和碳达峰以应对气候变化挑战成为近年来研究热点[1]。植被净初级生产力(NPP)是指植被在单位面积、单位时间内固定有机干物质的总量[2]，其作为陆地碳循环的重要组成部分，是判定生态系统碳循环的重要指标[3-4]。探究植被NPP 动态变化及其驱动因素，可以反映生态系统固碳能力的时空分异规律，揭示生态环境与气候变化之间的关系，为生态系统保护及“产碳端”和“固碳端”的动态平衡提供参考意见[5]。早在19 世纪80 年代就已经有学者开始研究植被NPP，但由于技术局限性，早期大多研究都停留在实地调查再构建回归模型计算NPP[6-8]，这种获取植被NPP 的传统方法不仅耗时、耗力且难以开展。随着遥感技术的发展，利用遥感数据与数学模型计算植被NPP 已经成为大尺度区域研究NPP的重要方法[9]。这种方法能获取高时密NPP 数据，避免了地面实测数据时间不连续、区域尺度小等缺点，并能够利用数据对其进行长时间序列监测及探索其驱动因子，为决策者制定气候政策达到碳达峰、碳中和等目标提供重要参考。近年来，许多学者利用MODIS 数据分析了不同时间、区域下的植被NPP 动态变化规律，研究发现近20 年间在我国植被NPP 在空间上呈现由西北向东南依次递增的格局，在时间上总体呈上升趋势[10-14]，此外研究表明NPP 的时空变化受气候、高程以及人类活动组合影响[15]，其中气候因子直接影响植被覆盖度变化对植被固碳趋势的贡献度即NPP 和GPP[16-18]。

成渝地区双城经济圈作为长江上游的重要城市群，其生态环境保护及可持续发展对整个长江流域的社会经济高质量发展有着举足轻重的作用。在“十四五”时期，成渝地区生态文明建设进入协同推进减污降碳、实现生态环境质量改善由量变到质变的重要阶段，监测成渝地区双城经济圈的植被NPP 时空变化及其气候因子影响，为长江上游生态文明建设及生态屏障构建提供决策性意见。现已有许多学者对川渝、西南地区的植被动态变化做了众多研究，例如茆杨、何宏昌等学者基于MODIS-NPP 数据对西南地区植被NPP 的时空变化及影响因素进行了研究[10,19]，郭斌、王耠熠等学者利用MODIS 及气象数据分析了近20 年来川西高原若尔盖地区的NPP 时空变化特征及气候响应[20-21]。基于MODIS-NPP 数据对2001—2020 年成渝地区双城经济圈植被NPP 时空变化特征进行研究，并利用ANUSPLIN 软件插值气象站点数据分析了植被NPP 对气候变化的响应，为长江上游城市群的经济发展与生态环境保护的动态平衡提供参考价值。

1 数据与方法

1.1 研究区概况

成渝地区双城经济圈，位于“一带一路”和长江经济带交汇处，以成都和重庆为中心，包含四川15 个市和重庆29 个区县（见图1），是西南地区经济水平最高，发展潜力最大的城镇化区域[22]。地形与四川盆地基本一致，中部多为平原，周围由低山、丘陵等环绕，气候适宜，自然资源丰富，生态禀赋优良[23]。

图1 研究区地理位置示意图Fig. 1 Schematic diagram of geographical location of the study area

1.2 数据来源

植被NPP 数据来源于美国NASA EOS/MODIS的2002—2021 年（实际计算植被NPP 年份为2001—2020 年）的MOD17A3HGF6.0 版本的全球NPP 数据，空间分辨率为500 m。基于GEE 平台对MODIS数据进行预处理，获取统一坐标的成渝地区双城经济圈20 年的植被NPP 数据。对预处理后的各年份NPP 数据乘以转换因子并剔除多年份NPP 平均值为0 的像元。

数字高程数据来源于地理空间数据云平台，分辨率为30 m。气象数据来自中国气象局数据信息中心的“中国地面气候资料日值数据集(V3.0)”，其数据时间尺度为2001 年1 月—2020 年12 月，包括日照时数、风速、降水、湿度、水面蒸发、日平均温度等6 个气象要素逐日数据。基于成渝地区双城经济圈内的气象站点日观测数据，采用ANUSPLIN 软件，以DME 为协变量插值生成空间分辨率为1 km的气象栅格数据，再经过重采样生成所用的气象数据，空间分辨率为500 m。对比其他插值方法，ANUSPLIN 插值精确度得到显著提高[15]。

1.3 研究方法

1.3.1 植被NPP 时空演变分析方法

采用一元线性回归模型来分析2001—2020 年成渝地区双城经济圈植被NPP 整体变化趋势及空间异质性。

式中， θ是NPP 年际变化趋势的斜率，n 为统计时段的累计年数，NPPi为第i 年植被npp 值，若 θ＞0 则表明植被NPP 呈上升趋势，若θ＜0 则表明植被NPP呈下降趋势，若θ=0 则表明植被NPP 几乎无变化。

采用F 检验法对植被NPP 的变化趋势进行显著性分析，其公式如下：

变异系数(Cv)表示为标准差和平均值的比值，是用来描述数据离散程度的指标[24]，采用Cv来反映近20 年植被NPP 的波动特征，其公式如下：

式中，NPPi为第i 年植被npp 值，是20 年内NPP 的平均值；Cv值越小，植被NPP 变化越小，反之亦然[25]。参考ArcGIS 软件中的自然分类法[26]，将Cv分为以下五个等级：稳定（Cv＜=0.10）,较稳定（0.10＜Cv＜=0.20），不稳定（0.20 ＜Cv＜=0.30），极不稳定（Cv＞0.30）。

1.3.2 相关性分析

基于像元偏相关分析方法计算植被NPP 与气温、降水影响因子的偏相关系数，其公式如下：

式中，Rxy,z表示为控制变量z 之后x 与y 的偏相关系数，x 表示NPP 值，y 表示降水量，z 表示气温，Rxy、Rxz、Ryz分别表示NPP 与降水量、NPP 与气温以及气温与降水量的线性相关系数。采用t 检验法对偏相关系数进行显著性检验。

复相关分析法是研究几个要素与某一要素的相关关系[27]。采用复相关分析法探究气温与降水对植被NPP 的共同影响作用，其公式如下：

式中，Rxy表示NPP 与降水量的线性相关系数，Rxz,y表示NPP 与气温的偏相关系数。采用F 检验法对复相关系数进行显著性检验。

1.3.3 驱动力分区准则

根据陈云浩等学者研究的植被变化驱动因子分区准则，将引起成渝地区双城经济圈NPP 变化的驱动力进行分区，根据偏相关系数与复相关系数的显著性检验结果，将驱动力分为下表五类（见表1）[27]。当复相关显著时：气温与NPP 偏相关的像元定义为气温驱动；降水与NPP 偏相关的像元定义为降水驱动，气温、降水分别与NPP 偏相关的像元定义为气温降水强驱动；气温、降水与NPP 偏相关不显著的像元定义为气温降水共同驱动。不满足复相关显著的像元定义为非气候因子驱动。

表1 NPP 变化驱动因素分区准则Tab. 1 Zoning criteria for NPP change drivers

2 结果与分析

2.1 植被NPP 时空演化分析结果

2.1.1 植被NPP 空间分布特征

研究区内2001—2020 年植被NPP 均值分布如图2所示，成渝地区双城经济圈植被NPP 平均值为563.14 g C·m-2·a-1，年均值分布在＜400 g C·m-2·a-1、400—500 g C·m-2·a-1、500—600 g C·m-2·a-1、600—700 g C·m-2·a-1、 700 —800 g C·m-2·a-1、 ＞800 g C·m-2·a-1的区域分别占总面积的5.49%、28.79%、34.31%、14.33%、8.95%、8.13%。从空间分布来看，研究区内植被NPP 呈现四周高中间低的分布格局，其与四川盆地地形特征有着密切联系。结合不同海拔下植被NPP 均值分布图来看，随着海拔的升高，植被NPP 均值也在增加，在海拔＞3 000 m 的情况下，均值开始显著降低，但通过整体趋势线可以看出植被NPP 呈增加趋势。在德阳、成都和雅安一带海拔较高的山脉植被NPP 均值高，在海拔较低的盆地内部的各城市中心植被NPP 均值最低，在重庆黔江区、彭水县等低山丘陵处植被NPP 均值较高。

图2 成渝地区双城经济圈2001—2020 年植被NPP 均值分布图及不同海拔下植被NPP 均值分布图Fig. 2 Mean NPP distribution of vegetation in Chengdu-Chongqing economic circle in 2001-2020 and mean NPP distribution of vegetation at different altitudes

2.1.2 植被NPP 时空变化分析

2.1.2.1 植被NPP 年际变化趋势分析

利用一元线性回归法及变异系数分析研究区植被NPP 的时空变化特征（见图3a 和3b），并对变化趋势做显著性检验，将结果分为：显著减少、减少、稳定、增加、显著增加五类（见表2 和图3c）。

表2 植被NPP 变化趋势显著性结果统计Tab. 2 Statistical results of significance of vegetation NPP change trend

图3 成渝地区双城经济圈2001—2020 植被NPP 变化趋势Fig. 3 Variation trend of vegetation NPP in Chengdu-Chongqing economic circle from 2001 to 2020

2001—2020 年植被NPP 变异系数如图3a 所示，可以看出成渝地区双城经济圈近20 年来植被NPP 较为稳定，大部分区域变异较弱，变异较强的地方存在于成都市和重庆市主城区周边以及沿长江一带。从斜率分布图(见图3b)及近20 年植被NPP 变化分类图（见图3c）来看，其变异较强的区域都呈现显著减少趋势；在植被NPP 均值较高的德阳、成都和雅安一带的山脉以及重庆市的黔江区、彭水县等低山丘陵处变化甚微，趋于稳定；在城市群内部盆地及丘陵地带植被NPP 呈现增加趋势；结合表2 统计结果来看，其植被NPP 呈现显著增加趋势的面积最大，占研究区总面积的70.30%；呈现减少趋势的面积最小，其面积占比为0.43%。从2001—2020 年植被NPP 年均值变化（见图3b 和3d）来看，整体呈现增加趋势，其中植被NPP 呈上升趋势的面积占比达95.08%，呈下降趋势的面积占比为4.92%；2005 年研究区植被NPP 年均值最小为471.93 g C·m-2·a-1，植被NPP 年均值最大出现在2020 年为644.99 g C·m-2·a-1。

2.1.2.2 不同海拔下植被NPP 动态变化趋势分析

将成渝地区双城经济圈海拔划分成不同梯度等级，统计各海拔梯度下植被NPP 的变化情况，从图4可以看出，在研究区内，随着海拔的增加，植被NPP的年变化趋势率不断下降，在海拔小于500 m 的区域植被NPP 增长趋势最显著，年变化趋势率达到8.72 g C·m-2·a-1；在海拔大于3 000 m 的区域植被NPP 最稳定，年变化趋势率为0.95 g C·m-2·a-1；在海拔500—1 000 m、1 000—2 000 m、2 000—2 500 m及2 500—3 000 m 高程区间内，年变化趋势率分别为7.24 g C·m-2·a-1、4.62 g C·m-2·a-1、2.99 g C·m-2·a-1、1.94 g C·m-2·a-1。从整体来看，成渝地区双城经济圈植被NPP 整体处于增长趋势，在海拔大于3 000 m的区域人类活动少、植被覆盖度高其植被NPP 均值较高且变化小，在海拔小于500 m 的区域，由于城市扩张、人类活动干扰其植被NPP 均值相对较低，近年来天然林保护、长江流域生态屏障构建等措施致使植被NPP 均值得到显著增加。

图4 成渝地区双城经济圈各海拔梯度植被NPP 年际变化图Fig. 4 Interannual variation of vegetation NPP at different elevations in Chengdu-Chongqing economic circle

2.2 植被NPP 对气候因子的响应

2.2.1 植被NPP 与气温、降水的相关性分析

通过一元线性趋势分析法，可知研究区20 年里气温和降水都呈现增长趋势，年变化趋势率分别为0.0067℃a-1、8.14 mma-1。

研究发现成渝地区双城经济圈植被NPP 与气温、降水的偏相关系数范围都在-0.83～0.87 之间，且相关性空间格局存在异质性（见图5）。分析表明研究区植被NPP 与气温呈显著负相关、负相关的区域面积占比分别为0.39%、25.96%，呈显著正相关、正相关的区域面积占比分别为6.22%、67.43%；与降水呈显著负相关、负相关的区域面积占比分别为0.21%、7.33%，呈显著正相关、正相关的区域面积占比分别为23.40%、69.06%。

图5 成渝地区双城经济圈植被NPP 与气温、降水的偏相关系数（a、b）及显著性（c、d）Fig. 5 Partial correlation coefficients (a, b) and significance (c, d) of vegetation NPP with air temperature and precipitation in Chengdu-Chongqing economic circle

从空间上来看，研究区植被NPP 与气温正负关系并存，在其西南部即眉山市、乐山市和雅安市等地植被NPP 与气温呈负相关关系，在重庆市黔江区及彭水县部分地区存在负相关关系且零星出现显著负相关关系；在研究区中部及东北部植被NPP 与气温正相关关系明显，尤其是在重庆城市群及东部城市群的部分地区存在显著正相关关系。研究区大部分区域植被NPP 与降水之间都呈现正相关关系，并且在成都城市群与东部城市群的区域正相关关系显著；只有在东部城市群零星存在负相关关系。其结果表明了气温在不同区域对植被NPP 变化的影响不同，而降水在研究区内几乎都起着积极影响。

2.2.2 植被NPP 变化气候因子驱动分析

植被NPP 的变化除了受气温、降水等单方面因子的作用，同时受气温与降水等因子的共同驱动作用，通过植被NPP 与气候因子的相关性及其显著性分析，将驱动因子进行分区，从图6 可知，植被NPP与气温、降水驱动因子的复相关系数范围为0.002～0.867，其高值区主要分布在东部城市群及成都市与德阳市的西部，在东部城市群的东北部、雅安市与宜宾市的屏山区系数值较低。研究发现成渝地区双城经济圈植被NPP 变化主要受非气候因子驱动，其面积占比达86.87%；受气候因子影响的区域面积占研究区总面积的13.13%，其中受气温驱动、降水驱动、气温降水共同驱动、气温降水强驱动的面积占比分别为1.58%、8.40%、0.16%、2.99%。从空间上来看，受气温驱动的区域主要分布在东部城市群的平行岭谷地带及重庆城市群的江津区、綦江区；受降水驱动的区域分布最广，包括了成都城市群的西北部及南部部分区域、南部城市群的泸州市部分区域、重庆市黔江区部分区域；受气温、降水共同驱动的区域十分有限，分布稀少；在南充仪陇、蓬安县一带受气温、降水强驱动的区域最为密集。

图6 复相关系数及驱动力分区图Fig. 6 Partition diagram of complex correlation coefficient and driving forces

3 结论与讨论

基于MOD17A 产品NPP 数据，利用一元线性趋势法分析2001—2020 年成渝地区双城经济圈植被NPP 变化特征，并结合气象数据，运用相关分析法，研究植被NPP 变化的气候因子驱动力。

研究结果表明成渝地区双城经济圈植被近20 年间NPP 在471.93—644.99 g C·m-2·a-1范围之间，均值为563.14g C·m-2·a-1，比全国平均值514.48 g C·m-2·a-1高48.66g C·m-2·a-1，其总体呈现增长趋势，与冯婉[12]、张凤英[28]等学者的关于成渝城市群地区的研究结果一致。从年际变化上来看，研究区内植被NPP 整体呈现缓慢增长的趋势，增长率为7.53 g C·m-2·a-1。在植被NPP 均值较高的德阳、成都和雅安一带的山脉以及重庆市的黔江区、彭水县等低山丘陵处趋于稳定；在城市群内部盆地及丘陵地带植被NPP 呈现增加趋势，这一现象得力于“退耕还林、还草政策”、“长防林工程”、“长江保护修复攻坚战行动计划”[29]。从空间上来看，研究区内植被NPP 呈现四周高中间低的分布格局，随着海拔的升高，植被NPP 值不断增加。海拔在2 500—3 000 m 的高程区间内植被NPP均值为805.21 g C·m-2·a-1，海拔小于500 m 的高程区间内植被NPP 均值为491.06 g C·m-2·a-1。这与四川盆地地形息息相关，这可能是因为研究区靠近川西高原的部分高海拔区域，由于气候影响，生态环境脆弱，植被覆盖率低，NPP 低且较为稳定，在盆地城市内部由于地势平坦，水热条件好，其城市不断扩张，NPP 值较低，在低山丘陵地区，因气候适宜而少受人类活动影响，植被覆盖广，NPP 值最高，该结果与胥晓[30]、张锐[31]等学者的研究结论相符合。

此外，基于气象站点数据，利用ANUSPLIN 软件对气象数据进行插值，分析气候因子对植被NPP变化的影响。结果表明研究区植被NPP 与气温、降水的相关性存在空间异质性，气温在眉山市、乐山市和雅安市等地以及重庆市黔江区、彭水县部分地区对植被NPP 起着抑制作用，在重庆城市群及东部城市群的部分地区促进了植被NPP 的增长，而降水与植被NPP 在整个研究区内都呈正相关的关系，这可能是因为水热条件决定了植被的生长，而降水的升高使植被获得更多的水分，促进植被生长，使得NPP 不断增加，但研究表明气温对植被NPP 的影响是个复杂的过程，气温升高能够为植被提供更多的热量促进植被光合作用从而导致NPP 增加，但过高的温带增长会对植被产生高温胁迫[10,32]。此外，研究表明研究区内受到非气候因子影响的区域面积占比达到86.87%，植被NPP 分异较大程度受到了人为因素干扰，但仅针对气候因子来说，气温、降水仍然是影响植被NPP 变化的主要因素。

以往大多数研究多选择长江流域、川西高原或横断山区这样的生态区内[4,12,33]，但对成渝这样的城市群体研究甚少，成渝城市群作为长江上游最大的城市群体，其植被NPP 变化趋势对长江上游生态屏障构建及生态修复有着举足轻重的影响。