祁連山地區生態系統服務間權衡的社會-生態環境響應機制研究

呂榮芳，趙文鵬，田曉磊，張建明

（蘭州大學 資源環境學院，甘肅蘭州 730000）

0 引言

生態系統服務指人類直接或間接從生態系統獲取的各項產品和服務［1］。它是連接自然生態環境和人類社會的橋梁，其可持續利用是實現生態保護和高質量發展的前提和基礎［2-3］。千年生態系統評估的結果顯示，全球60%的生態系統都處于退化狀態，評估氣候變化和人類活動對生態系統服務的影響，對生態系統的科學管理和政策制定具有重要意義［4-6］。生態系統服務與人類福祉息息相關，已成為地理學、生態學、經濟學等多學科交叉研究的前沿和熱點問題，國內外多位學者從概念定義［7］、分類體系［8］、定量評估［1］、權衡分析［9］、供需耦合［10］、驅動機制［11］、政策制定［12］等方面開展了大量研究。生態系統中，碳循環、水循環、污染降解等生態過程相互依賴、相互影響，生態系統服務依此而生，因此不同服務間通常存在復雜關聯，包括此消彼長的權衡關系和同升同降的協同關系［13］。服務間關聯分析可為土地利用管理和生態規劃提供關鍵信息，逐漸成為相關研究的熱點［14］。

生態系統服務間的關聯存在一定的非線性，會隨時間和空間分布而變化，因此難以通過簡單的統計分析方法來進行計算［15-17］。現有研究多采用相關性分析、主成分分析、聚類分析、疊置分析等方法來探索服務間的權衡和協同關系，簡單易行，但難以表征服務間的非線性關聯及其空間異質性［13，18］。Mouchet等［19］提出現行的分析方法僅適合提供關于服務的變化趨勢和相對高低的初步信息。相關性分析和主成分分析難以判別非線性關系，通常是在全區尺度上分析整體關聯，不足以反映服務間關聯的空間異質性［15］。疊置分析的結果存在很大的不確定性，會因個人觀點而變化，降低了不同地區研究的可對比性［13］。因此，迫切需要引入新的方法來探索服務間關聯的非線性和空間異質性。另一方面，服務的時空變化會受到降雨、氣溫、土壤、植被分布等自然生態環境和經濟發展、人口分布、城鎮擴展等社會經濟條件的綜合影響［5，20］。定量化探索生態系統服務及其關聯的社會-生態驅動機制，對協調服務間權衡、制定生態管理政策具有重要意義［21-22］。現有研究多對單項服務的社會-生態驅動機制進行分析，缺少對服務間關聯的形成機制的研究，難以對協調服務間權衡、增強服務間協同提供有效的管理建議［23］。同時，現有研究多在區域尺度上對服務的驅動機制進行整體研究，對服務間非線性關聯的空間異質性的研究較少，尤其是造成服務間權衡的關鍵環境因子及其閾值效應，不足以對土地利用管理和政策制定提供有效建議［15］。

山區占據全球陸地生態系統中大概1/4的地表面積，為全球大概12%的人口提供住所，是全球超過50%的物種的熱點分布區［24］。山區可為全球提供糧食、木材等原材料，可防護泥石流、風沙侵蝕等自然災害，能夠保存特有文化，并為多樣化物種提供生境，因此山區生態系統服務可全面表征山區為人類社會提供的服務和效益［25］。祁連山是我國西部重要的生態屏障，承擔著水源涵養、生物多樣性保護等生態重任，被譽為河西走廊的“母親山”［26］。它是我國干旱區與非干旱區的分界線，也是我國200 mm年等降水量線，其生態重要性不言而喻［27-28］。本研究對祁連山地區重要的生態系統服務進行研究，分析服務間非線性關聯的空間差異性及其對社會-生態環境的響應機制，為祁連山地區的生態管理和修復提供科學依據。

1 研究區概況

