川西高原植被时空格局及其对气候变化的响应

伍良旭, 王 晗, 邵怀勇, 周思彤

(成都理工大学 地球科学学院, 成都 610059)

植被是地球生态系统主要的组成部分，是连接气候与土壤的纽带，在全球物质能量交换和气候变化中具有不可替代的作用[1]。植被动态变化在高原生态环境中具有敏感的指示作用，合理利用区域资源和有效治理生态环境提供依据，已成为生态环境研究的热点[2]。NDVI是表征植被动态变化的一个重要指标，能较精确地反映地表植被的生长状态及演变特征，在区域植被退化、湿地景观格局等领域发挥着重要的作用[3]。川西高原位于我国青藏高原向川西南山地和四川盆地的过渡地区，是我国生态安全战略格局中青藏高原生态屏障和黄土高原—川滇生态屏障的重要组成部分，地形复杂，且受高原背风坡和西南季风等系统的共同影响，是典型的生态环境脆弱区和生态气候敏感区[2]，但目前针对该区域植被时空格局与气候因子的相关研究较少，因此对该地区的植被时空格局及其对气候变化响应的研究具有重要意义。

此外，现在的遥感数据通常不能同时具备较长的时间序列和较高的空间分辨率，如GIMMS NDVI3g是目前最长时间序列的NDVI数据集之一[4]，前人利用其进行了1982—2009年东北地区[5]、1982—2006年蒙古高原[6]和1982—2013年新疆地区[7]的植被变化研究，但8 km的空间分辨率应用在区域植被变化研究中过于粗糙。虽然如今高分辨率辐射计(AVHRR)、地球观测系统(SPOT)-VEGETATION，MODIS和ENVISAT中分辨率成像光谱仪(MERIS)等卫星传感器所采集的数据有较高的时间分辨率，通过合成一定时间内的图像，能够减少甚至消除云污染[4]，遥感数据质量不断提高，使得区域植被趋势分析在全球范围内得到广泛应用[8-10]，但如何提高数据的空间分辨率仍是亟待解决的问题。相对于GIMMS NDVI3g数据精度更高、更能表达植被时序变化的MODIS NDVI产品，却因其时间覆盖范围(2001年7月至今)无法针对长时间序列进行研究。因此，如何提高早期遥感数据的分辨率对于长时间序列的植被变化监测研究显得尤为重要。

针对遥感数据时空尺度的融合问题，前人做了大量的研究。如Gao等[11]提出的STARFM模型(spatial and temporal reflection fusion Model)可融合Landsat和MODIS地表反射率，生成具有高时空分辨率的合成反射率数据。Singh等[12]利用STARFM融合模型充分提取了农作物8 a生长季节的物候。Bhandari等[13]和Schmidt等[14]评估了STARFM合成长期时间序列数据的能力，并获取不同植被的物候特征进行研究。Meng等[15]和Zhu等[16]对STARFM算法改进，分别提出了时空融合自适应植被指数模型(STAVFM)和增强型自适应时空融合模型ESTARFM，以更准确地表征地表的时空变化。由于在地形复杂地区，ESTARFM的融合效果比STARFM较好，使得ESTARFM算法被频繁用于生成高时空分辨率的合成遥感图像[17-18]。然而，目前较少有学者针对川西高原使用ESTARFM模型将GIMMS NDVI3g与MODIS NDVI进行时空融合并开展1995年至今植被时空格局的相关研究，为此本文基于ESTARFM模型开展川西高原GIMMS NDVI3g和MODIS NDVI数据的时空融合，这对开展其长期的植被变化监测具有重要意义。

本研究尝试基于ESTARFM模型建立川西高原GIMMS NDVI3g与MODIS NDVI数据的函数关系，以提高GIMMS NDVI3g数据的空间分辨率。在此基础上分析川西高原长时间的植被变化特征及其对气候变化的响应，为当地的环境保护决策提供基础数据和科学参考。

