劉長雨,謝保鵬*,楊 潔,陳 英,裴婷婷
(1.甘肅農業大學管理學院,甘肅 蘭州 730070;2.甘肅農業大學草業學院,甘肅 蘭州 730070)
全球氣候變暖和人類活動的共同影響下,地球正以一定的速度綠化[1-2]。美國航天局衛星數據表明,全球從2000—2017 年新增的綠化面積中,約1/4來自中國,貢獻比例居全球首位[3]。青藏高原作為我國乃至亞洲的重要生態安全屏障[4],其獨特的地理位置和典型的高寒生態系統使其成為國內外學者研究的熱點區域。過去的50年里,青藏高原氣候變化的暖濕化特征明顯[5],溫度每10 年升高0.3℃~0.4℃,降水每10 年增加2.2%[6]。降水增加和溫度上升導致蒸散發(Evapotranspiration,ET)加劇,土壤水分嚴重下降并直接影響植被生長[1,7],同時,過度的蒸散發也會影響大氣循環、造成溫室效應,甚至出現極端氣候[7]。因此,準確識別青藏高原的蒸散發與植被變化的關系對陸地生態系統的穩定有十分重要的意義。
植被作為蒸散發的重要媒介,對氣候變化的響應明顯。近年來,青藏高原的植被呈現出退化趨勢,尤其是高原腹地植被退化趨勢明顯[8],因此,識別青藏高原植被退化梯度對蒸散發的研究有重要意義。當前國內外對青藏高原蒸散發的研究已取得較多成果,可以概括為以下幾個方面:一是對蒸散發模型的研發與應用[9-10],如張文旭等[11]分析了渦度相關法、蒸滲儀測定和FAO56Penman-Monteith公式3種方法計算青藏高原風火山的蒸散發,結果表明3種方法計算的蒸散發結果相似;二是通過站點實測數據或遙感數據探究青藏高原蒸散發的時空分布格局[12-14],如姚天次等[15]利用274個氣象站數據探究了1970—2017 年高原及周邊地區潛在蒸散發的空間格局,得出年潛在蒸散發呈現南北高、中部低的空間分布的特征;三是刻畫氣候變化的不同因子對蒸散發的影響研究[16-18],如Ma等[19]探究了青藏高原1982—2016 年的蒸散發及影響因素,研究結果表明降水增強是蒸散發加劇的主要驅動因素。這些研究均從計算蒸散發的方法和分析蒸散發的時間格局及其影響因子出發,較少關系植被退化與蒸散發的關系。因此,本文從植被退化出發,了解高寒地區植被退化與蒸散發的交互過程,這對當地生態安全和高寒植被退化的恢復治理有重要意義。
基于此,本研究選定青藏高原這一高寒干旱區,利用MODIS16A2數據,從時空尺度、空間尺度以及不同植被類型上分析蒸散發的變化特征,同時根據像元二分模型計算植被覆蓋度,識別青藏高原植被退化梯度和植被變化區,進而探究青藏高原不同退化梯度和不同植被變化區蒸散發的時空變化特征。
青藏高原是地球上最高的高原,地理范圍為68°~104° E,26°~40° N,總面積約250×104km2。該地區的地勢西高東低,年均氣溫-6℃~20℃,年降水量50~2 000 mm,是典型的高原高寒氣候,地域差異明顯,氣溫低,冬季長,日照充足,日較差大。從東南部向西北部,氣候由溫暖濕潤向寒冷干燥變化,水平帶上也先后出現了森林、灌木、草地等植被類型[20]。本文綜合考慮我國流域分區、省級行政區劃和“一帶一路”亞洲關鍵區域流域邊界[21],將青藏高原劃分為14 個子區域,分別是黃河流域、印度河流域、石羊河流域、疏勒河流域、黑河流域、塔里木河流域、阿姆河流域、雅魯藏布江流域、青海湖水系、羌塘高原區、長江流域、怒江流域、瀾滄江流域和柴達木盆地,以期從整體和局部的不同視角對青藏高原蒸散發的時空演變格局進行深入研究(圖1)。

