近40年青藏高原生态格局演变及其驱动因素

王子滢,李周园,*,董世魁,符曼琳,李泳珊,李生梅,武胜男,马春晖,马天啸,曹 越

1 北京林业大学草业与草原学院,北京 100083

2 中国科学院森林生态与管理重点实验室(中国科学院沈阳应用生态研究所),沈阳 110016

3 清华大学建筑学院景观学系,北京 100084

青藏高原是世界上海拔最高、分布面积最大、自然地理环境最复杂的高原,是长江、黄河、雅鲁藏布江、怒江、澜沧江以及恒河和印度河等亚洲主要河流的发源地,被誉为“世界屋脊”、“地球第三极”与“亚洲水塔”。青藏高原强烈隆升所带来的与周边区域的海拔差异影响着中国乃至亚洲地区的大气环流体系,形成了从热带到寒带,从湿润到干旱等多种气候类型[1—3],造就了藏高原地区丰富的生物区系及生态系统多样性,从东南向西北依次分布着森林、灌丛、草地、荒漠等自然生态系统,镶嵌分布着湿地生态系统、冰川生态系统,在农区与牧区之内包含着不断发展扩张的城镇与农田生态系统。青藏高原整体是中国乃至亚洲的重要生态安全屏障,在水源涵养、水土保持、气候调节、生物多样性保护、固碳增汇等方面发挥着基础而关键的作用[4]。

20世纪80年代至今的近40年间,在全球变化与人类活动加剧的作用下,青藏高原的生态系统发生着广泛而深刻的变化,敏感脆弱的高寒草甸、荒漠草地出现大规模退化,同时高原局地经历着不同程度的冰川消融、城镇化、农业垦殖、草地灌丛化以及生态修复保护等一系列空间格局的动态演变,不同生态系统类型此消彼长、相互转变。过去40年间也是遥感与地理信息技术迅猛发展与海量环境数据积累的阶段,在全球关键区的生态环境变化监测与宏观生态学机制研究中取得长足进展。利用卫星遥感定期观测和解译的土地覆被与分类数据,对青藏高原生态系统的宏观时空格局动态研究,主要集中于三个方面：(1)气候变化下以高寒草地为主自然植被生态系统面积与质量变化[5—9],(2)人为活动对高原生态系统的扰动与干预过程的监测与评估[10—12],(3)青藏高原生态系统变化驱动机制的量化分析。以往的研究针对某一类或几类相对独立的生态系统演变进行剖析,对于全域各类生态系统的横向比较和归纳甚少；对于外部驱动因素的分析指标,以往很多研究选取自然或人为活动中的几个或一个变量相对独立地进行统计分析,未见有选取综合的驱动因素开展多层次的机制分析。特别地,随着近年来时空数据分析与多元统计建模技术方法的发展与普及,对于生态系统格局演变驱动机制的理论和实证探讨朝向多元化和深层次发展,对于兼有时间和空间维度信息的关键变量的转换与分析,成为解析和阐发生态系统状态与过程之间因果耦合的有效途径[13]。

青藏高原之所以是我国生态较为脆弱的区域是由于该区植被对气候变化以及人类活动的影响较为敏感[14],且近年来青藏高原整体暖湿化特征明显,其中青藏高原大部分地区呈现降水量每10年增加2.2%的趋势[15—16]。在自然以及人为因素的双重影响下,青藏高原地区各类生态系统出现不可避免的相互转化[17],如气候变暖可导致冰川、草地、荒漠生态系统向湿地生态系统转变[18—19]。当前,相关学者基于对地观测与地理环境数据对青藏高原典型区域的生态系统时空变化及驱动因素开展了一系列研究[20]。Tong等[21]发现湿度差异和降水、蒸发变化引起的年际动态是促成三江源区域水体与湿地发生格局改变的主要驱动因子。而引起玛曲草地生态系统面积消减的主要驱动因子为人类活动的增加[22]。张镱锂等[23]认为,珠穆朗玛峰自然保护区内冰川等自然生态系统对于人类活动以及自然因子的影响均表现出较为迅速的响应。目前也有针对青藏高原全局尺度的土地利用景观格局[24]以及土地覆被变化驱动力[25]的相关研究,但深入解析青藏高原地区生态格局演变驱动因素,仍有待进一步采取综合的自然与人为活动指标,系统地构建多层次归因模型,量化外部环境对于具有不同特征的各类自然与人工生态系统变化的贡献比率。

