閆俊杰, 呂光輝, 徐海量, 徐新文, 凌紅波
(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空間分析技術,分析塔里木河干流植被覆蓋及蒸散發時空變化,探討植被覆蓋變化對蒸散發的影響,以期為掌握塔里木河生態恢復的耗水規律,以及生態水的合理規劃和配置提供相應參考。
塔里木河干流地處中緯度歐亞大陸腹地,遠離海洋,北靠天山,南鄰塔克拉瑪干沙漠,受大陸性暖溫帶極端干旱氣候控制,氣候干旱少雨,蒸發強烈,多年年平均降水僅50.4 mm,年平均蒸發高達1 880.0 mm,極端最高溫度39.4 ℃,多年平均氣溫10.7 ℃,干燥度12~19。該地區植被地帶性屬溫性灌木和半灌木,但受外來河水補給影響,其植被群落多以胡楊、檉柳、蘆葦等為建群種。
塔里木河干流通常指塔里木河阿拉爾斷面至臺特瑪湖的河段,其下游的中大西海子—臺特瑪湖河段植被覆蓋雖然經過多年輸水治理后有了顯著改善,但相對于阿拉爾—大西海子段,其植被仍非常稀疏,且植被帶寬度及面積也均非常有限,因此截取植被覆蓋相對較高的阿拉爾—大西海子段作為本文研究區域。
本文用到的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。
用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
以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。
由圖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年際變化
根據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%。
利用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變化率等級比例矩陣
張巧鳳等[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時空變化規律的認識。
本文利用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的年際變化與氣溫和降水相關性更好。