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

夏 龍,宋小寧,蔡碩豪,胡容海,郭 達(dá)

中國科學(xué)院大學(xué)資源與環(huán)境學(xué)院, 北京 100049

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

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

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

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

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

1 研究區(qū)概況

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

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

2 數(shù)據(jù)源和預(yù)處理

2.1 NDVI數(shù)據(jù)

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

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

2.2 LST數(shù)據(jù)

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

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

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

2.3 中國氣象局陸面數(shù)據(jù)同化系統(tǒng)(CLDAS-V2.0)近實(shí)時(shí)產(chǎn)品數(shù)據(jù)集

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

2.4 ERA5數(shù)據(jù)集

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

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

圖1為2017年日序數(shù)DOY(Day of year)=161重采樣后的ERA5 氣溫?cái)?shù)據(jù)與中國氣象局CLDAS氣溫?cái)?shù)據(jù)散點(diǎn)圖,R2=0.75,相關(guān)性較好,同時(shí)通過計(jì)算求得回歸估計(jì)標(biāo)準(zhǔn)差RMSE為2 K,說明重采樣后的氣溫?cái)?shù)據(jù)能夠較好的反應(yīng)研究區(qū)域內(nèi)的氣溫狀況,可以用于后續(xù)計(jì)算。

2.5 植被覆蓋度

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

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

(1)

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

3 研究方法

3.1 青藏高原草地退化遙感監(jiān)測和評(píng)價(jià)指標(biāo)體系的建立

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

表1 草地退化等級(jí)

3.2 區(qū)域退化指數(shù)

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

(2)

式中,GDI為區(qū)域草地退化指數(shù),Di為草地退化等級(jí)i的評(píng)分,Ai為草地退化等級(jí)i的分布面積,A為研究區(qū)與的總面積。依據(jù)退化分級(jí)標(biāo)準(zhǔn)分析草地退化情況,劃分標(biāo)準(zhǔn)如表2所示。

表2 區(qū)域草地退化指數(shù)

3.3 改進(jìn)溫度植被干旱指數(shù)

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

1994年,Carlson等人在研究不同分辨率的LST和NDVI數(shù)據(jù)時(shí),發(fā)現(xiàn)LST與NDVI之間呈現(xiàn)出明顯的負(fù)相關(guān)關(guān)系[34]。之后隨著不斷研究發(fā)現(xiàn),當(dāng)研究區(qū)域內(nèi)的植被覆蓋和土壤濕度變化較大時(shí),通過遙感影像獲得的NDVI和LST所構(gòu)成的散點(diǎn)圖呈梯形或三角形。

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

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

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

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

(3)

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

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

Tmax=a1+b1×NDVI

(4)

濕邊方程為：

Tmin=a2+b2×NDVI

(5)

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

(6)

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

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

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

在TVDI模型中,假設(shè)TVDI變化的主要來源是土壤水分,沒有考慮空氣溫度。當(dāng)TVDI進(jìn)行反演和應(yīng)用時(shí),通常假設(shè)空氣溫度在研究區(qū)域內(nèi)是固定的,這在大區(qū)域進(jìn)行應(yīng)用時(shí)明顯是錯(cuò)誤的[37]。

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

(7)

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

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

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

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

表3 TVDI和TVDIi的干濕邊方程

3.4 變化趨勢分析方法

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

(8)

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

3.5 多元線性回歸分析方法與貢獻(xiàn)率

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

(9)

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

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

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

4 研究結(jié)果與分析

4.1 青藏高原草地退化現(xiàn)狀分析

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

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

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

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

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

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

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

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

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

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

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

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

5 討論

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

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

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

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

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

6 結(jié)論

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

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

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

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