植被作为地表的重要覆盖类型, 是陆地生态系统的核心组成部分。植被不仅在生物圈、大气圈和全球碳循环中扮演着重要的角色, 而且也是连接土壤、水分和能量等自然要素的纽带^[1]。各种地表植被之间的差异可以反映出不同的水热组合, 所以监测植被之间的差异, 在地方、区域和全球范围内, 对生态系统的状态, 环境压力和景观变化的研究具有重要指示作用^[2]。草地植被是全球陆地生态系统中面积最大的可更新资源, 一方面, 草地植被对培育土壤肥力、防止水土流失、维持陆地生态系统平衡起着重要作用, 另一方面草地植被也可以为食草动物提供饲料, 为人类生产食物、药品和工业原料等^[3-4]。但是, 草地生态系统生境脆弱, 比较容易受到外界环境的干扰, 因而对全球气候环境变化十分敏感^[5-6]。

青藏高原总面积约为250万km², 约占中国陆地总面积的四分之一, 平均海拔在4000 m以上, 是世界上海拔最高的高原, 常被称为世界的“第三极”, 高海拔、地处中低纬度的地理条件以及特殊的下垫面条件也使得青藏高原成为全球对气候变化最为敏感的区域之一。近年来, 在全球气候变化和人为因素的共同作用下, 青藏高原生态环境发生了一系列的负面变化, 例如湿地和草地面积减小、冻土退化、土地沙漠化严重、水域减少、生物多样性减小及自然灾害增多等^[7-10]。青藏高原地区天然草地面积约为1.5×10⁸ hm², 约占全国草地总面积的三分之一。青藏高原地区特殊的地理环境加上草地生态系统本身易受到外界环境干扰的特点, 使得青藏高原草地区域具有典型的生态脆弱性。其脆弱性不仅表现在该地区的草地植被自我更新能力差, 草地的自我恢复能力和抵抗外界来自于人类或自然的扰动能力低, 而且当草地生态系统遭到破坏后, 在重建、保护、改善、治理区域生态环境的过程中所付出的经济成本和生态成本也明显高于低海拔地区^[11]。因此, 有必要对青藏高原地区草地植被进行长时间序列的动态监测。

传统的草地监测方法主要采用野外布设样方进行采样的方法, 这种方法不仅效率低、成本高, 而且青藏高原地区海拔高、部分地区生存环境恶劣, 人力往往很难到达, 所以无法快速、大范围对草地植被生长状况进行动态监测。遥感技术具有大范围、快速、连续监测等优点, 所以被广泛用于区域尺度的草地植被长时间序列监测^[12-14]。目前, 很多学者开展了针对青藏高原草地退化的遥感监测研究。曹旭娟等人基于NDVI3g数据集计算了青藏高原长时间序列(1986—2013)草地退化状况, 发现青藏高原草地退化存在明显的空间差异^[15]。李重阳等人基于归一化植被指数(normalized difference vegetation index, NDVI)数据研究了青藏高原牧户草场退化趋势, 以及导致草场变化的主要驱动力, 发现草场退化面积达到总面积的69%^[16]。刘晓枫等人利用遥感影像获取色达县植被净初级生产力(net primary productivity, NPP)及土壤有机质(soil organic matter, SOM), 进而展开草地退化的研究^[17]。在众多植被、土壤指数中, NDVI能够有效反应草地面积、植被高度和盖度以及生产力的变化, 是较优的草地退化监测指标^[18]。基于遥感数据获取植被监测参数, 是进行草地退化大尺度研究的重要手段, 也是目前研究草地退化的前沿趋势。但目前针对青藏高原草地退化的研究集中在一些子区域, 且大多停留在草地退化的时空特征分析, 缺乏对高原草地退化的归因分析。

