气候变暖问题是当前人类社会面临的最严峻的挑战之一, 深刻影响着全球的生态环境与人类的生存发展^[1]。以CO₂为代表的温室气体的排放已被公认为是导致全球气候变化的主要原因^[2—3]。在过去几十年里, 飞速的经济增长与高强度的土地开发导致了中国能源消费与碳排放的持续增长, 中国目前已经成为世界上最大的碳排放国, 年碳排放量约占全球总量的三分之一^[4—5]。为了应对日益严峻的资源与环境压力, 积极参与并推进全球气候治理工作, 中国政府于2020年9月提出将采取更加有力的政策与措施, 努力争取碳排放在2030年前达到峰值, 并在2060年之前实现“碳中和”目标^[6—7]。

生态系统服务是指通过生态系统的结构、过程和功能直接或间接得到的生命支持产品和服务, 其维持与改善是保障区域生态安全, 提升人类福祉的基础^[8]。然而, 由于自然生态环境在全球范围内的恶化, 全球生物多样性与生态系统服务水平正呈现迅速下降的趋势^[9]。过去的研究表明, 土地利用变化是引起区域生态系统服务功能和质量变化的重要因素^[10]。随着未来社会经济政策的转型与“双碳”目标的逐步实现, 土地利用变化的趋势也将发生转变。同时, 由于短期内我国以化石能源为主体的能源供给结构难以发生改变^[11], 低碳减排政策的实施不可避免地会造成一定的经济损失。在此背景下, 如何基于有限的土地资源, 合理协调低碳减排与经济效益之间的关系, 优化土地利用结构与空间布局, 并在此基础上模拟生态系统服务的未来演化趋势并有效规避生态系统服务退化风险, 是全球气候变化背景下推动经济社会低碳转型与生态系统服务价值持续提升的重大战略需求。

土地利用情景模拟是实现土地利用优化配置的重要方法^[12]。近年来, 研究者们从结构与空间布局两个方面开展了大量的土地利用模拟研究, 并开发出一系列的预测模型。结构预测模型主要包括马尔科夫(Markov)模型^[13]、系统动力学(SD)模型^[14]、多目标规划(MOP)模型^[15]等。其中, MOP模型不仅具有SD模型的正反馈特征, 而且突破了目标函数一维性的限制, 可以在兼顾多种预设目标的基础上进行土地利用优化配置, 已广泛应用于土地利用模拟研究之中^{[10, 15]}。然而, 过去的结构预测研究往往都是通过协调研究区的生态与经济效益, 实现土地利用结构的优化, 缺乏以低碳减排作为目标, 并与经济效益进行权衡的研究。空间布局预测模型中应用较为广泛的有CA模型^[16]、CLUE-S模型^[17]、FLUS模型^[18]等。然而, 这些模型在揭示土地利用变化的驱动因素及其影响方面有所欠缺, 难以灵活地模拟多类土地利用类型的斑块级变化。Xun L^[19]等开发的斑块生成土地利用模拟模型(Patch-generating Land Use Simulation Model, PLUS)既保留了FLUS模型中自适应惯性竞争和轮盘竞争机制的优点, 又基于两期土地利用数据间各类用地的相互转化样本进行训练, 以此得出各土地利用类型的发展潜力, 能够更加准确地模拟土地利用空间分布的动态变化。

北京市是中国的政治、文化、国际交往与科技创新中心, 在低碳发展的背景下, 对其未来土地利用变化进行多情景模拟预测, 不仅能稳步推进北京市的碳中和行动, 持续提高生态系统服务水平, 也能对其他城市起到一定的引领与示范作用。本文通过耦合MOP与PLUS模型, 在权衡经济效益与低碳减排的基础上, 构建自然演变、经济优先以及低碳发展3种情景, 分别对不同情景下北京市2030年的土地利用结构与空间布局进行预测。在此基础上, 通过量化不同情景下的生态系统服务价值, 进一步对比分析不同情景下生态系统服务的差异, 并据此对北京市低碳发展背景下的具有生态系统服务提升潜力的区域进行识别。研究结果可以为北京市未来的土地资源的配置以及生态环境的改善提供科学的决策依据。

