解天琪 李 龍 陳 鑫 周 翼 張國梁 陳龍乾
(1.中國農(nóng)業(yè)大學 土地科學與技術學院,北京 100193;2.中國礦業(yè)大學 環(huán)境與測繪學院,江蘇 徐州 221116;3.荷語布魯塞爾自由大學 地理系,比利時 布魯塞爾 1050)
土地利用變化能夠改變生物地理群落的結構和功能,從而影響碳循環(huán)過程,因此成為造成陸地生態(tài)系統(tǒng)碳儲量變化關鍵要素之一[1-2]。除了燃燒化石燃料之外,土地利用變化對大氣中CO2濃度增加的影響最大[3-4]。由于人類活動對土地利用的干擾加劇,我國土地利用模式發(fā)生了較大變動,使得陸地生態(tài)系統(tǒng)的結構和功能發(fā)生變化,進而影響碳儲量,引發(fā)了重大環(huán)境效應。因此,從土地利用的角度研究碳排放,不僅可促進生態(tài)保護,還有利于耕地合理利用[5]?;春=?jīng)濟區(qū)是我國主要農(nóng)副產(chǎn)品基地與糧倉,近年來由于該地區(qū)人口不斷增長,經(jīng)濟快速發(fā)展,對土地資源尤其是建設用地資源的需求不斷增加[6-7],導致了過度的土地非農(nóng)化、城市“攤大餅”蔓延等問題,影響到區(qū)域碳收支的平衡,而碳收支平衡對一個區(qū)域的生態(tài)水平起著至關重要的作用[8]。因此,進行淮海經(jīng)濟區(qū)土地利用時空變化特征分析及其碳儲量測算研究,對該區(qū)域土地的合理利用十分必要[9]。
早在20世紀50年代國際上就有研究估算了全球土壤有機碳庫的儲量[10]。Rubey等[11]、Bohn等[12-13]、Post等[14]隨后分別計算了全球土壤有機碳庫儲量。90年代以來,不同研究結果差別較小[15-18]。國內(nèi)對土壤碳儲量的研究相對較晚,從20世紀80年代開始,對土壤有機碳庫的研究才日益增多,基于不同空間尺度對土壤碳儲量進行了系列的估算與分析[19-29]。關于植被碳儲量的計算,初期主要通過野外實地調(diào)查生物量和面積統(tǒng)計資料,但在野外調(diào)查時易受調(diào)查人員主觀判斷的影響,從而導致較大誤差[30];近期的研究普遍采用建立生物量和蓄積量間關系的方法來估算植被碳儲量[31]。
已有研究側(cè)重于歷年碳匯數(shù)據(jù)預測碳儲量,往往忽略了外部因素的變化對碳儲量的影響[32]。因此,為掌握淮海經(jīng)濟區(qū)碳時空特征,本研究擬從土地利用的觀點出發(fā),基于遙感技術的應用和土地利用數(shù)據(jù),分析2006、2011和2017年3 個時期土地利用變化,并依據(jù)CA-Markov模型模擬預測淮海經(jīng)濟區(qū)2025年的土地利用分類情況,進而對因其變化引起的碳儲量變化進行分析,以期預測未來土地利用變化造成的碳收支變化,也為政策制定提供定量依據(jù)。
淮海經(jīng)濟區(qū)位于32°22′43″~36°33′42″ N,113°51′42″~120°55′44″ E,包括蘇豫魯皖4 個省20 個地級市147個縣級單位(2019年1月,國務院批復同意撤銷地級萊蕪市,將其所轄區(qū)域劃歸濟南市管轄,但由于本研究所涉及的數(shù)據(jù)未能根據(jù)國家相關規(guī)定實時更新,故仍將萊蕪市作為地級市進行研究),總面積達1.781×105km2,占全國土地總面積的1.86%(圖1)。研究區(qū)氣候以半濕潤季風氣候為主,溫和濕潤,雨熱同期,年降雨量為750~1 050 mm,年均氣溫在14 ℃左右,光熱條件好;植被以落葉闊葉林為主,土壤肥沃;主要作物包括小麥、稻谷、棉花、花生等,作物熟制為一年兩熟,輪作方式主要包括小麥-水稻、小麥-玉米等。
根據(jù)研究目的及研究區(qū)范圍,本研究中的遙感數(shù)據(jù)選擇2006、2011和2017年3 個時期MCD12Q1數(shù)據(jù),下載于美國地質(zhì)調(diào)查局網(wǎng)站(https:∥earthexplorer.usgs.gov)。
社會經(jīng)濟數(shù)據(jù)來源于2007、2012、2018年蘇豫魯皖統(tǒng)計年鑒;DEM數(shù)字高程數(shù)據(jù)從Google Earth網(wǎng)頁提??;碳密度數(shù)據(jù)來源于國家生態(tài)系統(tǒng)觀測研究網(wǎng)絡(http:www.cnern.org.cn)和參考文獻[23,33-37]的研究結果;道路矢量數(shù)據(jù)使用水經(jīng)注萬能地圖下載器獲取。
本研究結合《土地利用現(xiàn)狀分類》(GB/T 21010—2007)和淮海經(jīng)濟區(qū)實際情況,將研究區(qū)分為水域、林地、草地、耕地、建設用地、園地、裸地和其他用地等8 類(表1)。