地表温度和土壤湿度都是陆地生态系统关键过程研究中重要的因子, 一方面, 在接收太阳辐射后, 地表温度升高, 导致土壤水分蒸发进入大气中, 形成地表与大气之间能量的交换；另一方面, 土壤水分通过垂直运移, 与地下水产生联系, 从而供给地表植被的生存用水^[19]。地表温度和土壤湿度在土壤的物理和生化过程中都起着重要的作用, 同时二者对植物种子的萌芽和植株的生长也都有着重要影响^[20-22]。水热要素对于植被生长来说是至关重要的角色^[23], 在区域尺度上的相关研究主要是针对空气温度和降雨量^[24-27], 对于地表水热要素的研究主要是集中在样方尺度和点尺度^[28-29], 而少见区域尺度地表水热要素对草地植被生长影响的综合研究^[30-31]。相比较于气温和降水, 地表温度和土壤湿度在土壤的物理和生化过程中都起着重要的作用, 二者直接作用于植物的根部, 对植物种子的萌芽和植株的生长也都有着重要影响。

基于此, 本文利用遥感数据, 获取了地表温度(land surface temperature, LST)和温度植被干旱指数(temperature vegetation dryness index, TVDI)两个地表水热参数, 并采用植被盖度指标中的NDVI为主要遥感监测指标, 建立草地退化遥感监测和评价指标体系。首先, 对青藏高原草地退化状况进行了时空变化监测；之后分析了青藏高原草地退化的时空特征和变化规律；最后, 探究了地表水热要素对草地退化的影响, 从而加深对青藏高原地区气候和植被变化机理的认识, 具有一定的生态价值和现实意义。

1 研究区概况

青藏高原(73°15′E—104°47′E, 26°00′N—39°47′N)西起帕米尔高原, 东至横断山脉, 北达昆仑山、祁连山脉, 南接喜马拉雅山脉, 其东西长约2800 km, 南北宽约300—1500 km, 总面积约为250万km², 平均海拔在4000 m以上。

青藏高原区域内高山大川密布, 地势险峻多变, 地形复杂, 其平均海拔远远超过同纬度周边地区。青藏高原各处高山参差不齐, 落差极大。总体来说, 青藏高原地势呈西高东低的特点。青藏高原年平均气温由东南的20℃, 向西北递减至-6℃以下。由于南部海洋暖湿气流受多重高山阻留, 年降水量也相应由2000 mm递减至50 mm以下。青藏高原地区天然草地面积约为1.5×10⁸ hm², 并且在高原特殊的气候环境影响下, 形成了以高寒草甸和高寒草原为主草地类型。

2 数据源和预处理 2.1 NDVI数据

通过遥感数据获得的NDVI时间序列数据能够反映陆地生态系统植被的生长状态、季相和年际变化特征, 因此被广泛用于全球尺度和区域尺度的生态环境变化监测与模拟、植被覆盖动态变化研究、植被物候特征识别与信息提取等。本文使用的NDVI数据为2001—2017年青藏高原地区的MODIS影像数据MOD13A2产品, 该产品为16 d合成产品, 空间分辨率为1 km。数据预处理上, 首先使用MRT (MODIS Reprojection Tool) 工具对上述数据进行批量投影变换和拼接处理, 再使用ArcGIS将数据进行批量值转换、水体掩膜、裁剪等处理。

NDVI由于表征的是植被生长状态, 通常反应植被生长过程的时间序列曲线是连续平滑的, 但是由于传感器本身性能、数据传输过程失误、太阳光照角度、观测视角、地物双向性反射以及云、大气气溶胶等观测条件因时间不同存在差异, 此外还受到地表水、冰雪等因素干扰, 得到的观测值包含很多不可预测的噪声, 所以NDVI时间序列曲线呈锯齿状不规则波动, 反映季节变化趋势不明显, 对研究结果带来很大的干扰。为此本文参考在青藏高原地区的相关研究, 选取了适合该地区的NDVI时间序列重建方法^[32], 即非对称性高斯函数拟合法(A—G拟合法)对2001—2017年的MODIS NDVI(16天/1 km)数据集进行处理, 最后通过平滑连接各高斯拟合曲线实现时间序列重建。

2.2 LST数据

地表温度是大气和地表连接处的温度状况, 在土壤的物理和生化过程中都起着重要的作用, 同时对植物种子的萌芽和植株的生长也都有着重要影响。作为与植物生长息息相关的环境要素, 地表温度的变化在很大程度上影响着地表植被的变化。

