陆地生态系统在吸收CO₂，调节气候变化方面起着关键作用^[1]，当区域碳储量提高时，能有效缓解温室效应对气候的影响^[2]，鉴于此，陆地生态系统固碳作用成为各国政府和学者共同关注的研究热点^[3]。随着国家“双碳”目标任务的提出，如何提升陆地生态系统碳汇能力成为关注的重要方面，尤其是对大尺度长时序生态系统碳储量时空演变格局研究。有关研究表明，在土地利用变化驱动下陆地生态格局与生态系统发生改变进而会影响碳循环过程^[4]，是碳储量变化的重要原因^[5]。不同土地类型具有不同的固碳能力，正确的土地管理可补偿碳损失^[6]。不合理土地利用是导致陆地生态系统碳源/汇失衡进而引起温室气体浓度增加主因之一。因此，深入研究土地利用变化下的生态系统碳储量动态演变规律，将对促进区域“双碳”目标实现具有重要的参考价值和借鉴意义。

在生态系统碳储量测算中，传统方法主要包括样地调查、遥感反演和模型模拟、IPCC清单法等，虽然精度高但是采样周期长和工作繁琐等不足，无法有效适用大尺度研究。为了弥补传统估算方法的不足，综合集成的模型模拟法得到发展，其中生态系统服务和权衡的综合评估模型(integrated valuation of ecosystem services and trade-offs, InVEST)因较为适用于长时序大尺度研究，且操作简便、评估精度较高等优点，近年得到应用。同时为了进一步模拟预测碳储量变化趋势，基于土地扩张分析策略(land expansion analysis strategy, LEAS)和多类型随机种子元胞自动机模型集成发展了斑块级土地利用模拟(patch-generating landuse simulation, PLUS)模型，其保留基于轮盘赌的自适应惯性竞争机制获取土地利用变化综合概率^[7]，并具有模拟精度高、数据处理快、能有效挖掘土地利用变化的原因和模拟多地类复杂演变的过程的优点^[8]，进而实现对碳储量模拟预测。因此，综合InVEST与PLUS模型用于模拟各种地类复杂演变及估算碳储量具有较好的应用前景^[9-10]。目前，流域尺度下长序列碳储量研究成为热点，多侧重于开展基于InVEST模型的流域碳储量时空变化特征分析研究^[11-12]，并利用PLUS模型等开展多情景模拟预测碳储量变化趋势，有关研究认为碳储量与土地利用空间分布一致，碳储量空间分异主要受气候因素影响^[13]，尤其碳储量与降水呈明显正相关，降水量的增加有利于碳损失的缓解^[14]。认为生态保护情景能有效缓解碳损失^[15]，相较于自然发展情景更有利于碳的蓄积^[16-17]，生态工程建设是流域碳储量增加的主导因素^[18]，林地等增加将对区域碳汇能力产生正向影响^[19]。以上有关研究证明耦合InVEST与PLUS模型开展流域生态系统碳储量时空变化及模拟预测是比较有效的方法。

密云水库是首都重要的地表饮用水源地、水资源战略储备基地。密云水库流域作为密云水库的载体，是区域重要的生态屏障与水源涵养区，在维系区域生态平衡中起着极其重要的作用。目前，密云水库流域生态系统的研究主要关注不同影响因子与生态环境系统之间的响应机制，而对生态系统碳储量状况鲜有报道。因此，本文选取密云水库流域作为研究区域，探究1985—2020年密云水库流域生态系统碳储量时空演变格局，并模拟预测2035年不同情境下碳储量变化趋势，以期为密云水库流域“双碳”目标下的山水林田湖草沙一体化保护修复和流域生态补偿标准优化提供科学依据。

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

