耕地利用轉型碳排放時空分異特征與形成機理研究

蓋兆雪 詹汶羲 王洪彥 杜國明

(東北農業大學公共管理與法學院， 哈爾濱 150030)

0 引言

耕地是人類賴以生存與發展的寶貴資源，是保障國家糧食安全和治國安邦的根本條件[1]。近年來，隨著工業化、城鎮化的快速發展，各類用地通過競爭相互消長，導致城鄉空間布局和城鄉地域形態等迅速轉變，鄉村要素快速流動引發了耕地邊際化、非農化和破碎化等問題日益明顯；同時耕地轉變為其他用地類型，耕地的數量和空間結構發生改變[2]。耕地是碳循環的重要環節，耕地利用轉型不可避免地改變碳排放規模[3-4]。耕地在轉型過程中既會增加植被、土壤和水域的固碳能力，又會釋放出大量碳，進而導致氣候異常引發全球變暖，威脅人類的生存和發展[5]。當前中國正處于鄉村振興戰略初期以及城鎮化、工業化加速階段，耕地利用和生態環境問題愈加顯著。因此，探究耕地利用轉型過程中的碳排放強度及時空分異特征，對促進社會經濟低碳轉型和可持續發展具有重要意義。

耕地利用轉型是由土地利用轉型拓展而來，可分為理論研究和應用研究兩方面，前者主要包括耕地利用轉型概念與內涵[6-7]、理論研究框架[8]、轉型機理[9]和轉型路徑[10]等方面；后者主要包括耕地利用轉型的時空分異特征[11]、耦合關系[12-13]以及轉型效應[14]等研究。碳排放立足于從微觀—中觀—宏觀尺度研究土地利用碳排放效應，在微觀尺度上，基于土壤碳和植被碳計算土地利用變化的碳通量[15]；在中觀尺度上，研究土地利用變化碳排放的時空特征[16]、影響因素[17]和碳排放動態模擬[18]等方面；繼而從宏觀尺度提出土地利用減排政策和碳補償辦法[19]。綜合來看，一是已有研究側重挖掘土地利用面積變化的碳排放問題，而耕地利用轉型對碳源、碳匯的影響研究尚未涉及；二是缺乏對耕地利用轉型的碳排放影響因素研究。因此，本文基于1990、2000、2010、2020年4期土地利用數據，采用單元網格法和碳排放系數分析1990—2020年研究區耕地利用轉型碳排放強度，結合重心分析、探索性空間數據分析和冷熱點分析工具揭示耕地利用轉型碳排放的空間格局演變規律和趨勢，并借助地理探測器模型探索耕地利用轉型碳排放的形成機理，以期為緩解生態環境惡化和完善耕地保護政策提供參考。

1 研究區概況

松花江流域哈爾濱段位于125°42′～130°10′E，44°04′～46°40′N(圖1)，地處松嫩平原，位于黑龍江省南部，是黑龍江省重要的糧食生產基地，是東北城鎮化進程和資源環境變化較快的典型區域，同時也是“東北振興”和“一帶一路”戰略的重要地區。研究區土地利用類型多樣，土地面積706 913.32 hm2。耕地面積比重較大，2020年耕地面積為421 697.41 hm2，年平均氣溫為5.34℃，冬夏最大溫差77.10℃，年均日照時數2 180.80 h，降水主要集中在6—9月，占全年的70%[20]。據相關研究表明，研究區城鎮化進程快、土地開發強度大、耕地非農化現象明顯[20]。

圖1 研究區地理位置Fig.1 Geographical location of study area

2 數據來源及處理

土地利用數據來源于中國科學院資源環境科學數據中心(http:∥www. resdc.cn)，空間分辨率為30 m，根據土地利用現狀分類標準，參照前人的研究[20]以及地理特征，將研究區土地利用類型分為耕地、林地、草地、水域、建設用地、未利用地。利用ENVI軟件進行波段融合、幾何校正、圖像增強、鑲嵌、裁剪等預處理，同時根據野外實地調查對其進行人工修正，獲得1990、2000、2010、2020年4期土地利用數據，最終建立1990—2020年4期土地利用數據庫。

