縣域土地利用變化下生態系統服務價值損益預測——以重慶市巴南區為例

馮 源,朱建華,2,曾立雄,2,肖文發,2,*

1 中國林業科學研究院森林生態環境與保護研究所 國家林業和草原局森林生態環境重點實驗室, 北京 100091 2 南京林業大學南方現代林業協同創新中心, 南京 210037

生態系統服務指人類從生態系統結構、過程和功能直接或間接得到的產品或惠益[1]。以貨幣形式評估的生態系統服務價值(Ecological services value, ESV)是生態功能區劃、環境保護及生態補償的重要指標[2- 3]。土地利用/覆被變化(Land use and land use change, LUCC)是全球變化的重要驅動力,通過影響生態系統格局與過程,改變生態系統服務的提供能力,最終導致生態系統服務價值發生變化[4- 6]。量化土地利用變化對區域生態系統服務的影響有助于優化配置土地資源和宏觀調控、全面認識生態環境變化趨勢并促進區域可持續發展[7-8]。

國內外相關研究已從最初基于土地利用數量描述生態系統服務價值[9-10]發展到應用多種指標、多種模型定量評價土地利用強度變化對生態系統服務的影響[7,11],并且隨著精細尺度的研究需求不斷增強,研究者的關注范圍也從全國[12]、省級[5]或區域[6]尺度逐漸轉向較小的空間尺度[4,13-14]。Chen等[4]研究了1995—2015年全國縣域生態系統服務強度與土地利用變化之間的聯系；郭椿陽等[13]基于格網探究了1990—2015年伏牛山區土地利用變化對生態服務價值的影響；雷軍成等[14]基于Landsat TM/OLI影像及當量因子法評估了1995—2016年江西省尋烏縣生態系統服務價值的變化。縣域是我國最基本的行政單元[15],近年來針對縣域土地利用變化及生態系統服務價值方面的研究不斷增多[14- 17],但仍存在不足:(1)往往僅從數量上描述時間序列下土地利用變化造成的生態系統服務價值變化,缺乏對生態系統服務價值損益進行空間分析；(2)僅對一級地類進行分析,缺乏對同一地類內部生態系統服務質量差異進行描述；(3)多集中于生態系統服務價值現狀評估,不能預測未來的價值損益動態,導致人們對未來縣域精細尺度下生態系統服務價值損益的變化規律仍不清楚。

目前預測未來土地利用變化的方法主要有三種,分別是經驗統計模型、元胞自動機和系統動力學模型[16,18-19],其中CLUE-S模型(Conversion of land use and its effects at small regional extent)是發展較為成熟且應用最廣泛的經驗模型,它基于土地利用與驅動要素的關系、對各種地類競爭及空間分配進行系統模擬,適用于中小尺度研究,而且具有較高的解釋力與可信度[16,20- 21]。該模型對縣域土地利用變化的預測是量化未來生態系統服務價值損益時空動態的基礎。

作為國家中心城市,重慶市在城鎮化推進過程中面臨著日益嚴重的環境保護壓力,其中巴南區是重慶市主城九區中面積最大、城鎮與農村發展相互作用最劇烈的區域[22],選擇該區作為研究對象對了解縣域尺度土地利用變化下生態系統服務價值損益的時空動態具有典型意義。基于2000年以來連續5期遙感影像及CLUE-S模型所預測出的未來10年巴南區8種土地利用類型的變化,評估2000—2030年巴南區生態系統服務價值損益的時空動態、構建影響了強度模型度量土地利用變化對生態系統服務價值的影響,旨在評估精細尺度下生態系統服務損益時空動態、提供科學合理的縣域土地資源優化配置方案、推動城鄉區域協調發展,統籌中、小尺度國土空間規劃并促進社會經濟及生態的可持續發展。

1 研究區域與研究方法

1.1 研究區域

巴南位于重慶市西南部(106°26′2″—106°59′53″E、29°7′44″—29°45′43″N),東接重慶市涪陵區和南川區,南連綦江區,西鄰江津區、九龍坡區、大渡口區,北傍南岸區、江北區、渝北區和長壽區(圖1)。巴南區東西橫距46 km,南北縱距70 km,總面積1.83×105hm2。巴南區地貌復雜多樣,以丘陵和低山為主,平地極少；地形起伏較大,海拔范圍為154—1133 m；年平均氣溫為18.7 ℃,年均降水量為1000—1200 mm,日照1100—1300 h；無霜期大于300 d,屬于亞熱帶濕潤氣候。全區生態環境良好,森林覆蓋率達47.6%,擁有樹種227種,動物資源近百種,并具有豐富的水資源與礦產資源。

