環京津城市群生態系統服務的時空演變與驅動因素

王凱平,馮 悅,邱陳瀾,王欣言,馬 嘉,張云路

北京林業大學園林學院,北京 100083

生態系統服務被廣泛定義為直接或間接促進可持續人類福祉的生態特征、功能或過程,并按照生態系統對人類產生的直接影響分為供給、調節、文化服務,以及對維持生態系統的其他功能具有重要作用的支持服務[1]。當前,伴隨著生態系統服務對人類發展的重要性被廣泛認識,探求城市群生態系統服務的演變規律與驅動因子關系逐漸成為學界的關注點之一。相關研究主要將生態系統服務驅動因素劃分為自然因素和社會經濟因素兩方面。在自然環境方面,紀建萬、閆國振、蘇常紅等學者研究發現,氣候、土壤、高程等自然氣候因素變化與生態系統服務價值演變存在明顯相關[2—4]。在社會經濟方面,徐煖銀,耿甜偉、璩路路等學者發現社會因素對生態系統服務價值有較大的影響,并在空間上與生態系統服務具有高度相關性[5—8]。由此可見,自然環境因素和經濟社會因素變化對生態系統服務有顯著影響和密切聯系。因此,識別城市群生態系統服務的主要驅動因素,可以為未來規劃者在城市群發展過程中做出決策提供依據。

2018年11月中共中央、國務院印發《關于建立更加有效的區域協調發展新機制的意見》,提出“以北京、天津為中心引領京津冀城市群發展,帶動環渤海地區協同發展”的目標。環京津地區在自然資源的提供、廢氣廢物的處理等方面為京津發展提供支持,是京津地區生態安全的綠色屏障和重要生態緩沖地帶,因此對京津及環京津六市進行生態系統服務評價和驅動因素分析,可以為未來協調該地區生態保護與社會經濟發展提供支撐,并更好地促進京津冀一體化發展[9—11]。

同時,為避免當前廣泛使用的“基于單位服務功能價格的方法”和“基于單位面積價值當量因子的方法”分別存在的忽略生態系統服務在時空尺度上的內在豐富信息,不利于對生態系統保護的管理決策[12—13]和可以轉化為市場交易價格的生態系統服務指標較少,限制了對區域關鍵生態系統服務的評估[14—15]的問題。本次研究在系統梳理研究區域的發展脈絡的基礎上,選取長時間跨度中的多個時間節點,并以當前京津及環京津六市面臨的生態挑戰作為研究指標構建符合區域特征的生態系統服務評價體系,生成生態系統服務綜合指數。通過PCA-OLS-GWR空間統計分析,較為科學地分析研究時期內,生態系統服務演變的特征,辨析其發展的動態驅動機制。通過分析實現主要研究目標如下：(1)分析2000—2018年間生態系統服務的時空演變特征；(2)構建自然環境因素和社會經濟驅動因子指標體系,闡釋各項驅動因子對生態系統服務演變的驅動機制；(3)結合政策導向提出區域生態系統服務優化提升路徑。

1 實驗數據與方法

1.1 研究區域與研究數據

本次研究將北京、天津以及環京津六市(承德、唐山、張家口、保定、滄州、廊坊)8個城市劃定為研究區域,總面積16萬km2。該區域地處東經113°40′—119°19′,北緯37°29′—42°37′,位于黃河下游以北,北靠燕山山脈,南朝華北平原,西倚太行山,東臨渤海灣,屬于典型溫帶季風氣候。地勢由西北向東南逐漸傾斜,降水量自東南向西北遞減,地貌類型和生態系統類型多樣,從西北向東南依次分為壩上高原區、燕山和太行山區、京津保中心區、海岸海域區以及低平原區五大區域(圖1)。2018年常住人口為7486.8萬人,2018年GDP總值高達68081.46億元,是中國北方地區經濟發展、人口聚集的城市群,在中國經濟社會發展中具有重要的戰略地位。