耕地利用轉型碳排放是自然條件和社會經濟條件共同作用的結果，需綜合考慮自然條件和社會經濟因素，并兼顧指標選取的難易程度。因此，結合研究區實際情況，分別從自然、距離和社會經濟3方面選取影響因子，自然因子包括地形地貌因子和氣候因子，具體包含坡度X1、地形起伏度X2、高程X3、年平均降水量X4、年平均氣溫X5，其中坡度、地形起伏度和高程決定著土地利用類型的分布格局，坡度小、地形起伏度小、高程低的地帶，耕地極易與其他類型土地發生轉型。氣溫、降水條件決定著耕地、林地的地理分布和生產水平，是碳排放變化最敏感的因子。距離因子包括與城鎮中心距離X6、與鄉級以上道路距離X7、與水域距離X8，與城鎮中心距離越近，耕地越容易發生轉型，與鄉級以上道路、水域距離越近耕地發生轉型的機率也增加，如為運輸方便占用耕地而修建工業用地；為完善農田基礎設施，離水域越近的耕地易發生轉型。社會經濟因子包括地均GDPX9、人口密度X10和土地利用程度X11，表征人類活動對耕地利用轉型的干擾程度。具體指標見表1。

表1 形成機理指標體系Tab.1 Index system of formation mechanism analysis

3 研究方法

3.1 網格單元法

網格單元法可定量刻畫土地利用轉型碳排放時空的精細演變過程[26]。利用ArcGIS 10.2軟件，以正方形等積規則網格劃分樣本區域，綜合研究區面積、研究目的和計算機運算效率，經反復調試，以900 m×900 m的正方形網格為研究單元，獲得9 180個網格單元，并將所需研究數據轉入到對應的網格單元中。

3.2 碳排放測算

根據研究區土地利用類型，將耕地利用轉型分為耕地轉為林地、草地、水域、建設用地和未利用地5種類型。其中林地、草地、水域和未利用地作為碳匯，耕地、建設用地作為碳源。耕地利用轉型為相互轉型，需測算相互轉型的碳排放量，碳排放測算公式為

E=∑ei=∑Tiδi

(1)

式中E——碳總排放量

ei——耕地與其他土地利用類型相互轉型對應產生的碳排放量

Ti——耕地與其他土地利用類型相互轉型對應變化的土地面積

δi——耕地與其他土地利用類型相互轉型對應的碳排放系數

耕地碳排放系數的確定需要同時考慮碳排放和碳吸收兩方面(既產生CH4，又吸收CO2)，CAI等[27]研究表明，農作物碳排放系數為0.504 t/hm2；何勇等[28]研究表明，農作物碳吸收系數為0.007 t/hm2，由此得耕地凈碳排放系數為0.497 t/hm2。林地、草地的碳排放系數均來源于方精云等[29]的研究成果，中國森林碳匯效率的加權平均值為-0.581 t/hm2，草地排放系數為-0.021 t/hm2，考慮研究區的植被覆蓋和氣候情況，采用這一平均碳匯系數具有一定的合理性。水域碳排放系數也需要考慮碳源和碳匯兩方面(積水情況下是CO2的匯，被排干圍墾后是CO2的源)，賴力[30]的研究表明，中國水域的平均碳匯系數為-0.257 t/hm2，采取這一系數作為研究區水域碳排放系數。建設用地的碳排放系數既要考慮土壤的碳排放量又要考慮工業生產生活中的能源消耗等活動產生的碳排放，因此本文通過查閱統計年鑒、IPCC碳排放清單，結合研究區主要能源消耗(煤、石油、天然氣)結構和能源碳排放系數，確定建設用地碳排放系數為0.643 t/hm2。未利用地具有一定的碳吸收能力，碳排放系數參考相關研究[30]的公認值-0.005 t/hm2。

3.3 重心分析

