汶川縣生態系統服務對景觀格局及其變化的響應

徐建英,樊斐斐,劉焱序,陳吉星,魏建瑛

1 首都師范大學資源環境與旅游學院, 北京 100048

2 北京師范大學地理科學學部,地表過程與資源生態國家重點實驗室, 北京 100875

生態系統服務是指人類從生態系統中獲得的益處[1]。土地作為陸地生態系統的載體,其類型、格局和利用強度都會影響生態系統服務[2]。首先,土地類型直接決定了生態系統服務的主導類型,如農田生態系統糧食供給服務和自然生態系統各種調節服務[3]；其次,土地利用的變化引起生態過程發生相應的改變,進而導致生態系統服務發生量變或質變[4]；最后,不同的土地利用強度也會造成所提供生態系統服務的差異,這種對應關系往往需要較高精度的土地利用數據來反映[5]。景觀格局是由不同的土地利用類型組成并在空間上呈現不同的分布與組合形式[6- 7]。從生態系統服務供給與維持的角度出發,景觀格局的變化會影響生態系統物質流、能量流和信息流,引起生態過程的改變[8- 9],最終引起生態系統服務的改變。由于土地利用類型變化往往不能完全反映其布局形態的演化規律,景觀格局指標所展現的空間信息值得被特別關注。

目前,生態系統服務與景觀格局集成是區域環境政策和生態系統管理等跨學科領域的熱點問題[10- 11],但此類研究較多地直接以土地利用類型為基礎評估生態系統服務,通過相關和回歸分析探究生態系統服務與景觀格局的關系。例如：Su等分析了杭嘉湖4個生態區的景觀格局和生態系統服務價值變化,發現除了水源涵養和土壤保持外其他生態系統服務都與景觀指標顯著相關[12]；Ayinuer等發現新疆艾比湖流域的13個景觀指標中有8個顯示出與生態系統服務有顯著的相關性[13]；Zhang等分析廣西壯族自治區喀斯特地區生態系統服務價值變化對景觀格局的響應中發現,總面積、最大斑塊指數等指標與生態系統服務呈正相關,分離度和斑塊豐富度則與服務呈負相關[14]；Zhang等認為斑塊大小會造成景觀指標的明顯差異,從而影響景觀指標與生態系統服務間的關系[15]。然而,此類研究多采用土地利用類型賦值的價值當量法評估生態系統服務,進而解析生態系統服務與景觀格局的關系,忽視了在評估生態系統服務中使用土地利用類型,而其與景觀格局指標高度相關,導致相關或回歸分析在輸入數據上存在天然的因果關系。

考慮到區域的差異性和方法的局限性,仍有必要繼續探究重點區域的生態系統服務如何響應于景觀格局的變化。其中,通過生態遙感指標實現生態系統服務制圖,基于獨立于土地利用數據的生態系統服務評估結果,解釋景觀指標與生態系統服務之間的關聯,可以更有效刻畫出景觀格局與生態系統服務間的復雜關系。汶川縣是我國重要自然保護區——臥龍自然保護區的所在地,研究區域內保護區與非保護區并存,自然干擾(地震、滑坡和泥石流等自然災害)和人為干擾(汶川縣作為阿壩州經濟發展較快的區域)明顯,研究該區域內生態系統服務對于景觀格局及其變化的響應有助于了解保護政策的實施效果,自然和人為干擾對區域生態環境的影響,為區域生態保護、景觀規劃和生態系統服務提升提供理論依據。

1 研究區概況與研究方法

1.1 研究區概況

汶川縣位于四川省阿壩藏族羌族自治州的東南部(30°45′—31°43′N,102°51′—103°44′E),四川盆地的西北角,面積約4084 km2。境內地勢由西北向東南逐漸降低,西部海拔較高,其中著名的四姑娘山海拔約6250 m,東南部最低海拔約764 m,河流發育,屬于典型的高山峽谷地形。汶川縣屬于溫帶季風氣候,由于地勢起伏較大,研究區域氣候的垂直地帶性明顯,形成常綠闊葉林、針闊混交等不同的植被帶。該區域植被覆蓋度較高,生態系統碳固定作用明顯；由于地處降雨豐富的高山峽谷地形,高植被覆蓋可有效增強區域的水源涵養和土壤保持服務。縣轄8個鎮,4個鄉,為了保護珍稀動植物種,國家建立臥龍和草坡自然保護區,占全縣總面積的一半以上(如圖1)。

