陆地生态系统是碳循环中的关键碳库，具有调节区域气候和促进碳循环的重要作用^[1]。土地利用通过改变陆地生态系统的结构与功能来影响碳循环过程。“十四五”规划明确指出，以城市群为主体协调推进区域经济发展和新型城镇化，促进中国城市绿色低碳可持续发展^[2]。因此，开展城市群土地利用变化与碳储量变化耦合关系研究，对于未来城市土地可持续利用与低碳发展具有重要的参考价值和借鉴意义。

近年来，各国经济快速发展，导致碳排放量急剧增长，研究如何通过调整土地利用方式、促进陆地生态系统固碳减排受到高度关注。目前，碳储量的测算方法主要有实地调查法和基于模型估算法。实地调查法虽能获得较为真实的碳储量数据，但对于调研设备要求高，数据获取难度大且处理时间长，只能满足局部地区要求需求。基于模型估算方法主要有CASA(carnegie-ames-stanford approach)模型^[3]、植被固碳模型^[4]、簿记模型^[5]、FORCCHN(forest ecosystem carbon budge tmodel for China)^[6]、InVEST模型(integrated valuation of ecosystem services and trade-offs)等^[7-8]。其中，InVEST模型数据处理操作简单、评估精度较高，满足长时间序列和大尺度区域研究需求，在生态系统碳储量研究中得到广泛的应用^[9]。

为了得出2035年碳储量变化趋势，则需要模拟出2035年土地利用空间分布。土地利用模拟模型主要有SD(system dynamics)^[10]、ANN-CA(artificial neural network-cellular automata)^[11]、CA-Markov^[12]、FLUS(future land use simulation model software)等^[13]，但是以上模型不能深入挖掘土地利用变化的驱动机制，对于多种地类的模拟精度不高。PLUS(patch-generating landuse simulation)模型综合考虑自然环境、交通区位和社会经济因素对土地利用变化的影响，时空动态模拟区域斑块级地类的相互作用，并计算出较高精度的未来土地利用的分布^[14-16]。丁岳等基于InVEST—PLUS模型测算环杭州湾碳储量，发现建设用地的快速扩张导致碳储量降低^[17]；孙欣欣等^[18]以南京市为例，得出碳储量空间分布与土地利用分布保持一致；孙方虎等^[19]则发现安徽省林地、草地等高含碳量地类得到保护，碳储能力提升。耦合PLUS—InVEST模型能够基于土地利用数据得出研究区域碳储量的变化情况，能较好地探讨区域土地利用变化与碳储量的相互影响机制。

珠三角城市群是中国最具发展潜力城市群之一，是“特色金融+现代物流+先进制造”产业集群。它作为南方地区对外开放的门户，辐射带动华南、华中和西南的发展，也是建设国家级森林城市群和海洋强省的重要基地。珠三角城市群人口吸纳度高，常住人口高达7.80×10⁷人，对各类用地需求极大。因此，本文选定珠三角城市群作为研究区域，探究2005—2020年土地利用变化及其对碳储量的影响，并对2035年珠三角的土地利用格局和碳储量进行模拟预测，以期为珠三角地区“双碳”目标下的国土空间规划和生态决策提供科学依据，加快实现珠三角新时代高质量发展，同时为中国其他城市提供引领和示范作用。

1 材料与方法 1.1 研究区概况与数据来源 1.1.1 研究区概况

珠三角位于广东省中南部，在112°45′—113°50′E，21°31′—23°10′N之间，总面积约为54 000 km²，其中丘陵、山地和岛屿约占总面积30%。下辖广州、佛山、肇庆、深圳、东莞、惠州、珠海、中山、江门9个城市，毗邻港澳，与东南亚地区隔海相望。珠三角属于南亚热带季风气候，光照充足，四季多雨，终年温暖湿润，年平均气温为20~23 ℃。中部是平原，北部、东部和西部分布着丘陵、山地和岛屿，河口向南海延伸，形成相对闭合的“三面环山，一面临海”的三角区独特地形地貌。研究区概况如图 1所示。

1.1.2 数据来源

