植被生态系统是陆地生态系统的重要组成部分, 其时空分布格局的演替和变化对全球变化有响应和反馈作用^[1-3]。气温和降水等关键气候因子的时空变化, 通过影响对植物的生理生长和土壤湿度等生境条件^[4-6], 进而导致植被生态系统类型、结构及空间格局的演替变化^[7-11]。定量模拟和解析植被生态系统对全球变化的响应、适应及其反馈机制已成为全球变化及生态学研究领域的热点和前沿问题^[12-15]。因此, 如何利用气候观测数据并结合根据不同全球气候情景模式数据, 模拟预测全球变化影响下的植被生态系统长期演化趋势, 是深入揭示和厘清究生态系统结构、功能及其关键过程对人类活动影响和全球变化的响应与适应机制亟需研究的焦点问题之一^{[12, 16]}。另外, 植被生态系统模型的长足发展, 使在不同尺度上实现植被生态系统结构及其分布的时空模拟逐渐精细和深入^[16-17]。目前关于气候与植被生态系统相互作用的主流模型包括HLZ生态系统模型^[18]、BOX模型^[19]、DOLY模型^[20]、MAPSS模型^[21]、IBIS模型^[22]和LPJ-DGVM模型^[23]等。其中, HLZ生态系统模型可以依据年平均生物温度、年降水量和潜在蒸散比率3个关键的生物气候要素, 实现对不同植被生态系统类型的定量判别和时空模拟^[24-26]。由于其模型参数相对于其它模型更为简单, 基础数据容易获取, 被国内外学者广泛应用^[27-30]。

青藏高原被称之为“世界屋脊”和“地球第三极”, 是亚洲乃至北半球气候变化的“感应器”和“敏感区”^[31]。其独特而又复杂的地形地貌特征和自然环境格局^[32], 为不同生物区系的相互交汇与融合提供了特定的空间, 使青藏高原成为全球生物多样性保护的25个热点地区之一^{[31, 33-34]}。自20世纪以来, 国内外学者对青藏高原的地表过程开展了大量的持续性研究, 且在高原气候变化、冰川冻土变化、生态系统碳过程、土地覆被变化、生态退化区防治等方面取得了系列研究进展^[35-37]。在青藏高原植被生态系统方面, 现有研究主要集中于气候植被模型的研究进展^[38], 青藏高原局部区域的遥感植被参数提取分析^[39]、植被覆盖度的定量估算^[40]定量估算、以及年际植被变化规律^[41]等方面。

以上研究主要是对过去趋势以及青藏高原的局部区域的模拟分析, 缺乏对整个区域内的植被生态系统在气候变化驱动下的未来时空分布变化特征模拟的相关研究, 尤其很少从垂直分布的角度对整个青藏高原植被生态系统垂直分布的定量解析和模拟预测。因此, 该论文旨在基于Holdridge life zone(HLZ)模型, 结合数字高程模型(DEM)数据, 改变模型输入参数模式, 发展能够从垂直分布角度定量刻画和模拟青藏高原植被生态系统垂直分布变化的改进型HLZ生态系统模型, 进而结合1981—2010 (T0)时段的气候观测数据, 对青藏高原植被生态系统的空间分布格局进行模拟, 在采用CMIP5 RCP2.6、RCP4.5和RCP8.5三种情景数据, 从植被分布垂直地带性的角度, 对青藏高原2011—2040 (T1)、2041—2070 (T2)和2071—2100 (T3)三个未来时段的植被生态系统垂直分布的时空变化情景进行模拟。另外, 引入生态系统时空分析模型^[25]和生态多样性指数模型^[18], 定量求算青藏高原植被生态系统平均中心在未来不同气候情景下的时空演替趋势, 及其生态多样性变化的未来情景。

1 数据与方法 1.1 基础数据收集与处理