1 研究区域概况与数据来源 1.1 研究区域概况

北京市地处华北平原的北部, 东面与天津市毗连, 其余均与河北省相邻(图 1), 总面积约1.64万km²。近年来, 北京市高度重视应对气候变化工作, 持续推动产业结构优化和能源清洁转型、大力疏解非首都功能, 燃煤量大幅下降, 为全市碳达峰工作打下了扎实的基础。“十四五”时期, 北京将继续实行碳排放总量和强度“双控”机制, 并进一步开展碳减排专项行动, 实现碳排放稳中有降。

1.2 数据来源

本研究所使用的2期土地利用数据(2010年、2020年)来源于全球地表覆盖数据GlobeLand 30(http://www.globallandcover.com), 空间分辨率为30m, 为了便于后续的模拟分析, 参考前人研究^[20—21]将土地利用划分为6种类型：耕地、林地、草地、水域、建设用地、未利用地; DEM高程数据来源于地理空间数据云平台(http://www.gscloud.cn); 气象数据(包括年均温度、年均降水量、干燥度等)、土壤数据(土壤类型与侵蚀强度)以及社会经济数据(GDP与人口密度)来源于中国科学院资源环境科学数据中心(http://www.resdc.cn), 空间分辨率为1km; 夜间灯光数据来源于VIIRS_DNB_VNLV2数据集(https://eogdata.mines.edu), 空间分辨率为500m;现状道路、铁路、居民地以及河流水域矢量数据均来源于全国基础地理数据库(http://www.webmap.cn)。以上数据在Arcgis中进行投影变换、裁剪、重采样等一系列预处理工作, 并将栅格数据的分辨率统一转化为30m×30m。

2 研究方法 2.1 基于MOP模型的多情景土地利用结构预测

MOP即多目标规划模型, 是土地利用数量结构优化研究中常用的模型之一。该模型能够在一系列客观或主观约束条件下, 输出一个或多个目标达到最值时的决策^[22]。

2.1.1 多目标函数的构建

(1) 经济效益目标

经济效益目标函数公式如下：

(1)

式中：F₁为经济效益; x_i为各类土地利用类型(x_1—6分别代表耕地、林地、草地、水域、建设用地以及未利用地)的面积; E_i为各类土地利用类型的经济效益系数, 即单位面积的土地利用类型i所产生的经济效益。

参考前人的研究^{[10, 15]}, 通过对2010—2020年研究区内各土地利用类型对应的产值与其面积之比, 得到逐年的经济效益系数, 再借助GM(1, 1)灰色预测模型预测2030年研究区的经济效益系数E_i(表 1)。

(2) 低碳减排目标

低碳减排目标函数公式如下：

(2)

式中：F₂为碳排放量; x_i为各类土地利用类型的面积; C_i为各类土地利用类型的碳排放系数, 即单位面积的土地利用类型i所产生的碳排放量。

在碳排放量评估研究中, 常使用化石能源消费量间接评估建设用地的碳排放量^[23—24]。煤、石油、天然气这3类化石能源消费数据可以从《北京统计年鉴》中获取, 参照孙伟等^[25]的计算方法, 分别评估3类化石能源的燃烧释碳量。此外, 考虑到北京市人口众多, 本研究将人类呼吸产生的碳排放也纳入到评估范围内^[25]。化石能源与人类自身呼吸释碳量之和与建设用地面积的比值即为建设用地的碳排放系数。

本研究中, 耕地、林地、草地、水域与未利用地的碳排放系数均参考前人研究中对于北京市不同土地利用类型的碳排放/吸收能力评估结果进行设置^[26—29]。各土地利用类型的碳排放系数C_i详见表 1。

2.1.2 发展情景的设定及约束条件

本文共设定了3种未来发展情景, 分别为自然演变情景、经济优先情景以及低碳发展情景。

自然演变情景指研究区内各土地利用类型均按照现有的趋势进行自然演变, 本文借助马尔科夫(Markov)模型实现自然演变情景下各土地利用类型的面积预测; 经济优先情景是指在《北京城市总体规划(2016年—2035年)》(以下简称《规划》)的约束下, 以经济效益的最大化为优化目标的情景; 低碳发展情景是指在《规划》的约束下, 综合考虑经济效益与低碳减排两个优化目标的情景。由于两个目标函数的优化方向不同, 数量级与数值单位也有所差异。因此参考前人研究, 利用自然发展以自然演变情景下的函数解为参照对两个目标函数进行无量纲处理后再进行叠加计算, 并选择出综合效益最高的优化方案^[30]。其计算公式如下：

(3)

式中：F_C为低碳发展情景下的综合效益; F₁与F₂分别为该情景下的经济效益与碳排放量; F_N1与F_N2分别为自然发展情景下的经济效益与碳排放量。

为了保证多情境下的未来土地利用变化符合自然发展规律和政府的规划预期, 本文根据2030年自然发展情景下的土地利用需求数据以及《规划》来构建多目标优化规划的约束条件(表 2)。

2.2 基于PLUS模型的土地利用空间分布模拟

PLUS模型以2期土地利用变化数据为基础, 利用随机森林算法以及基于多类型随机种子机制的元胞自动机(CA)模型, 可以更好地挖掘不同土地利用类型扩张的驱动因素并精确地模拟未来多类土地利用的斑块级变化。Xun L.等^[19]的研究表明, PLUS模型的模拟精度高于目前相关研究中常用的其他土地利用模拟模型。

本研究以2010年与2020年2期土地利用数据为基础, 利用高程、坡度、气温、降水、干燥度、土壤类型、土壤侵蚀类型、人口密度、GDP、夜间灯光指数、距高速公路距离、距主要道路(国道、省道、环路等)距离、距其他道路距离、距铁路距离、距居民地距离以及距河流水域距离共计16个驱动因子, 评估北京市2010—2020年间的土地利用变化趋势与驱动因素。并以2030年土地利用结构优化结果为数量目标, 预测2030年的土地利用空间分布格局。此外, 为保证模拟精度达到研究需求, 在进行正式模拟前, 需要以2010年土地利用数据作为模拟基准图, 利用PLUS模型得到2020年土地利用模拟结果, 并将其与2020年真实土地利用数据进行对比验证, 评价模型的模拟精度。

2.3 生态系统服务价值评估

本文以250m*250m的的格网作为分析单元, 利用当量因子法对研究区生态系统服务价值进行量化。该方法最早由Costanza^[31]等提出, 谢高地等对该方法进行了改进并得出中国的生态系统服务价值当量^[32]。本研究以谢高地等计算的各土地利用类型生态系统服务价值当量为基础, 根据北京市的实际情况, 结合前人研究中的计算方法^[33—35], 最终确定研究区1个生态系统服务价值当量因子经济价值约为2375.05元/hm², 研究区生态系统服务价值系数如表 3所示。生态系统服务价值的计算公式如下：

(4)

式中, ESV为生态系统服务价值; i为土地利用类型; A_i为第i类土地利用类型的面积; VC_i为第i类土地利用类型单位面积的生态系统服务价值(即表 3所示的服务价值系数)。

3 研究结果 3.1 不同情景下土地利用结构预测

不同情景下的土地利用结构预测结果如表 4所示。自然发展情景下, 耕地面积显著减少, 建设用地面积显著增加。相较于2020年, 耕地面积减少了968.17km², 建设设用地面积增加了857.07km²。研究区林地、草地与水域均呈现缓慢增加的趋势。在经济优先情景下, 由于受到耕地保护、控制建设用地规模等规划政策的限制, 耕地面积减少的速度相较于自然演变情景有所减缓, 建设用地扩张的趋势也得到了一定程度的遏制。与2020年相比, 林地面积增加了315.37km², 草地与水域面积保持不变。在低碳发展情景下, 耕地面积的收缩得到了进一步的控制, 建设用地面积呈现减少趋势。林地面积增加了739.86km²。, 是3种情景中增长量最大的, 草地和水域面积与2020年相同。此外, 未利用地在研究区内面积占比较少, 且3种情景下差异不明显。

3.2 基于PLUS模型的土地利用空间分布模拟

本研究中, PLUS模型模拟的整体精度为84.13%, Kappa系数为0.77, Kappa值位于0.61-0.80区间内, 说明模拟结果与实际数据呈现高度一致性。模型模拟精度较高, 能够较好地模拟研究区土地利用动态变化。

将3种情景下的土地利用结构预测结果分别带入PLUS模型, 以2020年土地利用数据为基础, 对不同情景下研究区2030年土地利用空间分布进行模拟(图 2)。结果显示, 自然演变情景下, 北京东南部平原区大面积的耕地被建设用地侵占, 西北部山区的林地与草地面积得到了一定程度的增长, 密云水库与西侧官厅水库附近的水域面积得到了一定程度的扩张; 在经济优先情景下, 东南部平原区建设用地扩张的趋势得到了有效的遏制, 西北部山区内有一部分耕地转化成为林地与草地, 而一些原本的草地则进一步转变为林地; 在低碳发展情景下, 平原区的建设用地规模与分布与2020年基本一致, 西北部山区中的耕地与建设用地面积减少, 部分草地又转化升级为林地, 林地的面积显著增加, 整体的土地利用空间布局更为合理。

3.3 生态系统服务价值评估结果

利用土地利用空间分布模拟结果进行生态系统服务评估(表 5, 图 3), 在此基础上, 将3种未来情景下的生态系统服务价值与2020年的数据进行对比分析, 结果如图 4所示。

在自然演变情景下, 研究区生态系统服务价值从2020年的554.81亿元提升至567.68亿元, 共提升了12.87亿元。生态系统服务价值降低的区域主要位于东南部平原区, 而西北部山区的生态系统服务价值有一定程度的提升。

在经济优先情景下, 研究区2030年生态系统服务价值约为563.77亿元, 相较于2020年提升了8.96亿元, 服务价值总量略低于自然演变情景。生态系统服务价值降低的区域依旧主要分布在东南部平原区, 但规模得到了一定程度控制。同时, 得益于林地面积的提升, 西北部生态系统服务价值提升区域的面积显著增多。

在低碳发展情景下, 研究区生态系统服务价值约为590亿元, 与2020年相比提升了35.19亿元, 显著高于另外两种情景。东南部平原区的生态系统服务价值并没有呈现下降的趋势, 反而一些城市区域的服务价值有所提升。同时, 西北部山区的服务价值提升区域的面积也明显大于另外两种情景, 研究区整体的生态质量得到了显著提升。

3.4 低碳发展背景下生态系统服务提升空间的识别与分析

在本研究中, 生态系统服务价值提升的区域代表其具有生态质量改善的潜力。因此, 本文将低碳发展情景下服务价值增加的区域作为未来应重点建设的生态系统服务提升空间, 并按照服务价值提升的程度划定生态提升潜力(图 5), 为北京市规划政策的制定与生态工程的开展提供参考。

研究结果表明, 提升潜力较大的区域面积占比约为7.98%, 主要分布于研究区西北部的延庆区以及东北部的平谷区与密云区。提升潜力中等的区域面积占比约为26.57%, 主要分布于北部山区、门头沟区东部以及第二道绿色隔离地区(位于北京市五环路到六环路附近)。提升潜力较低的区域面积占比约为65.45%, 且在北部山区的分布较为密集。

4 讨论与结论 4.1 讨论

随着气候变化问题的日益严峻以及我国“双碳”目标的逐步落实, 如何合理协调经济效益与低碳减排之间的关系, 对有限的土地资源进行合理地开发利用, 并逐步提升生态系统服务价值, 将成为我国未来发展中的重要问题之一^[36]。刘慧灵等^[37]利用多目标规划模型, 从低碳经济的视角对福州市的土地利用结构进行优化。何海珊等^[38]以深圳市为例, 通过构建碳汇最大化与碳排放最小化两种低碳发展情景, 借助FLUS模型实现了低碳导向下的土地利用演变模拟。从长远来看, 通过土地利用多情景模拟研究, 能够科学准确地评估城市未来发展中不同区域的生态提升潜力, 并辅助决策者制定出更具有前瞻性与可操作性的规划与政策。然而, 过去的研究却缺乏与实际规划工作相互衔接的尝试。本研究通过耦合MOP与PLUS模型, 不仅提出了一种低碳发展背景下土地利用结构优化与空间布局模拟的新方法, 并且在此基础上通过评估生态系统服务价值, 进一步对北京市未来生态系统服务提升区域进行了识别, 将传统的多情景土地利用模拟研究进行了延伸与扩展, 使得其与实际规划工作可以更好地结合。

通过横向对比本研究中土地利用结构与空间分布模拟结果可以发现：东南部平原区的建设用地容易向耕地扩张, 从而导致耕地的退化; 西北部山区的部分耕地以及建设用地能够转化为林地与草地, 这与前人的研究结果类似^[39]; 在3种情景中, 只有自然发展情景实现了水域面积的增长, 这说明对于研究区而言, 盲目提升水域面积对于经济效益与低碳减排效益的提升效率并不高^[40]。通过对比各情景下土地利用模拟与生态系统服务价值评价结果不难发现：西北部林地面积的增加以及对耕地的有效控制是低碳发展情景下生态系统服务价值高于其他情景的主要原因; 自然发展情景与经济优先情景虽然都相较于2020年服务价值有所提升, 且总量近似, 但各类用地的贡献却有所差异, 前者主要得益于水域面积的增长, 而后者主要是由于林地面积的增加。

总体而言, 本研究可从低碳发展的视角为北京市未来发展提供如下建议：(1)坚守耕地和永久基本农田保护红线, 控制建设用地的规模, 尤其需要警惕东南部平原区建设用地对耕地的侵占。(2)林地面积的增加不仅能够降低净碳排放量, 还能显著提升北京市的生态系统服务价值, 政府可以通过人工造林、退耕还林等方式增加林地面积, 并逐步提升森林生态系统的质量。(3)延庆区、平谷区、密云区、门头沟区东部以及第二道绿色隔离地区的生态系统服务提升潜力较大, 未来可有针对性地实施相应的生态工程逐渐提升其生态质量。

值得注意的是, 本文参考前人研究, 使用土地利用类型来对经济效益与碳排放量进行量化, 与真实值存在一定的差距, 但具有相对意义^[38]。由于本文主要的目标是探讨两个目标权衡下的土地利用优化配置, 因此利用该方法进行量化可以满足研究的需求。

4.2 结论

本文不仅从低碳发展的视角多情景地模拟了北京市2030年的土地利用空间布局与生态系统服务价值。还进一步针对低碳发展情景, 划定了生态系统服务提升区域。研究结果表明：

(1) 本文提出的耦合MOP与PLUS模型的研究方法能够较好地权衡低碳减排与经济发展之间的关系, 对土地利用结构与空间布局进行优化与模拟。并且该方法在北京市的适用性较强, 土地利用模拟的整体精度为84.13%, Kappa系数为0.77, 为今后该区域的土地利用模拟提供了方法参考。

(2) 3种情景下的土地利用结构与空间布局存在差异。在自然演变情景下, 东南部平原区的大面积耕地被建设用地侵占, 林地、草地与水域面积有小幅度的增长; 在经济优先情景下, 平原区建设用地向耕地扩张的趋势得到了一定程度的遏制, 西北部山区中的林地面积也有所增长; 在低碳发展情景下, 城市建设用地规模有小幅度下降, 西北部山区中的耕地与建设用地逐渐向林地与草地转化, 林地面积显著增加。

(3) 3种情景下生态系统服务水平有所不同, 从总量上来看, 低碳发展情景最高(590亿元), 自然演变情景次之(563.77亿元), 经济优先情景最低(567.68亿元)。从空间分布上看, 低碳发展情景下服务价值提升区域的面积显著高于另外两种情景, 且服务价值下降的区域面积也显著低于另外两种情景。

(4) 通过对低碳发展情景下生态系统服务价值提升区域的识别, 为北京市规划政策的制定与生态工程的开展提供参考。其中, 提升潜力较大的区域面积占比约7.98%、潜力中等区域的面积占比约26.57%, 潜力较低区域的面积占比约65.45%。