不同土地利用情景下汾河上游地區碳儲量評估

張 燕,師學義,唐 倩

中國地質大學(北京),北京 100083

大氣中二氧化碳含量增加是國際社會廣泛關注的環境問題,如何減少碳排放、提升碳儲存能力成為了研究熱點和難點。陸地生態系統通過釋放和吸收CO2、N2O等溫室氣體調節區域氣候,是全球碳循環過程的重要環節[1]。陸地生態系統碳儲量提高能夠有效減少大氣中CO2的含量,是緩解溫室效應、減緩全球氣候變化的重要方式之一。人類活動改變土地利用格局,而土地利用變化進一步影響陸地生態系統結構、功能,改變陸地生態系統固碳能力,是陸地生態系統碳儲量變化的主導因素[2]。因此,預測不同情景下土地利用變化,定量評估土地利用變化所致的陸地生態系統碳儲量,為未來土地利用結構優化方向和陸地生態系統碳儲量提升途徑提供了指導作用,為區域可持續發展和氣候變化減緩做出貢獻[3]。

近年來,眾多學者基于土地利用格局變化對單一或綜合性陸地生態系統碳儲量進行了研究。研究結果表明,全球范圍內熱帶地區森林植被減少、森林生態系統功能退化造成了全球陸地生態系統碳儲量的下降[4-5]；陳耀亮等運用簿記(Bookkeeping)模型評估中亞地區森林生態系統碳儲量,指出森林砍伐造成了碳儲量下降且土地利用/覆蓋變化對碳儲量有顯著影響[6]；劉亞男等研究了中國主要濕地碳儲量并指出土地利用變化造成了濕地生態系統功能的退化,嚴重影響濕地的“碳匯”功能[7]；柯新利等人基于InVEST模型研究湖北省陸地生態系統碳儲量,認為其下降的主要原因是城市建設用地擴張和耕地保護[8]。這些研究大多通過InVEST模型、Bookkeeping模型評估碳儲量,表明了土地利用變化對陸地生態系統碳儲量有顯著影響,為土地利用結構優化提供了新的科學依據。一些學者運用模型模擬了未來土地利用格局與碳儲量變化,卻沒有考慮不同土地利用情景下的差異[9-11]。流域作為一個復雜系統,包括各種自然資源、生態系統等要素,具有明顯的層次結構和整體特征,從流域尺度分析并解決生態環境問題是切實可行的有效途徑[12]。然而,針對流域尺度,結合時間和空間尺度預測不同土地利用情景與生態系統碳儲量的研究并不多見,預測未來土地利用格局,模擬不同土地利用情景,評估生態系統碳儲量,有助于區域生態環境發展,促進土地利用優化配置。

本文選擇汾河上游作為研究區,該區域位于山西省中北部,汾河發源地,對山西省整體生態環境有重要影響,近年來人類活動對土地利用格局產生巨大影響,生態系統碳儲存功能下降。因此,本文以2007年、2011年、2017年汾河上游土地利用數據為基礎構建了SD_CLUE-S復合模型,預測2030年生態保護情景和自然增長情景下各土地利用類型面積,利用InVEST模型對2007—2017年及不同情景下2030年生態系統碳儲量進行估算并分析碳儲量空間分布情況,為汾河上游土地利用規劃提供科學依據,促進區域生態系統碳儲存功能提升,對山西省整體生態環境可持續發展具有重要意義。

1 研究區概況

汾河流域(38°7′57″—39°8′41″N,111°42′43″—112°36′55″E)位于山西省忻州市中部,東依原平市、忻府區,西隔管涔山與五寨、岢嵐縣為鄰、南連太原市的婁煩縣、古交市,北靠神池縣、朔州市的朔城區(圖1),總面積達到了399737.47hm2。空間上總體呈現“兩山一川”大格局,以蘆芽山、云中山為兩翼,以汾河川為中心從南到北貫穿整個區域,自然生態要素齊備,生態功能完整。該地區位于北溫帶大陸性季風氣候,多年平均氣溫6.30℃,最高和最低均出現在河川階地區；降水多集中在7—9月份,歷年平均降水量420—704mm,區域內平均年蒸發量為1812.60mm；境內山地多,整體海拔較高,地勢呈中部低、東西兩邊高,境內草地、林地、耕地所占面積較大,建設用地面積逐漸擴大,水體面積多年來持續下降。研究區是山西省重要的水源涵養區、華北平原重要生態屏障,全省汾河水源地保護與治理重點區域,生態功能系統碳儲存功能脆弱是制約經濟發展和社會進步的突出因素。