祁連山地區（93°56′～103°90′E，35°84′～39°97′N）地處內蒙古高原、青藏高原和黃土高原的交匯處，是疏勒河、黑河、石羊河三大內陸河流域的水源涵養地（圖1）。它橫跨甘肅、青海兩省，總面積約為19.47×104km2，海拔1 937～5 792 m，高程自東南向西北呈降低趨勢，是我國西部地區重要的生態屏障。內有雪山、冰川、森林、灌叢、草原、草甸、沼澤、濕地等多種土地利用和植被覆蓋類型，也是我國重要的野生動物遷徙廊道，社會-生態環境多樣［29-30］。

圖1 祁連山地區位置和高程圖Fig.1 The location of Qilian Mountains in China and its topography

2 研究方法

2.1 生態系統服務的計算

基于研究區的社會-生態環境特征和主要生態問題，本研究重點分析可以緩解全球變暖趨勢的碳服務（碳固定和碳儲存）［31］、水文相關服務（產水量、土壤保持和水質凈化）和物種多樣性保持的支持服務（生境質量）。主要數據來源如表1所示。

表1 主要數據來源Table 1 The main data sources

2.1.1 碳固定

碳固定服務利用InVEST模型的Carbon模塊進行計算，分為地下生物碳、土壤碳和死亡有機碳，以土地利用類型為評估單元，利用各類碳在不同土地利用類型內的平均碳密度進行評估［32，33，34］。其計算公式如下：

式中：Ctot為碳固定量，Cbelow為地下生物碳儲量，Csoil為土壤碳儲量，Cdead為死亡有機碳儲量。不同土地利用類型內各類碳的平均碳密度如表2所示。

表2 不同土地利用類型內各類碳庫的平均碳密度（單位：Mg·hm-2）Table 2 The average carbon densities of three pools in different land use types（unit：Mg·hm-2）

2.1.2 碳儲存

碳儲存通過植被凈初級生產力（Net Primary Productivity，NPP）進行表示，指的是地表植被在單位面積單位時間內積累的有機物總量，為植被通過光合作用產生的有機物量減去自養呼吸消耗量。本研究以朱文泉等提出的基于CASA模型的ENVI遙感計算模塊（zhu_npp_calc）進行計算［35］。其計算公式如下：

式中：NPP（x，t）為像元x在t月的植被凈初級生產力g C·m-2·a-1，APAR(x，t)為像元x在t月的光合有效輻射（g C·m-2），ε(x，t)為像元x在t月的實際光能利用率（g C·MJ-1），SOL(x，t)為像元x處在t月的太陽總輻射，FPAR(x，t)為像元x處在t月的光合有效輻射吸收率，常數0.5表示植被所能利用的太陽有效輻射占太陽總輻射的比例，Tε(x，t)、Wε(x，t)分別表示溫度和水分脅迫系數，εmax為在理想條件下植被的最大光能利用率。

2.1.3 產水量

產水量通過InVEST模型Water Yield模塊進行計算。該模塊基于水量平衡原理，綜合降雨、地表蒸發、植物蒸騰、植物最大根系深度以及可用水量等參數，假定各柵格所有產水量都經過地表或地下徑流到達流域出口，最終將各柵格的產水量在子流域尺度上進行加權或平均［36］。其計算公式如下：

式中：Yxj為土地覆被類型j上柵格單元x的年產水量（mm）；Px為柵格單元x的年降水量（mm）；AETxj為實際年平均蒸散發量（mm），ωx為氣候與土壤屬性關系的非物理參數；Rxj為潛在蒸散發與降水量的比值；Kxj為植被類型的蒸散系數；ET0x為柵格單元x的蒸散發量（mm），ET0-PM為計算參考作物蒸發蒸騰量的Penman-Monteith公式，其中T為計算時段內的平均氣溫（℃），Δ為飽和水氣壓（kPa·℃-1），Rn為太陽凈輻射（MJ·m-2·d-1），G為土壤熱通量（MJ·m-2·d-1），本文G＝0，γ為濕度計常數（kPa·℃-1），es為飽和水氣壓、ea為實際水氣壓（kPa），μ2是高出地面2 m處的平均風速（m·s-1）。模型基于水文站點的監測數據進行多次校準。

