京津冀城市群生态系统服务价值的时空动态演变

汪东川 ，孙志超，孙然好，陈俊合，张威，张仙，王潇

1. 天津城建大学地质与测绘学院，天津 300384；2. 天津城建大学/天津市土木建筑结构防护与加固重点实验室，天津 300384；3. 中国科学院生态环境研究中心/城市与区域生态国家重点实验室，北京 100085

生态系统作为地球生命与环境的统一，不仅是人类生存和发展的基础，而且为生态系统服务与功能的形成和维持提供了物质基础（李晓赛等，2015；虞依娜等，2010）。生态系统服务指生态系统形成和维持的人类赖以生存和发展的环境条件与效用（欧阳志云等，1999；Lin et al.，2018），并主要通过供给服务、调节服务、文化服务、支持服务等服务功能为人类提供直接或间接的利益（熊善高等，2016；陈能汪等，2009），其价值的定量评估结果对资源管理决策、维持生态平衡、促进区域经济与环境的协调发展具有重要意义（赵苗苗等，2017；陶宝先等，2010）。

自 20世纪末期以来，生态系统的可持续发展与生态安全日益受到重视，生态系统服务价值定量评估的研究已成为生态经济研究的热点（Sawut et al.，2013）。最初，Costanza（1998）、Daily et al.（2000）、Groot et al.（2002）国外学者对全球的生态系统服务进行了全面系统的研究，在此基础上，谢高地等（2001，2003）、欧阳志云等（1999a，1999b）、赵景柱等（2000）为代表的国内学者完善了生态系统服务的相关理论、明确了生态系统的主要功能、构建了适宜于中国生态系统服务价值的评价方法，并针对不同区域、不同尺度、不同类型的生态系统服务价值进行了评估，其结果对生态环境保护、生态功能分区等多个方面产生了一定的社会效益（Gong et al.，2017；Zhang et al.，2015；Zhang et al.，2015）。其中，谢高地于 2003年基于我国200多位生态学者的问卷调查内容提出的适合我国生态环境的“单位面积生态系统服务价值当量表”被广泛接受，并于 2008年从降低水域生态系统服务价值的角度修正了该当量表（谢高地等，2008）。修正后的方法本质依然是基于土地利用/覆被变化进行当量赋值。然后结合单位面积生态系统服务价值系数与各类生态系统面积实现对生态服务价值的定量评估，后来又于2015年构建了NPP、降水、土壤保持空间调节因子，探究了中国生态系统服务价值在省/直辖市尺度上的异质性（谢高地等，2015）。在此基础上，刘金雅等（2018）基于泰森多边形修正的市域边界分析了京津冀城市群在地级市尺度上的异质性，提高了生态服务价值的空间异质性精度。随着研究的深入，该方法依旧存在修正的空间。首先，该方法是问卷调查的结果，对区域内生态系统服务功能及价值的客观性反映明显低于主观性反映（杜加强等，2008；李伟等，2014）；其次，该方法适宜于较大尺度的整体区域评估，但对区域内生态服务价值空间差异性的反映精确度不足，尤其是同种土地利用类型的生态系统服务价值的空间差异性（Zhang et al.，2015）；最后，该方法以土地利用类型及面积属性作为核心评估标准，但是关于土地利用内部结构与整体格局对生态系统服务价值影响的研究较少（王航等，2017；赵江等，2016；Yirsaw et al.，2017），且多数研究采用的生态系统服务价值定量评估方法，以静态分析为主，对生态系统服务价值受自然和经济条件变化影响的动态演变研究不完善（Kreuter et al.，2001；Li et al.，2010）。

