基于多源遙感數據的青藏高原植被變化特征評價

趙耀華，彭小清，2，金浩東， 杜 冉，陳 聰，彭思佳

（1.蘭州大學資源環境學院，甘肅蘭州 730000；2.蘭州大學祁連山凍土生態環境野外科學觀測研究站，甘肅蘭州 730000）

0 引言

全球平均表面溫度持續上升，青藏高原氣候變暖更為顯著。近幾十年來，青藏高原升溫速率約為全球同期升溫速率的2倍，達到每10年升高0.3～0.4℃［1］?？焖俚脑鰷刈兣o作為我國重要生態屏障的青藏高原帶來巨大的生態安全挑戰［2］。一方面，植被是陸地生態系統的必要成分，是生態系統中最為活躍的因素，具有連接土壤、大氣、水分與其他生態環境因子的橋梁作用，其變化影響著地球重要圈層的物質循環與能量流動。另一方面，生態系統的變化可以被敏感的植被所反映［3］?？茖W地評估青藏高原植被變化有助于維護高原地區生態平衡以及恢復高原生態系統，以應對全球氣候變化帶來的負面影響。在生態文明加速發展的背景下，《生態文明體制改革總體方案》的部署也明確指出青藏高原生態屏障修復的重要性［4］。正確認識植被的變化情況將促進區域經濟與生態文明協調發展。

全球環境的變化正在迅速影響著陸地植被動態，但植被應對快速變化的環境的機制尚不明確。因此，植被生態對于氣候變化的反饋與響應受到專家學者持續性的關注。青藏高原屬于高寒生態系統，主要覆蓋著高寒草原與高寒草甸，相比其他生態系統更為脆弱，也由此成為相關研究的重點［5］。近些年來，對于青藏高原的生態系統的研究主要圍繞著植被物候、植被綠度、林線位置、生物多樣性、生態系統生產力與碳匯等主題［2，5-6］。其中植被綠度及物候變化是生態系統監測的重要內容，利用衛星遙感的衍生指標，比如歸一化植被指數（NDVI），評估青藏高原植被長期狀況，揭示植被變化的研究日漸增多。Yu等［7］通過GIMMS NDVI對青藏高原植被生長季動態進行研究，發現在氣溫的影響下，高原春季物候在20世紀90年代中期發生了延遲。Chen等［8］雖然表示并未在野外檢測中觀察到相關現象，但整合了草原退化導致的變冷效應及凍土凍結提前對春季物候延遲現象的解釋。Yi等［9］認為可能是青藏高原春季的較高氣溶膠指數影響了NDVI的監測值，造成春季物候延遲的現象。Piao等［10］也使用GIMMS NDVI得到了青藏高原春季物候在1982—1999年以0.88 d·a-1的速率提前，隨后春季物候出現了延遲，與Yu等［7］的研究結果一致。Zhang等［11］基于MODIS、SPOT與GIMMS NDVI數據集，揭示在1982—2011年間青藏高原春季物候以1.04 d·a-1的速度持續提前。并且，Zhang等［11］對三種NDVI數據集進行了不同時間尺度的年際變化分析，認為GIMMS NDVI存在質量問題。但Shen等［12］在隨后利用SPOT與MODIS NDVI指出，盡管青藏高原近10年春季溫度上升，但在區域尺度上植被的返青期并未顯著提前。關于青藏高原春季物候的提前、靜止或延遲，學界一直缺乏共識［13］，不同研究對于青藏高原植被評價的結論存在差異性，而這種不一致性很可能與研究所選取的衛星產品相關［11，14］。不同時間、空間分辨率的多源NDVI數據通常作為直接表征區域植被綠度的指標，用來探究植被變化，同時也是植被物候提取的數據來源，因此基于多源NDVI數據集對青藏高原地區植被變化特征進行評價并探討不同NDVI數據之間的差異性尤為重要。但是，現有研究對于多源遙感數據在植被變化特征差異性方面的探討缺少針對性：一方面，研究更多地探討不同衛星衍生指標之間的差異，而非聚焦同種指標［14-16］，如歸一化植被指數之間的差異；另一方面，即使基于歸一化植被指數進行研究，所使用的NDVI數據集種類并不全面［11，17］，研究范圍也往往并非針對青藏高原［18-21］，而光譜植被指數受到的區域限制因素眾多，存在極大的空間異質性。所以，本研究采用多種不同時空尺度的歸一化植被指數在青藏高原范圍內，對植被狀況及其長時間變化進行分析，探討和對比不同數據之間的差異性。