本文使用的LST数据为2001—2017年青藏高原地区的MOD11A2产品, 该数据为8 d合成产品, 空间分辨率为1 km。数据预处理上, 首先使用MRT(MODIS Reprojection Tool)工具对上述数据进行批量投影变换和拼接处理, 再使用ArcGIS将数据进行批量值转换、水体掩膜、裁剪等处理, 并将LST数据相邻两个时期的数据进行平均值计算处理, 获取与NDVI数据时间相匹配的16 d分辨率产品。

2.3 中国气象局陆面数据同化系统(CLDAS-V2.0)近实时产品数据集

CLDAS-V2.0产品为覆盖亚洲区域(0—65°N, 60—160°E)经纬度网格融合再分析产品, 其空间分辨率为0.0625°, 时间分辨率为1 h, 包括大气驱动场产品(2 m气温、2 m比湿、10 m风速、地面气压、降水、短波辐射6个要素)、土壤湿度产品(垂直分为5层：0—5、0—10、10—40、40—100、100—200 cm)产品。该数据集利用多种来源地面、卫星等观测资料, 采用多重网格变分同化、最优插值、概率密度函数匹配、物理反演、地形校正等技术研制而成, 在中国区域质量优于国际同类产品, 且时空分辨率更高。由于该数据集在中国地区仅有2017年至今的数据集, 时间范围上不满足本文需求, 但是空间分辨率较高, 所以本文采用该数据集中的2 m气温数据用作验证数据。

2.4 ERA5数据集

ERA5(The Fifth Generation of ECMWF Atmospheric Reanalyzes of the Global Climate)是欧洲中期天气预报中心(European Centre for Medium—Range Weather Forecasts, ECMWF)发布的第五代全球气候大气再分析网格数据集。再分析使用物理定律将模型数据与来自世界各地的观测结果组合成一个全局完整且一致的数据集, 以产生新的最佳估计大气状态, 该数据集包括从2000年至今的每小时0.25°×0.25°的气象因子数据, 本文选取其中的2 m高度的气温数据(2 m temperature)进行处理。

首先下载与Terre卫星过境时间(当地时间10：30)最接近的气温数据(NetCDF格式), 然后利用ArcGIS对数据进行格式转换、平均值合成、重采样等预处理获取于NDVI时空分辨率相匹配的2 m气温数据。其中, 重采样方法使用双线性内插法, 该插值方法适用于表示某种现象分布、地形表面的连续数据, 如DEM(Digital elevation model, 数字高程模型)、气温、降雨量分布、坡度等。本文将从中国气象局陆面数据同化系统(CLDAS-V2.0)中获取的气温作为基准值, 来研究重采样后ERA5数据的误差。

图 1为2017年日序数DOY(Day of year)=161重采样后的ERA5气温数据与中国气象局CLDAS气温数据散点图, R²=0.75, 相关性较好, 同时通过计算求得回归估计标准差RMSE为2 K, 说明重采样后的气温数据能够较好的反应研究区域内的气温状况, 可以用于后续计算。

2.5 植被覆盖度

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

(1)

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

3 研究方法 3.1 青藏高原草地退化遥感监测和评价指标体系的建立

根据国家标准GB19377—2003天然草地退化、沙化、盐渍化的分级指标, 20世纪80年代初期相同监测区域相同草地类型的草地植被特征可以作为退化草地的基准。因此, 以1982—1985年每个像元的年际最大植被覆盖度(AVHRR数据计算得到)作为基准, 将其退化程度分为未退化、轻度退化、中度退化、重度退化和极重度退化5个级别^[33], 具体见表 1。

3.2 区域退化指数

在草地退化等级划分的基础上, 采用Gao等人提出的草地退化指数(grassland degradation index, GDI)。草地退化指数的计算公式为：

(2)

式中, GDI为区域草地退化指数, D_i为草地退化等级i的评分, A_i为草地退化等级i的分布面积, A为研究区与的总面积。依据退化分级标准分析草地退化情况, 划分标准如表 2所示。