2.1.4 水質凈化

水質凈化通過InVEST模型NDR模塊進行計算，僅考慮面源污染，基于長期穩定的污染物移動計算各像元的污染物攔截量，綜合高程、根系限制層深度、降雨量、土地利用、植被可利用水含量、平均潛在蒸散發等參數進行計算。參考相關研究獲取不同土地利用類型的污染物攔截系數和過濾量［37］。其計算公式如下：

式中：NEi為像元i中的污染物截留量（kg·pixel-1）；fj為各土地利用類型j的污染物截留容量；ALVj為像元i內基于水文敏感系數的污染物載荷量調整值；poli為像元i的輸出系數；λi為像元i的輸出指數；為平均輸出指數。

2.1.5 土壤保持

土壤保持通過InVEST模型SDR模塊進行計算，該模塊基于通用土壤流失方程USLE進行計算，主要包括上坡泥沙來源下地塊自身攔截量與植被覆蓋和水土保持措施下減少的土壤侵蝕量，以潛在土壤侵蝕量減去實際土壤侵蝕量來表示［38］。其計算公式如下：

式中：SEDRx為柵格x自身攔截的沉積物保留量（單位：t）；SEx、USLEy以及SEz分別表示為柵格x的截留率、上坡柵格y產生的泥沙量（單位：t）以及上坡柵格的泥沙截留量（單位：t）；SEDRETxD為土壤保持量（單位：t）；RKLSx與USLEx表示為潛在土壤侵蝕量（單位：t）和實際土壤侵蝕量（單位：t）。

2.1.6 生境質量

生境質量通過InVEST模型Habitat Quality模塊進行計算，基于各土地利用類型的威脅因子敏感度、外界威脅強度以及威脅因子的空間權重和影響距離進行計算［39-40］。計算公式如下：

式中：Dxj為生境脅迫水平；Hj為土地利用類型j的生境適宜性；k為半飽和常數，取Dxj最大值的一半；z為歸一化常數，取2.5。

2.2 生態系統服務間權衡和協同關系分析

首先，創建樣本數據集，利用ArcGIS 10.4軟件中的“創建隨機點”和“值提取至點”工具創建生態系統服務時間變化的樣本數據集，內含100 000個樣本點，不同點之間距離超過300 m，以降低其空間自相關性。第二，在區域尺度上，采用相關性分析和主成分分析定量化研究六項服務間的關聯。對樣本數據集中各變量進行z-score標準化處理，以消除數據量綱的影響。Kolmogotov-Smirnov檢測的結果顯示變量存在非正態分布，因此采用Spearman相關系數來計算其相關性，并通過主成分分析探索不同服務間關聯的接近程度。以上分析主要通過R語言軟件“corrgram”和“vegan”包實現。

多元回歸樹分析（Multivariate Regression Tree，MRT）為機器學習算法中二元回歸樹分析的拓展，可同時處理多個因變量，用于分析復雜的生態數據，并探索多元生態數據和環境因子間的關系。該方法以自變量為分類節點，利用二分類法將因變量逐步分為不同類型，可有效模擬變量間的非線性關系，已廣泛應用于研究物種分布與環境因子之間的關系，并以此劃分植物群落類型［41］。基于研究區的社會-生態環境特征，以選擇因子的代表性和數據的可獲得性為原則，參考服務的計算方法和相關研究，選擇對六項服務可能存在顯著影響的社會-生態環境因子［42］（表3）。以六項服務為因變量，以13項環境因子為自變量，進行多元回歸樹分析，以此探索服務間權衡和協同關系的空間異質性及其對社會-生態環境的響應機制。

表3 生態系統服務的可能影響因子Table 3 The potential socio-ecological factorsof ecosystem services

3 研究結果

3.1 區域尺度上生態系統服務的空間分布格局和服務間關聯