鑒于青藏高原植被在地-氣能量交換、水文過程等的重要作用，而過去研究缺乏多源數據在青藏高原植被變化評估中的差異性和不確定性，為此本文圍繞不同時空分辨率的多源衛星遙感數據在青藏高原植被變化評估的差異性如何的科學問題開展相關研究?；贛ODIS（250 m、500 m、1 km、5.5 km）、GIMMS 8 km和SPOT 1 km NDVI數據集，開展不同時空尺度下青藏高原地區的植被變化特征評價。具體包括：使用最大值合成方法（MVC）探究青藏高原植被2000—2014年的空間分布；使用Theil-Sen趨勢估計和Mann-Kendall趨勢檢驗探究時間變化特征；最終評價不同時間空間分辨率的多源衛星遙感數據在青藏高原植被變化特征研究上的差異性。

1 研究區概況

青藏高原是繼南北極后的“第三極”，平均海拔高度超過4 000 m，面積約為2.5×106km2，是世界上海拔最高、面積最大的高原［22］。青藏高原范圍為73°29′56″～104°40′20″E，25°59′37″～39°49′33″N［23］，山脈縱橫，地域廣闊。青藏高原氣溫差異大，年平均氣溫為-10～15℃［24-25］［圖1（a）］。自1901年以來，氣溫增加趨勢最高可達0.15℃·（10a）-1以上［圖1（b）］，高原西北部增溫趨勢更為顯著。降水主要集中在夏季和秋季，而冬、春兩季相對干燥［26］。降水趨勢在整個高原尺度上表現出較小的季節和空間波動，但總體上略有增加［26］。年降水量由東南向西北變化，從高于1 000 mm減少至低于100 mm［27］。復雜的氣候特征造就了高原豐富的植被類型分布。相對更為溫暖濕潤的高原東南部與東部邊緣，主要分布有森林和灌木，向西植被逐漸發展為高寒草甸，干燥的高原中部與東北部區域覆蓋有溫帶與高寒草原，荒漠則分布在高原西北部的干旱地區［26］。

圖1 青藏高原1901—2018年青藏高原年平均氣溫（網格溫度數據來自CRU）（a）及其變化趨勢（b）Fig.1 Mean annual air temperature（grid temperature data is from CRU）of the Tibetan Plateau from 1901 to 2018（a）and its change trend（b）

2 數據與方法

2.1 數據來源

本文研究數據來自MODIS、SPOT-VGT、AVHRR衛星遙感源傳感器，涉及包括250 m、500 m、1 km、5.5 km、8 km共5種空間分辨率，跨越了從10-2km2至近102km2共5個數量級的空間尺度（表1）。

表1 衛星遙感數據集介紹Table 1 Description of satellite remote sensing datasets

本文使用搭載于EOS/Terra衛星上的中等分辨率成像光譜儀（MODIS）反演的植被指數產品，包括MOD13Q1、MOD13A1、MOD13A2、MOD13C1產品，空間分辨率為250 m、500 m、1 km、5.5 km，其提供了全球的植被狀況與具體時空變化的信息，可以用于觀測陸地植被光合作用的變化。MOD13Q1、MOD13A1與MOD13A2數據集來源于Google Earth Engine平臺，MOD13C1數據來源于美國宇航局（https：//ladsweb.modaps.eosdis.nasa.gov/）。GIMMS NDVI 3g數據（http：//ecocast.arc.nasa.gov/data/）由NOAA衛星搭載的AVHRR傳感器獲取，是目前在植被研究中使用最廣，涵蓋時間序列最長的NDVI產品。SPOT NDVI產品（http：//land.coperni?cus.eu/global/products）為GIOGL1_ATBD_ND?VI1km-V2.2數據集，數據由SPOT-VEGETATION提供（2014年后由PROBA-V提供），基于最大值合成法合成以10天為周期的最佳NDVI值，減少了噪聲的影響。