圖1 淮海經(jīng)濟區(qū)區(qū)位圖Fig.1 Location map of the Huaihai Economic Zone

表1 土地利用分類體系Table 1 Land use classification system
20世紀50年代初,現(xiàn)代計算機創(chuàng)始人Neumann提出標準元胞自動機(cellular automaton,CA)的雛形,其模型原理:元胞自動機的空間、時間和狀態(tài)均是離散的,每個元胞也都處于離散狀態(tài);各個元胞狀態(tài)按照相同的轉(zhuǎn)變規(guī)則進行同步更新,每個變量狀態(tài)是有限的,且僅在局部遵循轉(zhuǎn)變規(guī)則[38]。Markov模型是基于概率創(chuàng)建的一種隨機型的時序模型,根據(jù)某一變量的目前形態(tài)和趨向預測其未來狀態(tài)的方法,因此,Markov模型可根據(jù)不同時期土地利用分類情況計算得出的轉(zhuǎn)置機率來預測未來某一時期的土地利用情況[30]。然而,Markov模型只能預測土地覆蓋分類數(shù)量方面的變化,不能預測其在空間方面的變化,即無法預測土地利用類型空間上的分布,元胞自動機恰好具備模擬空間演變這一功能,因此將二者結合組成CA-Markov模型,能實現(xiàn)在數(shù)量和空間上對土地利用變化進行模擬預測。
根據(jù)CA-Markov模型,結合前期研究區(qū)不同時期遙感圖像分類結果,計算得到未來某時期土地利用分類情況。本研究通過計算2017年土地利用現(xiàn)狀柵格數(shù)量和模擬預測柵格數(shù)量的誤差絕對值來驗證模擬精度,計算公式如下[39]:
(1)
式中:E為第i種土地利用類型的數(shù)量的差值;Mix、Miy分別為第i種土地利用類型的實際柵格數(shù)量和模擬預測柵格數(shù)量,個。需要注意的是運用此方法比較的是E的絕對值,其絕對值越小表示模擬精度越高。
本研究所用的CA-Markov預測流程如圖2。