3.3 改进温度植被干旱指数

土壤湿度作为一个重要的环境参数, 是求解区域水量平衡、研究全球水循环过程的关键因子。一方面, 土壤水分通过接收太阳辐射, 地表温度升高, 使水分蒸发进入大气中, 形成地表与大气之间能量的交换；另一方面, 土壤水分通过垂直运移, 与地下水产生联系, 从而供给地表植被的生存用水, 将地表水与地下水联系在一起。此外, 土壤湿度的时空分布与动态变化不仅会影响其自身, 对地表反照率、蒸散发、生物化学循环、植被生长等也有着重要的作用。

1994年, Carlson等人在研究不同分辨率的LST和NDVI数据时, 发现LST与NDVI之间呈现出明显的负相关关系^[34]。之后随着不断研究发现, 当研究区域内的植被覆盖和土壤湿度变化较大时, 通过遥感影像获得的NDVI和LST所构成的散点图呈梯形或三角形。

实验中, 当研究区域足够大, 而且土地覆盖类型从裸土变化到完全植被覆盖, 土壤湿度从干旱变化到湿润, NDVI和LST构成的空间关系为三角形。如图 2所示, A点表示干燥裸土, 即NDVI小, LST高；B点表示湿润裸土, 即NDVI小, LST也小；C点表示植被完全覆盖, 同时土壤水分也很充足, 即NDVI大, LST低, 此时蒸散阻抗小(蒸散包括蒸发和植物的蒸腾, 蒸发又包括裸土的蒸发和植被冠层截获水分的直接蒸发)。AC边表示, 对于最高温度下, 不同植被覆盖类型土壤湿度最低的状况, 我们称这条边为“干边”。BC边表示, 在这一温度下, 土壤水分充足, 不会成为植被生长的限制因素, 我们称之为“湿边”。

当水分条件充足时, 植被生长迅速, NDVI值变高, 同时区域内植被生命活动旺盛, 蒸腾量变大, 整个像元内的蒸散阻抗变小, 潜热所占比例增大, 所以像元内的地表温度会降低；当水分不足时, 植物受到水分胁迫, 植物叶片气孔关闭, 降低了蒸腾所造成的水分损失, 进而导致地表潜热通量的降低, 根据能量平衡原理, 地表能量必须平衡, 所以地表感热通量会增加, 感热通量的增加最终会导致冠层温度的升高。

Sandholt等在研究利用地表温度—植被指数三角形特征空间研究土壤湿度时, 提出了利用温度植被干旱指数(TVDI)来表示土壤干湿状况^[35], 定义为：

(3)

式中, T表示任意像元的地表温度值；T_max表示某一NDVI对应的最大温度, 可以由NDVI与干边线性拟合得到；T_min表示某一NDVI对应的最小温度, 可以由NDVI与湿边线性拟合得到。TVDI的取值范围为0—1, TVDI值越大, 说明这一像元的土壤湿度越低, TVDI值越小, 说明这一像元的土壤湿度越高。

在三角形特征空间中对干湿边进行线性拟合, 可以得到干边方程为：

(4)

湿边方程为：

(5)

将式(4)和式(5)代入式(3)中, 可以得到TVDI的计算公式为：

(6)

式中, a₁, b₁和a₂, b₂分别是干边和湿边的线性拟合的系数。

如图 2所示, 在AC干边上, 不同植被指数下的温度最高, 土壤含水量低, NDVI低的点干燥迅速, 蒸发量少, 吸收较多的太阳能, 为干旱状态, TVDI为1；在BC湿边上, 土壤含水量高, 吸收的太阳能主要用来蒸发和蒸腾作用, 裸土表面和作物冠层温度差距不明显, BC边接近水平线, 为湿润状态, TVDI为0。

由于某个区域某一段时间内不同像元的NDVI与LST值都分布在干湿边界构成的特征空间内, 所以我们可以将这个特征空间看作是由土壤湿度等值线组成。地表温度—植被干旱指数TVDI可以用来表征区域尺度地表土壤含水量情况, 也可以称为温度植被干旱指数, 下文中将统一用温度植被干旱指数表示^[36]。

