北京市西北部生态涵养区未来土地利用及生态系统服务变化情景模拟

陈新云，王 甜，李宝健

(1.国家林业和草原局 调查规划设计院，北京 100714；2.常州工学院 计算机信息工程学院，江苏 常州 213032；3.中国石油天然气股份有限公司 管道分公司，河北 廊坊 065000)

生态系统服务是指人类从生态系统中获得的各种效益，包括有助于提升人类福祉或维持地球生命的产品和服务[1-2]。过去几十年里，因人口增长和城市化进程加快带来的如粮食、水、燃料等自然资源的短缺，迫使人类以前所未有的强度和速度改变着地球生态系统[1,3-4]。生态系统评估报告认为，全球超过60%的生态系统服务已经或正在发生退化，进而影响了当前和未来生态系统服务的供给[1]。人类活动导致的土地利用变化是影响生态系统服务供给的决定因素之一，其中，某些特定的土地利用类型与生态系统服务密切相关，如生态系统的气候调节和水土保持服务主要由森林提供，因此，探讨土地利用变化与生态系统服务之间的关系一直是各国科学家和决策者们关注的重点[5-7]。

近年来，土地利用变化对生态系统服务供给的影响研究已取得一定进展[8-10]。有研究表明，耕地的集约化会导致碳储量的减少和土壤侵蚀风险的增加，而以经济效益优先作为优化目标的城市化以及建设用地的增加会破坏地表水平衡并影响区域气候[11-12]。也有研究发现，土地利用变化对生态系统服务的影响因社会经济背景和时空尺度而存在较大差异，这主要是由于人类需求的多样化和时空异质性导致生态系统服务间出现更为复杂的相互作用[13]。此外，土地利用变化与生态系统服务之间的关系证明了生态系统服务评估在指导土地利用规划和生态系统管理方面的重要性[14]。具体而言，生态系统服务评估可用于评价、比较和选择不同土地利用情景下的土地利用规划方案，例如，D.A.Shoemakeretal[15]利用多层次的城市增长模型模拟了2030年不同情景下的土地利用，并与生态系统服务评估模型相结合探讨了城市格局与生态系统服务之间的权衡。Y.Wuetal[11]提出了自然演变、农田保护和生态恢复3种情景，并探讨了6种生态系统服务对土地利用时空变化的响应特征。尽管土地利用时空格局的模拟能有效地预测未来土地利用的发展轨迹并支持土地利用决策的制定，但多数研究只是模拟每个地类的演变过程，不同的土地利用/土地覆盖变化过程在多数情况下会相互作用[10]。因此，本研究提出了集成马尔可夫链(Markov chain)与FLUS模型的“交互耦合”模型，该模型不仅可以解决不同土地利用类型之间的相互竞争作用，也能使土地利用需求预测与局部变化分配相结合，有效地处理自然过程和人类活动影响下各种土地利用类型相互转化的不确定性和复杂性[16-17]。

北京市西北部生态涵养区具有丰富的景观结构、生境类型和自然资源，其土地利用的变化必然对首都的生态系统服务和居民福祉有着潜在的重要影响。因此，本研究以北京市西北部生态涵养区为例，结合土地利用预测模型FLUS和马尔可夫链(Markov chain)以及生态系统服务评估模型InVEST，探讨2015-2030年不同发展情景下的土地利用变化如何影响生态系统固碳和产水量服务。内容包括：1)基于FLUS模型和马尔可夫链预测2030年自然演变、生态控制和城市发展情景下的土地利用时空动态；2)评估2015-2030年3种发展情景下的固碳和产水量服务的时空变化特征；3)揭示土地利用变化对生态系统服务的潜在影响。

1 材料与方法

1.1 研究区概况

北京市西北部生态涵养区(39°31′-41°04′N，115°24′-117°29′E)涵盖北部和西部的山区，总面积约11 247 km2(图1)。最高海拔2 291 m，最低海拔11 m。该区属典型的温带季风气候区，年平均降水量达576.71 mm，主要发生在6-8月，年平均温度2.5℃～13.4℃。本区包括北京最重要的山区、水源涵养地和基本农田等核心的生态用地，生态系统类型多样，生物多样性丰富。自20世纪末以来，该地区经历了快速土地利用转型和土地利用政策，导致了生态系统服务的退化，如区域碳储存损失和供水和净化能力下降[18-19]。近年来，北京市政府将生态涵养区列为北京市可持续发展的重点保护区域。

图1 北京市西北部生态涵养区位置Fig.1 The map of the ecological conservation area in Beijing

1.2 数据来源及处理