依據重心模型構建耕地利用轉型碳排放重心模型[31]，揭示耕地利用轉型碳排放演化的空間軌跡。假定某區域由n個單元組成，其中，第i個子區單元內的中心坐標設定為(Xi,Yi)，Mi代表該區域單元在某種屬性特征意義中的數值，該屬性中的重心坐標表達式為

(2)

式中X——重心坐標的經度

Y——重心坐標的緯度

Xi——第i個網格中心的經度坐標值

Yi——第i個網格中心的緯度坐標值

Mi——第i個網格耕地利用轉型碳排放量

3.4 冷熱點分析

冷熱點分析用于確定空間聚集的高(熱點)/低(冷點)值區域[32]。通過該方法來反映空間上耕地利用轉型碳排放變化的聚類分布。冷熱點觀測值G用于描述冷熱分布，如果該指數為具有統計學意義的正值，則值越高熱點聚類的分布越集中，該區域為熱點區域，反之表明該區域為冷點區域。利用ArcGIS 10.2軟件中Getis-Ord General G工具進行冷熱點分析。冷熱點分布根據置信區間分為6個等級，即在99%、95%和90%置信區間的冷點與熱點分布。計算公式為

(3)

式中Mj——第j個網格的耕地利用轉型碳排放量

Wij——要素i和j的空間權重矩陣(值為1則表示空間相鄰，值為0則表示不相鄰)

3.5 地理探測器模型

地理探測器模型能夠定量評估影響研究區耕地利用轉型碳排放空間變化的因子及因子的影響強度[33 ]，計算公式為

(4)

式中q——形成機理解釋力，值域為[0,1]

L——因子的層數

Nh——第h層碳排放和因子對應的單元數

N——全區碳排放和因子對應的單元數

σ2——全區碳排放變化方差

交互作用探測主要為判斷各因子對因變量是獨立產生影響的還是相互作用后產生影響，影響的作用力是減弱還是增強，兩個因子之間的關系可分為以下幾類：①q(x1∩x2)Max(q(x1)，q(x2))，雙因子增強。④q(x1∩x2)=q(x1)+q(x2)，獨立。⑤q(x1∩x2)>q(x1)+q(x2)，非線性增強。

4 結果與分析

4.1 耕地利用轉型結構演化特征

利用ArcGIS 10.2平臺的疊加分析工具，分析1990—2020年研究區耕地利用轉型數量的分布特征及相互轉型的流向，揭示耕地利用轉型結構演化特征。1990—2000年耕地轉出面積為3 806.73 hm2，其中耕地轉為水域的面積最大，貢獻率為36.48%；其次是轉為建設用地，貢獻率為33.76%；另外，22.39%的耕地由于撂荒轉為未利用地。然而，耕地轉入面積為6 721.74 hm2，其中未利用地轉為耕地的面積最大，貢獻率為40.14%；其次是草地、林地，貢獻率分別為34.57%、19.46%。綜合來看，由于轉入量大于轉出量，1990—2000年耕地面積呈上升趨勢，增加了2 915.01 hm2，主要源于未利用地和草地的轉入(表2)。

表2 1990—2000年耕地利用轉型結果Tab.2 Transformation of utilization of cultivated land in 1990—2000

2000—2010年耕地轉型為其他類型土地的面積為44 131.50 hm2，其中轉為建設用地的面積最大，為21 381.48 hm2，貢獻率為48.45%；轉為林地、草地、水域和未利用地的面積相差不大，貢獻率在10%～16%之間。然而，其他類型土地轉為耕地的面積僅為27 506.61 hm2，主要來自建設用地(41.13%)和林地(24.75%)的轉入；其次為未利用地(17.45%)和草地(11.21%)的轉入。整體而言，由于轉出量大于轉入量，2000—2010年耕地面積呈下降趨勢，減少了16 624.89 hm2，主要流向建設用地和林地，可見耕地與建設用地、林地之間的轉型較為活躍(表3)。

表3 2000—2010年耕地利用轉型結果Tab.3 Transformation of utilization of cultivated land in 2000—2010

