雅鲁藏布江流域植被覆盖变化及其对气候变化的响应

韩先明，左德鹏，李佩君，徐宗学，高晓曦

(1. 北京师范大学水科学研究院城市水循环与海绵城市技术北京市重点实验室，北京 100875;2.水利部宣传教育中心，北京 100053)

植被作为重要的自然资源，是人类生存环境中不可或缺的组成部分。植被覆盖与气象要素密切相关，通过水汽输送的方式，植被适应着气候变化并在一定程度上表征气候变化。同时，植被参与并影响全球水循环、生化循环、能量平衡等过程，对生态系统变化产生深远的影响。植被覆盖是连接气象和水文的载体，通过生态系统的反馈可对其他生态要素产生不同程度的作用。植被指数在评价植被覆盖度、生长活力和生物量等植被状态信息时发挥着重要作用，其中归一化植被指数[1-3](normalized difference vegetation index,NDVI)是表征植被生长状态和植被空间分布密度的指示因子，能够很好地反映植被的动态变化，在气象、水文等领域的研究中被广泛应用。

雅鲁藏布江作为我国西南地区最大的外流河，地处青藏高原腹地，区位海拔跨度极大。流域内植被及生态环境十分脆弱，是我国重要的地理屏障，因此研究雅鲁藏布江流域植被变化及水文过程，提高对该地区植被覆盖变化的认识，厘清植被覆盖变化同气象要素的关系对我国的水能开发、环境保护、地缘政治有着重大战略意义[4-9]。许多学者[10-14]在不同地区对NDVI与气象要素的关系进行了探究，姜琳等[15]研究了降雨因子对雅鲁藏布江流域植被生长的影响；姜红涛等[16]利用NDVI、降水和气温数据，分析了艾比湖流域NDVI变化趋势及其对降水和气温的响应特征，指出NDVI总体上与年降水量呈正相关，而与平均气温呈负相关。朴世龙等[17]发现近年来我国大部分地区NDVI呈增长趋势，西北干旱区域和青藏高寒区域NDVI增长率在夏季达到最大，东部季风区域NDVI增长率在春季达到最大。大多研究表明NDVI对降水的响应存在一定滞后性[18-19]，同降水及气温存在正相关关系[20]；NDVI与海拔也存在一定关系，一般来说低海拔地区NDVI值相对较高[21]。就雅鲁藏布江流域而言，Li等[22]通过遥感手段对1999—2013年雅鲁藏布江流域NDVI年际变化及其沿海拔梯度变化特征进行了分析，发现以阔叶林和针叶林为主的低海拔地区植被覆盖增加明显。陈斌[23]利用SPOT数据集对雅鲁藏布江大峡谷自然保护区不同植被类型、海拔梯度下的NDVI时空变化规律进行了研究。吕洋等[24]分析了雅鲁藏布江流域植被覆盖时空变化规律，发现NDVI与降水之间存在较为明显的正相关关系。

已有研究对雅鲁藏布江流域植被覆盖变化进行了相关分析，然而不同植被类型下的植被生长状况及其持续性分析还未做充分探讨。考虑到研究区测站稀疏，属于缺资料地区，本文基于MODIS高精度遥感影像，采用趋势分析、Hurst指数、偏相关性分析对研究区2000—2016年生长季、年内和年际植被覆盖空间异质性做详细分析，探讨植被覆盖对气候变化的响应，以期为研究区生态环境保护和气候变化研究提供参考。

图1 研究区概况及气象站点分布

1 研究区概况和数据

1.1 研究区概况

雅鲁藏布江流域[25]位于北纬27°8′～31°2′，东经82°0′～97°1′之间，地处我国西藏自治区境内，是世界上平均海拔最高的流域。流域发源于西藏西南部喜马拉雅山脉北麓的杰马央宗冰川，上游称为马泉河[26]。流域海拔高度绝大部分在3 000 m以上，从源头杰马央宗冰川到下游河谷,由干旱、半干旱半湿润、湿润地区组成。流域东西狭长、南北窄短，南有喜马拉雅山脉，北有冈底斯山及念青唐古拉山脉，南北之间为藏南谷地，在米林以下的下游地区分布着著名的大拐弯地带[27]。流域降水年际变化不大，多年平均降水量300～500 mm，但空间分布差异极大且年内分配不均，部分地区夏季降水可达全年降水量80%。

通过咨询西藏水利专家并查阅相关文献[28]，将雅鲁藏布江流域分为上、中、下游，上游为日喀则市拉孜县以上宽谷区，植被覆盖多以草甸和草原为主；中游为河谷区，植被覆盖多以高山植被和草甸为主；米林县派镇以下为下游峡谷区，植被覆盖多以灌丛、针叶林和阔叶林为主。研究区水文气象站点主要位于流域中、下游且分布不均，属于缺资料流域，气象站分布、高程及流域概况如图1所示。

