地表水熱要素在青藏高原草地退化中的作用

夏 龍,宋小寧,蔡碩豪,胡容海,郭 達

中國科學院大學資源與環境學院, 北京 100049

植被作為地表的重要覆蓋類型,是陸地生態系統的核心組成部分。植被不僅在生物圈、大氣圈和全球碳循環中扮演著重要的角色,而且也是連接土壤、水分和能量等自然要素的紐帶[1]。各種地表植被之間的差異可以反映出不同的水熱組合,所以監測植被之間的差異,在地方、區域和全球范圍內,對生態系統的狀態,環境壓力和景觀變化的研究具有重要指示作用[2]。草地植被是全球陸地生態系統中面積最大的可更新資源,一方面,草地植被對培育土壤肥力、防止水土流失、維持陸地生態系統平衡起著重要作用,另一方面草地植被也可以為食草動物提供飼料,為人類生產食物、藥品和工業原料等[3-4]。但是,草地生態系統生境脆弱,比較容易受到外界環境的干擾,因而對全球氣候環境變化十分敏感[5-6]。

青藏高原總面積約為250萬km2,約占中國陸地總面積的四分之一,平均海拔在4000 m以上,是世界上海拔最高的高原,常被稱為世界的“第三極”,高海拔、地處中低緯度的地理條件以及特殊的下墊面條件也使得青藏高原成為全球對氣候變化最為敏感的區域之一。近年來,在全球氣候變化和人為因素的共同作用下,青藏高原生態環境發生了一系列的負面變化,例如濕地和草地面積減小、凍土退化、土地沙漠化嚴重、水域減少、生物多樣性減小及自然災害增多等[7-10]。青藏高原地區天然草地面積約為1.5×108hm2,約占全國草地總面積的三分之一。青藏高原地區特殊的地理環境加上草地生態系統本身易受到外界環境干擾的特點,使得青藏高原草地區域具有典型的生態脆弱性。其脆弱性不僅表現在該地區的草地植被自我更新能力差,草地的自我恢復能力和抵抗外界來自于人類或自然的擾動能力低,而且當草地生態系統遭到破壞后,在重建、保護、改善、治理區域生態環境的過程中所付出的經濟成本和生態成本也明顯高于低海拔地區[11]。因此,有必要對青藏高原地區草地植被進行長時間序列的動態監測。

傳統的草地監測方法主要采用野外布設樣方進行采樣的方法,這種方法不僅效率低、成本高,而且青藏高原地區海拔高、部分地區生存環境惡劣,人力往往很難到達,所以無法快速、大范圍對草地植被生長狀況進行動態監測。遙感技術具有大范圍、快速、連續監測等優點,所以被廣泛用于區域尺度的草地植被長時間序列監測[12-14]。目前,很多學者開展了針對青藏高原草地退化的遙感監測研究。曹旭娟等人基于NDVI3g數據集計算了青藏高原長時間序列(1986—2013)草地退化狀況,發現青藏高原草地退化存在明顯的空間差異[15]。李重陽等人基于歸一化植被指數(normalized difference vegetation index, NDVI)數據研究了青藏高原牧戶草場退化趨勢,以及導致草場變化的主要驅動力,發現草場退化面積達到總面積的69%[16]。劉曉楓等人利用遙感影像獲取色達縣植被凈初級生產力(net primary productivity,NPP)及土壤有機質(soil organic matter,SOM),進而展開草地退化的研究[17]。在眾多植被、土壤指數中,NDVI能夠有效反應草地面積、植被高度和蓋度以及生產力的變化,是較優的草地退化監測指標[18]。基于遙感數據獲取植被監測參數,是進行草地退化大尺度研究的重要手段,也是目前研究草地退化的前沿趨勢。但目前針對青藏高原草地退化的研究集中在一些子區域,且大多停留在草地退化的時空特征分析,缺乏對高原草地退化的歸因分析。