京津冀城市群作为推动华北地区社会经济发展的增长极，其快速的城市土地扩张模式改变了生态系统结构与功能，并严重制约了京津冀城市群的可持续发展（景永才等，2018）。2015年《京津冀协同发展规划纲要》的发布以及 2017年雄安新区的设立，不仅强调了京津冀的协同发展，而且把“坚持生态优先、绿色发展”放在雄安新区开发建设的首要位置（毛汉英，2017；徐涵秋等，2017）。鉴于此，本文从京津冀城市群生态系统服务一体化的角度出发，基于土地利用类型数据以及“中国生态系统单位面积生态服务价值当量”，利用3 km×3 km的格网分析方法，并结合县域尺度的粮食单产、NPP、降水量、土壤保持模拟量数据对传统ESV评估模型进行改进，在此基础上，利用热点分析方法与等量分类方法对生态系统服务价值的空间分布特征进行分析，并将空间格局指数、动态度、变化轨迹分析等方法相结合，探究生态系统服务功能与价值的时空动态演变，旨在为京津冀城市群区域生态系统服务功能增强、土地利用结构调整、生态景观建设提供科学依据。

1 研究区概况

图1 研究区地理位置Fig. 1 Geographical location of study area

京津冀城市群地处环渤海核心地带，位于113°34′－120°05′E，36°00′－42°40′N 之间（图 1），是由北京市、天津市以及河北省 11个地级市组成的城市群，土地总面积约216485.2 km2。境内地形以平原和山地为主，南北高程差异大。山地集中分布在京津冀城市群西北地区，以太行山脉、燕山等为主，土地利用类型主要为林地、草地；平原集中分布在京津冀城市群东南部，以华北平原为主，土地利用类型以耕地、人工表面为主。气候为温带大陆性气候，冬季严寒，夏季炎热，气温年较差大，故植被类型以温带落叶阔叶林为主，受气候影响较大。京津冀城市群作为中国的“首都经济圈”，其快速的城镇化及高强度的城市扩张对京津冀的土地利用类型改变强度大，使原本脆弱的生态环境更加脆弱，区域内部的生态服务功能明显下降，热岛、雾霾等生态环境问题日益严重，已严重威胁到京津冀城市群一体化的可持续发展（马程等，2013）。

2 数据与研究方法

2.1 数据来源及预处理

本研究在地理空间数据云网站（http://www.gscloud.cn/）中下载京津城市群2005年、2010年的Landsat ETM+以及2015年Landsat8 OLI遥感影像，影像均无云层干扰。根据“中科院土地利用覆盖分类体系”，利用目视解译的方法，将研究区分类为林地、草地、耕地、水域、人工表面、未利用地，空间分辨率为30 m×30 m。为了保证解译精度，结合同期更高分辨率的Google Earth历史影像及野外调查进行精度验证，3期影像的解译精度均在93%以上。

在美国Earth Date search网站（https://search.earthdata.nasa.gov/）中，下载 2005年、2010年、2015年的MOD17A3数据，空间分辨率为1 km×1 km，预处理后得到京津冀城市群的NPP年均数据；同时，下载以上3年7月内影像质量较好的MODIS植被指数16天合成产品（MOD13A2数据），空间分辨率为1 km×1 km，预处理后得到京津冀城市群的NDVI数据。

降水数据来源于中国国家级地面气象站基本气象要素日值数据集（V3.0），从中选择2005年、2010年、2015年京津冀150个气象站点的年平均降雨量数据，并利用克里金插值方法得到京津冀城市群的降水数据，空间分辨率为100 m×100 m。

土壤保持模拟量数据是基于降水数据、NDVI数据、土壤数据等其他数据，利用刘军会等（2018）研究中应用的RULSE模型计算得到，空间分辨率为1 km×1 km，其中土壤质地数据来源于国家地球系统科学数据共享服务平台土壤科学数据中心（http://soil.geodata.cn/），土壤深度数据来源于寒区旱区科学数据中心（http://westdc.westgis.ac.cn/），DEM数据来源于地理空间数据云。

为了保证数据图层间的运算更加科学，以及减小像元混合现象对研究结果带来的影响，本研究把所有数据全部依据最邻近法重采样为空间分辨率为100 m×100 m的数据，并利用ArcGIS 10.2软件中的Nibble工具使所有数据的边界重合。

