劉 洋,畢 軍,呂建樹
1 濟南大學商學院,濟南 250002 2 山東龍山綠色經濟研究中心,濟南 250002 3 南京大學環境學院,南京 210023 4山東師范大學地理與環境學院,濟南 250014
生態系統服務是指人類直接或間接從生態系統中獲取的各種惠益,涵蓋人類賴以生存的物質基礎和環境條件。多種生態系統服務在復雜的影響因素作用下往往表現出此消彼長的權衡關系(trade-off)和相互促進的協同關系(synergy)[1-2]。識別生態系統服務間的相互關系,厘清其時空變化特征及驅動機制,可協同優化多種服務,避免在提高一種服務的同時損害其他服務,對促進區域生態系統可持續管理、指導自然資源的合理開發、提高人類福祉具有重要意義[3],也是生態系統服務從理論研究走向管理實踐的關鍵環節[4-5]。
近年來,生態系統服務及相互關系研究已成地理學、生態學及生態經濟學等多學科的前沿領域[1,6-7]。國內外學者普遍采用統計分析法[8-9]、空間制圖法[10-11]、情景模擬法[12-13]及模型量化法[14-15]等方法,借助InVEST、ARIES、Envisio、EcoAIM等軟件[6,16],對生態系統服權衡/協同關系的時空表現特征、尺度效應、形成機制及不確定性等方面進行研究[8,17]。結果表明,在氣候、土地利用、社會偏好、激勵政策等因素[18- 20]的影響下,生態系統服務權衡與協同關系表現為空間異質性與時間動態性,并且具有尺度依賴性[21]和研究不確定性[22]。然而,當前研究尚處于起步階段[23],主要集中于供給服務和調節服務之間的權衡關系,且缺乏精確的數量模型,急需通過增加數據樣本、提高數據質量、優化研究方法等方式,進行生態系統服務關系的度量。同時,驅動機制的研究也較薄弱,僅局限于一種或幾種特定影響因素分析[12,24],缺少對驅動力的綜合解析及主導驅動力識別,從而限制生態系統服務關系研究在管理實踐中的應用。
流域作為復合生態系統,具有多樣的服務類型,在自然和人為因素的共同作用下,生態系統服務之間表現出復雜的時空關系特征。同時流域也是人類活動干擾最強烈的地區之一,面臨水質退化、水量不足、水土流失等問題,如何協調經濟發展與環境保護的關系、緩解上下游利益相關者之間的矛盾成為流域可持續發展的關鍵。生態系統服務是銜接自然環境與人類需求的重要紐帶,通過研究流域多種生態系統服務的權衡與協同關系,明確其驅動機制,有助于制定流域發展與生態保護“雙贏”的政策措施[1,25]。基于此,本文以太湖流域江蘇省為例,考慮到農業非點源污染是該區域的首要水環境問題,選擇氮、磷凈化作為主要服務,其次是水量供給與土壤保持服務;通過空間分布式的生物物理模型與經驗統計模型,計算研究區每種生態系統服務的數量值;在此基礎上,借助GIS空間分析和數理統計法,研究2種水體凈化服務與其他服務之間的權衡/協同關系。最后,結合驅動因子空間柵格化,利用Logistic回歸模型,定量分析影響生態系統服務關系的主導驅動力,以期為流域的水環境管理及生態環境保護奠定基礎。
太湖流域地處長江三角洲城市群,是中國經濟最發達的地區之一。研究區選擇太湖流域的江蘇省(圖1),因其位于相對完整的行政區和獨立的水利分區[26],便于研究數據收集及生態環境模擬。該區北靠長江,地勢西高東低,平均氣溫15—17 ℃,多年平均降雨量約為1181 mm,氣候溫暖濕潤,是中國重要的糧食產區;總面積為16984.48 km2,包括蘇州市、無錫市、常州市、鎮江市及南京市的部分地區,城市化水平較高。近年來,隨著社會經濟發展,區域人地矛盾加劇,生態環境日益退化,非點源污染尤為嚴重。研究與非點源污染物產生、排放、傳輸、凈化有關的生態系統服務類型,明確其相互關系的時空特征及驅動機制,從生態系統綜合管理的角度為區域水環境管理,尤其是非點源治理提供新思路。