地表溫度和土壤濕度都是陸地生態系統關鍵過程研究中重要的因子,一方面,在接收太陽輻射后,地表溫度升高,導致土壤水分蒸發進入大氣中,形成地表與大氣之間能量的交換；另一方面,土壤水分通過垂直運移,與地下水產生聯系,從而供給地表植被的生存用水[19]。地表溫度和土壤濕度在土壤的物理和生化過程中都起著重要的作用,同時二者對植物種子的萌芽和植株的生長也都有著重要影響[20-22]。水熱要素對于植被生長來說是至關重要的角色[23],在區域尺度上的相關研究主要是針對空氣溫度和降雨量[24-27],對于地表水熱要素的研究主要是集中在樣方尺度和點尺度[28-29],而少見區域尺度地表水熱要素對草地植被生長影響的綜合研究[30-31]。相比較于氣溫和降水,地表溫度和土壤濕度在土壤的物理和生化過程中都起著重要的作用,二者直接作用于植物的根部,對植物種子的萌芽和植株的生長也都有著重要影響。

基于此,本文利用遙感數據,獲取了地表溫度(land surface temperature, LST)和溫度植被干旱指數(temperature vegetation dryness index,TVDI)兩個地表水熱參數,并采用植被蓋度指標中的NDVI為主要遙感監測指標,建立草地退化遙感監測和評價指標體系。首先,對青藏高原草地退化狀況進行了時空變化監測；之后分析了青藏高原草地退化的時空特征和變化規律；最后,探究了地表水熱要素對草地退化的影響,從而加深對青藏高原地區氣候和植被變化機理的認識,具有一定的生態價值和現實意義。

1 研究區概況

青藏高原(73°15′E—104°47′E,26°00′N—39°47′N)西起帕米爾高原,東至橫斷山脈,北達昆侖山、祁連山脈,南接喜馬拉雅山脈,其東西長約2800 km,南北寬約300—1500 km,總面積約為250萬km2,平均海拔在4000 m以上。

青藏高原區域內高山大川密布,地勢險峻多變,地形復雜,其平均海拔遠遠超過同緯度周邊地區。青藏高原各處高山參差不齊,落差極大。總體來說,青藏高原地勢呈西高東低的特點。青藏高原年平均氣溫由東南的20℃,向西北遞減至-6℃以下。由于南部海洋暖濕氣流受多重高山阻留,年降水量也相應由2000 mm遞減至50 mm以下。青藏高原地區天然草地面積約為1.5×108hm2,并且在高原特殊的氣候環境影響下,形成了以高寒草甸和高寒草原為主草地類型。

2 數據源和預處理

2.1 NDVI數據

通過遙感數據獲得的NDVI時間序列數據能夠反映陸地生態系統植被的生長狀態、季相和年際變化特征,因此被廣泛用于全球尺度和區域尺度的生態環境變化監測與模擬、植被覆蓋動態變化研究、植被物候特征識別與信息提取等。本文使用的NDVI數據為2001—2017年青藏高原地區的MODIS影像數據MOD13A2產品,該產品為16 d合成產品,空間分辨率為1 km。數據預處理上,首先使用MRT (MODIS Reprojection Tool) 工具對上述數據進行批量投影變換和拼接處理,再使用ArcGIS將數據進行批量值轉換、水體掩膜、裁剪等處理。

NDVI由于表征的是植被生長狀態,通常反應植被生長過程的時間序列曲線是連續平滑的,但是由于傳感器本身性能、數據傳輸過程失誤、太陽光照角度、觀測視角、地物雙向性反射以及云、大氣氣溶膠等觀測條件因時間不同存在差異,此外還受到地表水、冰雪等因素干擾,得到的觀測值包含很多不可預測的噪聲,所以NDVI時間序列曲線呈鋸齒狀不規則波動,反映季節變化趨勢不明顯,對研究結果帶來很大的干擾。為此本文參考在青藏高原地區的相關研究,選取了適合該地區的NDVI時間序列重建方法[32],即非對稱性高斯函數擬合法(A—G擬合法)對2001—2017年的MODIS NDVI(16 天/1 km)數據集進行處理,最后通過平滑連接各高斯擬合曲線實現時間序列重建。

