耕地利用转型碳排放时空分异特征与形成机理研究

盖兆雪 詹汶羲 王洪彦 杜国明

(东北农业大学公共管理与法学院， 哈尔滨 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%，可见，水域转为耕地的面积最大。总体来看，由于转入量小于转出量，耕地面积呈下降趋势，减少了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。总体来看，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，说明研究区耕地利用转型碳排放空间集聚性增强。总体来看，研究区耕地利用转型碳排放非随机分布，具有较强的空间集聚性。

图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年耕地利用转型碳排放热点区呈增加趋势，出现了点状、带状和片状共存的空间格局，说明耕地利用转型碳排放量进一步增加；而冷点区呈现点状的分布格局，零星分布在研究区的东部和南部。总体来看，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)，因子交互作用呈现非线性增强，而与城镇中心距离与土地利用程度，因子交互作用呈现双因子增强效应，可见社会经济因子与自然因子之间的交互作用以非线性增强为主。总体来看，土地利用程度与其他因子之间的交互作用明显强于其他因子之间的交互作用，土地利用程度是人类活动作用的结果，对耕地利用转型碳排放的影响强度较强，因此土地利用程度与其他因子之间的作用强度最为复杂。社会经济因子与距离因子、自然因子之间的交互作用明显强于内部因子的交互作用，表明耕地利用转型碳排放空间分异特征并不是单独发生作用，而是呈现协同增强的作用效果(表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)单一因子对耕地利用转型碳排放空间分异的影响存在一定差异。与城镇中心距离是耕地利用转型碳排放的主控因子，其次是土地利用程度。双因子交互作用对耕地利用转型碳排放的影响强度明显高于单一因子。各因子之间交互作用以双因子增强为主，其中土地利用程度与年平均降水量、与城镇中心距离的相互作用的解释力相对较高，社会经济因子与距离因子、自然因子之间的交互作用明显强于内部因子的交互作用。