随着城市化进程的加快, 我国城镇化进入到“下半程”, 城市在经济高速发展的同时也经历着人口集聚、生态用地被侵占等一系列“城市化”问题^[1-4], 由此导致人类和自然面临越来越多严峻的考验, 这与当前社会主义生态文明观的要求有很大差距。在国家层面上, 将实行最严格的生态环境保护制度, 形成绿色发展理念和生活方式, 推动人与自然和谐发展, 形成现代化建设新格局。

生态风险评价, 作为一种重要的生态环境管理手段^[5], 可以为区域生态管理提供理论依据和技术支持^[6-8]。关于生态风险的评价方法有很多, 但基本上可以划分为两大类:基于风险源汇和基于景观格局^[9]。基于风险源汇的景观生态风险评价方法, 需要识别威胁区域生态系统健康的风险受体和风险源^[10-11]。基于景观格局的景观生态风险评价方法是以土地利用变化为诱因的生态风险评价^[12-13]。在缺乏生态监测资料的历史积累时, 可以利用相对较易获取的土地利用数据描述和评价土地利用方式的不同对生态系统的影响^[14]。土地利用是全球生态环境变化的重要组成部分, 对区域生态安全起着决定性作用。不同的土地利用方式对生态环境所产生的影响具有累积性和区域性^[15-17], 这可以直接反映在生态系统的结构和组成上, 例如景观要素及其功能的变化^[18-19]。很多学者从景观格局角度出发, 评价区域生态风险。然而景观尺度主要集中在流域和小尺度地区^[20-21], 针对城市化发展较快的大尺度城市群区域的研究相对较少。

权重赋值是生态风险评价的关键环节。评价指标的权重赋值可以分为定性和定量两种方法。定性评价根据指标之间的相对大小或相对优劣程度去对指标进行赋值, 然后评价生态风险的相对变化程度^[22-23]。定量评价则通过构建公式或模型, 对生态风险进行评价^[24-31]。在实际研究中, 需要根据具体的评价目标来确定合适的赋权方法。上述的定性或定量生态风险评价方法基本都是针对某一个时间点独立进行, 这已经不能满足探究长时间生态风险的动态转移要求, 因此需要寻求某种方法去探究它的转移方向。曾辉等人利用地统计学方法对整体的生态风险进行空间分析, 辅以半变异函数为工具进行生态风险指数的方向性分析^[32], 李明财等人利用空间重心转移模型对城市扩张的方向进行了分析^[33]。地统计学方法针对单个时间节点的方向性进行分析, 而未考虑多个时间节点的连续方向性的分析。而重心模型则能很好地表明某个变量的转移轨迹, 且经常被运用于经济、消费和人口等研究中^[34-36]。

京津冀地区作为我国面积较大的城市群区域, 在城市群发展过程中土地利用结构的破坏导致了生态环境质量下降, 阻碍了可持续发展的进程。鉴于以上背景, 本文利用生态风险评价方法去识别生态风险空间分布特征, 然后借助重心模型深入探讨综合高风险区域重心的时间演变过程及变化规律。

1 数据与方法 1.1 研究区概况及数据来源

京津冀地区区域总面积21.8万km², 截止2016年底总人口达1.12亿人, GDP总量7.56万亿元。行政区域范围涉及200区(县), 其中北京16个市辖区、天津16个市辖区、河北11个地级市以及168个区(县)。1982年《北京市建设总体规划方案》首次提出“首都圈”概念, 提出了京津冀三地协同发展的理念, 包括北京、天津两个直辖市以及河北的保定、唐山、廊坊、石家庄、秦皇岛、张家口、承德、沧州、邯郸、邢台、衡水等11个地级市。2015年4月30日, 中共中央政治局审议通过《京津冀协同发展规划纲要》, 并且制定了京津冀协同发展总体规划纲要和生态环境保护与交通一体化专项规划, 以促进京津冀地区的可持续发展^[37-38]。同时这些地区先后签订了以《环渤海区域合作框架协议》为代表的区域合作协议和诸多促进环境保护、经济合作、交通合作等方面的合作协议。

研究区气候属于暖温带大陆性季风型气候, 降水分配不均匀, 春秋短促且干旱多风, 夏季高温多雨, 冬季寒冷干燥。