气候数据包括观测数据和模式模拟情景数据。气候观测数据来源于1981—2010年位于青藏高原及周边区域的气象台站。气候情景数据采用IPCC CMIP5发布的能代表未来温室气体排放的高中低三种情景^[42], 即：RCP2.6(低排放情景)、RCP4.5(中间排放情景)、RCP8.5(高排放情景) (http://www.ipcc-data.org)。青藏高原的DEM数据采用SRTM数据, 数据来源于http://srtm.csi.cgiar.org, 空间分辨率为1km×1 km。

如何对气象站点的观测数据进行空间插值, 以及对气候情景数据进行空间降尺度, 进而获取高质量的气候要素空间数据, 对于植被生态系统分布特征模拟结果可靠性具有直接影响^{[18, 43]}。鉴于高精度曲面建模(HASM)方法^[44], 能够克服反距离加权模型(IDW)、三角网模型(TIN)、克里金模型(Kriging)、样条插值模型(Spline)等常用方法的理论缺陷并提升模型的模拟精度^[44-45]。因此, 在进行气候观测数据空间插值和未来气候情景数据空间降尺度的过程中, 采用HASM方法, 并结合经纬度及高程数据, 实现气候观测数据的空间插值和未来气候情景数据的空间降尺度^[45-46], 分别获得青藏高原1km×1km空间分辨率的T0、T1、T2和T3四个时段的年平均生物温度、年降水量和潜在蒸散比率数据。

1.2 植被生态系统垂直分布空间模拟模型

针对青藏高原的地形特征和植被分布的海拔垂直差异, 在对HLZ(Holdridge life zone) 生态系统模型^{[26-28, 47]}进行修正和拓展的基础上, 根据青藏高原的年平均生物温度、年降水量和潜在蒸散比率三个主要的生物气候要素的高精度空间模拟数据, 建立了青藏高原植被生态系统垂直分布的空间模拟模型, 其理论计算公式可表征为：

(1)

(2)

(3)

(4)

式中, MAB_Tibet(x, y, t)、TAP_Tibet(x, y, t)和PER_Tibet(x, y, t)分别为青藏高原t时刻位置(x, y)处的平均生物温度(℃)、年降水量(mm)和潜在蒸散比率；TEM_Tibet(j, x, y, t)和P_Tibet(j, x, y, t)分别为青藏高原第j天的积温和降水量；M(j, x, y, t)=log₂MAB_Tibet(x, y, t), T(j, x, y, t)=log₂MAB_Tibet(x, y, t), P(j, x, y, t)=log₂MAB_Tibet(x, y, t); M_i0、T_i0和P_i0分别是青藏高原第i种植被生态系统类型中心点的MAB、TAP和PER的判别标准值; HLZ_{Tibet, i}(x, y, t)为t时刻位置(x, y)处青藏高原第i种植被生态系统类型。

1.3 植被生态系统的斑块连通性和生态多样性指数

斑块连通性是一个可有效计算某一生态系统类型中的动植物扩散传播的平均效率的生态景观指数^[48]。因此, 引入斑块连通性指数模型对青藏高原植被生态系统景观单元的斑块连通性变化进行定量求算, 可以有效揭示在景观层次上青藏高原植被生态系统垂直分布对气候变化的响应情况, 理论公式可表达如下：

(6)

(7)

式中, t为时间变量；S_{i, j}(t), 为物种在斑块(i, j)内的迁移效率指数, A_{i, j}(t)和Pr_{i, j}(t), 分别为第i种植被生态系统类型第j个斑块的面积和周长, 为正六边形周长的平方与面积的比率；P_{i, j}(t)为第i种植被生态系统类型第j个斑块面积占整个研究区域的面积比率；Co(t)为t时段该生态类型的斑块连通性, 且0≤Co(t)≤1.1, 当所有斑块均为正六边形时, Co(t)=1.0。

在求解青藏高原植被生态系统斑块连通性指数的基础上, 引入度生态多样性指数模型^{[26, 48]}, 并对模型输入参数对象进行拓展, 即以青藏高原植被生态系统垂直分布的空间分布作为研究对象, 实现青藏高原植被生态系统多样性的定量计算, 进而揭示青藏高原植被生态系统类型丰富度对气候变化的响应机理。模型的理论公式可表达为

(5)

式中, t为时间变量；p_i(t)为第i种植被生态系统类型的概率；m(ε)为青藏高原植被生态系统类型的数量；, a研究区总面积, e取常数值2.71828；Div(t)即为t时段青藏高原植被生态系统的生态多样性。

1.4 植被生态系统平均中心偏移的时空分析模型

植被生态系统平均中心时空偏移分析模型, 是一个在求算植被生态系统平均中心分布的基础上, 根据其平均中心的时空偏移距离和方向, 对植被生态系统时空分布格局变化进行定量阐述的空间分析模型^[49]。通过植被生态系统平均中心偏移距离的大小可以有效揭示其对气候变化响应的强度, 而偏移方向的变化则可以反映其对气候变化影响的变化频次^[26], 理论公式可表达为:

(8)

(9)

式中, t为时间变量；I_j(t)第j种植被生态系统的斑块数；S_ij(t)为第j种植被生态系统第i个斑块的面积；S_j(t)为第j种植被生态系统的总面积；X_ij(t), Y_ij(t)分别为第j种植被生态系统第i个斑块几何中心的经度和纬度的坐标值；x_j(t), y_j(t)分别为第j种植被生态系统的平均中心经度和纬度的坐标值。第j种植被生态系统平均中心的偏移距离和方向可分别为:

(10)

(11)

式中, d_j为第j种植被生态系统平均中心从t到t+1时段的偏移距离；θ_j为第j种植被生态系统平均中心从t到t+1时段的偏移角度, (x_j(t), y_j(t))和(x_j(t+1), y_j(t+1))分别代表第j种植被生态系统t时段和t+1时段的平均中心坐标；当345° < θ_j≤15°时, 近似认为第j种植被生态系统的平均中心从t时段到t+1时段向东偏移；当15° < θ_j≤75°时, 近似认为第j种植被生态系统的平均中心从t时段到t+1时段向东北方向偏移；当75° < θ_j≤105°时, 近似认为第j种植被生态系统的平均中心从t时段到t+1时段向北偏移；当105° < θ_j≤165°时, 近似认为第j种植被生态系统的平均中心从t时段到t+1时段向西北方向偏移；当165° < θ_j≤195°时, 近似认为第j种植被生态系统的平均中心从t时段到t+1时段向西偏移；当195° < θ_j≤255°时, 近似认为第j种植被生态系统的平均中心从t时段到t+1时段向西南方向偏移；当255° < θ_j≤285°时, 近似认为第j种植被生态系统的平均中心从t时段到t+1时段向南偏移；当285° < θ_j≤345°时, 近似认为第j种植被生态系统的平均中心从t时段到t+1时段向东南方向偏移。

2 模拟结果 2.1 青藏高原植被生态系统的空间分布格局

运用1981—2010(T0)时段的气候观测数据对青藏高原植被生态系统垂直分布格局的模拟结果表明(图 1), 青藏高原共有17种植被生态系统类型, 其中, 低山有刺草原在RCP2.6和RCP4.5情景中均未出现, 仅在RCP8.5情景的T3时段出现, 其于16种植被生态系统类型垂直分布的总体格局为：荒漠生态系统类型主要分布在塔里木盆地的塔克拉玛干沙漠地带及其周围地区。从青藏高原的山前地带到高海拔地区依次为山前湿润森林、低山湿润森林、低山干旱生森林、低山荒漠灌丛、山地湿润森林、山地灌丛、山地草原、亚高山潮湿森林、亚高山湿润森林、亚高山干旱灌丛、高山雨苔原、高山潮湿苔原、高山湿润苔原、高山干苔原和冰雪/冰原等14种植被生态系统类型。其中, 高山潮湿苔原、亚高山湿润森林和冰雪/冰原是青藏高原主要的植被生态系统类型, 分别占总面积的23.69%、17.79%和14.78%。高山潮湿苔原主要分布在青藏高原的中部、西南部和东北部区域, 亚高山湿润森林主要分布在青藏高原的南部和东北部区域, 而冰雪/冰原则主要分布在青藏高原的中部、西北、西南和东北部的高海拔区域。

2.2 青藏高原植被生态系统面积的变化情景

对RCP2.6、RCP4.5和RCP8.5三种情景下的青藏高原植被生态系统模拟结果(图 2-4)进行统计分析表明, 在三种不同强度的气候变化情景下, 青藏高原植被生态系统类型的面积在2010—2100年间将呈现出不同的变化趋势和空间差异。

在RCP2.6情景下, 高山湿润苔原、荒漠和高山干苔原的面积在T0—T3将呈持续减少趋势, 平均每10年分别减少1.451×10⁴km²、1.343×10⁴km²和0.151×10⁴km², 其中高山干苔原和荒漠的面积的缩减幅度将分别高达92%和69%。亚高山潮湿森林、山地湿润森林和山地灌丛的面积将呈持续增加趋势, 平均每10年分别增加1.831×10⁴km²、1.286×10⁴km²和0.903×10⁴km², 其中山地灌丛和亚高山潮湿森林面积增加的幅度将分别接近170%和105%。此结果表明, 高山干苔原、荒漠、山地灌丛、亚高山潮湿森林等植被生态系统类型对气候变化的敏感性高于其它的植被生态系统类型。

在RCP4.5情景下, 冰雪/冰原、高山湿润苔原和荒漠三种植被生态系统类型面积在T0—T3时段将呈持续减少趋势, 平均每10年分别减少3.445×10⁴km²、1.920×10⁴km²和1.541×10⁴km²。亚高山潮湿森林、山地湿润森林和山地灌丛三种植被生态系统类型的面积增加最多, 在T0—T3时段平均每10年将分别增加3.266×10⁴km²、2.485×10⁴km²和1.240×10⁴km²。从T0到T3时段, 高山干苔原和荒漠的面积减少幅度分别为95%和79%, 而亚高山潮湿森林、山地灌丛、山地湿润森林和低山湿润森林面积的增加幅度将超过130%。此结果表明, 在RCP4.5情景驱动下的高山干苔原、荒漠、山地灌丛、亚高山潮湿森林、低山湿润森林、山地湿润森林等植被生态系统类型对气候变化的敏感性高于其它的植被生态系统类型。

在RCP8.5情景下, 冰雪/冰原、高山湿润苔原和荒漠三种植被生态系统类型的面积持续减少最多, 在T0—T3时段平均每10年将分别减少4.034×10⁴km²、2.498×10⁴km²和1.866×10⁴km²。山地湿润森林、亚高山潮湿森林和山地灌丛三种植被生态系统类型面积的面积持续增加最多, 在T0—T3时段平均每10年将分别增加5.179×10⁴km²、5.153×10⁴km²和1.430×10⁴km²。从面积增加比例来看, 在T0—T3时段内, 亚高山潮湿森林、山地灌丛、山地湿润森林、低山干旱森林、低山湿润森林面积的增长比率均大于290%, 表明这几种植被生态系统类型对气候变化响应的敏感性高于其它的植被生态系统类型。

2.3 植被生态系统类型多样性及斑块连通性分析

对青藏高原植被生态系统的生态多样性及斑块连通性的模拟结果(表 1)表明：在未来气候变化驱动下, 青藏高原植被生态系的生态多样性整体上呈减少趋势, 斑块连通性则呈持续减少趋势。其中, RCP8.5情景下的生态多样性呈持续减少趋势, 在T0—T3时段内平均每10年将减少0.227%；RCP2.6和RCP4.5情景下的生态多样性呈波动减少趋势, 在T0—T3时段内平均每10年将分别减少0.056%和0.041%；RCP2.6、RCP4.5和RCP8.5三种情景下的青藏高原植被生态系统斑块连通性均呈持续减少趋势, 其平均每10年将分别减少0.225%、0.226%和0.419%。这一研究结果显示, 未来气候变化将导致青藏高原植被生态系统的生态多样性和斑块连通性的整体呈减少趋势, 且气候变化越快(RCP8.5), 其生态多样性和斑块连通性减少速度越快。这意味着未来气候变化强度将直接影响青藏高原植被生态系统生态多样性和斑块连通性的变化趋势和强度。

2.4 青藏高原植被生态系统平均中心的时空偏移情景

运用生态系统平均中心时空分析模型对青藏高原植被生态系统平均中心的时空偏移情景进行模拟的分析结果(表 2-4, 图 5)表明, 在RCP2.6、RCP4.5和RCP8.5三种情景的驱动下, 青藏高原不同植被生态系统的平均中心将呈现出不同的时空偏移趋势。

在RCP2.6情景下(表 2, 图 5), 低山荒漠灌丛的平均中心在T1—T2时段将向东偏移0.158 km, 然后在T2—T3时段, 其平均中心将转向南偏移0.542 km。高山干苔原平均中心的偏移幅度最大, 其在T0—T1时段将向东北为偏移4.974 km, 在T1—T2时段转向西南偏移2.199 km, 在T2—T3时段继续向西南偏移2.351 km。山地湿润森林的平均中心在T0—T3时段内将持续向西南方向偏移。冰雪/冰原、高山湿润苔原、高山潮湿苔原、亚高山湿润森林和亚高山潮湿森林等植被生态系统类型的平均中心在T0—T3时段内将整体呈向西南偏移趋势。低山干旱森林的平均中心偏移距离最大, 在T0—T3时段内平均每10年将偏移1.33km。低山湿润森林的平均中心偏移距离最小, 在T0—T3时段内平均每10年仅偏移0.05km。在RCP 2.6情景下, 青藏高原所有植被生态系统类型的平均中心在T0—T3时段内平均每10年的偏移距为0.38km。

在RCP4.5情景下(表 3, 图 5), 低山荒漠灌丛植被生态系统仅出现在T1—T3段出现, 其平均中心将向持续向北分别偏移1.053 km和0.617 km。亚高山干旱灌丛、亚高山潮湿森林以和山地草原三种植被生态系统的平均中心偏移距离大于其它的生态系统类型。高山湿润苔原和高山雨苔原的平均中心在T0—T3时段将持续向南偏移, 而山地湿润森林和荒漠生态系统的平均中心在T0—T3时段将持续向西南方向偏移。冰雪/冰原、亚高山湿润森林和山地湿润森林的平均中心在T0—T3时段将整体向西南方向偏移。低山湿润森林的平均中心偏移距离最小, 在T0—T1时段内向西南偏移0.329 km, 继而在T1—T2时段转向东北方向偏移0.286 km, 在T2—T3时段则继续向东北方向偏移0.138km。

在RCP8.5情景下(表 4, 图 5), 低山荒漠灌丛从T1时段开始出现, 其平均中心在T1—T3时段内先向西北偏移1.197 km, 然后转向西南偏移0.713 km。青藏高原高山湿润苔原、山地草原和荒漠生态系统类型的平均中心偏移强度高于其它的植被生态系统类型, 在T1—T3时段内平均每10年将分别偏移1.012km、1.361km和1.724km。山地湿润森林和荒漠生态系统类型的平均中心在T1—T3时段内将持续向西南方向偏移。冰雪/冰原、高山干苔原、高山苔原、亚高山湿润森林和山地灌丛等植被生态系统类型的平均中心将整体上向西南方向偏移。低山湿润森林的平均中心偏移强度最小, 在T1—T3时段内平均每10年仅偏移0.105km。

2.5 青藏高原不同海拔植被生态系统变化情景

为了更好地分析不同海拔高度上青藏高原植被生态系统在不同气候情景驱动下的变化特征, 结合青藏高原植被生态系统情景模拟结果, 以1000m海拔高程差为间隔, 分别对RCP2.6、RCP4.5和RCP8.5三种情景下的青藏高原不同海拔上植被生态系统变化情况进行空间统计和对比分析(表 5)显示: 在青藏高原所有的海拔高度带上, RCP8.5情景下的植被生态系统的变化强度最大, 在RCP4.5情景下次之, RCP2.6情景下变化强度最小。在RCP2.6、RCP4.5和RCP8.5三种情景下, 青藏高原海拔大于8000 m的植被生态系统变化强度整体上均高于其它海拔高度带上的植被生态系统变化强度, 在T0—T3时段内三种情景下每10年的变化面积分别占该梯度总面积的7.41%、7.41%和18.52%。另外, 在RCP2.6和RCP4.5两种情景下分布于青藏高原海拔小于1000m的植被生态系统变化最为缓慢, 在T0—T3时段内三种情景下每10年的变化面积分别占该梯度总面积的1.64%和5.08%。在RCP45情景下, 分布于青藏高原海拔1000—2000m之间的植被生态系统变化最为缓慢, 在T0—T3时段内每10年的变化面积占该梯度总面积的2.68%。总之, 在三种情景下, 青藏高原的植被生态系统变化强度从低海拔到高海拔区域整体呈上升趋势。

2.6 不同气候情景驱动下青藏高原植被生态系统变化的对比分析

在RCP2.6、RCP4.5和RCP8.5三种气候情景驱动下的青藏高原植被生态系统将呈现处不同的差异特征。整体上表现为, 青藏高原植被生态系统在RCP8.5情景驱动下的面积变化、平均中心偏移、生态多样性变化及不同海拔高度分布的变化强度整体上高于在RCP2.6和RCP4.5情景驱动下的变化情景。其中, 冰雪/冰原和荒漠两类生态系统类型的面积在三种情景下的T0—T3时段内均呈持续减少趋势, 且RCP8.5情景下的减少速度高于RCP2.6和RCP4.5情景下的减少速度；三种情景下的青藏高原植被生态系统的生态多样性和斑块连通性整体上呈减少趋势, RCP2.6情景下的减少速度最慢, 而RCP8.5情景下的减少速度最快；三种情景下的冰雪/冰原、亚高山湿润森林和山地湿润森林的平均中心在T0—T3时段均将向西南方向偏移, 且在RCP8.5情景下的偏移幅度最大, RCP4.5情景次之, RCP2.6情景下的偏移幅度最小；从青藏高原的低海拔到高海拔区域, 分布于不同海拔高度上的植被生态系统类型的变化强度呈逐渐提升趋势。

3 讨论与结论

该论文针对现有HLZ生态系统模型主要运用于从维度和经度二维的角度模拟和解释植被生态系统在不同尺度上的变化情景^{[17-18, 26-28, 30]}, 而缺乏从垂直角度模拟植被生态系统情景变化的应用不足, 结合气候观测数据和气候情景数据, 从植被生态系统分布垂直地带性的角度, 实现了RCP2.6、RCP4.5和RCP8.5三种气候情景驱动下的青藏高原植被生态系统垂直分布的时空变化情景模拟。结果显示, 有效集成DEM数据的改进型HLZ生态系统模型, 能够有效地对不同气候情景驱动下的青藏高原植被生态系统垂直分布带的时空分布格局变化及未来情景进行模拟。

模拟结果表明, 青藏高原共有16种植被生态系统类型, 其中高山潮湿苔原、亚高山湿润森林和冰雪/冰原为青藏高原主要的植被生态系统类型, 分别占总面积的23.69%、17.79%和14.78%。在T0—T3期间, 青藏高原的高山湿润苔原、高山干苔原、荒漠呈持续减少趋势, 平均每10年将分别减少1.96×10⁴ km²、0.15×10⁴ km²和1.58×10⁴ km²；亚高山潮湿森林、山地湿润森林和山地灌丛呈持续增加趋势, 平均每10年将分别增加3.42×10⁴ km²、2.98×10⁴ km²和1.19×10⁴ km²；RCP8.5情景下青藏高原的植被生态系统平均中心的偏移幅度最大, RCP4.5情景下的偏移幅度次之, 而RCP2.6情景下的偏移幅度最小。另外, 在三种气候变化情景驱动下, 青藏高原植被生态系统的生态多样性呈减少趋势。青藏高原区域内的高山干苔原、亚高山潮湿森林、山地灌丛、山地湿润森林和荒漠等植被生态系统类型对气候变化的敏感性高于其它的植被生态系统类。总之, 气候变化强度将直接影响青藏高原植被生态系统的面积、平均中心及其生态多样性的时空变化幅度, 在T0—T3时段内, RCP8.5情景下的变化幅度最大, 而RCP2.6情景下的变化幅度最小。气温和降水作为植被生态系统分布及其多样性的直接驱动因子, 在未来气候变化情景下, 青藏高原植被生态系统在不同的垂直带上将呈现处不同的差异。随着海拔的升高, 青藏高原植被生态系统空间分布的动态差异对气候变化的响应强度逐渐增大^{[25, 49]}。尤其是分布于青藏高原海拔8000m以上的植被生态系统类型对气候变化的响应最为敏感。

因此, 在青藏高原生态系统时空变化趋势及未来情景的综合评估, 及其对自然气候变化与人类活动耦合驱动效应研究的过程中, 不仅要关注低海拔地区及荒漠化区域的生态系统效应分析, 而且需要更加关注高海拔区域的生态系统的时空演化机理分析, 这对深入揭示气候变化对青藏高原生态系统的影响效应及驱动机理, 以及如何部署青藏高原生态安全动态监测体系, 提升国家生态安全屏障的总体功能具有重要意义^[32]。