2.2 研究方法

2.2.1 最大值合成法與標準差

使用最大值合成法分別將不同時空分辨率的MODIS、SPOT、GIMMS NDVI數 據 合 成2000—2014年多年季節性NDVI［28］，以探究青藏高原多年季節性植被空間分布狀況。為了保證不同數據源NDVI之間的可比性，反映不同數據來源NDVI的差異性，使用最近鄰插值將所有數據集的空間分辨率轉換為250 m［29］，計算6種不同NDVI數據集之間的標準差［14］。

2.2.2 Theil-Sen趨勢估計和Mann-Kendall趨勢檢驗

由于NDVI時間序列通常不滿足數據獨立且正態分布的參數趨勢檢驗要求［15］，在探究青藏高原植被變化時，使用了兩種非參數方法——Theil-Sen趨勢估計和Mann-Kendall趨勢檢驗，以探究NDVI序列變化的趨勢［30］。

研究使用Theil-Sen斜率估計［31］作為NDVI隨時間變化程度的衡量指標，更為準確地定量化描述植被變化程度［32］。Theil-Sen斜率估計是一種穩健的趨勢統計方法［30］，雖然與線性最小二乘回歸方法相關，但該方法基于非參數統計，通過中值計算斜率，與前者相比能夠抵抗噪聲的干擾，對數據中的離群值有著良好的規避能力［6，18］。

NDVI時間序列趨勢的顯著性使用非參數的Mann-Kendall趨勢檢驗，Mann-Kendall趨勢檢驗比較樣本數據的相對大小并計算得到統計量S及其方差Var（S），然后計算檢驗統計量Z，來確定數據變化的趨勢顯著性［15，33］。

Theil-Sen趨勢估計和Mann-Kendall趨勢檢驗相結合能更加有效地探究NDVI的時間序列變化，這兩種方法被廣泛應用于地球科學的許多領域［34］。

3 結果與分析

3.1 青藏高原NDVI空間分布特征

通過對不同遙感數據集的青藏高原2000—2014年季節性NDVI最大值空間分布的分析，結果表明，青藏高原區域植被分布情況有著明顯的空間異質性。各季節青藏高原NDVI空間分布均呈現東高西低的特征，NDVI由東南到西北逐漸降低，空間分布格局明顯（圖2）。青藏高原植被分布格局的形成可能與植被對于溫度的敏感性相關，高原年均氣溫的分布同樣是東南部相對較高。雖然植被生長與氣溫之間的關系十分復雜且具有顯著的區域差異［35］，但對于高緯度地區與青藏高原來說，氣溫依然被廣泛認為是植被生長的重要限制因素。

圖2 青藏高原生長季NDVI的空間分布Fig.2 Spatial distribution of NDVI in the growing season on the Tibetan Plateau

研究區NDVI在0～1的范圍內均有分布（圖2）。不同季節間，NDVI的主要分布范圍有所差異。青藏高原夏季NDVI、秋季NDVI與生長季NDVI均呈現雙峰式集中分布，夏季NDVI與生長季NDVI在0～0.3與0.7～0.9處集中。秋季NDVI主要在0～0.3與0.6～0.8相對集中，但集中在0.6～0.8的像元占比遠小于前者。春季NDVI整體較低，約58.06%的像元NDVI集中在0～0.2范圍。

不同NDVI數據集之間的標準差由北向南逐漸增加（圖3）。秋季的高原整體標準差最小。春季，6種NDVI數據集差異最為明顯，較大的標準差出現在高原東南部及南部。在夏季與生長季，不同ND?VI數據集主要在高原南部出現較大的差異。SPOT NDVI分布更多的集中在0.1～0.4之間（圖4）。SPOT NDVI在0.1～0.4之間像元占比最大，平均多于其他數據集4.11%。在0.4～0.7的NDVI分布區間內，GIMMS NDVI較其他數據平均多約1.83%。而MODIS數據普遍在0.7～1區間的NDVI像元占比更大。6種NDVI數據集中，250 m分辨率的MODIS NDVI在0.7以上分布占比最具優勢。而隨著空間分辨率的降低，MODIS NDVI數據在NDVI>0.7的面積逐漸減少。