圖1 研究區地理位置概況

2 材料與方法

2.1 數據來源及處理

研究區2007、2011、2017年遙感影像來自地理數據空間云(http://www.gscloud.cn),運用ENVI軟件進行遙感解譯且經過野外驗證及Google Earth Pro驗證精度達到了89.20%、91.45%,90.87%,根據《土地利用現狀分類(GB/T21010—2007)》將研究區3期土地利用數據分為草地、耕地、建設用地、林地、未利用地、水體6種類型。

研究區年降雨量及年均氣溫來自《忻州市統計年鑒》；DEM數據來源于ASTER GDEM,驅動因子中海拔因子和坡度因子數據來自DEM數據的ArcGIS處理；土壤質地因子來源于聯合國糧農組織(FAO)和維也納國際應用系統研究所(IIASA)所構建的世界土壤數據庫以及第二次全國土地調查南京土壤所所提供的1∶100萬土壤數據；距離因子中省道、縣道、鄉道、河流、鐵路、農村居民點以及縣域中心的數據均來自全國地理信息資源目錄服務系統,運用歐式距離法得到研究區每個柵格到某空間驅動因子的距離；人口因子數據來源于地理國情監測云平臺中全國1km×1km的人口柵格數據。

2.2 研究方法

2.2.1系統動力學

系統動力學基于信息反饋系統認識并解決系統問題,認為系統的行為模式與特性主要取決于系統內部的動態結構與反饋機制[13],一般用系統結構框圖、因果關系圖以及流圖來構建特定模型進行長期預測,以滿足實際應用需要,在研究數據不夠充分的情況下尤其適用。

本文基于2007—2017年研究區域土地利用轉移矩陣信息構建了土地利用變化的系統動力學模型(圖2),各土地利用類型設置為水平變量,土地轉移面積設為流速變量[14],轉移概率設為影響轉換面積的常量,并根據實際情況對轉移概率進行修正和調整,以便更精確地預測2030年土地利用格局。

圖2 研究區土地利用變化SD模型

2.2.2CLUE-S模型

CLUE-S模型分為非空間需求和空間分配模塊,其中,非空間需求模塊確定土地利用變化驅動因子及其與土地利用變化分布概率之間的關系[15],本文通過二元Logistic回歸模型及ROC曲線驗證了驅動因子選取的精度,同時根據土地利用的SD模型預測2030年研究區各土地利用類型面積。空間分配模塊分為輸入和輸出兩個階段,輸入階段包括土地利用轉換規則,限制性區域、空間特征以及非空間土地需求模塊求得的土地需求結果[16],通過迭代將2030年土地利用類型的空間分配結果輸出。

根據《忻州市土地利用總體規劃(2006—2020年)調整方案》將各生態修復區域及基本農田區域選定為限制區域(圖3),設定生態保護情景下草地、耕地、建設用地、林地、未利用地、水體的轉移彈性值為0.5、0.6、0.7、0.3、0.7、0.6,并設置土地利用變化轉移矩陣(表1)；自然增長情景下各土地利用類型ELSA值設為0.5、0.5、0.7、0.3、0.9、0.9,同時將土地利用變化轉移矩陣中水體轉移為未利用地設為0,其余皆為1,并且不設置限制區域。

圖3 2030年生態保護控制區

表1 生態保護情境下土地利用變化轉移矩陣

2.2.3InVEST模型

InVEST模型根據不同土地利用類型對應的地上生物碳庫、地下生物碳庫、土壤碳庫和凋落物碳庫估算碳儲量[17],并生成平均碳密度空間分布圖。根據研究區具體實際情況以及數據獲取難度,本文選取對碳儲量評估影響較大的地上生物碳庫、地下生物碳庫以及土壤碳庫進行估算,計算方法如下：

Ckt=Ckabove+Ckbelow+Cksoil

(1)

(2)

式中,Ckt表示土地利用類型k的碳密度總量(t/hm2)；Ckabove表示土地利用類型k的地上生物碳密度(t/hm2)；Ckbelow表示土地利用類型k的地下生物碳密度(t/hm2)；Cksoil表示土地利用類型k的土壤碳密度(t/hm2)；C表示研究區碳儲量(t);Sk表示土地利用類型k的面積(hm2)。

根據李克讓、解憲麗等人的研究得到全國六種地類的碳密度數據[18-19],基于S.A.Alam、陳光水等人提出的公式[20-21]對碳密度進行修正(表2),由于死亡有機物碳密度難以獲得且所占比例不大,因此將其填為0。

表2 研究區四大碳庫碳密度

CSP=3.3968×MAP×+3996.1

(3)

CBP=6.798×e0.0054×MAP

(4)

CBT=28×MAT+398

(5)

式中,CSP表示年降雨量修正下得到的土壤碳密度(kg/m2)；MAP表示年均降雨量(mm)；CBP和CBT分別表示年降雨量和年均氣溫修正下得到的生物量碳密度(kg/m2)；MAT指的是年均氣溫(℃)。

(6)

(7)

KB=KBP×KBT

(8)

(9)

式中,KBP表示生物量碳密度年降雨量因子修正系數；KBT表示生物量碳密度年均氣溫因子修正系數；KB表示生物量碳密度修正系數；KS表示土壤碳密度修正系數；C′和C″分別表示研究區及全國相應數據。

3 研究結果

3.1 2007—2017年汾河上游地區土地利用變化分析

3.1.1土地利用變化幅度分析

2007年與2011年草地是研究區主要土地利用類型,2011年與2007年相比草地面積下降且降幅為3%,耕地、建設用地、林地以及未利用地面積增加,其中林地增加幅度最大達到了2.29%,其次是耕地增加幅度2.03%,最小的是未利用地。此外,2011年水體面積減少了7209.89 hm2(如圖4、表3所示)。

2011—2017年草地仍然是土地利用的主要類型,但是面積仍舊在下降且降幅達到了2.49%,此外未利用地與水體面積存在不同程度下降,其中未利用地面積減少了4080.85 hm2,主要原因是人們加大了對未利用地的開發程度；水體面積下降了0.24%,近年來水體面積的持續下降對研究區域的生態環境產生了惡劣影響。耕地、建設用地、林地繼續保持上升趨勢,其中耕地面積增加了9468.66 hm2,建設用地增加幅度達到了0.49%,林地面積增加了3547.33 hm2(如圖4、表3所示)。

圖4 研究區土地利用數據

表3 研究區土地利用面積及占比

2007—2017年,耕地、建設用地、林地面積持續上漲,增加幅度分別為4.40%、0.74%、2.18%,草地、未利用地、水體面積下降,其中草地面積降幅最大為5.49%,其次是水體為2.04%,減少面積達到了8144.77 hm2,未利用地面積則下降了3137.40 hm2(如圖5所示)。

圖5 土地利用變化幅度分析

3.1.2土地利用類型轉移矩陣分析

草地轉出主要流向耕地,其次是林地,建設用地、水體與未利用地所占比重較小；耕地主要轉移成了草地,林地、建設用地、水體、未利用地次之；建設用地主要流向耕地,轉出比例達到了41.11%,其次是草地；林地轉化為其他土地利用類型的程度最小,其中轉為草地的比重最大；未利用地轉移成林地比例最大,其次是草地、耕地與建設用地；對水體的占用主要是草地與耕地,兩者的轉移概率總和將近80%,其次是建設用地、林地,未利用地沒有占用水體(如表4、表5所示)。

表4 2007—2017年土地利用轉移矩陣

表5 2007—2017年土地利用轉移概率/%

3.2 土地利用模擬預測結果及分析

3.2.1基于SD模型的土地利用仿真及預測

根據系統動力學模型得到2007—2017年每一年各個土地利用類型面積的仿真結果(表6),將其與2017年解譯結果進行對比,確定SD模型構建的可行性和準確性。

表6 2007—2017年各土地利用類型SD模型預測面積

根據仿真值與解譯值對比可以看出,所有土地利用類型仿真結果的相對誤差都在1%以下(表7),說明SD模型仿真結果與實際情況基本一致,可以進行下一步的預測,即輸入2017年各土地利用類型面積作為初始值,設置步長為13年,預測2030年研究區土地利用情況。

表7 2017年仿真與解譯值對比

根據系統動力學模型預測結果可以看到(表8),研究區域在未來13年間各土地利用類型變化面積較大,其中耕地變化面積最大,減少了1117.00 hm2；林地增加量次之,達到了1045.46 hm2；建設用地和未利用地面積分別減少了181.10 hm2和101.69 hm2；草地和水體分別增加了234.00 hm2和118.16 hm2。

表8 2017—2030年各土地利用類型SD模型預測面積

3.2.2土地利用變化驅動因子分析

本文選取了11個驅動因子進行分析,其中距離因子包括距省道的距離、距縣道的距離、距鄉道的距離、距鐵路的距離、距農村居民點的距離、距河流水系的距離和距縣中心的距離；地形因子包括高程和坡度；土壤因子與社會經濟因子分別為土壤類型和人口密度因子。

本文通過空間數據分析、二元Logistic回歸分析和ROC曲線測算,得到驅動因子的B值和Exp(B)(如表9所示),其中B值是各驅動因子的回歸系數,Exp(B)表示發生比率,即驅動因子每變化一個單位,各土地利用類型發生比的變化；ROC曲線能夠有效驗證土地利用模型擬合度,一般來看ROC值大于0.5即可表明驅動因子對土地利用類型空間分布具有較好的解釋能力,值越大表明擬合度越高,根據計算得到ROC值均大于0.65,擬合度較高,符合要求。

表9 驅動因子分析

研究結果表明,研究區草地分布與縣道、鄉道、農村居民點、河流水系距離以及坡度呈正相關,與海拔因子呈負相關,其中坡度對其影響最大,坡度每增加1°,草地分布概率就增加2.9%；耕地與距農村居民點的距離、與縣中心的距離、高程、坡度呈負相關,與距河流水系的距離及人口密度因子無關,與其他因子為正相關,其中縣域中心距離對耕地分布影響最大,距離每增加1m,耕地分布概率就減少62.5%,；建設用地分布與縣域中心距離相關性最強；土壤類型對林地分布影響最大,將土壤類型設置為壤土-0,砂土-1,壤砂-2,砂壤-3,粉粘-4,因此土壤類型對應編碼每增加一個單位,林地分布概率就減少14.96%；未利用地分布僅受到海拔與坡度的正相關影響,其中海拔對其分布概率影響最大,海拔每增加1m,未利用地分布概率就增加36.89%；水體分布概率與縣域中心距離相關性最強。

3.2.3基于生態保護情景的土地利用變化模擬

將2007年數據作為起始年份數據,仿真模擬2017年土地利用數據驗證Kappa精度,保證2030年土地利用數據預測準確性,如圖6所示,模擬正確的柵格有59973個,總柵格數有63704個,則P0為0.94,PC為0.30,PP為1,得到Kappa系數為91.58%,即2017年研究區域模擬精度較高,完全滿足進一步的研究需要,輸入自然情景以及生態保護下預測2030年土地利用數據的相關基礎數據,運行模型并輸出結果(圖7)。

圖7 兩種情景下2030年研究區土地利用數據

自然增長情景下,林地面積大幅下降,減少了近52908.26hm2,相較于生態保護情景減少了近54925.56hm2；自然增長情景下2030年草地面積達到了198100.00hm2,而生態保護情景下僅有145325.00hm2,主要原因是ELSA值、轉移矩陣及空間分布概率影響；生態保護情景下退耕還林、還草等生態保護措施實施耕地面積下降,與之相比自然增長情景下耕地面積增加了467.73hm2；建設用地在生態保護情景下發生增長而在自然增長情景下減少,兩種情景下面積相差較小；生態保護情景下大部分未利用地轉化為林地、草地及水體,且利用面積相較于自然增長情景多了1230.00hm2；研究區作為汾河的重要源頭,水體占比具有重要意義,2007—2017年水體面積不斷下降,水資源環境持續惡化,在自然增長情景下水體仍舊呈現下降趨勢,而生態保護情景下實現了142.16hm2的增長。

3.3 汾河上游地區碳儲量評估

生態保護情景下,2007—2030年研究區碳儲量呈穩步上升趨勢(圖8),碳密度持續增長,其中2007—2011年碳儲量增幅最大,達到了868423.72t,2011—2017年碳儲量增量下降為368719.31t,2017—2030年研究區碳儲量增加了341217.04t；2007年碳密度為144.45t/hm2,2030年碳密度達到了147.95t/hm2,且2007—2011年的碳密度變化最大；從土地利用類型來看,草地碳儲量在2007—2017年間持續下降,2030年與2017年相比發生了51714.05t的增長,2007—2017年耕地碳儲量持續增長,2017—2030年下降且降幅達到了160758.96t；建設用地碳儲量在2011—2017年間增加量最大,其中增幅最小為8550.58t；未利用地在2007—2011年間有小幅增長,在2011—2030年持續下降,且2011—2017年下降幅度較大；水體在2007—2017年間碳儲量持續下降,在2017—2030年發生了增長,其中2007年下降幅度最大,達到了11455.27t(圖9)。

圖8 2007—2030年研究區碳儲量及碳密度變化

圖9 2007—2030年各土地利用類型碳儲量

自然增長情景下,2017—2030年碳儲量下降176917.86t,碳密度變化0.41t/hm2,相較于生態保護情景碳儲量減少341217.04t,碳密度下降0.85t/hm2,從土地利用類型分析,草地碳儲量發生極大增長,比生態保護情景多8176915.15t,同時耕地與未利用地分別增長228058.52t和47338.88t；反之,林地大幅度下降且減少了8771373.06 t,建設用地下降21579.67t,水體碳儲量變化較小,有576.87t下降量。

本文將碳儲量變化量在≤-150t/hm2范圍設定為下降,處于-150t/hm2到150t/hm2設定為基本不變,≥150t/hm2表示增加。由此可知,2007—2011年碳儲量下降區域主要以帶狀分布在靜樂縣汾河流域周圍以及寧武縣陽方口工礦鎮區域,其余較為零散地分布在研究區林地范圍內,增長主要發生在靜樂縣縣城、寧武縣蘆芽山自然保護區、縣內汾河及汾河支流附近。2011—2017年靜樂縣碳儲量增長區域集聚分布在以汾河流域為中心的范圍內,而下降則發生在汾河流域內部；寧武縣碳儲量沿蘆芽山自然保護區邊界、恢河流域及縣城周圍村莊下降,碳儲量增長主要在縣城中心范圍、東馬坊鄉及荷葉坪區域(圖10)。

圖10 2007—2017年碳儲量變化

生態保護情景下2030年碳儲量下降區域集中在研究區汾河流域,增長區域零散分布在寧武縣蘆芽山自然保護區、恢河流域以及東馬坊鄉；自然增長情景下碳儲量基本沒有變化,僅在寧武縣恢河流域附近發生增長(圖11)。

圖11 兩種情景下2017—2030年碳儲量變化

4 結論和討論

4.1 討論

(1)總體來看,2007—2017年汾河上游社會經濟發展,生產活動頻繁,人們對耕地、建設用地的需求增加,導致耕地、建設用地面積持續增長；寧武縣工礦企業增加造成了水體污染,同時人們生產生活中不合理利用水資源導致了汾河水量、水質不斷下降,水體面積逐年萎縮；汾河上游擁有大面積的原始森林,隨著政府實行林地管控及退耕還林政策,林地面積發生增長。

(2)系統動力學(SD)能夠反饋系統內部要素相互作用,借助模型對復雜問題實行運算,能夠有效預測不同情況下土地利用變化、土地利用需求等,系統仿真提高了SD模型預測精度,但是SD無法獲得空間信息和空間分布預測結果[22]。CLUE-S模型的空間模塊能夠模擬土地利用空間分布,但是無法獲得土地利用需求量,在應用中人們往往會結合其他模型進行多領域的研究,高志強等和將其與Dinamica EGO模型結合,分析了中國未來土地利用格局和土地利用變化驅動力因子[23]；周銳等結合Markov模型預測了辛莊鎮在3種不同情景下未來土地利用格局的變化情況[24]；劉靜怡等結合灰色線性規劃方法預測了嘉興北部地區不同情景下的土地利用格局,并進行土地利用結構的優化[25]。然而這些方法主要運用數學方法進行土地利用需求預測,忽略了人為決策因素,或者需要大量數據支撐,這在實際研究中可能難以獲得,SD模型可以很好地解決這些問題。因此本文運用SD_CLUE-S復合模型消除了單一模型的不足,避免了其他模型的劣勢,兩者相輔相成,提高了土地利用格局預測精度,為進一步研究提供技術支撐,但是不足之處在于SD模型構建中缺少了社會經濟因素、政策因素考量,這些影響因素的量化和空間表達是未來研究重點；此外空間尺度對土地利用格局影響較大,本文分析一種空間尺度下研究區土地利用變化,今后應當設置不同空間尺度對比分析。

(3)目前陸地生態系統碳儲量評估大多通過采樣或植被、土壤清查數據計算獲得,具體方法可分為地表測量、模型模擬和遙感估算三類[26],其中地表測量大多是通過采樣、實地勘測等,結果精確性高,但實施難度較大；遙感估算是當前研究的新方向和新思路,通過遙感影像人們不僅可以獲得地表信息,還能對檢測地下生物碳儲量,但遙感影像精度會影響研究結果；模型模擬方法中人們普遍采用InVEST模型,數據獲取方便,結果準確性較高,但是模型忽略了植被種類、植被生長情況等對碳密度的影響,因此在研究成果基礎上對研究區各土地利用類型碳密度進行修正,這比以往直接引用全國碳密度值準確性高,不足之處是缺乏實地測量數據對修正結果進行驗證,只能通過相近地區研究結果進行對比,在今后研究中可以采用實驗和預測模型或實地調研和預測模型相結合的方法,提高碳儲量估算精度。

4.2 結論

通過SD_CLUE-S模型模擬汾河上游生態保護情景和自然增長情景下未來土地利用格局,并運用InVEST模型估算生態系統碳儲量。結果表明,2007年汾河上游碳儲量和碳密度分別為57740767.95t和144.45t/hm2,2007—2017年汾河上游碳儲量持續增長,但增長幅度明顯減緩。2030年汾河上游自然增長情景下碳儲量和碳密度出現較為明顯的下降,與2017年相比分別減少176917.86t和0.44t/hm2。生態保護情景改變了這種趨勢,碳儲量和碳密度顯著增加,增幅達到了2164299.18t和0.41t/hm2,增強了汾河上游固碳能力。未來對汾河上游土地利用結構優化,應基于生態保護情景,對汾河進行河道治理,提高河流岸坡生態環境,保證汾河水質、水量,積極保護林地,合理規劃蘆芽山保護區的旅游資源,限制工礦用地過度擴張,確保基本農田數量,這樣在提高碳儲量的同時也可以保證汾河上游地區社會經濟可持續發展。