2.2 LST數據

地表溫度是大氣和地表連接處的溫度狀況,在土壤的物理和生化過程中都起著重要的作用,同時對植物種子的萌芽和植株的生長也都有著重要影響。作為與植物生長息息相關的環境要素,地表溫度的變化在很大程度上影響著地表植被的變化。

本文使用的LST數據為2001—2017年青藏高原地區的MOD11A2產品,該數據為8 d合成產品,空間分辨率為1 km。數據預處理上,首先使用MRT(MODIS Reprojection Tool)工具對上述數據進行批量投影變換和拼接處理,再使用ArcGIS將數據進行批量值轉換、水體掩膜、裁剪等處理,并將LST數據相鄰兩個時期的數據進行平均值計算處理,獲取與NDVI數據時間相匹配的16 d分辨率產品。

圖1 重采樣后的ERA5氣溫數據與中國氣象局CLDAS-V2.0氣溫數據對比Fig.1 Comparison between resampled ERA5 air temperature data and CLDAS-V2.0 air temperature data CLDAS：中國氣象局陸面數據同化系統 CMA Land Data Assimilation System；ERA5：第五代ECMWF大氣重新分析全球氣候 The Fifth Generation of ECMWF Atmospheric Reanalyzes of the Global Climate

2.3 中國氣象局陸面數據同化系統(CLDAS-V2.0)近實時產品數據集

CLDAS-V2.0產品為覆蓋亞洲區域(0—65°N, 60—160°E)經緯度網格融合再分析產品,其空間分辨率為0.0625°,時間分辨率為1 h,包括大氣驅動場產品(2 m氣溫、2 m比濕、10 m風速、地面氣壓、降水、短波輻射6個要素)、土壤濕度產品(垂直分為5層：0—5、0—10、10—40、40—100、100—200 cm)產品。該數據集利用多種來源地面、衛星等觀測資料,采用多重網格變分同化、最優插值、概率密度函數匹配、物理反演、地形校正等技術研制而成,在中國區域質量優于國際同類產品,且時空分辨率更高。由于該數據集在中國地區僅有2017年至今的數據集,時間范圍上不滿足本文需求,但是空間分辨率較高,所以本文采用該數據集中的2 m氣溫數據用作驗證數據。

2.4 ERA5數據集

ERA5(The Fifth Generation of ECMWF Atmospheric Reanalyzes of the Global Climate)是歐洲中期天氣預報中心(European Centre for Medium—Range Weather Forecasts, ECMWF)發布的第五代全球氣候大氣再分析網格數據集。再分析使用物理定律將模型數據與來自世界各地的觀測結果組合成一個全局完整且一致的數據集,以產生新的最佳估計大氣狀態,該數據集包括從2000年至今的每小時0.25°×0.25°的氣象因子數據,本文選取其中的2 m高度的氣溫數據(2 m temperature)進行處理。

首先下載與Terre衛星過境時間(當地時間10：30)最接近的氣溫數據(NetCDF格式),然后利用ArcGIS對數據進行格式轉換、平均值合成、重采樣等預處理獲取于NDVI時空分辨率相匹配的2 m氣溫數據。其中,重采樣方法使用雙線性內插法,該插值方法適用于表示某種現象分布、地形表面的連續數據,如 DEM(Digital elevation model,數字高程模型)、氣溫、降雨量分布、坡度等。本文將從中國氣象局陸面數據同化系統(CLDAS-V2.0)中獲取的氣溫作為基準值,來研究重采樣后ERA5數據的誤差。

圖1為2017年日序數DOY(Day of year)=161重采樣后的ERA5 氣溫數據與中國氣象局CLDAS氣溫數據散點圖,R2=0.75,相關性較好,同時通過計算求得回歸估計標準差RMSE為2 K,說明重采樣后的氣溫數據能夠較好的反應研究區域內的氣溫狀況,可以用于后續計算。

2.5 植被覆蓋度

