植被是生态系统的重要组成部分，也是连结土壤、大气和水分的重要纽带，在调节气候、保持水土以及物质循环中均扮演着重要的角色^[1-3]。遥感以其空间覆盖广、获取速度快以及周期性观测等优势，已广泛应用于大面积观测植被覆盖及其变化情况的研究中。归一化植被指数(normalized difference vegetation index, NDVI) 与叶面积指数和生物量等均有较高的相关性，可以反映地表植被的茂盛程度，在一定程度上能够表示地表植被覆盖的变化^[4-7]。

青藏高原被称为地球的“第三级”，其生态环境的变化对全球生态^[8]系统具有深刻的影响。因此，对青藏高原地区植被的分布规律及其动态演变进行研究，具有重大的生态价值与现实意义。赵紫薇^[8]使用1982—2012年GIMMS-NDVI和2001—2013年MODIS-NDVI数据发现，30 a来青藏高原植被整体呈改善趋势；卓嘎等^[9]利用2000—2016年MODIS NDVI数据对青藏高原植被时空分布进行探究，发现东南部地区植被状况明显好于西北部；孔冬冬等^[10]将青藏高原划分成8个植被分区，使用GIMMS NDVI 3g数据，发现物候变化最剧烈的植被分区是高寒草甸和高寒灌木草甸；梁大林等^[11]以阿里高寒草原和三江源高寒草甸为研究对象，发现2000—2017年两地超过50%的植被NPP均呈上升趋势，且高寒草甸植被NPP增长趋势更明显。

近年来，对青藏高原的植被覆盖研究多基于不同的气候分区，较少有研究从整个研究区的尺度上对植被覆盖整体进行等级划分，进而分析其时间的变化和空间的转移，且研究对近几年的时间尺度涉及较少。本文使用MOD13Q1-NDVI数据，利用NDVI值将青藏高原地区植被覆盖程度划分为不同的等级，使用一元线性趋势回归和马尔科夫矩阵模型分析整个研究区植被时间和空间的变化情况。对植被的影响因子研究中，多数研究选取气象站点气温数据分析气温对植被生长的影响。地表温度同样会对植被的生长产生一定的影响，而基于地表温度对植被影响的研究较少，因此，本文选取MOD11A2-LST地表温度，并结合TRMM 3B43降水数据，使用偏相关分析对研究区2000—2019年20 a来植被与地表温度和降水的相关性进行探究。

1 研究区概况

青藏高原是世界上海拔最高的高原，区域范围在北纬26°00′~39°47′，东经73°19′~104°47′之间，是反映全球气候变化的敏感地带。在中国境内，青藏高原西起帕米尔高原，东接横断山脉，北靠昆仑-祁连山北部，南抵喜马拉雅山脉南缘，跨越青海、西藏等6个省区^[12]。青藏高原地区地形复杂，包括高山、盆地、湖泊、冰川等多种地貌类型。气候复杂多样，全年地表温度偏低，不同地区间气候差异明显，干湿分明，西北地区寒冷干旱，东南部分地区较暖热湿润。青藏高原的植被类型丰富多样，主要类型包括高寒灌丛、高寒草原、高寒荒漠等，东南地区也分布有阔叶林和针叶林^[13]。

2 数据来源与研究方法 2.1 数据来源与处理 2.1.1 植被指数产品MOD13Q1数据

MODIS植被产品数据MOD13Q1来自美国航空航天局数据中心，该数据的空间分辨率为250 m，时间分辨率为16 d。该产品经过了气溶胶、水、云等处理，数据质量高，因此已被广泛运用于植被覆盖变化的研究中^[14]。本研究提取2000—2019年连续20 a的NDVI数据作为植被数据源。

2.1.2 地表温度产品MOD11A2数据

MODIS-MOD11A2数据来自美国航空航天数据局中心，空间分辨率为1 km，由每日地表温度产品MOD11A1合成，存储了8 d中晴朗天气下的地表温度平均值，消除了部分云的干扰，提高了LST数据的有效性。该数据采用分裂窗算法，使用MODIS的31、32波段反演得到。由于其像元值是依据LST的真值转化得到的易于存储的数值，因此，需要将其转化为真实的地表温度值。计算公式如下：

$ \theta = 0.{\rm{ }}02p - 273.{\rm{ }}15 $

(1)

