基于格網的南四湖流域土地利用碳排放與其生態系統服務價值時空關系分析

王 丹,荊延德,*,韓善梅,高明秀

1 曲阜師范大學地理與旅游學院, 日照 276826 2 日照市國土空間規劃與生態建設重點實驗室, 日照 276826 3 山東農業大學資源與環境學院, 泰安 271018 4 土肥資源高效利用國家工程實驗室, 泰安 271018

生態系統服務定義為：由生態系統提供的,人類從中收獲的福祉和利益[1]。通常使用經濟價值來評價生態系統服務,以確定我們從自然中獲得的利益,使不同區域具備可比性。在全球變暖的現狀下,碳減排已經變為一個全球性的緊要話題[2]。兩者因土地利用類型與陸地表面特征和功能而有很大區別。經濟發展通過頻繁改變土地利用方式間接影響陸地生態系統和人為的碳排放。土地利用變化對碳排放不只是簡單的線性關系,而是更為深遠的影響。同時土地利用方式的改變也極大地影響了自然生態系統,從而影響了它們所提供的服務。因此土地利用變化不僅影響著碳排放的空間分布、范圍和強度,而且還通過改變生物多樣性、區域資源和生態系統類型的空間分布間接影響ESV[3]。隨著工業化的快速發展,中國成為全球碳排放量最高的國家[4]。特別是1950—2005年因土地利用變化引起的累積碳排放量為106×108t,是同期世界土地利用變化碳排放量的12%[5]。為了減緩碳排放和全球變暖,中國制定了不同的環境政策并采取了有效措施[6]。特別是,2020年習近平總書記在第七十五屆聯合國大會一般性辯論上提出了碳達峰、碳中和的目標。2021年國家林草局為落實這一決定提出“十四五”期間將扎實開展林業和草原碳匯行動,提升林草生態系統碳匯能力。因此,在碳減排政策和改善生態環境的雙重背景下,研究兩者間的空間關系模式,對解決社會發展過程中既要保持經濟增長,又要兼顧生態環境這一難題有重大意義。

自20個世紀90年代開始,一方面,國內外學者圍繞土地利用碳排放問題從國家[7]、省域[8—9]、城市[10]、區域[11—12]等不同尺度對碳排放與經濟增長的脫鉤效應分析[13—15]、單一地類碳排放機理剖釋[16—17]、影響因素[7,18—19]、變化規律[10,20]、效率[21]、地類變化與碳源/匯的關系[11]、碳排放空間差異與關聯[22]等角度展開了大量的研究。另一方面,在評估土地利用變化造成ESV波動方面,學者們還進行了許多研究[23—26],如ESV和土地利用的時空變化特征[27—28]、相互作用關系[25,29—30]等。綜上可知,以往都是從不同角度和尺度對兩者單獨開展了研究,但是關于土地利用碳排放和ESV之間空間關系的研究較少[31—32],并且南四湖流域還未開展土地利用碳排放的相關研究。

南四湖流域位于魯南經濟圈,是我國北方最大的淡水湖泊,肩負重要糧食生產的重任,同時在京杭大運河和國家戰略性工程中發揮重要作用。未來如何推動流域內各縣市在經濟持續發展的前提下,實現節能減排、保護生態環境是值得思考的問題。鑒于此,本文以南四湖流域為例,利用土地利用數據和統計資料,對2000—2018年間的土地利用轉移中碳排放量和ESV的時空演變開展研究,并基于格網對碳排放強度和生態系統服務強度進行空間關聯分析,以期為流域提出碳減排政策,建立低碳土地利用結構,努力改善生態環境總體質量,從而走向綠色低碳發展的道路提供支撐。

1 研究區概況

圖1 研究區地理位置Fig.1 Location of the study area

南四湖流域處于中國東部淮河流域的北部,涉及山東、江蘇兩省5市,面積約3.17×104km2(圖1)。地處魯中南山區西側的山前沖擊洪積平原,中、西部為平坦的平原,東部為山地丘陵地形。海拔在- 42—630 m之間,絕大部分地區低于70 m。耕地分布廣泛,占全區面積的68.5%以上。自21世紀以來由于城鎮化進程的推進,流域國內生產總值由1174.91×108元變為12078.12×108元,翻了9.28倍；建設用地面積從4635.15 km2增長到6035.81 km2,增長了30.22%。近年來,該區受社會發展的直接影響,土地利用方式劇烈變化,能源消耗需求增加,對當地的生態和碳排放系統帶來一系列問題。同時由于其所處地理位置及其作用,各級政府和相關部門都十分重視生態保護,啟動了《南四湖流域水污染綜合整治三年行動方案(2021—2023年)》等一系列工作,以防止環境惡化。