本研究用到的数据主要包括：高程数据(DEM)来源于全国地理信息资源目录服务系统(http://www.webmap.cn/)。气象数据，来源于国家地球系统科学数据共享服务平台，主要有年降水量、月降水量、温度和日照时数等(http://www.geodata.cn)。土壤属性、根系限制层深度和200 cm深度植物枯萎点的有效土壤含水量数据，来源于世界土壤数据库(Harmonized World Soil Database，http://webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database/)。植被蒸散发数据，主要来源于联合国粮农组织FAO的在线目录(http://www.fao.org/3/X0490E/x0490e0b.htm)。土地利用数据是2000年和2015年生态涵养区土地利用数据，来源于中国科学院资源环境科学数据中心(http://www.resdc.cn/)，采用欧阳志云等[20]提出的生态系统分类体系，对数据进行整合和重分类，最终将土地利用类型分为7类：草地、水域、农田、建设用地、未利用地、林地和灌木林地。利用Arcgis10.3软件将所有空间数据的分辨率重采样为30 m，投影坐标系统一为(Beijing_1954_3_Degree_GK_CM_ 114E)。

1.3 情景设定与土地利用预测

1.3.1 情景设定 基于对2000-2015年土地利用变化的驱动因素和北京市政府实施的土地利用政策的分析，提出了3种土地利用情景：自然演变情景、生态控制情景和城市快速发展情景，并将2015年的土地利用数据作为基础数据，在其基础上进行土地利用的模拟。

1.3.1.1 自然演变情景 主要是基于过去研究区土地利用的转变速率持续保持不变的假设而构建的，其人口、经济和技术革新的发展趋势将继续与目前的状况保持一致，因此，土地利用变化率被认为与2000-2015年的年度变化一致，选择马尔可夫模型来模拟土地利用需求量。此外，本情景假设气候变量的年均温和降水量与当前保持一致。

1.3.1.2 生态控制情景 主要是保障研究区的生态功能，保护环境质量安全，以生态效益优先作为优化目标，强化生态控制区的保护。根据《北京市总体规划(2016-2035年)》的要求，主要以保护密云水库、自然保护区、基本农田、公园绿地等生态区域为主，控制建设用地规模，减少经济发展对环境的影响。

1.3.1.3 城市快速发展情景 将充分发挥经济产出潜力较高的土地利用效益，以研究区的经济效益最大化为目标。在这种背景下，人口、技术和经济将快速增长，交通水利建设用地面积迅速扩大。此外，由于人口对粮食的需求不断增长，农田面积的下降幅度较其他情景小。

1.3.2 土地利用预测模型 首先利用FLUS模型的神经网络算法(ANN)获取各类用地的适宜性概率，通过元胞自动机(CA)模型提高模型的适用性，其次通过马尔可夫链(Markov chain)获取2030年各用地类型需求量，并设置邻域权重参数、成本矩阵及限制发展区等相关参数完成对2030年土地利用情景的模拟。从现有文献中选取土地利用变化的驱动因素作为预测变量进行检验(图2)[21-23]，主要计算模块如下：

1.3.2.1 基于神经网络的适宜性概率计算 神经网络算法(ANN)包括预测和训练阶段，由输入层、隐含层和输出层组成，其计算公式为:

(1)

式中，sp(p,k,t)为t时间和p栅格下第k种土地利用类型适宜性概率；wjk为输出层与隐含层之间的权值；Sigmoid()是从隐含层到输出层的激励函数；netj表示第j个隐藏栅格p在t时刻接收到的信号，神经网络算法输出的各土地利用类型适宜性概率之和均为1，即：

∑ksq(p,k,t)=1

(2)

(3)

(4)

(5)

1.3.3 模型验证 通过模拟2000-2015年的土地利用变化，测试和比较FLUS模型的性能，其他输入包括模拟参数、空间限制区域和适宜性概率数据。本研究以2000年土地利用空间分布图为模拟基准图，在FLUS模型中运行得到2015年土地利用模拟分布图，使用Kappa系数和FoM系数来评估模型的可靠性，以衡量模拟结果与实际土地利用数据之间的一致性[24-25]。该模型的总体精度为83.25%，Kappa值和FoM值分别为0.756和0.314，说明该模型参数符合精度要求，可以用于预测未来的土地利用格局。因此，使用2015年的有效参数和土地利用分布图，预测2030年3种情景下土地利用分布情况。

1.4 生态系统服务的量化

选取InVEST 3.7.0版本的固碳和产水量服务模块，基于2015年和2030年不同情景下的土地利用数据，定量评估北京市生态涵养区的生态系统服务[26-27]。每个参数的数据基于已发布的共享数据库。

1.4.1 固碳服务 根据土地利用当前时间段内的碳储量数据，利用InVEST 3.7.0的“碳储存和封存：气候调节”工具，对每个情景的碳储量进行计算。该模型使用土地利用地图和特定分类汇总了4个主要碳库中的碳储量：地上生物量、地下生物量、土壤和死亡有机物[28]。土地利用类型j中栅格单元(x，y)的总碳CSjxy公式为：

CSjxy= Ax (Cajxy+Cbjxy+Csjxy+Cdjxy)

(6)

式中，Cajxy、Cbjxy、Csjxy、和Cdjxy分别是土地利用类型j中栅格单元(x，y)的地上碳密度/(Mg·hm-2)、地下碳密度/(Mg·hm-2)、土壤有机碳密度/(Mg·hm-2)和死亡有机物碳密度/(Mg·hm-2)。本研究假设碳在土地覆盖变化期间会立即变化。

1.4.2 产水量服务

1.4.2.1 评估模型 利用InVEST 3.7.0模型的“季节产水量”模块，计算每个网格单元的年降水量(mm)、土壤最大根系埋藏深度(mm)、植物可利用水含量(0%～1%)、年平均潜在蒸散发量(mm)、土地利用[29-31]。土地利用类型x上每个栅格单元的年产水Yx的计算公式如下：

Yx=(1-AETx/Px)·Px

AETx/Px=(1+PETx/Px)-[(1+(PETx/Px)ω]1/ω

PETx=Kx·ET0/Px

ωx=Z·AWCx/Px+1.25

AWCx=Min(SoilDepth，RootDepth)·PAWC

(7)

式中，AETx表示栅格x的年实际蒸散量，Px为栅格x的年降水量，AETx/Px为土地利用类型的植被蒸散发，PETx为潜在蒸散量，ωx为自然气候-土壤性质的非物理参数，ET0为栅格单元x的参考作物蒸散，Kx为栅格单元x中特定土地利用/覆被类型的植物蒸散系数，AWCx为土壤有效含水量，PAWC为植物可利用含水率(0%～1%)。

1.4.2.2 模型参数的选择Z值为经验常数，又称“季节常数”，能够表征区域降雨分布及其他水文地质特征，参考北京地区的相关文献，取Z值为4.02[31-32]。根据中国植被类型分布图，研究区的土地利用类型主要有草地、水域农田、建设用地、未利用地、林地和灌木林地等类型，参考相关文献和InVEST模型使用手册中的研究结果，确定各土地利用类型的植被蒸散指数Kc(表1)[33-35]。

表1 不同土地利用类型的植被蒸散指数KcTable 1 The plant evapotranspiration coefficient (Kc) of different land use types

2 结果与分析

2.1 土地利用变化特征

生态涵养区在2015-2030年不同情景下的土地利用变化特征有较大差异(表2和图3)。林地和灌木林地是研究区最主要的土地利用类型，在2015年，其面积分别占研究区总面积的44.14%和32.10%，其次是农田(8.36%)以及建设用地(6.38%)，主要分布在研究区的西部和东部边缘地区。此外，水域和未利用地面积所占比例均较小，均占5%以下。在2015-2030年城市快速发展情景下，土地利用变化的主要特征是建设用地的快速扩张，相对于2015年，城市快速发展情景的建设用地增加了39.57%，远高于自然演变情景和生态控制情景的15.10%和2.86%。农田和水域在2015-2030年自然演变情景的下降幅度最大，均高于20%。在2030年的生态控制情景下，建设用地和林地面积的变化较为稳定，相对于2015年，分别增加了2.86%和3.10%，农田的减少率低于自然演变情景和城市快速发展情景(-8.74%)。

表2 研究区2015-2030年3种情景下土地利用变化Table 2 Land use changes from 2015 to 2030 under three scenarios in ecological conservation area

注：(a)：2015年；(b)：自然演变情景；(c)：生态控制情景；(d)：城市快速发展情景。下同。图3 2015和2030年3种情景下的土地利用覆盖Fig.3 Land use maps from 2015 to 2030 underthree scenarios

2.2 生态系统服务变化特征

2.2.1 固碳服务 研究区碳储量的最大值是生态控制情景下的99.53×106Mg，最小值是城市快速发展情景下的96.31×106Mg，相对于2015年的98.36×106Mg，分别变化了1.17×106Mg和-2.05×106Mg(表3)。从土地利用类型的碳储量及其变化特征可以看出，林地拥有最大的碳储量，其中，生态控制情景下林地的碳储量最大(66.80×106Mg)，相对于2015年增加了1.27×106Mg，其余如水域、建设用地和未利用地的碳储量均为0。就空间分布特征来看，生态涵养区中部主要分布的是森林和灌木林地，其碳储量较为丰富，西侧的延庆城区和东部边缘地区主要分布着建设用地和农田，碳储量相对较少，因此可以看出，碳储量呈现出中部多东西偏少的基本分布格局(图4)。结合研究区地形地貌特征来看，高水平的碳储量主要分布在海拔相对较高的区域，而地势较为平缓的区域碳储量也相对较少(图1)。

图4 2015年-2030年3种情景下的碳储量空间分布格局Fig.4 Spatial distribution of carbon storage from 2015 to 2030 under three scenarios

表3 研究区2015-2030年3种情景下的碳储量及其变化情况Table 3 The variation ofcarbon storage from 2015 to 2030 under three scenarios ×106Mg

2.2.2 产水量服务 研究区产水量的最大值是城市快速发展情景下的362.50×107m3，最小值是生态控制情景下的361.5×107m3，相对于2015年的361.08×107m3，分别变化了1.42×107m3和0.42×107m3(表4)。从土地利用类型的产水量及其变化特征来看，林地拥有最大的产水量，其中，在生态控制情景下的林地产水量最大(161.13×107m3)，相对于2015年增加了5.27×107m3。城市快速发展情景下的建设用地产水量为36.81×107m3，相对于2015年增加量最大(9.98×107m3)。从产水量的分布特征来看，研究区产水量呈现出南高北低以及中部低东部和西部偏高的基本分布格局(图5)。结合研究区的地形地貌特征来看，高产水量区域主要分布在海拔相对较低的区域，而地势较高的区域产水量也相对较少(图1)。此外，城市和农田被确定为提供产水量服务的热点地区。

图5 2015-2030年3种情景下的产水量空间分布Fig.5 Spatial distribution of water yield from 2015 to 2030 under three scenarios

表4 研究区2015-2030年3种情景下的产水量及其变化情况Table 4 The variation of water yieldfrom 2015 to 2030 under three scenarios ×107 m3

2.3 土地利用变化对生态系统服务的影响

从生态涵养区2015-2030年土地利用与碳储量和产水量的变化特征来看，碳储量与产水量变化同各用地类型的面积变化具有一定关系(表5)。从单位面积变化对碳储量影响的绝对量来看，灌木林地的影响强度在城市快速发展情景中达到最大，其面积每减少1 km2，碳储量损失3.91×104Mg。从单位面积变化对产水量影响的绝对量而言，未利用地的影响强度在自然演变情景中最大，其面积每减少1 km2，产水量减少4.93×105m3。对比各用地类型面积增减幅度以及与碳储量和产水量的变化关系发现，尽管面积同碳储量和产水量呈正相关关系，但面积增减仅意味着碳储量和产水量的增减，其增减幅度并不能表征碳储量和产水量的增减幅度。

表5 2015-2030年3种情景下土地利用变化对生态系统服务的影响Table 5 Impacts of land use changes on ecosystem services under three scenarios from 2015 to 2030

3 结论与讨论

3.1 结论

利用FLUS模型耦合马尔可夫链(Markov chain)模拟了2030年自然演变、生态控制和城市快速发展情景下的土地利用分布，并探讨了2015-2030年3种土地利用3种情景下的土地利用变化如何影响固碳和产水量服务。

北京市西北部生态涵养区2015-2030年的土地利用变化程度因不同情景而出现较大差异，其中，相对于2015年，城市快速发展情景的建设用地面积增加了39.57%，生态控制情景下的建设用地和林地面积的变化较为稳定，分别增加了2.86%和3.10%，农田的减少率低于自然演变情景和城市快速发展情景(-8.74%)。

研究区内碳储量和产水量的最大值分别出现在生态控制情景(99.53×106Mg)和城市快速发展情景(362.50×107m3)。林地拥有最大的碳储量和产水量，其中，生态控制情景下的碳储量(66.80×106Mg)和产水量(161.13×107m3)最大。碳储量在生态涵养区的中部较为丰富，在西侧的延庆城区和东部边缘地区的建设用地和农田较少，而产水量呈现出东部和西部偏高其余地区偏低的分布格局。高水平的碳储量主要分布在海拔相对较高的区域，而高水平的产水量区域主要分布在相对较低的区域。

碳储量与产水量变化同各用地类型的面积变化具有一定关系，其中，灌木林地单位面积变化对碳储量的影响强度在城市快速发展情景中最大，未利用地单位面积变化对产水量的影响强度在自然演变情景中最大。

3.2 讨论

3.2.1 土地利用变化对生态系统服务的影响 研究结果显示，在2015和2030年的3种情景中，生态控制情景下的林地面积占比最大，其单位面积的碳储量最高，建设用地、水域和未利用地的碳储量为0，说明林地作为固碳服务的重要供给区，对生态涵养区至关重要。此外，城市快速发展情景下的碳储存总量最小，主要是由于建设用地并不是有效的碳库，其面积的扩张会导致诸如农田和林地等面积的减少[36]。

林地、灌木林地和草地的单位面积产水量较低，而建设用地和农田生态系统的单位面积产水量较高，这与窦攀烽等[37]的研究结果相一致。有研究表明，产水量主要受植被的蒸散作用、根系深度和下渗过程的影响，而建设用地属于不透水地表，其蒸散量和下渗水量都较小，降水到达地表后易迅速形成径流，其产水量高于其他土地利用地类型，因此，建设用地面积在城市快速发展情景下最大，进而增加了产水总量[38]。尽管如此，降水到达不透水地表后易快速形成径流进入城市下水道，难以被人类利用，因此，水资源的有效供给也会受技术因素的影响[38]。农田的植被蒸散量与森林和草地类似，但由于根系深度和植被密度等原因，农田的下渗水量明显小于森林和草地，因此，农田的单位面积产水量高于森林和草地。

土地利用类型面积的扩大或缩小会导致区域碳储量和产水量不同程度的增加或减少，该结论是基于每个用地类型内部均质的假定，仅从各用地类型量变的角度进行了研究，并未考虑质的差别，如果将用地类型内部的差别考虑在内，土地用地类型面积的变化不一定能够产生同面积等量的碳储量和产水量的变化，如林地每增加1 km2，碳储量不一定新增0.14×104Mg。

3.2.2 模型优点与局限性 本研究方法主要包括土地利用的情景模拟和生态系统服务评估，其优点是为探讨土地利用变化对3种生态系统服务的影响提供了一种可供借鉴的方法，该方法对于其他地区的研究也是有效的。此外，InVEST模型被认为是适用于多个尺度的生态系统服务评估，可较为直观地表达评估结果[26,28]。同时，InVEST模型也存在一些局限性，首先，生态系统服务的评估结果过于依赖土地利用分类，如本研究将土地利用类型分为7个大类，每个土地利用类型的内部均被假定为均质的，并没有考虑内部的异质性[39]。例如，森林的碳储量受温度、海拔、降雨和森林年龄结构的影响，而在划分较为粗糙的土地利用类型中，碳储量的异质性特征被部分掩盖了。其次，尽管土地利用数据是由30 m分辨率的遥感影像解译而来，但土壤、降水、温度等模型参数数据均来源于1 000 m分辨率的数据源，因此，数据的不匹配增加了生态系统服务评估的不确定性[9]。此外，“季节性产水量”模块是基于年平均数据，忽略了极端情况，没有考虑季节性的供水时间模式[40]。

3.2.3 政策建议 根据2000-2015年的土地变化趋势预测了2030年自然演变情景下的土地利用状况，并以此为基础设计了生态控制和城市快速发展2种土地利用情景。针对《北京城市总体规划(2016-2035年)》中对森林、水域和农田的保护规划，本研究在城市快速发展和生态控制情景的设置中禁止将水库、自然保护区和国家森林公园转换为其他用地。结果表明，与自然演变情景相比，建设用地的扩张在生态控制情景下将得到有效控制。此外，生态控制情景是具有更可持续性的发展策略，因为其将2030年人类为满足生产、生活和其他需求所需的土地考虑在内，因此，加强自然资源保护规划和调控的实施，是缓解生态系统退化的有效途径。

此外，本研究识别了固碳和产水量服务及其响应土地利用变化的热点区域，为维持和加强生态系统服务保护提供了有效的参考意见。例如，土地利用变化结果表明，如按自然演变情景和城市快速发展情景，农田和水域在接下来的几十年将大面积的减少，而建设用地将大幅度的增加，这也为如何防止潜在的城市无序扩张以及保护基本农田提供了依据[41]。西部的延庆城区以及东部边缘地带适当的增加林地和农田将促进固碳和产水量服务的区域空间平衡，因此，多种管理策略指导土地利用和其他资源(如粮食，石油和技术)的分配显得尤为必要。此外，尽管城市的产水量较大，但大部分降水进入城市下水道后难以被人类利用，因此，如何有效利用降水和提高水资源的有效供应是关键。最后，建议政策制定者应采取多样化的土地管理策略来指导不同情景下的土地资源的配置，从而实现生态系统服务的可持续利用。