圖1 巴南區位置及鄉鎮區劃Fig.1 Location and township division of Banan District

1.2 數據來源

本文以2000—2020年連續5期土地利用/覆被數據為主要數據源。該數據來源于Landsat TM/ETM+/OLI遙感影像(http://glovis.usgs.gov/),空間分辨率為30 m,利用Erdas 9.2軟件進行波段組合、幾何校正、輻射校正、裁剪等預處理,采用人機交互式目視解譯。根據IGBP、FAO及馬里蘭大學的全球土地利用分類體系,劃分為水域、建設用地、未利用地、耕地和林地5種土地利用類型；同時結合巴南區土地利用/覆被特點和研究需求,又將林地劃分出針葉林、針闊混交林、灌木林和常綠闊葉林這4種二級地類。基于森林資源調查、固定樣地等實測數據,應用隨機樣點進行精度驗證,所有影像的總體分類精度均高于85%,能夠滿足研究需求。海拔及地形數據來自地理空間數據云的SRTMDEM數據(http://www.gscloud.cn),空間分辨率為30 m。該區城鎮居民點、道路數據來源于國家基礎地理信息中心(http://www.ngcc.cn/),將這些作為預測未來土地利用格局變化的驅動要素。巴南區主要農作物類型、總產量、播種面積及價格等數據來源于《重慶市統計年鑒》[23],并結合2000—2015年《全國農產品成本收益資料匯編》[24]估算生態系統服務價值當量。

1.3 CLUE-S模型應用

1.3.1模型介紹

CLUE-S由荷蘭瓦赫寧根大學研究團隊開發,基于每種土地利用類型與多個驅動因素的經驗關系,對研究區內所有地類進行明確的空間分配,適宜于精細尺度土地利用變化的研究[20- 21],主要步驟如下。

①基于2000—2020年巴南區土地利用歷史變化趨勢,應用馬爾科夫鏈模型估算2025—2030年各類土地利用類型的面積。馬爾科夫鏈模型是基于時間序列的隨機過程,廣泛應用于地學數值模擬[18]。根據巴南區現有土地利用類型集合X(t)預測未來各地類集合狀態X(t+n),基本形式為:

(1)

式中n為時間步長；Pij為土地利用轉移概率矩陣。

②設置土地利用轉化規則及轉化彈性。根據巴南區現有土地利用變化規律設置土地利用轉化規則,將歷史時段時內不同土地類型間發生過轉化情況賦值為1,將不同地類間未發生過轉化情況賦值為0[16]。彈性系數描述了不同地類轉化的穩定性,參考重慶市主城區已有的研究成果[25]并經過多次模擬調試,確定巴南區水域、建設用地、未利用地、耕地、林地(含所有二級地類)的彈性系數分別為0.9、0.7、0.5、0.5、0.6。

③選擇地形因子(海拔、坡度、坡向)和區位因子(距道路距離、距城鎮距離)作為驅動力對各地類進行Logistic回歸(表1),得到各地類在空間某點的分布概率(Pi,u)。

表1 巴南區土地利用類型與驅動因子的Logistic回歸方程及ROC值

④根據分布概率、彈性系數以及迭代系數得到的總概率對各地類進行空間分配,公式如下。

TPROPi,u=Pi,u+ELSAu+ITERu

(2)

式中,TPROPi,u為柵格i是土地利用類型u的總概率；Pi,u為Logistic回歸求得的空間分布概率；ELASu為土地利用類型u的彈性系數；ITERu為迭代系數。

巴南區8種土地利用類型與驅動因素的逐步回歸結果表明所有土地利用類型均與坡度的關系最密切,同時所有土地利用類型回歸模型的ROC(Relative operating characteristics)曲線面積接近或大于0.7,表明逐步回歸模型能夠較好地解釋巴南區不同地類的發生概率。

1.3.2馬爾科夫鏈模擬檢驗及CLUE-S模擬結果檢驗

應用卡方檢驗對馬爾科夫鏈所預測的未來各種土地利用面積結果進行驗證。使用歷史土地利用轉移概率矩陣估算2020年巴南區各地類面積,并與其真實土地利用格局進行比較(圖2),應用Kappa系數檢驗各種地類空間分布格局預測的準確性。卡方檢驗結果和Kappa系數分別由公式(3)和公式(4)計算得到。

圖2 CLUE-S預測結果及巴南區實際土地利用格局對比Fig.2 Comparison of CLUE-S simulated result with the real land-use pattern of Banan District

(3)

Kappa=(Pa-Pc)/(Pp-Pc)

(4)

1.4 生態系統服務價值評估

參考謝高地等[3,26]所提出中國陸地各種生態系統類型單位面積生態系統服務價值當量及生態經濟學方法[27](公式5—6),構建了巴南區不同土地利用類型的各項生態系統服務價值當量表(表2)。建設用地可提供的生態系統服務價值很低,不參與計算[28]。生態系統服務價值當量因子的經濟價值為當年全國平均糧食市場價值的1/7[26]。由《重慶市統計年鑒》[23]和《全國農產品成本收益資料匯編》[24]得到2000—2015年全國糧食平均出售價格范圍為1098.80—2325.60 元/t,此期間全國與巴南區單位面積糧食產量范圍分別為4.16—7.01 t/hm2、3.68—3.89 t/hm2,獲得巴南區生物量修正因子為0.56—0.88。由于目前尚未能獲得2020年糧食產量及價格數據,因此在2020—2030年也以2015年的當量因子價值進行估算,確保未來生態系統服務價值評估與歷史時段結果具有可比性。同時為消除物價變動對生態系統服務價值評估的影響,本文使用居民消費價格指數(Consumer price index,CPI)以2015年為基期,對2000—2010年當量因子價值進行修正[11],得到2000—2030年巴南區生態系統服務價值的當量因子價值范圍為818.53—1335.56 元(圖3)。

圖3 巴南區生態系統服務價值當量因子價值 Fig.3 Per unit value of equivalent of ecosystem service values (ESV) of Banan District

表2 巴南區不同土地利用類型生態系統服務價值當量

1.5 土地利用變化對生態系統服務價值的影響

使用單一動態度和綜合動態度分別描述各種土地利用類型及研究區整體土地利用變化的速率趨勢及劇烈程度[29]。基于土地利用轉移矩陣核算得到的生態系統服務價值損益矩陣可揭示服務價值的數量變化特征、增減方向及來源[7]；應用劉亞茹等[7]所構建的影響程度指標來評價土地利用變化對服務價值的影響程度(公式7—10)。為從精細尺度分析2000—2030年巴南區生態系統服務價值損益的空間分布格局,本文使用格網作為研究單元對巴南區生態系統服務價值損益做探索性空間分析,了解損益的空間自相關特性。常見的格網尺度有500 m×500 m、1 km×1 km、3 km×3 km、5 km×5 km及10 km×10 km等[13,16]。考慮到本文研究尺度精細以及最小可塑性單元對結果的影響,故而增加120 m×120 m和250 m×250 m這兩種尺度,通過分析不同格網尺度下生態系統服務價值損益時空動態及差異,發現以250 m×250 m網格尺度的分析結果能夠充分體現巴南區生態系統服務價值空間自相關特征,是現有研究尺度、數據量和空間斑塊情況下最合適的方案。因此本文將巴南區按250 m×250 m大小進行魚網劃分。按照自然斷點法將生態系統服務價值損益從高到低分為5級,使用Geoda軟件進行價值損益的全局自相關分析分別計算全局、局部Moran′sI指數及局部Gi指數。

(5)

ESV=∑(Ai×Vi,k)

(6)

LPi=∑Ai×(VCi,k-VCi,k′)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

2 結果分析

2.1 巴南區土地利用變化

2000—2020年巴南區耕地面積由1.39×105hm2減少了38.77%(圖4)。減少的耕地76.32%轉化為林地、17.97%轉化為建設用地；轉化為水域和未利用地的比例最低,分別為4.82%和0.88%。由于耕地的轉化,林地面積增長最多(4.1×104hm2),未利用地增加最少(472.0 hm2)。CLUE-S模型預測上述土地利用變化趨勢將在未來10年內延續。2020—2030年耕地繼續減少1.85×105hm2,其中63.01%轉化為林地,29.62%轉化為建設用地。2000—2030年耕地占巴南區面積的比例由75.99%降低到36.41%,而林地占比增長了28.86%并將在2025年后取代耕地成為最主要的土地利用類型,其中針葉林、灌木林、針闊混交林和常綠闊葉林的面積占巴南區面積比例增長依次降低,分別增長了13.42%、7.15%、7.07%和1.22%(圖4)；建設用地和水域的比重分別增長8.29%和2.25%；未利用地占比最小,其所占比重僅增長0.18%。

圖4 2000—2030年巴南區不同土地利用類型面積變化Fig.4 Area changes of different land use type of Banan District during 2000 to 2030

2000—2010年是巴南區土地利用變化最為劇烈的時期,土地利用綜合動態度由0.67%增加到1.08%(圖5),主要來源于建設用地和林地的迅速增長,其中建設用地單一動態度最高,由12.16%增加到29.29%(圖5),林地單一動態度也較高(4.94%—6.40%),說明建設用地和林地是轉化最顯著的兩種土地利用類型。此后巴南區土地利用變化速率放緩。2020年以后所有土地利用的單一動態度均趨于平緩,綜合動態度在2025—2030年最低(0.48%)。2000—2030巴南區建設用地、水域、林地、未利用地和耕地的單一動態度依次降低,分別為41.58%、5.26%、4.55%、0.83%和-1.74%；耕地的單一動態度始終小于0,即耕地持續轉化為其他土地利用類型。

2.2 巴南區生態系統服務價值變化

2000—2020年巴南區生態系統服務價值由11.69億元增長到27.64億元(表3),此后仍逐漸增長,到2030年達到30.19億元。耕地生態系統服務價值在2010年達到最大值5.71億元,此后逐漸降低至3.38億元。2000—2030年林地生態系統服務價值共增長16.19億元,其中針葉林和針闊混交林是巴南區生態系統服務價值增長最多的兩種地類,分別增加7.61億元和5.06億元；水體生態系統服務價值增長3.40億元,未利用地生態系統服務價值變化極小。2000—2030年耕地及林地的面積變化對生態系統服務價值影響最明顯(圖5),平均影響程度分別為32.24%和30.18%；水域面積增長對服務價值的影響逐漸增強,而未利用地則基本無影響,二者影響程度分別為5.60%和1.29%；建設用地對服務價值的影響程度為-2.39%,表明其面積增加將導致生態系統服務價值的降低。2000—2030年除食物生產的服務價值將減少外(減少量為0.08億元),其他服務價值均增加(表4)。氣候調節價值增量最高,為4.82億元,而維持養分循環價值增量最低,為0.11億元。耕地面積減少引起了食物生產價值持續下降,而林地、水域面積的增長則是其他生態系統服務價值增長的主要原因。

表3 巴南區不同土地利用類型生態系統服務價值動態/億元

表4 巴南區不同類型生態系統服務價值動態/億元

2.3 生態系統服務價值損益時空動態

2000—2030年巴南區生態系統服務價值損失0.58億元、總增益19.09億元,凈增益為18.50億元。耕地、林地轉化為建設用地是生態系統服務價值損失的主要來源,分別占全區損失的80.64%和8.42%(表5)；而耕地轉化為林地是服務價值增益的最主要方式,貢獻全區總增益的57.84%。土地利用結構變化和當量因子價值升高共同導致了生態系統服務價值增長,二者分別貢獻了凈增益的70.96%和29.04%。

表5 2000—2030年巴南區生態系統服務價值損益矩陣/萬元

2000年巴南區耕地占絕對優勢,使全區生態系統服務價值呈現較為均勻的空間分布格局,僅北部水系及海拔較高的桃子蕩山的生態系統服務價值較高。隨著時間推移,城鎮化發展的壓力則使地勢平坦的中西部地區耕地轉化為建設用地,使生態系統服務低值區由最初零星分布于西部邊緣逐漸向中部及北部擴張。同時該區中南部、東南部海拔較高的山區耕地大面積轉化為以針葉林為主的林地,使生態系統服務價值明顯增加,使全區生態系統服務價值分布格局差異逐漸顯現(圖6)。探索性空間分析結果顯示(表6),巴南區生態系統服務價值損益的全局Moran′sI指數及P值分別為0.6077、0.001,說明該區生態系統服務價值損益表現為顯著的正向空間自相關。若將巴南區由東北到西南劃一條對角線,其西側海拔較低、地勢平坦而東側海拔較高。生態系統服務價值冷點集中在對角線西側海拔較高的李家沱、花溪街道、南泉街道等地區,表現為“低值-低值”集聚特征和冷點(圖7)；而生態系統服務價值增益熱點主要分布于對角線東側海拔較高的石龍鎮、石灘鎮、接龍鎮、跳石鎮南部及安瀾鎮西側,呈現“高值-高值”集聚特征。

表6 巴南區生態系統服務價值損益全局自相關

圖6 2000—2030年巴南區地均生態系統服務價值空間分布格局及損益Fig.6 Spatial distribution pattern and the profit and loss of ecosystem service values in Banan District

圖7 巴南區生態系統服務價值損益的局域空間自相關分析Fig.7 Local spatial autocorrelation of the profit and loss of ecosystem service values in Banan District

3 討論

本文基于土地利用變化評估并預測了縣域生態系統服務價值損益時空動態,為未來區域土地資源的合理規劃及有效利用提供了科學依據,對保護生態環境和維持區域社會經濟發展具有現實的指導意義。過去20年間巴南區最明顯的變化是耕地持續轉化為其他土地利用類型,基于該趨勢馬爾科夫鏈模型預測出2030年巴南區耕地面積和林地面積分別為6.64×104hm2、9.13×104hm2,符合《巴南區及龍洲灣街道等22個鎮街土地利用總體規劃(2006—2020年)調整方案》對永久基本農田保護紅線(4.66×104hm2)和生態保護(1.97×104hm2)的面積規定,表明按現有趨勢發展未來巴南區土地利用變化預測合理。

基于CLUE-S模型對精細尺度土地利用格局的準確模擬,不僅能夠獲得未來生態系統服務價值的數量變化,還可直觀顯示出不同土地利用類型及同一地類內部二級地類的生態系統服務價值損益的空間分異。2000—2030年巴南區土地利用結構變化導致的價值增長為13.13億元,占凈增益的70.96%；而當量因子價值升高則貢獻了凈增益的29.04%。巴南區生態系統服務價值損益與地形密切相關。生態系統服務價值在西部低海拔地區隨著城鎮化加劇以及耕地轉化為建設用地而損失逐漸加劇,而在中南部及東南高海拔地區則隨著耕地轉化為林地而明顯升高,使價值損益在全局和局域尺度上均呈現顯著的正向空間自相關特性。與我國同期大部分城市發展導致生態系統服務價值持續減少的情況不同[7,11],盡管巴南區建設用地逐漸擴張,但林地面積的大幅增長使該區生態系統服務價值增益遠高于其損失。這不僅反映出過去十幾年間天然林資源保護工程及退耕還林工程的實施能夠有效緩解城鎮化與生態環境保護的沖突,而且也說明現有土地利用發展模式能夠兼顧經濟發展及生態環境保護兩方面的需求。

縣域是介于省市及鄉鎮之間的中間尺度,優化縣域空間布局可有效推動大中城市和小城鎮的協調發展。近期在重慶市“一區兩群”發展格局規劃中巴南區是主城都市區之一,同時也是“國家城鄉融合發展試驗區”之一,依據本文所得結果可從以下3方面優化該區未來空間布局,以構建科學合理的城鎮體系:(1)地勢較平坦西部及北部地區應提高已有建設用地的利用強度、合理規劃城市用地邊界以避免減損水域及林地等生態系統服務價值高值區、促進產業結構調整和升級以增強中心城市綜合承載能力；(2)南部海拔較高的地區繼續進行退耕還林等生態保護工程、推行造林綠化、加強相鄰區縣間生態環境保護的聯防聯治、促進區域生境質量的全面提升；(3)由于耕地面積持續減少會引起食物生產服務的減弱,因此需要嚴守耕地紅線、使用先進技術提高糧食產量和保障糧食安全、搭建城鄉產業協同發展平臺,推動城鄉區域協調發展。

本文局限性主要包含以下3個方面:(1)CLUE-S難以反映政策對土地利用的直接影響,所使用的Logistic逐步回歸方程所擬合的土地利用類型與驅動因素的關系是單向作用過程,不能反映土地利用類型變化的反饋機制與過程[30]；(2)由于可獲得數據及方法的限制,本文未能預測未來生態系統服務當量因子經濟價值的變化,仍延用2015年的經濟價值進行計算會造成未來生態系統服務價值低估或高估；(3)僅按照當前發展趨勢對未來土地利用格局進行預測,未能設置多種土地利用情景進行對比評估。今后需要應用更為準確的生態系統服務價值核算體系進行精細評估,同時需依托各類生態系統定位觀測研究站獲取連續的自然驅動要素與社會經濟數據,設置多種發展模式進行模擬,加強對微觀土地利用變化過程與機理的研究、增強對生態系統服務流以及多種生態系統服務權衡與協同關系的認識。

4 結論

基于土地利用變化及CLUE-S模型可以直觀評價并預測未來縣域尺度生態系統服務價值損益時空動態。2000—2030年巴南區生態系統服務價值損失0.58億元、總增益19.09億元,凈增益18.50億元；土地利用結構變化和當量因子價值升高分別貢獻了凈增益的70.96%和29.04%。耕地、林地轉化為建設用地是導致巴南區生態系統服務價值損失的主要途徑,而耕地轉化為林地則促進生態系統服務價值增益。巴南區生態系統服務價值損益與地形密切相關,呈現顯著的正向空間自相關特性。巴南區現有土地利用發展模式能夠兼顧經濟發展及生態環境保護兩方面的需求,但未來仍需實施生態工程、合理規劃城市用地邊界、嚴守耕地紅線以實現區域長期可持續發展。