2.2 单位面积生态系统服务价值的当量确定

本研究基于谢高地等（2008）在 2008年提出的“中国生态系统单位面积生态服务价值当量”（表1），结合京津冀城市群土地利用一级分类结果，实现了对该地区生态系统服务价值的评估。该表中将1 hm2全国平均产量的农田每年自然粮食产量的经济价值定义为 1，其他生态系统服务价值当量因子是指该生态系统产生的生态服务相对于农田食物生产服务的贡献大小（谢高地等，2003），在此不考虑人工表面的生态系统服务价值（姜翠红等，2016；向悟生等，2009）。

2.3 单位面积生态系统服务价值系数计算

参考《2006中国粮食年鉴》、《2011中国粮食年鉴》、《2016中国粮食年鉴》中的北京市、天津市和河北省3年中各粮食作物的播种面积、粮食单产以及对应的全国平均价格，利用式（1）（赵江等，2016）确定京津冀城市群单位面积农田食物生产服务经济价值：

式中，Eat为单位面积农田生态系统在第t年提供的食物生产服务功能的经济价值（yuan·hm-2）；i为研究区粮食作物的种类，京津冀城市群的主要作物有稻谷、小麦、玉米和大豆；n为研究区的粮食作物种类和，在此n=4；Sit为第i种粮食作物研究区在第t年的种植面积（hm2）；St为n种粮食作物在第 t年总的种植面积（hm2）；pit为第 i种粮食作物在第 t年的全国平均价格（yuan·kg-2）；qit为第 i种粮食作物在第t年的单产（kg·hm-2）。1/7是指在没有人力投入的自然生态系统提供的经济价值是现有单位面积农田提供的食物生产服务经济价值的1/7。

表1 生态系统单位面积生态服务价值当量Table 1 Equivalent value per unit area of ecological service

根据式（1）可知，单位面积农田食物生产服务经济价值与粮食单产存在正相关，故本研究基于京津冀城市群县域水平的粮食单产，构建单位面积农田食物生产服务经济价值的调节因子，实现对单位面积农田食物生产服务经济价值的调节。计算公式如下：

式中，Bat为单位面积农田食物生产服务经济价值在第t年的调节因子，xit为第i个县在第t年的粮食单产，为n个县在第t年粮食单产的平均值，此时n=166。由于在北京市、天津市等城市中心，不存在或存在的农田较小，没有粮食单产数据，在此区域的Bat统一设定为Bat的最小值0.5。

利用式（3）（赵江等，2016）确定京津冀城市群不同县域内其他生态系统单位面积生态服务价值，从而确定京津冀城市群不同县域的单位面积生态系统服务价值系数，计算公式如下：

式中，Eijt为单位面积在第t年第i种生态系统对应的第 j种生态系统服务价值（yuan·hm-2）；fij为第i种生态系统对应的第j种生态系统服务价值当量因子；i为6种生态系统类型；j为生态系统服务功能类型，包括食物生产、原料生产、气体调节、气候调节、水文调节、废物处理、土壤保持、维持生物多样性和美学景观共9种。

2.4 单位面积生态系统服务价值时空动态调节因子

由于生态系统在不同单位区域、不同时间段的内部结构与外部形态不一致，因而其所具有的生态服务功能及其价值量存在时空异质性（李晓赛等，2015），而这种时空异质性与各区域的生物量、降水量、植被覆盖度、地形坡度、土壤性质等因素息息相关（李士美等，2010）。谢高地等（2015）通过构建NPP、降水、土壤保持时空动态调节因子探索了生态系统服务价值在中国省（直辖市）域尺度上的空间异质性，本研究基于格网的尺度探究京津冀城市群生态系统服务价值的空间异质性（唐秀美等，2017），若引用谢高地的调节因子可能会过分的夸大或降低某些格网内的生态系统服务价值，故本部分基于归一化的思想修正了谢高地的时空动态调节因子。此外，选取的格网尺度过大会影响要素异质性的充分体现，格网过小会导致数据量过大，故本研究基于研究数据的最小分辨率（1 km×1 km），选取 3 km×3 km 的网格，将研究区划分为24767个相互独立的地理空间单元，实现了对京津冀城市群单位面积生态系统服务价值的时空调节。