2 數據來源與研究方法

2.1 數據來源及預處理

本文使用地理國情監測云平臺的2000—2015年的土地利用數據,其分辨率為30 m。2018年土地利用數據來源于中國科學院資源環境科學數據中心,是通過人工目視解譯Landsat 8遙感影像獲得的。依據中科院土地利用分類標準,將研究區的地類分為6個一級類型,依次為耕地、林地、草地、水體、建設用地和未利用地。能源消耗和社會經濟數據分別來源于對應年份的《中國能源統計年鑒》、《統計年鑒》和《全國農產品成本收益資料匯編》等統計資料。利用網格法按照5 km×5 km大小的網格對研究區進行劃分,共得到1387個評價單元。由于數據可獲取性,本文中2000—2018年南四湖流域能源消耗數據,根據研究區內各區縣GDP與所在省份GDP的比值來進行推算。

2.2 研究方法

2.2.1 土地利用類型轉移矩陣

土地利用轉移矩陣以二維矩陣的形式反映19年間各用地類型的數量轉移情況,清楚掌握用地類型在2000年的流動方向及2018年的來源構成。其公式為：

(1)

式中：Sij為土地面積；n為6；i、j分別為2000年與2018年的地類序號。

2.2.2 碳排放量的估算

研究區碳排放量計算：即各用地類型的面積和相應的碳排放系數相乘再求和。碳排放計算公式為：

Eκ=∑ei=∑Si×δi

(2)

式中：EΚ表示直接碳排放量；ei表示各用地類型的碳排放量；Si和δi分別表示用地類型i的面積和碳排放系數,參考相關文獻[33—34]和研究區的現實情況確定各用地類型的碳排放系數,見表1所示。

表1 用地類型碳排放系數

基于2000—2018年南四湖流域化石能源消耗量計算建設用地碳排量,然后除以其面積得到建設用地碳排放強度。建設用地碳排放量的計算采取間接估算法,結合《中國能源統計年鑒》和2006年IPCC《國家溫室氣體排放清單指南》(表2)[38],通過各種能源的消耗量和碳排放系數估算得到。

表2 各類能源折標準煤系數和碳排放系數

建設用地碳排放計算公式：

Eη=∑ei=∑Ei×μi×εi

(3)

式中,E為建設用地碳排放量；ei為能源i消耗產生的碳排放量；Ei為能源i的消耗量；μi為能源i消耗量轉換為標準煤的系數；εi為能源i的碳排放系數。

土地利用碳排放強度由格網內各用地類型的面積和碳排放系數決定,即碳排放強度越大說明單元格網的碳排放量越大。計算公式如下：

(4)

式中,CRI為土地利用碳排放強度；Si和Pi分別為格網內用地類型i的面積和碳排放系數；S為格網總面積。

2.2.3 生態系統服務價值

參照謝高地等[39—40]的研究,1個生態系統服務價值當量的價值量由每公頃農田年平均自然糧食產量的經濟價值1/7[40—41]來確定,已廣泛在流域[39—40]等尺度開展深入研究。為避免各年份和不同省份糧食價格的差異對價值計算的影響,本文僅以山東省2018年小麥、玉米和稻谷3種主要農作物的播種面積、單產和糧食價格為基礎數據,計算出南四湖流域ESV的單位當量為2085.53 元,其計算模型如下：

(5)

式中,Ea為南四湖流域ESV當量因子的經濟量；n為3；pi為農作物i在山東省的價格；qi為農作物i單產；mi為農作物i面積；M為農作物總面積。

由于我國疆域遼闊,各地的生態系統各具特點,為降低直接使用2015年謝高地完善的“單位面積生態系統服務價值當量表”[39]產生的誤差,本文參考研究區已有研究[42]和現實情況修正了各地類的生態系統服務價值系數(表3),其中不考慮建設用地的ESV。

表3 南四湖流域各地類生態系統服務價值系數/(元 hm-2 a-1)

生態系統服務價值：

(6)

(7)

生態系統服務價值強度：

(8)

式中,Aij為網格j中地類i的面積；Ci為地類j的價值系數；S為網格面積；n為6種土地利用類型,m為1387。

2.2.4 雙變量空間自相關性分析