2010—2020年耕地轉為其他類型土地的面積為58 097.62 hm2，其中耕地轉為建設用地的面積最大，達到36 434.47 hm2，貢獻率為62.71%；其次是轉為林地和草地，貢獻率分別為16.65%、15.35%。而其他類型轉為耕地的面積為37 872.87 hm2，其中8 373.24 hm2的建設用地、8 695.29 hm2的水域、7 776.92 hm2的林地轉為耕地，貢獻率分別為22.11%、22.96%、20.53%，可見，水域轉為耕地的面積最大?？傮w來看，由于轉入量小于轉出量，耕地面積呈下降趨勢，減少了20 224.75 hm2，耕地主要轉為建設用地(表4)。

表4 2010—2020年耕地利用轉型結果Tab.4 Transformation of utilization of cultivated land in 2010—2020

4.2 耕地利用轉型碳排放時間變化特征

1990—2000年研究區耕地利用轉型碳排放量為3 704.12 t，此階段轉型過程中碳排放量大于碳吸收量，導致耕地利用轉型碳排放最終呈現為碳源形式。其中耕地與未利用地之間相互轉型所產生的碳源最大，碳排放量為1 336.70 t；其次是與草地之間的轉型，碳排放量為1 152.55 t。而耕地與水域之間的轉型最終以碳匯形式展現，碳排放量為-163.06 t。2000—2010年耕地利用轉型碳排放量是上一時間段的近5.87倍，碳排放量達到21 743.53 t，主要來自與建設用地之間的轉型，碳排放量達到19 370.90 t，占總碳排放量的89.09%；而林地、水域與耕地之間的轉型產生碳匯，碳排放量分別為-674.64、-731.53 t。2010—2020年耕地利用轉型碳排放量為35 656.29 t，是上一時間段的近1.64，碳源同樣主要來自與建設用地之間的相互轉型，碳排放量為27 588.86 t，占總碳排放量的77.37%；而耕地與林地之間的轉型最終形成碳匯，碳排放量為-1 751.00 t?？傮w來看，1990—2020年研究區耕地利用轉型碳排放量呈上升趨勢，由1990—2000年的3 704.12 t增加到2010—2020年的35 656.29 t，增加了近8.63倍，這與耕地非農化、寂寞化、邊際化等有關(表5)。

表5 1990—2020年耕地利用轉型碳排放量Tab.5 Carbon emissions from cultivated land use transformation in 1990—2020 t

4.3 耕地利用轉型碳排放空間分異特征

4.3.1耕地利用轉型碳排放重心分析

1990—2020年研究區耕地利用轉型碳排放重心基本保持穩定，向東移動了15.17 km(圖2)。1990—2020年耕地利用轉型碳排放重心向東遷移，其中1990—2010年重心移動距離最大，向東北方向移動了11.49 km，說明東北地區耕地利用轉型碳排放量顯著增加，而西南地區耕地利用轉型碳排放有所改善。2010—2020年重心發生改變，向東南方向遷移，但偏移距離不大，僅移動了3.68 km，可見東南地區碳排放惡化，而西北地區碳排放明顯改善。從市轄區來看，松北區耕地利用轉型碳排放量在減少，而道外區耕地利用轉型碳排放量增加，主要由于大量耕地轉為建設用地，導致碳排放量增加。

圖2 碳排放重心轉移軌跡Fig.2 Trajectory of carbon emission center of gravity transfer

4.3.2耕地利用轉型碳排放冷熱點分析

基于GeoDa平臺中的全局自相關工具探測研究區耕地利用轉型碳排放的空間集聚性，結果如圖3所示。1990—2000年、2000—2010年和2010—2020年研究區耕地利用轉型碳排放的Global Moran’sI分別為0.524、0.507、0.578，并通過P<0.001顯著性檢驗。Global Moran’sI整體呈上升趨勢，其中1990—2000年、2000—2010年Global Moran’sI基本保持穩定，略呈下降趨勢，減少了0.017，表明耕地利用轉型碳排放空間集聚性略有所下降，但下降幅度可以忽略不計。2000—2010年和2010—2020年Global Moran’sI變化較為明顯，增加了0.071，說明研究區耕地利用轉型碳排放空間集聚性增強?？傮w來看，研究區耕地利用轉型碳排放非隨機分布，具有較強的空間集聚性。

