劉長雨,謝保鵬*,楊 潔,陳 英,裴婷婷
(1.甘肅農(nóng)業(yè)大學管理學院,甘肅 蘭州 730070;2.甘肅農(nóng)業(yè)大學草業(yè)學院,甘肅 蘭州 730070)
全球氣候變暖和人類活動的共同影響下,地球正以一定的速度綠化[1-2]。美國航天局衛(wèi)星數(shù)據(jù)表明,全球從2000—2017 年新增的綠化面積中,約1/4來自中國,貢獻比例居全球首位[3]。青藏高原作為我國乃至亞洲的重要生態(tài)安全屏障[4],其獨特的地理位置和典型的高寒生態(tài)系統(tǒng)使其成為國內外學者研究的熱點區(qū)域。過去的50年里,青藏高原氣候變化的暖濕化特征明顯[5],溫度每10 年升高0.3℃~0.4℃,降水每10 年增加2.2%[6]。降水增加和溫度上升導致蒸散發(fā)(Evapotranspiration,ET)加劇,土壤水分嚴重下降并直接影響植被生長[1,7],同時,過度的蒸散發(fā)也會影響大氣循環(huán)、造成溫室效應,甚至出現(xiàn)極端氣候[7]。因此,準確識別青藏高原的蒸散發(fā)與植被變化的關系對陸地生態(tài)系統(tǒng)的穩(wěn)定有十分重要的意義。
植被作為蒸散發(fā)的重要媒介,對氣候變化的響應明顯。近年來,青藏高原的植被呈現(xiàn)出退化趨勢,尤其是高原腹地植被退化趨勢明顯[8],因此,識別青藏高原植被退化梯度對蒸散發(fā)的研究有重要意義。當前國內外對青藏高原蒸散發(fā)的研究已取得較多成果,可以概括為以下幾個方面:一是對蒸散發(fā)模型的研發(fā)與應用[9-10],如張文旭等[11]分析了渦度相關法、蒸滲儀測定和FAO56Penman-Monteith公式3種方法計算青藏高原風火山的蒸散發(fā),結果表明3種方法計算的蒸散發(fā)結果相似;二是通過站點實測數(shù)據(jù)或遙感數(shù)據(jù)探究青藏高原蒸散發(fā)的時空分布格局[12-14],如姚天次等[15]利用274個氣象站數(shù)據(jù)探究了1970—2017 年高原及周邊地區(qū)潛在蒸散發(fā)的空間格局,得出年潛在蒸散發(fā)呈現(xiàn)南北高、中部低的空間分布的特征;三是刻畫氣候變化的不同因子對蒸散發(fā)的影響研究[16-18],如Ma等[19]探究了青藏高原1982—2016 年的蒸散發(fā)及影響因素,研究結果表明降水增強是蒸散發(fā)加劇的主要驅動因素。這些研究均從計算蒸散發(fā)的方法和分析蒸散發(fā)的時間格局及其影響因子出發(fā),較少關系植被退化與蒸散發(fā)的關系。因此,本文從植被退化出發(fā),了解高寒地區(qū)植被退化與蒸散發(fā)的交互過程,這對當?shù)厣鷳B(tài)安全和高寒植被退化的恢復治理有重要意義。
基于此,本研究選定青藏高原這一高寒干旱區(qū),利用MODIS16A2數(shù)據(jù),從時空尺度、空間尺度以及不同植被類型上分析蒸散發(fā)的變化特征,同時根據(jù)像元二分模型計算植被覆蓋度,識別青藏高原植被退化梯度和植被變化區(qū),進而探究青藏高原不同退化梯度和不同植被變化區(qū)蒸散發(fā)的時空變化特征。
青藏高原是地球上最高的高原,地理范圍為68°~104° E,26°~40° N,總面積約250×104km2。該地區(qū)的地勢西高東低,年均氣溫-6℃~20℃,年降水量50~2 000 mm,是典型的高原高寒氣候,地域差異明顯,氣溫低,冬季長,日照充足,日較差大。從東南部向西北部,氣候由溫暖濕潤向寒冷干燥變化,水平帶上也先后出現(xiàn)了森林、灌木、草地等植被類型[20]。本文綜合考慮我國流域分區(qū)、省級行政區(qū)劃和“一帶一路”亞洲關鍵區(qū)域流域邊界[21],將青藏高原劃分為14 個子區(qū)域,分別是黃河流域、印度河流域、石羊河流域、疏勒河流域、黑河流域、塔里木河流域、阿姆河流域、雅魯藏布江流域、青海湖水系、羌塘高原區(qū)、長江流域、怒江流域、瀾滄江流域和柴達木盆地,以期從整體和局部的不同視角對青藏高原蒸散發(fā)的時空演變格局進行深入研究(圖1)。
圖1 青藏高原地理位置示意圖
1.2.1MODIS蒸散發(fā)數(shù)據(jù) MODIS16A2是NASA發(fā)布的全球陸地蒸散發(fā)產(chǎn)品數(shù)據(jù),基于Penman-Monteith公式計算得到,時間分辨率16天,空間分辨率500 m。通過全球通量塔臺站的檢驗,MOD16A2數(shù)據(jù)擬精度可達86%[22]。MOD16A2數(shù)據(jù)在時空上具有高分辨率、連續(xù)覆蓋等特點,在全球得到了較為廣泛的應用[23]。目前,MOD16A2數(shù)據(jù)已用于我國地表蒸散發(fā)的時空變化特征分析[24],總體上具有適用性,在站點資料稀缺的青藏高原區(qū),MOD16A2數(shù)據(jù)能夠較準確地呈現(xiàn)蒸散發(fā)的空間分布規(guī)律[25]。本研究選用2001—2020 年的MOD16A2產(chǎn)品數(shù)據(jù)(來源:https://search.earthdata.nasa.gov/),通過Arc GIS10.2軟件的像元統(tǒng)計工具計算年際蒸散發(fā)數(shù)據(jù)。
1.2.2植被數(shù)據(jù) 本研究使用的中國植被圖比例尺為1:1 000 000,由中國科學院資源環(huán)境科學與數(shù)據(jù)中心(https://www.resdc.cn/)提供。該數(shù)據(jù)將其劃分為13 大類,可歸納為6類,分別是森林(針葉林、闊葉林、高山植被)、灌木(灌叢)、草地(草原、草甸、草叢)、荒漠、栽培植被和其他(沼澤和其他),本文選取了森林、灌木和草地進行研究。
1.2.3植被覆蓋度 植被覆蓋度是基于像元二分模型進行計算。像元二分模型是一種簡單實用的遙感估算模型,它假設一個像元的地表由有植被覆蓋部分地表與無植被覆蓋部分地表組成,而遙感傳感器觀測到的光譜信息也由這2個組分因子線性加權合成,各因子的權重是各自的面積在像元中所占的比率,如植被覆蓋度可以看作是植被的權重。
(1)
式中,NDVIsoil作為裸土的NDVI值,理論上應該趨近于0,但是由于遙感影像受到大氣環(huán)境和粗糙地表、土壤顏色、土壤地表濕度等環(huán)境的影響,NDVIsoil值不是固定值,會在一定范圍內變化,一般為-0.1~0.2。NDVIveg為研究區(qū)與像元最大NDVI值。鑒于此,對研究區(qū)內像元NDVI的灰度值進行統(tǒng)計分析,截取置信區(qū)間累計頻率在5%~95%對應的NDVI值分別作為NDVI最大值和最小值,從而計算得到植被覆蓋度FVC。
1.3.1青藏高原植被退化遙感監(jiān)測和評價指標體系建立 根據(jù)國家標準GB19377-2003天然植被退化、沙化、鹽漬化的分級指標,以1982—1985 年每個像元的年際最大植被覆蓋度(AVHRR數(shù)據(jù)計算得到)作為基準,將其退化程度分為未退化、輕度退化、中度退化、重度退化和極重度退化5個級別[26-27](表1)。
表1 植被退化等級
1.3.2趨勢分析法 本文采用Theil-Sen中位數(shù)趨勢分析和Mann-Kendall檢驗來研究青藏高原蒸散發(fā)的空間趨勢特征。Theil-Sen中位數(shù)趨勢分析是一種穩(wěn)健的非參數(shù)統(tǒng)計的趨勢計算方法,對測量誤差和離群數(shù)據(jù)不敏感,可用于ET的時間序列趨勢分析。而Mann-Kendall趨勢檢驗作為非參數(shù)統(tǒng)計檢驗的方法,用于評估Theil-Sen斜率估計的顯著性,即檢驗青藏高原ET趨勢的顯著性[28]。具體公式如下:
式中:β為斜率;i和j代表年份;如果β>0,ET數(shù)據(jù)集時間序列具有正趨勢;如果β<0,ET數(shù)據(jù)集序列具有負趨勢。
(3)
(4)
(5)
式中:xj和xi為ET時間序列數(shù)據(jù);n為ET序列中數(shù)據(jù)個數(shù);sgn()為符號函數(shù);如果|Zc|≥Z1-α/2,則原假設被拒絕,ET序列呈顯著變化趨勢;如果Zc>0,則表明ET序列呈上升趨勢;如果Zc<0,則表明ET序列呈下降趨勢;α為給定的置信水平(顯著性檢驗水平),|Zc|≥1.96時表示ET序列通過了置信度為95%的顯著性檢驗。
1.3.3變異系數(shù) 變異系數(shù)用來表示地理數(shù)據(jù)的波動程度[29],在一定程度上可表明陸地蒸散量的變化情況。其計算公式如下:
(6)
1.3.4Hurst指數(shù) 基于重標極差(R/S)分析方法的Hurst指數(shù)(H值)是非線性時間序列分析的一種方法。常用于分析長時間序列的分形特征和長期記憶過程,也稱為重新標度極差分析[30]。H值從0到1。當H=0.5表明時間序列為隨機序列,前后變化無關聯(lián);當H>0.5時,未來趨勢與前期保持一致(即繼續(xù)上升或下降),H越高,持續(xù)性越強。當H<0.5時,未來趨勢很可能是下一個時期的反持續(xù)性(即從上升到下降,反之亦然),H越低,反持續(xù)性越大[3]。
2.1.1青藏高原ET的年際變化特征 如圖2A所示,2001—2020 年青藏高原年蒸散發(fā)呈波動上升趨勢且始終保持在600~800 mm之間,平均值為730.48 mm。20年間,2001年研究區(qū)的ET最低,為643.58 mm,最高值出現(xiàn)在2013年,為780.49 mm,其次是2020年的779.83 mm。
青藏高原的植被類型包括森林、灌木和草地,20年里各植被類型的年均ET呈波動上升的趨勢,ET大小表現(xiàn)為:灌木>森林>草地(圖2B)。2006 年和2013 年,森林ET高于灌木ET,且在2012—2013年3種植被類型的ET均有大幅上升。由圖2C可以看出,2001—2020年青藏高原包含的14 個子流域的ET均呈波動上升的趨勢,其中黃河流域區(qū)ET最高,各年份均在800 mm以上,長江流域次之,均在750 mm以上。羌塘高原區(qū)ET最低,年均ET在620 mm以下。四大內陸河流域中石羊河流域ET>黑河流域ET>疏勒河流域ET>塔里木河流域ET,其中石羊河流域的ET在20年里的變化量最大,為376.81 mm。
圖2 青藏高原蒸散發(fā)(A)、不同植被類型蒸散發(fā)(B)以及子流域蒸散發(fā)變化圖(C)
2.1.2青藏高原ET的空間變化特征 青藏高原年平均蒸散發(fā)在空間上呈中部低,東南高的態(tài)勢(圖3A)。據(jù)統(tǒng)計,年均ET范圍為800~1 000 mm的面積最大,占總面積的31.14%,其次為1 000~1 200 mm,占總面積的30.74%,0~200 mm面積最小,占總面積的0.03%。從空間上來看,位于西藏自治區(qū)的拉薩市、日喀則地區(qū)、那曲地區(qū)以及青海省的西寧市的ET較低,這可能與拉薩、西寧作為省會城市,城市化的發(fā)展速度較快導致植被覆蓋度降低有關。
2001—2020年青藏高原ET的變化趨勢如圖3B所示,斜率介于62.9 至54.3 之間,80.02%的研究區(qū)域的ET呈增加趨勢。青藏高原東部地區(qū)的ET變化相對強于西部地區(qū),中部和南部的ET顯著下降,如西藏自治區(qū)的林芝地區(qū),這可能與氣候變化和高山高原植被退化有關。2001—2020年青藏高原ET的變異系數(shù)在0.02~1.125之間(圖3C),表明ET的變化具有顯著的空間異質性。青藏高原大部分地區(qū)ET波動較低,是由于該區(qū)域人口相對較少,高山植被覆蓋穩(wěn)定,變異系數(shù)高于0.5的區(qū)域集中在區(qū)域南部,占全域面積的0.29%。青藏高原蒸散發(fā)的Hurst指數(shù)在0.06~0.904之間(圖3D),平均值為0.408。Hurst指數(shù)<0.5的區(qū)域占總面積的86.52%,這表明ET在青藏高原的未來趨勢中具有顯著的反持續(xù)性,其中那曲地區(qū)的東南部、海西蒙古族自治區(qū)和玉樹藏族自治區(qū)的交接處ET的未來變化趨勢與當前變化趨勢相反,即未來將呈增加趨勢。
青藏高原各植被類型存在顯著的空間異質性,草地面積最大且廣泛分布于研究區(qū)的中部和西部,而森林和灌木呈片狀或帶狀集中分布在青藏高原的東南部(圖3E)。2001—2020年,森林、灌木和草地年均ET分別為802.00 mm,832.21 mm,705.68 mm。從圖3F可以看出,青藏高原子流域的年均ET在空間上呈“西北低,東南高”的態(tài)勢,這與本研究區(qū)的海拔呈顯著負相關關系,即海拔高的地區(qū)ET較低。青藏高原14個子流域中羌塘高原區(qū)的年均ET最低,為515.15 mm,這可能與羌塘高原區(qū)位于青藏高原腹地,海拔高且四周高山重圍、氣候寒冷干燥,植被難以生長有關。黃河流域的年均ET最高,為909.69 mm,其次為長江流域,為834.13 mm。
2.2.1青藏高原植被退化時空變化格局 青藏高原植被退化自西北向東南呈先減輕后加劇的態(tài)勢(圖4A),極重度退化的區(qū)域集中在青藏高原西北部,如柴達木盆地和昆侖山脈周邊。植被未退化及輕度退化區(qū)域普遍位于黃河流域上游、青海湖水系以及長江流域、瀾滄江流域的部分地區(qū)。同時,青藏高原東南部的森林地區(qū)也出現(xiàn)一定程度的退化。
從時間上來看(圖4B),2001—2020年青藏高原各植被退化區(qū)面積一直保持著:極重度退化區(qū)>重度退化區(qū)>未退化區(qū)>輕度退化區(qū)>重度退化區(qū)的態(tài)勢,其中輕度退化和重度退化在20年里的變化不大,且面積基本持平。20年間,極重度退化區(qū)面積呈下降趨勢,減少了9.88×104km2,重度退化區(qū)面積呈波動上升趨勢,增加了3.05×104km2,未退化區(qū)面積總體呈大幅度下降繼而顯著上升趨勢,共增加了10.91×104km2,其中降幅最大的時間在2013年,降低了11.16×104km2,增幅最大在2020年,增加了7.71×104km2。這說明20年間青藏高原植被退化程度高的區(qū)域在不斷縮小,未退化區(qū)域不斷增加且增幅在不斷增強。
圖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年極重度退化區(qū)域面積最大,2020年則為未退化區(qū)域面積最大,表明青藏高原的植被退化在不斷減弱。
表2 青藏高原植被退化面積占比轉移矩陣
2.3.1青藏高原不同退化梯度內ET變化特征 青藏高原2001—2020年各退化梯度下的ET總體呈波動上升趨勢,其中各退化區(qū)ET的均值表現(xiàn)為:未退化區(qū)>輕度退化區(qū)>極重度退化區(qū)>中度退化區(qū)>重度退化區(qū),分別是869.89 mm,738.77 mm,731.14 mm,684.40 mm,610.83 mm,這說明植被退化在最低程度或最高程度時蒸散發(fā)較高(圖5)。
五種退化梯度中未退化區(qū)ET的增幅最大,年均斜率為6.56,其次為極重度退化、輕度退化、中度退化和重度退化,其斜率分別為:4.74,4.37,3.84,2.92,同時各退化梯度ET均在2012年大幅下降后在2013年顯著上升。
2.3.2青藏高原植被變化區(qū)內ET的時空變化特征 圖6A可以看出,青藏高原植被不變區(qū)、恢復區(qū)以及惡化區(qū)面積分別是:126.77×104km2,52.27×104km2,26.32×104km2。在空間上,植被不變區(qū)廣泛分布于青藏高原全域,占全域面積的61.7%,植被恢復區(qū)集中在青海湖周邊、柴達木盆地外圍以及黃河流域的西北部地區(qū),植被惡化區(qū)則呈點狀或片狀分布于青藏高原的西部和南部。識別出青藏高原20年間植被變化區(qū)對該地區(qū)生態(tài)修復政策制定和推廣具有一定意義。
2001—2020年青藏高原各退化區(qū)的ET整體呈波動上升趨勢(圖6B),且植被恢復區(qū)ET>植被不變區(qū)ET>植被惡化區(qū)ET,3個退化區(qū)均在2012 年大幅下跌后在2013 年急劇上升。20年間,植被恢復區(qū)的ET增量最多,共增加了184.43 mm,植被不變區(qū)和植被惡化區(qū)的ET分別增加了118.69 mm和118.03 mm。20年間,青藏高原植被恢復相較于植被惡化更能促使地區(qū)蒸散發(fā)的增加。
圖6 青藏高原植被變化區(qū)的蒸散發(fā)時空變化圖
在全球氣候變暖和綠化趨勢增強的背景下,蒸散發(fā)的增強是主要趨勢[31-32]。Pascolini-Campbell等[1]證明全球蒸散發(fā)呈顯著的正線性增長趨勢,增長的部分來源于降水量的增加和徑流的減少,Zeng等[33]證明1982—2011 年,全球蒸散發(fā)呈現(xiàn)出顯著的增長趨勢,其中超過一半的增長歸因于地球的綠化。氣象因素對蒸散發(fā)也具有顯著的驅動作用,如降水量被認為是影響中國ET變化的主導因素,風速、溫度、日照時數(shù)和相對濕度等均對蒸散發(fā)有影響[7]。本研究表明,2001—2020 年青藏高原蒸散發(fā)呈波動上升趨勢蒸散發(fā)的增強會導致干旱半干旱地區(qū)更多的水分流失[1],從而抑制植物生長,破壞生態(tài)系統(tǒng)水分循環(huán),最終可能影響整個陸地生態(tài)系統(tǒng)的穩(wěn)定性。此外,青藏高原年均蒸散發(fā)在全域尺度、子流域尺度以及不同植被類型中均表現(xiàn)為2012 年大幅下降后2013 年大幅上升,這種現(xiàn)象可能與該年份出現(xiàn)比正常年份更高的溫度、降水以及更低的相對濕度等,也可能與我國生態(tài)恢復政策有關[34-35]。裴婷婷等[36]證明2012 年以后全國的降水明顯升高,同時降水引起的蒸散發(fā)也大幅度增加,這在一定程度上可以解釋2012 年以后蒸散發(fā)的顯著升高。青藏高原草地面積占比最大,其次是森林和灌木,其中灌木的年均ET和森林的年均ET相差不多,但略大于森林。張永強等[37]在探究2003—2017年植被變化對全球陸面蒸散發(fā)的影響時,得出2012 年以后植被變化對蒸散發(fā)的影響顯著加強,灌木和耕地對蒸散發(fā)的影響最為顯著,這與本文的研究結果較為一致。青藏高原14 個子流域中陸地蒸散量均處于波動上升的趨勢,這與Li[38]對子流域的ET不斷升高的研究結果相同。因此,本文的研究結果對認識青藏高原不同空間維度的蒸散發(fā)格局提供了參考。
2001—2020 年,青藏高原植被退化在不斷減弱,這與部分學者的研究結果相似。如夏龍等[39]在對青藏高原草地退化的研究中證明青藏高原草地整體上在朝著改善的方向發(fā)展。Wang等[40]在探究青藏高原植被變化的主導因素中得出,青藏高原植被在2000—2015 年有改善趨勢。同時,全球植被綠化程度的不斷增強也能從側面證明植被退化在一定程度上減輕[2]。本文的研究結果表明蒸散發(fā)在未退化梯度上最高,其次為極重度退化梯度,這可能是由于蒸散發(fā)的途徑包括植被蒸騰、植被冠層截留降水以及露水蒸發(fā),從而使得植被覆蓋度高的地區(qū)蒸散發(fā)較高。極重度退化區(qū)的蒸散發(fā)也較高可能與降水和地表徑流直接被劃分為蒸散發(fā)有關。
青藏高原的青海湖和柴達木盆地周邊地區(qū)的植被在20 年里呈逐步恢復態(tài)勢且恢復程度顯著,這可能與我國出臺的生態(tài)修復政策有關,如針對青海湖和柴達木周邊地區(qū)的生態(tài)問題,集中采用退耕還草、輪牧、灌草結合、圍欄封育等措施來促進植被恢復[41]。柴達木盆地周邊地區(qū)的主要植被類型為草地,1980—2018 年草地面積呈增加趨勢,其中高中覆蓋度草地面積增加趨勢明顯,說明柴達木盆地植被退化程度在逐漸減輕[42-43]。劉秋漫[44]對柴達木盆地的研究表明,近30 年間該地區(qū)生態(tài)環(huán)境質量整體上呈現(xiàn)出先下降后上升的趨勢,生態(tài)環(huán)境質量總體上呈現(xiàn)出向好的方向發(fā)展的態(tài)勢,這與本文的研究結果相似。本研究仍存在一些不足,首先在產(chǎn)品選擇上,MODIS產(chǎn)品數(shù)據(jù)集的反演模型在全球通用,但在研究區(qū)的適用性有待提高,同時產(chǎn)品數(shù)據(jù)的精度相對不高,研究結果存在不確定性,其次本文采用中國科學院資源環(huán)境科學與數(shù)據(jù)中心的植被分類體系將植被劃分為林地、草地和灌木,該分類體系會對研究結果產(chǎn)生影響。
本文基于2001—2020年MODIS16A2數(shù)據(jù),探究了青藏高原蒸散發(fā)的時空變化格局同時識別了青藏高原植被退化梯度,進而探究五種退化梯度下青藏高原蒸散發(fā)的變化特征。
2001—2020年青藏高原年均ET呈波動上升的趨勢,其中2001年ET最低,為643.58 mm,2013年最高,為780.49 mm。從空間上來看年均ET呈中部低,東南高的態(tài)勢,低值集中在羌塘高原區(qū)且86.52%的區(qū)域未來趨勢具有顯著的反持續(xù)性,即青藏高原大部分區(qū)域的蒸散發(fā)將呈現(xiàn)下降趨勢。
青藏高原各自然植被類型ET中:灌木>森林>草地,且在2012年大幅下降繼而在2013年大幅上升。青藏高原范圍內的子流域ET均呈波動上升趨勢,其中黃河流域的ET最大,羌塘高原區(qū)的ET最小,四大內陸河流域中石羊河流域ET>黑河流域ET>疏勒河流域ET>塔里木河流域ET。
2001—2020年青藏高原植被退化不斷減弱,在空間上呈自西向東呈退化先減輕后逐步加強的趨勢,各退化區(qū)面積一直保持著:極重度退化區(qū)>重度退化區(qū)>未退化區(qū)>輕度退化區(qū)>重度退化區(qū)。青藏高原各植被變化區(qū)的面積為:不變區(qū)>恢復區(qū)>惡化區(qū),在空間上表現(xiàn)為不變區(qū)廣泛分布于青藏高原內部,恢復區(qū)集中分布在青藏高原東北部,如青海湖和環(huán)柴達木盆地周邊區(qū)域,惡化區(qū)呈點片狀分布于青藏高原西南部。20年里各植被變化區(qū)的年均ET呈波動上升態(tài)勢,表現(xiàn)為:恢復區(qū)>不變區(qū)>惡化區(qū)。