在區域尺度上，對祁連山地區六項生態系統服務的空間分布進行Spearman相關性分析和主成分分析，定量化探索六項服務的空間分布格局和服務間的整體關聯（圖2、圖3）。結果顯示，2019年，六項服務的空間分布存在一定的相似性，東南部相對較高，西北部相對較低，各項服務就空間分布而言均呈正相關性。六項服務中，碳固定、碳儲存和生境質量的空間分布較為相似，兩兩之間相關性均高于0.82，且在第一主成分上均有較大載荷，在東南部地區較高，在西南部地區較低，其中碳儲存的空間異質性相對最高。水質凈化與前三項服務也存在較高的正相關性，相關系數為0.65～0.85，同樣呈西低東高的分布格局，但在西部地區隨高程出現較大差異。此外，土壤保持和產水量的相關性為0.46，與其他服務間的相關性均小于0.42，顯示兩服務與其他服務的空間分布存在較大差異，其較高值主要出現在西北部的石羊河流域。總的來說，在區域尺度上，就空間分布而言，除了產水量與水質凈化間之外，六項服務間均存在顯著的協同關系。六項服務可分為三組：第一組為碳固定、碳儲存和生境質量，第二組為土壤保持和產水量，第三組為水質凈化，同一組內服務關聯更為緊密，空間分布較為相似，不同組間關聯性較低，空間分布差異較大。

圖2 2019年祁連山地區六項生態系統服務的空間分布圖Fig.2 The spatial distribution of six ecosystem services in Qilian Mountains in 2019

圖3 基于空間分布的六種生態系統服務的關聯圖及原理分析結果Fig.3 The corrgram and principle analysisresult of six ecosystem services

3.2 生態系統服務間關聯的空間差異及其環境響應機制

3.2.1 空間差異性

通過多元回歸樹分析，定量化研究六項服務及服務間關聯的空間差異性，并探索其環境響應機制。基于CVRE值最小表示預測能力最佳為原則，經交叉驗證，最終的多元回歸樹由5片葉子組成，模擬精度（R2）達0.904，模擬結果良好（圖4）。研究結果顯示，基于六項服務及服務間關聯的空間差異性，整個研究區可被分為5個聚類，回歸樹的每片葉子即為一個聚類。然后對不同聚類內各項服務的值域范圍進行對比（表4）。結果顯示，聚類1中，產水量和土壤保持處于全區平均值水平，其他服務均低于平均值，其中碳固定、產水量和生境質量遠低于其他聚類，而碳儲存、土壤保持和水質凈化處于較低水平。因此，六項服務均呈協同關系。聚類2中，產水量遠高于平均值，土壤保持略高于平均值，水質凈化略低于平均值，碳固定、碳儲存和生境質量遠低于平均值。其中，產水量和土壤保持要高于其他聚類，碳固定和水質凈化相對較低，而碳儲存和生境質量遠低于其他葉子。因此，聚類2中產水量和土壤保持與其他服務間呈權衡關系。聚類3中，碳固定和水質凈化處于全區平均值水平，產水量低于平均值，三者相較于其他聚類均處于中等水平，而碳儲存、土壤保持和生境質量略高于全區平均值，三者呈協同關系。聚類4中，碳固定和產水量略低于全區平均值，碳儲存和水質凈化遠高于平均值，土壤保持位于平均值水平，而生境質量遠低于平均值。其中水質凈化遠高于其他聚類，碳固定相對較高，產水量相對較低，而其他服務處于中等水平。因此，聚類4中水質凈化和碳固定呈協同關系，與產水量均呈權衡關系。聚類5中，碳固定和水質凈化處于全區平均值水平，碳儲存遠高于平均值，產水量、土壤保持和生境質量高于平均值，其中碳固定和生境質量遠高于其他葉子，水質凈化處于平均水平，而其他三項服務均處于略高水平。因此，除水質凈化外的五項服務在聚類5內均呈協同關系。

表4 多元回歸樹結果的定量統計Table 4 Qualitative interpretation of the resultsof Multivariate Regression Tree

3.2.2 社會-生態環境響應機制