圖2 CA-Markov預測流程圖Fig.2 CA-Markov prediction flow chart
考慮到碳密度計算所需的數(shù)據(jù)的可得性,本研究中的碳密度數(shù)據(jù)參考徐釗等[33]、賈松偉[34-35]、李偉等[36]、揣小偉等[37]已發(fā)表的研究結果;土壤碳密度計算,直接參考揣小偉等[37]研究結果,該研究有區(qū)域有5 個在淮海經(jīng)濟區(qū);植被碳密度的計算根據(jù)生物量碳密度估算法,即基于生物量推算出不同植被碳密度,通過對淮海經(jīng)濟區(qū)所涉及的蘇、豫、魯、皖四省的主要作物面積及產(chǎn)量進行統(tǒng)計分析發(fā)現(xiàn)其結構相似,因此,河南、山東、安徽各土地利用類型碳密度可參考江蘇省的研究結果[37],各地類碳密度數(shù)據(jù)見表2。
利用碳密度計算不同地類的土壤碳儲量,計算公式如下[40]:
TOC=SOC+ZOC
(2)
(3)
(4)
式中:TOC為淮海經(jīng)濟區(qū)總碳儲量,Tg;SOC為土壤總碳儲量,Tg;Xi為第i種土地利用類型的面積,km2;SOCDi為第i種土地利用類型的土壤碳密度,kg/m2;ZOC為植被總碳儲量,Tg;ZOCDi為第i種土地利用類型的植被碳密度,kg/m2。

表2 淮海經(jīng)濟區(qū)不同地類碳密度Table 2 Carbon density of different land types in the Huaihai Economic Zone kg/m2
研究發(fā)現(xiàn)淮海經(jīng)濟區(qū)土地利用類型以耕地為主,在2006、2011、2017年所占比例分別為87.51%、87.10%和86.21%,呈現(xiàn)出緩慢減少的趨勢;建設用地面積持續(xù)增加,由2006年的9 697.94 km2增加到2017年的10 962.10 km2;草地不斷擴張,占地面積由5 786.50 km2擴張到6 177.75 km2;水域、林地和園地面積所占比例較少,呈緩慢增加趨勢;在此期間,其他土地利用類型面積不斷減少,2017年總面積僅為142.94 km2,所占比例不足1%。
為顯示不同土地類型之間的轉(zhuǎn)化的細節(jié),利用Origin2020制作?;鶊D(Sankey diagram)(圖3)。由圖3可知:研究區(qū)內(nèi)主要是耕地、草地、建設用地和園地之間的轉(zhuǎn)化;草地的轉(zhuǎn)出面積與轉(zhuǎn)入面積都是最大的,轉(zhuǎn)入面積大于轉(zhuǎn)出面積導致草地面積總體上有所增加,說明草地主要是空間分布的變化;從轉(zhuǎn)出方向看,草地轉(zhuǎn)為耕地的面積最大,占草地的31.74%,水域次之;從轉(zhuǎn)入方向看,耕地和水域是草地面積增加的主要原因,其中9.86%的耕地轉(zhuǎn)換為草地;建設用地轉(zhuǎn)出方向較多且面積相當,但由于2.49%耕地的轉(zhuǎn)入,建設用地面積在總量上有所增加,表明了淮海經(jīng)濟區(qū)正處于經(jīng)濟發(fā)展階段,同時也映射出在人口持續(xù)增長的壓力下對建設用地的需求不斷增加;裸地面積有所增加,原因在于研究區(qū)礦產(chǎn)資源豐富,煤、石油、天熱氣等儲量豐富,開采量不斷增加,造成區(qū)域內(nèi)裸地面積不斷增加;水域呈現(xiàn)持續(xù)擴張的趨勢,其原因是淮海經(jīng)濟區(qū)對生態(tài)環(huán)境保護的日益重視。

圖3 2006—2017年土地利用面積變化轉(zhuǎn)移Fig.3 Change and transfer of land use area from 2006 to 2017
3.2.1CA-Markov模型模擬精度驗證
首先對2017年土地利用變化進行模擬,并將其與2017年土地利用現(xiàn)狀進行比較以驗證模擬精度(表3)。
由表3可知,除了其他用地之外,各個地類面積預測誤差均小于5%,表明該模擬的準確程度較高,能夠較為精確的預測未來某一時間的土地利用情況。其他用地誤差較大的原因可能是其所占面積較小,2017年底其他用地面積僅占研究區(qū)域總面積的0.08%,基線值較低,且呈遞減趨勢,在后續(xù)2025年模擬預測中可忽略不計。

