2000-2014年塔里木河干流的植被覆盖与蒸散发时空变化及其关系
2018-07-26闫俊杰吕光辉徐海量徐新文凌红波
闫俊杰, 吕光辉, 徐海量, 徐新文, 凌红波
(1.伊犁师范学院, 新疆 伊宁 835000; 2.中国科学院 新疆生态与地理研究所, 新疆 乌鲁木齐 830011; 3.中国科学院大学, 北京 100049; 4.新疆大学 资源与环境科学学院, 新疆 乌鲁木齐 830046)
干旱与半干旱区绿洲系统所具有的适宜小气候及稳定的植被群落是保证其内部人类及其他生物活动所需环境的基础[1]，充足的水源供给则是维持绿洲环境稳定的关键[2]，然而，随着绿洲内部社会及经济的发展，人类生产生活用水不断挤占自然生态用水[3-4]，致使社会发展与环境稳定之间的矛盾日趋紧张[3-5]，威胁到绿洲系统良性发展，合理配置有限的水资源对绿洲自然与社会协同发展意义重大。蒸散发(evapotranspiration， ET)代表了植被及地面向大气输送的水汽总通量[6]，在绿洲系统耗水中占有重要比例，量化绿洲系统蒸发与植被变化及其相互关系对指导绿洲系统水资源调控与分配具有重要指导价值。
遥感方法是开展全球或区域尺度蒸散发与植被动态监测的有效方法[7-9]。目前已有多个世界机构和科研团队公开发布了不同时间及空间尺度的蒸散发与植被遥感的成品数据，其中蒸散发数据有欧洲气象卫星应用组织发布的LSA-SAF MSG ET数据集、Jung等[9]制作的1982—2008年全球陆地蒸散发数据集以及美国航空航天局(NASA)的MODIS MOD16全球ET数据集；由于与植被生产力具有良好相关性，归一化植被指数(normalized difference vegetation index，NDVI)是陆表植被动态监测的有效指标[9-10]，目前主要的成品数据有起始于1982年的NOAA-AVHRR NDVI，起始于1998年的SPOT-VGT数据以及起始于2000年的MODIS NDVI数据[10]。对于3种蒸散发数据，多位学者已对MOD16数据在亚洲、中国西北及新疆地区的精度进行验证[11-13]，证明了其有效性及广泛适用性。而MODIS NDVI数据则是植被动态监测中被广泛采用的数据[14-15]。
塔里木河干流位于新疆塔里木盆地北缘，地处塔克拉玛干沙漠边缘，干旱少雨，自然生态主要依靠河水及地下水来维持，但长期以来流域内自然生态用水与农业用水矛盾突出，耕地的持续扩张促使河水及地下水被过度引用和开采，生态用水被挤占，致使塔里木河水量锐减甚至断流[16]，生态环境受到严重损害[17]，自2000年起，为改善塔里木河流域的生态环境问题，国家及当地有关部门制定和实施了生态输水、统一管理等一系列生态恢复措施[18]，经过多年持续治理，塔里木河生态环境得到了明显改善[18]，目前，针对生态恢复措施对塔里木河生态环境的影响，大量学者就生态需水、植被变化等方面进行了深入而广泛的研究[14,19-20]，但关于生态恢复过程中蒸散发的变化规律及其与植被变化协同关系尚缺乏深入研究。基于此，本文拟利用2000—2014年MODIS NDVI及MOD16 ET数据产品，综合利用遥感及GIS空间分析技术，分析塔里木河干流植被覆盖及蒸散发时空变化，探讨植被覆盖变化对蒸散发的影响，以期为掌握塔里木河生态恢复的耗水规律，以及生态水的合理规划和配置提供相应参考。
1 研究区概况
塔里木河干流地处中纬度欧亚大陆腹地，远离海洋，北靠天山，南邻塔克拉玛干沙漠，受大陆性暖温带极端干旱气候控制，气候干旱少雨，蒸发强烈，多年年平均降水仅50.4 mm，年平均蒸发高达1 880.0 mm，极端最高温度39.4 ℃，多年平均气温10.7 ℃，干燥度12～19。该地区植被地带性属温性灌木和半灌木，但受外来河水补给影响，其植被群落多以胡杨、柽柳、芦苇等为建群种。
塔里木河干流通常指塔里木河阿拉尔断面至台特玛湖的河段，其下游的中大西海子—台特玛湖河段植被覆盖虽然经过多年输水治理后有了显著改善，但相对于阿拉尔—大西海子段，其植被仍非常稀疏，且植被带宽度及面积也均非常有限，因此截取植被覆盖相对较高的阿拉尔—大西海子段作为本文研究区域。
2 材料与方法
2.1 数据来源及预处理
本文用到的NDVI数据和ET数据分别为MODIS MOD13 Q1产品和MOD16 A3产品，两者空间分辨率分别为250和1 000 m，时间分辨率分别为16和8 d，时间序列为2000—2014年。MODIS ET数据由美国蒙大拿大学森林学院工作组在Penman-Monteith公式基础上融合MODIS遥感数据制作而成[21]，计算过程考虑了符合绿洲和荒漠下垫面特征的土壤表面蒸发、冠层截流水分蒸发和植物蒸腾的差异，适用于干旱与半干旱区地表蒸散发研究[9]。MODIS NDVI数据是的NOAA/AVHRR NDVI的延续和发展[22]在红光波段(RED)和近红外波段(NIR)波段的波幅以及辐射定标技术上均有较大改进[23]。
对获得的遥感数据除进行了数据格式转换、镶嵌、投影转换及研究区提取等预处理外，对NDVI数据进行了Savitzky-Golay滤波和MVC合成处理，以降低噪音信息对数据影像，并获得代表植被生长最好状况的年NDVI数据；对于ET数据，由于塔里木河流域距离河道较远的区域植被覆盖较低，容易造成研究区内ET数据出现空缺值，因此对ET数据中的空缺值进行了8领域均值插值。最后为保证ET数据与NDVI数据的空间匹配，将两种数据像元重采样为250 m×250 m。
2.2 趋势检验
用Mann-Kendall检验与Theil-Sen median趋势分析判断NDVI及ET时间序列数据变化趋势，并量化其变化率[24-25]。Mann Kendall检验统计量计算过程为：
(1)
(2)
(3)
(4)
式中：xk，xi——NDVI或ET样本时序数据集合；n——数据集合长度； sign——符号函数。在给定显著性水平α下，当|Zc|>U1-α/2时，表示时间序列在α水平上变化趋势显著。若统计量Zc值为正，表示变化趋势上升，若Zc值为负，则表示变化趋势下降。
Theil-Sen median趋势分析用于量化变化趋势，表示单位时间内的变化量，计算公式为：
(5)
式中：1
3 结果与分析
3.1 NDVI及ET整体特征
以2000—2014年多年平均NDVI和ET代表研究区植被覆盖和ET的总体特征，并根据研究区NDVI及ET直方图的特征，分别将NDVI和ET划分为了5个等级(表1—2)。此外，还根据研究区ET和NDVI的空间分布特征，将研究区以十四团和恰拉2个断面划分为了上(阿拉尔—十四团)、中(十四团—恰拉)和下(恰拉—大西海子)3段，便于更为详尽的分析研究区NDVI及ET的空间分异特征。
表1 研究区各河段及全区平均NDVI及NDVI等级比例 %
表2 研究区各河段及全区平均ET及ET等级比例 %
研究区NDVI与ET的空间分异一致，阿拉尔—十四团的上段均为NDVI和ET高值的集中分布区，中段中的十四团—新渠满及英巴扎—恰拉河段则均为NDVI和ET低值的集中分布区，恰拉—大西海子的下段及中段的新渠满—英巴扎为NDVI和ET中值的集中分布区。研究区NDVI与ET空间分异的一致性在一定程度上也表明了ET的大小与植被覆盖的高低具有较高的相关性。
根据表1和表2中的统计结果，研究区全区多年平均NDVI为0.33，69.34%的区域的NDVI<0.40，仅有0.72%的区域NDVI>0.80，而NDVI介于0.60～0.80的面积比例也仅达到11.50%，剩余18.44%的区域的NDVI介于0.40～0.60。研究区全区多年平均ET仅为118.41 mm，与其干旱缺水的环境特性一致；全区69.50%的区域ET<120 mm，ET>200 mm的面积比例仅为3.88%，剩余7.76%和18.86%的区域的ET分别介于160～200 mm和120～160 mm。
对于不同河段，上段的NDVI及ET均最高，分别为0.53，157.05 mm，分别高出其全区平均值60.61%和32.63%，该河段47.32%区域的NDVI>0.60，48.97%区域的ET>160 mm。中段因79.91%面积的NDVI<0.40，80.72%面积的ET<120 mm，其平均NDVI和ET均最低，分别为0.28，108.94 mm;下段NDVI和ET分别为0.36，121.05 mm，均高于全区平均值，虽然其大面积的NDVI<0.40及
ET<120 mm，但其39.01%面积的NDVI>0.40，40.40%面积的ET>120 mm。
3.2 NDVI及ET变化趋势
由图1和表3可知，除了2007年外，NDVI与ET呈现相同的年际波动，均于2003，2008及2012年出现相对峰值，于2001，2009年出现相对谷值。但对研究区全区平均NDVI和平均ET的趋势分析显示，NDVI呈显著增加趋势，而ET虽然变化趋势未达到显著水平但却呈下降趋势，两者变化率分别为0.003单位/a和-0.03 mm/a (表3)。
图1 研究区全区NDVI和ET年际变化
表3 研究区全区及不同河段NDVI和ET变化趋势
对于不同河段，上段NDVI和ET一直为最高，其次为下段，中段一直最低(图2—3)。不同河段NDVI和ET大小虽然差异较大，但其年际波动则基本相同，且与全区平均NDVI和ET年际波动一致。对不同河段NDVI和ET的变化趋势分析(表3)发现，3个河段之中上段不仅植被覆盖最高，其NDVI和ET变化速率也最大，分别为0.007单位/a，1.91 mm/a，均远高于其他2个河段，其中NDVI呈显著增加变化趋势，ET也呈增加趋势，但增加未通过α=0.05显著检验；中段植被覆盖最低，其NDVI和ET变化速率也最小，分别为0.002单位/a和-0.48 mm/a，其NDVI为非显著增加，ET则为非显著减小；下段NDVI和ET变化趋势分别为显著增加和非显著减小，其变化量分别为0.004单位/a和-0.52 mm/a。
图2 研究区各河段NDVI年际变化
从上述分析可以发现，无论是全区平均还是各河段，其NDVI与ET呈现几乎一致的年际波动，但除植被覆盖较高的上段NDVI与ET呈现相同变化趋势外，其他2个河段却呈现相反的年际变化趋势，在一定程度上表明虽然年ET总量与NDVI存有较高的相关性，但NDVI年际变化并非是决定ET年际变化的控制因素。
图3 研究区各河段ET年际变化
3.3 NDVI及ET变化的空间特征
根据Mann-Kendall检验与Theil-Sen median趋势分析计算公式，对研究区NDVI和ET逐像元计算，获得NDVI和ET Mann-Kendall检验统计量Zc值及年际变化速率β值空间分布数据。依据此数据，对NDVI和ET的变化趋势及变化速率进行分类和划等，统计其比例特征(表4—6)。
2000—2014年NDVI和ET的变化趋势存在较大的差异。根据表4中统计结果，15 a内全区内51.17%区域的NDVI有所增加，而ET发生增加的区域仅为29.43%，NDVI发生减少的比例为48.83%，ET发生减少的比例则高达70.57%。但同时也应该看到，全区绝大部分区域NDVI和ET的增加或减少并未达到显著水平，其中NDVI的比例为58.82%，而ET的比例则为76.39%。
对于不同河段，各个河段NDVI与ET的变化趋势的空间分布均存在较大差异。上段是NDVI和ET增加的集中区域，根据表4中的统计结果，其NDVI与ET增加和减少的面积比例虽相对一致，分别为74.15%和73.35%，25.85%和26.65%，但NDVI显著增加的比例高达41.27%，而ET的该比例却为24.62%，远低于NDVI；同时NDVI与ET非显著增加的比例也有较大差别，分别为32.88%和48.73%。中段和下段NDVI与ET变化趋势空间分布的差异均远高于上段，两河段NDVI显著增加的比例分别达到17.29%和24.49%，而ET显著增加的比例的则仅为2.21%和1.61%，同时，两个河段NDVI和ET非显著减少和显著减少的差异也很大，而显著减小的差异则均较为一致。
表4 研究区全区及各河段NDVI和ET不同变化趋势比例 %
表5 研究区全区及各河段NDVI不同变化速率等级比例 %
表6 研究区全区及各河段ET不同变化速率等级比例 %
对于研究区NDVI和ET变化率的空间特征，从表5—6可知，研究区NDVI的变化率主要介于-0.005～0.005单位/a，其比例为55.34%，变化率<-0.005单位/a及变化率>0.005单位/a的比例分别为17.42%和27.23%。研究区ET的变化率主要介于-0.20～0 mm/a，其比例为66.58%，变化率介于0～0.20 mm/a及变化率>0.20 mm/a的比例分别达到了19.05%和10.68%，而ET变化率<-0.20 mm/a的比例仅为3.69%；不同河段，上段是NDVI及ET增加速率最快的集中区域，其NDVI变化速率>0.005 mm/a和ET变化速率>0.20 mm/a的比例分别达到了49.84%和41.70%。中段NDVI和ET减少的集中区域，其NDVI变化速率<-0.005单位/a的比例达到了20.05%，其ET变化速率主要介于-0.20～0 mm/a，比例为75.33%。下段NDVI主体表现为增加，呈增加的比例达到了65.95%，ET主体表现为减少，减少速率主要介于-0.20～0 mm/a，其比例为72.40%。
4 讨 论
4.1 植被覆盖与蒸散发关系
利用2001—2015年研究区NDVI和ET的15 a平均空间数据以及其变化率空间数据建立其散点图(图4—5)，并将NDVI和ET变化率按0.01个单位及0.002 5单位/a的间隔划分区间，计算每个区间内所有像素点ET平均值和ET变化率的平均值，绘制ET及ET变化率随NDVI及NDVI变化率的变化曲线和其趋势线，分析NDVI与ET及两者变化的相关关系。
图4 研究区15 a平均NDVI和平均ET
对于研究区ET与NDVI的关系，由图4可知，研究区ET随NDVI的增加逐步增加，两者拟合曲线的R2达到了0.71，同时经相关分析，两者相关系数也达到0.81(p=0.001)，即植被覆盖高的区域也是ET的高值区域，这也与NDVI和ET的空间分布格局一致，也证明了植被覆盖是决定研究区ET空间分布的重要因素。此外，袁国富等[26]通过对塔里木河干流下游(大西海子—台特玛湖段)涡度相关观测数据的分析也表明，蒸散的空间格局受到植被叶面积指数的控制，植被盖度越大，蒸散量越大，与本文的研究结果一致。
图5 研究区NDVI和ET变化率
而对于NDVI与ET年际变化的关系，从图5中可以看出，ET变化率随NDVI变化率的增加而逐步增加，但两者拟合关系远低于图4中NDVI与ET的拟合关系；同时相关分析表明，研究区2001—2014年全区平均年NDVI和年ET的相关系数为0.49(p=0.067)，NDVI变化速率与ET变化速率的相关系数为0.59(p=0.001)，均低于图4中NDVI与ET的相关系数。此外，根据表7中统计结果，研究区NDVI和ET变化趋势一致的面积比例为69.08%，变化趋势相反的面积比例为30.92%，其中26.73%的区域为共同增加，42.35%为共同减少，26.35%区域NDVI呈增加趋势而ET则呈减少趋势，4.57%区域NDVI呈降低趋势而ET呈增加趋势。可见，研究区植被覆盖并非是ET年际变化的决定因素。
表7 研究区ET变化率等级与NDVI变化率等级比例矩阵
4.2 ET年际变化影响因素
张巧凤等[27]通过对锡林郭勒草原ET年际变化分析表明，ET的年际变动态是气候及植被覆盖等多种因素共同作用的结果，ET变化与降水量、NDVI和水汽压呈极显著正相关关系，而与气温呈负相关关系；赵燊等[28]对山东省ET的研究也表明ET的时空变化与诸多气象因子相关，其中与降水及温度的关系最为密切；代超[29]，蹇东南等[30]通过对塔里木河整个流域长时间序列的数据分析表明，除了气候变化因素外，径流量及耕地面积的增加也是该区域ET时空变化的重要因素。基于前人研究，本文对ET与年平均气温、年平均水汽压、年降水量及年径流量进行了相关分析，发现相对其它要素，ET年际变化与降水和气温相关性更好(表8)。
表8 塔里木河ET与各要素的相关系数及其p值
刘波等[31]的研究认为，对于干湿条件不同的气候区，影响ET变化的主导因子存在明显差异，而对于中国西北干旱区，影响ET变化的主要因子是供水条，蹇东南等[30]的研究也发现出山口径流量减少是致使1997—2013年塔里木河流域ET减少的重要因素，然而代超[29]的研究表明塔里木河流域ET变化是灌溉引水及气候变化共同作用的结果。同时，本文分析发现ET与径流量相关系数却仅为0.34(p=0.209)，相关性差，除此之外，自2000年开始，塔里木河实施综合规划管理[32]，在耕地扩张、生态闸截水/调水及补给地下水等人为干扰因素的影响下，塔里木河下垫面特征及水资源的时空分布发生较大改变[33]，必将在很大程度上改变ET与植被、径流及气象等因素的时空匹配特征，因此量化人为因子的变化特征，分析其对ET的贡献率将一进加深对塔里木河干流ET时空变化规律的认识。
5 结 论
本文利用MODIS的ET和NDVI数据，及Mann-Kendall检验与Theil-Sen median趋势分析方法，选择阿拉尔—大西海子段作为代表区域，分析了2000—2014年塔里木河干流植被覆盖及ET时空变化规律，并探讨了两者的相互关系及ET变化的影响因素，得出如下结论：
(1) 研究区多年平均NDVI及ET分别为0.33，118.41 mm；全区69.34%的面积的NDVI<0.40，69.50%的面积的ET<120 mm。空间上ET及NDVI空间分异总体一致，均表现为上段>下段>中段。
(2) 研究区NDVI与ET平均变化趋势相反，NDV显著增加(ZC>1.96)，而ET则表现为非显著下降(-1.96 (3) 研究区NDVI与ET的变化趋势及变化速率均存在较大空间差异，全区NDVI发生降低的比例为48.83%，而ET的该比例则高达70.57%；全区NDVI变化速率主要介于-0.005～0.005单位/a，ET的变化速率主要介于-0.20～0 mm/a；NDVI及ET增加的区域集中分布在上段，而NDVI和ET减少的区域则集中分布在中段。 (4) ET的空间分布受植被覆盖控制，但在年际变化上，ET与NDVI相关性较低，且全区30.92%的面积NDVI和ET呈现相反的变化变化趋势。相对于NDVI，径流及水汽压，ET的年际变化与气温和降水相关性更好。