圖1 研究區位置

1.2 數據來源與預處理

本文所采用的土地利用數據(2000、2005、2010、2015年4期)來源于中國科學院資源環境科學數據中心(http://www.resdc.cn/),分辨率為30m,其土地利用類型分為林地、草地、耕地、水域、建設用地和未利用地六大類；DEM和EVI數據源于地理空間數據云(http://www.gscloud.cn/search/)提供的GDEMDEM和MODEV1M產品,分辨率分別為30m和250m；ET和NPP數據是源自美國地質勘察局(USGS)網站提供的MODIS影像MOD16A3產品和MOD17A3H產品,分辨率500 m,由于ET與NPP產品在網站中分別缺少2000年和2015年的數據,在分析中分別2001年作為ET的起始年,2014年作為NPP變化的終止年；月降水數據來源于國家氣象信息中心(http://data.cma.cn/),采用全國站點在Anusplin 4.2中完成基于地形的薄板樣條插值,并裁切得到汶川縣2000—2015年月均降水量；土壤數據源自FAO的HWSD 1.2全球土壤同化數據庫(http://westdc.westgis.ac.cn/)。為了保證數據分辨率的統一性,本文將所有數據的空間分辨率通過三次卷積內插重采樣到30 m分辨率。

2 研究方法

2.1 景觀格局變化

為了明晰2000—2015年汶川縣土地利用變化過程,本研究采用土地類型變化趨勢指數(Ps)為描述指標。該指數的取值在[-1,1]內,若Ps>0,則表示此類土地類型面積在增加,越接近1則表示其他地類主要轉為該地類；Ps<0,則相反,詳見以下公式[16]：

(1)

式中,Ps為某一類土地類型的變化趨勢指數；ΔUin為其他類型轉變為該類型的總面積；ΔUout為研究T期間內某一類土地類型轉變為其他類型的總面積；T為研究時間。

景觀格局指數是反映其結構構成和空間特征的簡單量化指標,分為斑塊水平、類型水平、景觀水平[17]。本文選擇類型和景觀2個水平進行分析,從景觀聚集程度、破碎程度、形狀復雜程度和多樣性等方面描述汶川縣的景觀格局特征,在類型水平選取了最大斑塊指數(Largest Patch Index, LPI)、斑塊形狀指數(Landscape Shape Index, LSI)、斑塊密度(Patch Density, PD)、平均斑塊面積(Mean Patch Area, AREA_MN)、邊緣密度(Edge Density, ED)、聚集度指數(Aggregation Index, AI)6個指標；在景觀水平選取了最大斑塊指數(Largest Patch Index, LPI)、斑塊形狀指數(Landscape Shape Index, LSI)、斑塊密度(Patch Density, PD)、平均斑塊面積(Mean Patch Area, AREA_MN)、邊緣密度(Edge Density, ED)、聚集度指數(Aggregation Index, AI)、香濃多樣性指數(Shannon′s Diversity Index, SHDI)和蔓延度指數(Contagion Index, CONTAG)形成共8個指標。考慮到部分土地類型的覆蓋面積較大,較小的統計尺度會形成大量均質景觀單元,因而將景觀格局指數統計到常見的千米尺度進行分析,包含像元一千余個。利用ArcGIS中Create Fishnet工具得到3 km邊長格網,輸入各年份土地利用類型數據,通過Fragstats 4.2軟件計算得到每個格網內的各項景觀格局指數。計算景觀格局指數使用了林地、草地、耕地、水域、建設用地和未利用地6類土地類型,但由于研究區水域變化不明顯以及建設用地、耕地和未利用地的總面積僅占全域的3.9%,因此本文在敘述類型水平景觀指數時只展示了林地和草地的景觀格局指數。

2.2 生態系統服務評估

考慮到生態系統調節服務是自然生態系統最為重要的生態系統服務類型,對生態系統調節服務的空間分析一直是學界研究的熱點議題[18]。尤其是汶川縣大面積處于自然保護區中,生態系統供給服務和文化服務相對而言不突出,因而生態系統調節服務對景觀格局及其變化的響應成為本研究的關注重點。

基于生態遙感指標,本研究評估了最為常見的三種生態系統調節服務：碳固定、水源涵養、土壤保持[19]。碳固定采用MOD17A3H產品中的各年NPP取值予以指代,單位為g C/m2；水源涵養主要基于地形和土壤指標對產水量進行修正[20- 21],土壤保持主要依循通用水土流失方程(The Revised Universal Soil Loss Equation, RUSLE)進行計算[22],具體過程如下。

2.2.1水源涵養

在水源涵養量計算中,首先通過水量平衡原理得到汶川縣的產水量,然后根據土壤的厚度、滲透性、地形和流速系數等因素計算得到涵養量,具體的計算公式如下[20]：

Y=P-ET

(2)

式中,Y、P和ET分別為產水量、降水量和蒸發量。

WR=min(1,249/V)×min(1,0.9×D/3)×min(1,Ksoil/300)×Y

(3)

式中,WR為多年平均水源涵養量(mm)；V為流速系數(采用模型參數表數據[21])；D為地形指數,無量綱；Ksoil為土壤飽和導水率(cm/d)。

(4)

式中,Watershed pixel count(Wpc)為集水區柵格數,無量綱；Soil depth(Sd)為土壤深度(mm)；Percent slope(Ps)為百分比坡度。

2.2.2土壤保持

采用RUSLE方程計算土壤侵蝕量,考慮到大尺度人類工程措施無法表征,默認土壤保持措施因子為1,方程的公式如下[22]：

A=R?K·L·S·C·P

(5)

式中,A為單位時間單位面積的平均土壤流失量(t hm-2a-1)；R為降雨侵蝕力因子(MJ mm hm-2h-1a-1)；K為土壤可蝕性因子(t h MJ-1mm-1))；L為坡長因子；S為坡度因子；C為作物覆蓋與管理因子；P 為土壤保持措施因子。