植被覆蓋度基于像元二分模型進行計算。像元二分模型是一種簡單實用的遙感估算模型,它假設一個像元的地表由有植被覆蓋部分地表與無植被覆蓋部分地表組成,而遙感傳感器觀測到的光譜信息也由這2個組分因子線性加權合成,各因子的權重是各自的面積在像元中所占的比率,如其中植被覆蓋度可以看作是植被的權重。

FVC=(NDVI-NDVIsoil)/(NDVIveg-NDVIsoil)

(1)

式中,NDVIsoil作為裸土的NDVI值,理論上應該趨近于0,但是由于遙感影像受到大氣環境和粗糙地表、土壤顏色、土壤地表濕度等環境的影響,NDVIsoil值不是固定值,會在一定范圍內變化,一般為-0.1—0.2。NDVIveg為研究區與像元最大NDVI值。鑒于此,基于對研究區域內像元NDVI灰度值的統計分析,截取置信區間累計頻率在5%—95%對應的NDVI值分別作為NDVI最大值和最小值,從而計算得到植被覆蓋度FVC。

3 研究方法

3.1 青藏高原草地退化遙感監測和評價指標體系的建立

根據國家標準GB19377—2003天然草地退化、沙化、鹽漬化的分級指標,20世紀80年代初期相同監測區域相同草地類型的草地植被特征可以作為退化草地的基準。因此,以1982—1985年每個像元的年際最大植被覆蓋度(AVHRR數據計算得到)作為基準,將其退化程度分為未退化、輕度退化、中度退化、重度退化和極重度退化5個級別[33],具體見表1。

表1 草地退化等級

3.2 區域退化指數

在草地退化等級劃分的基礎上,采用Gao等人提出的草地退化指數(grassland degradation index,GDI)。草地退化指數的計算公式為：

(2)

式中,GDI為區域草地退化指數,Di為草地退化等級i的評分,Ai為草地退化等級i的分布面積,A為研究區與的總面積。依據退化分級標準分析草地退化情況,劃分標準如表2所示。

表2 區域草地退化指數

3.3 改進溫度植被干旱指數

土壤濕度作為一個重要的環境參數,是求解區域水量平衡、研究全球水循環過程的關鍵因子。一方面,土壤水分通過接收太陽輻射,地表溫度升高,使水分蒸發進入大氣中,形成地表與大氣之間能量的交換；另一方面,土壤水分通過垂直運移,與地下水產生聯系,從而供給地表植被的生存用水,將地表水與地下水聯系在一起。此外,土壤濕度的時空分布與動態變化不僅會影響其自身,對地表反照率、蒸散發、生物化學循環、植被生長等也有著重要的作用。

1994年,Carlson等人在研究不同分辨率的LST和NDVI數據時,發現LST與NDVI之間呈現出明顯的負相關關系[34]。之后隨著不斷研究發現,當研究區域內的植被覆蓋和土壤濕度變化較大時,通過遙感影像獲得的NDVI和LST所構成的散點圖呈梯形或三角形。

實驗中,當研究區域足夠大,而且土地覆蓋類型從裸土變化到完全植被覆蓋,土壤濕度從干旱變化到濕潤,NDVI和LST構成的空間關系為三角形。如圖2所示,A點表示干燥裸土,即NDVI小,LST高；B點表示濕潤裸土,即NDVI小,LST也小；C點表示植被完全覆蓋,同時土壤水分也很充足,即NDVI大,LST低,此時蒸散阻抗小(蒸散包括蒸發和植物的蒸騰,蒸發又包括裸土的蒸發和植被冠層截獲水分的直接蒸發)。AC邊表示,對于最高溫度下,不同植被覆蓋類型土壤濕度最低的狀況,我們稱這條邊為“干邊”。BC邊表示,在這一溫度下,土壤水分充足,不會成為植被生長的限制因素,我們稱之為“濕邊”。

圖2 溫度植被干旱指數模型Fig.2 Temperature vegetation dryness index modelA點：裸土且土壤濕度最小；B點：裸土且土壤濕度最大；C點：全植被覆蓋；AC：理論干邊；BC：理論濕邊