研究數據包括基礎地理信息數據、土地利用數據、遙感(NNP、NDVI)數據、數字高程(DEM)數據、氣象數據、土壤屬性數據、社會經濟統計數據7類。基礎地理信息數據及行政區劃、土地利用數據、溫度、降雨量等氣象數據和土壤屬性數據均來源于中國科學院資源環境科學數據中心(www.resdc.cn)；遙感數據和高程數據來自國家地理信息空間數據云平臺(http://www.gscloud.cn/)；人口等社會數據來源于2019年北京、天津及河北各市縣的統計年鑒；基礎交通設施POI矢量數據來自于世界數據地圖(www.openstreetmap.org),以上數據分辨率均為1km×1km。通過對比1998—2018年NDVI的平均值,根據突變點并結合研究區域的發展歷程篩選出具有代表性的4個時間點——2000年、2005年2015年和2018年,劃分3個時期進行研究。

1.2 研究方法

相關學者的研究表明,研究區域近年來面臨著植被減少、水土流失、水源匱乏、城市熱島、動物多樣性等一系列問題[16—19]。因此,本研究根據“千年生態系統評估”框架并依據當前研究區域生態系統的特征和面臨的生態壓力,從供給服務、調節服務、支持服務和文化服務4種服務類型出發,綜合凈初級生產力(NPP)、水土保持、固碳釋氧、水源涵養、熱島緩解、生物多樣性保護、游憩利用7種指標表征區域生態系統服務整體水平(表1)。通過ENVI遙感解譯、地理信息疊加分析等技術進行生態系統服務評估。

1.2.1凈初級生產力(NPP)

凈初級生產力(NPP)是綠色植物呼吸后所剩下的單位面積單位時間內所固定的能量或所生產的有機物質。公式如下：

NPP=GPP-Ra

式中,NPP單位為：g m-2a-1,Ra表示自養呼吸的消耗量,GPP表示總初級生產量。

表1 生態系統服務評價指標

1.2.2水土保持

水土保持是生態系統通過其結構與過程減少由于水蝕所導致的土壤侵蝕的作用。公式如下：

Spro=NPPmean×(1-K)×(1-Fslo)

式中,Spro為水土保持服務能力指數,NPPmean為植被凈初級生產力平均值,Fslo為坡度因子,K為土壤可蝕性因子。

1.2.3固碳釋氧

固碳釋氧是生態系統通過綠色植物的光合作用進行氧氣和二氧化碳的物質交換過程,生態空間每生產1g的植物干物質,可以固定1.63g的CO2,同時會釋放出1.2g的O2。公式如下：

Gv=1.63Rc×A×NPP；G0=1.19A×NPP

式中,Gv為植被年固碳量,單位為g/d,G0為生態空間年釋氧量,Rc為CO2中碳的含量,取值27.27%,A為生態空間面積(hm2),NPP為植被凈初級生產力平均值。

1.2.4水源涵養

水源涵養是生態系統通過其特有的結構與水相互作用,對降水進行截留、滲透、蓄積,并通過蒸散發,實現對水流、水循環的調控。公式如下：

WR= NPPmean×Fsic×Fpre×(1-Fslo)

式中,WR為生態空間水源涵養服務能力指數, NPPmean為植被凈初級生產力平均值,Fsic為土壤滲透因子,Fpre為平均降水量因子,Fslo為坡度因子。

1.2.5熱島緩解

熱島緩解是通過增加植被、水體及濕地等生態系統,減輕由于人為原因導致的城市溫度高于周邊地區的熱島效應。公式如下：

Ts=Ti+A(Ti-Tj)×B

式中,Ts為地表溫度,單位為℃/a,Ti和Tj為熱通道i和j的亮度溫度,A和B為系數,由大氣透過率和地表比輻射率等因子決定。

1.2.6生物多樣性維護

生物多樣性維護是生態系統支撐服務的重要內容之一,旨在維持物種、基因多樣性中發揮的作用。通過該項指標的評價,識別現狀和未來不同區域承擔生物多樣性保護的能力強弱。公式如下：

Sblo=NPPmean×Fpre×Ftem×(1-Falt)

式中,Sblo為生物多樣性維護服務能力指數,NPPmean為植被凈初級生產力平均值,Fpre為平均降雨量,Ftem為平均氣溫,Falt為海拔因子。

1.2.7游憩利用

游憩娛樂為人類生活提供了休閑、娛樂與美學享受。通過對游憩資源、游憩設施、游憩區位三方面多個指標進行評價,判斷生態系統文化服務的能力水平[12]。公式如下:

式中,S為游憩利用綜合評價指數,Wi為游憩因子等級值,Xi為不同游憩因子權重值,n為評價因子的個數。

1.3 生態系統服務綜合指數

為避免同一類型服務因指標因子多而造成數值過高的問題,本次計算首先對上述7項服務指標進行最大值-最小值歸一化,然后按照生態系統服務4大分類等權疊加后再次進行歸一化,最終將四項服務指標累加獲取生態系統服務綜合指數,表征區域生態系統服務的總體水平。具體公式如下：

式中,yij、xij為生態系統服務因子i的標準化數據和原始數據;xmax和xmin分別為狀態指數j的原始數據的最大值和最小值。

CIES為生態系統服務綜合指數;ES為第i類生態系統服務的標準化結果;n為4。

1.4 因子篩選與相關性分析

1.4.1普通最小二乘法(Ordinary Least Squares,OLS)

OLS模型作為全局回歸模型,常用于解釋單一因變量與多個自變量之間的驅動關聯。本研究將其作為篩選影響城市群生態系統服務相關因子的評價途徑,計算生態系統服務價值與各項驅動因子基礎回歸關系,并進行因子篩選。計算公式如下：

yi=β0+∑iβixi+εi

式中,β0是常數項；βi為回歸系數；εi為隨機誤差項

1.4.2地理加權回歸(Geographical Weighted Regression, GWR)

在建立OLS模型的基礎上,通過GWR模型建立嵌入空間加權函數的線性回歸模型,用以探索研究對象在某一尺度下的空間變化關系及驅動因素。相比使用OLS模型進行回歸分析,GWR模型能夠反映參數在不同空間的非平穩性,使變量間關系可以隨空間位置變化,同時考慮空間對象的局部效應,其結果具有更高的準確性[20—21]。其評價公式如下：

βk(Ui,Vi)=(XTW(Ui,Vi)X)-1XTW(Ui,Vi)y

式中,β0(Ui,Vi)是(Ui,Vi)空間位置的地理加權回歸截距,βk(Ui,Vi)是k項自變量(驅動因子)在(Ui,Vi)空間位置的加權回歸系數,xk(Ui,Vi)為第k項自變量(驅動因子)在(Ui,Vi)空間位置的取值,εi為算法殘差,XT是自變量(驅動因子)的轉置,W(Ui,Vi)是距離權重矩陣,h為AIC準則的帶寬,為空間i位置與空間j位置的距離。

2 研究結果

2.1 生態系統服務時空演變特征

對研究區域生態系統服務指數進行具體分析發現,4類服務指數在各個時期表現出的特征各異(圖2),具體情況如下：1)供給服務指數在2000—2018年間總體上呈現從上升向平穩過渡的特征,并在2015年達到研究時期最大值0.316。2)調節服務指數在2000—2018年間呈現波動逐漸增強的變化特征,在2015年達到研究時期最高值0.333后并在2018年驟降至0.292。3)文化服務指數在2000—2018年間在總體上呈現先上升后下降末期數值高于初期的特征,在2005年達到研究時期最高值0.057。4)支持服務指數在2000—2018年間變化幅度較大,總體呈現先上升后下降再上升的特征。在2005年達到研究時期最高值0.057并在2015年降至研究時期最低值0.020,整體波動較大。