1 研究区概况

川西高原(27.59°—34.31°N，97.36°—104.62°E)地属青藏高原与四川过渡地带，覆盖阿坝州和甘孜州等高原藏区，区域面积约为25.38万km2。研究区内地势高，地形起伏大，平均海拔约4 000 m。境内属亚热带高原季风气候，地域差异明显，多年降水量约435.6～1 584.1 mm。近年来，由于全球气候变化以及人类活动的增加，当地植被覆盖逐渐变化，植被的退化对境内生态系统安全会造成极大威胁[19]，因此，研究分析川西高原植被覆盖时空分布特征，并掌握主要气候因素对植被变化的影响非常必要[20]。

2 数据与方法

2.1 数据来源与处理

(1) 遥感数据：MODIS 13Q1植被指数数据集产品，空间分辨率500 m×500 m，时间分辨率30 d。首先进行Savitzky-Golay滤波处理降低异常值等噪声影响，随后利用最大值合成法合成2005年和2015年MODIS NDVI数据。该方法能减少云雾等大气效应对数据的影响，提高数据精度[21]。

1995年、2005年和2015年的GIMMS NDVI3g数据空间分辨率为8 km×8 km，时间分辨率为15 d。采用最大值合成法合成3个时间段的年GIMMS NDVI3g数据，并将其空间分辨率重采样为500 m×500 m，为后续使用ESTARFM模型进行数据融合做准备。

(2) 气象数据：年均气温、年降水量空间插值数据集产品。薄板光顺样条法采用平滑样条函数对多变量数据进行分析和插值，即使用函数逼近曲面，是一种较好的气候要素插值的方法[22]。因此，该数据集采用基于薄板光顺样条法开发的ANUSPLIN插值软件(Hutchinson，2001)对气候要素进行插值，其中将气象站点的经度、纬度和海拔高程作为年降水量空间插值的自变量，而年均温度空间插值则考虑了温度随垂直高度的变化[22]。

MODIS NDVI数据来源于中国科学院计算机网络信息中心地理空间数据云平台(http:∥www.gscloud.cn/)；GIMMS NDVI3g数据来源于ECOCAST网站(https:∥ecocast.arc.nasa.gov/data/pub/gimms/3 g.v1/)；气象数据和高程数据来源于中国科学院资源环境科学数据中心(http:∥www.resdc.cn)；为提高试验结果精度，转换所有数据到统一的投影坐标系WGS_1984_Zone_48N。

2.2 研究方法

2.2.1 ESTARFM算法 为生成高时空分辨率的遥感图像，使用ESTARFM模型融合MODIS NDVI和GIMMS NDVI3g数据。ESTARFM模型是在STARFM模型[11]基础上根据像元的异质性调整赋权方法，通过设置转换系数来改进预测结果，异质性越大，预测精度越高，可以保存更多的空间特征细节[23]。

本文根据第1时期m时刻MODIS NDVI与tp时刻的GIMMS NDVI3g得到融合后tp时刻的GIMMS NDVI3g数据，记为mm=(Xw/2,yw/2,tp,B)，或利用第2时期n时刻的观测数据得到融合后的tp时刻GIMMS NDVI3g数据，记为mn=(Xw/2,yw/2,tp,B)。通过tm和tn与时刻tp间检测到的GIMMS NDVI3g变化幅度来计算权重，如公式(1)所示[16]：

(k=m,n)

由公式(2)计算中心像元的值：

M(xw/2,yw/2,tp,B)=Tm·Mm(xw/2,yw/2,tp,B)+

Tn·Mn(xw/2,yw/2,tp,B)

(2)

由于篇幅限制，具体融合原理请见参考文献[16]。统计2005年和2015年融合后的GIMMS NDVI3g和MODIS NDVI对应的像元值绘制散点图(图1)。融合的GIMMS NDVI3g数据与MODIS NDVI数据具有高度相似性，各像元值大致相同，两者基本分布在y=x的趋势线上(p<0.01)，预测效果较好。

