碳储量是生态系统服务功能中的重要环节，植被和土壤是陆地生态系统最重要的两大碳库，其固碳功能在缓解气候危机上发挥着重要作用^[1-3]。碳储量变化受到国际科学联合会(ICSU)、世界气象组织(WMO)和联合国环境规划署(UNEP)等多个组织的高度关注。国外学者对碳储量研究较早，20世纪末Simpson^[4]核算出北美东部和北部陆地针叶林碳储量以及不同植被类型的固碳能力，Hurtt等^[5]采用森林资源调查的方法研究区域碳储量变化特点。近年来对于大区域的碳储量研究更多的使用模型核算，如李传华^[6]采用改进的CASA模型计算石羊河流域的植被净第一生产力(NPP)，揭示影响植被碳储量变化的驱动因子。相较于CASA模型只能计算植被碳储量，InVEST(integrated valuation of ecosystem services and tradeoffs)模型可评估包含植被地上地下、土壤和死亡有机质碳储量^[7]。该模型是由美国斯坦福大学、世界基金会(WWF)和大自然保护协会(TNC)共同开发，2010年后国内学者使用InVEST模型研究区域碳储量变化，在白龙江、太行山等区域都取得较好的科研成果^[8-9]。荣月静^[10]基于“全国生态环境遥感调查”的数据在InVEST模型中核算太湖流域碳储量演变，指出林地和湿地转入建设用地是太湖生态系统碳储量下降的主要原因。

土地利用变化是造成生态系统碳储量变化的重要原因，导致大量的碳从陆地生态系统流向大气生态系统^[11]。城市建设用地扩张通过侵占林地、耕地和草地等引起的土地利用变化，使得城市生态系统面临严峻的碳流失问题。城市土地利用变化对生态系统的影响是近年来生态学领域研究的热点，陆地碳储量作为生态系统的重要环节也是众多学者关注的目标^[12]。目前关于碳储量研究大多集中在碳储量总量变化上的分析，鲜有模拟未来城市土地利用与碳储量的相关研究。土地利用变化的预测可以分为数量预测和空间预测，现有模型在模拟中大多未考虑多种土地类型在转换中的相互影响，忽略了土地类型的竞争关系。FLUS(future land use simulation model)模型基于神经网络训练的土地适应性概率分布，提出自适应惯性竞争机制，能较好的模拟未来多种土地类型的分布^[13-14]。已有多位学者使用FLUS模型在预测未来城市土地变化方面的研究，如刘小平等^[15]采用该模型预测珠三角城市群的开发边界，林伊琳^[16]预测昆明市未来建设用地的扩张和景观格局变化。最近的一些研究中也试图在将未来土地变化与碳储量模型链接起来作为评估土地变化对未来生态系统的影响，如Zhao^[17]使用CA-Markov耦合InVEST模型评估生态工程对中国西北黑河流域碳储量影响，预测实施生态工程后，2015-2029年该区域碳储量可增加1.00×10⁷ t，证明耦合模型在评估未来碳储量研究中具有一定科学价值。对于未来城市土地预测已有相关研究，但对于链接未来城市土地变化与碳储量的研究尚少。探讨城市建设用地扩张下土地利用变化对碳储量的影响，揭示碳储量时空演变和未来空间分布趋势，能为城市规划和生态脆弱区实施精准保护提供科学参考。

1 材料与方法 1.1 研究区概况与数据介绍

广州市地处珠江下游，濒临南海，位于东经112°-114°，北纬22°-24°之间。属于丘陵地貌，东北高西南低，北部植被覆盖密集，南部为冲击平原。受海洋性亚热带季风气候影响，光热充足，降水丰富，境内河流众多，水域面积广阔。广州是粤港澳大湾区和一带一路枢纽城市，我国重要的经贸中心。该市下辖11个区，总面积约7 400 km²，2018年常住人口达到1 500万，城镇化快速发展使得建设用地逐年扩增。研究所需要的数据包括土地利用数据和土地驱动因子数据，土地利用数据：经过遥感解译得到1990-2018年5期土地利用数据，分为林地、耕地、草地、水域、建设用地和未利用地。用于土地预测的驱动因子数据有：气象数据，包括降水和气温；地形数据，包括DEM(数字高程)和坡度；社会经济数据，包括行政区划矢量、交通网络、人口密度和GDP；NDVI(归一化植被指数)数据。

1.2 研究方法 1.2.1 InVEST模型

