青藏高原保護區(qū)植被與生態(tài)系統服務功能變化的一致性分析

程琳琳, 劉 華, 劉焱序

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

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

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

1 數據與方法

1.1 研究區(qū)與數據來源

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

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

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

表1 不同生態(tài)系統服務功能的生態(tài)系統容量評估方法

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年內平均值合成結果中的青海湖區(qū)域呈現集中的EVI顯著下降。值得注意的是，EVI非顯著增加的像元和非顯著減少的像元往往相互交疊，并不呈集中連片分布。因此，并不能有充足的證據將大面積的植被增加或減少直接歸因于某一局部區(qū)域的生態(tài)恢復或破壞過程。將EVI統計至保護區(qū)尺度上，EVI年內平均值合成與EVI年內最大值合成的趨勢提取結果總體上具有一致性。其中，EVI年內平均值合成結果中有3個區(qū)塊出現EVI顯著減少，此外有非顯著減少區(qū)塊37個，非顯著增加區(qū)塊64個。而在EVI年內最大值合成結果中，不存在EVI顯著變化的樣本，非顯著減少和增加的區(qū)塊分別為35個和69個。其中，面積最大的四個區(qū)塊，即羌塘、阿爾金山、三江源、可可西里，全部呈現出EVI的非顯著增加。總體而言，EVI變化趨勢在像元和保護區(qū)2個尺度上均反映出青藏高原植被整體好轉、局部退化的特征。

2.2 保護區(qū)生態(tài)系統服務功能變化

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

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

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

圖2 青藏高原2000-2010年生態(tài)系統服務功能保護區(qū)平均值

3 討 論

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

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

圖3 青藏高原保護區(qū)2000-2010年生態(tài)系統服務功能變化趨勢

4 結 論

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

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