當水分條件充足時,植被生長迅速,NDVI值變高,同時區域內植被生命活動旺盛,蒸騰量變大,整個像元內的蒸散阻抗變小,潛熱所占比例增大,所以像元內的地表溫度會降低；當水分不足時,植物受到水分脅迫,植物葉片氣孔關閉,降低了蒸騰所造成的水分損失,進而導致地表潛熱通量的降低,根據能量平衡原理,地表能量必須平衡,所以地表感熱通量會增加,感熱通量的增加最終會導致冠層溫度的升高。

Sandholt等在研究利用地表溫度—植被指數三角形特征空間研究土壤濕度時,提出了利用溫度植被干旱指數(TVDI)來表示土壤干濕狀況[35],定義為：

(3)

式中,T表示任意像元的地表溫度值；Tmax表示某一NDVI對應的最大溫度,可以由NDVI與干邊線性擬合得到；Tmin表示某一NDVI對應的最小溫度,可以由NDVI與濕邊線性擬合得到。TVDI的取值范圍為0—1,TVDI值越大,說明這一像元的土壤濕度越低,TVDI值越小,說明這一像元的土壤濕度越高。

在三角形特征空間中對干濕邊進行線性擬合,可以得到干邊方程為：

Tmax=a1+b1×NDVI

(4)

濕邊方程為：

Tmin=a2+b2×NDVI

(5)

將式(4)和式(5)代入式(3)中,可以得到TVDI的計算公式為：

(6)

式中,a1,b1和a2,b2分別是干邊和濕邊的線性擬合的系數。

如圖2所示,在AC干邊上,不同植被指數下的溫度最高,土壤含水量低,NDVI低的點干燥迅速,蒸發量少,吸收較多的太陽能,為干旱狀態,TVDI為1；在BC濕邊上,土壤含水量高,吸收的太陽能主要用來蒸發和蒸騰作用,裸土表面和作物冠層溫度差距不明顯,BC邊接近水平線,為濕潤狀態,TVDI為0。

由于某個區域某一段時間內不同像元的NDVI與LST值都分布在干濕邊界構成的特征空間內,所以我們可以將這個特征空間看作是由土壤濕度等值線組成。地表溫度—植被干旱指數TVDI可以用來表征區域尺度地表土壤含水量情況,也可以稱為溫度植被干旱指數,下文中將統一用溫度植被干旱指數表示[36]。

在TVDI模型中,假設TVDI變化的主要來源是土壤水分,沒有考慮空氣溫度。當TVDI進行反演和應用時,通常假設空氣溫度在研究區域內是固定的,這在大區域進行應用時明顯是錯誤的[37]。

在本文中,青藏高原地區面積較大,且地面異質性較大,所以考慮將空氣溫度引入到TVDI方程中,參照Parinaz等的相關研究[38],可以得到改進后的溫度植被干旱指數計算公式為：

(7)

構建TVDI特征空間時,本文采用的算法將NDVI每間隔0.025劃分步長,提取同一NDVI數值下的最大最小溫度(置信區間95%),根據最大最小溫度提取特征空間的干邊和濕邊,考慮到NDVI在0.1以下和0.9以上對植被生長狀況指示存在不靈敏或飽和的情況,本文在擬合干濕邊界時,只提取NDVI在0.1—0.9之間的值。

如圖3中所示,由于選取驗證的數據的時間處于6月份,此時青藏高原地表溫度較高,地表土壤含水量較低,對應的TVDI較高,所以散點集中分布在右下部分。另外,改進后的TVDIi比TVDI在1：1反比線左右更加聚合。綜上可得改進后的溫度植被干旱指數TVDIi在研究區域更加適用,更能夠反應土壤水分的變化。

圖3 TVDI與改進TVDI干濕邊擬合Fig.3 The fitting of dry and wet edge of the TVDI and improved TVDI modelTVDI：溫度植被干旱指數,temperature vegetation dryness index；TVDIi：改進后的溫度植被干旱指數,improved temperature vegetation dryness index