InVEST模型旨在权衡土地利用与生态系统服务功能的关系，目前已经发展到包含水源涵养、生境质量、碳存储等多个模块，形成多种生态系统服务评估功能为一体的重要模型^[18]，其中碳储存服务模块在美国、南美洲、非洲、东南亚等全球多区域成功运用到实际研究中。生态系统的碳储量包括植被地上碳储量、植被地下碳储量、土壤碳储量和死亡有机质碳储量^[19]，即：

$ {C_{{\text{total}}}} = {C_{{\text{above}}}} + {C_{{\text{below}}}} + {C_{{\text{soil}}}} + {C_{{\text{dead}}}} $

(1)

式中：C_toatl表示总体碳储量；C_above表示植被地上碳储量；C_below表示植被地下碳储量；C_soil表示土壤碳储量；C_dead表示死亡有机质碳储量。由于死亡有机质碳储量含量非常低, 一般不考虑^[20]。

InVEST模型需要输入研究区土地利用数据和各土地类型对应的碳密度值(即“碳密度表”)，通常是借助文献查阅得到研究区的碳密度数据。前人研究中认为处于同一气候带的土地类型碳密度差异较小^[21-22]，以荣检^[23]、朱鹏飞^[24]研究中得到的广西地区土地类型碳密度成果为主要参考，在土地类型碳密度选择上，采用降水和气温因子通过公式修正本地化土地类型碳密度数据，反演得到广州市的碳储量。修正因子计算公式^[25-26]为：

$ {K_{{\text{BP}}}} = \frac{{{{C'}_{{\text{BP}}}}}}{{{{C''}_{{\text{BP}}}}}};{K_{{\text{BT}}}} = \frac{{{{C'}_{{\text{BT}}}}}}{{{{C''}_{{\text{BT}}}}}} $

(2)

$ {K_{\text{B}}} = {K_{{\text{BP}}}}\cdot{K_{{\text{BT}}}} = \frac{{{{C'}_{{\text{BP}}}}}}{{{{C''}_{{\text{BP}}}}}}*\frac{{{{C'}_{{\text{BT}}}}}}{{{{C''}_{{\text{BT}}}}}} $

(3)