NPP时空动态调节因子公式：

式中，Pkt指第 k个地理空间单元在第 t年的NPP时空动态调节因子，Bkt指第k个地理空间单元在第 t年内的 NPP（t·hm-2），Bmt指京津冀城市群所有地理单元在第t年内的NPP平均值。该因子用于调节区域内食物生产、原材料生产、气体调节、气候调节、废物处理、维持生物多样性和提供美学景观的服务功能价值。

降水时空动态调节因子公式：

式中，Rkt指第k个地理空间单元在第t年的降水空间动态调节因子，Wkt指第k个地理空间单元在第 t年内年降水量（t·hm-2），Wmt指京津冀城市群所有地理空间单元在第t年的平均降水量。该因子用于调节区域内水文调节服务功能价值。

土壤保持时空动态调节因子公式：

式中，Skt指第k个地理空间单元在第t年的土壤保持空间调节因子，Lkt指第k个地理空间单元在第 t年的土壤保持模拟量（t·hm-2），Lmt指京津冀城市群所有地理空间单元在第t年的平均土壤保持模拟量。该空间动态调节因子用于调节区域内土壤保持服务功能价值，其中土壤保持模拟量是利用修正通用的水土流失方程RULSE模型（刘军会等，2018）计算得到，公式为：

式中，Ac为土壤保持模拟量，Ap、Ar依次为潜在土壤侵蚀量与实际土壤侵蚀量，R、K依次为降雨侵蚀力因子与土壤可蚀性因子，R值、K值参考王万忠等（1996）研究，分别基于降水资料计算得到与基于土壤质地赋予不同的 K值；L、S依次为坡长、坡度，基于DEM数据提取后得到；C为植被覆盖度因子，在此引用NDVI值。

2.5 生态系统服务价值计算

为反映京津冀城市群生态系统服务价值的异质性，本研究基于式（8）利用空间调节因子对单位面积生态系统服务价值系数进行调节，并结合 3 km×3 km的地理空间单元，利用式（9）统计分析每一地理空间单元的生态系统服务价值。

式中，Ekt为第k个单位面积在第t年的生态系统服务价值（yuan·hm-2）。

空间地理单元内生态系统服务价值计算公式：

式中，ESVt表示第 t个地理空间单元的生态系统服务价值总量（元）；n表示第 t个地理空间单元内包含的单位面积总数；1代表1 hm2，为1单位面积。

2.6 生态系统服务价值的空间异质性表达

热点分析方法能够识别具有统计显著性的高值（热点）和低值（冷点）的空间聚类区域并把空间聚类的位置进行可视化表达。本研究利用热点分析方法对京津冀城市群生态系统服务价值进行分析，得到生态系统服务价值低值及高值发生聚类的位置及值，值越趋近于0时，说明生态系统服务价值没有聚集情况，呈随机分布；的绝对值越大，说明生态系统服务价值低值及高值聚集的程度就越高。当值为正时，该聚集区域为“热点区域”，是高值聚集的区域；当值为负时，该聚集区域为“负热点区域”，即冷点区域，是低值聚集的区域。公式（王蓓等，2016）为：

式中：n表示地理空间单元数；wij为第 i与 j地理空间单元的空间权重系数，反映第i与j区域在空间上的关系，本研究中 wij由反距离方法自动生成空间权重文件；xi是第i单元某现象的发生数，本研究中为地理空间单元内生态系统服务价值总量；为所有地理空间单元生态系统服务价值总量的平均数。

2.7 生态系统服务价值的时空动态演变