密云水库流域位于东经115.44°—117.55°，北纬40.34°—41.61°之间，流域总面积15 396.3 km²，东西长约184 km，南北宽约130 km。流域地势整体呈现出西北高、东南低特征，地貌以山地为主，流域内西部为白河水系，东部为潮河水系，气候区属于半干旱与半湿润过渡带，具有典型大陆性季风气候特点，多年平均气温范围6~12 ℃，多年平均降水量约488 mm，降水主要集中在6~8月雨季。流域地跨河北省张家口市赤城县、沽源县，承德市兴隆县、滦平县、丰宁满族自治县以及北京市延庆区、怀柔区、密云区。

1.1.2 数据来源

(1) 土地利用数据。本文选取密云水库流域1985，1990，1995，2000，2005，2010，2015和2020年的土地利用类型数据，数据来源于中国科学院资源环境科学数据中心(http://www.resdc.cn/data)，空间分辨率30 m×30 m。按照《土地利用现状分类(GB/T 21010—2017)》对土地类型进行了部分整合和重分类，采用ArcMap 10.2重分类划分耕地、林地、草地、水域、建设用地和未利用土地6大类。

(2) 土地利用变化驱动因子。为了提升PLUS模型对未来土地利用空间格局的模拟精度，主要从影响土地利用变化的自然驱动、社会经济发展驱动及交通可达性因素3个主要方面出发，综合参考已有研究^{[9-10, 13]}选取自然、社会、经济方面16项驱动因子，包括1个限制转换因素、7个自然环境因素和8个社会经济因素，具体指标包括限制转换区域数据、自然环境数据(包括DEM、坡度、坡向、年均温度、年均降水量、土壤类型、到水域距离)、社会经济因素(包括人口密度、GDP值、距居民点距离、距铁路距离、距高速公路距离、距省道距离、距二级公路距离、距三级公路距离)。

对社会驱动因子(到铁路交通节点、公路交通线、水体、居民点的距离等)在进行数据预处理时，均通过ArcGIS 10.2的欧氏距离(Euclidean distance)分析工具计算获得。根据模型需要，将土地利用驱动因子栅格化后统一成与土地利用数据相同的投影坐标、行列号和空间分辨率。研究区土地利用变化驱动因子与数据来源详见表 1。

1.2 研究方法 1.2.1 不同情景设置

由于密云水库流域横跨京冀，不是完整的行政单元，其未来土地利用需求预测有一定难度。鉴于此，本文采用情景分析法，通过修改Markov模型的过程转移概率估算研究区2035年不同情景各用地类型的需求面积。具体情景包括: ①自然发展情景。利用PLUS模型Markov Chain模块，基于2010—2015年土地利用转移规律，根据研究区的实际情况，不设置各类用地之间的转换概率，以5 a为间隔设置条件下，预测2035年土地利用状况，并作为其余情景模拟的基础^[19]。②耕地保护情景。耕地是密云水库流域产业发展的重要载体。参考《北京市国土空间总体规划(2020—2035年)》《北京市土地资源保护利用规划》对耕地保护的目标要求，以《河北省国土空间规划(2021—2035年)》对承德、张家口耕地保护任务要求，本文设定耕地保护情景为减缓耕地向其他地类的转移速率并遏制建设用地的扩张速度。为此，设置耕地向林地、草地及水域的转换概率降低30%，向建设用地的转换概率降低80%，向未利用地的转换概率降低100%，以5 a为间隔预测2035年该情景下的土地利用情况。③生态保护情景。参照《京冀密云水库水源保护共同行动方案》《北京市密云水库上游地区空间保护规划(2021—2035年)》《张家口首都水源涵养功能区和生态环境支撑区建设规划(2019—2035年)》生态建设任务，加强对林地、草地和水域等生态用地的保护并设定其他用地向生态用地的转化。本情景遏制建设用地扩张，将林地、草地等生态用地向建设用地转换概率降低90%，向未利用地转换概率降低100%，假设开展退耕还林等措施，将耕地向林地的转换概率增加20%，将建设用地转换概率降低60%，以5 a为间隔预测2035年该情景下的土地利用情况。

1.2.2 基于PLUS模型的不同情景下的土地利用变化模拟

本研究中，初期数据用2010年土地利用数据，基准数据用2015年的，基于PLUS模型神经网络算法输入土地利用数据结合多种驱动力因子获取各地类的适宜性概率，得到2020年土地利用模拟结果，并与现实2020年土地利用情况进行匹配比较，表明耕地模拟精度达到86.59%，林地模拟精度达到96.59%，草地模拟精度达到90.15%，水域模拟精度达到73.39%，未利用地模拟精度为74.12%，建设用地模拟精度为85.33%。借助kappa系数验证结果显示，kappa值为0.877，总精度为0.937，说明在本流域尺度上模型模拟效果较好，可用于预测模拟2035年土地利用情况。

1.2.3 基于InVEST模型的碳储量估算

InVEST模型的碳储量模块把生态系统的碳储量划分为4个基本碳库^[20]：碳库可分为地上生物碳库、地下生物碳库、土壤有机质碳库和死亡有机质碳库^[21]。基于InVEST模型对不同土地利用的平均碳密度进行统计，然后用各土地利用的面积与其对应碳密度相乘并求和，得出研究区的总碳储量。

碳密度的计算公式为：

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

(1)

式中：C为地类的总碳密度(t/hm²); C_above为地类的地上生物的碳密度(t/hm²); C_below为地类的地下生物的碳密度(t/hm²); C_dead为地类的死亡有机质碳密度(t/hm²); C_soil为地类的土壤有机质碳密度(t/hm²)。

碳储量计算公式为：

$ C_{\text {total }}=C \cdot A $

(2)

式中：C_total为区域内地类的总碳储量(t); C为地类的总碳密度(t/hm²); A为地类的面积(hm²)。

碳密度数据参考李克让等^[22]、陈利军等^[23]、任玺锦等^[24]有关全国碳密度研究成果数据并结合密云水库流域实际进行修正处理。由于死亡有机质碳库数据来源缺乏且占碳库比例较小，本文不予考虑。由于地上、地下生物量碳密度与土壤碳密度与气温和降水密切相关^[25]，本文采用Alam等^[26]研究中的公式作为修正因子的公式。

碳密度修正公式为：

$ C_{\mathrm{SP}}=3.3968 \times P+3996.1\left(R^2=0.11\right) $

(3)

$ C_{\mathrm{BP}}=6.798 \times \mathrm{e}^{0.0054 P} \quad\left(R^2=0.70\right) $

(4)

$ C_{\mathrm{BT}}=28 \times T+398 \quad\left(R^2=0.47, p <0.01\right) $

(5)

式中：C_SP为根据年降水量得到的土壤有机质碳密度(kg/m²)；C_BP、C_BT分别为根据年降水量和年均气温得到的生物量碳密度(kg/m²); P为年降水量(mm); T为年均气温(℃)。

$ K_{\mathrm{BP}}=\frac{C_{\mathrm{BP}}^1}{C_{\mathrm{BP}}^2} $

(6)

$ K_{\mathrm{BT}}=\frac{C_{\mathrm{BT}}^1}{C_{\mathrm{BT}}^2} $

(7)

$ K_B=K_{\mathrm{BP}} \cdot K_{\mathrm{BT}} $

(8)

$ K_S=\frac{C_{\mathrm{SP}}^1}{C_{\mathrm{Sp}}^2} $

(9)

式中：K_BP为降水因子修正系数；K_BT为气温因子修正系数; K_B为生物量碳密度修正系数; K_S为土壤碳密度修正系数; C¹为研究区域碳密度(kg/m²); C²为全国碳密度(kg/m²)。

研究区域各地类碳密度详见表 2。

2 结果与分析 2.1 土地利用时空变化特征

由图 1可以看出，1985—2020年密云水库流域主要土地利用类型为林地、草地和耕地，面积总和约占全流域面积的99%;水域、建设用地、未利用地面积比例较小，仅占总面积的1%。1985—2020年，林地、建设用地、水域分别增加1 873.82, 124.89和20.32 km²，其中，林地比例由52.80%上升至64.97%，建设用地比例由0.44%上升至1.25%，水域比例由0.76%上升至0.90%。耕地、草地、未利用地分别减少970.82, 1 047.91和0.32 km²，其中，耕地比例由16.90%降至10.60%，草地比例由29.09%降至22.28%，未利用地比例由0.005%降至0.003%。

基于8期土地利用数据，利用PLUS模型Markov chain模块得到土地利用转移矩阵(表 3)。1985—2020年，密云水库流域21.84%的土地发生转移，主要发生在草地、耕地、林地和建设用地间(图 2)。其中，转出面积最大的土地利用类型为草地和耕地，分别占总流转面积的54.64%和35.66%，而转入面积最大的土地利用为林地，占总流转面积的64.68%。草地有89.97%转向林地，耕地有46.85%的转向草地，转移主要出现在流域西部赤城县、西北部沽源县及中北部丰宁县等地区; 耕地有42.93%转向林地，耕地有8.61%的转向建设用地，分散分布于流域各个城乡建成区周边。究其原因，一方面为保护密云水库流域实施退耕还林、京津风沙源治理、山水工程等生态建设，促进林地增加，耕地减少; 另一方面社会经济活动加强，城镇化进程加快，对土地开发利用需求增大，建设用地的增加，耕地、草地减少。

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

运用InVEST模型中的Carbon模块估算密云水库流域总碳储量。1985，1990，1995，2000，2005，2010，2015和2020年密云水库流域总碳储量分别为4.94×10⁸ t, 4.96×10⁸ t, 5.03×10⁸ t, 5.05×10⁸ t, 5.10×10⁸ t, 5.11×10⁸ t, 5.15×10⁸ t, 5.18×10⁸ t。1985—2020年，总碳储量持续增加，35 a间碳储量增加2.35×10⁷ t。从不同土地类型碳储量来看，林地碳储量最高，其次为草地、耕地，水域、建设用地，未利用地碳储量最少。1985—2020年，各土地利用类型中对碳储量的贡献程度最大的为林地、耕地和草地，共占总碳储量的98%(图 3)。

研究期间，林地、建设用地碳储量持续上升，耕地持续下降，草地碳储量波动下降，未利用地和水域碳储量波动上升。从碳储量空间分布来看，碳储量高值主要分布在延庆、怀柔、密云、滦平、兴隆等林地稠密地区，沽源、赤城、丰宁碳储量相对偏低(图 4)。

2.3 土地利用变化模拟与碳储量预测 2.3.1 驱动因子分析

通过PLUS模型的用地扩张分析模块得出主要用地类型的驱动因子贡献度(图 5)。耕地变化驱动力主要受高程、坡度、人口密度影响，耕地主要增加区域位于高程较低、坡度较平缓、人口密度相对较少的地区; 林地驱动力主要受坡度、年均降水量、高程、距居民点距离等因素影响，流域内中部、南部、东部地区有广袤地区适宜林地建设，这些地区是林地增加潜在区域; 草地变化驱动力主要受年均降水量、坡度、高程、人口密度、年均温度等因素影响，流域西北部适宜的降水量，较高高程等促进草地面积增加; 水域变化驱动力主要受高程、距水域距离、人口密度、坡度等因素影响，距现有水源越近、适宜地理环境等条件促进水域面积增加; 未利用地变化驱动因素主要是年均降水量、距三级公路距离; 建设用地变化驱动因素主要是人口密度、高程、坡度等，尤其是人口密度较高地区，随着城镇化建设，建设用地表现出较快的增长趋势。

2.3.2 密云水库流域2035年多情景模式下土地利用预测

从表 4可以看出，与2020年各土地利用面积相比，3种情景下2035年林地、建设用地均呈增长趋势，而耕地、草地、水域、未利用地均呈减少趋势，主要的土地利用转移为耕地、草地转向林地，耕地转向建设用地(图 6)。其中在生态保护情景下，与2020年相比林地增加最多，达到911.01 km²，建设用地扩张在3种情景下最小，说明通过生态保护限制了其他土地利用类型向建设用地的转换，生态保护下更有利于流域生态环境的持续改善。

2.3.3 2035年碳储量预测

从表 5看出，与2020年相比，3种情景下研究区碳储量总量均增加。在自然发展情景下，碳储量增加1.08×10⁷ t。在耕地保护情景下，碳储量增加1.24×10⁷ t。在生态保护情景下，生态保护修复力度加大，林地面积持续增长，进而带来碳储量的增长，期间增长1.31×10⁷ t。在3种情景下，林地碳储量贡献率最大，是所有土地利用类型中最主要的碳储量来源(图 7)。总体上，采取生态保护措施能够取得较高的碳储量水平，其他土地利用类型向林地转化的概率增加，进一步增强森林碳汇功能，流域整体固碳能力得到巩固和提升。

3 讨论

(1) 在土地利用变化对碳储量的响应方面，土地利用变化是影响陆地生态系统中碳储存变化的关键因素，无论是草地到林地、耕地到林地的转移，还是建设用地的开发扩张，这些变化对陆地生态系统碳储量时空格局变化产生深远的影响。从研究结果来看，各用地类型面积的增减将会影响碳储量的增减，这与有关研究结果相似，碳储量的空间分布与区域土地利用空间布局有显著一致性^[27]，说明土地利用的空间布局和结构变化在很大程度上影响着区域碳储量，可通过生态建设增加林草地、控制林地转为其他用地及合理控制建设用地扩张来实现土地结构的优化，进而促进生态系统碳储量巩固与提升。

(2) 在模型及碳密度选取准确性与局限性方面，InVEST模型对碳储量测算是基于大尺度土地利用变化，对土地利用内部结构和植被年龄的固碳能力差异考虑不足^[28]。InVEST模型中的Carbon模块对碳循环过程进行了简化，用静态碳密度来进行估算，而且研究中土地利用类型碳密度综合参考已有全国碳密度研究成果，并根据区域实际情况进行修正处理，但仍与研究区真实的碳密度数据存在一定的误差，影响碳储量测算的准确性，下一步结合样地调查数据对碳密度数据进行修正，以提高测算精度。

(3) 2020年林地面积占流域总面积的64.97%，相关研究表明，不同的树种对碳储量影响也存在差异，下一步需要对不同树种碳汇能力及分布状况进行研究，以期在增加林地和控制林地转为其他地类的同时注重树种结构的优化。本研究重点探讨长时序土地利用变化引起碳储量时空演替特征，反过来碳储量对土地利用变化脆弱性响应的研究还有待开展。碳储存服务作为生态系统服务的重要组成之一^[29]，研究其变化及脆弱性对缓解生态系统服务损失及退化有着重要意义^[30]，下一步将深入探讨碳储量脆弱性等问题。

4 结论

(1) 密云水库流域主要土地利用类型为林地、草地和耕地，面积总和约占全流域总面积的99%。1985—2020年，密云水库流域21.84%的土地发生转移，其中，转出面积最大的土地利用类型为草地和耕地，而转入面积最大的土地利用为林地，林地持续增加，林地占全流域面积比由52.80%上升至64.97%。

(2) 1985—2020年，总碳储量持续增加，由1985年4.94×10⁸ t增加到2020年5.18×10⁸ t, 35 a间碳储量增加2.35×10⁷ t。从不同土地类型碳储量来看，林地碳储量最高，其次为草地、耕地，水域、建设用地，未利用地碳储量最少。期间，林地、建设用地碳储量持续上升，耕地持续下降，草地碳储量波动下降，未利用地和水域碳储量波动上升。

(3) 流域碳储量呈现出“中南部高，西部、东北低”的空间格局，中南部对应流域山区地貌，为水源涵养核心区域，森林覆盖率较高，而西部、东北部主要分布草地、耕地及城镇建设用地等，相对碳储量偏低。

(4) 预测2035年生态保护情景更有利于碳储量储存。总体上，采取生态保护措施能够取得较高的碳储量水平，其他土地利用类型向林地转化的概率增加，流域整体固碳能力得到巩固和提升。