R通過采用Wischmeier提出的相關公式[23]計算得出,如下：

(6)

式中,R為降水可持續性因子,P為年降水量(mm),Pi為第i月的降水量(mm)。

K采用EPIC模型[24]得到,計算公式如下：

(7)

SNI=1-SAN/100

(8)

式中,K為土壤可侵蝕性；SAN、SIL、CLA為砂粒、粉粒、粘粒含量；SNI為非砂粒所占百分比；C為土壤有機含量(%)。

坡度因子S根據坡度θ[25]可計算得到：

(9)

坡長L根據以下公式計算得到[26- 28]：

L=(λ/22.1)m

(10)

m=β/(1+β)

(11)

β=(sinθ/0.0896)/(3.0(sinθ)0.8+0.56)

(12)

式中,λ、m、β分別為水平坡長、坡長指數、細溝侵蝕。

考慮到EVI比常用的NDVI不易飽和,更適用于植被密集地區,因而根據EVI換算植被覆蓋度c,從而得到植被因子C[29]：

(13)

其中EVI為12月的平均值。

土壤保持量應為潛在土壤侵蝕量與現實侵蝕量的差值,其中潛在土壤侵蝕量應排除地表覆被類型,即C=1：

Ar=Am-A

(14)

Ar為土壤保持量,Am為潛在土壤侵蝕量,A為現實土壤侵蝕量。

2.3 Mann-Kendall (M-K)檢驗

利用M-K顯著性檢驗分析生態系統服務在連續時間內的變化趨勢,該方法常用于連續時間序列的趨勢分析[30- 31],同時采用Theil-Sen (TS)過程輔助,其計算過程如下[32- 33]：

(15)

式中,Slpoe > 0表示增加,Slpoe<0表示減小；n是時間序列長度,xi和xj分別是時間i和j的生態系統服務量；σ是標準偏差,| Z | > 1.64則認為變化趨勢顯著。

3 結果與分析

3.1 景觀格局

