以温室气体CO₂浓度增长为主要驱动的全球变暖现象, 不仅对植物物候变化产生重要影响, 而且可能导致陆地植被分布及生态质量发生明显变化^[1-2]。植被总初级生产力(Gross primary productivity, GPP)是指单位时间内绿色植物通过光合作用途径所固定的有机碳总量, 表征了进入陆地生态系统的初始物质和能量^[3]。作为陆地生态系统吸收CO₂的重要指标, GPP在全球碳循环研究中具有极其重要的作用^[4-6]。科学研究地表植被GPP时空变化规律及其对气候变化的响应, 对深入理解陆地生态系统与大气之间的碳交换机理、评价陆地生态系统的环境质量以及预估未来气候变化的生态效应等具有重大意义^[7-8]。近年来, 随着对全球变化以及对碳循环方面的深入研究, 以GPP和NPP(Net primary productivity)为代表的陆地植被生产力及其对气候的响应以及与自然因素的关系逐渐成为研究重点。特别是在遥感技术的快速发展和支持下, 越来越多的遥感数据被用于开展大尺度区域性植被生产力时空变化研究。MODIS平台的GPP/NPP数据产品因其时空连续性好、数据精度高和获取成本低等优点, 在陆地生态系统碳循环研究中被广泛应用^[9-16]。众多研究结果表明, 气候变化通过改变植物生长发育进程中光热和水分的匹配而影响生态系统生产力^[17-18], 气候因子对植被生产力的影响存在显著的时空差异^[19-20], 而不同植被类型生产力对气候因子的响应作用也不一致^[21]。

长江流域是我国重要的工农业生产区和生态安全屏障, 横贯东、中、西三大经济区。作为一个独特而完整的自然生态系统, 长江流域具有强大的水土保持、生物育种、释氧固碳、环境净化等功能, 对维护周边地区乃至全国的生态平衡和安全具有重要作用。关于长江流域的植被及其生产力动态变化, 不少学者已开展过相关研究。如柯金虎等^[22]基于生态过程模型-CASA模型研究分析了长江流域植被NPP的动态变化。贾松伟^[23]基于4次全国森林资源清查数据并结合生物量估算模型和植被含碳系数, 研究了长江流域森林植被碳储量、碳密度分布特征及动态变化。Qu等^[24]研究了1982—2015年间以NDVI为表征的长江流域生长季植被覆盖度变化, 结果显示受益于生态修复工程引起的土地利用变化, 长江流域自1994年以来植被覆盖度显著上升。就气候驱动因子来看, 大多研究表明温度是决定长江流域植被覆盖度的一个控制因子, 而降水的影响相对较低^[24-26]。这些研究对于揭示长江流域陆地生态系统及其生产力时空格局特点、变化过程及驱动机制等具有重要意义。GPP研究能够量化植被的生长状况, 是全球碳循环研究的重要组成部分, 但目前对长江流域陆地生态系统GPP时空变异及控制因素的系统研究尚不多见。此外, 长江流域地形地貌、气候条件空间分异显著, 植被GPP变化与气候因子、土地利用类型之间关系及其空间差异如何还不清楚。为此, 本文拟利用MODIS系列产品, 并结合土地利用和气象观测数据, 研究近年来长江流域GPP空间分异特征和时间变化规律, 探讨不同土地利用类型GPP的变化差异, 客观评价主要气候因子对陆地生态系统植被GPP的影响及其空间差异, 以期为了解变化环境下区域植被生长状况和生物固碳能力、掌握生态环境质量提供科学基础和参考依据。

1 研究区概况

长江是我国第一大河, 世界第三长河, 发源于青藏高原唐古拉山各拉丹冬雪山。长江干流自西向东流经青海、西藏、云南、四川、重庆、湖北、湖南、江西、安徽、江苏等地, 于上海市注入东海, 干流全长近6400km。长江流域总面积约180×10⁴ km², 占全国总面积的18.75%(图 1)。流域地势西高东低, 由河源至河口, 总落差约5400m, 跨越中国地势的三大阶梯。流域大部分区域属于亚热带季风气候, 在青藏高原的部分属于高山高原气候。流域多年平均气温呈东高西低、南高北低的分布趋势。位于江源附近气温极低, 年平均气温在-4℃上下, 而中下游大部分地区年平均气温在16—18℃之间。流域多年平均降雨量约1100mm, 年降水量的时空分布很不均匀, 空间上呈由西北向东南递增的趋势。