以2009年DOY=161為例,構建改進后的溫度植被干旱指數的特征空間時,考慮到NDVI在0.1以下和0.9以上對植被生長狀況指示存在不靈敏或飽和的情況,本文在擬合干濕邊界時,只提取NDVI在0.1—0.9之間的值。通過圖3和表3得知改進后的TVDI對干濕邊界的擬合效果更好。

表3 TVDI和TVDIi的干濕邊方程

3.4 變化趨勢分析方法

Slope定義為在一定時間范圍內采用最小二乘法擬合年度變量均值的斜率,它能夠反映出每個柵格的變化趨勢。計算公式如下：

(8)

式中,i為年份序號,取值范圍為1—17,n為研究的時間序列長度,為17,Xi為第i年的年度平均值(例如當X為NDVI時,這里表示的則為NDVI的年度平均值)。若Slope大于0則說明該像元的變量X在這17年間的變化趨勢是增加的,反之則是減少的。

3.5 多元線性回歸分析方法與貢獻率

Yi=β0+β1·X1i+…+βk·Xki+μi

(9)

式中,Yi為因變量；Xki為自變量；i=1, 2, …,n；k為解釋變量的個數,βk為回歸系數。

本章節將以地表溫度LST和溫度植被干旱指數TVDI(改進后的TVDI,下同)兩個因子作為回歸模型的自變量,然后將NDVI作為因變量,來建立回歸模型。參考Zhu等人的研究[39],如果在建立回歸模型之前,將遙感柵格數據均進行歸一化處理,然后用17年的年尺度歸一化數據建立回歸模型,最后求得的每個自變量的系數則為該像元該自變量對因變量的定量貢獻率,并且可以通過對比不同自變量定量貢獻率的大小來決定主要影響因子。

首先對年尺度的數據進行歸一化處理。之后基于MATLAB建立多元線性回歸模型,逐像元計算LST和TVDI二元變量的系數及常數項,將系數和常數項分別出圖,得到LST和TVDI對因變量的貢獻率空間分布圖像及反應其他因子貢獻率的殘差圖像。

4 研究結果與分析

4.1 青藏高原草地退化現狀分析

從圖4,2001和2017年兩年的青藏高原地區草地退化等級空間分布可以看出,青藏高原植被退化程度空間差異明顯,柴達木盆地和青海湖附近退化較為嚴重,喜馬拉雅山脈北部、昆侖山脈南部、岡底斯山脈北部交匯的地區退化也較嚴重。但是2017年總體退化程度較2001年輕。

圖4 2001和2017年草地退化空間格局Fig.4 Spatial pattern of grassland degradation degrees in 2001 and 2017

從2001和2017年兩年的青藏高原地區植被退化等級占比隨時間變化圖可以看出來(圖5,圖6),2001—2017年每年植被未退化面積基本處于50%—60%之間,且不同退化等級的面積隨著時間的變化也均有變化,其中未退化面積從50.6%上升到59%,說明青藏高原植被退化整體上在朝著改善的方向發展。

圖5 青藏高原地區草地各退化等級面積占比Fig.5 The area proportion of different grassland degradation degrees on the Tibetan Plateau

圖6 青藏高原地區未退化草地面積占比Fig.6 The proportion of non-degraded grassland area on the Tibetan Plateau

GDI指數主要用來計算某個區域整體退化狀況。圖7為青藏高原草地區域2001—2017年的每年草地區域的GDI指數。整體上,青藏高原草地區域GDI指數大部分處于1.5—2之間,說明大部分時間處于輕度退化狀況,而2001年和2015年GDI均達到了2,趨向中度退化,說明這兩年草地植被生長出現了惡化。但是整體上,從2001—2017年,草地區域GDI指數呈減小的趨勢,說明青藏高原草地生長這些年之間退化情況在減輕。

圖7 青藏高原地區草地GDI年度變化Fig.7 Annual change of grassland GDI on the Tibetan Plateau

4.2 青藏高原地表水熱要素對草地退化的影響

