青藏高原保護區植被與生態系統服務功能變化的一致性分析

程琳琳, 劉 華, 劉焱序

〔1.中國礦業大學(北京) 地球科學與測繪工程學院, 北京 100083;2.北京師范大學 地理科學學部 地表過程與資源生態國家重點實驗室, 北京100875〕

植被是陸地生態系統的主體，是調節全球和區域物質和能量平衡的關鍵要素，也是生態系統提供服務的重要載體[1-2]。在近年研究中，遙感植被指數已被廣泛用于陸地植被時空變化的刻畫，其中歸一化差值植被指數(NDVI， the normalized difference vegetation index)是衡量植被變化的最常用遙感監測參數[3-6]。不同數據源的NDVI產品交互驗證表明，近30 a來中國植被變化普遍呈現增加趨勢，但不同數據產品和趨勢提取方法所呈現的區域植被變化狀態略有差異[7-9]。然而，NDVI在高覆蓋植被區的取值往往存在飽和現象，導致不利于監測高密度植被的時空演化。為解決這一問題，增強型植被指數(EVI， enhanced vegetation index)被逐步應用于植被變化研究中，并已成為監測大尺度植被變化特征的重要工具[10-12]。

青藏高原被譽為亞洲水塔，是世界多條大江大河的發源地。青藏高原生態系統的碳固定、土壤保持、水源涵養、棲息地保護等生態系統服務功能對中國和亞洲人類生活具有重要的影響。作為評估生態系統服務功能的關鍵指標，青藏高原植被變化一直備受關注[13-15]。1982—2011年，青藏高原生態系統碳匯的增加占全國增加碳匯的10%左右，但也有部分地區出現了植被退化[16]。對應了生態十年評估，是前人研究往往采用的時段，并且2000年以來，中分辨率成像光譜儀(MODIS， moderate-resolution imaging spectroradiometer)系列產品的發布為更高精度的刻畫植被變化、評估生態系統服務功能提供了新的數據支持。在此背景下，本研究聚焦青藏高原各類保護區，試圖回答以下兩項科學問題： ①基于MODIS EVI產品的植被變化制圖是否反映了青藏高原植被整體好轉局部退化的特征； ②青藏高原保護區內2000—2010年的生態系統服務變化是否與同期的植被變化一致。研究結果將有助于進一步明晰生態系統服務功能和遙感植被指數的關系，為區域生態系統服務功能評估與生態保護提供科學參考。

1 數據與方法

1.1 研究區與數據來源

本研究所采用的青藏高原邊界源自張鐿鋰等的劃定結果，源自全球變化科學研究數據出版系統[17-18]。青藏高原的保護區邊界源自世界保護區數據庫(WDPA, world database on protected areas)(https:∥www.protectedplanet.net/)?？紤]到數據庫中保護區矢量圖層相互交疊，以取最大面積為規則，對圖層拓撲關系進行調整，共提取保護區矢量區塊105個。2000—2010年生態系統服務功能制圖數據源自Zhang等的評估結果，包括碳固定、土壤保持、水源涵養、棲息地保護以及總生態系統服務功能5項[19]。所有生態系統服務功能圖層標準化至0～1，并統計各保護區的平均值。

MODIS EVI數據源自地理空間數據云中國月合成產品(http:∥www.gscloud.cn)。在月合成影像向年值轉化中，最大值合成法較為常見[20-23]。同時，也有研究采用生長季的平均NDVI作為各個季節的植被變化[24-26]。本研究擬同時采用以上兩種方法，進行相互對比。由于青藏高原東南與西北的生長季時段差異較大，本研究將EVI小于0.1的部分作為非植被區域予以去除，而將其余部分作為有植被覆蓋的生長季EVI，從而得到基于EVI年最大值和年平均值。

1.2 變化趨勢識別

為避免最小二乘法方法對異常點擬合而產生的參數誤差，Theil-Sen(TS)中位數斜率與Mann-Kendall(M-K)單調性檢驗作為一種非參數趨勢識別方法逐漸被引入植被變化分析中[8,27-28]。TS-MK算法對植被變化的具體形式可表征為[29-30]：