本研究基于土地利用覆被变化数据、气候数据,对青藏高原地区全域9类主要生态系统,结合人为活动强度、土壤侵蚀度和生物丰度等综合化的驱动因素空间数据集,对该区域的生态系统宏观演变特征、环境变化格局进行研究,利用结构方程模型定量计算驱动作用关系,旨在探究如下科学问题：(1)青藏高原1980—2018年间生态系统的分布格局在时空上是如何演变的？(2)其变化的主要驱动机制如何？本研究揭示的高原生态系统演变及驱动因素的作用路径,为系统剖析过去近40年间青藏高原生态系统的演变与驱动因素以及发展保护该地区生态屏障功能对策提供实证理论依据与科学参考。

1 材料与方法

1.1 研究区域

青藏高原位于中国的西南部,范围为26°00′12″—39°46′50″ N,73°18′52″—104°46′59″ E,总面积为2.6×106km2,占全国陆地总面积的26.8%[26],主要包括青海省和西藏自治区全域,以及甘肃省西南部地区、四川北部地区、云南西北部、新疆维吾尔自治区南部地区[27]。青藏高原地处高原气候区,全域最冷月平均气温低达-15—10℃,大部分地区最暖月均温<10℃[28]。从东南部海拔3500 m以上的横断山脉到中部海拔5000 m以上的唐古拉山,再到西北缘海拔6000 m以上的昆仑山[29],依次分布着森林、灌丛、草地、荒漠生态系统等,水体与湿地、荒漠、冰川、农田以及城镇生态系统镶嵌其中(图1)。

图1 青藏高原海拔图与生态系统类型分布

1.2 数据来源与处理方法