式中：θ为像元的真实地表温度值(℃)，p为像元的灰度值(K)。

2.1.3 TRMM降水数据

所使用的TRMM 3B43第7版数据是TRMM卫星和其他卫星以及地面观测数据联合反演出的月降水产品，其空间分辨率为0. 25°×0. 25°。该数据覆盖全球50°S~50°N的地区。李函璐^[15]在对TRMM数据在降水研究现状中提到，TRMM卫星资料已广泛应用于喜马拉雅山、青藏高原、太平洋等地降水的研究，且在热带和亚热带地区研究成熟，在青藏高原等中纬度地区具有较好的适用性。

2.1.4 数据预处理

由于各数据的投影、空间分辨率不同，为便于计算，将NDVI、地表温度以及降水数据转换为统一的Albers投影，并将空间分辨率重采样为1 km。使用月最大值合成法将16 d NDVI数据合成月尺度数据，再将月尺度数据通过年最大值合成法合成年尺度数据。对地表温度数据进行年平均值合成法处理，将TRMM逐月降水数据合成为年降水量，最后对各数据进行研究区的裁剪得到2000—2019逐年数据集。本研究所采用的青藏高原边界数据总集^[16-19]来源于“国家青藏高原科学数据中心” (http://data.tpdc.ac.cn)。

2.2 研究方法 2.2.1 一元线性回归趋势分析

采用一元线性回归模型，对青藏高原地区2000—2019年20 a的植被覆盖程度变化趋势进行定量分析，计算公式为：

$ S = \frac{{t\sum\limits_{i = 1}^t {\left( {i{N_i}} \right)} - \sum\limits_{i = 1}^t i \sum\limits_{i = 1}^t {{N_i}} }}{{t\sum\limits_{i = 1}^t {{i^2}} - {{(\sum\limits_{i = 1}^t i )}^2}}} $

(2)

式中: S表示趋势线的斜率，N_i表示像元第i年的NDVI值，t为数据的监测时间(t=20 a)；S＞0，表示NDVI值增加，即植被覆盖呈改善趋势，反之表示植被覆盖呈退化趋势。在趋势分析的基础上，使用F检验，判断植被NDVI变化趋势的显著性。

2.2.2 马尔科夫模型

空间转移矩阵可以反映土地利用的结构特征以及各类型之间的转换情况与方向, 表征各个类型的转入转出情况。本研究利用马尔可夫模型，定量表征青藏高原地区不同等级的植被覆盖程度变化情况，计算公式如下：

$ {\mathit{\boldsymbol{V}}_{ij}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{V}}_{11}}{\mathit{\boldsymbol{V}}_{12}} \ldots {\mathit{\boldsymbol{V}}_{1h}}}\\ {{\mathit{\boldsymbol{V}}_{21}}{\mathit{\boldsymbol{V}}_{22}} \ldots {\mathit{\boldsymbol{V}}_{2h}}}\\ { \vdots \;\;\; \vdots \;\;\; \ddots \;\;\; \vdots }\\ {{\mathit{\boldsymbol{V}}_{h1}}{\mathit{\boldsymbol{V}}_{h2}} \ldots {\mathit{\boldsymbol{V}}_{hh}}} \end{array}} \right] $

(3)

式中：V_ij表示第i等级植被覆盖程度转换成第j等级植被覆盖程度的面积(本研究采用面积百分比代替)，h为植被覆盖程度等级的总类别，i和j分别为转移前后的植被覆盖程度等级。

2.2.3 偏相关分析

偏相关系数的大小可以反映地表温度和降水对植被的影响，该方法被广泛应用于研究地表温度和降水对植被生长的影响。其计算公式如下：

$ {r_{xy, z}} = \frac{{{r_{xy}} - {r_{xz}}{r_{yz}}}}{{\sqrt {1 - {r^2}_{xz}} \sqrt {1 - {r^2}_{yz}} }} $

(4)

式中：变量x表示NDVI的像元值，y和z分别表示地表温度和降水的像元值，r_{xy, z}表示控制降水的情况下NDVI与地表温度的偏相关系数。r_{xy, z}的取值范围为-1~1，当r_{xy, z}＞0时，表示NDVI与地表温度呈正相关，反之，则二者呈负相关；且偏相关系数的绝对值越接近于1，表示二者的相关性越强。

