白馬雪山國家級自然保護區景觀格局與生態系統服務關系研究*

陳星宇 趙筱青 楊永貴 普軍偉 王 茜 李 琴 趙祖軍

(1.云南大學地球科學學院，云南 昆明 650500；2.云南大學國際河流與生態安全研究院，云南 昆明 650500；3.云南省生態環境監測中心，云南 昆明650034)

頻繁的人為活動導致自然保護區有效保護范圍縮小、土壤退化、景觀破裂程度增加、景觀多樣性改變等問題[1-2]，同時引起生態系統結構和功能的變化，從而導致生態系統服務功能變化。為促進自然保護區持續發展，給保護區內動植物提供良好的生態環境，有必要對保護區內景觀格局與生態系統服務關系進行研究。許多學者對景觀格局和生態系統服務關系進行了探討。BENNETT等[3]定量研究了生態系統服務受景觀異質性的影響程度。張明陽等[4]、童晨等[5]的研究均發現，生態系統服務與一些景觀格局指數密切相關。但是對于景觀格局與生態系統服務的空間分異規律及兩者的空間耦合關系還未進行深入探討。本研究擬利用GeoDA軟件來探討這一問題。

白馬雪山國家級自然保護區(以下簡稱白馬雪山保護區)是全球24個重點生物地理區之一，也是世界上生物多樣性保護熱點地區之一[6]。近年來，白馬雪山保護區也出現了生物資源和自然環境管理不完善、農民過度放牧和薪材消耗量大、保護區周邊基礎設施建設大增、旅游業快速發展等問題，導致整體景觀格局發生變化，生態系統服務受到不同程度的影響。因此，本研究定量計算了白馬雪山保護區的景觀格局指數和生態系統服務功能，并基于空間自相關分析方法探討景觀格局與生態系統服務的空間關系，以期為白馬雪山保護區生態環境保護以及區域生態紅線劃定提供參考依據。

1 數據與方法

1.1 數據來源