結果顯示（圖4），12個社會-生態環境因子中三個因子的影響最為顯著，分別是土地利用格局（第一層和第三層的分割點）、降雨量（第二層的分割點）和植被覆蓋度（第二層的分割點），其他環境因子的影響被前三個因子所掩蓋或代替。在回歸樹的第一級分層中，依照土地利用類型為無植被覆蓋用地（水體、建筑用地、裸地和冰川/積雪）和有植被覆蓋用地（耕地、林地、草地、灌叢和濕地）而分割，前者依照年均降雨量高低被分為葉子1（＜440.2 mm）和葉子2（≥440.2 mm），而后者則依照植被覆蓋度高低進行第二層分級，植被覆蓋度較低（＜0.559）的為葉子3，植被覆蓋度較高的繼續進行第三級分層，其中土地利用類型為耕地的區域為葉子4，用地類型為林地、草地、灌叢和濕地的區域為葉子5。每片葉子即為一個聚類。聚類1主要分布在中西部地區和東部和青海湖地區，以水體、建筑用地、裸地和冰川/積雪等無植被覆蓋用地為主，且降雨相對較低。聚類2主要分布在東北部地區高程較高的地區，同樣無植被覆蓋，但降雨量相對較高。聚類3遍布整個研究區，以耕地、林地、草地、灌叢和濕地等有植被覆蓋用地為主，且植被覆蓋度相對較低。聚類4面積最小，僅占分布在黑河和石羊河流域內部分地區，以耕地為主，且植被覆蓋度相對較高。聚類5主要分布在研究區東部地區，以林地、草地、灌叢和濕地為主，且植被覆蓋度相對較高。

4 討論

4.1 生態系統服務間關聯的空間異質性

生態系統中碳循環、水循環、植被生長等生態過程相互依賴、相互影響，生態系統服務依此而生，因此不同服務間通常會存在一定關聯［13］。地形、地質、地貌、水文、氣候等生態環境因子和道路網、人口分布等社會經濟因子的時空異質性共同造就了服務間關聯的復雜性。該復雜性不僅存在于空間尺度上，即服務間關聯可能隨社會-生態環境不同而變化，還存在于時間尺度上，即服務間關聯可能從早期的協同關系變為后期的權衡關系［43］。因此，在制定相關的管理政策過程和土地利用規劃過程中，需要對服務間關聯的時空異質性進行詳加考慮。國內外多位學者對此展開了大量研究，現有結果顯示山區的調節服務間多呈協同關系［24］，這與本研究中全區尺度上的結果相一致，但與聚類2和聚類4內調節服務間多呈權衡關系的結果相反。這與兩區域內特殊的生態環境存在很大關系。喬木、灌木、草本等地表植被在生態系統中具有重要作用，它可通過光合作用從空氣中吸收大量碳，也可對降雨和地表污染物進行截留，同時固定土壤阻止土壤流失，在生態系統服務具有重要的生態功能和作用［20］。聚類2內地表無植被覆蓋且降雨量較高，因此產水量較高而碳儲存較低，二者因而呈權衡關系。聚類4的情況與此相反，植被覆蓋度較高，且多為耕地，對降雨的攔截作用較強，因而產水量較低而水質凈化較高，二者呈權衡關系。因此，山區內調節服務間關聯并不一定為協同，降雨量過高或單一的土地利用類型（耕地）會使服務間關聯從協同變為權衡，在降雨較高的山地區可鼓勵植被種植和恢復，在耕地區推行滴灌、水肥一體化等節水降污措施，對服務提高和服務間權衡關系的協調十分重要。

4.2 社會-生態環境因子對生態系統服務及服務間關聯的影響

研究結果顯示，在干旱區半干旱區的山地生態系統中，就碳固定、碳儲存、產水、水質凈化、土壤保持和生境質量六項服務而言，土地利用類型、降雨量和植被覆蓋度是生態系統服務供給的重要決定因素，其影響要超過其他環境因子。降雨量的高低更多的受到全球尺度上氣候的宏觀影響，在區域尺度上難以控制和改變［44-45］。相較而言，土地利用格局和植被覆蓋度為可控管理因子。研究結果可為土地利用管理和空間規劃提供有效參考，有助于政策制定者和管理者在不同環境條件下判別主要的保護對象和關鍵生態問題［46-48］。如在黑河流域的耕地區域，服務供給以高水質凈化和低產水量為主要特征。該地區農業生產以灌溉為主，依靠上游地區降水和冰川融水，其自身的產水服務不足以滿足日常生產和生活所需。而農業生產在追求經濟效益的同時，農藥化肥的大量使用帶來了嚴重的面源污染問題［49］，耕地地表植被覆蓋度較高，可截留部分污染，提供相對較高的水質凈化服務，但仍有一定的污染流入河流［50］。因此，在該地區農業生產過程中推行滴灌、水肥一體化等節水降污措施對提高生態質量具有重要意義。