在TVDI模型中, 假设TVDI变化的主要来源是土壤水分, 没有考虑空气温度。当TVDI进行反演和应用时, 通常假设空气温度在研究区域内是固定的, 这在大区域进行应用时明显是错误的^[37]。

在本文中, 青藏高原地区面积较大, 且地面异质性较大, 所以考虑将空气温度引入到TVDI方程中, 参照Parinaz等的相关研究^[38], 可以得到改进后的温度植被干旱指数计算公式为：

(7)

构建TVDI特征空间时, 本文采用的算法将NDVI每间隔0.025划分步长, 提取同一NDVI数值下的最大最小温度(置信区间95%), 根据最大最小温度提取特征空间的干边和湿边, 考虑到NDVI在0.1以下和0.9以上对植被生长状况指示存在不灵敏或饱和的情况, 本文在拟合干湿边界时, 只提取NDVI在0.1—0.9之间的值。

如图 3中所示, 由于选取验证的数据的时间处于6月份, 此时青藏高原地表温度较高, 地表土壤含水量较低, 对应的TVDI较高, 所以散点集中分布在右下部分。另外, 改进后的TVDI_i比TVDI在1：1反比线左右更加聚合。综上可得改进后的温度植被干旱指数TVDI_i在研究区域更加适用, 更能够反应土壤水分的变化。

以2009年DOY=161为例, 构建改进后的温度植被干旱指数的特征空间时, 考虑到NDVI在0.1以下和0.9以上对植被生长状况指示存在不灵敏或饱和的情况, 本文在拟合干湿边界时, 只提取NDVI在0.1—0.9之间的值。通过图 3和表 3得知改进后的TVDI对干湿边界的拟合效果更好。

3.4 变化趋势分析方法

Slope定义为在一定时间范围内采用最小二乘法拟合年度变量均值的斜率, 它能够反映出每个栅格的变化趋势。计算公式如下：

(8)

式中, i为年份序号, 取值范围为1—17, n为研究的时间序列长度, 为17, X_i为第i年的年度平均值(例如当X为NDVI时, 这里表示的则为NDVI的年度平均值)。若Slope大于0则说明该像元的变量X在这17年间的变化趋势是增加的, 反之则是减少的。

3.5 多元线性回归分析方法与贡献率

(9)

式中, Y_i为因变量；X_ki为自变量；i=1, 2, …, n；k为解释变量的个数, β_k为回归系数。

本章节将以地表温度LST和温度植被干旱指数TVDI(改进后的TVDI, 下同)两个因子作为回归模型的自变量, 然后将NDVI作为因变量, 来建立回归模型。参考Zhu等人的研究^[39], 如果在建立回归模型之前, 将遥感栅格数据均进行归一化处理, 然后用17年的年尺度归一化数据建立回归模型, 最后求得的每个自变量的系数则为该像元该自变量对因变量的定量贡献率, 并且可以通过对比不同自变量定量贡献率的大小来决定主要影响因子。

首先对年尺度的数据进行归一化处理。之后基于MATLAB建立多元线性回归模型, 逐像元计算LST和TVDI二元变量的系数及常数项, 将系数和常数项分别出图, 得到LST和TVDI对因变量的贡献率空间分布图像及反应其他因子贡献率的残差图像。

4 研究结果与分析 4.1 青藏高原草地退化现状分析

从图 4, 2001和2017年两年的青藏高原地区草地退化等级空间分布可以看出, 青藏高原植被退化程度空间差异明显, 柴达木盆地和青海湖附近退化较为严重, 喜马拉雅山脉北部、昆仑山脉南部、冈底斯山脉北部交汇的地区退化也较严重。但是2017年总体退化程度较2001年轻。

从2001和2017年两年的青藏高原地区植被退化等级占比随时间变化图可以看出来(图 5, 图 6), 2001—2017年每年植被未退化面积基本处于50%—60%之间, 且不同退化等级的面积随着时间的变化也均有变化, 其中未退化面积从50.6%上升到59%, 说明青藏高原植被退化整体上在朝着改善的方向发展。