汶川縣土地類型主要是在林地、草地和耕地之間轉換其中在2000—2015年間,林地增加了112.21 km2、草地減少了139.58 km2、耕地增加了20.34 km2, 而水域、建設用地和未利用地的土地變化相對較小,分別為0、3.23 km2和0。為了強調研究區景觀格局的變化效應,考慮到研究區的水域、建設用地和未利用地的占地面積以及變化面積較小,因此在后續分析中主要關注耕地、林地和草地變化對生態系統服務的影響。具體而言,汶川縣2000—2015年內耕地和林地面積處于緩慢增長的趨勢,草地面積變化趨勢則相反。其中,不同時間序列變化趨勢指數Ps(表1)表明,2000—2005年間耕地和林地的變化趨勢指數Ps均小于0,即其面積處于減少趨勢,而草地的變化趨勢指數Ps大于0,面積處于增長趨勢；2005—2010年耕地、林地的變化趨勢指數均大于0,面積處于增長趨勢,尤其是耕地增長量較大；草地相反,Ps低至-0.42；2010—2015年林地、草地和耕地的變化趨勢和動態度指數Ps波動較小(-0.01、0.02、0.11),變化均不明顯。

表1 不同土地類型2000—2015年的變化趨勢指數(Ps)

2000—2015年景觀指數的變化趨勢見圖2,分析結果表明研究時段內汶川縣整體景觀更加破碎化,斑塊形狀復雜化,斑塊異質性增加。在景觀水平,整體景觀的AREA_MN從452.57減少到424.37、AI從97.91減少到97.79,而LSI從23.82增加到25.13、LPI從50.84增加到52.76、CONTAG從71.35增加到73.67、SHDI從0.81增加到0.83。草地的景觀變化也顯示出與整體景觀較為一致趨勢,其AREA_MN 、AI和LPI均表現出下降趨勢,PD、ED和LSI有所增加,表明景觀破碎度增加,聚集程度減弱,最大斑塊面積在逐漸下降,形狀趨于復雜。但是,林地景觀的變化規律則與上述兩者不同,林地的AREA_MN從2000年的1116.61增加到2015年的1870.84,LPI從50.81增加到52.76、AI從98.21增加到98.33,但PD從0.05下降到0.03,說明林地破碎化減弱,團聚度增強。

圖2 不同時期景觀指數的變化趨勢

3.2 生態系統服務變化

由圖3中可知,研究區2000—2015年絕大部分地區水源涵養和土壤保持量呈增長趨勢,而NPP呈下降趨勢。其中,2000—2015年內汶川縣水源涵養量波動上升,除臥龍自然保護區的西部外,大部分地區水源涵養呈增長趨勢；另外,草坡和臥龍保護區的水源涵養量主要在300mm左右波動上升,而非保護區的水源涵養在250mm處波動上升,即保護區的水源涵養要明顯高于非保護區,表明保護區的生物多樣性保護起到一定的積極作用(圖3)。汶川縣土壤保持量在15年間波動較大,出現多處低谷(2000、2004、2007、2015年)和峰值(2001、2006、2008、2013年)；大部分地區土壤保持處于增長趨勢,但在城鎮比較聚集的地區土壤保持量明顯減少(圖3)。碳固定服務在2000—2015年呈減少趨勢,研究發現三個區域中NPP變化趨勢一致,呈現同升同降的規律,可推測碳固定服務的變化受更宏觀的氣候條件影響較大(圖3)。三個區域的生態系統服務比較分析表明,臥龍和草坡保護區內的水源涵養和土壤保持有明顯的提高,但在臥龍保護區與縣域的交界處,水源涵養和碳固定服務呈明顯的下降趨勢。

圖3 2000—2015年生態系統服務時空變化趨勢

3.3 生態系統服務對景觀格局的響應

將2015年生態系統服務評估結果以分位數法等分為低、較低、中等、較高、高5個等級,根據每個等級上景觀指數的波動范圍,得到不同生態系統服務等級所對應的景觀指數(圖4和圖5)。

圖4 不同生態系統服務對景觀水平景觀指數的響應

3.3.1景觀水平