通過對研究區域在2000年、2005年、2015年、2018年生態系統服務綜合指數進行綜合量化分析,發現研究區域生態系統服務綜合指數呈現出先上升后下降、研究期末生態系統服務綜合指數略高于研究初期的總體特征。

圖2 2000—2018年生態系統服務指數變化圖Fig.2 Map of changes in ecosystem service index from 2000 to 2018

通過對4類生態系統服務指數在空間變化上進行研究后發現(圖3)：1)供給服務指數在研究期間內呈現以太行山燕山山脈為界東南、西北兩側分別以“增強-減弱-增強”和“減弱-增強-減弱”的間隔規律差異發展的總體特征。2)調節服務指數在京津保中心區和燕山區域內逐漸增強,但在壩上高原和西部太行山區呈現相反演變特征,而海岸海域區、低平原區以2015年為界,先增強后降低。3)文化服務指數在研究期間,以2005年為界,燕山和太行山區、京津保中心區、低平原區呈現先增強后減弱變化趨勢。海岸海域區則呈現相反的態勢。壩上高原區在研究期間變化趨勢不明顯。4)支持服務指數,除了在太行山區呈現先增強后減弱再增強的變化波動現象外,整體變化不明顯。

圖3 生態系統服務指數空間分布圖Fig.3 Spatial distribution of ecosystem service index