景观指数可以高度浓缩景观格局信息，能够反映景观结构组成和空间配置等方面特征。在基于等量分类的方法对京津冀城市群生态系统服务价值从低到高依次分为5个等级（I，II，III，IV，V，对应的栅格值依次为1，2，3，4，5），选取类型层次上的景观面积（Patch Density，PD）、蔓延度指数（Aggregation Index，AI）表达不同等级的生态系统服务价值在空间上的团聚程度或延展趋势（彭建等，2006），并结合动态度、变化轨迹分析方法反映不同等级的生态系统服务价值在全域、局部的变化情况。

动态度（刘盛和等，2002）模型的表达式为：

式中：K为研究时段内某服务价值等级的动态度；Ua和 Ub分别为研究期初及研究期末生态服务价值等级的斑块密度或聚合度指数；t为研究时段长，当t的时段设定为年时，K的值就是该研究区某服务价值等级的年变化率，本研究中t=5。

变化轨迹分析方法（黄勇等，2015）是利用轨迹代码记录空间内某一属性的动态变化方法。首先利用式（12）计算两个时间段内生态系统服务价值的变化情况，其次利用式（13）生成变化轨迹代码：

3 分析与讨论

3.1 调节前后单位面积生态系统服务价值的对比分析

根据表2可知，从2005－2015年，在未被空间调节因子调节之前，单位面积生态系统服务价值的最大值依次为38244.109、61315.926、80092.641 yuan·hm-2，平均值依次为 12925.517、20572.184、26381.263 yuan·hm-2；在利用空间调节因子调整后，单位面积生态系统服务价值的最大值依次为62053.742、102983.422、126188.164 yuan·hm-2，平均值依次为 13081.535、20407.107、27056.767 yuan·hm-2。发现单位面积生态系统服务价值处于不断增加状态，造成这种现象的原因不是生态系统的服务功能提高了，而是随着科学技术的进步以及社会经济的发展，粮食单产不断增加，粮食价格逐渐增高的原因。

表2 单位面积生态服务价值统计表Table 2 Statistical table of ecological service value per unit area

调节前后，单位面积生态服务价值的平均值比较接近，只是调整后的最大值明显大于调整前的最大值，研究结果与李晓赛等（2015）的研究结果相似。除此之外，调节后的结果也表明了同一种土地利用类型的单位面积生态系统服务价值表现出了明显的空间差异性（图2），例如：林地在承德、秦皇岛、唐山市的生态服务价值较高，在张家口市的生态服务价值较低；内陆湖泊的单位面生态服务价值明显高于环渤海湿地的单位面积生态服务价值。这种空间差异与区域降水、净初级生产力、土壤保持能力等因素息息相关（谢高地等，2015），更符合生态系统受环境和生物条件影响而表现出的复杂性特点（李晓赛等，2015）。

3.2 单位面积生态服务价值空间特征分析

图3表明，京津冀城市群单位面积生态系统服务价值在2005年、2010年、2015年的空间分布特征相似，热点区域集中分布在京津冀城市群北部、太行山及燕山等山地区域、环渤海湿地、白洋淀及衡水湖等内陆湖泊地区，冷点区域集中分布在平原地区及张家口市等高原区域。由表1可知，林地的各项生态服务功能均较强，故林地的综合生态系统服务功能相对较强（赵江等，2016），在京津冀北部以及太行山、燕山附近，景观类型以林地为主，植被覆盖度高，人为干扰较少，故其生态系统服务价值较高，生态系统服务价值的等级多为IV、V类；另外，由于河湖的水文调节功能及废物处理功能最强，且该两项生态服务功能远大于其他土地利用类型，故其生态系统服务价值较高，服务价值等级多为第V类（王蓓等，2016）；京津冀城市群的东南地区以及西北地区，生态系统服务价值等级多为第II、III类，单位面积生态服务价值较低，受土地利用类型影响较大，其土地利用类型以耕地为主，水域、林地等生态用地的面积较少，耕地的各类生态系统服务功能与林地、水域比较相对较弱（Lin et al.，2018）；由于本研究假定人工表面的生态系统服务价值为0（姜翠红等，2016），所以生态系统服务价值的第 I类区域的土地利用类型多为人工表面，故生态系统服务价值最低，在各城市建成区表现出明显的低值聚集现象。