圖3 1990—2020年耕地利用轉型碳排放Moran散點圖Fig.3 Moran scatter charts of carbon emission from cultivated land use transformation in 1990—2020

為了進一步探測1990—2020年耕地利用轉型碳排放的空間異質性，利用冷熱點分析工具進行研究區耕地利用轉型碳排放局部空間自相關分析，以此揭示研究區耕地利用轉型碳排放的空間分布格局(圖4)。1990—2000年耕地利用轉型碳排放熱點區呈現點狀的分布格局，其中道里區、松北區點狀熱點區較多，表明這些區域碳排放量有所增加；冷點區域零星分布在阿城區，主要由于阿城區林地資源較多，退耕還林政策的實施改善了當地的碳排放量。2000—2010年耕地利用轉型碳排放熱點區呈現帶狀的分布格局，形成以中心區為圓心向四周擴散的趨勢，熱點區主要圍繞南崗區向周圍邊界擴散，這與研究區國民社會經濟發展方向相一致；冷點區仍然分布在東南部地區，但冷點區面積比1990—2000年增加，從點狀向片狀轉變。2010—2020年耕地利用轉型碳排放熱點區呈增加趨勢，出現了點狀、帶狀和片狀共存的空間格局，說明耕地利用轉型碳排放量進一步增加；而冷點區呈現點狀的分布格局，零星分布在研究區的東部和南部?？傮w來看，1990—2020年耕地利用轉型碳排放的冷熱點空間分布格局與社會經濟發展趨勢、三北防護林工程、退耕還林政策的實施密不可分。

圖4 耕地利用轉型碳排放冷熱點分布圖Fig.4 Cold and hot spot distribution diagrams of carbon emission from cultivated land use transformation

4.4 耕地利用轉型碳排放形成機理

由于不同時期影響因子對耕地利用轉型碳排放的影響具有相似性，因此以2010—2020年為目標年進行耕地利用轉型碳排放形成機理分析。利用地理探測器模型中的因子探測法定量評估影響研究區耕地利用轉型碳排放的因子及因子的影響強度(圖5)。2010—2020年耕地利用轉型碳排放的影響因子作用強度從大到小依次為與城鎮中心距離(0.107 8)、土地利用程度(0.103 0)、年平均降水量(0.045 7)、高程(0.039 3)、年平均氣溫(0.033 3)、與鄉級以上道路距離(0.029 4)、地均GDP(0.028 1)、人口密度(0.027 2)、地形起伏度(0.021 0)、坡度(0.018 9)、與水域距離(0.014 1)，其中與城鎮中心距離是耕地利用轉型碳排放的主控因子，土地利用程度也具有較強的解釋力；而地形起伏度、坡度、與水域距離對耕地利用轉型碳排放的影響強度較低，說明與城鎮中心距離、土地利用程度對耕地利用轉型的碳排放量有較大的影響。隨著研究區城鎮化、工業化進程的加快，建設用地的經濟收益明顯高于耕地的經濟收益，耕地在適宜的條件下極易發生轉型，同時人類生活及生產活動需求的增加，不斷促進耕地非農化，進而導致耕地利用轉型碳排放增加。

圖5 2010—2020年影響因子的形成機理解釋力Fig.5 Explanatory power of formation mechanism of impact factor in 2010—2020

