青藏高原不同退化梯度下植被蒸散发的时空格局研究

刘长雨，谢保鹏*，杨 洁，陈 英，裴婷婷

(1.甘肃农业大学管理学院，甘肃 兰州 730070；2.甘肃农业大学草业学院，甘肃 兰州 730070)

全球气候变暖和人类活动的共同影响下，地球正以一定的速度绿化[1-2]。美国航天局卫星数据表明，全球从2000—2017 年新增的绿化面积中，约1/4来自中国，贡献比例居全球首位[3]。青藏高原作为我国乃至亚洲的重要生态安全屏障[4]，其独特的地理位置和典型的高寒生态系统使其成为国内外学者研究的热点区域。过去的50年里，青藏高原气候变化的暖湿化特征明显[5]，温度每10 年升高0.3℃～0.4℃，降水每10 年增加2.2%[6]。降水增加和温度上升导致蒸散发(Evapotranspiration,ET)加剧，土壤水分严重下降并直接影响植被生长[1,7]，同时，过度的蒸散发也会影响大气循环、造成温室效应，甚至出现极端气候[7]。因此，准确识别青藏高原的蒸散发与植被变化的关系对陆地生态系统的稳定有十分重要的意义。

植被作为蒸散发的重要媒介，对气候变化的响应明显。近年来，青藏高原的植被呈现出退化趋势,尤其是高原腹地植被退化趋势明显[8]，因此,识别青藏高原植被退化梯度对蒸散发的研究有重要意义。当前国内外对青藏高原蒸散发的研究已取得较多成果，可以概括为以下几个方面：一是对蒸散发模型的研发与应用[9-10]，如张文旭等[11]分析了涡度相关法、蒸渗仪测定和FAO56Penman-Monteith公式3种方法计算青藏高原风火山的蒸散发，结果表明3种方法计算的蒸散发结果相似；二是通过站点实测数据或遥感数据探究青藏高原蒸散发的时空分布格局[12-14]，如姚天次等[15]利用274个气象站数据探究了1970—2017 年高原及周边地区潜在蒸散发的空间格局，得出年潜在蒸散发呈现南北高、中部低的空间分布的特征；三是刻画气候变化的不同因子对蒸散发的影响研究[16-18]，如Ma等[19]探究了青藏高原1982—2016 年的蒸散发及影响因素，研究结果表明降水增强是蒸散发加剧的主要驱动因素。这些研究均从计算蒸散发的方法和分析蒸散发的时间格局及其影响因子出发，较少关系植被退化与蒸散发的关系。因此，本文从植被退化出发，了解高寒地区植被退化与蒸散发的交互过程，这对当地生态安全和高寒植被退化的恢复治理有重要意义。

基于此，本研究选定青藏高原这一高寒干旱区，利用MODIS16A2数据，从时空尺度、空间尺度以及不同植被类型上分析蒸散发的变化特征，同时根据像元二分模型计算植被覆盖度，识别青藏高原植被退化梯度和植被变化区，进而探究青藏高原不同退化梯度和不同植被变化区蒸散发的时空变化特征。

1 材料与方法

1.1 研究区概况

青藏高原是地球上最高的高原，地理范围为68°～104° E，26°～40° N，总面积约250×104km2。该地区的地势西高东低，年均气温-6℃～20℃，年降水量50～2 000 mm，是典型的高原高寒气候，地域差异明显，气温低，冬季长，日照充足，日较差大。从东南部向西北部，气候由温暖湿润向寒冷干燥变化，水平带上也先后出现了森林、灌木、草地等植被类型[20]。本文综合考虑我国流域分区、省级行政区划和“一带一路”亚洲关键区域流域边界[21]，将青藏高原划分为14 个子区域，分别是黄河流域、印度河流域、石羊河流域、疏勒河流域、黑河流域、塔里木河流域、阿姆河流域、雅鲁藏布江流域、青海湖水系、羌塘高原区、长江流域、怒江流域、澜沧江流域和柴达木盆地，以期从整体和局部的不同视角对青藏高原蒸散发的时空演变格局进行深入研究(图1)。

图1 青藏高原地理位置示意图

1.2 数据来源与处理