圖1 研究區地理位置及基本狀況圖Fig.1 Location and basic situation of the study area
生態系統服務及驅動力空間量化所需的數據包括降雨、溫度、太陽輻射等氣象數據,土壤類型、理化性質、含水量、深度等土壤數據,用于植被覆蓋度計算的遙感數據,土地利用及道路的空間數據,人口密度、GDP、農業等社會經濟數據,以及用于模型校正的水文、水質數據等,具體來源如表1所示。
1.3.1水體凈化服務模型
生態系統的植被、土壤等要素通過儲存、轉換等方式移除徑流中的營養鹽污染物,從而減少排放到外界的污染物,達到凈化水體的作用。生態系統凈化功能越強,營養鹽輸出量越小。因此,利用氮輸出量(Nitrogen Export, NE)和磷輸出量(Phosphorus Export, PE)表征生態系統氮、磷凈化服務,這兩個反向指標值越大表明水體凈化服務越低,反之亦然。本文利用InVEST軟件的“營養鹽凈化”模型,結合經驗統計模型、專家咨詢及文獻調研等方法修正參數,以實現對研究區2種水體凈化服務的空間量化。由于該模型考慮了不同地類氮、磷營養鹽參數的差異性,在一定程度上降低了服務指標之間的相關性,可用于后續的相互關系分析。模型原理如下:

表1 本研究的數據來源
(1)
式中,Expx為上游柵格x的氮、磷輸出值(kg/hm2);Ey表示每個下游柵格y過濾氮、磷營養鹽的效率(%);X代表氮、磷輸出量從上游產生到下游水體的運輸路徑;ALVx為校正的柵格x的氮、磷輸出值(kg/hm2),計算公式如下:
ALVx=HSSx×polx
(2)
(3)
(4)

1.3.2水量供給服務模型
生態系統通過水循環產生可供人類蓄水、發電等利用的水量供給服務(Water Supply, WS),并且水循環過程會溶解、析出氮、磷營養鹽,進而間接影響區域非點源污染。本文采用InVEST軟件的“水產量”模型來計算研究區的地表產水量,實現水量供給服務的空間化。模型原理如下:
WS(x)=Prex-AETx
(5)
式中,WS(x)為柵格x的年產水量(mm);Prex為柵格x的年降雨量(mm);AETx為柵格x的實際蒸散量(mm)。該公式可進一步表示為:
(6)
式中,AETx/Prex為基于Budyko曲線的表達方式,其計算公式如下:
(7)
式中,ωx為非物理參數,Rx為柵格x的干燥指數,各自計算公式如下:
(8)
(9)
式中,Z為降雨時間分布和水文生態特征的經驗常數;AWCx為柵格x的植被有效含水量(mm),由公式(10)求得。kj為不同地類的蒸散系數;ET0為潛在蒸散量(mm/d),通過Modified-Hargreaves法進行計算。
AWCx=Min(soil.depth, root.depth) ×PAWC
(10)
式中,soil.depth為限制根系生長的土壤最大深度(mm),root.depth為根系長度(mm),PAWC為植物有效持水量(%)。
1.3.3土壤保持服務模型
生態系統在預防土壤侵蝕和保持土壤過濾污染物功能等方面發揮重要作用,進而影響區域非點源污染。本文使用RUSLE模型估算土壤保持服務(Soil Retention, SR),并基于小流域的實驗數據對其進行校正,以得到本地化的計算模型,如下:
SR(x)=R×K×LS×(1-C×P×A)
(11)
式中,R、K、LS、C分別為降雨侵蝕力因子、土壤可蝕性因子、坡長坡度因子、植被覆蓋因子,計算公式如下所示。P表示工程管理措施因子,為0—1的無量綱數,各地類賦值參考呂建樹等的研究[27]。A為本地化修正系數,用來調整土地利用類型和地表光滑度對土壤侵蝕的影響,根據研究區實測數據確定值為1.3。
(12)
式中,p為年平均降雨量(mm),pi為月降雨量(mm)。


(13)
式中,Sa、Si、Sn分別為砂粒、粉砂與粘粒的含量值(%);C表示有機碳含量值(%)。
(14)
式中,Fa為匯水累積閾值,Cs為柵格大小,s為坡度,n為坡度相關參數,β為流向相關參數。