长江流域支流水系众多, 不同水系子流域内地形地貌和气象条件差异显著。根据全国水资源管理区划, 长江流域划分为12个二级水资源区: 金沙江石鼓以上流域(JSJ-1)、金沙江石鼓以下流域(JSJ-2)、岷沱江流域(MTJ)、嘉陵江流域(JLJ)、乌江流域(WJ)、宜宾至宜昌段长江干流区间(UM)、洞庭湖流域(DL)、汉江流域(HJ)、鄱阳湖流域(PL)、宜昌至湖口段长江干流区间(MM)、湖口以下长江干流区间(LM)、太湖流域(TL)等(图 1)。本文后期将以此分区为依据, 具体分析长江流域不同区域GPP变化的空间差异。

2 数据及方法 2.1 数据来源及处理

本文所用的GPP数据来自于美国蒙大拿大学密苏拉分校地球动态数值模拟研究组(Numerical Terra dynamic Simulation Group-NTSG)发布的2000—2015年的MOD17A3全球GPP数据产品(Version 55)(下载网址为: http://www.ntsg.umt.edu/project/modis/mod17.php), 空间分辨率为30弧秒(0.0083°), 数据格式为GeoTIFF, 投影格式为WGS84投影, 比例系数为0.1。该数据的GPP估算模型参考利用了BIOME-BGC模型与光能利用率模型, 模型模拟得到陆地生态系统日GPP, 再进一步推算出8 d GPP值与年GPP值^[27-28]。与美国国家航空航天局(NASA)提供的8天和年GPP数据(V5)产品相比, 该数据很好的消除了云对植被叶面积指数和光合有效辐射值的影响, 提高了数值精度^[27-29]。目前, 该数据已与全球多个地区的通量站点数据进行了对比验证, 在全球和区域植被生物量估算、碳循环和全球变化等研究中得到广泛应用。本研究通过对该数据进行格式转换、影像裁剪等处理, 得到长江流域2000—2015年GPP年值数据集。

土地利用数据来源于中国科学院西部生态环境数据中心提供的2010年全国土地利用1km栅格数据。按照二级分类体系, 将长江流域土地利用类型分为水田、旱地、有林地、灌木林、疏林地、其他林地、高覆盖度草地、中覆盖度草地和低覆盖度草地等9种, 水域、裸地以及城镇等其他类型均归纳为无植被覆盖区类型(图 1)。

论文所有气象数据由中国气象科学数据共享服务网(http://cdc.nmic.cn)提供, 包括整个长江流域内共计175个气象站点的蒸发皿蒸发量、气温、降水、相对湿度、日照时数和风速等月值数据, 数据时段为2000—2015年。各站点气候因子年平均数据由月数据汇总所得。长江流域及各二级水资源区各气候因子年平均值均由各自区域内所有气象站点实测值通过算术平均法求得。

2.2 研究方法

采用简单的一元线性回归模型(y=αt+β, α为回归系数, β为常数, t为年份, y为待分析变量)来分析2000—2015年16年期间长江流域GPP整体变化趋势及其空间异质性。该模型中的回归系数α即为线性趋势斜率, 其中, α < 0表示变量在计算时段内呈下降趋势, 反之α>0表示增加趋势。对于线性变化趋势的显著性, 采用Mann-Kendall(MK)统计检验法^[30]进行检验。根据该方法, 在0.1、0.05和0.01显著性水平下, 如果标准化统计量| Z |≥1.64、| Z |≥1.96和| Z |≥2.32, 则不接受无趋势的零假设。此外, 研究中通过计算标准偏差来反映植被GPP的空间变异程度。

本研究通过计算偏相关系数来分析主要气候因子与植被GPP年际变化之间的本质联系。因为在多元相关分析中, 变量之间的相互关系很复杂, 它们可能受到不止一个变量的影响, 简单的线性相关系数不能真实的反映出两变量之间的本质联系。偏相关分析也称净相关分析, 它是在控制其他变量的线性影响的条件下分析两变量间的线性相关性, 所采用的工具是偏相关系数(净相关系数)。根据偏相关系数, 从而可以判断自变量对因变量的影响程度。

3 结果分析 3.1 GPP空间分布特征

图 2显示了2000—2015年长江流域植被覆盖区多年平均GPP的空间分布。由图可知, GPP总体上东南高西北低, 由东南向西北递减。GPP较大的区域主要位于鄱阳湖流域的南部和东北部, 此外金沙江石鼓以下沿长江一定范围内的局部区域GPP也相对较大。空间上, GPP最大值为2765 gC m^-2 a^-1, 最小值为0.3 gC m^-2 a^-1, 流域整体均值约990.46 gC m^-2 a^-1。统计结果显示, 累计约65%的区域GPP集中在800—1500 gC m^-2 a^-1之间, 25%的区域GPP均值小于800 gC m^-2 a^-1, GPP均值大于1500 gC m^-2 a^-1的区域不足10%(图 2)。

图 3显示了海拔每隔100m统计的长江流域GPP的变化规律。由图可知, 当海拔低于2300m的时候, 植被GPP整体较高, 在1000—1300 gC m^-2 a^-1之间, GPP随海拔升高呈现出先增大、后降低、然后再增大的过程。在海拔300m和2300m左右, 出现GPP的两个峰值。当海拔高于2300m的时候, GPP整体上随海拔升高呈现急剧减小的趋势。海拔5000m以上, 有植被的地方GPP已经降到很小, 不足80 gC m^-2 a^-1。

就长江流域12个二级水资源区来看(表 1), 鄱阳湖流域单位面积GPP最高(1361.83 gC m^-2 a^-1), 其次是洞庭湖流域(1255.99 gC m^-2 a^-1)和宜宾-宜昌段长江干流区间(1178.41 gC/m² a^-1), 最低的是金沙江石鼓以上流域(231.63 gC m^-2 a^-1)。值得注意的是, 宜昌至湖口段长江干流区间、湖口以下长江干流区间和太湖流域单位面积GPP大体相当；金沙江石鼓以下流域因南部沿长江一定范围内的局部区域GPP较高, 使得整体单位面积GPP明显高于金沙江石鼓以上流域和岷沱江流域；岷沱江、嘉陵江和汉江等长江干流以北流域单位面积GPP小于长江干流以南流域。从标准差来看, 金沙江石鼓以下流域、岷沱江流域和鄱阳湖流域较大, 分别514.86 gC m^-2 a^-1、385.85 gC m^-2 a^-1和376.50 gC m^-2 a^-1。表明这三个流域植被GPP的离散程度大, 空间变异剧烈。考虑到不同区域的面积大小, 各区域GPP总量差异更为突出。计算结果表明, 整个长江流域多年平均GPP总量为1.735 P gC, 其中洞庭湖流域GPP总量为0.325 P gC, 占比约18.76%。其次, 金沙江石鼓以下流域和鄱阳湖流域, 其GPP总量分别为0.265 P gC和0.216 P gC, 占比约15.28%和12.42%。太湖流域因面积小, 其GPP总量为0.029 P gC, 仅占长江流域GPP总量的1.69%。

3.2 GPP年际变化特征

2000—2015年间长江流域植被覆盖区单位面积GPP总体上呈波动增加趋势(图 4)。其中, 2000—2001年间GPP相对偏低, 2002—2008年间GPP整体较高, 2009—2012年间GPP出现一个低谷, 之后的2013—2015年间GPP又恢复较高水平。趋势分析可知2000—2015年间长江流域GPP年际变化线性趋势率为2.39 gC m^-2 a^-1, 上升趋势不显著(P>0.1)。

空间分析结果表明(图 5), 长江流域累计约68%的区域GPP线性变化趋势率大于0, 反映整个长江流域绝大部分地方GPP呈上升趋势, 其他约32%的区域GPP呈下降趋势。具体来看, 嘉陵江和汉江流域北部区域、金沙江流域南部区域GPP上升趋势最为明显, 太湖流域、洞庭湖流域以及四川盆地大部区域GPP下降最为突出。

各二级水资源区中, 除了洞庭湖流域和太湖流域单位面积GPP呈下降趋势外, 其他各区GPP均呈上升趋势(表 1)。其中, 太湖流域GPP年际变化显著(P < 0.05), 线性趋势率为-8.81 gC m^-2 a^-1。在所有GPP上升的区域中, 金沙江石鼓以下流域和汉江流域GPP年际变化显著(P < 0.05), 线性趋势率分别为5.49 gC m^-2 a^-1和5.82 gC m^-2 a^-1。另外, 金沙江石鼓以上流域单位面积GPP也有相对显著的上升趋势(P < 0.1)。

3.3 不同土地利用类型GPP变化比较

长江流域面积大, 不同土地利用类型的GPP分布情况不同, 本文依据二级分类体系, 对研究区不同土地利用类型的地表植被进行了分类研究。图 6所示的统计结果表明, 耕地中的水田和旱地植被多年平均GPP大体相当, 分别为1075.80 gC m^-2 a^-1和1064.10 gC m^-2 a^-1。耕地GPP整体上略低于林地。林地各类型GPP介于1122.80—1276.90 gC m^-2 a^-1之间, 其中有林地最高, 然后依次为疏林地、其他林地和灌木林。草地各类型GPP明显低于耕地和林地, 而且不同覆盖度的草地GPP相差很大。高、中、低覆盖度草地GPP分别为896.70 gC m^-2 a^-1、671.80 gC m^-2 a^-1和198.50 gC m^-2 a^-1。标准差结果反映高覆盖度草地和中覆盖度草地GPP离散程度最大, 空间变化程度最为突出。水田、旱地GPP标准差最小, 说明这两种人工植被土地利用类型GPP空间变异最小。

图 6显示的2000—2015年间不同土地利用类型的单位面积GPP趋势变化中, 除水田GPP呈微弱下降外, 其他各土地利用类型的GPP均呈上升趋势。其中, 线性变化趋势率以高、中覆盖度草地为最, 分别为4.63 gC m^-2 a^-1和6.41gC m^-2 a^-1, 随后为灌木林(3.35 gC m^-2 a^-1)、旱地(3.19 gC m^-2 a^-1)、疏林地(3.07 gC m^-2 a^-1)、有林地(2.15 gC m^-2 a^-1)、其他林地(1.62 gC m^-2 a^-1)和低覆盖度草地(1.47 gC m^-2 a^-1)。显著性检验表明, 高覆盖度草地和中覆盖度草地GPP上升趋势显著(P < 0.05), 低覆盖度草GPP上升趋势也较大(P < 0.1)。

3.4 GPP时空变化与气候因子的关系

气候条件是影响陆地生态系统GPP变化的重要环境因素。本文选取蒸发皿蒸发量、降水、相对湿度、日照时数、风速和气温等6个变量作为长江流域GPP变化的潜在气候驱动因子。表 2统计了2000—2015年长江流域及各二级水资源区GPP年际变化与主要气候因子的偏相关系数。总体上, 蒸发对植被GPP变化起到一定的负作用(除太湖流域外), 温度和日照主要起正作用, 降水、相对湿度和风速对GPP变化的影响有正有负。不同气候因子对植被GPP年际变化的影响程度存在明显的区域差异。长江上游源头的金沙江石鼓以上流域和岷沱江流域GPP与温度呈显著的正偏相关关系, 其他气候因子与GPP的偏相关关系均不显著。长江上游金沙江石鼓以上流域和嘉陵江流域的蒸发对GPP年际变化的负作用更为突出, 其偏相关关系均较显著。另外, 嘉陵江流域降水与GPP之间的偏相关关系也达到显著性水平, 且其对GPP年际变化的正作用超过蒸发的负作用。乌江流域和长江上游干流区间, 气候因子对GPP年际变化的影响作用十分复杂。其中, 乌江流域相对湿度对GPP年际变化起显著负作用, 其次是温度和降水, 尽管两者与GPP之间的偏相关系数为正, 但其影响程度均不及相对湿度, 其他因素影响不显著。长江上游干流区间蒸发、降水、风速和温度与GPP之间的偏相关系数都较显著, 但是蒸发和风速对GPP年际变化起负作用, 而降水和温度起正作用, 并且其影响程度要强于前者。长江中游洞庭湖流域, 温度和风速与GPP的偏相关关系均较为显著, 其中温度影响的正作用要强于风速的负作用。汉江流域GPP年际变化主要受蒸发的负作用影响, 两者之间呈显著的负偏相关关系, 其他气候因子的影响均不明显。鄱阳湖流域各气候因子与GPP之间偏相关系数均不显著, 但总体上温度的影响强于其他因子。长江中游干流区间的降水、温度与GPP之间呈显著的偏相关关系, 两者对GPP年际变化起到较强的正作用。长江下游干流区间和太湖流域各气候因子与GPP之间偏相关系数均不显著。对长江流域整体而言, GPP年际变化主要受温度影响, 两者之间呈显著的正偏相关关系, 其他气候因子的作用均不明显。

对于不同土地利用类型, 植被GPP年际变化与主要气候因子的偏相关系数见表 3。由表可知, 不同土地利用类型GPP年际变化同样主要受温度和蒸发影响, 其中温度对植被GPP变化起正作用, 蒸发起到一定的负作用(除低覆盖度草地外)。林地各主要类型以及高、中覆盖度草地的GPP与温度之间的偏相关关系达到显著性水平, 而耕地(水田与旱地)与低覆盖度草地的GPP与温度之间的偏相关关系不显著。蒸发的负作用对高覆盖度草地GPP的影响十分显著, 而对其他土地利用类型的影响并不突出。对所有土地利用类型而言, 降水、相对湿地、日照和风速对GPP变化的影响均不显著。

4 讨论

水分和温度及其相互配合构成的水热条件是影响植被生长和分布的主要因素^[31]。总体上, 长江流域陆地植被GPP的空间分布特征与流域水热条件的分布情况高度一致。然而, 值得注意的是, 长江流域植被单位面积GPP随高程变化, 在海拔300m和2300m左右存在两个峰值。对于这一现象, 主要是受以下因素影响: (1)长江流域海拔300m左右的高程区, 主要分布于长江中下游和四川盆地的低山丘陵, 地表植被以亚热带常绿阔叶林为主。长江中下游和四川盆地海拔300m以下地区大多的以水田和旱地为主, 而耕地GPP整体低于林地。300m以上区域, 气温受垂直地带性影响随海拔升高逐渐降低, 从而不利于植被生长；(2)海拔2300m左右的高程区, 主要分布在我国地势一、二级阶梯的交界处, 集中分布在四川盆地和云贵高原向青藏高原过度的横断山区外围地带。这里山高谷深, 地表植被主要分布于河谷两岸的山坡地带, 以高山阔叶林为主。由于河谷地带大多气候干热, 局地“焚风”气候效应形成了独特的垂直植被带倒置现象: 河谷低海拔地带因气候干热植被稀少, 而海拔较高的山坡地带随降水增多, 植被生长逐渐转好。从植被GPP随高程的变化结果来看, 西部干热河谷区海拔2300m处整体水热条件最好。张文江等^[32]对岷江上游干热河谷植被分布的研究也表明, 岷江上游以2000—3000m高程带水热组合条件最好, 分布了川滇高山栎等阔叶林；3000m以上逐渐受热量制约, 植被依次为针叶林、灌木、草甸和高山植被；2000m以下受干热河谷效应及人类活动影响, 植被以灌木类为主。

以往的研究中, 大多注重气温和降水对植被生产力的影响^{[10, 11, 13-14, 16]}。本研究除气温和降水外, 还系统探讨了蒸发、湿度、日照和风速等气候因子的协同作用。对长江流域整体而言, GPP年际变化主要受温度影响, 这与众多相关研究的结论一致, 因为温度对植被光合作用和生产力的促进作用最为突出^[33]。此外, 由于长江流域地处亚热带湿润气候, 流域水分条件好, 能够充分满足植被生长的需要, 使得降水和相对湿度的变化对GPP的影响作用较弱, 而日照、风速等因素可以通过影响蒸发过程从而对植被GPP产生作用^[34]。总体上, 陆地生态系统的组成、结构和功能受气候干湿状况的影响较为突出, 特别是干旱对植被生长极为不利。杜文丽等^[35]的研究表明, 在时间尺度上, 1980—2013年中国陆地生态系统GPP与SPEI年际变化格局基本吻合, 干旱年份的GPP显著降低；在空间尺度上, 北方大部分地区的GPP与SPEI呈正相关, 南方大部分地区呈负相关。Zhao和Running^[29]的研究也表明, 2000—2009年间全球陆地植被NPP的减少是干旱引起的。需要指出的是, 受长江流域内部地形地貌及气象条件分异显著的影响, 不同气候因子对GPP年际变化的影响存在明显的空间差异。本文的研究结果表明, 在不同子流域内, 蒸发、降水、相对湿度、风速和温度等都有可能成为GPP年际变化的主要气候因子(表 2)。对不同土地利用类型而言, GPP年际变化主要受温度和蒸发的影响, 降水等其他因素的影响较小(表 3)。因此, 有关气候变化对植被生产力的影响研究中, 不应简单的分析气温和降水的影响, 而应考虑流域所处气候条件, 分析其他气候因子或流域干湿状况的综合作用。

在时间上, 长江流域陆地植被GPP呈波动增加趋势, 流域累计约68%的区域GPP线性变化趋势率大于0。就不同土地利用类型来看, 除水田GPP呈微弱下降外, 其他均呈上升趋势, 尤其是草地GPP上升趋势较为显著。这些结果表明长江流域整体的固碳能力在逐步增强, 生态系统环境质量有逐渐改善的趋势, 这与近年来长江流域植被覆盖度显著上升有较好的对应^[24]。从流域气候变化的角度来看, 近几十年来长江流域气温上升明显, 局部地区降水也呈现出一定的增加趋势^[36], 这对提高流域植被GPP, 促进生态系统质量改善具有积极作用, 但应注意气候过程的极端化对生态系统可持续发展的潜在风险。此外, 流域内人类活动对植被GPP的影响也不容忽视。特别是近年来, 随着城市化的快速扩张, 部分区域地表覆被状况发生了显著改变, 从而影响地表植被生产力和生态系统服务功能的时空变化^[37-38]。另一方面, 在国家生态文明建设的倡导下, 长江流域积极开展生态修复工程, 导致植被覆盖度显著上升, 流域生态系统质量得到显著改善^[24]。为客观认识变化环境下陆地生态系统固碳能力和生态服务功能的变化, 积极开展气候和土地利用变化对植被生产力和生态系统状况的影响研究具有重要的现实意义。

5 结论

本文以MODIS GPP数据为基础, 系统研究了长江流域2000—2015年植被GPP的时空变化特征, 并探讨了植被GPP年际变化与主要气候因子的关系, 获得结论如下:

(1) 长江流域植被覆盖区GPP在0.3—2765 gC m^-2 a^-1之间, 均值约990.46 gC m^-2 a^-1。植被GPP在海拔300m和2300m左右存在两个峰值, 分别为1269.30 gC m^-2 a^-1和1251.80 gC m^-2 a^-1。当海拔高于2300m的时候, GPP随海拔升高急剧减小。不同水资源区内, 单位面积GPP以鄱阳湖流域最高, 金沙江石鼓以上流域最低。长江流域多年平均GPP总量为1.735 P gC, 其中洞庭湖流域总量最高, 约占整个长江流域总量的18.76%, 太湖流域最小, 占比仅1.69%。

(2) 长江流域GPP在2000—2015年间呈缓慢的波动上升趋势, 趋势率为2.39 gC m^-2 a^-1。空间上, GPP上升区和下降区分别占总流域面积的68%和32%。对于各二级水资源区, 除洞庭湖流域和太湖流域GPP呈下降趋势外, 其他各区均呈上升趋势, 尤其是金沙江流域和汉江流域的上升趋势较为著性。

(3) 长江流域不同土地利用类型GPP平均值在198.50—1276.90 gC m^-2 a^-1之间。其中, 草地各类型空间变化程度最为突出, 水田和旱地空间变异程度最小。各土地利用类型中除水田GPP呈微弱下降外, 其他均呈上升趋势, 尤其是高、中、低覆盖度草地GPP上升趋势显著。

(4) 总体上, 主要气候因子中蒸发对植被GPP变化起到一定的负作用(除太湖流域外), 温度和日照主要起正作用, 降水、相对湿度和风速对GPP变化的影响有正有负。不同气候因子对植被GPP变化的影响程度在不同区域、不同土地利用类型间均存在一定差异, 但就长江流域整体而言, GPP年际变化受温度影响显著, 其次是蒸发, 而降水等其他气候因子的影响不大。