城市边缘区作为城市与乡村的过渡地带和未来城市的扩展区域^[1], 是城市化进程中生态环境比较敏感的地区, 在城市发展中占有重要地位。很多学者研究认为城镇化对生态环境带来了负面影响, 使生态系统服务价值降低^[2], 例如, 城镇扩张、道路建设、工业发展等破坏自然植被、打断水文系统^[3-4]。此外, 城市边缘地区集体产业虽然发展快, 但规模小且数量多, 导致很多农业用地^[5]和生态用地被违法占用, 造成整个生态系统的结构和功能改变。但是, 人为干扰下造成土地利用结构和功能的改变是否一定会使边缘地区生态系统服务价值降低, 目前尚无定论。例如：Pickett等^[6]认为扩大城市空间有利于提高生物多样性。因此, 在国家大力推进生态文明建设和美丽乡村建设的背景下, 探索城市边缘地区景观格局变化及其生态系统服务功能价值影响具有重要现实意义。

定量评估生态系统服务价值已成为国际可持续发展研究的热点之一^[7], 是当前生态学、生态经济学、环境科学以及土地管理科学的交叉前沿领域。对于生态系统服务价值的核算, 在国内主要是以谢高地等^[8]的研究为基础, 根据不同区域实际情况调整生态服务当量因子表, 从而估算研究区生态系统服务价值, 能够比较科学地反映研究区生态系统服务价值情况, 为区域生态环境建设提供依据^[9]。在研究尺度上, 当前对于生态系统服务价值的研究主要集中于大尺度、多区域土地利用变化对生态系统服务价值的研究, 例如：刘金勇等^[10]探讨了济南市土地利用变化及其对生态系统服务价值的影响；赵丹等^[11]以淮北市为例, 阐述了1987—2007间城市土地利用变化对生态系统服务的影响；刘桂林等^[12]分析了整个长三角地区1980—2010年的土地利用变化对生态系统服务价值的影响。但将土地利用格局变化运用到量化城市边缘区生态系统服务价值鲜有报道。相比于大尺度研究中复杂的政策和土地利用变化, 小尺度研究更容易发现产生变化的驱动机制^[13]。

根据《2004—2020年北京城市总体规划》, 顺义新城是东部发展带的重要节点和重点发展的新城之一, 带动区域发展的临空产业中心和先进制造业基地, 是首都国际航空中心核心区。其中牛栏山-马坡镇是顺义新城建设的核心区, 又是未来顺义发展的产业核心区。在持续的快速经济发展背景下, 两地土地利用格局已发生并有可能持续发生剧烈的土地利用格局变化。因此, 研究生态系统服务价值对当地可持续发展具有重要意义。本研究旨在分析与探讨社会经济快速发展的背景下, 牛栏山-马坡镇景观格局的变化趋势及其引起的生态服务价值的动态变化和驱动机制。本研究主要创新之处是分析了在快速发展地区, 时间变化造成的土地利用空间差异和结构变化带来的生态服务价值变化, 并判断整个区域生态服务价值增减的主要土地利用类型及原因。因此, 本研究不仅为北京边缘地区土地资源优化配置、土地利用效益评价和可持续发展提供支撑, 而且通过分析土地利用类型景观指数与生态系统服务价值的关系, 可为其他类似城市边缘地区提高生态服务价值提出针对性对策提供理论依据。

1 研究区与研究方法 1.1 研究区概况

牛栏山-马坡镇(40°00′—40°18′ N, 116°28′—116°58′ E)位于北京市外边缘顺义区^[14](图 1), 共有23个行政村, 占地面积7029 hm²。该区域属暖温带半湿润大陆性季风性气候, 年平均气温为11.5℃, 最低气温零下19.1℃, 最高气温达40.5℃。年日照2750 h, 无霜期195 d左右。年均相对湿度50%, 年均降雨量约625 mm, 是华北地区降水量较均衡的地区之一, 全年降水的75%集中在夏季。其中, 属地财税收入马坡乡由2004年的0.61亿元增长到2014的11.73亿元。牛栏山镇由2009年的3.53亿元增长到2013年的11.26亿元, 经济跨越式增长明显。

1.2 数据来源