1.2 数据来源及处理

1.2.1植被数据

NDVI数据由美国航空航天局NASA提供的MODIS(moderate resolution image spectroradiometer)遥感产品(MOD13A3)(https://modis.gsfc.nasa.gov/)处理所得，已经过质量控制，该产品时长为2000年2月至2016年12月。MOD13A3是MODIS的三级产品，是在二级产品的基础上，对由传感器成像过程产生的边缘畸变(Bowtie效应)进行校正所得。该产品空间分辨率为1 km×1 km，时间分辨率为月，原始数据为EOS-HDF格式，有效值域范围为-0.2～1.0，数值越大说明植被覆盖状况越好。数据通过MRT软件批量处理成GEO TIFF格式，并通过ArcGIS投影拼接裁剪处理，其中合成方法选用最大合成法。植被类型数据来源于国家自然科学基金委员会中国西部环境与生态科学数据中心(http://westdc.westgis.ac.cn)，空间分辨率为1∶100万。

1.2.2气象数据

气象数据包括逐月降水和平均气温，均来自中国气象科学数据共享服务网(http://data.cma.cn/)的中国地面气候资料月值数据集。该数据集已通过质量控制，月年值各要素无错误数据。主要选取雅鲁藏布江流域内15个气象站点(当雄、拉孜、南木林、日喀则、尼木、拉萨、墨竹贡卡、琼结、泽当、江孜、嘉黎、波密、加查、林芝、米林)自2000年2月至2016年12月的逐月降水和平均气温数据，经去趋势处理后用于和NDVI进行偏相关分析。研究区气象站点分布如图1所示。

2 研究方法

2.1 偏相关系数法

偏相关系数[29]是考察多变量相关程度时采取的研究方法，是在控制其他变量的条件下，衡量其中某两个变量直接线性相关关系的指标。其中，一阶偏相关系数的计算公式为

(1)

式中：r12-3为控制第3个变量的情况下第1个变量和第2个变量的相关系数；r12为第1个变量和第2个变量的相关系数；r13为第1个变量和第3个变量的相关系数；r23为第2个变量和第3个变量的相关系数。r12-3越接近1，表明偏相关性越高。

图2 2000—2016年生长季平均NDVI空间分布

2.2 Hurst指数法

基于重标极差(R/S)分析方法的Hurst指数[30]是定量描述时间序列信息长期依赖性的有效方法，该方法是由英国水文学家Hurst对尼罗河进行的长期水文观测归纳发明，其基本原理是：给定一时间序列{ξt}，对于任意正整数τ，定义均值序列：

(2)

(3)

(4)

式中：〈ξ〉τ为给定时间序列的均值序列；Xt,τ为时间序列与其均值序列的累计离差；Rτ为累计离差的极差；Sτ为时间序列的标准差。若存在R/S∝τH，则说明时间序列存在Hurst现象，H值称为Hurst指数，其值可在双对数坐标系(lnτ, ln(R/S))中用最小二乘法拟合得到。Hurst指数主要有3种形式：①当0.5

2.3 一元线性回归法

一元线性回归法是根据自变量和因变量的相关关系，建立两者的线性回归方程对因变量进行分析的方法[31-32]。一元线性回归法的模型为

Yt=a+bΧt

(6)

(7)

(8)

式中：Χt为自变量时间序列；Yt为因变量时间序列；a和b为一元线性回归方程的参数；n为系列长度。当b>0时，自变量与因变量呈正线性相关；当b<0时，自变量与因变量呈负线性相关；当b=0时，自变量与因变量无线性相关关系。

3 雅鲁藏布江流域NDVI特征分析及其对气候变化的响应

3.1 植被覆盖空间分布特征

3.1.1NDVI空间分布特征

从基于最大值合成法计算的2000—2016年雅鲁藏布江流域生长季NDVI空间分布(图2)可以看出，NDVI高值区主要分布于流域中游和下游地区，尤其是下游地区，米林、波密等站附近NDVI值大多在0.8以上；而拉孜以上上游地区由于海拔较高，不少地区冰雪终年覆盖，罕有植被生长，故多为NDVI低值区，基本小于0.4。和图1对比可以发现高海拔地区NDVI值普遍较低，低海拔地区NDVI值保持稳定高值，说明植被覆盖分布特征与海拔存在有一定相关关系。

3.1.2NDVI沿高程变化特征

以250 m为分度值对雅鲁藏布江流域NDVI值沿高程变化特征进行分析(图3)。可以看出，NDVI值呈现一定程度的海拔分异性，随海拔升高，NDVI值显著下降。中低海拔地区为NDVI高值区，海拔高度为1 000～2 500 m范围内NDVI值基本保持同一水平，均值为0.91。2 500 m以上地区NDVI值逐步下降，在3 000 m附近呈微弱上升趋势之后急剧下降。6 000 m以上NDVI值较低，均值为0.03，且随海拔高度增加NDVI值基本保持不变。

图3 NDVI沿高程分布变化

图4 2000—2016年NDVI 10 a变化趋势空间分布

3.2 植被覆盖时间趋势特征

3.2.1NDVI年际变化分析

雅鲁藏布江流域狭长且空间差异性较大，不同地区NDVI年际变化差别显著。对研究时段内流域NDVI年际变化趋势进行空间分析(图4)，可以看出，流域内NDVI年际变化空间异质性显著，流域上游和下游地区年际NDVI变化趋势中增加比例分别占56.82%和58.60%，而流域中游地区NDVI变化趋势中增加比例仅占42.07%。流域上游和下游地区NDVI 10 a变化趋势在-1%～1%之间的比例分别为73.17%和58.70%，而流域中游地区NDVI 10 a变化趋势在-1%～1%之间的比例为51.72%。结果表明流域上游和下游地区植被覆盖以改善为主且绝大部分区域变化轻微，流域中游地区植被覆盖以恶化为主且部分区域恶化程度较重。雅鲁藏布江流域人类活动对下垫面植被变化的影响不容小觑，下一步可以联系土地利用变化等因素进一步探讨其驱动机制，在更长时间尺度上进行植被覆盖对气温和降水的时间滞后性分析，探究下垫面变化对湿度和蒸散量等更多气象要素的响应关系，并分析研究结果的不确定性，并继续补充数据源，例如采用SPOT VGT和AVHRR遥感数据，应用多源遥感数据时空融合技术研究下垫面长期时空演变特征。

分析2000—2016年雅鲁藏布江流域上游、中游、下游地区NDVI年均值变化(图5)，可以看出，NDVI总体呈上升趋势，增长率为0.002%/10a(显著性水平p<0.001)，表明近17 a流域植被覆盖总体呈缓慢好转趋势，尤其是中游和下游地区NDVI自2015年后出现显著增长。NDVI自上游、中游至下游呈依次增大趋势，NDVI均值分别为0.21、0.36和0.55，全流域NDVI均值为0.41，这与上中下游地区海拔差异有一定关系。上游大部分地区位于雪线以上或山谷边缘，坡度变化剧烈，植被生长稀疏；而下游地区海拔较低，气温较高，降水丰沛，在适宜的雨热条件下植被生长旺盛。因此，可以看出雅鲁藏布江流域植被分布因地形和海拔分布不均，且具有很大空间差异性。

图5 2000—2016年NDVI年际变化

3.2.2NDVI年内变化分析

对2000—2016年各月NDVI均值进行计算后的空间分布如图6所示，为了更加清楚地展现2000—2016年各月NDVI均值的最大值和最小值空间分布情况，制作最值统计图如图7所示。从图6和图7可以看出：2月NDVI低值区面积最大，中上游NDVI值均小于0.4，仅有流域下游植被覆盖面积较大，NDVI值维持在0.9左右；8月NDVI高值区面积最大，中上游NDVI值均在0.4至0.5之间，流域中游植被覆盖面积达年内最大值，维持在0.7左右。

图6 2000—2016年NDVI年内变化空间分布

图7 2000—2016年NDVI年内最值空间分布

分别统计研究区上游、中游、下游和全流域NDVI年内变化如图8所示。可以看出NDVI值随月份变化呈先增大、后减小趋势，全流域1—3月NDVI值维持在0.22左右，在3月达最小值0.21，之后稳步上升在7月达最大值0.48，9月之后下降较为明显，在12月下降为0.25。流域内各分区和全流域年内变化趋势基本一致，但下游地区NDVI值在1—3月呈明显下降趋势。

图8 2000—2016年NDVI年内变化

3.2.3NDVI季节变化分析

为了准确反映各季节植被生长情况，按照春季(3—5月)、夏季(6—8月)、秋季(9—11月)和冬季(12—2月)对雅鲁藏布江流域各季NDVI采用最大值合成法进行统计分析，结果见图9。从图9可以看出，2000—2016年间流域夏季NDVI值最大，均值达0.40，其中2013年的0.42和2010年的0.37分别为最大和最小值；其次为秋季，NDVI均值为0.36，在2000—2005年间呈缓慢上升趋势，并在2005年达最大值0.39，此后波动下降，在2008年达最小值0.33，2009年后基本维持在0.35左右；春季和冬季NDVI值相对较小，均值均为0.21。夏季和秋季NDVI呈下降趋势，10 a变幅均为-0.000 3%(夏季p=0.876，秋季p=0.897)，而春季和冬季NDVI呈上升趋势，10 a变幅分别为0.003%(p=0.131)和0.002%(p=0.097)。

图9 2000—2016年各季NDVI变化过程

3.2.4NDVI变化持续性分析

由上文分析可以发现雅鲁藏布江流域NDVI变化具有明显的周期性，因此有必要对该流域NDVI变化的持续性进行研究。通过对NDVI月值逐格点进行Hurst指数计算，结果空间分布如图10所示，相关统计信息如表1所示。

图10 2000—2016年逐月NDVI的Hurst指数空间分布

表1 雅鲁藏布江流域逐月NDVI持续性统计

从图10可以发现，NDVI的Hurst指数分布具有一定的空间差异性，在流域中游墨竹贡卡附近拉萨河流域及流域上游地区空间持续性较强，表明未来与过去变化趋势趋同概率较高。全流域Hurst指数均值为0.51，说明雅鲁藏布江流域NDVI呈弱持续性，未来与过去的变化趋势趋同概率较低。如表1所示，研究区呈弱持续性比例最大，达43.93%，流域面积为105 633 m2；其次是弱反持续性，面积为85 826 m2，占35.69%；强持续性与强反持续性所占比例最小，分别为12.69%和7.69%，所占面积分别为30 526 m2和18 498 m2。

3.3 不同植被类型NDVI时空变化特征

雅鲁藏布江流域植被类型数据来源于寒区旱区数据中心提供的1∶100万中国植被图集，通过对不同植被类型年际及年内NDVI进行叠加并采用最大值合成法进行统计分析，结果如图11所示。由图11可知，不同植被类型年际和年内分布从大到小依次为阔叶林、草丛、针叶林、栽培植被、灌丛、草甸、草原、高山植被和其他，均值分别为0.88、0.86、0.76、0.56、0.53、0.37、0.33、0.27和0.13。其中，2009年栽培植被和草原NDVI值较其他年份小，灌丛和高山植被NDVI值较其他年份大。植被非生长期栽培植被和草原NDVI值趋同，生长期开始后，4—7月草原和阔叶林NDVI值趋同，7—10月栽培植被和灌丛NDVI值趋同，这体现生长期草原和阔叶林、栽培植被和灌木分布差距较小的特点，草丛、灌丛相较于阔叶林、栽培植被等植被类型受植物节律影响，年内变化显著的特点。研究表明[20]高海拔地区不同高程植被覆盖与植被类型有关，通常在3 600 m左右植被覆盖开始出现急剧下降。对流层平均气温直减率为0.65℃/100 m，3 600 m以上绝大部分地区平均温度在0℃以下，植被类型从阔叶林针叶林向草原荒漠过渡可能是不同海拔高度植被覆盖度差异的原因。

图11 2000—2016年不同植被类型NDVI年际和年内变化过程

图12 NDVI与气象要素相关性分布

3.4 植被覆盖对气候变化的响应关系

在雅鲁藏布江流域植被与气象要素的相关性分析中，为了探讨区域NDVI与站点气象要素的关系，根据研究区气象站点位置，参考已有研究[33-37]，利用ArcGIS提取以各站点为中心，半径为1.5 km的区域NDVI平均值，并对2000—2016年NDVI、气温和降水数据进行去除趋势处理，然后分别计算NDVI与气温和降水的偏相关系数，结果如图12所示。NDVI与气温的偏相关系数在0.27～0.85之间，均值为0.47；而NDVI与降水的偏相关系数在-0.28～0.57之间，均值为0.28，且通过了显著水平为0.05的显著性检验。下游同中上游相比，NDVI和气温的偏相关系数较降水大，说明流域下游地区植被覆盖受气温影响更大。全流域NDVI与气温的平均偏相关系数大于降水，说明整体来说雅鲁藏布江流域植被生长受气温影响要强于降水。

4 结 论

a. 流域植被覆盖自上游、中游至下游呈递增趋势，海拔2 500 m以上植被覆盖开始下降，至6 000 m以上后保持在极低水平。流域上游和下游植被覆盖以改善为主且绝大部分区域变化轻微，流域中游植被覆盖以恶化为主且部分区域恶化程度较重。

b. 植被覆盖在2月低值区面积最大，8月高值区面积最大。NDVI夏季最大，秋季次之，春季和冬季相对较小。夏季和秋季NDVI呈下降趋势，而春季和冬季呈上升趋势。

c. 流域Hurst指数空间差异性显著，均值为0.51,说明流域植被覆盖整体呈弱可持续性，墨竹贡卡附近拉萨河流域及流域上游空间持续性较强。整体来说流域植被覆盖受气温影响要强于降水，下游最为明显。

d. 不同植被类型NDVI年际和年内分布从大到小依次为阔叶林、草丛、针叶林、栽培植被、灌丛、草甸、草原、高山植被和其他。生长期草原和阔叶林、栽培植被和灌木分布差距较小，非生长期栽培植被和草原分布较为相近。