(15)
式中,C為0—1的無量綱數,c為植被覆蓋度,通過歸一化植被指數(NDVI)求得。
1.3.4驅動力分析方法
考慮研究數據的可得性,在多因子相關分析的基礎上,結合專家咨詢方法,最終確定氣候、土壤、經濟、用地等10類31種驅動力,具體如表2所示。
在生態系統服務關系和驅動力空間量化的基礎上,將各自柵格值賦予GIS生成的10000個隨機點,形成生態系統服務關系為因變量,驅動力為自變量的數據集,再運用多元Logistic回歸模型進行驅動機制研究。模型原理如下:
(16)
式中,pi為生態系統服務權衡/協同關系發生的概率,xi為驅動力;α為截距,即模型的常數項;βi為斜率,即模型中每個驅動力的偏回歸系數,其值為正則表示因變量與自變量的變化方向相同,負值則相反。
為更好的測量變量之間的關聯性,明確自變量對因變量發生概率的作用程度,一般應用p/(1-p)來解釋回歸系數,即發生比率(odds ratio, OR),其計算公式如下:
OR(p)=exp(α+β1X1+β2X2+,…,+βnXn)
(17)
本文選擇Homsmer-Lemeshow(HL)指標進行模型擬合度檢驗,當HL值不顯著表示模型擬合好,而顯著則說明模型擬合不好。通過Wald統計量檢驗模型的回歸系數,取0.05為顯著水平;當自變量的Wald所對應的P值小于0.05時,保留在方程中,而大于0.05則因未通過顯著性檢驗而被去除。

表2 生態系統服務權衡與協同關系的驅動力
本文計算了2000、2005及2010年研究區4種生態系統服務指標的物理量,如圖2。其中,綠、黃、紅色分別表示服務供給的高、中、低值區,藍色代表水體。結果顯示,氮凈化服務高值區(氮輸出量低值區)主要分布在研究區南部、西部及湖體周圍,低值區(氮輸出量高值區)則零散分布于北部及中部,中值區分布廣泛。從2000到2010年,氮輸出量單位面積年均值呈現先增加后略有下降的趨勢,增幅由8.89%降低為8.08%。磷凈化服務高值區(磷輸出量低值區)主要分布于研究區周邊的林地、草地、濕地等生態用地,低值區(磷輸出量高值區)在農村居民點及城鎮建設用地附近突出,并零散分布于部分農耕地。2000—2010年,磷輸出量年均值逐漸增加,增幅為24.69%。過量的化肥施用、畜禽養殖等農業活動會增加氮、磷的排放,而農村生活污染和城鎮污水也會顯著影響磷素富集[28]。因此磷凈化服務分布與農村居民點、城鎮建設用地存在明顯關系,進而導致氮、磷指標時空變化的差異性。已有研究也表明磷負荷與城市用地、農村建設用地密切相關[29-30]。水量供給服務的空間分布呈東西走向,且變化顯著,高值區主要分布于研究區周邊,中值和低值區由研究區西部逐漸轉向中、東部。2000—2010年,水量供給服務年均值呈先下降后上升的趨勢,降幅由-27.88%變為-12.76%。研究區土壤保持服務普遍較低,低值區分布廣泛,高值區零星分布。2000—2010年,土壤保持服務年均值呈下降趨勢,降幅為-21.61%,主要受降雨、日照等自然因子[31]以及農田、林地、城鎮等用地變化[32-33]的影響。

圖2 2000—2010年研究區4種生態系統服務時空分布Fig.2 The temporal and spatial distribution of four ecosystem services in study area from 2000 to 2010
在ArcGIS軟件中將2010年和2000年每種生態系統服務相減,得到各自的空間變化圖,并將變化分為正向、負向及未變化3類。利用GIS的空間疊加分析功能得到兩種服務變化的空間分布格局(圖3)。其中,正向和負向變化區為兩種服務均增加或均減少的區域,表明兩者具有相同變化方向的協同關系;相反變化區為兩種服務變化方向相反,表明兩者具有權衡關系;未變化/正(負)向變化區為一種服務沒變化而另一種服務增加(減少),表明兩者關系不顯著;未變化區是兩種服務均未變化,其關系也不顯著。統計每類變化區占區域總面積的百分比,結果見表3。本文以水體凈化服務為主,分析該類服務指標與其他服務的相互關系。