遥感影像数据来源于地理空间数据云(http://www.gscloud.cn/), 其中包括1984年的Landsat 4 MSS数据、1990年的Landsat 5 TM数据、2000年, 2005年和2010年的Landsat ETM数据以及2015年的Landsat 8 OLI数据。选择6月至8月期间没有云或少云(云量低于10%)的影像数据以确保地面特征识别和解译的高精度。

参考国家土地利用现状分类标准(GB/T 21010-2017)和中国科学院地理科学与资源研究所刘纪远等在建设“中国20世纪LUCC时空平台”建立的LUCC(土地利用/土地覆被变化)分类系统, 并结合京津冀地区的实际土地分类情况将土地利用类型归结为7个一级分类:林地、草地、湿地、耕地、人工表面、其他用地和海域。本文的海域指的是2015年的国土最大轮廓线内的每个时间节点海水面积。遥感影像解译主要是在ArcGIS平台上进行的。利用波段组合、图像校正和图像增强技术, 对空间分辨率为80m(MSS)、30m(TM)和15m(ETM + / OLI)的遥感影像进行预处理。在目视解译中使用来自调查网站的多边形样本和来自Google Earth的高分辨率图像。先解译较高分辨率的2015年遥感影像, 然后将解译后的2015年LUCC数据作为基准图, 将不同时间节点的影像进行叠加, 生成各节点的土地利用/覆被数据。依据现场采样点和谷歌地球样本进行精度验证, 解译结果的准确性被证实高于95%。

1.2 生态风险指数 1.2.1 生态风险指数的构建

生态风险指自然或人为因素影响下景观格局与生态过程相互作用可能产生的不利后果^[9]。生态风险指数是指不同景观对生态环境所产生的风险值大小, 用来表征不同景观所导致生态风险的贡献度。本文从两个方面建立土地利用类型与基于景观格局的生态风险之间的联系, 一是利用土地利用生态风险的主观权重, 耦合相应土地利用类型面积比例^{[21, 31, 39]}, 二是利用景观破碎度。景观破碎化主要表现两个方面, 一是景观斑块面积由大变小的过程, 二是随着斑块面积的减少和地理空间隔离程度的增加产生的景观结构破碎化, 可分别用斑块密度指标和蔓延度指标来表示^{[7, 40]}。斑块密度指标表达的是单位面积上的斑块数, 其值越大说明斑块数量增加, 斑块形状趋于不规则, 内部斑块面积缩小, 斑块彼此隔离越来越大, 反映了景观受干扰较强烈, 景观破碎化程度高; 蔓延度是景观连通度的指示指标, 是景观中不同斑块类型的团聚程度或延展趋势, 是对斑块连通水平的衡量, 蔓延度越高, 表明景观中的优势斑块形成了良好的连接, 说明景观中起主导作用的斑块完整性和连通性好, 景观低破碎; 反之则说明景观具有多种要素的散布格局, 斑块间连通性差、隔离度较高、破碎度高。景观生态风险指数的计算公式如下:

(1)

式中, ERI为生态风险指数; n为土地利用类型的数量; A_i为研究区内第i种土地利用类型的总面积(km²); A为研究区内土地利用类型的总面积(km²); W_i为研究区内第i种土地利用类型所反映的生态风险强度参数, 也就是不同土地利用类型生态风险的权重; F斑块密度(个/100hm²); M为蔓延度(%)。

无论选用哪种景观生态风险评价方法, 对于评价指标的权重进行赋值都是景观生态风险评价中至关重要的环节。根据层次分析法^{[21, 39, 41]}的原理和方法确定土地利用生态风险权重。首先对土地利用类型的相对强度指标进行两两比较、判断, 构造出判断矩阵; 其次采用加权法计算指标的标准化权重, 并进行一致性检验。根据一致性检验结果判断所构造的矩阵是否具有满意的一致性。本文采用AHP软件完成如上分析, 一致性比例的计算结果为0.0022(< 0.01), 表明判断矩阵具有满意的一致性, 从而获得6种土地利用类型生态风险的主观权重分别为:林地0.0708、水域0.0961、草地0.1290、耕地0.1691、未利用地0.2238、人工表面0.3112。

1.2.2 基于格网的生态风险小区划分

公式(1)表示景观格局指数一方面由土地利用类型的生态风险贡献度及其面积比构成, 另一方面由斑块密度与景观蔓延度构成, 景观斑块与景观蔓延度由Fragstats 4.2软件中的移动窗口法获得。因此, 为了能够将土地利用类型的生态风险贡献度及其面积比的乘积空间化, 在考虑研究区空间异质性和斑块大小的基础上, 将研究区划分5km×5km的采样网格, 共计9053个生态风险小区(图 1)。在GIS软件中计算生态风险小区的生态风险值, 并将每个网格的值赋予该网格的中心点, 最后将所有中心点值采用克里金插值法进行空间插值得到生态风险的空间分布图^{[12-13, 20-21]}。

1.3 生态风险重心

“重心”的概念源于力学中重心的概念。重心是指在某一区域之内, 如果把该区域看作是一个质地均匀的平面, 那么, 这个平面上便存在一个支点, 使平面保持平衡, 该支点就是该区域的重心点^[39]。在本文中, 生态风险重心是指用生态风险分区中的综合高风险区替代区域质量时形成的重心, 即在生态风险发展过程中, 在综合高风险区域空间上存在一点, 在该点前后左右各个方向上的生态风险力量对比保持相对均衡, 该点就是生态综合高风险重心。

本文采用几何重心法来衡量某一属性的空间分布状况^[41], 计算公式如下:

(2)

式中, X、Y分别表示生态风险重心的行号、列号, 准确表达了重心的地理位置; X_i分别为第i个综合高风险次级区域(即每一个栅格单元)中心的行号, Y_i分别为第i个综合高风险次级区域(即每一个栅格单元)中心的列号, M_i为第i个综合高风险次级区域的属性值。由式(2)可知, 影响重心迁移的因素是各栅格单元综合高风险区的空间坐标和属性值的大小。

2 结果与分析 2.1 生态风险的时空格局分析

生态风险的时空格局及演变特征是区域景观生态环境质量的直接反映, 也代表生态系统受外界胁迫与抗胁迫能力的强弱。根据公式(1)计算获得1984-2015年土地利用生态风险图(图 2), 并提取各类生态风险区域分布面积并计算其所占比例。

如图 2和表 1所示, 京津冀城市群的生态风险越来越严重。其中, 较高风险区和高风险区都呈现出不断扩张趋势, 较高风险区的上升趋势较明显, 高风险区所占比例从0.069%增加至0.262%, 较高风险区所占比例从0.679%增加至1.825%, 高风险区和较高风险区与城市建设用地的吻合度较高, 对应面积也呈逐渐增加的趋势。表明建设区域周边生态风险正在不断地加剧。低风险区主要分布在沧州市、衡水市、北京市北部、保定市西部、邢台市东部和邯郸市东部, 并且所占比例呈现出下降趋势, 从80.744%降低至68.252%, 而较低风险区面积有所增加, 所占比例从14.843%降低至23.370%, 较小程度上加剧了生态风险。中等风险区呈现出上升趋势, 所占比例由3.665%上升至6.292%。景观高风险区呈点状分布于北京北部、天津中部、保定北部、石家庄西部以及张家口西部以及承德市西北部。京津冀城市群区域的较高风险区由呈面状零散分布变化至面状连片分布。其中高风险区和较高风险区多位于各大城市的中心城区附近, 佐证了城市化的快速发展对生态风险影响显著。

2.2 13个地级市综合高风险区比例分析

如表 2所示, 在研究时间内, 各地级市综合高风险区的变化各不相同。其中, 北京市、天津市面积增长百分点分别达到了13.15%和14.85%;唐山市、廊坊市大于零, 分别为2.27%和0.56%;沧州市和衡水市几乎没变化; 邢台市、邯郸市、秦皇岛市、石家庄市、保定市、张家口市和承德市的面积变化百分比小于零, 其中, 邢台市和邯郸市面积降低幅度较小, 承德市综合高风险区的面积降低百分比最大, 达到了10.40%, 其次是张家口市。在1984年, 承德市和张家口市综合高风险区所占比例较高, 分别为27.85%、和20.48%, 而在廊坊市、张家口沧州市和衡水市所占比例均为0.00%;在2015年, 北京综合高风险区所占的比例达到了最高, 达到了22.64%, 其次是承德市和天津市。

2.3 生态综合高风险区的重心转移

由于廊坊市、沧州市和衡水市的综合高风险区所占比例较小, 所以选择其余10个地级市进行综合高风险区重心的计算。根据公式(2)获取1984-1990年、1990-2000年、2000-2005年、2005-2010年和2010-2015年5个阶段内研究区内综合高风险区的重心转移(图 3)。

如图 3所示, 综合高风险区重心转移方向可以大致分为3类:向首都转移、向东部沿海区域变化和按照自身要求转移。向转移发展的地级市有保定市, 其中, 北京既要控制人口, 又要减少城乡建设用地的规模, 保定优越的地理位置势必就将成北京人口和产业疏解的重点方向; 向东部海域发展的地级市有天津市和唐山市, 其中, 天津市是由于滨海新区的快速发展, 天津的城市规划可以被形象的概括为“一个扁担两头挑”, 即以海河为带状中心和发展轴线, 一头是中心城区, 一头是滨海新区。随着滨海新区在2005年开始被写入“十一五”规划并纳入国家发展战略, 成为国家重点支持开发开放的国家级新区。2010年, 滨海新区生产总值达到5030.11亿元, 超过浦东新区, 增速为25.1%;唐山加大了各类新区的开发力度, 南堡经济开发区, 曹妃甸新区、凤凰新城、空港工业园区等新区的建设极大促进了GDP的增长及产业空间与结构调整^[42]。按照自身要求发展的地级市有北京市、邯郸市、石家庄市、邢台市、秦皇岛市、承德市和张家口市, 北京市积极推进京津冀一体化发展、加强生态环境保护、推动产业升级转移, 且由于北京北部山区的阻拦作用, 重心开始向东偏南方向转移; 邯郸在加快邯郸西部矿区污染源治理、深化大气污染防治的同时建立了邯郸经济开发区和冀南新区等发展新区^[43], 一方面邯郸西部的生态环境得到了一定的改善, 另一方面邯郸市东部地区的生态风险有了一定的加剧; 石家庄市的地理位置决定了城市发展方向, 因为其西部地处太行山中段, 东部为滹沱河冲积平原; 邢台市的《邢台市城市总体规划》 (2005-2020)指出:未来城市将主要向东南和西北方向发展, 从2008年的邢台重点发展西北新区到现在的重点发展邢东新区, 实现了中心城区建设方向的战略东移; 秦皇岛既是全国综合交通枢纽城市, 又是旅游城市, 将建设成为一流国际旅游城市为目标, 推进完成新一轮城市总体规划编制、土地利用总体规划和秦皇岛港口总体规划调编, 高标准制定国际旅游城市总体规划, 制定落实建设国际旅游城市的实施意见, 加快形成层次完善、结构合理、体现国际一流要求的旅游规划体系; 承德市的生态风险转移方向分两个阶段, 在1984-1990年向东南方向发展, 1990-2015年向西北方向发展, 一方面是因为城市方向, 另一方面则是水土流失问题所导致的生态风险转移方向发生了变化; 张家口市相比较于1984年, 2015年的风险位置位于其西部, 张家口市由于其草原的破碎化(道路切割)对生态系统产生了负面的影响, 从而导致了生态风险重心的转移。考虑研究区尺度范围较大, 加强城市之间降低生态风险的能力是京津冀区域生态一体化的重点。人类活动所导致的城市扩张是生态风险的主要来源, 所以在土地利用规划和城市规划布局时, 应注意合理的功能分区。

3 结论

本文基于土地利用数据, 利用生态风险评价分析和重心方法, 分析了1984-2015年京津冀地区土地利用时生态风险时空分布特征, 并在此基础上对土地利用综合高风险区的转移过程进行分析, 得到以下结论:

(1) 从京津冀城市群区域格局把握土地利用生态风险, 生态风险总体上呈上升趋势, 并且生态风险格局存在明显的时空分异, 综合高风险区域主要分布在地级市中心城区周边区域, 土地利用生态风险调控的重点是北京市、天津市。

(2) 京津冀城市群各地级市综合高风险的变化各不相同, 但是各个地级生态风险变化有3种情况:上升、下降和基本不变。对于生态风险升高和几乎不变的区域, 要加强城市内部生态的建设, 对于生态风险下降的区域在加强内部建设的同时, 加强与邻近地区的联系。从区域的角度考虑生态风险的影响, 加强各个城市之间的生态联系, 积极促进京津冀生态一体化。

(3) 京津冀城市群综合高风险对生态环境的影响具有一定的区域性, 各个地级市之间的生态环境的变化既有相同之处, 也有不同之处, 主要表现在各地级市综合高风险区的重心转移方向上。根据地级市的生态风险重心, 积极调整土地利用结构, 加强城市之间的生态凝聚力, 发展成为生态环境更加良好的城市群。

本文利用斑块密度、景观蔓延度及土地利用类型面积比及其主观权重对土地利用生态风险定量化评价方法进行了探索与尝试, 还存在一定的不足。土地利用生态风险的影响因素、表征状态及其影响结果都是十分复杂的, 以土地利用数据为基础, 以人类活动所导致的生态风险为主要研究对象的生态风险指数的构建有待进一步研究。