景觀水平上生態系統服務主要與LSI、SHDI、AI、 LPI等指數密切相關(圖4)。以景觀格局指數的中位數為對象,當LSI分別在小于1.95、小于1.88時,所對應的土壤保持、水源涵養量處于較高水平,當LSI范圍小于1.87時,碳固定服務處于較低水平；當SHDI小于0.55時,水源涵養和土壤保持均處于較高水平,而碳固定則反之；AI在98.36以上時,土壤保持和水源涵養量相對較高,碳固定值較低；LPI大于85時,土壤保持和水源涵養量較高,而碳固定則反之。但三種生態系統服務對CONTAG并不敏感,CONTAG值一直在60左右波動。由此可知,當LSI和SHDI較低、AI和LPI較高時,土壤保持和水源涵養量較高,而碳固定相對較低。此外,箱圖中顯示的奇異值被界定為大于上四分位數1.5倍四分位數差的值,或者小于下四分位數1.5倍四分位數差的值,結果中奇異值較多存在與LSI、AI、LPI、SHDI等景觀指數,體現了格網中的這些指數變化幅度較大。說明汶川縣景觀格局具有較強的異質性,而這種異質性特征很可能導致了生態系統服務空間分布的異質性。

3.3.2類型水平

以景觀指數的中位數作為參考,生態系統服務對類型水平景觀格局的響應與景觀水平不同。由圖5可知,對于土壤保持,林地和草地斑塊形狀較為簡單時,會使得景觀的土壤保持能力提高,主要表現為林地和草地的LSI對土壤保持服務的響應與景觀水平較為一致, LSI分別小于1.90和1.79時,林地和草地的土壤保持量處于較高水平(景觀水平為1.88)。增大草地的最大斑塊面積是提高景觀土壤保持能力的有效途徑,主要體現在LPI分別在17.83和81.85時,林地和草地土壤保持達到較高水平；說明草地LPI對土壤保持服務的響應與景觀水平(82.40)較為一致,而林地則不一致。草地斑塊破碎度較低時,土壤保持量較高,而林地土壤保持對PD的響應并不明顯,表現為PD在0.2上下波動。

圖5 不同生態系統服務對類型水平景觀指數的響應

水源涵養對景觀指數的響應與土壤保持大多一致,其明顯區別在于,增大林地的最大斑塊面積可有效提高景觀水源涵養能力的有效因素。其主要表現在林地與草地LPI在81.13和12.14時,水源涵養處于較高水平,而林地LPI與景觀水平(82.40)范圍大致相當。對于碳固定服務,林地和草地斑塊形狀較為復雜時的碳固定相對較高。主要表現在林地與草地LSI分別低于1.89和1.42時,碳固定服務處于低水平,與景觀水平LSI(1.47)對碳固定的響應大致相當。增大草地的最大斑塊面積是提高景觀碳固定能力的有效方式,主要表現在林地和草地LPI分別在18.19、90.04時,碳固定處于低水平,其對草地LPI的響應與景觀水平(89.97)一致。

3.4 生態系統服務變化對景觀格局變化的響應

2000—2015年生態系統服務與景觀格局變化量的對應關系表明,只有碳固定變化明顯響應于景觀格局變化,景觀變化會引起研究區大面積范圍碳固定值的降低,而水源涵養和土壤保持在研究區中并沒有明顯響應于景觀格局的變化。不同地類的景觀指數變化對碳固定服務變化的敏感性存在差異(圖6)。以碳固定服務顯著變化像元所占格網面積百分比的中位數為對象,在景觀水平上,當LSI、LPI和AI發生明顯變化時,會造成約50%以上的碳固定服務顯著減少,而在LPI、LSI和PD穩定區間內該面積百分比中位數僅為29.98%、33.25%和32.36%。

圖6 碳固定顯著減少比例對景觀指數變化的響應(2000—2015年)

在類型水平上,生態系統服務響應景觀格局現象更具有針對性。例如,林地最大斑塊面積顯著減少、斑塊形狀趨向簡化、聚集度減弱時,引起碳固定在大面積范圍內顯著減少。主要體現在LPI、PD在減少區間內,對應的碳固定顯著減少所占面積百分比分均在約55%以上,高于其他區間。草地的最大斑塊面積增加、形狀趨于復雜、聚集程度增加時會造成碳固定在較大范圍內顯著減少,表現在LPI、LSI、PD在增加區間內,碳固定顯著減少所占面積百分比中位數均在60%以上,而其穩定區間的百分比中位數低于40%。