圖3 6種NDVI標準差的空間分布Fig.3 Spatial distribution of standard deviation of six NDVI：spring（a），summer（b），autumn（c）and growing season（d）

圖4 青藏高原各季節NDVI的頻率分布Fig.4 Frequency distribution of NDVI in each season on the Tibetan Plateau

3.2 青藏高原NDVI時間變化特征

3.2.1 區域平均NDVI變化趨勢

通過對6種NDVI數據集在青藏高原范圍內進行區域平均計算和Theil-Sen趨勢估計分析（圖5）。結果表明，在青藏高原的區域尺度上，SPOT NDVI與MODIS 250 m、500 m、1 km NDVI數據集除秋季外，均表現顯著的綠化趨勢（表2），MODIS 5.5 km與GIMMS NDVI數據集分別在春季、夏季、生長季呈現顯著綠化。

表2 青藏高原平均NDVI變化趨勢［單位：（10a）-1］Table 2 Trends in average NDVI changes on the Tibetan Plateau［unit：（10a）-1］

圖5 青藏高原平均NDVI變化Fig.5 Average NDVI changes on the Tibetan Plateau

對于不同數據源而言，NDVI的顯著變化趨勢相差較大，從0.0092（10a）-1到0.0184（10a）-1。即使季節相同，不同NDVI數據集之間的顯著綠化趨勢差異最大可以達到0.0090（10a）-1（生長季）。SPOT NDVI在四個季節均表現顯著的綠化趨勢，且平均綠化趨勢達到0.0162（10a）-1，高于其他數據平均變化趨勢0.0085（10a）-1，綠化顯著且迅速。GIMMS NDVI僅在2000—2014年間的生長季變化分析中具有0.0092（10a）-1的顯著綠化趨勢。具有最高空間分辨率的MODIS 250 m NDVI數據集往往表現出最小的顯著變化趨勢值。

就季節而言，除秋季外其他季節均具有顯著的變化趨勢。各個NDVI數據集在15年間春季的顯著變化趨勢介于0.0098～0.0184（10a）-1之間。多數數據集夏季顯著綠化趨勢明顯高于其他季節，夏季平均顯著綠化趨勢達到0.0140（10a）-1，而青藏高原夏季的NDVI往往深刻影響整個生長季的植被變化。

3.2.2 青藏高原多年植被變化趨勢分析

青藏高原植被的長期變化存在空間差異，并非完全遵循在區域平均尺度上分析得到的規律均勻發展［35］，受到復雜的生物地球物理與生物地球化學過程影響，而在氣候變化的背景下，這些循環過程變得尤為復雜。此外，植被受到氣溫、降水、海拔、經緯度等環境變量的影響，植被分布狀況具有明顯的空間異質性，因此對于植被變化的評估十分困難。但其是探究陸地生態系統對氣候變化響應的關鍵［36］。對于青藏高原這一巨大的地理單元，在探究區域植被的具體變化趨勢時，應考慮廣闊的高原因自然條件的不同而導致的差異性與復雜性。因此，進一步基于年際季節性最大NDVI數據逐像元進行Theil-Sen趨勢估計并結合Mann-Kendall趨勢檢驗，對青藏高原植被隨時間變化趨勢及其相應空間格局進行探究。

在2000—2014年15年間，青藏高原大部分區域的NDVI呈現上升趨勢，通過Theil-Sen趨勢估計確定的植被變化斜率在平均74.04%的范圍呈現正值，反映高原植被逐漸綠化，其中29.29%的像元表現為顯著的綠化趨勢（圖6～9）。但是絕大部分區域NDVI變化趨勢較小，集中在0～0.02（10a）-1。綠化趨勢較大的區域基本分布于高原的東部，具體的分布格局因數據源與研究季節而異。

圖6 青藏高原春季NDVI變化趨勢的空間分布Fig.6 Spatial distribution of trends in NDVI changes in spring on the Tibetan Plateau