GDI指数主要用来计算某个区域整体退化状况。图 7为青藏高原草地区域2001—2017年的每年草地区域的GDI指数。整体上, 青藏高原草地区域GDI指数大部分处于1.5—2之间, 说明大部分时间处于轻度退化状况, 而2001年和2015年GDI均达到了2, 趋向中度退化, 说明这两年草地植被生长出现了恶化。但是整体上, 从2001—2017年, 草地区域GDI指数呈减小的趋势, 说明青藏高原草地生长这些年之间退化情况在减轻。

4.2 青藏高原地表水热要素对草地退化的影响

参考Zhu等人的研究^[39], 如果在建立回归模型之前, 将遥感栅格数据均进行归一化处理, 然后用17年的年尺度归一化数据建立回归模型, 最后求得的每个自变量的系数则为该自变量对因变量的定量贡献率, 并且可以通过对比不同自变量定量贡献率的大小来决定主要影响因子。通过逐像元计算方式, 可以获得各因变量的贡献率在空间上的分布情况。首先计算TVDI与LST对NDVI的二元回归关系, 其回归系数可以反映因子对植被影响的大小(图 8)。当回归系数在-0.5—0.5之间, 可以认为因子对植被几乎没有影响；当回归系数为-1.0至-0.5, 0.5至1, 因子对植被有一般程度的影响；当回归系数在小于-1或大于1因子对植被具有较显著的影响。地表温度对NDVI的回归系数主要为负数, 说明大部分区域地表温度与植被NDVI的关系为反向关系, 而对TVDI, 与植被生长状况主要呈现正向关系。

通过比较TVDI、LST回归系数与常数项的大小, 确定青藏高原区域各像元的主导因子, 结果如图 9所示。其中, 在青藏高原东南部, 水热因子对植被的影响相对较小。水分主导的区域主要分布在青藏高原北部, 相对较干旱的区域。其次, 地表温度在大部分地区都存在与植被的显著相关关系。从草地植被类型角度, 统计了水热因子在各类型草地中主导区域的面积比例, 其结果如图 10所示。地表温度在各个植被中均占据较大的主导面积。地表温度在温性草原、高寒草原、温性草甸、高寒草甸中占据主导面积的比例分别为32.73%, 42.24%, 34.59%, 26.70%。TVDI在温性草原、高寒草原、温性草甸、高寒草甸中占据主导面积的比例分别为7.83%, 16.92%, 8.65%, 10.80%。

5 讨论

从计算结果来看, 地表温度对植被NDVI具有较显著的回归关系, 但是很难说明地表温度对植被具有很高的影响。地表温度与地表湍流能量分配有关, 而地表湍流能量分配主要取决于植被覆盖度和土壤湿度^[40]。因此, 在特定的土壤水分含水量下, 植被NDVI和LST之间存在很强的线性关系, 这一关系也被用于发展计算土壤水分亏缺状况的方法, 如三角形特征空间模型^[41]。为了更深一步认识LST与NDVI之间的强相关关系, 本文选取了2014年的NDVI与LST数据, 通过对NDVI按照0.005间隔划分, 计算相应区间内LST均值, 得到图 11统计结果。从图 11可以看出, 总体上, NDVI越大, 对应的LST越小。由于遥感观测的LST包括土壤组分和植被组分两部分, 相对于裸土, 在相同的辐射下, 植被组分的温度更低。因此, 尽管回归模型中, LST对NDVI的回归关系较优, LST对植被生长的关系也更为复杂。在低植被覆盖区域, LST主要组分来自裸土的情况下, 地表温度的增高有利于植被的生长, 同时, 在高植被覆盖区域, 植被覆盖度会反过来影响地表温度, 随着覆盖度的增大, 地表温度降低。从回归系数来看(图 8), LST对NDVI的回归系数也偏向于负数, 说明了LST与NDVI之间较强的负相关关系。进一步地, 建议获取土壤组分温度, 从而深入研究区域尺度上土壤温度对草地植被生长的影响。