在偏相关分析的基础上，采用T检验法判断偏相关系数的显著性，其公式如下：

$ w = \frac{{{r_{xy, z}}}}{{1 - {r^2}_{xy, z}}}\sqrt {n - m - 1} $

(5)

式中：n表示样本数，m表示自变量的个数。根据T检验临界值分布表，可以得到在α显著水平上的临界值w_α，若w>w_α，则表示显著偏相关，反之，则表示相关性不显著。

3 结果与分析 3.1 归一化植被指数时间变化特征

青藏高原地区2000—2019年NDVI最大值的面平均值变化趋势如图 1所示。可以看出，20 a来研究区NDVI值的变化呈波动上升趋势，2012—2015年间，植被有退化趋势，自2016年逐渐改善，年均NDVI值为0. 35~0. 40。总体来看，研究区NDVI值增长速度较缓，增速为0. 001 2 a^-1 (R²=0.607 4)，植被覆盖呈轻微改善趋势。

3.2 归一化植被指数空间变化特征

根据青藏高原地区NDVI值的大小，将研究区植被覆盖状况划分为5个等级，如表 1。

从空间上看(图 2)，青藏高原地区NDVI呈现由北到南、由西向东的递增趋势。低植被覆盖区主要分布在研究区的西部和北部，包括新疆维吾尔自治区的南部、西藏自治区北部以及青海省的西北部，高植被覆盖区主要分布在研究区的东部和南部地区，主要位于四川省、云南省境内。根据统计结果显示，植被覆盖面积所占比例最大的是高植被覆盖区域，占研究区总面积的28. 22%，其次为低覆盖区域，占比24. 32%，较高植被覆盖等级所占比例较少，仅为10. 26%。

根据马尔科夫模型构建青藏高原地区植被覆盖程度等级转移矩阵，计算结果如表 2, 可以看出，2000—2019年，研究区裸地区域明显减少，主要向低覆盖区转换，从2000年的21. 05%减至2019年的15. 79%，高植被覆盖区有所增加，主要由中覆盖与较高覆盖等级转入，到2019年，增至28. 88%。低覆盖、中覆盖以及较高覆盖植被面积均发生轻微变化。总体来看，植被转换的趋势呈正向演进，植被覆盖等级保持不变的区域占76. 28%，有所改善的区域占总面积的18. 37%，退化区域占比仅5. 35%。

3.3 归一化植被指数变化趋势分析

根据一元线性回归与F检验的结果，将青藏高原地区NDVI的变化情况划分为5个等级，如表 3。

根据经过显著性检验的NDVI趋势结果(图 3)，青藏高原地区NDVI显著增大的面积占研究区的27. 08%，较显著增大面积占比6. 95%；NDVI显著减小和较显著减小的面积分别占1. 81%和1. 25%。由此可见，近20 a来青藏高原地区植被状况整体趋于好转，植被退化较少。从空间上看，NDVI显著增大的区域主要分布在研究区的北部，例如新疆维吾尔自治区的南部、青海省西北部以及西藏自治区的北部等，南部分布较为零散，呈破碎斑块状，部分海拔较高地区NDVI显著增大，如西藏的阿里、日喀则等地区。研究区南部与东部地区NDVI总体趋于稳定，变化并不显著。

3.4 归一化植被指数变化与温湿度因子的相关性 3.4.1 地表温度和降水的时间和空间分布特征

地表温度、降水以及光照等气候因子在植被生长和发育的过程中是不可或缺的，其变化会造成植被生长环境的改变，从而影响植被的生长状况。本研究选取地表温度和降水2个因子，分析青藏高原地区NDVI变化与二者的相关性。

青藏高原地区总体地表温度较低，从时间上看(图 4)，研究区2000—2019年平均地表温度介于11. 5~14. 0 ℃，最低为11. 7℃，在2016年地表温度达到最高值，为14. 0 ℃。20 a来青藏高原地区年降水量介于400~510 mm，整体呈波动上升的趋势，降水量最低值在2006年，为402 mm，2018年，降水量达到最高值，为510 mm。