不同數據集對探究青藏高原NDVI時間變化趨勢所造成的差異難以忽視。數據的變化趨勢值標準差說明了6種NDVI普遍在青藏高原南部及東南部差異較大（圖10），在春季差異最為明顯，且廣泛分布在高原東南部。在其他季節中，不同數據集的變化趨勢值標準差相對較小，趨勢斜率的差異主要分布在青藏高原南部邊界。與區域平均NDVI長時間序列的變化趨勢結果一致（表2），SPOT NDVI反映的植被綠化趨勢相對于其他數據更加顯著，91.18%的像元在生長季中呈現綠化，29.01%的像元呈現顯著綠化，分別超出MODIS 250 m、500 m、1 km NDVI、MODIS 5.5 km NDVI與GIMMS ND?VI 7.16%、5.53%、4.65%、2.87%與15.89%。顯著綠化趨勢達到0.02（10a）-1以上有23.86%，超出趨勢分布具有較高一致性的MODIS數據11.25%，超出綠化像元占比最少的GIMMS數據16.76%（圖11）。其次，春季SPOT NDVI呈現顯著增長的像元占比最大，達到33.86%，大范圍分布在青藏高原東南部；GIMMS NDVI呈現顯著增長的像元則在春季占比最小，而顯著下降的像元達到了7.88%，占比最大，分布在高原南部邊界。這可能是春季相比其他季節標準差更大的原因。在其他季節中，GIMMS NDVI也具有呈現植被綠化趨勢較小，褐化像元占比明顯高于其他5種數據集的特點。GIMMS NDVI顯著褐化的像元占比平均達到5.83%。MODIS NDVI隨著分辨率的提高，顯著綠化趨勢的像元占比與綠化趨勢卻往往逐漸降低。總體上，逐像元的NDVI長時間變化趨勢與青藏高原區域平均的多年植被變化相一致，GIMMS NDVI表現高原植被綠化范圍與綠化趨勢最小。綠化像元最多且趨勢最大的SPOT NDVI，在區域平均NDVI變化斜率也最大。

圖10 6種NDVI變化趨勢標準差的空間分布Fig.10 Spatial distribution of standard deviation of six trends in NDVI changes：spring（a），summer（b），autumn（c）and growing season（d）

圖11 青藏高原各季節NDVI變化趨勢的頻率分布Fig.11 Frequency distribution of trends in NDVI changes in each season on the Tibetan Plateau

4 討論

MODIS NDVI與SPOT NDVI數據集是目前評估長時間序列植被變化的可靠數據來源，有著廣泛的一致性。在青藏高原區域尺度NDVI時間變化分析與像元尺度NDVI時間變化分析中，MODIS 250 m、MODIS 500 m、MODIS 1 km與SPOT NDVI數據集顯示青藏高原植被整體呈現綠化趨勢，MODIS 5.5 km NDVI也在春、夏兩季呈現顯著的增長趨勢。這與青藏高原正在變得溫暖的認知相符合［1］，氣溫的上升往往被認為有利于植被的生長［6］。氣候變暖主導了青藏高原的綠化趨勢，可以通過促進植物光合作用及延長生長季來增加寒冷地區的植被綠度，對植被變化產生積極的影響［37-38］。但另一方面，在青藏高原的干旱與半干旱區域，高寒草甸與高寒草原對水分條件十分敏感，變暖引起的蒸散發增加可能致使植被受到干旱脅迫［39-40］，從而導致植被退化。冬季的增溫過強可能使得植物休眠延遲，從而推遲春季物候，縮短生長季［7］，導致植被綠度降低。這能夠解釋在青藏高原的部分地區，NDVI呈現下降趨勢的現象。

圖7 青藏高原夏季NDVI變化趨勢的空間分布Fig.7 Spatial distribution of trends in NDVI changes in summer on the Tibetan Plateau