青藏高原在2001—2017年间, 极重度退化和重度退化面积比例变化很小, 说明重度及以上退化一般呈现不可逆的退化趋势, 短时间内很难恢复。同时, 未退化草地面积呈现逐年上升趋势, 说明青藏高原草地在2001—2017年间的退化状况正在改善, 且主要是从轻度、中度退化过渡到未退化。如图 4所示, 青藏高原东北部地区(主要分布在海西蒙古族藏族自治州, 柴达木沙漠周围地区)呈现出较为明显的草地退化改善趋势。从图 8, 这些地区的TVDI与NDVI具有较高的相关系数, 土壤湿度的增大可能是导致地区草地退化改善的主要原因, 其次, 地表温度与植被NDVI较高的负相关关系, 也体现出草地植被NDVI的增大趋势。草地退化同时受多种因素影响, 除了自然因素, 人类活动在玉树藏族自治州、甘南藏族自治州等以牧场草地为主的地区也有重要的影响^{[15-16, 42]}。由于近年来有规划的合理放牧, 2000年至2016年, 牧场草地的退化面积也呈现明显的上升趋势^[16]。从植被类型角度, 受地表水热影响较大的草地类型为高寒草原, 其次为温性草甸、温性草原和高寒草甸。高寒草原在青藏高原的分布面积最大, 相对于其他草地类型, 受人类活动影响的面积比例也较低, 因此, 地表水热等自然因素对该类型草地影响较大。总体来看, 地表水对草地植被的恢复具有积极影响, 地表含水量同时也反映了土壤的持水能力和土壤质地状况^[43], 而草地退化的根本原因之一可能为土壤的退化。

同时, 在研究过程中仍存在需要完善和值得深入思考的地方, 比如本文主要研究内容为青藏高原植被退化对地表水热要素的响应(LST、TVDI), 其他的影响因子如太阳辐射、人类活动等对植被生长均有影响, 本文也将除地表水热因子以外的影响因子统称为“其他影响因子”进行了分析, 在今后的研究中可以考虑将更多的影响因子单独纳入分析中。同时, 由于青藏高原地区范围广、气象监测站点少、自然环境恶劣, 其他因素比如气象因子或人为因子数据源难以和遥感数据进行空间尺度匹配, 因此在今后的研究中应该加强各种区域尺度数据的获取, 开发适合青藏高原复杂下垫面的气象要素插值方法, 研究将县、市、省级的人为影响因子进行尺度转化的方法, 从而与常见的遥感数据尺度进行匹配。

6 结论

(1) 研究结果表明从2001到2017年间, 青藏高原植被退化程度空间差异明显, 柴达木盆地和青海湖附近退化较为严重, 喜马拉雅山脉北部、昆仑山脉南部、冈底斯山脉北部交汇的地区退化也较严重。

(2) 从2001—2017年间, 青藏高原草地未退化面积从50.6%上升到59%, 说明青藏高原植被退化整体上在朝着改善的方向发展。2001—2017年内, 大部分年份GDI指数在1.5—2之间, 说明青藏高原草地大部分时间处于轻度退化状态, 其中, 2001年和2015年, GDI达到了2以上, 这两个年份青藏高原草地退化达到中等退化水平。整体上看, GDI呈现减小趋势, 青藏高原草地退化在减轻。

(3) 青藏高原草地植被14.04%的区域主要受到土壤水分含量影响, 从回归关系来看, 36.61%的区域NDVI与地表温度有较强的相关关系, 但主要反映在植被覆盖影响地表温度。综合来看, 青藏高原地表水热对草地植被退化具有显著影响。

致谢: 国家气象科学数据中心(http://data.cma.cn/site/index.html)提供气象数据产品；美国国家航空航天局土地流程分布式活动档案中心提供MODIS影像数据(https://search.earthdata.nasa.gov/search)；欧洲中期天气预报中心提供再分析气象数据产品(https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels)；中国科学院空天信息创新研究院遥感科学国家重点实验室生态系统遥感研究室提供青藏高原覆被数据, 特此致谢。