3.3 生态系统服务价值动态变化研究

3.3.1 景观格局指数动态度分析

由表3可知，从2005－2010年，京津冀城市群生态系统服务价值的第I、II、V类的PD动态度依次为 0.020、0042、0.029，AI的动态度依次为-0.007、-0.046、-0.025。这 3种等级的斑块密度均增加了，聚合度均降低了，与图4反映的信息一致，且第二等级的变化程度更剧烈，侧面反应了这段时期内这3种等级的破碎度增加了，且第二等级的破碎状况更严重。第III、IV类的PD动态度依次为-0.015、-0.014，AI的动态度依次为 0.027、0.022，与第I、II、V类的正负性相反，且两者之间的变化程度接近，说明这两个等级的破碎度降低了，聚合度升高了。从 2010－2015年，生态系统服务价值的第I、II、V类的 PD动态度依次为-0.002、-0.003、-0.002，AI的动态度依次为0.001、0.002、0.002；第 III类的 PD动态度为 0.002，AI的动态度为-0.003；第IV类的PD动态度与AI动态度均为0。

图2 单位面积生态系统服务价值Fig. 2 Ecological service value per unit area

综合以上分析，后5年内PD、AI的动态度与前5年相反，且动态度的绝对值明显较小。在前5年内，城市化发展强度高，城市扩张、人口向城市聚集的现象明显，侵占了林地、草地、湿地等生态用地，改变了景观的内部结构，使得景观破碎化严重。同时人们的环保意识不足，为了追求高的GDP，重点发展工业，三废排放严重，使得雾霾、高温、气候异常、河流污染等环境问题日益严重（严珅等，2018），这些原因使得水文调节、气候调节、气体调节、维持生物多样性等生态服务功能明显下降，从而造成该时段内生态系统服务价值破碎度增加较多，蔓延度降低幅度大（Su et al.，2012）；而在后5年内，随着人们的环保意识加强，第三产业比例一直处于上升状态，环境污染得到了一定的控制，另外国家于 2014年强调的京津冀一体化政策以及京津冀协同发展重大国家战略的实施更是强调了区域的协同发展，强调了人口经济资源环境相协调的重要性（苑清敏等，2017），对各区的盲目城市扩张现象进行了限制，降低了对生态服务功能的干扰程度（吕金霞等，2018）。

3.3.2 变化轨迹分析

图3 生态系统服务价值等级及热点分析结果Fig. 3 Ecosystem service value level and hotspot analysis results

图4 服务价值等级的PD、AI指数统计Fig. 4 PD, AI index statistics of service value level

由图5可知，从2005－2015年，第I类生态服务价值的变化区域集中分布在天津市、沧州市、衡水市等平原区域以及张家口市。在天津市周边，变化情况多为前5年不变，后5年升高；在沧州市、衡水市，变化情况多为前5年等级提高，后5年等级不变。在平原地区，第II类生态服务价值的变化区域分布较多，先升后不变或下降的区域集中分布在衡水市、邢台市，先降后不变或上升的区域集中分布在沧州市。前两类的生态服务价值较低，且集中分布在以耕地为主要景观的平原地区以及以草地、耕地为主要景观的高原地区，这些区域人类活动强度大，景观受人为干扰程度严重，景观变化往往会削弱景观变化较为密集的城乡结合部的生态系统服务功能和景观价值（Wu et al.，2015），由于前5年内的城市活动比后5年内的城市活动强，故有此现象。第 III类生态服务价值的持续下降区域或先下降后不变的区域集中分布在平原与山地丘陵的交界处以及邯郸市，先升后不变或下降的区域多分布于张家口市。因为该区域是生态服务功能强弱的过渡区，耕地、林地、草地面积均较多，毁林造田、毁山造田等人类活动破坏了原有的植被，降低了植被覆盖度，增加了水土流失的风险指数，导致该区域内的水土保持、水文调节等服务功能下降，从而降低该区域内的生态服务功能（严珅等，2018）。由于生态服务价值的第IV、V类区域多集中分布在燕山、太行山等山地丘陵地区以及环渤海区域，故其变化区域也集中分布于此，第 III类生态服务价值的变化以先升后降或不变为主，多分布于山地丘陵的西部；第 III类生态服务价值的变化以持续下降或先下降后不变为主，集中分布于承德市东部。山地丘陵西部的平均海拔较高，人为干预少，而承德市东部的海拔较低，是人类活动相对频繁的山地区域。