參考Zhu等人的研究[39],如果在建立回歸模型之前,將遙感柵格數據均進行歸一化處理,然后用17年的年尺度歸一化數據建立回歸模型,最后求得的每個自變量的系數則為該自變量對因變量的定量貢獻率,并且可以通過對比不同自變量定量貢獻率的大小來決定主要影響因子。通過逐像元計算方式,可以獲得各因變量的貢獻率在空間上的分布情況。首先計算TVDI與LST對NDVI的二元回歸關系,其回歸系數可以反映因子對植被影響的大小(圖8)。當回歸系數在-0.5—0.5之間,可以認為因子對植被幾乎沒有影響；當回歸系數為-1.0至-0.5,0.5至1,因子對植被有一般程度的影響；當回歸系數在小于-1或大于1因子對植被具有較顯著的影響。地表溫度對NDVI的回歸系數主要為負數,說明大部分區域地表溫度與植被NDVI的關系為反向關系,而對TVDI,與植被生長狀況主要呈現正向關系。

圖8 LST、TVDI對植被NDVI的回歸系數空間分布 Fig 8 The spatial distribution of the regression coefficients of LST and TVDI to NDVILST：地表溫度land surface temperature；TVDI：溫度植被干旱指數 temperature vegetation dryness index

圖9 青藏高原草地區域受地表水熱因子影響的空間分布Fig.9 Spatial distribution of grassland affected by surface hydrothermal factors in Tibetan Plateau

通過比較TVDI、LST回歸系數與常數項的大小,確定青藏高原區域各像元的主導因子,結果如圖9所示。其中,在青藏高原東南部,水熱因子對植被的影響相對較小。水分主導的區域主要分布在青藏高原北部,相對較干旱的區域。其次,地表溫度在大部分地區都存在與植被的顯著相關關系。從草地植被類型角度,統計了水熱因子在各類型草地中主導區域的面積比例,其結果如圖10所示。地表溫度在各個植被中均占據較大的主導面積。地表溫度在溫性草原、高寒草原、溫性草甸、高寒草甸中占據主導面積的比例分別為32.73%,42.24%,34.59%,26.70%。TVDI在溫性草原、高寒草原、溫性草甸、高寒草甸中占據主導面積的比例分別為7.83%,16.92%,8.65%,10.80%。

5 討論

從計算結果來看,地表溫度對植被NDVI具有較顯著的回歸關系,但是很難說明地表溫度對植被具有很高的影響。地表溫度與地表湍流能量分配有關,而地表湍流能量分配主要取決于植被覆蓋度和土壤濕度[40]。因此,在特定的土壤水分含水量下,植被NDVI和LST之間存在很強的線性關系,這一關系也被用于發展計算土壤水分虧缺狀況的方法,如三角形特征空間模型[41]。為了更深一步認識LST與NDVI之間的強相關關系,本文選取了2014年的NDVI與LST數據,通過對NDVI按照0.005間隔劃分,計算相應區間內LST均值,得到圖11統計結果。從圖11可以看出,總體上,NDVI越大,對應的LST越小。由于遙感觀測的LST包括土壤組分和植被組分兩部分,相對于裸土,在相同的輻射下,植被組分的溫度更低。因此,盡管回歸模型中,LST對NDVI的回歸關系較優,LST對植被生長的關系也更為復雜。在低植被覆蓋區域,LST主要組分來自裸土的情況下,地表溫度的增高有利于植被的生長,同時,在高植被覆蓋區域,植被覆蓋度會反過來影響地表溫度,隨著覆蓋度的增大,地表溫度降低。從回歸系數來看(圖8),LST對NDVI的回歸系數也偏向于負數,說明了LST與NDVI之間較強的負相關關系。進一步地,建議獲取土壤組分溫度,從而深入研究區域尺度上土壤溫度對草地植被生長的影響。

圖10 青藏高原草地區域受地表水熱因子影響的面積占比Fig.10 The proportion of grassland affected by surface hydrothermal factors in Tibetan Plateau

圖11 2014年青藏高原NDVI與地表溫度的關系 Fig.11 The relationship between NDVI and land surface temperature over Tibetan Plateau in 2014NDVI：歸一化植被指數,Normalized difference vegetation index