我們使用雙變量Moran′sI研究了碳排放強度與ESV強度之間的空間聚類和離散。本文使用了兩種類型的雙變量Moran′sI,即全局雙變量Moran′sI和局部雙變量Moran′sI。全局雙變量Moran′sI研究了整個區域內碳排放強度與ESV強度之間是否存在以及空間相關性程度,而局部雙變量Moran′sI則在具體的空間格網上顯示空間相關性。計算公式如下：

(9)

(10)

3 結果與分析

3.1 土地利用變化分析

如表4所示,2000—2018年南四湖流域各土地利用類型間均進行了復雜的轉移,發生轉移的總面積為4468.86 km2,是流域總面積的14.10%。耕地對轉移面積的貢獻最大,為2579.12 km2,其中有2084.55 km2變為建設用地,占比80.82%,15.72%的轉出量轉為水體,雖占比較低,但轉出面積大,為405.40 km2；其次是建設用地,轉出量為854.38 km2,95.87%的土地轉為耕地；其它大部分地類的轉出量轉為耕地。出現此情況的原因是流域經濟快速發展促使城鎮外圍的耕地被大量占用,但為保證糧食安全,故改變其它生態用地的使用方式,從而確保耕地總量,之后又由于實施“兩退三還”(退耕、退池、還林、還濕、還湖)工程和居民點整治等原因導致土地轉換復雜且多樣。建設用地是轉入最多的土地類型,轉入2255.05 km2,其次是轉入面積為1581.84 km2的耕地。其他地類均發生轉入轉出,情況各異。

表4 南四湖流域土地利用轉移情況

3.2 土地利用碳排放分析

3.2.1 南四湖流域2000—2018年土地利用碳排放測算與分析

研究區土地利用碳排放量的改變受地類和土地利用程度的影響。本節主要以各地類的碳排放量、碳源/匯和凈碳排放量等指標測算南四湖流域土地利用碳排放,分析2000—2018年間各指標的動態變化特征。利用上述計算公式和相關收集數據計算南四湖流域在研究期間的相應碳排放指標,結果見表5。

通過計算可知,南四湖流域的凈碳排放量表現為碳排放(表5)。2000—2018年碳排放能力顯現穩定增長趨勢,由1262.51×104t增長至6438.41×104t,年均增長率為22.78%,變化趨勢與碳源量相似。其中,建設用地碳排放量增長迅速,從2000年的1176.36×104t增長至2018年的6355.35×104t,但增長量具有“先快速,后放緩”的特點,這是由于流域內城鎮化水平得到大幅度提高,伴隨著大量的能源消耗和建設用地增加,之后由于能源結構的優化和節能減排的實施；耕地的碳排放量由于本身面積的減少呈現緩慢下降的趨勢,由2000年的96.02×104t下降至2018年的91.81×104t,年均下降率為0.24%。碳匯量在19年間呈現先增加后減小并趨于穩定的態勢。其中林地和草地的碳吸收量在2005年后出現下降,水體基本保持不變,未利用地波動較為明顯,呈先減少后增加的特點。

從碳源/匯構成上分析,碳源方面,建設用地的碳排放是最主要碳源,占總碳源量的92.45%—98.58%,并呈逐年增長態勢。同時耕地作為碳源之一,碳排放貢獻率低,但面積占比高。地均建設用地碳排放量表示各年建設用地的碳排放強度以及區域的能源消費結構,自2000年以來表現為快速增長趨勢,各年的地均建設用地碳排放量依次為25.38、62.16、80.28、96.02、105.29 t/hm2,其中2018是2000年的4.15倍,因此建設用地在碳減排方面具有較大潛力；碳匯方面,林地和水體對碳吸收的貢獻率相當,兩者之和占總碳匯的96.11%以上,說明林地和水體具有重要的碳匯能力,而在流域內兩者占總面積的5.89%—8%,凸顯了其重要的碳吸收能力,另外草地和未利用地的碳吸收量較小。由于碳源量與碳匯量相差懸殊,因此該流域今后的工作重心是碳源的控制與碳匯的增加。

表5 南四湖流域土地利用碳排放量核算結果

3.2.2 土地利用碳排放強度的時空特征