(2)

(3)

(4)

(5)

式中：slope>0表示植被增加; slope<0表示植被減少；ti，j——年份; EVIi，j——第i和j年的觀測像元EVI值； var(S)表示方差；n——時間長度;Z的大小表示趨勢的顯著性，基于S指數判定，若|Z|>1.64則趨勢是顯著的。

1.3 生態系統服務功能評估

當前，生態系統服務功能的評估方式尚未統一，一些評估模型在運算中采用土地利用參數賦值的方式，掩蓋了同種土地利用類型中不同植被覆蓋條件的生態系統服務功能差異。Zhang等改進了一種針對生態系統服務供給的生態系統容量評估方法，其以CASA(carnegie-ames-stanford approach)模型評估逐月凈初級生產力(net primary productivity, NPP)，并以年NPP為基礎指標，將不同生態系統服務功能表示為NPP及其對應的環境因子加權(表1)。由于NPP與植被覆蓋程度直接關聯，這種生態系統評估方式可以更有效的體現出年際變化，更有助于觀測生態系統服務功能對植被變化的響應情況，并服務于生態系統動態評估。由于荒漠、水體難以估算其NPP，本研究生態系統服務估算主要以森林、草地、耕地等有植被覆蓋的生態系統為對象。

表1 不同生態系統服務功能的生態系統容量評估方法

2 研究結果

2.1 青藏高原EVI變化

對比2000—2010年青藏高原EVI變化如圖1所示。在柵格尺度上，EVI年內平均值合成與EVI年內最大值合成的趨勢提取結果具有差異性。具體而言，EVI年內平均值合成結果中，7.44%的像元顯著增加，63.49%的像元呈現非顯著增加趨勢；而在EVI年內最大值合成結果中，顯著增加像元僅占3.38%，非顯著增加像元也降至54.14%。相比之下，EVI年內平均值合成結果中27.98%的像元呈現非顯著減少，而在EVI年內最大值合成結果中該比例提升至40.83%。同時，兩組提取結果中呈現顯著減少的像元極少，僅有EVI年內平均值合成結果中的青海湖區域呈現集中的EVI顯著下降。值得注意的是，EVI非顯著增加的像元和非顯著減少的像元往往相互交疊，并不呈集中連片分布。因此，并不能有充足的證據將大面積的植被增加或減少直接歸因于某一局部區域的生態恢復或破壞過程。將EVI統計至保護區尺度上，EVI年內平均值合成與EVI年內最大值合成的趨勢提取結果總體上具有一致性。其中，EVI年內平均值合成結果中有3個區塊出現EVI顯著減少，此外有非顯著減少區塊37個，非顯著增加區塊64個。而在EVI年內最大值合成結果中，不存在EVI顯著變化的樣本，非顯著減少和增加的區塊分別為35個和69個。其中，面積最大的四個區塊，即羌塘、阿爾金山、三江源、可可西里，全部呈現出EVI的非顯著增加?？傮w而言，EVI變化趨勢在像元和保護區2個尺度上均反映出青藏高原植被整體好轉、局部退化的特征。

2.2 保護區生態系統服務功能變化

圖2為青藏高原2000—2010年生態系統服務功能保護區平均值。如圖2所示，青藏高原生態系統服務功能的保護區平均值體現為碳固定最高，土壤保持次之，棲息地保護最低。這是由于青藏高原大部分保護區如羌塘、阿爾金山等地盡管面積較大，卻大部分屬于荒漠生態系統，不能提供高質量的棲息地。水源涵養相對低也是由于荒漠生態系統廣布所致。就時間變化而言，水源涵養和棲息地保護服務略有增加，而碳固定和土壤保持服務功能的變化不明顯。