圖3 2000—2010年生態系統服務間關系的空間格局Fig.3 Spatial pattern of ecosystem services relationships between 2000 and 2010
結果表明,氮輸出指標與水量供給關系中,相同變化區(包括正、負向變化區,表示氮凈化與水量供給服務的權衡關系區)面積占統計面積的33.16%,其中正向變化區占11.97%,主要分布于研究區西北部、西南角及東南角地區,負向變化區占21.19%,與城鎮建設用地分布較一致;相反變化區(表示氮凈化與水量供給服務的協同關系區)占統計面積的29.98%,主要分布于研究區中、東部。磷輸出指標與水量供給關系中,相同變化區和相反變化區的面積占比分別為27.21%和31.92%,在研究區西北部、西南角及中部地區集中分布;其次是未變化/正向變化區,占統計面積的26.43%,主要位于研究區的中西部及北部一段沿江地區。氮、磷輸出指標與土壤保持關系中,未變化/正向變化區分布范圍均最廣,分別占統計面積的48.06%和55.37%,而正向、負向及相反變化區所占比例較小,說明水體凈化服務與土壤保持服務的權衡/協同關系不明顯。氮與磷指標關系中,正向變化區廣泛分布,占統計面積的58.86%,說明研究區水體凈化服務主要表現為協同關系。

表3 2000—2010年生態系統服務間關系的面積比
該統計面積沒考慮研究區未變化的水域面積;單位為%
選擇顯著的生態系統服務關系進行驅動力分析,主要為氮凈化服務與水量供給服務的權衡關系、兩種水體凈化服務的協同關系。通過多元Logistic回歸模型進行定量分析,模型中因變量是兩種服務關系,取值為0和1值,0表示兩者關系不顯著,即一種變化而另一種未變化或者兩者均未變化;1表示兩者具有協同關系(均正向變化或負向變化)或權衡關系(相反變化)。據此,將生態系統服務關系的空間格局圖進行重分類,得到0和1的因變量圖層。同時,在多期驅動因子空間柵格化的基礎上計算其平均值,得到自變量圖層。將因變量和自變量的圖層值賦予隨機點,基于點的屬性值確定Logistic回歸模型的輸入值。由于數據來源于點的隨機抽樣,避免了空間自相關性;同時在隨機點選取時,盡量保證因變量的0和1值數量相等,以避免對模型常數項的影響。模擬結果如表4、5所示。
表4模型的HL指標值為0.817,統計不顯著,反映模型的擬合度較好;模型預測正確率為81.32%,表明模型較穩定。根據Wald統計量可判斷,影響氮凈化和水量供給服務關系的最重要因素是城鎮建設用地密度和植被覆蓋度,其次是離農村距離、農村居民點密度、水網密度、離城鎮距離及人口密度等。其中,城鎮建設用地密度和農村居民點密度的回歸系數為正,說明每平方公里增加1hm2的城鎮或農村建設用地,將使氮凈化和水量供給服務關系權衡變化關系的發生概率分別增加29.100倍和10.282倍。而植被覆蓋度和水網密度的回歸系數為負,表明此權衡關系將隨著這些因素的增加而減小。此外,年相對濕度、年均氣溫等氣候因素的Wald統計變量也比較顯著,說明氣候因素對兩者權衡關系也具有一定的影響,尤其是年平均氣溫對其貢獻率達到28.979。土壤因素中,總磷含量是最重要的解釋因變量,且其發生比率也最大,達到37.490。社會經濟因素中,除人口密度的貢獻率最大外,農業GDP也是較為重要的影響因素。農業GDP增加可視為當地農業發展水平提高,此時將帶來兩者權衡關系的下降,其發生概率為0.997倍。因此,當氣候變化帶來區域水量供給增加時,由于農業的發展,可采取適應性措施,減弱對區域生態系統水體凈化服務的影響。