表3 2017年土地利用模擬精度Table 3 Land use simulation accuracy in 2017
3.2.22025年土地利用模擬
利用CA-Markov模塊制作淮海經(jīng)濟區(qū)2025年土地利用類型的適宜性圖集從而得到2025年土地利用預測數(shù)據(jù)(圖4)。
對預測結果進行統(tǒng)計,2025年淮海經(jīng)濟區(qū)水域、林地、草地、耕地、建設用地、園地、裸地、其他用地面積分別為2 019.30、195.09、6 209.72、151 958.73、11 863.31、1 566.64、4 804.71、87.04 km2,分別占研究區(qū)總面積的1.13%、0.11%、3.47%、85.02%、6.66%、0.88%、2.69%、0.05%。
3.3.1碳儲量測算
按照2.3土地利用數(shù)據(jù)和碳儲量計算公式,對淮海經(jīng)濟區(qū)2006、2011、2017年碳儲量計算結果如表4所示,對2025年碳儲量預測結果如表5所示。
3.3.2碳儲量時空演變特征
1) 碳密度和碳儲量時間變化特征
利用Origin繪圖功能制作2006—2025年研究區(qū)各土地利用類型碳密度和碳儲量隨著時間的變化情況,如圖5所示。

圖4 2025年土地利用預測數(shù)據(jù)Fig.4 Land use forecast data for 2025

表4 各土地利用類型碳儲量Table 4 Carbon density and carbon storage of land use types Tg

表5 2025年各土地利用類型預測碳儲量Table 5 Carbon density and carbon storage of land use types in 2025 Tg
由圖5可見:除耕地和其他用地碳儲量在研究期呈遞減趨勢外,其余地類呈不同程度遞增趨勢。2006—2025年,林地碳密度最大,變化也最為明顯,總體上呈勻速增加趨勢,結合這一期間林地面積數(shù)據(jù)得出林地碳密度變化主要原因在于林地面積的不斷擴張;水域、草地、耕地、建設用地碳密度幾乎沒有變化;園地僅次于林地的碳密度,總體上呈緩慢增加趨勢。除耕地外,建設用地碳儲量最大且在整個研究時期內(nèi)總體呈現(xiàn)不斷遞增趨勢;草地在2006—2011年內(nèi)增長幅度較為明顯,在2011—2025年期間內(nèi)增長速度有所減緩;耕地碳儲量總體呈現(xiàn)遞減趨勢,但減少速度極慢,是因為耕地面積在一定范圍內(nèi)有縮減;水域、林地、園地碳儲量變化較?。宦愕卦?006—2011年碳儲量增長幅度較大,在2011—2025年內(nèi)仍然呈現(xiàn)遞增趨勢,但增長速度大大減??;其他用地碳儲量明顯減少,原因在于其面積不斷被其他土地利用類型侵占。
2) 碳密度和碳儲量空間變化特征
基于ArcGIS繪制整個研究期間淮海經(jīng)濟區(qū)碳密度、碳儲量空間分布情況(圖6)。由圖(a)可知:在碳密度方面,2006—2017年碳密度增加區(qū)域在研究區(qū)內(nèi)零散分布,其變化范圍大多在0.99~3.97 kg/m2,原因在于草地向林地的轉(zhuǎn)換;在2017—2025年,碳密度變化總體上呈現(xiàn)出“多零散,少集聚”的分布,究其原因是濟寧市2018年開展實施建設國家森林城市,林地面積大幅度增加,和棗莊市、徐州市接壤處碳密度變化較大并且集中。由圖(b)可知:在碳儲量方面,2006—2011年,碳儲量總體上變化不明顯,臨沂南部增加0~9.03 Tg,濟寧市東南部、宿遷市和淮安市交界處碳儲量有所減少;2011—2017年碳儲量總體上減少,變化值處在9.08~107.26 Tg,這是由于城市化進程中建設用地需求不斷增加,侵占其他地類,導致該區(qū)域內(nèi)碳密度降低從而碳儲量減少;2017—2025年,碳儲量延續(xù)前期減少趨勢,但濟寧市與棗莊市、徐州市交界處碳儲量基本不變,原因在于該區(qū)域碳密度有所增加。
在上述結果的基礎上,利用ArcGIS的空間分析功能的標準差橢圓對2006—2025年碳儲量方向分布變化進行分析,生成橢圓面積表示范圍,橢圓扁率表示方向性的明顯程度。結果表示:整個研究期間,生成的橢圓集中在淮海經(jīng)濟區(qū)東部,表示碳儲量變化多集中在研究區(qū)中心城市徐州及其周邊地區(qū);橢圓面積呈現(xiàn)出先增后減的趨勢,究其原因在于2011—2017年淮海經(jīng)濟區(qū)經(jīng)濟高速發(fā)展,城鎮(zhèn)化速度不斷加快,地類之間的相互轉(zhuǎn)換涉及范圍較廣,導致碳儲量變化區(qū)域增大;2006—2025年橢圓扁率不斷減小,意味著碳儲量變化方向性越來越模糊,呈現(xiàn)出由棗莊、徐州、宿遷和淮安一帶向四周發(fā)散的趨勢。
本研究以土地利用為出發(fā)點,模擬計算了淮海經(jīng)濟區(qū)歷史時期的碳儲量變化,精度較高,并以此預測了未來一段時期的碳儲量變化,對其變化情況進行分析,得到的主要研究結論如下:
1)分析土地利用歷史時空變化可知,耕地為淮海經(jīng)濟區(qū)主要的土地利用類型,應對耕地進行二次劃分以提高研究精度。研究發(fā)現(xiàn)耕地非農(nóng)化趨勢明顯,究其原因是為了追求最大利潤,發(fā)展鄉(xiāng)鎮(zhèn)企業(yè)與“開發(fā)區(qū)熱”,耕地不斷被占用導致[41];建設用地不斷擴張,原因在于淮海經(jīng)濟區(qū)各市經(jīng)濟水平不斷提高,城鎮(zhèn)化進程不斷發(fā)展,大規(guī)模開發(fā),導致了建設用地面積不斷增加;水域、林地、草地、園地、裸地面積所占比例較少,總體上呈遞增趨勢;其他用地面積不斷減少。