雖然遙感植被指數是估算生態系統服務功能的重要參數，但后者更多的考慮了區域實際的氣候狀況，二者不能等同。圖3所反映的青藏高原保護區生態系統服務功能變化表明，植被增加并不一定意味著生態系統服務功能增加。具體而言，在保護區尺度上，青藏高原生態系統服務功能減少的區塊數量大于服務功能增加的區塊數量，與植被整體好轉、局部退化的總體現象相反。總生態系統服務功能呈現顯著減少的區塊共9個，非顯著減少的區塊共66個，而非顯著增加的區塊僅29個。面積較大的生態系統服務功能顯著減少區包括色林錯和雅魯藏布江中游河谷黑頸鶴保護區等。就各類生態系統服務功能而言，碳固定服務功能與土壤保持服務功能的變化規律比較相似，顯著減少區塊均為7個，非顯著減少區塊分別為61和60個，而非顯著增加區塊分別為36和37個。水源涵養服務功能的顯著減少區塊高達12個，非顯著減少區塊78個，而非顯著增加區塊僅為14個；與之對應，棲息地保護服務功能的顯著減少區塊僅為4個，非顯著減少區塊72個，而非顯著增加區塊為28個。同時，羌塘、阿爾金山、三江源、可可西里4個面積最大的區塊均呈現出棲息地保護服務功能的增加，該服務功能與EVI的對應關系相比其他服務功能更為緊密，但在東南部其他區塊仍顯現出與EVI相反的變化趨勢特征。該結果證明，2000—2010年青藏高原保護區的各項生態系統服務功能并未隨植被增加有明顯的增加趨勢。

圖1 青藏高原2000-2010年EVI變化

圖2 青藏高原2000-2010年生態系統服務功能保護區平均值

3 討 論

目前，遙感植被指數已經成為表述植被變化的常用手段，雖然不同的遙感產品、不同的植被指數、不同的年內合成方法勢必影響最終結果的表征，但全球和中國近年來的植被總體增加已經在學界達成共識[21,8,28]。本研究也證實，盡管EVI與NDVI在取值上確有差異，不同的年內合成方法也導致了一定幅度的結果改變，但不影響青藏高原植被總體增加的結論。然而，與以往認知有所差異的是，本研究認為植被增加并不一定導致除水生態系統以外的陸地生態系統服務的增加。該現象一方面可能源于生態系統服務自身評估方法的不確定性，而另一方面并更有可能是由于生態系統服務功能評估需要考慮當前的氣候狀況，未經過模型處理的遙感參數信息并不應直接作為生態過程、功能的表征途徑。

近年來，InVEST(integrated valuation of ecosystem services and tradeoffs)等生態系統服務功能評估模型在生態系統服務評估中被廣泛采用。然而，此類模型的部分模塊如碳固定、生境質量等對土地利用數據的依賴較強，并不能細致刻畫植被的疏密特征，以及氣候的年際變化。本研究采用以NPP為基礎數據的生態系統服務供給的生態系統容量評估方法，雖然在權重系數的刻畫上較為簡單，但可以有效反映氣候變化下幾項關鍵生態系統服務功能的時空變化格局。尤其是NPP本身估算采用CASA模型，該模型已相對成熟并在近年研究中取得了大量應用。NPP所反映的碳固定生態系統服務功能呈現非顯著退化，即反映了碳固定能力的變化并不能等同于植被變化，這也是符合植被年內季節生長的碳循環規律的。在未來研究中，可以考慮將生態模型的計算結果與遙感參數觀測結果進行進一步對比，從而深入驗證遙感植被指數的實際應用能力。

圖3 青藏高原保護區2000-2010年生態系統服務功能變化趨勢

4 結 論

(1) 青藏高原EVI在2000—2010年以非顯著變化為主，在柵格和保護區兩個尺度上均表現為非顯著增加樣本多于非顯著減少樣本，即青藏高原植被整體好轉局部退化。同時，青藏高原EVI變化并沒有呈現明顯的集中連片特征，推測氣候對青藏高原植被變化的影響大于人類活動對其的直接影響。

(2) 青藏高原保護區各項生態系統服務功能在2000—2010年也以非顯著變化為主，105個區塊中生態系統服務功能顯著減少樣本在4～12個之間，生態系統服務功能非顯著減少樣本明顯多于非顯著增加樣本，與植被變化趨勢不一致。推測年際氣候的差異可能是導致該不一致性的主要因素。