(1) 土地利用数据。本文所需珠三角土地利用数据来源于中国科学院资源环境科学与数据中心(http://www.resdc.cn)，空间分辨率为30 m，并依据土地利用分类国标将地类划分为耕地、林地、草地、水域、建设用地和未利用地6大类。

(2) 土地利用变化驱动因子。结合已有研究和数据的可获得性，土地利用扩张主要是受到自然地理环境、交通路线发展情况和社会经济发展的综合影响。因此，本文从自然环境、社会经济和交通区位3方面来确定13个土地利用驱动因子。自然环境包括DEM、坡向、坡度、年均降水、年均气温。DEM来源于地理空间数据云平台(http://www.gscloud.cn/)，坡度和坡向数据经DEM提取而得，年均降水和年均气温来源于中国气象数据网(http://www.data.cam.cn/)。社会经济数据的人口和GDP数据来源于中国科学院资源环境科学与数据中心(http://www.resdc.cn)；交通区位数据包括到水域、高速、铁路、国道、省道、县道距离，数据来源于openstreetmap(http://www.openstreetmap.org)；选取的到河流、高速、铁路等距离通过ArcGIS10.8的欧氏距离分析计算而得。根据模型需要，将土地利用驱动因子栅格化后统一成与土地利用数据相同的投影坐标、行列号和空间分辨率。

1.2 研究方法 1.2.1 PLUS模型

PLUS模型是基于FLUS，SD和CA等模拟模型发展起来的未来土地利用模拟模型，能够动态地模拟多类土地斑块的生成和时空演化，深入挖掘各类土地利用变化的潜在机理^[20]。模型包括土地扩张分析策略(land expansion analysis strategy，LEAS)和基于多类随机斑块种子机制的CA模型(CA based on multiple random patch seeds，CARS)^[21]，它较以往的模型具有更高的模拟精度和更真实土地景观格局。

模型通过提取两期土地利用数据之间的扩张部分作为样本进行训练，采用随机森林算法对土地扩张部分和土地利用驱动因子进行计算，得到各土地利用类型的发展概率以及驱动因子在该时段内对各类用地扩张的贡献度。基于Markov得到未来各类用地需求，进而结合随机斑块的生成和转移矩阵的设置来确定未来土地利用空间分布格局。

(1) 土地扩张分析策略。土地扩张分析策略融合已有转化分析策略(transformation analysis strategy，TAS)和格局分析策略(pattern analysis strategy, PAS)的优势，克服了随着用地类型的增加而导致模型复杂性增加、灵活度降低等问题，使模型更深入挖掘未来诱发土地利用变化的因素，具有更好的解释性^[22]。该策略提取两期土地利用数据间各类用地的扩张部分，采用随机森林算法得到不同用地类型的增长概率和各驱动因子对未来不同用地类型扩张的贡献度。其表达式为^[23]：

$ P_{i, k}^d(x)=\frac{\sum\limits_{i=1}^n I\left[h_n(x)=d\right\rfloor}{M} $

(1)

式中：P^d_{i, k}(x)为空间单元i处k类土地利用类型增长概率; x为由驱动因子组成的向量; M为随机森林决策树数量; d取值1或0, 1表示地类i可以转至地类k, 0表示地类i不可转为地类k; h_n(x)为计算得到的土地利用预测类型; I〔h_n(x)=d〕为决策树的函数。

(2) 基于多类随机斑块种子机制的CA模型。基于多类随机斑块种子的机制的CA模型包含了基于多重土地利用的随机种子的斑块生成机制。在模拟过程中，土地利用需求将会通过自适应惯性机制去影响土地用地竞争，驱使土地利用的数量达到未来需求^[24]。为预测多类用地，CARS模块在传统的CA模型的基础上，结合随机种子生成和阈值递减机制来执行未来用地模拟的操作^[25]。其表达式为^[23]：

$ \mathrm{OP}_{i, k}^{d=1, t}=P_{i, k}^{d=1} \times \mathit{\varOmega}_{i, k}^t \times D_k^t $

(2)

式中：OP_{i, k}^{d=1, t}为用地类型k在空间单元i上生成的随机种子；P_{i, k}^d=1为用地类型k在空间单元i上的增长概率; Ω^t_{i, k}为空间单元i的领域效应; D_k^t为未来用地类型k需求的影响。

1.2.2 InVEST模型与碳密度数据

InVEST模型中碳储量模块基于各时期土地利用数据和各地类的碳密度数据，测算各时期的碳储量。碳储量模块将陆地生态系统的碳储量划分为地上、地下、土壤有机质和死亡有机质碳密度4个基本碳库，并且各地类的碳密度是一个常量^[26]。其表达式为：

$ C_i=C_{i-\text{above }}+C_{i-\text{below }}+C_{i-\text{dead }}+C_{i-\text{soil }} $

(3)

$ C_{i-\text{total }}=\sum\limits_{i=1}^n C_i \times A_i $

(4)

式中：C_i为地类i的碳储量; C_i-above为地类i的地上碳密度; C_i-below为地类i的地下碳密度; C_i-dead为地类i的死亡有机质碳密度; C_i-soil为地类i的土壤有机质碳密度; C_i-total为地类i的总碳储量; A_i为地类i的面积。因为死亡有机质碳储量含量极低且数据难以获取, 此处不予考虑^[27]。

基于已有研究结果可知，处于同一气候带的土地类型碳密度差异较小^[28]。因此，使用相似气候区用地类型碳密度进行碳汇情况分析，通过年均气温和年均降水对碳密度进行修正，经检验最终参照刘萍等^[29]、朱志强等^[30]测定的碳密度，以表 1作为碳密度标准。

2 结果与分析 2.1 土地利用变化分析

2005—2020年珠三角各时期土地利用类型面积及变化见表 2。研究期间，耕地、草地和水域分别减少1 403.90，84.03，916.93 km²，建设用地、林地和未利用地分别增加2 225.56，167.83，11.48 km²。其中，耕地面积减少最多，所占比例从2005年的30.21%下降至2020年的27.62%；建设用地增加最多，从2005年的8.45%增加至2020年的12.55%，15 a间面积增长明显。

2005年以来，珠三角土地利用类型转移较为明显(表 3)。研究期间，耕地的转出量最多，转出面积为3 838.53 km²，其中有45.99%转化为建设用地和45.15%转化为林地。林地主要转移为耕地和建设用地，表现出耕地和林地的双向转化和建设用地扩张逐渐侵占城镇周边林地的现象。建设用地转入量大，转出量小，转入量为2 258.88 km²。15 a内，新增建设用地中有78.16%的面积由耕地转化而来，可见其侵占耕地严重。建设用地侵占耕地的能力逐年减弱，侵占的速度也逐年放缓，说明城镇化进程中限制耕地转化为建设用地初见成效。水域的转出量位居第三，主要有919.76 km²转化为耕地和321.12 km²转化为建设用地。草地和未利用地的面积较小，转移为其他用地较少，变化不明显。

2005—2020年珠三角土地利用空间分布特征明显，空间变化显著，主要表现为耕地空间范围缩减和建设用地向外扩张(图 2)。研究期间，建设用地主要集聚在佛山的禅城区和顺德区，广州的越秀区、荔湾区、海珠区，东莞市的东城、南城区，深圳的宝安、龙华、龙岗、福田区等。2005年，这4个市的建设用地在空间上分布较为零散，到了2020年趋于片块分布，空间扩张较为明显，另外其余地区的建设用地也有不同程度的增长。说明随着珠三角经济建设的发展和社会人口的增长，城市间的交流联系加强，建设用地激增。耕地则主要分布在惠州市、江门市以及广州、佛山建设用地外围区域，耕地因城市建设的需要受到了明显的侵占。

林地空间分布广，主要分布在肇庆、惠州、江门等地，随着封山育林、退化林生态修复以及山水林田湖一体化治理等工程项目的实施，林地面积有所增加，东南部的林地增长最为明显。水域主要由珠江、西江、北江等水系和零散分布的坑塘水面构成。到2020年，佛山、中山以及西江、珠江沿岸的水域空间分布减少。一方面，受到围湖造田的影响，部分水域整理复垦为耕地；另一方面，由于城市建设和产业发展的需要，部分的坑塘水面和近海沿岸的水面转换为工业园区和港口码头等建设用地，以缓解建设用地占用耕地的紧张局面。

2.2 碳储量时空变化特征分析

基于InVEST模型的碳储量模块得出珠三角2005，2010，2015年和2020年的碳储量分别为3.32×10⁸，3.34×10⁸，3.36×10⁸，3.32×10⁸ t。2005—2015年，碳储量呈缓慢增长趋势，到了2020年，碳储量由前期的缓慢增长转变为下降趋势。在6类地类中，林地、耕地的碳储量所占比例较高，2020年林地碳储量所占比例高达81.85%；建设用地及其他地类碳储量较低。

由表 4可知，土地利用面积变化与碳储量变化表现出显著一致性。2005—2020年，耕地、草地、水域减少导致碳储量减少5.10×10⁷ t，林地、建设用地、未利用地的增多促使碳储量增长了4.82×10⁷ t。碳储量减少的主要原因是高含碳量的耕地面积减少，而碳储量增加的主要原因则是高含碳量的林地和较低含碳量的建设用地面积增长。水域面积减少相对较多，但碳储量变化不大，主要是因为水域的碳密度较低，水体固碳能力较弱。土地利用与碳储量耦合变化表明土地利用的结构变化对碳储量的含量与分布具有明显的影响。

由图 3可知，2005—2020年，碳储量的空间分布表现出一定的聚集性，呈现出“四周高中部低”的空间分布格局。碳储量高值区主要分布在珠三角东北部、西北部以及东南部，这些区域以高海拔高覆盖度的山区林地为主，集中了珠三角主要的林业资源，地形环境复杂，建设用地的扩张受到了限制。珠三角中部地区地势较为平坦，以建设用地为主，城市化程度高，人口密度大，水系丰富，交通发达，有利于人们进行社会生产活动，因而固碳能力相对较弱。因此，碳储量高值区和低值区的空间分布与土地利用类型在空间上的结构和布局具有较为明显的相关性。

2005—2020年期间，珠三角碳储量空间变化主要表现为市镇周围碳储量减少，郊区地带碳储量增加(图 4)。碳储量增加的区域主要集聚在佛山与广州交界处、广州与珠海交界处以及中山、惠州、江门等郊区，境内分布着大稠顶、北岭山、南昆山、北峰山国家森林公园和象头山国家级自然保护区等山林，是珠三角主要的林业资源。在退耕还林、封山育林等生态保护政策的驱动下，人们环境保护意识也随之增强，使得高碳密度的林地得到保护与增长，这是珠三角碳储量增加的主要原因。

碳储量减少区域零散分布在珠三角境内。随着三角工业化城镇化加快推进，城市规模不断扩大，城镇人口不断增加，各地建设用地均有不同程度的增加，占用耕地和水域的趋势增强，使得高碳密度地类减少，低碳密度增加，导致研究区碳密度总体上稍有下降。

2.3 土地利用变化模拟与碳储量预测 2.3.1 精度验证

本文以2010年土地利用数据作为训练样本，预测2020年土地空间分布。将实际的2020年土地利用数据和模拟而得的2020年土地利用数据一同输入模型精度验证模块，得出模拟的总体精度为0.87，kappa系数为0.80，表明构建的预测模型模拟精度较高，满足本文研究需求，可以用于模拟预测2035年珠三角土地利用。

2.3.2 驱动因子分析

通过PLUS模型的用地扩张分析模块得出各用地类型的驱动因子贡献度(图 5)。未来各地类的变化受到自然环境因素、社会经济因素和交通区位因素的综合影响，但影响程度各有不同。耕地、林地和建设用地主要受到GDP、人口、距离省道的距离和DEM等的影响；草地和未利用地主要受到DEM、坡度、气温降水等的影响；水域则受DEM，距水域的距离，人口以及气温的影响。

2.3.3 2035年土地利用模拟

土地利用变化受到人类活动的干扰，同时也受到自身自然资源条件以及地理空间位置的限制。基于土地利用扩张驱动因子贡献度，模拟得到珠三角2035年土地利用空间分布，通过ArcGIS10.8的栅格计算器分析得到未来15 a的用地扩张情况(图 6)。

未来15 a，珠三角还是以建设用地和林地的增长以及耕地的减少为主。到2035年，建设用地面积较2020年预计增加1 493.23 km²，增幅为18.68%，相比2005—2020年的扩张强度有所减缓。建设用地以佛山、广州和深圳中心城区的新增与扩展。林地面积增加152.64 km²，耕地面积减少1 531.29 km²。水域和未利用地减少，草地面积增加。

总体来看，未来珠三角依托《粤港澳大湾区发展规划纲要》，经济建设高质量发展，GDP随之增长，吸引外来人口参与到城市建设中，着力构建“六纵四横”的高速铁路网和“四纵两横”的普速高铁网，建设用地势必增长，区位优势明显的耕地作为新增建设用地的主要来源，耕地则相应的减少。在地势险要区域，开发利用成本高，限制因子难以消除，林地和草地则增长。

2.3.4 2035年碳储量预测

基于2035年土地利用空间分布格局预测2035年碳储量分布与变化情况(图 7)。2035年珠三角碳储量为3.90×10⁸ t，较2020年增加了5.75×10⁷ t，增长较为明显。结合2035年土地利用空间格局，碳储量空间分布与2005和2020年的分布差异不大，碳储量低值区还是集聚在珠三角中部地区，碳储量低值区分布在珠三角外围地带，与当年土地利用分布相对应。碳储量减少区域主要集聚在广州、佛山、深圳、中山等地的中心城区及其外围扩展区域；增加区域主要是珠三角边界城市，比如惠州、江门和肇庆等城市以及水系分布区域，与2035年新增地类表现出高度一致性，即建设用地增长区、耕地减少区为碳储量下降区，林地和水域增长区为碳储量增加区域。

未来15 a，随着粤港澳大湾区建设的协同推进，对建设用地需求持续加大，导致低碳储量用地增加，高含碳量用地进一步减少。但是，就目前形势看，珠三角积极响应国家号召，合力推进珠三角国家森林城市群建设进程，落实蓝色海湾战略部署，加快建成中心主城区人工林和森林公园，对污染退化天然林的修复保护力度加大，促进珠三角林地资源向好发展，增加了高含碳量的林地面积，提高森林覆盖率，依托生物、工程和技术等生态保护修复手段，提高生态系统固碳能力。

3 讨论与结论 3.1 结论

本研究通过PLUS和InVEST模型测算2005—2020年珠三角的土地利用变化及其碳储量空间格局分布，并预测分析2035年珠三角土地利用结构布局与碳储量变化趋势，得出以下结论。

(1) 2005—2020年珠三角碳储量与土地利用类型面积同增同减。15 a间，土地利用转移主要表现为耕地转化为建设用和林地。耕地和水域分别减少1 403.90，84.03 km²，导致碳储量减少4.71×10⁷ t，2.73×10⁶ t；建设用地和林地面积增长了2 225.56，167.83 km²，带动碳储量增长3.25×10⁷，1.56×10⁷ t；15 a间，珠三角碳储量呈缓慢上升后缓慢下降趋势，其中林地、耕地等为高含碳量地类，建设用地、水域和未利用地为低含碳量地类。

(2) 碳储量的空间分布与土地利用空间布局表现出显著一致性。2005—2020年珠三角碳储量空间分布有明显聚集性，呈现出“四周高中部低”的空间分布格局。碳储量高值区分布在珠三角东北部、西北部和东南部的林地和草地；低值区主要集中在广州、佛山、深圳等以建设用地为主的中部地区。

(3) 预计2035年，珠三角碳储量为3.90×10⁸ t，比2020年增加了5.75×10⁷ t。碳储量增加地区主要在未来新增林地和新增耕地，比如肇庆北部、江门西南部以及惠州东部地区；碳储量减少区域主要以广州、佛山、深圳等地的城市周边及其扩展区域。随着珠三角逐步建成国家森林城市群以及对林地保护修复力度加大，生态环境向好发展。

3.2 讨论

(1) PLUS模型和InVEST模型的优势与不足。本文虽选取13个驱动因子，但忽略了“三区三线”、国家自然保护区等政策限制因子对土地利用变化的约束，今后为提高模型模拟精度，应将限制性因素纳入驱动因子体系。碳密度数据主要从已有结果中来获取，与研究区真实的碳密度数据存在一定的差异，影响实际碳储量计算，今后的研究应将加入实地调查数据对碳密度数据进行修正。

(2) 土地利用变化对碳储量的影响机制。各用地类型面积的增减将会带动碳储量的增减，碳储量的空间分布与区域土地利用空间布局有显著一致性，说明土地利用的空间布局和结构变化在很大程度上影响着区域碳储量，可通过对土地利用类型的管控实现区域碳储量的走向。

(3) 未来的发展规划建议。在保证未来用地需求的前提下，加大对林地草地的修复保护力度，进一步挖掘林业资源固碳潜力。积极稳妥推进海岸线占补平衡制度，保证未来海洋强省的建设。严格实施耕地占补平衡制度和用途管控制度，严禁非农建设占用永久基本农田，还可以通过与邻近省份的耕地指标交易来缓解本区域的用地问题。