Sentinel-Ⅱ(2019)遙感影像數據來源于美國地質勘探局網站(http://glovis.usgs.com/)；歸一化植被指數(NDVI)遙感影像數據來源于2019年美國國家航空航天局(https://earthdata.nasa.gov/)的MODIS13Q1產品，空間分辨率為250 m，時間分辨率為16 d；數字高程模型(DEM)數據，來源于地理空間數據云(http://www.gscloud.cn/)；氣象數據(包括氣溫、降水量、日照時數、太陽輻射等)來源于中國氣象數據網(http://data.cma.cn/)；土壤數據(包括土壤深度，根系深度和土壤砂粒、粉粒、黏粒及有機質質量分數等)來源于世界土壤數據庫(HWSD)的中國土壤數據集(v1.1)，空間分辨率為1 km。所有空間數據均統一轉換為WGS_1984_UTM_47N坐標系。

1.2 遙感數據處理

1.2.1 土地利用類型數據提取

用大氣校正、波段融合及影像裁剪等技術對Sentinel-Ⅱ(2019)遙感影像數據進行預處理，在野外考察并掌握遙感影像總體特征的基礎上，建立解譯標志，通過人機交互解譯的方式提取白馬雪山保護區2019年的土地利用類型數據。參考《土地利用現狀分類標準》(GB/T 21010—2017)，并結合白馬雪山保護區實際情況，將土地利用類型劃分為林地、草地、耕地、水域及水利設施用地、住宅用地、交通運輸用地、其他用地7類一級用地，并進一步劃分為喬木林地、灌木林地、疏林地、中覆蓋草地、高覆蓋草地、河流水面、湖泊水面、冰川及永久性積雪、內陸灘涂、農村宅基地、公路用地、裸土地、裸巖石礫地、水田、旱地15類二級用地。然后，根據野外調查時間應與遙感數據同期或相近原則，選取各土地利用類型中具有代表性的地塊再次進行野外實地驗證，進一步完善遙感解譯結果。最終解譯結果的Kappa系數大于0.85，解譯精度為88.18%。

1.2.2 NDVI的月和年數據合成

NDVI遙感影像數據通過最大合成法消除云、大氣和太陽高度角等干擾，獲取NDVI時間序列數據集。利用MRT軟件對NDVI遙感影像數據進行投影批處理和平滑處理，每月兩個16 d的NDVI遙感影像數據用最大合成法取最大值作為該月的數據，然后用最大合成法取每月的最大值作為該年的數據。

1.3 景觀格局指數

斑塊數、邊緣密度、蔓延度和香農多樣性4個景觀格局指數是研究區域景觀破裂程度的常用指標[7-8]，每個柵格(30 m×30 m)的各景觀格局指數的計算參見文獻[9]。

通過Fragstats4.2軟件移動窗口法來分析景觀格局指數的空間分布特征。使用移動窗口法時，窗口半徑的選擇十分重要[10]。比較了300、500、1 000 m為窗口半徑時的景觀格局指數穩定性，發現當窗口半徑為500 m時，景觀格局指數的空間變異特征最穩定，可以有效反映景觀格局的空間變異特征，因此選擇500 m為窗口半徑。

1.4 生態系統服務功能

(1) 水源供給服務功能計算

(1)

式中：ETn為n月參考蒸散發量，mm；Tn為n月氣溫，℃。

AWC=SR×(54.509-0.132×SA-0.003×SA2-0.055×SI-0.060×SI2-0.738×CL+0.007×CL2-2.688×OC+0.501×OC2)

(2)

式中：SR為土壤深度或根系深度的較小值，mm；SA、SI、CL、OC分別為土壤砂粒、粉粒、黏粒和有機質質量分數，%。

(2) 土壤保持服務功能計算

每個柵格的土壤保持服務功能通過InVEST模型的SDR模塊計算定量。導入該模塊的相關數據有降水侵蝕因子(R，通過式(3)計算得到)、土壤可蝕性因子(K，通過式(4)計算得到)，植被覆蓋、作物管理和土壤保持措施因子參考《InVEST模型用戶手冊》設置，地形因子從DEM數據中提取。

R=0.444 88·p0.969 82

(3)

式中：p為降水量，mm。

(4)

(3) 植被凈初級生產力(NPP)服務功能計算

每個柵格的NPP服務功能是基于改進的CASA模型形成的植被凈初級生產力(NPP)插件在ENVI軟件中計算定量的。需要導入的相關數據有降水量、太陽輻量、氣溫、土地利用類型、NDVI等數據。通過上述數據導入模型求出月NPP數據，并將一年12個月數據進行累加獲取2019年NPP服務功能。

(4) 生境支持服務功能計算

每個柵格的生境支持服務功能通過InVEST模型的Habitat Quality模塊計算定量。導入該模塊的相關數據有土地利用類型數據、威脅因子數據(指土地利用類型數據中的農村宅基地、公路用地、水田、旱地數據)，對威脅因子的影響距離和土地利用類型對各威脅因子的敏感度參考《InVEST模型用戶手冊》設置。

1.5 雙變量空間自相關分析

空間自相關分析能看出兩個變量是否具有聚集性(包括全局空間自相關性和局部空間自相關性)[11]。本研究運用GeoDA軟件來分析景觀格局指數和生態系統服務功能的空間自相關性[12]。

2 結果與分析

2.1 景觀格局指數空間分布特征

白馬雪山保護區的景觀格局指數空間分布見圖1。

圖1 景觀格局指數空間分布

(1) 斑塊數

白馬雪山保護區景觀的破裂程度較高，斑塊數總數達到8 684個，其空間分布總體呈北高南低的趨勢，說明白馬雪山保護區的北部地區景觀破裂程度更高。

(2) 邊緣密度

白馬雪山保護區景觀的邊緣密度平均值為31.37 m/hm2，其空間分布也呈北高南低的趨勢。

(3) 蔓延度

白馬雪山保護區景觀的蔓延度平均值為72.42%，其空間分布除呈南高北低的趨勢外還表現出四周較高中間較低的趨勢，同樣說明北部地區景觀破裂程度更高，同時邊緣地區連通性較好，而中間地區連通性較差，說明中間存在更多的小斑塊，破碎化程度更高。

(4) 香農多樣性

白馬雪山保護區景觀的香農多樣性平均值為1.26，其空間分布呈北高南低的趨勢，高值區在北部及邊緣地區，景觀元素較多，景觀多樣性較大；低值區在南部地區，景觀元素較少，景觀多樣性較小。這種特殊情況可能是斑塊數造成的。

2.2 生態系統服務功能空間分布特征

白馬雪山保護區的生態系統服務功能空間分布見圖2。

圖2 生態系統服務功能空間分布

(1) 水源供給服務功能

白馬雪山保護區水源供給服務功能為124.94～1 094.56 mm/(hm2·a)，其空間分布呈西南高、東北低的趨勢，這是因為西南地區植被覆蓋度較高，而東北地區的土地利用類型主要為水田、旱地、交通運輸用地、水域及水利設施用地以及其他用地等，植被覆蓋度較低，受人類活動干擾較大。

(2) 土壤保持服務功能

白馬雪山保護區土壤保持服務功能為0～3 689.15t/(hm2·a)，總體較低，其空間分布呈南高北低的趨勢，南部地區土地利用類型主要是林地和草地，植被覆蓋度較高且土壤侵蝕量較小，水土保持情況較好，而北部地區人類活動影響強度較大，因此水土流失比較嚴重，土壤保持能力相對較弱。

(3) NPP服務功能

白馬雪山保護區NPP服務功能為169.27～431.68 g/m2，其空間分布呈四周高中間低的趨勢，特別是西北部、東北部、東南部區域都較高，因為四周的蔓延度較高，植被資源豐富。

(4) 生境支持服務功能

白馬雪山保護區生境支持服務功能為0～1，其空間分布呈北低南高的趨勢。

2.3 景觀格局與生態系統服務的相關性分析

2.3.1 景觀格局指數與生態系統服務功能的全局空間自相關性

由表1可見，白馬雪山保護區NPP服務功能、水源供給服務功能、生境支持服務功能3種生態系統服務功能都與蔓延度呈空間正相關關系，而與斑塊數、邊緣密度、香農多樣性呈空間負相關關系；土壤保持服務功能與蔓延度、斑塊數、邊緣密度呈空間正相關關系，與香農多樣性呈空間負相關關系。

表1 景觀格局指數與生態系統服務功能的莫蘭指數

2.3.2 景觀格局指數與生態系統服務功能的局部空間自相關性

白馬雪山保護區水源供給服務功能與斑塊數、邊緣密度、香農多樣性的高-高聚集區在中部偏西北地區，說明該區域景觀斑塊類型復雜多樣且水源供給服務功能好；高-低聚集區在南部偏東地區，該區域水源供給服務較好但景觀斑塊類型較為單一；低-高聚集區集中在北部地區，該區域水源供給服務不好但景觀斑塊類型較為復雜多樣；水源供給服務功能與蔓延度在空間上聚集區較少且分散?？偟膩碚f，北部地區斑塊數、邊緣密度、香農多樣性對水源供給產生正向影響；南部地區斑塊數、邊緣密度、香農多樣性對水源供給產生負向影響。

NPP服務功能與斑塊數、邊緣密度、香農多樣性的高-高聚集區分布在北部地區，景觀類型豐富，景觀邊緣形狀相對不規整，景觀破裂程度較高，說明斑塊數、邊緣密度、香農多樣性對NPP服務功能產生了正向影響；高-低聚集區分布在南部地區，景觀破裂程度較小，斑塊類型較為單一而NPP服務功能則較高；NPP服務功能與蔓延度的高-高聚集區分布面積很少，主要分布于東部和東北部的部分地區。總的來說，北部地區斑塊數、邊緣密度、香農多樣性對NPP服務產生正向影響；南部地區斑塊數、邊緣密度、香農多樣性對NPP服務產生負向影響。

土壤保持服務功能與蔓延度、斑塊數、香農多樣性幾乎沒有高-高聚集區；低-高聚集區分布在北部地區，這3個景觀格局指數對土壤保持服務產生負向影響；低-低聚集區分布在南部地區，這3個景觀格局指數對土壤保持服務產生正向影響；土壤保持服務功能與蔓延度的空間聚集關系不明顯?？偟膩碚f，北部地區蔓延度、斑塊數、香農多樣性對土壤保持服務產生負向影響；南部地區蔓延度、斑塊數、香農多樣性對土壤保持服務產生正向影響。

生境支持服務功能與景觀格局指數的高-高聚集區較少且比較分散；生境支持服務功能與蔓延度、斑塊數、香農多樣性的低-高聚集區分布在北部地區，這3個景觀格局指數對生境支持服務產生負向影響；高-低聚集區分布在南部地區，生境支持服務功能較好，3個景觀格局指數對生境支持服務產生負向影響?？偟膩碚f，景觀格局指數對生境支持服務產生負向影響。

3 討 論

3.1 景觀格局與生態系統服務及其空間關系研究方法探討

以往研究一般采用標準法計算景觀格局指數，雖然能夠定量反映景觀結構組成和空間配置特征的簡單量化指標，但不能從空間上反映景觀格局的特征[13]。本研究使用Fragstats4.2軟件移動窗口法計算景觀格局指數，分析景觀格局空間分布特征，可以實現區域尺度下景觀格局指數空間分布的可視化，形成景觀格局指數柵格分布圖，使景觀格局的信息更加明朗[14]，并為進一步研究景觀格局指數與生態系統服務功能空間關系提供數據基礎。對于生態系統服務功能計算，以往一般使用物質量、價值量和能值量方法等，結果的空間可視化表達和分析比較困難[15]，而本研究通過空間數據量化生態系統服務功能，并在地圖上以可視化的形式定量表達生態系統服務功能。景觀格局和生態系統服務的空間關系分析上，目前國際上較為常用的空間關系統計分析軟件有R軟件和GeoDA軟件，R軟件只能以散點圖的形式顯示，不能以地圖的形式顯示局部聚集性[16]，本研究將GeoDA軟件運用于生態系統服務及景觀格局關系研究，可以生成局部空間自相關的顯著圖和聚集圖，從而可以更加直觀的看出生態系統服務和景觀格局的聚集關系，在景觀生態空間分析中有很好的適用性。

3.2 景觀格局與生態系統服務及其空間關系結果探討

對于景觀格局指數，白馬雪山保護區北部和邊緣地區景觀格局趨于破碎化，這是因為這些地區受人類活動干擾較多，用地類型主要為農村宅基地、水田、旱地，人口密度較大，人類活動頻繁；南部和中部地區景觀格局破裂程度較小，這是因為這些區域主要用地類型為林地和草地，人口密度較小，受人類活動干擾較少，而且保護區成立后更是加強了對林區的保護。

對于生態系統服務功能，本研究利用InVEST模型計算的土壤保持服務功能為0～3 689.15 t/(hm2·a)，水源供給服務功能為124.94～1 094.56 mm/(hm2·a)，生境支持服務功能為0～1，與孫興齊[17]用InVEST模型計算的香格里拉生態系統服務功能較為接近。這是因為這兩個地區生態環境比較類似，且都使用的是InVEST模型，同時也說明了用InVEST模型計算生態系統服務功能的可行性。

在景觀格局指數和生態系統服務功能空間關系分析上，與張金茜等[18]、王航等[19]、董敏等[20]的研究結果相似，本研究也表明生態系統服務與景觀格局存在相關性，研究結果合理、科學。不過，不同景觀格局指數和生態系統服務功能的影響比較復雜，既有正向的影響也有負向的影響[21]。因此，景觀格局與生態系統服務的相互作用還需要進一步研究。

4 結 論

(1) 景觀格局指數斑塊數、邊緣密度、香農多樣性在保護區北高南低，而蔓延度除呈南高北低的趨勢外還表現出四周高中間低的趨勢，景觀破裂程度高導致景觀多樣性增加。

(2) 生態系統服務功能中水源供給服務功能為124.94～1 094.56 mm/(hm2·a)，土壤保持服務功能為0～3 689.15 t/(hm2·a)，NPP服務功能為169.27～431.68 g/m2，生境支持服務功能為0～1。

(3) 白馬雪山保護區NPP服務功能、水源供給服務功能、生境支持服務功能3種生態系統服務功能都與蔓延度呈空間正相關關系，而與斑塊數、邊緣密度、香農多樣性呈空間負相關關系；土壤保持服務功能與蔓延度、斑塊數、邊緣密度呈空間正相關關系，與香農多樣性呈空間負相關關系。不同景觀格局指數和生態系統服務功能的影響比較復雜，既有正向的影響也有負向的影響。