青藏高原在2001—2017年間,極重度退化和重度退化面積比例變化很小,說明重度及以上退化一般呈現不可逆的退化趨勢,短時間內很難恢復。同時,未退化草地面積呈現逐年上升趨勢,說明青藏高原草地在2001—2017年間的退化狀況正在改善,且主要是從輕度、中度退化過渡到未退化。如圖4所示,青藏高原東北部地區(主要分布在海西蒙古族藏族自治州,柴達木沙漠周圍地區)呈現出較為明顯的草地退化改善趨勢。從圖8,這些地區的TVDI與NDVI具有較高的相關系數,土壤濕度的增大可能是導致地區草地退化改善的主要原因,其次,地表溫度與植被NDVI較高的負相關關系,也體現出草地植被NDVI的增大趨勢。草地退化同時受多種因素影響,除了自然因素,人類活動在玉樹藏族自治州、甘南藏族自治州等以牧場草地為主的地區也有重要的影響[15-16,42]。由于近年來有規劃的合理放牧,2000年至2016年,牧場草地的退化面積也呈現明顯的上升趨勢[16]。從植被類型角度,受地表水熱影響較大的草地類型為高寒草原,其次為溫性草甸、溫性草原和高寒草甸。高寒草原在青藏高原的分布面積最大,相對于其他草地類型,受人類活動影響的面積比例也較低,因此,地表水熱等自然因素對該類型草地影響較大。總體來看,地表水對草地植被的恢復具有積極影響,地表含水量同時也反映了土壤的持水能力和土壤質地狀況[43],而草地退化的根本原因之一可能為土壤的退化。

同時,在研究過程中仍存在需要完善和值得深入思考的地方,比如本文主要研究內容為青藏高原植被退化對地表水熱要素的響應(LST、TVDI),其他的影響因子如太陽輻射、人類活動等對植被生長均有影響,本文也將除地表水熱因子以外的影響因子統稱為“其他影響因子”進行了分析,在今后的研究中可以考慮將更多的影響因子單獨納入分析中。同時,由于青藏高原地區范圍廣、氣象監測站點少、自然環境惡劣,其他因素比如氣象因子或人為因子數據源難以和遙感數據進行空間尺度匹配,因此在今后的研究中應該加強各種區域尺度數據的獲取,開發適合青藏高原復雜下墊面的氣象要素插值方法,研究將縣、市、省級的人為影響因子進行尺度轉化的方法,從而與常見的遙感數據尺度進行匹配。

6 結論

(1)研究結果表明從2001到2017年間,青藏高原植被退化程度空間差異明顯,柴達木盆地和青海湖附近退化較為嚴重,喜馬拉雅山脈北部、昆侖山脈南部、岡底斯山脈北部交匯的地區退化也較嚴重。

(2)從2001—2017年間,青藏高原草地未退化面積從50.6%上升到59%,說明青藏高原植被退化整體上在朝著改善的方向發展。2001—2017年內,大部分年份GDI指數在1.5—2之間,說明青藏高原草地大部分時間處于輕度退化狀態,其中,2001年和2015年,GDI達到了2以上,這兩個年份青藏高原草地退化達到中等退化水平。整體上看,GDI呈現減小趨勢,青藏高原草地退化在減輕。

(3)青藏高原草地植被14.04%的區域主要受到土壤水分含量影響,從回歸關系來看,36.61%的區域NDVI與地表溫度有較強的相關關系,但主要反映在植被覆蓋影響地表溫度。綜合來看,青藏高原地表水熱對草地植被退化具有顯著影響。

致謝：國家氣象科學數據中心(http://data.cma.cn/site/index.html)提供氣象數據產品；美國國家航空航天局土地流程分布式活動檔案中心提供MODIS影像數據(https://search.earthdata.nasa.gov/search)；歐洲中期天氣預報中心提供再分析氣象數據產品(https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels)；中國科學院空天信息創新研究院遙感科學國家重點實驗室生態系統遙感研究室提供青藏高原覆被數據,特此致謝。