图1 真实MODIS NDVI与融合后的GIMMS NDVI3g数据散点图

对比2005年和2015年GIMMS NDVI3g、融合后的GIMMS NDVI3g和真实MODIS NDVI数据细节(图2)可以看出，融合后的GIMMS NDVI3g与真实MODIS NDVI数据像元分辨率较好，地物纹理清晰，部分特征地区的细节能够凸现出来，能较好地再现区域植被覆盖细节特征。

注：A为2005年GIMMS NDVI3g；B为2005年融合后的GIMMS NDVI3g；C为2005年MODIS NDVI；D为2015年GIMMS NDVI3g；E为2015年融合后的GIMMS NDVI3g；F为2015年MODIS NDVI。图2 融合数据对比图

2.2.2 一元线性回归模型 基于融合处理后的GIMMS NDVI3g数据(值域为[-1,1])，使用像元二分模型计算植被覆盖度，由于篇幅限制，具体的植被覆盖度计算公式和原理请见参考文献[24]。

植被覆盖度变化主要用于分析植被覆盖度随时间变化的关系，他们的关系满足一元线性回归模型的条件，以局部像元时间变化特征表示整个空间的变化规律[24]，其计算方式见公式(3)：

(3)

式中:n为总监测年数;fi为第i年植被覆盖度;Slope为植被覆盖度线性拟合斜率。利用两者的相关关系来判断植被覆盖度时间变化的显著性,斜率为负表示植被覆盖度下降,反之则表示植被覆盖度上升,拟合斜率绝对值越大,变化越明显[24-25]。

2.2.3 偏相关性分析法 地理系统是一个由多要素组成的复杂系统，偏相关分析可剔除其他变量的影响，只研究一个因素与另一个因素的相关程度[25]。本文基于像元对川西高原1995年、2000年、2015年植被覆盖度与年降水量、年均气温两个气候因素进行偏相关性研究，见公式(4)[25]：

(4)

式中:rab,c表示控制变量c的值不变,变量a，b间的偏相关系数。rab，rbc，rbc分别表示变量a和b，b和c，a和c的简单相关系数。其中a，b和c分别为植被覆盖度、年均温度和年降水量。rab,c>0时表示两变量呈正相关,rab,c<0时表示两变量呈负相关,rab,c的绝对值越大,两变量的相关性越大。

简单相关系数计算见公式(5)[25]：

(5)

3 结果与分析

3.1 植被覆盖特征分析

3.1.1 植被覆盖空间分布特征 参照前人研究[26]将川西高原植被覆盖度分为极低覆盖度(0～0.2)、低覆盖度(0.2～0.4)、中覆盖度(0.4～0.6)、中高覆盖度(0.6～0.8)和高覆盖度(0.8～1.0)5个等级，得到川西高原1995年、2005年和2015年植被覆盖度空间分布图(图3)。

图3 不同等级植被覆盖度空间分布

分析可得，1995年、2005年和2015年研究区极低、低和中覆盖度区域主要分布在西部高海拔的白玉、巴塘、理塘和东部低海拔的泸定、康定、九龙、理县等县区，南部的稻城、和西北部的玛多县主要为中高覆盖度，其余地区多属于高覆盖度区域。西部地区海拔较高，多为高山植被，西部德格、甘孜一带和东部泸定县区有荒漠和草甸区，其余地区水热条件良好，主要为针叶林和阔叶林。2005年川西高原高覆盖度区域面积相对1995年有所减少，中植被覆盖和中高植被覆盖区域面积增加。2015年高植被覆盖区域覆盖面积继续降低，但趋势减缓，中、中高植被覆盖区域面积增加，且后者相对前者增加明显。