生態系統服務綜合指數的分級方法是在地理信息系統軟件中將4類服務做疊加處理,得到生態系統服務綜合指數空間分布。同時,進一步采用自然間斷點分類對其進行5級分類。通過對在研究期間內的生態系統服務綜合指數進行評估得出(圖4)：研究期間內,總體上燕山山脈及太行山山脈北部區域為研究區域生態系統服務綜合指數的高值區域,壩上區域為低值區。生態系統服務綜合指數波動主要集中在太行山和燕山區,燕山區域生態系統服務綜合指數在2000—2015年間呈下降的趨勢,2015—2018年間呈上升的趨勢,而太行山區域在整個研究時期呈現下降趨勢。

圖4 生態系統服務綜合指數空間分布及變化Fig.4 Spatial distribution and change of ecosystem service index

2.2 構建驅動模型

為了探索自然社會因素與生態系統服務之間的關系,本研究在可獲取的數據范圍內,結合研究區域實際情況[22—24]和相關研究成果[25—28],最終選擇涵蓋氣候、地形土壤、人類活動及土地利用5類共8項因子代表區域自然地理特征和社會經濟發展的特征(表2)。這些因子作為研究區域生態系統服務的預選驅動變量,并確定最鄰近狀態的2018年進行驅動要素分析評估。為避免自變量出現明顯的共線性問題,導致因子對模型的解釋力出現偏差,研究以擬合優度大于0.5、P值邊界0.05、冗余檢驗小于7.5作為篩選條件,運用OLS模型篩選出5個驅動因子,即降雨量、溫度、土壤含沙量、人口密度和基礎設施密度。在對比OLS模型與GWR模型的基礎上(表3),將驅動因子導入GWR模型中進行地理空間加權回歸分析。模擬結果在ArcGIS 10.6軟件中進行空間可視化,獲得的GWR模型擬合結果顯示圖(圖5)。在擬合優度(R2)模擬結果上,除保定南部極少部分區域外整體擬合優度均高于0.6,平均值為0.655,模型效率可信；在擬合t值模擬結果上,整體區域模型t值范圍均在-2.58—2.58內,模型結構可信。

圖5 地理加權回歸模型擬合常量R2和t值分布圖Fig.5 Geographically weighted regression model fitting constants R2 and t-value distribution

表2 驅動因子統計分析表

2.3 驅動結果分析