表4 氮凈化與水量供給服務關系的驅動力模型估計結果
表5模型的HL指標值為1.017,統計不顯著,反映模型擬合度較好,同時模型也較穩定,其預測正確率為89.49%。根據Wald統計量可判斷,影響兩種水體凈化服務協同關系的最重要驅動力是植被覆蓋度,其次是耕地比例、林地比例、離農村距離、水網密度,并且這些解釋變量的回歸系數均為正,表明這些因素的增加將帶來兩種服務協同關系的提高。從發生比率來看,驅動因素的貢獻率均不大,普遍在1.000左右,僅有農村居民點密度的概率達2.406。總體來看,所有驅動力中植被用地的影響最明顯,通過增加耕地、林地及草地的面積或建立植被緩沖帶將顯著促進兩種水體凈化服務的協同變化;而增加林地、草地等生態用地可有效提高生態系統功能,凈化多種水質污染物,達到協同減排的目的。此外,土壤砂粒、粘粒等顆粒物組成也會促進兩者的協同關系。因此,通過調節土壤屬性或改良土壤結構可同時改變兩種水體凈化服務,為區域非點源污染物的協同治理提供有效途徑。

表5 氮、磷凈化服務關系的驅動力模型估計結果
復雜的內外界環境因素作用于生態系統的結構、過程及功能,進而影響生態系統服務供給,表現為多種服務在時空尺度上的同增同減、一增一減或無明顯變化關系。生態系統服務間相互關系研究有助于管理者識別多種服務的變化特征,明確區域生態系統的變化程度,以采取生態空間管制、環境污染防治等針對性政策措施進行調控,促進區域生態系統的健康、持續發展。流域是自然環境和人類活動交互作用強烈地區,面臨諸多水環境問題,結合流域的實際情況,選擇關鍵的生態系統服務,進行多種服務之間的權衡與協同關系研究,為實現流域生態系統可持續管理奠定基礎。
本文以太湖流域江蘇省為研究區,運用空間顯示的生物物理模型和經驗統計模型計算區域4種生態系統服務指標,在此基礎上借助GIS空間分析手段進行服務權衡與協同關系研究,并利用多元Logistic模型定量識別影響生態系統服務關系的主導驅動力。研究結果表明:
(1) 從2000到2010年,研究區4種生態系統服務指標的單位面積年均值呈現不同程度的變化,氮輸出指標先增加后略有下降的趨勢,磷輸出指標逐漸增加,水量供給先下降后上升,土壤保持則呈現逐漸下降趨勢。研究區4種服務指標的空間分布及變化格局具有較明顯的空間異質性。
(2) 氮、磷輸出指標與水量供給關系的空間分布格局較為相似,表現為廣泛分布的相同變化區和相反變化區;氮、磷輸出指標與土壤保持的正向、負向及相反變化區占比均較小,表明水體凈化服務與土壤保持服務的權衡/協同關系不明顯;氮、磷輸出指標的正向變化區廣泛分布,說明研究區水體凈化服務主要表現為協同關系。
(3) 氮凈化和水量供給服務權衡關系的正向驅動力為城鎮建設用地密度和農村居民點密度,且每平方公里增加1公頃的城鎮或農村建設用地,其權衡關系的發生概率將增加29.100倍8和10.282倍;而植被覆蓋度和水網密度則具有顯著的負作用。氮、磷服務協同關系的首要影響因素是植被覆蓋度,其次是耕地、林地比例等因素,且這些因子均表現為正向促進作用。
盡管生態系統服務的研究已逐步深入,然而從理論走向實踐仍面臨一系列挑戰[1,3, 5],尤其是如何更精準的計算服務量,如何科學識別多種服務之間相互關系,如何定量分析服務關系的驅動機制等問題。本研究綜合利用生態、地理、統計、計量等方法手段,從時空多維度進行生態系統服務關系的表征及驅動力分析,以期拓展相關研究,為進一步理解生態系統服務之間復雜關系及其驅動機制提供借鑒,并促進生態系統服務研究在流域水環境管理中的應用。生態系統服務具有時間積累性,而不同類型服務的變化速度和閾值也存在差異,影響服務關系的時間表現特征。由于數據的有限性,本文僅進行了2000—2010年的研究,未能充分考慮每種生態系統服務變化的周期。未來可結合生態系統變化的閾值分析進一步探究不同時間周期及更長時間范圍上生態系統服務之間的相互關系,并探索驅動力變化程度對生態系統服務關系的影響。