川西高原整体的植被覆盖较高，呈现较为明显的区域差异。高植被覆盖区主要分布在北部和南部地区，低值区分布在中部地区，此外沿海拔梯度，植被覆盖呈先升高后降低的垂直分布特点。

3.1.2 植被覆盖变化趋势特征 为研究1995—2015年川西高原植被覆盖的时间变化规律，对植被覆盖与时间进行一元线性回归分析(图4)，根据趋势变化的方向和程度，划分为5个等级(表1)[24]。植被覆盖呈增长趋势的面积占17.63%，其中明显增加面积占15.53%，轻度增加面积占2.10%，主要分布在高海拔地区的玛多、石渠、白玉等县区，以玛多县分布尤为密集。植被覆盖呈减少趋势的面积占27.36%，其中明显减少面积占18.29%，轻度减少面积占9.07%，减少区域广泛分布在各个区域，以理塘、康定县区分布较广，结合前人研究分析认为[3]，理塘县自20世纪80年代开始，因过度挖采、放牧等人为因素、鼠虫害、气候变化、暴雨冲刷等多种因素综合并长期作用使理塘县草原退化、沙源裸露导致植被覆盖减少；康定县由于乱砍滥伐，开荒种地，旅游开发等人为干扰因素和自然灾害等多因素引起大面积土地退化，山坡地水土流失十分严重，植被减少趋势明显。植被覆盖基本稳定的区域面积占55.01%，分布广泛，主要在研究区的东北部的阿坝、若尔盖县区，该地区主要为沼泽、草甸等，植被生长状况稳定，植被变化不明显。总体而言，1995—2015年川西高原植被覆盖呈增长趋势面积占17.63%，呈增长趋势面积占27.36%，基本稳定面积占55.01%，植被整体变化趋势稳定。

图4 川西高原1995-2015年植被覆盖变化趋势

表1 1995-2015年川西高原植被覆盖变化面积及比例

3.2 植被对气候变化的响应

3.2.1 年均气温与年降水量空间分布特征 植被的变化受到多方面因素的影响，其中年降水量和年均气温是对植被变化影响较大的两个气候因素[25]。计算1995年、2005年和2015年的年均气温和年降水量(图5)。

图5 1995-2015年年均气温和年降水量分布

分析可得，1995年、2005年和2015年川西高原年均气温为3.3℃，最高温度为17.9℃，位于最南部的盐源县，最低温度为-13.1℃，位于西北部的玛多县；年均降水量为848.6 mm，最大降水量为1 584.1 mm，位于西南部的泸定县；最小降水量为435.7 mm，位于西北部的玛多县。研究区内年均气温和年降水量的分布具有空间异质性，玛多县所在的西北部和泸定县所在的东南部区域气温整体低于中部地区，泸定县所在的东南地区，位于青藏高原和四川盆地过渡区，受东南季风和西南季风影响，降水充沛，气温不高；玛多县地区靠近青藏高原，从研究区中部到西北部地势逐渐增高，气温逐渐下降。

3.2.2 植被覆盖对气候变化的响应 计算1995—2015年年均气温和年降水量与植被覆盖的偏相关系数(图6A—B)。年植被覆盖与年均气温的偏相关系数变化范围为[-0.71，0.71]，年植被覆盖与年降水量的偏相关系数变化范围为[-0.71，0.71]。统计分析可得年植被覆盖与年均气温呈正相关的区域面积占14.9%，年植被覆盖与年降水量呈正相关的区域面积占35.4%。根据年均气温和年均降水量与年植被覆盖的偏相关系数方向，对研究区域进行分区讨论[27]，年植被覆盖与年均气温负相关，与年降水量负相关，归于分区一；年植被覆盖与年均气温负相关，与年降水量正相关，归于分区二；年植被覆盖与年均气温正相关，与年降水量负相关，归于分区三；年植被覆盖与年均气温正相关，与年降水量正相关，归于分区四；年植被覆盖与年均气温和年降水量的偏相关系数较低(都为0)，归于分区五(图6C)。