圖5 各土地利用類型碳密度變化和碳儲量時間變化Fig.5 Time variation of carbon density and carbon storage by land use type

圖6 各土地利用類型碳密度變化和碳儲量空間變化Fig.6 Spatial changes in carbon density and carbon storage by land use type
2)利用CA-Markov模型模擬2025年淮海經(jīng)濟區(qū)土地利用取得了較好的效果;2017—2025年期間土地利用變化趨勢與2006—2017年保持一致:耕地持續(xù)非農(nóng)化;建設用地持續(xù)擴張;水域、林地、草地、園地、裸地總體上仍呈遞增趨勢;其他用地面積持續(xù)減少。
3)土地利用碳儲量測算結果顯示:時序上,2006、2011和2017年碳儲量分別為1 720.018、1 720.020、1 716.531 Tg,預計到2025年碳儲量會達到1 720.726 Tg,較2017年增加4.195 Tg。由于耕地有利于減緩碳密度和碳儲量的下降,雖然研究期間耕地面積不斷減少但其仍占很大比例;空間上,碳儲量變化總體上呈現(xiàn)集聚特征,其減少區(qū)域成片分布但其減少量較小,增加區(qū)域分布零散但增加量較大。
因此,淮海經(jīng)濟區(qū)土地利用方式的轉(zhuǎn)變導致碳儲量的變化,耕地在一定程度上減緩了碳儲量的減少,而且研究區(qū)碳儲量隨著林地面積的增加而增加。因此,應在切實保護耕地的基礎上,積極實施退耕還林,保護林地資源,既保證了糧食安全又有利于提高研究區(qū)生態(tài)環(huán)境質(zhì)量。