从空间上看(图 5)，青藏高原地区地表温度分布受海拔以及纬度影响较大，海拔和纬度越高，地表温度越低。地表温度低值区主要分布在研究区西北部的昆仑山地区以及东北部的祁连山地区，高值中心出现在东北部的柴达木盆地地区，年均地表温度达到20. 0 ℃以上，东南部雅鲁藏布江谷底地区，受印度洋暖湿气流的影响，其地表温度也较高。研究区降水自西北向东南逐渐增多，与NDVI的空间分布具有较高的一致性。西北地区降水量较少，普遍低于300 mm，NDVI值也较低，以小于0. 20分布为主，东南地区NDVI值较高，大多高于0. 60，降水量也较多，部分地区超过1 500 mm。

3.4.2 植被与温湿度因子的偏相关分析

从像元角度对地表温度和降水与青藏高原地区植被之间的关系进行偏相关分析，计算每个像元的NDVI与地表温度和降水因子的偏相关系数和显著性(如图 6和图 7所示)。

如图 6，青藏高原植被NDVI与地表温度的偏相关分析结果显示，青藏高原地区NDVI与地表温度的偏相关系数介于-0. 96~0. 95之间，NDVI与地表温度呈正相关的区域占研究区的36. 68%，主要分布在西藏的西北部、青海的南部等地；NDVI与地表温度呈负相关的区域占总面积的63. 14%，主要分布于青海的西南部以及西藏的中部地区。其中，NDVI与地表温度呈显著正相关(p < 0.01) 的区域为2. 08%，呈显著负相关的区域为8. 42%。

如图 7，青藏高原植被NDVI与降水的偏相关分析结果显示，青藏高原地区NDVI与降水的偏相关系数介于-0. 92~0. 95，平均值为0. 09。降水与植被呈正相关的地区占总面积的64. 63%，主要分布在西藏的中部地区和青海的中南部，其中，NDVI与降水呈显著正相关(p < 0.01) 的地区占9. 34%，主要分布在研究区的中北部地区，这些地区年降水量相对不高，低于600 m，降水量的增加会对植被的生长产生促进作用；降水与NDVI呈负相关的地区占总面积的35. 37%，主要分布在西藏东南部；其中呈显著负相关(p < 0.01) 的地区占总面积的1. 69%，分布较为零散。

4 讨论与结论

本研究基于2000—2019年连续的MOD13Q1数据，计算了青藏高原20 a来NDVI变化情况，利用一元线性回归模型分析了区域NDVI变化趋势，同时结合地表温度和降水两个温湿度因子，探究了青藏高原植被变化与二者之间的关系，得到的结论如下：

从时间上，2000—2019年20 a间青藏高原地区植被整体处于改善的趋势，具体表现为研究区北部NDVI值显著增加，南部NDVI值变化不显著，退化区域面积较小；从空间上，青藏高原地区NDVI变化表现为从西北低，东南高，北部干旱地区NDVI值较小，大部分小于0. 20，而东南地区林地较多，整体NDVI值较高，大于0. 60，此结论与韩炳宏等^[20]、卓嘎等^[9]研究结果一致。

青藏高原地区植被与地表温度和降水存在一定的相关性，在西藏的中部以及青海的中南部地区，NDVI值与降水呈现正相关性，偏相关系数较高；在西藏的西北部、青海的南部等地，NDVI值与地表温度呈正相关关系，在西藏中部等地，NDVI与地表温度呈负相关关系。此结论与薛宇轩等^[21]研究结果基本一致，同时，其在文中也探究了NDVI与气象站点获取温度数据的相关性，研究结果表明，在青藏高原地区，地表温度与气温的变化趋势相似，而地表温度的升幅一般高于气温。因此，相比于地表温度，气温可能更易成为植被生长的限制因子，即地表温度更容易与NDVI呈现负相关关系。

本研究结果可对青藏高原地区的生态环境保护以及规划提供一定的参考，但还存在一些不足。由于气候变化存在着随机波动，需要在大的时间尺度上才能显示出来，而由于数据的时间尺度有限，本研究的时间范围为2000—2019年20 a间，在今后将采取更大范围的时间尺度进行更深一步的探究。同时，地区的植被除了与地表温度、降水等气候因素相关，和人类活动等各种因素均存在一定的相关性，但由于缺乏相关数据，无法进一步从更深的层次作出解释。今后将收集更多相关材料，对青藏高原地区的植被及其生长相关因子做出进一步研究。

致谢: 文中青藏高原边界数据来源于“国家青藏高原科学数据中心” (http://data.tpdc.ac.cn)，感谢数据中心为本文提供的数据支持！