4 结论

本研究基于3 km×3 km的网格，基于格网尺度改进了单位面积生态系统服务价值的调节因子，并结合热点分析方法、景观格局指数、动态度模型、变化轨迹分析方法分析了京津冀城市群生态系统服务价值的空间异质性及其动态变化特征，主要有以下结论：

表3 服务价值等级的CA、AI指数动态度Table 3 PD, AI index dynamics of service value level

（1）京津冀城市群单位面积生态系统服务价值的热点区域集中分布在京津冀城市群北部、太行山及燕山等山地区域、环渤海湿地、白洋淀及衡水湖等内陆湖泊地区，服务价值等级以 V类为主，IV类为辅；冷点区域集中分布在平原地区及张家口市等高原区域（图3），服务价值等级以I、II类为主，III类为辅。这种差异与京津冀城市群的景观类型、植被覆盖度等因子相关密切。

（2）从 2005－2010年，京津冀城市群生态系统服务价值的第I、II、V类的PD动态度均为正值，AI动态度均为负值，说明斑块密度均增加了，聚合度均降低了；第III、IV类的PD动态度均为负值，AI的动态度均为正值，说明斑块密度降低了，聚合度升高了。从 2010－2015年，京津冀城市群生态系统服务价值的第I、II、III、V类的PD、AI动态度的正负性与前5年相反，第IV类的PD、AI动态度均为0值，说明该类的斑块密度与聚集度变化不明显。

图5 服务价值等级的变化轨迹分析结果Fig. 5 Analysis results of the change trajectory of service value level

（3）在张家口市、衡水市、邢台市，生态系统服务价值的等级多为先升后不变的区域，沧州市生态系统服务价值的等级既存在先升后不变的区域，又存在先降后不变的区域；第 III类生态服务价值的持续下降区域或先下降后不变的区域集中分布在平原与山地丘陵的交界处以及邯郸市；第 III类的变化以先升后降或不变为主，多分布于山地丘陵的西部；第 III类的变化以持续下降或先下降后不变为主，集中分布于承德市东部。

本研究选取 3 km×3 km的网格来缩小研究尺度，并改进了谢高地等（2015）基于市域尺度上提出的生态系统服务价值的空间调节因子，改进的模型降低了基于网格统计信息可能会存在的异常值对研究结果的影响程度，并从县域的尺度上修正了单位面积农田食物生产服务价值。其中调节因子是基于格网尺度构建的，其受网格尺度的影响，故后续可利用不同尺度的网格做对比分析，以求可以表达研究结果的更优尺度，从而降低网格尺度大小对生态系统服务价值空间差异准确表达的影响程度。同时，该方法不仅适用于城市群的尺度，而且通过改变网格的大小适用于任何区域的生态系统服务价值的空间异质性探索，包括市域、县域尺度。

致谢：感谢京津冀国家重点研发计划课题（2016YFC0503001）提供了这次写作的机会，感谢地理空间数据云网站(http://www.gscloud.cn/)以及美国Earth Date search网站（https://search.earthdata.nasa.gov/）提供的遥感影像资料，感谢中国国家级地面气象站提供的气象数据资料，感谢国家地球系统科学数据共享服务平台土壤科学数据中心提供的土壤质地数据，感谢寒区旱区科学数据中心提供的土壤深度数据。