山區地勢陡峭，其地表植被可有效緩沖和保護周邊低地免于滑坡、雪崩等自然災害的危害，也可截留大氣降雨和地表徑流以降低洪水威脅，保護土壤免遭侵蝕并增強土壤發育過程［51-52］。結果顯示，有植被覆蓋區內相較于無植被覆蓋區，除產水量可能較低外，其他五項服務都較高，表明了植被覆蓋在山區環境質量提高和生態效益供給方面的重要性。另一方面，在有植被覆蓋的地區，植被覆蓋度對生態系統服務的影響要高于其他地形、氣候、水文和人類活動因子，植被覆蓋度較低的地區碳儲存服務較低，這與秦嶺-大巴山區的研究結果相一致［40］。其次為土地利用類型，耕地相較于其他用地，因其地表植被較為密集且為氮磷等污染物的重要來源而提供較高的水質凈化服務，同時受人類活動的干擾較強，其環境和生存物種較為單一，因此提供的生境質量服務較低。

4.3 多元回歸樹分析的優劣

本研究通過多元回歸樹分析來判別生態系統服務間關聯的空間異質性和驅動機制，該方法可有效考慮服務間關聯的非線性和閾值效應。例如，研究區在年降雨量為440.2 mm處存在閾值效應，高于該值產水量與碳儲存間的關系從協同變為權衡，在生態管理政策制定過程中需加以考慮。第二，數理統計方法通常對輸入數據有所要求［53-54］，如皮爾遜相關系統分析要求數據呈正態分布［15］。相較于其他方法，多元回歸樹分析對輸入數據的要求較低，便于對多源數據進行處理，結果易于解讀，便于管理者和研究者就政策制定進行探討，有利于研究結果的推廣。第三，多元回歸樹分析可將社會-生態驅動因子和生態系統服務納入同一平臺上，判別產生服務聚類和服務間關聯的環境條件。但該方法也存在缺點，其中最主要的是未考慮服務的空間集聚性，相較于應用更為廣泛的自組織迭代聚類法，獲得的服務聚類在空間上分散性較高，特別是聚類3，這也是未來研究待解決的重要問題。

5 結論

本研究基于多元回歸樹分析法定量化分析2019年祁連山地區六項重要生態系統服務（碳固定、碳儲存、產水量、土壤保持、水質凈化和生境質量）間的非線性權衡和協同關系，揭示服務分布和服務間關聯的空間異質性，探索其社會-生態驅動機制和關鍵閾值。六項服務的空間分布主要受土地利用類型、降雨和植被覆蓋度的影響，整體均呈東南高、西北低的格局，其中碳固定、碳儲存和生境質量的空間分布更為相似，土壤保持和產水量更為相似，而水質凈化與其他服務的分布格局差異較大。在全區尺度上，六項服務整體呈協同關系，但該關聯在不同社會-生態環境下存在顯著的空間異質性。

基于六項服務空間分布和服務間關聯的差異性，祁連山地區可被分為五個聚類，各聚類內六項服務的供給特征和服務間關聯相對一致。無植被覆蓋度區內服務供給和服務間關聯受降雨量的影響更大，年均降雨量高于440.2 mm時，整體服務相對較高，產水量和土壤保持與其他服務間呈權衡關系。在植被覆蓋區內，服務供給受植被覆蓋度的影響更大，植被覆蓋度低于0.559或耕地區內服務間關聯減弱。

謹以此文，紀念恩師李吉均先生！