圖1 青藏高原地理位置示意圖
1.2.1MODIS蒸散發數據 MODIS16A2是NASA發布的全球陸地蒸散發產品數據,基于Penman-Monteith公式計算得到,時間分辨率16天,空間分辨率500 m。通過全球通量塔臺站的檢驗,MOD16A2數據擬精度可達86%[22]。MOD16A2數據在時空上具有高分辨率、連續覆蓋等特點,在全球得到了較為廣泛的應用[23]。目前,MOD16A2數據已用于我國地表蒸散發的時空變化特征分析[24],總體上具有適用性,在站點資料稀缺的青藏高原區,MOD16A2數據能夠較準確地呈現蒸散發的空間分布規律[25]。本研究選用2001—2020 年的MOD16A2產品數據(來源:https://search.earthdata.nasa.gov/),通過Arc GIS10.2軟件的像元統計工具計算年際蒸散發數據。
1.2.2植被數據 本研究使用的中國植被圖比例尺為1:1 000 000,由中國科學院資源環境科學與數據中心(https://www.resdc.cn/)提供。該數據將其劃分為13 大類,可歸納為6類,分別是森林(針葉林、闊葉林、高山植被)、灌木(灌叢)、草地(草原、草甸、草叢)、荒漠、栽培植被和其他(沼澤和其他),本文選取了森林、灌木和草地進行研究。
1.2.3植被覆蓋度 植被覆蓋度是基于像元二分模型進行計算。像元二分模型是一種簡單實用的遙感估算模型,它假設一個像元的地表由有植被覆蓋部分地表與無植被覆蓋部分地表組成,而遙感傳感器觀測到的光譜信息也由這2個組分因子線性加權合成,各因子的權重是各自的面積在像元中所占的比率,如植被覆蓋度可以看作是植被的權重。
(1)
式中,NDVIsoil作為裸土的NDVI值,理論上應該趨近于0,但是由于遙感影像受到大氣環境和粗糙地表、土壤顏色、土壤地表濕度等環境的影響,NDVIsoil值不是固定值,會在一定范圍內變化,一般為-0.1~0.2。NDVIveg為研究區與像元最大NDVI值。鑒于此,對研究區內像元NDVI的灰度值進行統計分析,截取置信區間累計頻率在5%~95%對應的NDVI值分別作為NDVI最大值和最小值,從而計算得到植被覆蓋度FVC。
1.3.1青藏高原植被退化遙感監測和評價指標體系建立 根據國家標準GB19377-2003天然植被退化、沙化、鹽漬化的分級指標,以1982—1985 年每個像元的年際最大植被覆蓋度(AVHRR數據計算得到)作為基準,將其退化程度分為未退化、輕度退化、中度退化、重度退化和極重度退化5個級別[26-27](表1)。

表1 植被退化等級
1.3.2趨勢分析法 本文采用Theil-Sen中位數趨勢分析和Mann-Kendall檢驗來研究青藏高原蒸散發的空間趨勢特征。Theil-Sen中位數趨勢分析是一種穩健的非參數統計的趨勢計算方法,對測量誤差和離群數據不敏感,可用于ET的時間序列趨勢分析。而Mann-Kendall趨勢檢驗作為非參數統計檢驗的方法,用于評估Theil-Sen斜率估計的顯著性,即檢驗青藏高原ET趨勢的顯著性[28]。具體公式如下:
式中:β為斜率;i和j代表年份;如果β>0,ET數據集時間序列具有正趨勢;如果β<0,ET數據集序列具有負趨勢。
(3)
(4)
(5)
式中:xj和xi為ET時間序列數據;n為ET序列中數據個數;sgn()為符號函數;如果|Zc|≥Z1-α/2,則原假設被拒絕,ET序列呈顯著變化趨勢;如果Zc>0,則表明ET序列呈上升趨勢;如果Zc<0,則表明ET序列呈下降趨勢;α為給定的置信水平(顯著性檢驗水平),|Zc|≥1.96時表示ET序列通過了置信度為95%的顯著性檢驗。
1.3.3變異系數 變異系數用來表示地理數據的波動程度[29],在一定程度上可表明陸地蒸散量的變化情況。其計算公式如下:
(6)

1.3.4Hurst指數 基于重標極差(R/S)分析方法的Hurst指數(H值)是非線性時間序列分析的一種方法。常用于分析長時間序列的分形特征和長期記憶過程,也稱為重新標度極差分析[30]。H值從0到1。當H=0.5表明時間序列為隨機序列,前后變化無關聯;當H>0.5時,未來趨勢與前期保持一致(即繼續上升或下降),H越高,持續性越強。當H<0.5時,未來趨勢很可能是下一個時期的反持續性(即從上升到下降,反之亦然),H越低,反持續性越大[3]。
2.1.1青藏高原ET的年際變化特征 如圖2A所示,2001—2020 年青藏高原年蒸散發呈波動上升趨勢且始終保持在600~800 mm之間,平均值為730.48 mm。20年間,2001年研究區的ET最低,為643.58 mm,最高值出現在2013年,為780.49 mm,其次是2020年的779.83 mm。
青藏高原的植被類型包括森林、灌木和草地,20年里各植被類型的年均ET呈波動上升的趨勢,ET大小表現為:灌木>森林>草地(圖2B)。2006 年和2013 年,森林ET高于灌木ET,且在2012—2013年3種植被類型的ET均有大幅上升。由圖2C可以看出,2001—2020年青藏高原包含的14 個子流域的ET均呈波動上升的趨勢,其中黃河流域區ET最高,各年份均在800 mm以上,長江流域次之,均在750 mm以上。羌塘高原區ET最低,年均ET在620 mm以下。四大內陸河流域中石羊河流域ET>黑河流域ET>疏勒河流域ET>塔里木河流域ET,其中石羊河流域的ET在20年里的變化量最大,為376.81 mm。

圖2 青藏高原蒸散發(A)、不同植被類型蒸散發(B)以及子流域蒸散發變化圖(C)
2.1.2青藏高原ET的空間變化特征 青藏高原年平均蒸散發在空間上呈中部低,東南高的態勢(圖3A)。據統計,年均ET范圍為800~1 000 mm的面積最大,占總面積的31.14%,其次為1 000~1 200 mm,占總面積的30.74%,0~200 mm面積最小,占總面積的0.03%。從空間上來看,位于西藏自治區的拉薩市、日喀則地區、那曲地區以及青海省的西寧市的ET較低,這可能與拉薩、西寧作為省會城市,城市化的發展速度較快導致植被覆蓋度降低有關。
2001—2020年青藏高原ET的變化趨勢如圖3B所示,斜率介于62.9 至54.3 之間,80.02%的研究區域的ET呈增加趨勢。青藏高原東部地區的ET變化相對強于西部地區,中部和南部的ET顯著下降,如西藏自治區的林芝地區,這可能與氣候變化和高山高原植被退化有關。2001—2020年青藏高原ET的變異系數在0.02~1.125之間(圖3C),表明ET的變化具有顯著的空間異質性。青藏高原大部分地區ET波動較低,是由于該區域人口相對較少,高山植被覆蓋穩定,變異系數高于0.5的區域集中在區域南部,占全域面積的0.29%。青藏高原蒸散發的Hurst指數在0.06~0.904之間(圖3D),平均值為0.408。Hurst指數<0.5的區域占總面積的86.52%,這表明ET在青藏高原的未來趨勢中具有顯著的反持續性,其中那曲地區的東南部、海西蒙古族自治區和玉樹藏族自治區的交接處ET的未來變化趨勢與當前變化趨勢相反,即未來將呈增加趨勢。
青藏高原各植被類型存在顯著的空間異質性,草地面積最大且廣泛分布于研究區的中部和西部,而森林和灌木呈片狀或帶狀集中分布在青藏高原的東南部(圖3E)。2001—2020年,森林、灌木和草地年均ET分別為802.00 mm,832.21 mm,705.68 mm。從圖3F可以看出,青藏高原子流域的年均ET在空間上呈“西北低,東南高”的態勢,這與本研究區的海拔呈顯著負相關關系,即海拔高的地區ET較低。青藏高原14個子流域中羌塘高原區的年均ET最低,為515.15 mm,這可能與羌塘高原區位于青藏高原腹地,海拔高且四周高山重圍、氣候寒冷干燥,植被難以生長有關。黃河流域的年均ET最高,為909.69 mm,其次為長江流域,為834.13 mm。
2.2.1青藏高原植被退化時空變化格局 青藏高原植被退化自西北向東南呈先減輕后加劇的態勢(圖4A),極重度退化的區域集中在青藏高原西北部,如柴達木盆地和昆侖山脈周邊。植被未退化及輕度退化區域普遍位于黃河流域上游、青海湖水系以及長江流域、瀾滄江流域的部分地區。同時,青藏高原東南部的森林地區也出現一定程度的退化。
從時間上來看(圖4B),2001—2020年青藏高原各植被退化區面積一直保持著:極重度退化區>重度退化區>未退化區>輕度退化區>重度退化區的態勢,其中輕度退化和重度退化在20年里的變化不大,且面積基本持平。20年間,極重度退化區面積呈下降趨勢,減少了9.88×104km2,重度退化區面積呈波動上升趨勢,增加了3.05×104km2,未退化區面積總體呈大幅度下降繼而顯著上升趨勢,共增加了10.91×104km2,其中降幅最大的時間在2013年,降低了11.16×104km2,增幅最大在2020年,增加了7.71×104km2。這說明20年間青藏高原植被退化程度高的區域在不斷縮小,未退化區域不斷增加且增幅在不斷增強。

圖4 青藏高原植被退化時空變化圖
2.2.2青藏高原植被退化轉移矩陣 表2可知,青藏高原2001年未退化、輕度退化、中度退化、重度退化、極重度退化面積占比分別為:23.46%,10.37%,9.97%,27.33%,28.88%;2020年面積占比分別為:28.32%,10.60%,10.26%,27.41%,23.41%。2001年極重度退化區域面積最大,2020年則為未退化區域面積最大,表明青藏高原的植被退化在不斷減弱。

表2 青藏高原植被退化面積占比轉移矩陣
2.3.1青藏高原不同退化梯度內ET變化特征 青藏高原2001—2020年各退化梯度下的ET總體呈波動上升趨勢,其中各退化區ET的均值表現為:未退化區>輕度退化區>極重度退化區>中度退化區>重度退化區,分別是869.89 mm,738.77 mm,731.14 mm,684.40 mm,610.83 mm,這說明植被退化在最低程度或最高程度時蒸散發較高(圖5)。
五種退化梯度中未退化區ET的增幅最大,年均斜率為6.56,其次為極重度退化、輕度退化、中度退化和重度退化,其斜率分別為:4.74,4.37,3.84,2.92,同時各退化梯度ET均在2012年大幅下降后在2013年顯著上升。
2.3.2青藏高原植被變化區內ET的時空變化特征 圖6A可以看出,青藏高原植被不變區、恢復區以及惡化區面積分別是:126.77×104km2,52.27×104km2,26.32×104km2。在空間上,植被不變區廣泛分布于青藏高原全域,占全域面積的61.7%,植被恢復區集中在青海湖周邊、柴達木盆地外圍以及黃河流域的西北部地區,植被惡化區則呈點狀或片狀分布于青藏高原的西部和南部。識別出青藏高原20年間植被變化區對該地區生態修復政策制定和推廣具有一定意義。
2001—2020年青藏高原各退化區的ET整體呈波動上升趨勢(圖6B),且植被恢復區ET>植被不變區ET>植被惡化區ET,3個退化區均在2012 年大幅下跌后在2013 年急劇上升。20年間,植被恢復區的ET增量最多,共增加了184.43 mm,植被不變區和植被惡化區的ET分別增加了118.69 mm和118.03 mm。20年間,青藏高原植被恢復相較于植被惡化更能促使地區蒸散發的增加。

圖6 青藏高原植被變化區的蒸散發時空變化圖
在全球氣候變暖和綠化趨勢增強的背景下,蒸散發的增強是主要趨勢[31-32]。Pascolini-Campbell等[1]證明全球蒸散發呈顯著的正線性增長趨勢,增長的部分來源于降水量的增加和徑流的減少,Zeng等[33]證明1982—2011 年,全球蒸散發呈現出顯著的增長趨勢,其中超過一半的增長歸因于地球的綠化。氣象因素對蒸散發也具有顯著的驅動作用,如降水量被認為是影響中國ET變化的主導因素,風速、溫度、日照時數和相對濕度等均對蒸散發有影響[7]。本研究表明,2001—2020 年青藏高原蒸散發呈波動上升趨勢蒸散發的增強會導致干旱半干旱地區更多的水分流失[1],從而抑制植物生長,破壞生態系統水分循環,最終可能影響整個陸地生態系統的穩定性。此外,青藏高原年均蒸散發在全域尺度、子流域尺度以及不同植被類型中均表現為2012 年大幅下降后2013 年大幅上升,這種現象可能與該年份出現比正常年份更高的溫度、降水以及更低的相對濕度等,也可能與我國生態恢復政策有關[34-35]。裴婷婷等[36]證明2012 年以后全國的降水明顯升高,同時降水引起的蒸散發也大幅度增加,這在一定程度上可以解釋2012 年以后蒸散發的顯著升高。青藏高原草地面積占比最大,其次是森林和灌木,其中灌木的年均ET和森林的年均ET相差不多,但略大于森林。張永強等[37]在探究2003—2017年植被變化對全球陸面蒸散發的影響時,得出2012 年以后植被變化對蒸散發的影響顯著加強,灌木和耕地對蒸散發的影響最為顯著,這與本文的研究結果較為一致。青藏高原14 個子流域中陸地蒸散量均處于波動上升的趨勢,這與Li[38]對子流域的ET不斷升高的研究結果相同。因此,本文的研究結果對認識青藏高原不同空間維度的蒸散發格局提供了參考。
2001—2020 年,青藏高原植被退化在不斷減弱,這與部分學者的研究結果相似。如夏龍等[39]在對青藏高原草地退化的研究中證明青藏高原草地整體上在朝著改善的方向發展。Wang等[40]在探究青藏高原植被變化的主導因素中得出,青藏高原植被在2000—2015 年有改善趨勢。同時,全球植被綠化程度的不斷增強也能從側面證明植被退化在一定程度上減輕[2]。本文的研究結果表明蒸散發在未退化梯度上最高,其次為極重度退化梯度,這可能是由于蒸散發的途徑包括植被蒸騰、植被冠層截留降水以及露水蒸發,從而使得植被覆蓋度高的地區蒸散發較高。極重度退化區的蒸散發也較高可能與降水和地表徑流直接被劃分為蒸散發有關。
青藏高原的青海湖和柴達木盆地周邊地區的植被在20 年里呈逐步恢復態勢且恢復程度顯著,這可能與我國出臺的生態修復政策有關,如針對青海湖和柴達木周邊地區的生態問題,集中采用退耕還草、輪牧、灌草結合、圍欄封育等措施來促進植被恢復[41]。柴達木盆地周邊地區的主要植被類型為草地,1980—2018 年草地面積呈增加趨勢,其中高中覆蓋度草地面積增加趨勢明顯,說明柴達木盆地植被退化程度在逐漸減輕[42-43]。劉秋漫[44]對柴達木盆地的研究表明,近30 年間該地區生態環境質量整體上呈現出先下降后上升的趨勢,生態環境質量總體上呈現出向好的方向發展的態勢,這與本文的研究結果相似。本研究仍存在一些不足,首先在產品選擇上,MODIS產品數據集的反演模型在全球通用,但在研究區的適用性有待提高,同時產品數據的精度相對不高,研究結果存在不確定性,其次本文采用中國科學院資源環境科學與數據中心的植被分類體系將植被劃分為林地、草地和灌木,該分類體系會對研究結果產生影響。
本文基于2001—2020年MODIS16A2數據,探究了青藏高原蒸散發的時空變化格局同時識別了青藏高原植被退化梯度,進而探究五種退化梯度下青藏高原蒸散發的變化特征。
2001—2020年青藏高原年均ET呈波動上升的趨勢,其中2001年ET最低,為643.58 mm,2013年最高,為780.49 mm。從空間上來看年均ET呈中部低,東南高的態勢,低值集中在羌塘高原區且86.52%的區域未來趨勢具有顯著的反持續性,即青藏高原大部分區域的蒸散發將呈現下降趨勢。
青藏高原各自然植被類型ET中:灌木>森林>草地,且在2012年大幅下降繼而在2013年大幅上升。青藏高原范圍內的子流域ET均呈波動上升趨勢,其中黃河流域的ET最大,羌塘高原區的ET最小,四大內陸河流域中石羊河流域ET>黑河流域ET>疏勒河流域ET>塔里木河流域ET。
2001—2020年青藏高原植被退化不斷減弱,在空間上呈自西向東呈退化先減輕后逐步加強的趨勢,各退化區面積一直保持著:極重度退化區>重度退化區>未退化區>輕度退化區>重度退化區。青藏高原各植被變化區的面積為:不變區>恢復區>惡化區,在空間上表現為不變區廣泛分布于青藏高原內部,恢復區集中分布在青藏高原東北部,如青海湖和環柴達木盆地周邊區域,惡化區呈點片狀分布于青藏高原西南部。20年里各植被變化區的年均ET呈波動上升態勢,表現為:恢復區>不變區>惡化區。