全局回歸結果表明(表4),各因子的驅動強度排序如下：降雨量(5.244405)>土壤含沙(-2.185510)>人口密度(-1.009130)>溫度(0.788722)>基礎設施密度(-0.649611)。由此可見,與社會經濟因素相比,自然環境因素對研究區域生態系統服務能力影響更大。GWR回歸結果如圖6所示。在自然驅動因子中,降雨量對生態系統服務的驅動影響呈現明顯的正相關性,驅動作用強度從承德西北部和張家口東南部山地向南呈現遞減的趨勢。溫度對生態系統服務驅動影響呈現由北向南從正相關向負相關過渡的趨勢,正向驅動作用在保定及廊坊以北逐漸增強,在承德則出現最為明顯的正驅動作用,而在以南地區呈現負相關。土壤含沙量與研究區域生態系統服務整體呈負相關性,土壤含沙量的增加會給生態系統服務帶來負面影響,且從北向南逐漸增強。南部地區保定、滄州市較中北部城市呈現的負驅動作用最為明顯。人口密度和基礎設施密度2項社會驅動因素都對生態系統服務驅動呈現負相關。人口的增長會給生態系統服務能力帶來負面影響,且從東北向西南逐漸增強,南部地區保定、滄州市較中北部城市呈現出更為明顯的負驅動作用。基礎設施密度對生態系統服務的負向驅動作用由太行山脈部向東南和東北方向逐漸增強,以壩上高原區及海岸海域區負向驅動作用最為明顯。

表3 基于探索性回歸的最小二乘法(OLS)模型參數因子篩選結果表

表4 全局回歸結果的計算結果

圖6 各因子與生態系統服務回歸系數空間分布圖Fig.6 Spatial distribution diagram of regression coefficients between various factors and ecosystem service

3 討論

3.1 自然-社會因子對生態系統服務的驅動作用分析

人為因素和自然因素是生態系統變化的主要驅動力因素,在長時間跨度上,自然因素對生態系統服務的影響是巨大的,而在短時間內,生態系統服務則更易受人為因素影響[29]。京津及環京津六市是京津冀城市群中人口最為稠密、土地利用方式復雜、經濟生態社會脅迫效應顯著的區域。在研究期間內,該區域的生態系統服務呈現明顯且明確的空間分異,即在空間維度上,生態系統服務的高值區域集中在燕山山脈區域,而生態系統服務低值區域集中在壩上高原地區。本研究運用“OLS降維模型+GWR回歸模型”探索了驅動因素對生態系統服務的空間異質性影響[30—31],使得研究結果更加全面且符合區域實際。GWR回歸結果表明,由于研究區域位于屬于中國半干旱向半濕潤過渡的地帶北方農牧交錯帶,水資源較為匱乏,生態本底較差[32]。因此,降雨量對該區域生態系統服務的正向演化具有積極作用,并對較為干旱的壩上高原區域表現出更強的驅動力。同時溫度也是限制植物生長的主要因素[33],由于太行山區及壩上高原區氣候寒冷,而南部平原地區較為溫暖,所以溫度對生態系統服務的正向驅動作用由南向北增強。其次,土壤的顆粒組成及含沙量是影響土壤可蝕性差異的主要因素。土壤中其含沙量越大,土壤可蝕性越高,易發生侵蝕[34]。而在研究區域內,南部地區用地性質以農田為主,植被單一加之人類不合理的生產活動,導致土壤養分流失和土地退化以及水土流失等生態問題,從而影響生態系統的穩定性[35—37]。社會因子對生態系統服務的影響也較為明顯,人口密度和人類活動強度的增加對生態系統產生較強的干擾和破壞作用[38—39]。人口增長會增加對周邊自然資源需求和改變土地覆蓋變化,對周邊生態環境造成不利影響[31]。本研究區域北部山地城市人口密度較小,基礎設施建設量也相對較少,對自然系統的破壞較小。而中部南部平原城市由于北京首都功能疏解導致的人口外溢和人口的高流動性,使該區域人口聚集,對生態系統服務的負驅動作用更為顯著。由此可見,人類活動多集中在易于建設的平原地區,社會因素對平原的生態系統服務負向驅動力明顯高于山區與高原地區。

3.2 生態系統服務時空演變分析