图6 1995-2015年年均气温和年降水量与植被覆盖度偏相关分析及偏相关性分区

分区一占总面积的6.4%，主要分布在石渠县、甘孜县和色达县的北部，植被类型主要为高寒草甸，分析认为该区域鼠灾泛滥和人为过度放牧导致植被覆盖稀疏，气温和降水反而加速了水土流失，抑制植被生长。分区二占总面积的30.8%，主要分布在西部的白玉县、理塘县和东北部的黑水县等地，植被类型主要为草甸、常绿阔叶林、针叶林和灌丛，地形复杂，气温较高而降水量低，水土蒸发量大，水土流失较严重[28]，气温升高减弱其光合作用，抑制植被的生长，适当的降水可增加土壤含水量，促进植被生长。分区三占总面积的10.3%，主要分布在马尔康县、道孚县、九龙县一带，该区域气温相对较低而降水量高，降水持续增加引发泥石流、滑坡等地质灾害，并使植被光照条件下降，减弱其光合作用，抑制植被生长，从而使植被覆盖降低[29-30]。分区四占总面积的4.6%，主要分布在西北部的玛多县及东南部的汶川县和康定县地区。汶川地区自2008年汶川地震后，该区域的地质活动较为频繁，生态系统也在逐渐恢复，气温和降水作为植被生长的两个重要因素，促进该区域生态重建[31]。玛多县位于高海拔地区，气温低、降水量少、人类活动影响程度小，气温和降水成为植被覆盖变化的主要诱因[32]。分区五占总面积的43.8%，主要分布在若尔盖县和红原县等东北部高寒高纬度地区，土地沙化严重，植被变化的主导驱动因素是人类活动，因此，该区域的年植被覆盖变化与年均气温和年降水量的相关性不大[33]。

在研究区不同的区域，年植被覆盖与年均气温的偏相关性和年植被覆盖与年均降水量的偏相关性的方向与程度不同，整体而言，植被覆盖对降水的变化相较于气温更加敏感，此外研究区内年植被覆盖与年降水量呈负相关的区域，往往与年均气温呈正相关，这与仙巍等[20]研究川西高山高原过渡带植被对气候变化响应的结果相符。

4 结 论

(1) 运用ESTARFM模型融合的2005年和2015年的GIMMS NDVI3g数据与对应的真实MODIS NDVI数据吻合度较高，能凸显地区特征，融合结果可用于地形复杂的川西高原植被变化的研究，为早期的大尺度植被变化研究提供了一种新的思路。

(2) 1995年、2005年和2015年川西高原植被覆盖结构整体较好，空间上呈现沿海拔先升高后降低的垂直分布特点，时间上植被变化趋势稳定。植被覆盖度较低的区域主要分布在低海拔区，近20 a来研究区植被覆盖度整体变化趋势比较稳定，其中以高覆盖度区域面积的逐年减少和较低覆盖度区域面积的逐年增加趋势较为明显。

(3) 1995年、2005年和2015年川西高原年均气温和年降水量存在空间异质性，气温空间差异表现为研究区北部气温低于南部、西部气温低于东部，降水空间差异表现为研究区西部降水量低于东部。

(4) 年植被覆盖与年均气温呈正相关的区域面积占14.9%，年植被覆盖与年均降水量呈正相关的区域面积占35.4%，此外研究区内年植被覆盖与年降水量呈负相关的区域，往往年植被覆盖与年均气温呈正相关，植被覆盖度对降水的变化相较于气温的变化更加敏感。

本研究中使用的数据时间间隔较长，可能会忽略一些细节。此外，由于数据获取等原因，仅选取了年均气温和年降水量两个主要的气候要素与植被进行响应分析，未考虑地表因素和人为因素等其他因素对植被的响应。未来的研究将会考虑更短的时间间隔和更多因素对植被变化的影响研究。