由于流域內地類結構和發展水平的差異,因此為了具有可比性以及直觀地分析土地利用碳排放量的變化情況,故本文利用碳排放強度指標,依次分為Ⅰ、Ⅱ、Ⅲ、Ⅳ和Ⅴ 5個區間(見圖2)。19年間南四湖流域碳排放強度時空變化明顯,其特點如下：流域碳排放強度低值多分布于東部,整體上呈現增加的特點,碳排放強度最大值從2000年的21.61 t/hm2增長到2018年的101.42 t/hm2,增長了4.69倍。落在Ⅰ區間的格網數逐年減少,尤其是2000—2005年間減少了78.3%,出現了向Ⅱ區間轉移的現象,而同時南四湖湖區始終屬于Ⅰ區間；落在Ⅱ區間的格網呈現先增加后減少的特征,2010年以來土地利用類型主要為林地和草地的東部和南四湖湖區的周邊穩定的處于該區間；Ⅲ區間的格網數的增長率最高,在2018年占總格網數的46.36%,絕大部分與耕地的分布范圍相吻合；Ⅳ和Ⅴ區間的格網從無到有并逐年增加,呈團塊狀,主要集中在人類活動活躍的建設用地。

圖2 2000—2018年土地利用碳排放強度空間演變Fig.2 Spatial patterns of carbon emission intensity on land use during 2000—2018

3.3 生態系統服務價值的時空演變

2000—2018年南四湖流域生態系統服務總價值呈波動變化,總體上呈增加趨勢,總價值704.07×108元增加到743.37×108元,增加了5.58%(表6)。從各地類的ESV來看,耕地和水體是總ESV的貢獻主體,依次占24%和61%以上。19年間水體的ESV占比穩定增長,由61.85%到67.79%；而耕地的ESV占比逐年下降,由27.03%到24.48%；其它地類的ESV總體上呈下降態勢。從土地利用構成中得到,19年間建設用地的面積大幅度擴張,林地、耕地、草地和未利用地均發生不同程度的減少,而水體面積的增加,反而使總ESV增加29.30×108元,說明水體的ESV強度最高,具有極其重要的生態地位,同時各地類間的相互轉換是導致研究期間內總ESV波動變化的關鍵原因。

表6 2000—2018年南四湖流域各土地利用類型ESV統計

為了直觀地分析流域ESV的空間分布特征,同時避免流域邊緣面積的影響,本節利用ESV強度表示各格網,并對其重分類依次分為Ⅰ、Ⅱ、Ⅲ、Ⅳ 4個區間(見圖3)。南四湖流域的ESV強度具有明顯的空間差異,西部地區的ESV強度整體高于東部地區。其中,處于Ⅰ區間的ESV強度的分布范圍廣,其格網占總格網的72.45%以上,主要由于耕地的空間分布相對均一且占比高；處于Ⅳ區間的格網占總格網的4%左右,土地利用類型為水體,主要分布在南四湖湖區。與其他年份相比,2018年處于Ⅰ區間的格網變少,出現向Ⅱ區間轉變的現象,同時Ⅱ區間的分布范圍主要位于流域的東部,集中分布著草地和林地；Ⅳ區間的分布范圍由原先的條帶狀出現斷裂；Ⅲ區間的分布范圍變化不明顯。

圖3 2000—2018年南四湖流域生態系統服務價值時空變化Fig.3 Spatio-temporal changes of ecosystem service value in Nansi Lake Basin during 2000—2018

3.4 碳排放強度和ESV強度的關系

為了深入地了解碳排放強度與ESV強度的相互交互作用,首先利用SPSS軟件得到2000—2018年兩者的相關系數依次為-0.453、-0.413、-0.347、-0.335和-0.341；然后借助空間自相關揭示兩者的空間關聯規律(表7),各年份的Moran′sI均為負,P值小于0.001,表明研究區內土地利用碳排放強度與ESV強度之間具有空間負相關性,即隨著碳排放強度的增強,ESV強度呈下降的態勢,并且Moran′sI具有逐年下降的特點,表現出整體較為分散、小范圍聚集的空間布局格局。

表7 2000—2018年南四湖流域雙變量全局Moran′s I統計表

通過LISA集聚圖可以看出(圖4),2000—2018年間,低高聚集區穩定地分布在南四湖湖區,表示該區碳排放強度低值與ESV強度高值形成集聚現象；2000和2005年的高低聚集區集中分布在耕地面積占比高,建設用地點狀零散分布的西部和北部,而在2005年之后出現了大面積向低低聚集區轉變的現象,高低聚集區與呈片狀分布的主要城鎮用地相吻合,主要是由于2005—2010年間地類之間的變化劇烈,總轉移面積為3581.62 km2,其中有1687.82 km2的其他類型土地轉化為建設用地,說明在該時間段流域的城鎮化水平得到快速發展,能源消耗增加,導致局部地區的碳排放強度發生大幅度的上升；雖然建設用地和耕地的價值系數低且都屬于碳源,但是前者的碳排放系數要遠高于耕地,從而使耕地面積占比高的格網與周圍的低ESV強度形成集聚現象；低低聚集區在2000—2005年零散分布且變化不大,之后分布范圍擴大,主要分布在西部和北部,地類為具有面積優勢的耕地；高高聚集區呈零星分布,這是由于周圍是水體,屬于碳匯。