本研究通過總結研究區域當前面臨的關鍵挑戰,選取確定了凈初級生產力(NPP)、水土保持、固碳釋氧、水源涵養、熱島緩解、生物多樣性保護、游憩利用7種評價指標,生成該區域生態系統服務綜合指數,與相關學者對該區域的生態系統服務時空分布的研究結果相似[28],但由于在指標選取方面存在差異,生態系統服務的空間分布也存在一定的差異,同時相比與運用生態系統服務價值時空分布,該實驗結果可以更加準確地反應時空尺度上的內在豐富信息。通過對2018年研究區域生態系統服務驅動因子的分析,結合研究區2000—2018年生態系統服務指標時空演變的特征,對近20年的生態系統服務的驅動因子作用機制進行分析：2000—2015年間,退耕還林等生態治理工程、京津風沙源治理工程等生態治理工程在研究區域內不斷推進。這些生態保護項目注重地表覆蓋的綠化,植被覆蓋的明顯提升,能夠降低土壤侵蝕,促進土壤保持服務、固碳等服務能力的提升[40—41]。2015年,中共中央政治局召開會議,審議通過《京津冀協同發展規劃綱要》,為京津冀地區注入了新的發展活力。在這4年間,中國政府開始大規模推進京津冀城市群一體化發展戰略,開展一系列的基礎設施建設活動,改變了區域土地利用的方式,造成了綠地面積的減少和生態空間的破碎化。這無疑加大對生態資源的消耗,對當地的降雨、溫度等自然環境因子產生了影響,造成生態系統服務能力的退化[42—43],也使得研究區域生態系統服務總體則呈現下降的趨勢。因此在未來該區域發展的過程中,應將與人口密度、基礎設施建設相關的政策作為必要的輔助行動,來維護區域的生態系統健康。

3.3 生態資源與空間管控優化策略

社會與生態是相輔相成的,二者只有相互協調達到良好的耦合狀態,才能實現研究區域生態保護和社會發展的雙重目標。因此,未來應該依據分區規劃,因地制宜的采取有針對性的生態資源與空間管控優化策略,并積極推動生態補償體系建設,以實現區域協調發展[44—45]。在京津保中心區,通過疏解非必要的首都功能,控制人口聚集、城市擴張來減弱對周邊生態系統的影響,并適當增加生態用地比重。在燕山和太行山區內應繼續開展一系列的生態修復工程。在保護生物多樣性的前提下,可以適當促進文旅產業發展,增加當地收入。在壩上高原區,受到降雨等影響較大,且生物多樣性低,應通過增加植被覆蓋、保護河流,增強該區域的水源涵養和地表蓄水能力,提升區域生物多樣性[46]。在低平原區應關注土壤沙化問題,合理開發水資源,保護生態系統免受退化[47]。

4 結論

由于城市化進程的加快,人類活動極大地改變了城市群的生態系統服務能力。本研究選擇以2000—2018年近20年的長時間跨度,探索生態系統對研究區域的功能反饋效益。并探討了生態-社會驅動因子與生態系統服務的關系,揭示其在空間尺度上對生態系統服務的影響程度與作用機制。結果表明：1)2000—2018年,研究區域的生態系統服務綜合指數呈現出先上升后下降,在2015年達到最高值。2)2000—2018年生態系統服務綜合指數呈現明顯的空間分異特征,以燕山和太行山為界,燕山和太行山區域呈現先減弱后提升的趨勢,西北區域則整體呈現減弱趨勢,東南區域從2015年起逐漸降低。3)自然-社會因子對生態系統服務的驅動作用從強到弱排列依次是：降雨量>土壤含沙>人口>溫度>基礎設施4)。在驅動因子中,分析結果表明：人口密度、土壤含沙量、基礎設施密度呈較為明顯的負相關驅動特征,降雨量呈較為明顯的正相關驅動特征。溫度呈南部負相關,北部正相關。同時,本文還依據實驗結果提出了相應的優化策略,為研究區域實現生態系統可持續管控提供有價值的政策指引。