1.2.1MODIS蒸散发数据 MODIS16A2是NASA发布的全球陆地蒸散发产品数据，基于Penman-Monteith公式计算得到，时间分辨率16天，空间分辨率500 m。通过全球通量塔台站的检验，MOD16A2数据拟精度可达86%[22]。MOD16A2数据在时空上具有高分辨率、连续覆盖等特点，在全球得到了较为广泛的应用[23]。目前，MOD16A2数据已用于我国地表蒸散发的时空变化特征分析[24]，总体上具有适用性，在站点资料稀缺的青藏高原区，MOD16A2数据能够较准确地呈现蒸散发的空间分布规律[25]。本研究选用2001—2020 年的MOD16A2产品数据(来源：https://search.earthdata.nasa.gov/)，通过Arc GIS10.2软件的像元统计工具计算年际蒸散发数据。

1.2.2植被数据 本研究使用的中国植被图比例尺为1：1 000 000，由中国科学院资源环境科学与数据中心(https://www.resdc.cn/)提供。该数据将其划分为13 大类，可归纳为6类，分别是森林(针叶林、阔叶林、高山植被)、灌木(灌丛)、草地(草原、草甸、草丛)、荒漠、栽培植被和其他(沼泽和其他)，本文选取了森林、灌木和草地进行研究。

1.2.3植被覆盖度 植被覆盖度是基于像元二分模型进行计算。像元二分模型是一种简单实用的遥感估算模型，它假设一个像元的地表由有植被覆盖部分地表与无植被覆盖部分地表组成，而遥感传感器观测到的光谱信息也由这2个组分因子线性加权合成，各因子的权重是各自的面积在像元中所占的比率，如植被覆盖度可以看作是植被的权重。

(1)

式中，NDVIsoil作为裸土的NDVI值，理论上应该趋近于0，但是由于遥感影像受到大气环境和粗糙地表、土壤颜色、土壤地表湿度等环境的影响，NDVIsoil值不是固定值，会在一定范围内变化，一般为-0.1～0.2。NDVIveg为研究区与像元最大NDVI值。鉴于此，对研究区内像元NDVI的灰度值进行统计分析，截取置信区间累计频率在5%～95%对应的NDVI值分别作为NDVI最大值和最小值，从而计算得到植被覆盖度FVC。

1.3 研究方法

1.3.1青藏高原植被退化遥感监测和评价指标体系建立 根据国家标准GB19377-2003天然植被退化、沙化、盐渍化的分级指标，以1982—1985 年每个像元的年际最大植被覆盖度(AVHRR数据计算得到)作为基准，将其退化程度分为未退化、轻度退化、中度退化、重度退化和极重度退化5个级别[26-27](表1)。

表1 植被退化等级

1.3.2趋势分析法 本文采用Theil-Sen中位数趋势分析和Mann-Kendall检验来研究青藏高原蒸散发的空间趋势特征。Theil-Sen中位数趋势分析是一种稳健的非参数统计的趋势计算方法，对测量误差和离群数据不敏感，可用于ET的时间序列趋势分析。而Mann-Kendall趋势检验作为非参数统计检验的方法，用于评估Theil-Sen斜率估计的显著性，即检验青藏高原ET趋势的显著性[28]。具体公式如下：

式中：β为斜率；i和j代表年份；如果β>0，ET数据集时间序列具有正趋势；如果β<0，ET数据集序列具有负趋势。

(3)

(4)

(5)

式中：xj和xi为ET时间序列数据；n为ET序列中数据个数；sgn()为符号函数；如果|Zc|≥Z1-α/2，则原假设被拒绝，ET序列呈显著变化趋势；如果Zc>0，则表明ET序列呈上升趋势；如果Zc<0，则表明ET序列呈下降趋势；α为给定的置信水平(显著性检验水平)，|Zc|≥1.96时表示ET序列通过了置信度为95%的显著性检验。

1.3.3变异系数 变异系数用来表示地理数据的波动程度[29]，在一定程度上可表明陆地蒸散量的变化情况。其计算公式如下：

(6)

1.3.4Hurst指数 基于重标极差(R/S)分析方法的Hurst指数(H值)是非线性时间序列分析的一种方法。常用于分析长时间序列的分形特征和长期记忆过程，也称为重新标度极差分析[30]。H值从0到1。当H=0.5表明时间序列为随机序列，前后变化无关联；当H>0.5时，未来趋势与前期保持一致(即继续上升或下降)，H越高，持续性越强。当H<0.5时，未来趋势很可能是下一个时期的反持续性(即从上升到下降，反之亦然)，H越低，反持续性越大[3]。

2 结果与分析

2.1 青藏高原ET的时空变化特征

2.1.1青藏高原ET的年际变化特征 如图2A所示，2001—2020 年青藏高原年蒸散发呈波动上升趋势且始终保持在600～800 mm之间，平均值为730.48 mm。20年间，2001年研究区的ET最低，为643.58 mm，最高值出现在2013年，为780.49 mm，其次是2020年的779.83 mm。

青藏高原的植被类型包括森林、灌木和草地，20年里各植被类型的年均ET呈波动上升的趋势，ET大小表现为：灌木>森林>草地(图2B)。2006 年和2013 年，森林ET高于灌木ET，且在2012—2013年3种植被类型的ET均有大幅上升。由图2C可以看出，2001—2020年青藏高原包含的14 个子流域的ET均呈波动上升的趋势，其中黄河流域区ET最高，各年份均在800 mm以上，长江流域次之，均在750 mm以上。羌塘高原区ET最低，年均ET在620 mm以下。四大内陆河流域中石羊河流域ET>黑河流域ET>疏勒河流域ET>塔里木河流域ET，其中石羊河流域的ET在20年里的变化量最大，为376.81 mm。

图2 青藏高原蒸散发(A)、不同植被类型蒸散发(B)以及子流域蒸散发变化图(C)

2.1.2青藏高原ET的空间变化特征 青藏高原年平均蒸散发在空间上呈中部低，东南高的态势(图3A)。据统计，年均ET范围为800～1 000 mm的面积最大，占总面积的31.14%，其次为1 000～1 200 mm，占总面积的30.74%，0～200 mm面积最小，占总面积的0.03%。从空间上来看，位于西藏自治区的拉萨市、日喀则地区、那曲地区以及青海省的西宁市的ET较低，这可能与拉萨、西宁作为省会城市，城市化的发展速度较快导致植被覆盖度降低有关。

2001—2020年青藏高原ET的变化趋势如图3B所示，斜率介于62.9 至54.3 之间，80.02%的研究区域的ET呈增加趋势。青藏高原东部地区的ET变化相对强于西部地区，中部和南部的ET显著下降，如西藏自治区的林芝地区，这可能与气候变化和高山高原植被退化有关。2001—2020年青藏高原ET的变异系数在0.02～1.125之间(图3C)，表明ET的变化具有显著的空间异质性。青藏高原大部分地区ET波动较低，是由于该区域人口相对较少，高山植被覆盖稳定，变异系数高于0.5的区域集中在区域南部，占全域面积的0.29%。青藏高原蒸散发的Hurst指数在0.06～0.904之间(图3D)，平均值为0.408。Hurst指数<0.5的区域占总面积的86.52%，这表明ET在青藏高原的未来趋势中具有显著的反持续性，其中那曲地区的东南部、海西蒙古族自治区和玉树藏族自治区的交接处ET的未来变化趋势与当前变化趋势相反，即未来将呈增加趋势。

青藏高原各植被类型存在显著的空间异质性，草地面积最大且广泛分布于研究区的中部和西部，而森林和灌木呈片状或带状集中分布在青藏高原的东南部(图3E)。2001—2020年，森林、灌木和草地年均ET分别为802.00 mm,832.21 mm,705.68 mm。从图3F可以看出，青藏高原子流域的年均ET在空间上呈“西北低，东南高”的态势，这与本研究区的海拔呈显著负相关关系，即海拔高的地区ET较低。青藏高原14个子流域中羌塘高原区的年均ET最低，为515.15 mm，这可能与羌塘高原区位于青藏高原腹地，海拔高且四周高山重围、气候寒冷干燥，植被难以生长有关。黄河流域的年均ET最高，为909.69 mm，其次为长江流域，为834.13 mm。

2.2 青藏高原植被退化时空变化特征

2.2.1青藏高原植被退化时空变化格局 青藏高原植被退化自西北向东南呈先减轻后加剧的态势(图4A)，极重度退化的区域集中在青藏高原西北部，如柴达木盆地和昆仑山脉周边。植被未退化及轻度退化区域普遍位于黄河流域上游、青海湖水系以及长江流域、澜沧江流域的部分地区。同时，青藏高原东南部的森林地区也出现一定程度的退化。

从时间上来看(图4B)，2001—2020年青藏高原各植被退化区面积一直保持着：极重度退化区>重度退化区>未退化区>轻度退化区>重度退化区的态势，其中轻度退化和重度退化在20年里的变化不大，且面积基本持平。20年间，极重度退化区面积呈下降趋势，减少了9.88×104km2，重度退化区面积呈波动上升趋势，增加了3.05×104km2，未退化区面积总体呈大幅度下降继而显著上升趋势，共增加了10.91×104km2，其中降幅最大的时间在2013年，降低了11.16×104km2，增幅最大在2020年，增加了7.71×104km2。这说明20年间青藏高原植被退化程度高的区域在不断缩小，未退化区域不断增加且增幅在不断增强。

图4 青藏高原植被退化时空变化图

2.2.2青藏高原植被退化转移矩阵 表2可知，青藏高原2001年未退化、轻度退化、中度退化、重度退化、极重度退化面积占比分别为：23.46%,10.37%,9.97%,27.33%,28.88%；2020年面积占比分别为：28.32%,10.60%,10.26%,27.41%,23.41%。2001年极重度退化区域面积最大，2020年则为未退化区域面积最大，表明青藏高原的植被退化在不断减弱。

表2 青藏高原植被退化面积占比转移矩阵

2.3 青藏高原不同退化梯度下ET的时空变化特征

2.3.1青藏高原不同退化梯度内ET变化特征 青藏高原2001—2020年各退化梯度下的ET总体呈波动上升趋势，其中各退化区ET的均值表现为：未退化区>轻度退化区>极重度退化区>中度退化区>重度退化区，分别是869.89 mm,738.77 mm,731.14 mm,684.40 mm,610.83 mm，这说明植被退化在最低程度或最高程度时蒸散发较高(图5)。

五种退化梯度中未退化区ET的增幅最大，年均斜率为6.56，其次为极重度退化、轻度退化、中度退化和重度退化，其斜率分别为：4.74,4.37,3.84,2.92，同时各退化梯度ET均在2012年大幅下降后在2013年显著上升。

2.3.2青藏高原植被变化区内ET的时空变化特征 图6A可以看出，青藏高原植被不变区、恢复区以及恶化区面积分别是：126.77×104km2,52.27×104km2,26.32×104km2。在空间上，植被不变区广泛分布于青藏高原全域，占全域面积的61.7%，植被恢复区集中在青海湖周边、柴达木盆地外围以及黄河流域的西北部地区，植被恶化区则呈点状或片状分布于青藏高原的西部和南部。识别出青藏高原20年间植被变化区对该地区生态修复政策制定和推广具有一定意义。

2001—2020年青藏高原各退化区的ET整体呈波动上升趋势(图6B)，且植被恢复区ET>植被不变区ET>植被恶化区ET，3个退化区均在2012 年大幅下跌后在2013 年急剧上升。20年间，植被恢复区的ET增量最多，共增加了184.43 mm，植被不变区和植被恶化区的ET分别增加了118.69 mm和118.03 mm。20年间，青藏高原植被恢复相较于植被恶化更能促使地区蒸散发的增加。

图6 青藏高原植被变化区的蒸散发时空变化图

3 讨论

在全球气候变暖和绿化趋势增强的背景下，蒸散发的增强是主要趋势[31-32]。Pascolini-Campbell等[1]证明全球蒸散发呈显著的正线性增长趋势，增长的部分来源于降水量的增加和径流的减少，Zeng等[33]证明1982—2011 年，全球蒸散发呈现出显著的增长趋势，其中超过一半的增长归因于地球的绿化。气象因素对蒸散发也具有显著的驱动作用，如降水量被认为是影响中国ET变化的主导因素，风速、温度、日照时数和相对湿度等均对蒸散发有影响[7]。本研究表明，2001—2020 年青藏高原蒸散发呈波动上升趋势蒸散发的增强会导致干旱半干旱地区更多的水分流失[1]，从而抑制植物生长，破坏生态系统水分循环，最终可能影响整个陆地生态系统的稳定性。此外，青藏高原年均蒸散发在全域尺度、子流域尺度以及不同植被类型中均表现为2012 年大幅下降后2013 年大幅上升，这种现象可能与该年份出现比正常年份更高的温度、降水以及更低的相对湿度等，也可能与我国生态恢复政策有关[34-35]。裴婷婷等[36]证明2012 年以后全国的降水明显升高，同时降水引起的蒸散发也大幅度增加，这在一定程度上可以解释2012 年以后蒸散发的显著升高。青藏高原草地面积占比最大，其次是森林和灌木，其中灌木的年均ET和森林的年均ET相差不多，但略大于森林。张永强等[37]在探究2003—2017年植被变化对全球陆面蒸散发的影响时，得出2012 年以后植被变化对蒸散发的影响显著加强，灌木和耕地对蒸散发的影响最为显著，这与本文的研究结果较为一致。青藏高原14 个子流域中陆地蒸散量均处于波动上升的趋势，这与Li[38]对子流域的ET不断升高的研究结果相同。因此,本文的研究结果对认识青藏高原不同空间维度的蒸散发格局提供了参考。

2001—2020 年，青藏高原植被退化在不断减弱，这与部分学者的研究结果相似。如夏龙等[39]在对青藏高原草地退化的研究中证明青藏高原草地整体上在朝着改善的方向发展。Wang等[40]在探究青藏高原植被变化的主导因素中得出，青藏高原植被在2000—2015 年有改善趋势。同时，全球植被绿化程度的不断增强也能从侧面证明植被退化在一定程度上减轻[2]。本文的研究结果表明蒸散发在未退化梯度上最高，其次为极重度退化梯度，这可能是由于蒸散发的途径包括植被蒸腾、植被冠层截留降水以及露水蒸发，从而使得植被覆盖度高的地区蒸散发较高。极重度退化区的蒸散发也较高可能与降水和地表径流直接被划分为蒸散发有关。

青藏高原的青海湖和柴达木盆地周边地区的植被在20 年里呈逐步恢复态势且恢复程度显著，这可能与我国出台的生态修复政策有关，如针对青海湖和柴达木周边地区的生态问题，集中采用退耕还草、轮牧、灌草结合、围栏封育等措施来促进植被恢复[41]。柴达木盆地周边地区的主要植被类型为草地，1980—2018 年草地面积呈增加趋势，其中高中覆盖度草地面积增加趋势明显，说明柴达木盆地植被退化程度在逐渐减轻[42-43]。刘秋漫[44]对柴达木盆地的研究表明，近30 年间该地区生态环境质量整体上呈现出先下降后上升的趋势，生态环境质量总体上呈现出向好的方向发展的态势，这与本文的研究结果相似。本研究仍存在一些不足，首先在产品选择上，MODIS产品数据集的反演模型在全球通用，但在研究区的适用性有待提高，同时产品数据的精度相对不高，研究结果存在不确定性，其次本文采用中国科学院资源环境科学与数据中心的植被分类体系将植被划分为林地、草地和灌木，该分类体系会对研究结果产生影响。

4 结论

本文基于2001—2020年MODIS16A2数据，探究了青藏高原蒸散发的时空变化格局同时识别了青藏高原植被退化梯度，进而探究五种退化梯度下青藏高原蒸散发的变化特征。

2001—2020年青藏高原年均ET呈波动上升的趋势，其中2001年ET最低，为643.58 mm，2013年最高，为780.49 mm。从空间上来看年均ET呈中部低，东南高的态势，低值集中在羌塘高原区且86.52%的区域未来趋势具有显著的反持续性，即青藏高原大部分区域的蒸散发将呈现下降趋势。

青藏高原各自然植被类型ET中：灌木>森林>草地，且在2012年大幅下降继而在2013年大幅上升。青藏高原范围内的子流域ET均呈波动上升趋势，其中黄河流域的ET最大，羌塘高原区的ET最小，四大内陆河流域中石羊河流域ET>黑河流域ET>疏勒河流域ET>塔里木河流域ET。

2001—2020年青藏高原植被退化不断减弱，在空间上呈自西向东呈退化先减轻后逐步加强的趋势，各退化区面积一直保持着：极重度退化区>重度退化区>未退化区>轻度退化区>重度退化区。青藏高原各植被变化区的面积为：不变区>恢复区>恶化区，在空间上表现为不变区广泛分布于青藏高原内部，恢复区集中分布在青藏高原东北部，如青海湖和环柴达木盆地周边区域，恶化区呈点片状分布于青藏高原西南部。20年里各植被变化区的年均ET呈波动上升态势，表现为：恢复区>不变区>恶化区。