4 討論

本研究結果表明,并不是所有景觀格局變化均會引起生態系統服務變化。以顯著變化的面積占到格網的百分比的中位數為對象,土壤保持服務的變化量對景觀格局的變化呈現出震蕩波動的狀態(圖7)。例如,景觀水平LSI、LPI、PD的不同變化區間內,土壤保持顯著增加所占面積百分比均在15%上下波動。在類型水平也是如此,如草地也是在LSI、LPI、PD的不同變化區間內,面積百分比也在15%上下波動。造成這種現象的原因可能是：雖然景觀格局的改變影響土壤流失方程中的C值,但在30m分辨率植被覆蓋度遙感監測中這種變化不易被識別,尤其是對高植被覆蓋的保護區而言林地和草地均呈現較高的植被覆蓋度,導致遙感評估得到的土壤保持服務變化往往不依賴于景觀格局變化。

圖7 土壤保持顯著增加比例對景觀指數變化的響應(2000—2015年)

不同尺度也會影響生態系統服務變化對景觀格局變化的響應程度。在全縣范圍內,水源涵養變化與景觀格局變化并沒有明顯的響應,但是,在局部區域景觀劇烈變化時,水源涵養變化隨景觀格局變化的響應依然有明顯響應。如在臥龍保護區的東南部選擇34個格網單元的景觀指數如LPI、LSI、PD等景觀指數變化相對明顯：中位數分別為18.29、-0.61、-0.35,其水源涵養量顯著增加面積中位數達56.34%；而臥龍保護區的東北部25個格網單元的景觀格局變化較小：LPI、LSI、PD中位數分別為-0.89、0.07、0.11,對應的水源涵養顯著增加面積中位數僅為1.41%。

就研究方法而言,本文揭示了生態系統服務與景觀格局關系的測度直接依賴于生態系統服務評估手段,高精度的生態過程地面觀測數據補充將有助于進一步揭示生態系統服務與景觀格局之間的復雜聯系。此外,臥龍、草坡保護區與其他地區的交界處往往是生態系統服務易于變化的區域,因此,保護區與非保護區生態系統服務供應存在明顯差異。草地對研究區生態系統服務的影響較大,優化區域區內合理的林地、草地規模配比與空間布局將更有效的提升研究區的生態系統服務水平。

5 結論

針對重點區域的生態系統服務如何響應于景觀格局變化的科學問題,本文識別了汶川縣土壤保持、水源涵養和碳固定對景觀格局及其變化的響應,得到以下三點發現：

2000—2015年汶川縣林地、草地和耕地之間存在轉變,整體景觀更加破碎化、斑塊形狀復雜化,斑塊異質性增加,整體景觀格局的變化受草地影響較大。研究區碳固定服務有所減少,而大部分區域水源涵養和土壤保持呈增長趨勢,在臥龍和草坡保護區內的水源涵養和土壤保持有明顯提高,但在臥龍保護區與縣域交界處的生態系統服務呈下降趨勢。

不同生態系統服務對景觀格局的響應存在差異,其不僅表現在指標上,也體現程度上。在景觀水平,土壤保持、水源涵養和碳固定服務主要對LSI、SHDI、LPI等指數較為敏感,而對CONTAG響應并不明顯；在類型水平,生態系統服務對類型水平景觀格局的響應與景觀水平有所不同,林地類型水平與景觀水平差異較大,而草地類型水平與景觀水平大致相當。例如土壤保持高值區間對應林地和草地LPI中位數分別為17.83和81.85,而景觀水平LPI中位數82.40。

汶川縣生態系統服務與景觀格局變化量的對應關系發現,只有碳固定變化明顯響應于景觀格局變化,景觀的顯著變化會引起研究區大面積范圍碳固定值的顯著降低,而水源涵養和土壤保持在研究區中并沒有觀測到這一響應方式。在景觀水平,景觀指數變化直接影響碳固定服務變化；在類型水平上,林地LPI和PD的減少和草地LPI、LSI和PD的增加,會導致碳固定服務的顯著減少。