圖4 2000—2018年雙變量空間聯系的局部指標(LISA)集聚圖Fig.4 Bivariate local indicators of spatial association(LISA) cluster graph during 2000—2018

4 討論

以往的研究大多集中在獨立的行政單元(如縣、市等)的經濟和社會等方面的碳排放效應分析[7—9,18]。本文將跨越行政區的流域作為一個研究整體,基于格網探究因土地利用變化引起的碳排放和ESV以及兩者的空間聚集規律。研究發現流域的建設用地面積和凈碳排放量在這19年內呈增長的趨勢,與Chen等[43]以成渝城市群為研究區、李玉玲等[14]以陜西省為研究區的研究結果一致,說明土地利用的改變是影響碳排放量和碳源/匯格局的關鍵原因,不當的土地開發會影響碳平衡；同時與中原城市群[31]的空間關系一樣,碳排放強度與ESV強度之間具有空間負相關性,碳排放強度對ESV強度有空間溢出效應。結合研究結果,通過在流域內實施能源消費結構的升級、提高利用效率、引進低碳生產技術、防范耕地被濫用和提高建設用地集約度,可以實現碳排放的降低和有效保護生態系統。

盡管本研究展示了19年間碳排放和ESV的時空特征,但仍存在一些潛在的局限性需進一步改進。首先,本文忽略了空間異質性,只是借鑒前人研究和資料得到各地類地碳排放/吸收系數[44],如耕地碳排放受種植作物種類、土壤的性質和耕作特點等多種因素的影響而有所差異。本文與國內絕大部分研究類似將耕地視為碳源,而馬濤[45]的研究表明上海的農業活動表現為碳匯。各能源碳排放系數參考了2006年IPCC國家溫室氣體清單指南,這可能并不完全適用于研究區。因此,今后根據區域的實際情況科學合理地確定碳排放系數,或者形成統一的測算方法；其次,本文使用了高精度的土地利用數據,但基于5 km×5 km網格對流域進行劃分。由于生態和地理過程伴隨著復雜的尺度效應[46],不同的尺度研究的側重點不同。因此未來有必要從多尺度探究碳排放強度與ESV強度關系特征；最后,由于受數據收集的限制,本文通過GDP占比推算流域的能源消耗,同時沒有將工業生產、生活和廢棄物排放加入建設用地的碳排放核算中,所以計算結果有誤差,但不影響區域碳排放的時空分析。因此,今后需要尋找全面且受數據限制小的計算方法。

5 結論

(1)研究區具有面積優勢的地類是耕地和建設用地,分別占總面積的68.5%和14.5%以上。在19年內各地類之間進行了不同程度的轉移,其中耕地和建設用地是變化量最大的類型。產生該變化的原因與城鎮化、人類生產生活和退耕還湖等相關政策有關。

(2)19年內流域的總價值隨地類之間的轉化呈波動變化,整體上是增加,水體面積的增加是導致其的決定性原因。通過ESV強度圖可知,東部的ESV強度明顯高于西部,同時Ⅲ、Ⅳ區間分布在南四湖湖區,Ⅱ區間分布在海拔較高的東部,這與地類、價值系數和地形有關。

(3)南四湖流域的凈碳排放量表現為碳排放,且穩定增長。其中,碳源量要遠大于碳匯量,地均建設用地碳排放量19年間增長了4.15倍,這與經濟社會發展有密切的關聯。通過碳排放強度的時空演變分析可知,Ⅰ區間的格網轉變為Ⅱ和Ⅲ區間的格網,后者的分布范圍迅速擴張,同時Ⅳ和Ⅴ區間的格網從無到有且逐年增加,與主要建制鎮的位置相吻合。整體上呈現東部高、西部低的特點,與ESV強度的空間分布有一定的關聯。

(4)對南四湖流域土地利用碳排放強度和ESV強度進行數量和空間統計分析,得知兩者具有負相關性,且通過了P值檢驗。19年內兩者的空間關聯規律發生較大的變化。