2010—2020年耕地利用轉型碳排放的空間異質性源于多種因子的共同作用，單一因子并不能完全解釋耕地利用轉型碳排放的空間分異特征。根據地理探測器模型的交互探測法可知，不同因子之間的相互作用會產生更強的作用效果，其作用方式包括非線性增強和雙因子增強兩種。耕地利用轉型碳排放的影響因子交互作用以雙因子增強為主。其中土地利用程度與年平均降水量(0.166 3)的相互作用最大；其次為與城鎮中心距離與土地利用程度(0.161 4)、土地利用程度與高程(0.148 7)。其中土地利用程度與年平均降水量(0.166 3)、土地利用程度與高程(0.148 7)，因子交互作用呈現非線性增強，而與城鎮中心距離與土地利用程度，因子交互作用呈現雙因子增強效應，可見社會經濟因子與自然因子之間的交互作用以非線性增強為主?？傮w來看，土地利用程度與其他因子之間的交互作用明顯強于其他因子之間的交互作用，土地利用程度是人類活動作用的結果，對耕地利用轉型碳排放的影響強度較強，因此土地利用程度與其他因子之間的作用強度最為復雜。社會經濟因子與距離因子、自然因子之間的交互作用明顯強于內部因子的交互作用，表明耕地利用轉型碳排放空間分異特征并不是單獨發生作用，而是呈現協同增強的作用效果(表6)。

表6 2010—2020年各影響因子的交互作用(q值)Tab.6 Interaction of influencing factors in 2010—2020

5 討論

本文利用網格單元法，基于長時間序列土地利用現狀數據揭示了耕地利用轉型碳排放的時空分異特征，從自然因子、社會經濟因子和距離因子3方面挖掘耕地利用轉型碳排放的影響因子，探測不同因子以及因子之間交互作用的影響強度，探索耕地利用轉型碳排放空間分布格局的形成機理，在一定程度上為耕地可持續利用奠定了基礎。

通過研究發現，1990—2020年松花江流域哈爾濱段耕地利用轉型碳排放量明顯增加，歸因于研究區耕地利用轉型比較頻繁，且主要轉型為碳排放較高的建設用地。2010—2020年耕地利用轉型最為劇烈，耕地利用轉型碳排放量是1990—2000年的9.63倍，碳源量遠大于碳匯量。與城鎮中心距離、土地利用程度對耕地利用轉型碳排放的影響較強，與城鎮中心越近，人類活動越活躍，耕地發生轉型的概率就會增加，且極易轉為建設用地[34]，進而導致碳排放量增加。從因子交互作用結果來看，土地利用程度與年平均降水量、與城鎮中心距離的交互作用對耕地利用轉型碳排放影響程度明顯增強。因此，碳排放的控制需考慮社會經濟因素和自然因素的共同作用。

此外，本文僅從自然、社會經濟和距離因子3方面探究耕地利用轉型碳排放形成機理，但缺乏對政策因子、土壤類型、新型經營主體等方面的考慮。下一步擬結合土地利用政策、土壤類型等方面開展耕地利用轉型碳排放形成機理分析，并開展碳分區補償研究。

6 結論

(1)1990—2000年研究區耕地面積呈上升趨勢，源于未利用地和草地的轉入；2000—2010年、2010—2020年耕地面積均呈下降趨勢，耕地與建設用地、林地之間的轉型最為頻繁，耕地主要轉型為建設用地和林地。

(2)1990—2020年研究區耕地利用轉型碳排放量呈上升趨勢，其中2010—2020年耕地利用轉型碳排放最為劇烈，耕地利用轉型碳排放最終呈現為碳源形式。耕地利用轉型碳排放重心持續向東移動，其中1990—2010年重心移動距離最大，呈現東北地區碳排放惡化，而西南地區碳排放明顯改善的特點。

(3)1990—2020年研究區耕地利用轉型碳排放非隨機分布，具有較強的空間集聚性。熱點區主要圍繞南崗區向周圍邊界擴散，冷點區零星點狀分布在東南部地區，耕地利用轉型碳排放的冷熱點空間分布格局與社會經濟發展趨勢、三北防護林工程、退耕還林政策的實施密不可分。

(4)單一因子對耕地利用轉型碳排放空間分異的影響存在一定差異。與城鎮中心距離是耕地利用轉型碳排放的主控因子，其次是土地利用程度。雙因子交互作用對耕地利用轉型碳排放的影響強度明顯高于單一因子。各因子之間交互作用以雙因子增強為主，其中土地利用程度與年平均降水量、與城鎮中心距離的相互作用的解釋力相對較高，社會經濟因子與距離因子、自然因子之間的交互作用明顯強于內部因子的交互作用。