本研究中8期(1980、1990、1995、2000、2005、2010、2015、2018年)栅格土地利用数据由中国科学院资源环境科学与数据中心提供(http://www.resdc.cn/),分辨率为30 m。该数据包含6个一级类型以及25个二级类型,对原始不同土地利用与覆被类型依生态系统属性进行划分,得到9类青藏高原的生态系统类型,即森林、灌丛、草地、农田、城镇、冰川、水体与湿地、裸地、荒漠。将该9类生态系统中的森林、灌丛、草地、水体与湿地、冰川、裸地、荒漠,依据其主体物质构成划分为以植被为主体和以无机环境要素为主体的两类自然生态系统,前者地表形态以植物覆盖的有机生命体为主组成,而后者以不同形态的无机物或有机非生命体组成为主,即水、冰、二氧化硅、碳酸钙、腐殖质等。另外的农田与城镇作为人工生态系统单作一类(表1)。

表1 青藏高原土地覆被类型及其生态系统类型归类

年均气温、年均降水数据(1980—2015年,共36期)来源于中国科学院资源环境科学与数据中心平台提供的全国气象插值数据产品,空间分辨率为1 km。数字高程数据来源于北京大学城市与环境学院地理数据平台(http://geodata.pku.edu.cn),空间分辨率为250 m,利用数字高程图在QGIS 3.16里进行坡度计算,生成以百分比为单位的坡度栅格图层,以反映区域对生态系统形成与变化影响的海拔、坡度因素。土壤侵蚀度数据来源于中国科学院资源环境科学与数据中心平台,空间分辨率为1 km,用于定量反映土地要素对区域生态系统格局演变的作用。根据生物多样性与生态系统功能理论[30—31],考虑到生物丰度格局对于生态系统功能变化过程有直接效应,选取以2005年土地覆被数据为基础、依据生态环境质量评价技术规范所设置的转化参数,对不同生态系统赋予不同权重,生成生物丰度指数空间分布栅格图层,该变量指单位面积上不同生态系统类型在生物物种数量上的差异,间接地反映被评价区域内生物的丰贫程度,空间分辨率为1 km[32]。人类活动相对强度数据来源于荒野度指数(Wilderness Quality Index,WQI)数据产品,已公开发表[33],空间分辨率为1 km。该图层基于来自2015—2018年间的相关数据集,通过测度土地利用自然度、人口密度、距居民点遥远度、距道路遥远度、居民点密度和道路密度6项指标并进行加权线性叠加得到荒野度。利用该数据集,取差值方式计算得到量化人类活动相对强度的综合指数(Human Activity Index,HAI),

HAI=1-WQI

(1)

以上图层以第二次青藏高原综合科学考察所划定的矢量边界进行统一裁剪并叠置图层,以高斯克吕格统一投影坐标系,以便后续计算和提取数值分析。

1.3 统计分析

1.3.1生态系统格局演变量化

由于相邻年份之间各类生态系统变化幅度不明显,本研究选择2000年为间隔,分别计算1980—2000年、2000—2018年以及1980—2018年期间不同类型生态系统面积相比较起始年份的变化比率。同时,青藏高原生态系统面积建立分辨率为1 km的空间网格,对39年间每1 km网格内9类典型生态系统所占百分比的动态趋势进行线性回归,以每10 a百分比(%/10a)为单位,用来反映不同生态系统宏观动态消长变化速率,利用Python语言编程处理得到生态系统动态变化率栅格图层。

1.3.2驱动因素变化速率计算

驱动因素数据空间分辨率可综合时间序列长度以及空间尺度特征,结合研究需要进行选择,且分辨率的大小对最后结果无显著影响[34],因此本文选择将气温、降水、海拔、坡度、人为活动指数信息以及土壤侵蚀度、生物丰度指数等地理信息数据产品重采样,建立分辨率为50 km的空间网格。对于动态驱动因素(气温、降水),分别计算每个网格点平均气温和降水的动态变化速率,即以年均值作为因变量、时间序列的年份作为自变量进行线性回归,取得气候因子斜率值分布,作为气候驱动因子的量化结果。

1.3.3驱动因素综合分析

本研究将上述处理得到的生态系统格局变化速率、气候变化速率、海拔、坡度、人为活动指数以及土壤侵蚀度、生物丰度指数的图层数据汇总整合,作为驱动因素,即自变量,以三大类生态系统的变化率加和作为量化生态格局演变的作用效益,即因变量,分别在R语言进行结构方程模型计算。其中,结构方程模型构造中,根据以往研究基础和知识提出先验假设,以生态系统变化率作为解释因变量的终点目标,气候变化、地形地貌与人为活动既对青藏高原生态系统动态演变有直接作用,又可以通过影响植被地下的土壤侵蚀程度和植被本身的生物多样性分布来间接作用生态系统的演变速率,即构造一套以气候变化、地形地貌以及人为活动强度为外生变量、其他因子为内生中间变量的多层次结构方程模型(图2),利用分步回归的方式在‘piecewiseSEM’程序包里,利用采集得到的栅格图层数据进行统计建模、检验并计算变量之间的路径系数,完成驱动因素分析。

图2 青藏高原生态系统动态变化与驱动因子关系的概念图

2 结果与分析

2.1 青藏高原生态系统分布特征

青藏高原地区1980—2018年间的生态系统分布总体变化情况如表2所示。由表可知,由自然植被要素为主的生态系统(森林、草地、灌丛)是青藏高原最主要的自然生态系统,1980年三者占青藏高原总面积的61.7%,2018年时三者面积占比变化不大,约为61.9%。以无机环境要素为主的生态系统(水体与湿地、冰川、荒漠、裸地)1980年整体占青藏高原全区总面积的36.0%,2018年占青藏高原全区面积的35.8%。对于人工生态系统(城镇、农田)而言,1980年代城镇与农田生态系统面积仅占青藏高原全域的2.2%,2018年二者总面积占青藏高原全域的2.3%。

表2 1980—2018年青藏高原生态系统面积/(×104 km2)

2.2 青藏高原生态格局时空变化特征

青藏高原生态格局时空变化的量化计算结果显示,从表3可知,1980—2000年期间,森林、冰川以及荒漠生态系统面积呈现负增长,其中冰川生态系统减少最为明显,局部缩减达到28.1%,其余生态系统面积均呈现增加趋势,其中人工生态系统(城镇、农田)增长率分别达到了7.9%和1.5%。2000—2018年间,森林、草地、冰川、荒漠以及农田生态系统呈现负增长趋势,其中农田生态系统减少了1.8%,而城镇生态系统依旧呈现较高的增加幅度,达到30.2%,同时水体与湿地生态系统面积增加也较多,达到1.7%。总体说来,以自然植被为主的生态系统在近40年中,森林生态系统面积较1980年减少了0.3%,灌丛与草地生态系统则各增加了0.2%和0.4%；以无机环境要素为主的生态系统整体变化幅度大于自然植被为主的生态系统,其中冰川生态系统面积减少幅度最大,减少了28.2%,水体与湿地与裸地生态系统则分别增加了2.4%和1.8%。人工生态系统中城镇生态系统增幅较大,40年间增加了近40.4%,而农田生态系统则呈现先增后降的趋势。

表3 青藏高原生态系统面积平均变化比率/%

图3展示了各类生态系统变化速率的空间分布,在1980—2018年间,青藏高原东北缘祁连山脉至横断山脉以及冈底斯山脉附近的森林生态系统面积变化呈现增加趋势,西藏高原东南缘至念青唐古拉山脉的森林生态系统面积减少,同时该区域也是灌丛生态系统变化集中发生的地区。青藏高原昆仑山以西至喜马拉雅山以北地区草地生态系统面积的减少率高达33%/10a,其余地区大部分草地生态系统的变化率则呈现正增长。水体与湿地生态系统的时空动态变化则与草地生态系统呈现出相反的态势。冰川生态系统消减速率较快,部分地区可达-25.0%/10a,但有仍有区域冰川生态系统处于扩张趋势,增长速率达到8%/10a。西藏西部裸地与荒漠生态系统增加速率较快,速率分别可达32%/10a以及15%/10a。在青藏高原东部省份内,人工生态系统如农田以及城镇生态系统扩张明显,呈现正增长的趋势。

图3 1980—2018年青藏高原九类主要生态系统变化速率空间分布格局

2.3 青藏高原生态格局演变的潜在驱动机制

图4显示了1980—2018年青藏高原的气温变化速率格局,全域气温变化速率介于-0.01—0.22 0.1℃/10a之间。温度每10年的变化速率在空间上呈现明显的区域分布,西北边缘以及东南地区呈现气温降低的趋势,中部地区气温呈现增加的趋势。图4显示了1980年至2018年青藏高原地区的降水变化速率格局,全域降水变化速率介于-15—25 mm/10a。同样的,降水变化速率也呈现出一定的空间规律,从西部到东部逐渐降低。图4还显示了青藏高原地区的土壤侵蚀度以及人为活动指数的空间格局分布,土壤侵蚀程度较高的地区集中在柴达木盆地以及青藏高原西南部,人为活动指数则呈现出明显的区域分异特征,青藏高原东南边缘人为活动指数超过0.5。

图4 青藏高原生态系统关键自然与人为驱动因子空间格局

结构方程模型分析结果可知(图5),气温、降水、海拔、土壤侵蚀度以及生物丰度对以自然植被为主的生态系统演变的速率的影响是负向显著的,其路径系数分别为-0.21、-0.13、-0.20、-0.34和-0.56,只有坡度对其演变速率的影响是正向显著的,其路径系数为0.25。而对于以无机环境要素为主的生态系统而言,气温、降水、海拔、土壤侵蚀度以及生物丰度对其演变的速率的影响是正向显著的,其路径系数分别为0.20、0.13、0.24、0.34和0.58,只有坡度对其演变速率的影响是负向显著,路径系数为-0.27。仅有海拔显著负向影响人工生态系统的演变速率,其路径系数为-0.13。

气温、降水、海拔、坡度以及人为活动的变化能直接对以自然植被为主的生态系统、以无机环境要素为主的生态系统以及人工生态系统的演变速率产生影响,也可通过改变土壤侵蚀度以及生物多样性间接对三大类生态系统的演变速率产生影响。通过直接连接和间接连接路径乘积并按不同驱动因素加和计算可知(图5),以自然植被为主的生态系统的演变速率与除坡度以外的环境因子变化之间的耦合关系呈现负反馈和保守性的变化,以无机环境要素为主的生态系统的演变速率与环境变化之间的耦合关系则呈现相反的态势。综合来看,直接驱动因子中气温变化因子以及间接驱动因子中的生物丰度因子对不同脆弱性的生态系统演变速率有着较大的影响。而对人工生态系统而言,只有海拔因子对于其演变速率有着相对较大的影响。

图5 青藏高原生态格局演变驱动因素分析的结构方程模型量化结果

3 讨论

引起青藏高原地区生态系统格局演变的因素错综复杂,自然以及人为因素则是造成其变化的主要原因。通过本研究发现,气温是本研究分析的驱动因素中,对青藏高原生态格局变化速率影响最大的因子。近40年来青藏高原地区整体增温现象明显,表现出高于全球增温幅度的暖化现象[35],青藏高原地区的降水情况则存在区域分异的情况,大部分地区降水量每10年增加2.2%,局部地区如青藏高原东部以及东南部地区年降水量则产生下降趋势[15—16, 36],气候整体呈现暖湿以及冷干的组合。本研究显示气温以及降水可以直接或者通过改变土壤侵蚀程度或是生物丰度间接影响不同类型生态系统格局的演变速率,同时气温、降水还可以通过影响土壤的矿化程度[37]、湿度[38],或是植被物候[39]、植被指数[40]来进一步驱动生态系统格局的演变。本研究中以自然植被为主的生态系统(森林、灌丛、草地)演变速率随着气温、降水以及海拔变化幅度的增加而减慢,抵抗力较强,表现出相对保守的演变特征。造成该现象的原因主要有以下三方面：(1)该类生态系统具有相对复杂的群落结构、多样的物种以及多层次的营养级,在面对外界扰动时,具有更强的抗干扰能力和稳定性[41—43]。在考虑降水量的条件下,森林、灌丛生态系统所在地区拥有更低的土壤侵蚀强度,较高植被覆盖度的草地生态系统也可以缓解一部分由降水所带来的侵蚀力[44]。(2)气候的改变创造了更有利于森林、灌丛、草地生态系统面积增加的条件。由于青藏高原地区森林、灌丛与草地生态系统对于热量及水分的需求处于亏欠的状态,气温与降水的增加大大缓解了它们对于水分及热量的需求,降低灌丛因低温而导致死亡的概率[45],增加草地净初级生产力[15],使其能够保持原有的状态正常生长,并不会出现极大速率的改变。(3)暖湿以及冷干的气候分异导致青藏高原生态系统呈现总体趋好,局部变差的态势。随之可能带来的就是部分地区生态系统分布面积增加,部分地区生态系统退化,面积减少[15],再加之大部分森林、灌丛、草地生态系统在这三者内相互转移,因此从青藏高原这一大范围来看,虽局部地区变化明显,但整体依旧呈现出相对稳定的格局状态。

本研究中以无机要素环境为主的生态系统(水体与湿地、冰川、荒漠、裸地)表现出的韧性有限,随环境变化而改变迅速,对环境扰动敏感,具有较高的脆弱性。造成该趋势的原因可能有以下两方面：(1)以水为组成主体的冰川以及水体与湿地生态系统对气温变化极为敏感。在近几十年内,81%的冰川生态系统处于退缩状态[46],冰川融水以及地表径流量变化将导致水体与湿地生态系统分布呈现局部扩张的格局[47]。(2)该大类生态系统植被覆盖较为稀疏,生态系统稳定性较差。青藏高原整体“暖湿化”的趋势可能会在青藏高原局部地区带来集中降水,冰雹与暴雪的产生将会对拥有稀疏植被覆盖的荒漠以及裸地生态系统表面产生严重影响,如在地势陡峭区域则会引发泥石流、滑坡等自然灾害,加剧土壤侵蚀程度,对生态系统产生不可逆的破坏。本研究中人工生态系统(城镇、农田)的变化速率除随着坡度增加而减慢以外,受其他自然因素的影响较小。这可能是因为人工生态系统受较强的人为意识的控制,人类选择的低海拔的平缓地区更有利于农田以及城镇的扩张。

总地来讲,可以归纳指出青藏高原以植被为主体的自然生态系统,空间上发生的消长变化反映其对环境的响应主要由地表植被与环境互作,对外部环境改变具有生态幅的适应度,反映无机非生命体所没有的特质,有机生命体组成的植被对气候响应通过有机组分和生理生态过程,其呈现的负反馈关系表示了青藏高原的植被具有典型的自适应和滞后性的复杂系统的特征[48—49]；而对于以无机环境要素为主体的自然生态系统,其宏观动态变化最主要反映的是非有机生命体对外部环境变化的响应,如冰川消融、水流侵蚀、水体和风沙输运等过程,该生态系统下所包含的生命有机体总量低,整体所体现韧性幅度窄,对环境改变的响应反馈即时迅速[50—52]。

人为活动因子作为驱动因素之一放大了自然因子对于生态系统格局演变的作用。本研究得出的结果表明,从青藏高原全域的尺度来看,近40年内人为活动强度对各大类生态系统格局的演变的影响尚不明显,但有研究表明在人口密集区域,人为活动相较于气候影响,对植被施加了更大的压力[53—54]。对生态资源利用的不恰当行为如森林采伐、过度放牧、煤矿开采、基础设施的建设等都对各类生态系统产生了一系列负面影响[47, 55],其中过度放牧等人类活动是导致青藏高原高寒草地加速退化的重要因素[56—57],但是积极的人为活动如生态保护与建设工程、退牧还草工程以及国家公园建设等措施都对生态系统分布格局、生态系统质量起到了正向的调节作用[58—59]。

4 结论

本文通过对1980—2018年归纳后的青藏高原生态系统时空变化特征以及造成其变化的驱动因素分析可得：

(1)近40年期间,青藏高原全域内各类生态系统保持相对稳定,其中以自然植被为主的生态系统(森林、灌丛、草地)面积占比大,总计可达61.9%,以无机要素环境为主的生态系统(水体与湿地、冰川、荒漠、裸地)平均占比分别约为4.6%、1.5%、12.2%、0.2%,人工生态系统(城镇、农田)平均占比分别约为0.1%以及2.1%。以自然植被为主的生态系统(森林、灌丛、草地)面积随年份的增加变化幅度不大,分别变化了-0.3%、0.2%和0.4%,森林生态系统变化速率发生较快地区为青藏高原东部祁连山脉、横断山脉以及念青唐古拉山脉地区,青藏高原整体灌丛增加速率显著大于灌丛面积减小速率,在东南横断山脉附近局部扩张最高可达约7%/10a。草地生态系统波动速率最大,局部可达30%/10a。无机要素环境为主的生态系统(水体与湿地、冰川、荒漠、裸地)中冰川生态系统面积变化率最高,达到了-28.2%,其中喜马拉雅山脉的冰川消减最为明显。人工生态系统(城镇、农田)中城镇生态系统变化幅度最大,增加了约40.4%,其中城镇生态系统在省会地区变化较为明显,青海省东部、四川省、云南省以及西藏自治区东南部的农田生态系统变化较为剧烈。

(2)本研究中气温、降水、海拔、土壤侵蚀度以及生物丰度因子对以自然植被为主的生态系统(森林、灌丛、草地)和无机要素环境为主的生态系统(水体与湿地、冰川、荒漠、裸地)的演变速率有着显著影响,从影响程度来看,气温变化这一直接驱动因素的总路径系数最大。以自然植被为主的生态系统在外部环境驱动因素的影响下反映出具有相对高的抵抗力,具有一定的保守性；而无机要素环境为主的生态系统随环境变化而改变迅速,对环境扰动敏感,具有一定脆弱性。因此根据上述关键发现和实证量化结果,保护面积最大、韧性最强、具有较强恢复潜力的森林、灌丛和草地生态系统将更有利于发挥青藏高原的生态屏障功能。