$ {K_S} = {{C'}_{{\text{SP}}}}/{{C''}_{{\text{SP}}}} $

(4)

式中：K_BP，K_BT分别表示植被碳密度降水因子和气温因子修正系数；K_B表示地上地下植被碳密度修正系数；K_S表示土壤碳密度修正系数；C′，C″分别表示广州和广西地区的碳密度, 由年均温和年降水量代入公式计算得到：

$ {C_{{\text{BP}}}} = 6.798{{\text{e}}^{0.005{\text{ }}4{\text{MAP}}}} $

(5)

$ {C_{{\text{BT}}}} = 28{\text{MAT}} + 398 $

(6)

$ {C_{{\text{SP}}}} = 3.396{\text{ }}8{\text{MAP}} + 3{\text{ }}996.1 $

(7)

式中：MAP表示降水；MAT表示气温；C_BP，C_BT分别表示根据降水量和气温得到的植被碳密度；C_SP表示根据降水量得到的土壤碳密度。综合所得的各土地类型的碳密度值详见表 1。

1.2.2 FLUS未来土地预测模型

FLUS模型是在元胞自动机原理的基础上耦合马尔科夫链和神经网络模型预测未来土地利用，在预测中考虑到人文因素和自然因素对土地变化的多重影响，可以模拟多类用地的演变过程，对土地类型变化的预测结果更加精确。模型主要分为两部分：第1部分基于ANN网络的样本训练得到土地适宜性概率，ANN包括输入层、隐藏层和输出层，通过反向传播算法不断更新各层之间的权重系数来达到最优估计目的，第2部分为自适应惯性和竞争机制，核心是在迭代中自适应调整惯性系数〔公式(8)〕，最终模拟结果不仅取决于神经网络得到的土地分布概率，还受到邻域和转换成本等限制，体现出土地类型变化中相互作用的竞争模式^[27]。

$ I_k^t = \left\{ \begin{gathered} I_k^{^{t - 1}}\left( {if|{D^{t - 2}}\left| \leqslant \right|{D^{t - 1}}|} \right) \hfill \\ I_k^{^{t - 1}} \times \frac{{{D^{t - 2}}}}{{{D^{t - 1}}}}(if0 > {D^{t - 2}} > {D^{t - 1}}) \hfill \\ I_k^{^{t - 1}} \times \frac{{{D^{t - 1}}}}{{{D^{t - 2}}}}(if{D^t} - 1 > {D^{t - 2}} > 0) \hfill \\ \end{gathered} \right. $

(8)

式中：I_k^t为迭代次数为t时的惯性系数；D_k^t-1为到t-1时预测的某土地类型与真实土地类型的像元个数差值。土地利用变化受到人文和自然多重因素的共同影响，结合研究区特点和数据的可获取性，本次研究中选择坡度，DEM，年均降水，NDVI，年均气温，到市中心距离，交通网络，GDP和人口密度作为驱动因子，具体数据来源详见表 2。

模型中设置的成本矩阵表示当前土地类型转为其他用地的难易程度，参考自然发展条件下城市土地转移的规则，通过探索性试验获得成本矩阵。根据研究区的实际情况及土地转移概率矩阵，设置邻域因子权重参数范围为0~1，越接近1表示该地类的扩张能力越强^[28]，参考历史土地转移特征，经过调试和验证，得到模拟精度较高的邻域因子参数表(表 3)。依据上述步骤，分别计算出每个栅格的总概率，将预测的土地类型分配到栅格，计算公式为：

$ T_{_{p,k}}^t = s{p_{p,k}}\cdot\Omega _{_{p,k}}^t \cdot I_k^t\cdot(1 - sc) $

(9)

式中：T_{p, k}^t表示像元p在迭代t时转成土地类型k的概率；sp_{p, k}是土地适宜性概率分布；Ω_{p, k}^t为邻域作用；(1-sc)表示土地类型转换的难易程度。

将所有驱动因子重采样到相同分辨率，进行标准化处理，神经网络隐藏层设为12，先通过神经网络样本训练得到研究区土地类型适宜性概率分布，标准偏差为0.226，再基于土地类型适宜性概率分布，在自适应惯性和竞争机制模型下预测未来土地分布。模型的验证以2010年土地类型作为训练集，预测2018年土地类型分布，总体精度94.89%，Kappa系数为93.93%，表明模型预测和真实地物分布有较高的一致性。此外，通过Precision(精确率)、Recall(召回率)和F₁综合分数进一步评价模型的预测表现，指标的计算见公式(10)-(11)，各土地类型精度的综合评价结果详见表 4。综合以上评价结果，表明模型预测精度达到要求，可作为预测研究区未来土地利用变化。

$ {\text{Precision}} = \frac{{{t_p}}}{{{t_p} + {f_p}}}{\text{, Recall}} = \frac{{{t_p}}}{{{t_p} + {f_n}}} $

(10)

$ {F_1} = 2 \times \frac{{{\rm{Preision}} \times {\rm{Recall}}}}{{{\rm{Preision + Recall}}}} $

(11)

式中：t_p表示真阳性记录百分比；f_n表示假阴性记录百分比；f_p表示假阳性记录百分比。

2 结果与分析 2.1 广州市1990-2018年土地利用与碳储量变化 2.1.1 土地利用变化

表 5-6分别为研究期内广州市土地类型面积和土地类型动态度。从表 5可知，林地是广州市的最大景观类型，面积所占比例在42%以上，是研究区的优势土地景观类型。其他土地类型面积所占比例依次是耕地、建设用地、水域、草地、未利用地。1990-2018年土地类型变化中，耕地面积的减少量最高，1990年耕地面积所占比例接近40%，2018年所占比例下降到30.15%，总共减少681.0 km²。其次是林地，减少158.5 km²，草地减少11.1 km²，未利用地减少0.9 km²，水域增加59.84 km²。面积增幅最大的是建设用地，增加791.8 km²，年均增加40.02 km²，面积所占比例从8.60%上升到19.64%，城市扩张特征明显。建设用地是变化量最高的土地类型，其次是耕地和林地，其他地类变化量相对较小。

土地利用动态度表示研究区某一定时间范围的土地类型变化情况，用来反映区域土地变化的幅度和速度，通过各时期动态度变化来研究土地利用变化特点。

表 6是土地类型面积和动态度变化，1990-2005年土地利用动态度较高，土地变化剧烈，2005-2018年土地利用的动态度逐期下降，土地变化相对缓和。耕地在2005-2010年和2010-2018年减少量分别占2000-2005年减少量的36.69%和15.67%，研究时期的动态度分别是-0.82%，-2.25%，-0.93%和-0.42%，2005年后变化有明显下降趋势。1990-2018年林地动态度为-0.18%，是所有土地类型动态度最低的地类，由于林地面积所占比例最高，面积减少量仍然较高，仅次于耕地。建设用地各时期的动态度最高，分别是3.11%，9.20%，2.64%和1.08%，1990-2010年建设用地的增量较大，2010-2018年扩张趋于稳定，增量减少，增加量趋于稳定，只占2000-2005年的19.36%，建设用地迅速扩张且扩张能力明显下降。草地在2000-2010年减少7.78 km²，受到退耕还草政策的实施效果，2010-2018年草地面积增加1.27 km²。1990-2000年水域变化明显，面积增加81.27 km²，未利用地的面积变化较小。1990-2018年的土地变化情况来看，建设用地的动态度远高于其他地类，研究时间内广州市土地类型变化的主要特征是建设用地面积的快速增加，耕地和林地大量转入建设用地。

2.1.2 碳储量变化

通过查阅相关文献以及公式修正获得广州市土地类型碳密度数据，结合当年土地利用数据导入InVEST模型碳模块下运行，得到1990-2018年5期碳储量数据及其变化。图 1是根据模型运算得到碳储量结果。结果表明1990-2018年广州市碳储量减少2.48×10⁶ t，降幅6.2%，年均下降8.86×10⁴ t。其中2000-2005年的下降幅度最大，年均下降2.18×10⁵ t，2010-2018年下降幅度为0.8%，年均下降3.74×10⁴ t。1990-2000年广州市建设用地扩张能力较强，进入21世纪后，广州市城市建设进入一个新的高峰，随着珠江新城的开发和周边区域的城市扩张，2000-2005年碳流失达到高峰。总体上1990到2010年区域碳储量的变化比较剧烈，该时期广州市经济快速增长，城镇化速度加快，对于土地开发需求也较为强烈。2010年后建设用地扩张趋于缓和，土地变化逐渐趋于平稳，这一时期广州市严重碳流失逐步得到缓解。

2.1.3 碳储量空间变化特征

从碳储量空间分布来看，广州市碳储量空间分布具有显著的空间异质性，图 2分别是1990，2005和2018年陆地碳储量分布情况。由图 2可知，高密度碳储量主要分布在东北部一带，该区域海拔较高，拥有高覆盖的森林面积，植被覆盖率高；低密度碳储量主要分布在珠江下游，该区域主要为平原，城市化程度较高；南部主要覆被农田等，受人类活动影响较大，碳储量维持在较低水平。东北部海拔较高，主要是山区林地覆被，碳储量较高。高海拔决定了坡度和坡向，限制了城市用地的扩张和开垦活动，在一定程度上影响植被类型的分布和土壤的性质，且在水源涵养、森林资源等方面更容易占据优势的生态位；南部属于冲击平原，地势平坦，水系丰富，适合人类社会生产活动，因而固碳能力相对较弱。

图 3为1990-2005年和2005-2018年陆地碳储量空间变化情况。从碳储量空间变化来看，变化的区域具有大集聚和零星分布的特点。1990-2005年碳储量显著下降的区域分布在花都区、番禺区和黄埔区，主要集中在东南和西部区域，这一时期建设用地扩张剧烈，番禺区、海珠区、黄埔区和花都区扩张明显，大量耕地和林地转为建设用地，靠近城市中心附近的碳储量显著减少。2005年后建设用地扩张能力相对下降，土地类型转移趋于稳定，区域碳储量下降也逐渐缓和，受城市东扩影响，碳储量显著大面积下降的区域主要是黄埔区和增城区西南，城市中心区域的周边也有零星的减少，如白云区和花都区等。东北部区域的从化区和增城区北部碳储量较为稳定，该区域覆被大量森林资源，城镇扩张能力较低，碳储量相对稳定，珠江下游区域受城市发展影响，碳储量变动较大。总体上广州市碳储量在2005年之前有明显的下降，2005年后下降幅度低于前期，2005-2018年碳储量总体变化较之前缓和。

2.2 未来建设用地扩张与陆地碳储量变化 2.2.1 未来土地预测与碳储量潜力

从预测的2034年土地分布来看，建设用地会进一步扩张。从2018年1 409.11 km²增加到2034年的1 683.48 km²，面积增加274.37 km²，年均增长17.14 km²，动态度1.22%，耕地与林地面积持续下降，预计分别下降173.6和91.8 km²，低于1990-2018年土地类型变化程度，总体上未来广州市土地变化趋于稳定。预测结果表明建设用地仍然会有一定程度的扩张，分析未来建设用地扩张对碳储量变化影响，揭示在自然发展条件下未来广州市碳储量分布变化。从广州市2018-2034年土地利用变化预测情况分析(图 4)，未来白云区和花都区的交界处建设用地的扩张能力仍较强，交界区域的城市扩张明显，中心城区与花都区和白云区具有较高可能性连成一个整体，建城区的斑块聚集度提高，散落在中心城区外围和周边区域的耕地在未来有较大的可能性会转为建设用地，聚集成更大的斑块。且扩张侵占的土地上，可能会进一步侵占更高碳密度的土地类型，如分布在从化区、白云区和增城区的林地等也有转为建设用地的趋势。2034年广州市碳储量仍会进一步降低，预计将会减少1.20×10⁶ t，降幅3.2%，年均下降7.50×10⁴ t，低于1990-2010年的年均下降1.09×10⁵ t，高于2010-2018年的年均下降6.00×10⁴ t，表明在自然发展条件下，广州市碳储量流失整体上已经没有1990-2010年剧烈。预测2034年从化区碳储量所占比例将上升到40%，广州市北部和中部森林资源丰富，北部受制于地形因素，城镇建设较为缓慢，未来仍是广州市乃至珠三角重要的碳汇区域，中部林地靠近市区，受到城市建设用地扩张风险较大，碳储量比重进一步降低。预测碳储量减少的区域在花都区和白云区的交界较为显著，未来这一区域的碳储量最易流失，增城区西部和黄埔区也会存在较大的碳储量流失风险，主要原因可能是该区域的林地有较大的概率转为建设用地。碳储量下降最高的区域分别是增城区、白云区和花都区，预计分别下降2.95×10⁵ t，2.10×10⁵ t和1.83×10⁵ t。黄埔区和白云区交界处的白云山是临近中心城区的较大森林覆盖区域，拥有丰富的林地和生物资源，对于保护中心城区的碳平衡和生态系统的稳定具有非常重要的价值。

此外其他区域也会发生不同程度的碳储量变化，从化区南部有一定碳流失风险，北部碳储量微量增加，南沙区整体较为稳定。未来部分区域碳储量流失仍然较为严峻。

2.2.2 建设用地与低密度碳储量重心迁移

城市建设用地扩张与碳储量生态系统服务是碳储量研究中的重要内容，城市建设用地扩张造成耕地林地等高密度固碳土地类型转成低密度固碳土地类型，加重碳储量生态系统功能的危机。在计算得到的碳储量分布图上划分碳密度等级，选择低密度碳储量空间分布，通过不同时期的重心迁移分析建设用地与碳储量的关系。图 5为广州市1990-2034年不同时期建设用地与低密度碳储量的重心变化，从重心迁移结果看，1990-2000年建设用地扩张的整体朝向东南部，2000-2010年重心向东北方向迁移，2010到2018年城市持续向东北扩张，但迁移距离减小，随着西部花都区等区域副中心城镇的扩张仍充满潜力，向西迁移趋势有所增加。总体上城市重心整体在向东北偏移，且逐渐由东北偏向正北，迁移距离先增后减。预测2018-2034年建设用地重心朝向西部，迁移距离在减少，说明城市扩张幅度在波动中趋于稳定，低密度碳储量空间重心迁移与建设用地扩张方向基本吻合，表明建设用地的扩张与低密度碳储量空间演变具有显著的相关性。

3 讨论与结论

本文在预测未来土地的驱动参数上选择有限，对人文因子仅考虑人口、GDP和交通分布，人文因子是一个复杂的指标，包括工业产值、工厂分布、政策制定、开发红线等，今后在驱动因子选择中可以构建科学的人文因子模型作为预测条件。未来城市土地变化中受到耕地红线、生态红线的限制，本次预测中暂不考虑政策因素，只考虑了自然情景下未来土地类型的预测，没有针对生态保护、耕地保护等限制条件下的情景模拟，以后的研究中可以考虑在受到“三生空间”政策影响下不同情景的未来城市生态系统功能变化。

(1) 1990-2018年广州市土地利用变化的特征表现为建设用地的快速扩张，增幅128.3%，变化量最高，其他用地均有不同程度变化，其中耕地面积减少量最高，其次是林地。1990-2005年是土地利用变化动态度较大的时期，2010年后土地利用变化动态度显著下降。

(2) 广州市碳储量分布具有空间差异规律，总体上表现为高值集聚在北部，低值集聚在西南。北部森林覆盖面积大，海拔高，是碳储量分布的主要位置，从化区、花都区和增城区是广州市重要的碳汇区域；西南部城市化程度高，城市扩张活动频繁，城市中心碳储量比较低。1990-2005年是碳流失较严重时期，与城市扩张密切相关，2010-2018年碳流失相对缓和。

(3) 预测未来建设用地面积仍会继续增加，扩张幅度将会下降，扩张重心朝向西部花都区和白云区等。预测未来碳储量会进一步下降，碳储量减少的区域主要位于广州市西北和主城区东部附近。西北区域建设用地未来可能与主城区形成大面积的低碳储量区域，对于城市生态系统稳定是较大的威胁，未来广州市部分区域仍会存在一定的碳流失问题。