研究区范围为牛栏山-马坡镇, 包括1992年、2001年、2015年3期土地利用数据, 其中1992年和2001年土地利用数据来源于中国科学院遥感与数字地球研究所(http://www.radi.ac.cn/)和北京市国土局1 : 1万土地利用现状图, 2015年土地利用数据根据worldview 2高分影像目视解译得到。参考我国GB/T21010—2007《土地利用现状分类》和生态用地相关分类体系^[15], 将研究区景观类型分为2个一级分类, 分别为建设用地和生态用地, 其中生态用地又包括农田、林地、水体、城市绿地、未利用地。具体分类见表 1。

1.3 研究方法 1.3.1 景观格局分析

景观格局的变化直接影响栖息地质量^[16]以及营养元素食物迁移路径^[17-18], 间接改变生态系统中的能量流动和物质循环过程。因此, 对景观格局进行定量分析并研究其动态特征, 是理解生态过程的基础^[19]。由于景观指数种类较多, 且往往是高度相关的, 所以本文以Fragstat 4.2软件作为技术平台获取景观格局变化动态指数, 根据不同功能及属性选取斑块个数(NP, Number of Patches)、斑块密度(PD, Patch Density)、最大斑块指数(LPI, Largest Patch Index)、景观形状指数(LSI, Landscape Shape Index)、香农多样性指数(SHDI, Shannon′s Diversity Index)、散布与并列指数(IJI, Interspersion Juxtaposition Index)和蔓延度指数(CONTAG, Contagion Index)进行分析。其中, 斑块个数(NP)反映景观的空间格局, 经常用来描述整个景观的异质性, 其值的大小与景观的破碎度有很好的正相关性；斑块密度(PD)反映景观总体斑块的分化程度或破碎化程度, 斑块密度高表示单位面积斑块规模小, 景观异质性高；最大斑块指数(LPI)是景观中最大斑块面积与总面积之比, 有助于确定景观优势类型；景观形状指数(LSI)是景观边界总长度与总面积平方根之比, 值越大, 表明类型越复杂, 斑块越分离；香农多样性指数(SHDI)反映景观类型丰富和复杂程度, 值越大表明景观中没有优势类型且各景观类型趋于均匀分布；散布与并列指数(IJI)可以用来表征景观总体分散情况, 指标越大不同斑块交替出现的规律越明显；蔓延度指数(CONTAG)一般用来表征景观中不同斑块类型的团聚程度或蔓延趋势, 高蔓延值表明景观类型中的某种优势斑块类型形成了良好的连接性, 反之则表明景观是具有多种要素的密集格局, 景观的破碎化程度较高^[20]。

1.3.2 土地生态系统服务价值估算

(1) 生态系统服务价值系数

在国内, 对于生态系统服务价值的评估主要采用谢高地等^[8]的中国生态系统服务价值当量因子表。由于本研究中土地利用类型划分与谢高地等的划分不完全一致, 需要将两者的生态系统类型进行对比, 并根据其他生态系统进行调整。其中, 由于建设用地是一个复杂的生态系统, 包括乔木、灌木、草地等组成的绿地系统和水体组成的流域生态系统等, 与程琳等^[20]的研究中土地利用分类相似, 且都是基于城市生态服务价值的评估, 因此, 参考程琳等^[21]的方法调整建设用地生态系统的当量因子, 并将其他用地因转化为建设用地而减少的损失纳入建设用地生态系统服务价值的进行计算。由于没有1992年建设用地绿地面积, 本研究设定当年建设用地生态服务价值为0。

(2) 生态系统服务价值估算

根据各土地利用类型的生态系统服务价值系数, 基于谢高地等^[8]的研究, 生态服务价值系数=生态系统服务价值当量×全国平均粮食单产市场价值/7, 分别计算1992—2015年牛栏山-马坡镇各用地类型的生态系统服务价值, 计算公式如下：

式中, ESV_n为某类生态系统服务的价值量, V_in为第i类用地类型的第n种生态系统服务的价值系数, A_i为一种用地类型的面积。

1.3.3 景观指数与生态系统服务价值相关性分析

建立景观指数与生态系统服务价值的耦合关系, 能够较好反映景观指数变化对生态系统服务价值的影响。顾泽贤等^[22]以澜沧县为研究对象, 利用2000年、2005年和2010年影像的解译结果, 建立并分析了景观格局变化与其ESV的相关性关系。基于顾泽贤的研究, 本研究利用SPSS软件, 选用Spearman相关系数, 采用1992年、2001年和2015年3期影像解译出的结果, 建立景观指数与ESV的关系, 分析研究区景观指数变化对ESV的影响。

2 结果与分析 2.1 牛栏山-马坡镇1992—2015年景观格局空间变化特征

根据牛栏山-马坡镇1992、2001、2015年3期土地利用图(图 2), 并参照上表中土地利用一级分类, 得到2002—2015年牛栏山-马坡镇不同景观类型面积变化及空间格局差异。总体来看, 1992—2015年期间, 牛栏山-马坡镇建设用地呈现向四周扩张的趋势, 农田大面积收缩, 林地面积向四周显著扩张；水域面空间分布由分散到集中, 从分散在北部和东部到集中到东南部。

2.2 景观格局指数变化分析

1992—2015年牛栏山-马坡镇景观斑块数量(NP)、斑块密度(PD)、香农多样性指数(SHDI)呈增长趋势, 且景观斑块数量(NP)年增长均值有减弱趋势；说明在1992—2015年虽然景观破碎化程度加剧, 异质性不断增高, 但有减弱趋势, 景观类型趋于丰富多元化, 斑块类型趋于均匀分布。景观蔓延度指数(CONTAG)、散布与并列指数(IJI)呈递减趋势, 且景观蔓延度指数(CONTAG)递减趋势减弱；说明研究区的景观破碎化程度升高, 但升高趋势有所减弱, 不同斑块类型的邻接概率降低；最大斑块指数(LPI)先减小后增大, 但景观形状指数(LSI)却正好相反, 表明牛栏山-马坡镇1992—2015年优势景观类型发生了改变, 且景观形状发生剧烈改变。可知研究区的景观格局受到人类活动的严重干扰(表 2)。

2.3 不同景观类型面积动态变化

图 3表明1992—2015年农田面积不断减小, 建设面积、林地不断增大, 尤其在2001—2015年变化趋势更加明显, 在2015年建设用地已成为优势景观类型, 水体和未利用地面积变化不大, 建设用地和林地主要由不断减小的农田转化而来。1992—2001年北京边缘区牛栏山-马坡镇城镇用地面积从1988.7 hm²增加到2078.54 hm², 增长率为4.52%, 城镇化并不明显, 但2001—2015年牛栏山-马坡镇受人为干扰严重, 城镇用地面积从2078.54 hm²增加到3091 hm², 增长率达32.76%, 城镇化速度明显加快。

2.4 不同景观的生态系统服务功能价值分析

依据牛栏山-马坡镇3期土地利用变化数据, 结合不同土地利用类型的生态服务价值系数(表 3), 估算牛栏山-马坡镇景观格局变化引起的生态系统服务价值(ESV)变化(图 4、图 5)。结果表明：较1992—2001年, 2001—2015年生态系统服务总价值增加了2528.39万元；其中, 林地ESV增加明显, 从1992年的1109.67万元增加到2015年的8940.98万元。水域、未利用地ESV变化不明显, 农田和建设用地ESV逐年减小。在1992—2001年水体提供的ESV最多, 是牛栏山-马坡镇提供生态系统服务的主要土地类型；但在2001—2015年, 林地所提供的ESV超过水域, 成为提供研究区ESV的主要土地类型。1992—2015年随着景观格局改变, 牛栏山-马坡镇食物生产和废物处理的ESV减少, 其他各项生态服务功能价值均呈增加趋势。由于农田主要提供食物材料等价值, 从侧面说明农田面积不断较少。

2.5 不同景观类型对生态系统服务价值贡献率分析

从不同时期各景观类型对生态系统服务总价值的贡献率来看(图 6), 1992—2015年, 林地对生态系统服务贡献率最大, 农田和建设用地贡献率最小, 水体和未利用地对生态系统服务总价值贡献率并没有显著变化。这表明：研究区ESV提高的主要原因是由于林地的增加, 农田面积的减少和建设用地的增加是造成生态系统服务总价值减少的原因。对于建设用地, 2001—2015年相比1992—2015年对ESV的贡献率明显增大, 表明随着居民生活水平的提高, 牛栏镇-马坡镇生态环境有所改善。

2.6 景观指数与生态系统服务价值相关性分析

根据牛栏山-马坡镇1992—2015年不同景观指数与生态服务功能做关联度分析, 得到相关系数(表 4), 由表 4可知：在选择的6类景观指数中, 生态系统服务总价值与香农多样性指数(SHDI)、斑块密度(PD)、最大斑块指数(LPI)、散布与并列指数(IJI)、蔓延度指数(CONTAG)相关性较高, 说明景观类型越丰富, 破碎度越大有利于提高整体服务价值, 景观类型越分散可能不利于ESV的提升。在研究区虽然连通性逐渐降低, 生态系统服务价值却逐年升高, 林地生态系统服务总价值与斑块密度、最大斑块指数、景观形状指数均呈高度正相关, 与散布与并列指数呈高度负相关, 表明提高林地斑块密度、最大斑块指数、景观形状指数有利于提高林地ESV, 进而有利于提高总ESV；香农多样性指数、斑块密度、散布与并列指数、蔓延度指数与各项生态系统服务功能之间存在较高相关性, 提高香农多样性指数、斑块密度对提高原材料生产、气体调节、气候调节、水文调节、保持土壤、维持生物多样性和提供美学景观价值有较大影响, 提高散布与并列指数、蔓延度指数有可能提高食物生产和废物处理价值。

3 讨论 3.1 生态系统服务价值变化分析

通常认为, 建设用地增加等人为扰动会造成ESV降低。例如, 王佼佼^[22]在北京ESV研究指出, 由于建设用地急剧增长, 耕地和水域面积大幅度减少, 造成北京顺义区2000—2008年ESV降低；张轶秀等^[23]对广州城市边缘区ESV变化进行研究, 发现在2000—2008年ESV随着经济的发展也逐渐下降；马凤等^[24]对同属于一线城市的上海临港新城ESV变化的研究发现, 由于人为干扰和土地利用的改变, 2000—2006年ESV逐渐降低。总体来看, 由于经济的发展, 在2000—2008年期间一线城市边缘地区ESV普遍降低。与此形成对比的是：牛栏山-马坡镇1992—2015年建设用地虽然不断侵占农田和生态用地。但ESV并没有减少, 反而增加3382.73万元, 且2002—2015年增幅加大。这是由于虽同处城市边缘区, 但经济发展和政策导向存在城市间差异。本研究中牛栏山-马坡镇ESV增大主要原因包括：(1)牛栏山-马坡镇大量农田转化为林地, 例如, 苗圃面积增加1394 hm², 占林业总增长量的81.25%, 而林地的生态价值系数相对较大, 林地对ESV起决定性作用；(2)研究区植被覆盖率显著提高。由于2001年实施的退耕还林工程、水上奥运场馆及配套绿化工程、公路绿色交通体系建设、生态廊道建设工程、百村万户绿化美化工程、首都平原百万亩造林工程等使牛栏山-马坡镇植被覆盖率增大, 林地面积增加, 使得植被固碳释氧、吸热降温、空气和水质净化等生态功能, 因此对气体调节、气候调节、保持土壤价值有显著提高作用^[25]。同时增大了截留降水和补充地下水的能力, 对水源涵养, 水文调节功能有较大贡献；(3)在经济拉力上, 由于果业在收益上高于粮食作物, 所以传统耕地更多用于种植经济林, 直接增加了植被覆盖度^[26], 从而提高ESV。整体而言：牛栏山-马坡镇林地和园地带来的生态经济效益弥补了由于建设用地的扩大造成的ESV的损失；加强林业建设、扩大园地面积、加大建设用地绿地覆盖率是提高牛栏山-马坡镇生态服务价值的主要措施。因此, 牛栏山-马坡镇生态环境建设模式可为其他边缘地区城市发展提供参考。

3.2 景观格局对生态系统服务价值影响分析

景观格局演变与生态系统服务变化之间存在着密切相关性^[27]。第一, 景观功能决定生态过程, 由于气候、天气状况等会影响景观功能, 景观功能的改变影响能量的流动和营养物质的循环, 从而影响生态系统服务价值。例如, 农田耕作、施肥等改变土壤物理化学性质, 影响土壤对有机物的储存, 增加非点源污染的发生。最终影响生态系统服务价值^[28]。第二, 景观格局的改变影响景观指数的变化, 从而影响生态系统服务价值。王云等^[29]以西安市为例, 利用1999年、2004年和2011年影像数据, 研究了都市农业景观破碎化过程对ESV的影响, 表明最大斑块指数与各项ESV均呈较强的正相关关系, 但斑块密度和分离度呈负相关。张明阳等^[30]分析了桂西北典型喀斯特区ESV与景观格局的关系, 表明最大斑块指数、蔓延度指数与ESV呈曲线正相关, 随着连通性的增强, ESV有所增加；分离度指数、斑块丰富度与ESV呈曲线负相关, 随着斑块的破碎与分离的加大生态服务价值有所减少。岑晓腾^[25]以杭州湾南岸区域为例, 研究了景观格局与ESV的关系, 表明增加景观多样性、提高破碎度和分散度有利于提高整体服务价值。从目前研究来看, 不同区域, 由于景观格局的不同, 景观指数与ESV的关系不同, 每种景观指数对生态系统服务影响的程度有待于深入研究。本研究结果表明牛栏山-马坡镇受城镇化影响, 斑块密度增大, 最大斑块指数先减后增, 蔓延度不断减小, 破碎度增大, ESV不断增大。

鉴于以上分析, 提出以下几点对策建议, 对未来北京边缘区发展规划提供参考：

(1) 潮白河与牛栏山镇紧密相连, 且位于牛栏山镇内, 潮白河是顺义新城发展的生态屏障。因此, 对于毗邻水系的城市边缘区, 要加强对林地的保护力度, 要严格控制耕地及其他地类对林地、水域的侵占, 合理利用国家政策, 提升牛栏山镇潮白河作为生态屏障的服务功能。

(2) 牛栏山-马坡镇在完成耕地保护指标的同时, 要保护提高全域林地和水域面积, 注重维持现有苗圃面积和经济林面积, 通过各种方式提升牛栏山-马坡镇生态服务价值。

(3) 牛栏山-马坡镇要控制建设用地比例, 并且提高建设用地绿地覆盖率, 缓解由于建设用地面积增加, 生态服务价值剧烈减小的趋势。

(4) 对于北京市边缘区要通过提高斑块面积及数量, 提高整个区域景观格局斑块密度指数、最大斑块指数、景观形状指数和多样性指数, 注重提高林地密度指数、最大斑块指数、景观形状指数, 以提高生态服务价值, 改善生态环境。

4 结论

本文以牛栏山-马坡镇1992年、2001年和2015年3期土地利用数据为基础开展城市边缘区景观格局变化及其生态服务价值影响分析, 年份跨度大, 尺度小, 能够较好的反映城市边缘地区小尺度上景观格局变化及其对ESV空间变化影响, 并且构建了“土地利用景观指数—生态系统服务价值”相关性分析, 找出了决定牛栏山-马坡镇生态系统服务价值变化的主要土地利用类型以及相关指数, 可以较好地分析自然生态系统与人类社会系统之间的驱动和反馈机制, 结果表明：

(1) 1992—2015年牛栏山-马坡镇受到人类强烈干扰, 景观格局发生显著变化, 破碎化程度加剧, 异质性不断提高, 且景观类型趋于丰富多元化, 斑块类型趋于均匀分布, 优势景观类型由农田转变为建设用地。

(2) 1992—2015年牛栏山-马坡镇ESV总量不断增大, 其中1992—2001年牛栏山-马坡镇的生态系统服务总价值变化不大, 仅增加427.17万元, 但2001—2015年生态系统服务总价值增加了2955.56万元, 增幅明显。除食物生产和废物处理价值减少外, 其他各项生态服务功能价值均呈增加趋势。主要提供ESV的优势景观类型由水域变为林地。

(3) 从景观格局对ESV影响分析得知：不同区域景观指数与生态系统服务价值相关性不同, 在城市边缘区林业对ESV起着决定性作用, 提高林地斑块数量及面积有利于提高林地生态服务功能, 进而提高生态系统服务总价值。加强林业建设, 提高建设用地绿地覆盖率是克服建设面积扩张而导致生态服务功能降低的有效措施。