多源NDVI數據集在青藏高原植被的空間分布與長時間序列變化中都存在著難以忽視的差異。盡管MODIS NDVI與SPOT NDVI在青藏高原區域范圍內出現了較為清晰一致的變化模式，但具體的綠化幅度及綠化范圍因分析的數據集而異。在探究2000—2014年的植被變化時，GIMMS呈現褐化趨勢的像元占比顯然多于不同分辨率的MODIS NDVI與SPOT NDVI，區域尺度上的GIMMS NDVI的變化斜率也明顯小于二者。前人研究指出，青藏高原地區GIMMS NDVI要小于SPOT NDVI與MODIS NDVI，而GIMMS NDVI與SPOT NDVI、MODIS NDVI的差異大于SPOT NDVI與MODIS NDVI的差異［19］。SPOT NDVI與GIMMS NDVI數據在反映青藏高原1998—2006年植被變化時，呈現了較大的異質性，SPOT NDVI所反映的綠化范圍與綠化趨勢明顯大于GIMMS NDVI［11］。此外，許多關于青藏高原植被動態的研究也得出了基于GIMMS NDVI所確定的春季物候變化趨勢與基于MODIS NDVI、SPOT NDVI、MODIS增強型植被指數（EVI）或日光誘導葉綠素熒光（SIF）等其他遙感數據集所確定春季物候變化趨勢不一致的結論［1，12，14，41-42］。部分學者認為這種物候變化趨勢的不一致可能與GIMMS NDVI的降低有關［11，41］，并由此推斷在青藏高原區域使用GIMMS NDVI探究植被變化會帶來研究結果的不確定性。

圖8 青藏高原秋季NDVI變化趨勢的空間分布Fig.8 Spatial distribution of trends in NDVI changes in autumn on the Tibetan Plateau

不同數據集對植被變化分析的差異性可能由許多因素導致，例如空間分辨率、時間分辨率、傳感器設計、校正方法、投影系統、地理定位誤差之間的差異［15，17-18，21］。首先，從8 km到250 m的空間分辨率，隨著觀測尺度的變化，尺度依賴性和尺度效應的影響隨之體現［19，43］。其次，衛星平臺及其所搭載的傳感器對NDVI產品有著重要影響。GIMMS NDVI、SPOT NDVI、MODIS NDVI數據集均來自于多傳感器系統［21］，但相比AVHRR，其他二者的傳感器擁有更穩定的軌道與改進植被監測的光譜配置［20］。MODIS傳感器與SPOT-VGT傳感器在紅外與近紅外的波段設計更為接近［19］，這可能也是其獲取的NDVI更為接近的原因之一。來自于AVHRR傳感器的GIMMS NDVI是目前時間序列最長的NDVI產品，但AVHRR系列傳感器最初并非為植被研究設計［18］，缺少機載校準、光譜規格不連續［44］且光譜結構也難以進行更為精準的大氣校正［18］。但大氣成分的散射與吸收是NDVI監測的重要誤差來源。有研究指出GIMMS NDVI與其他衛星遙感數據的差異可能受青藏高原春季氣溶膠濃度增加的影響［11］。

5 結論

研究結果表明，SPOT NDVI所反映的植被綠化呈現顯著且迅速的特征，27.44%的像元呈現顯著綠化趨勢，分別超出MODIS NDVI與GIMMS ND?VI 4.10%、15.89%；在生長季，SPOT顯著綠化趨勢達到0.0182（10a）-1，高于其他顯著變化的數據0.0078～0.0090（10a）-1。MODIS NDVI隨著分辨率的提高，顯著綠化的像元占比與綠化趨勢卻逐漸降低；MODIS數據之間顯著綠化的像元占比相差不足2.80%。GIMMS NDVI顯著褐化的像元占比平均達到5.83%，超出其他數據集3.37%～5.51%。在春季，GIMMS NDVI表現顯著褐化像元占比最大，達到7.88%，與顯著綠化的像元占比相當。GIMMS的區域平均NDVI具有最小的顯著綠化趨勢，為0.0092（10a）-1，并且在春、夏兩季呈現不顯著的褐化趨勢。基于GIMMS NDVI探究青藏高原植被變化特征時，特別是對春季物候研究，可能會造成結果具有較大的不確定性。而SPOT與MO?DIS NDVI體現了較高的一致性，可以互為補充，探究高原植被變化。

圖9 青藏高原生長季NDVI變化趨勢的空間分布Fig.9 Spatial distribution of trends in NDVI changes in growing season on the Tibetan Plateau