基于MCR与Markov-CA模型的城镇增长边界划定研究——以江苏省新沂市为例

■ 李全宝

（江苏省新沂市国土资源局，新沂 221400）

0 引言

当前中国正处于新型城镇化快速发展时期，人口将持续向城镇集中，经济增长与产业发展需要更多空间支持，城镇边界也将不断扩张。据统计，2009—2016年，我国城镇化率从48.34%提高到57.35%，年均提高1.29%，城市建设用地面积从38726.92 km2增加到52761.30km2，年均增长4.52%，土地城镇化速度远高于人口城镇化速度[1]。自改革开放以来，我国城镇扩张蔓延粗放，“摊大饼”式的发展带来了一系列城市问题，如生态和粮食安全、环境污染、人口过度拥挤、占用耕地等，严重制约社会经济可持续发展[2]。因此，如何通过划定用地增长边界来严格控制城镇无序蔓延及引导合理增长，已成为城镇用地增长管控的核心任务。

当前我国城镇空间增长已步入优化转型时期，开始表现为从追求规模到追求质量[3-4]。 第三轮《全国土地利用总体规划纲要（2006—2020年）》提出研究“中心城区空间增长边界”的要求[5]。2016年中共中央办公厅、国务院办公厅印发的《省级空间规划试点方案》中明确提出以主体功能区规划为基础，划定城镇、农业、生态空间以及生态保护红线、永久基本农田、城镇开发边界[6]。城镇用地增长边界划定已成为当前学术界研究的热点问题，但目前该项任务仍缺乏统一的划定方法。目前在研究和探讨划定城镇用地增长边界（UGB）的方法中，大体上可以分为如下三类：一是从“反规划”角度，通过明确禁止建设侵入地区，以框定城镇增长边界[7]。丛佃敏等基于“反规划”的理念，以天水市为例，选取了DEM、水土流失、地表水源保护区、植被覆盖率、基本农田等17个生态因子，运用构建综合生态安全格局，划定不适宜城市扩张的区域，从而划定城市扩展的边界，促进城市与生态环境的协调发展[8]。二是基于GIS技术的土地适宜性评价划定方法，通常采用地理信息系统（GIS）技术，对自然资源、交通等图层叠加分区划定UGB。李旭锋以哈尔滨市为例[9]，祝仲文等以防城港市中心城区为例[10]，均运用此方法划定UGB，并提出了空间管制对策。三是构建数学模型划定UGB。主要是元胞自动机模型（CA）、约束性CA、基于BP人工神经网络(ANN)模型等[11-18]。如龙瀛等运用地理信息系统软件与约束性CA模拟了北京城市空间的增长，划定了中心城、新城和乡镇三个层次的UGB[18]。

本研究以江苏省新沂市作为研究区域，在借鉴前人研究的基础上，将最小累积阻力模型（Minimal Cumulative Resistance，MCR）与Markov-CA模型相结合，既弥补了以往“千层饼”评价对景观水平过程的忽略，又能弥补MCR无法进行时间模拟预测的缺陷，同时考虑了绝对阻力对城镇扩张的刚性约束，从而使划定城镇用地增长边界更符合实际，以期为研究区城镇用地增长边界划定提供决策参考。

1 研究方法

本文通过结合阻力分析与动力分析两方面机制，耦合MCR模型与Markov-CA模型划定城镇增长边界。一方面通过MCR模型对城镇扩张的阻力面进行分析，重在探究城镇“基底”的发展潜质，从而评估城镇扩展的可行性及限制性，进而获得研究区的刚性增长边界；另一方面，以Markov-CA模型作为动力模型模拟预测城镇的扩张边界，该模型弥补了单纯的MCR模型无法进行时间模拟预测的缺陷。MCR与Markov-CA模型的结合，旨在从城镇扩展的外在阻力约束与内在动力机制两个方面更好地模拟真实城镇扩展，以此划定城镇增长边界。

1.1 MCR模型的构建方法

1.1.1“源”的确定

源（Source）是事物或事件向外扩散的起点和基地，具有内部同质性和向四周扩张或吸引的能力，扩张或吸引能力的大小与“源”的性质和四周传播媒体的性质有关。选择现状城镇用地作为扩张“源”，现状城镇用地作为人口集聚地和经济活动中心，对周围的土地具有较强的辐射能力，辐射能力随距城镇距离的增加而减小，即城镇扩张的范围越大，它所要克服的阻力越大。

1.1.2 基面阻力及系数

因基面属性的不同，城镇扩张所面临的阻力也不同。阻力系数（Resistance Coefficient）表示城镇用地在扩张过程中通过某一点的难易程度。基面阻力分为绝对阻力和相对阻力两类。绝对阻力即城市发展不可突破的生态控制底线，是城镇增长规模的上限，包括自然保护区、风景名胜区、洪水调蓄区、生态公益林等，这些区域作为城镇扩张所遇到的“壁垒”，一般受政策严格调控，一旦确定则对城镇用地扩张阻力非常大。相对阻力可以从地形坡度、植被覆盖、土地利用现状、区位条件、建设密度、交通可达性等6个方面建立。通常情况下，坡度越大、植被覆盖度越高、区位条件越落后、建设密度越低、交通越不发达，城镇扩张的阻力越大，其阻力系数越高；不同土地利用现状类型对城镇扩张的生态阻力差异较大。

1.1.3 阻力面模型

MCR最早应用于物种扩散过程研究，后来被广泛用于生态空间、土地利用格局等研究。该模型主要考虑三个方面的因素，即源、距离和成本面（或称阻力面，阻力越大成本越高）。基本公式[19-20]为

其中，f是一个未知的负函数，反映扩张过程与最小累积阻力值之间的负相关关系；Dij是地块i到源j的空间距离；Ri是地块i对源j扩张的阻力值；min∑表示从地块i到源j的最小累积阻力值。因此，阻力面反映了源i在空间中扩张的潜在可能性及趋势。

1.1.4 模型理论分析

根据相对阻力评价指标体系，计算生成各个因子的阻力值评价图，并对其进行综合加权叠加分析，得到垂直过程的各个因子对于城市扩张的阻力大小，运用ArcGIS中的Cost-distance模块最终得到城镇扩张在水平运动过程中的最小累积阻力。

1.2 Markov-CA模型

城镇扩张的阻力分析是主要考虑城镇用地外部条件对城镇扩张的影响，但城镇用地本身也有其内在的发展规律，即自身存在扩张的动力，因此选用Markov-CA模型对城镇扩张进行动力分析。

Markov理论是指随机过程中，在有限的时序t1、t2…tn中，任意时刻tn的状态an只与其前一时刻tn-1的状态an-1有关，这种无后效性的过程即Markov过程。在城市扩展研究中，各类型土地的相互转化可视为Markov过程，它只与其前一时刻的土地利用类型相关。各类土地之间相互转化的面积或比例即为状态转移概率，公式为[21-24]：

式（2）中，Pij为状态转移矩阵，S（t）、S（t+1）分别表示t、t+1时刻的土地利用状态，状态转移概率矩阵预测了事件发生的状态及其发展变化的趋势。

元胞自动机CA（Cellular Automata)具有强大的空间运算能力，常用于自组织系统演变过程的研究。它是一种时间、空间、状态都离散，空间相互作用和时间因果关系都为局部的网格动力学模型。CA有四个基本要素：元胞、状态、邻域和转换规则[25-26]。公式为：

Markov-CA模型综合了Markov长期总量预测的优势和CA模拟复杂空间演化的能力，提高了模型模拟的精度。

2 实证分析

新沂市位于江苏省北部，是徐州都市圈的重要节点，也是沿东陇海线产业带的中心节点城市，具有得天独厚的交通区位优势。2015年，新沂实现地区生产总值(GDP)507.63亿元，其中，三次产业比例为12.2∶41.8∶46.0，在中国百强县位居第64位，竞争力等级为A+级。新沂市区所在城镇与其余城镇之间的过渡不自然，城镇体系分布较松散，城镇规模等级结构尚不完善。随着未来新沂市城镇化进程加快，其城镇内部的空间结构势必发生重要变化，不同地类的使用在其过程中遇到的冲突与矛盾也会愈演愈烈。因此，划定该区域城镇用地增长边界，对引导城镇用地合理有序扩张显得非常重要与紧迫。

2.1 城镇增长阻力分析结果

绝对阻力的选择是依据新沂市生态红线划定的内容，包括马陵山风景名胜区、马陵山生态公益林；沂河、沐河、新沂河洪水调蓄区、新沂骆马湖湿地自然保护区等，总面积167.84km2，占全市面积的10.54%，该部分区域即为城市扩张的禁止建设区。相对阻力的分析通过MCR模型实现，现状城镇用地作为城镇扩张源，总面积55.95km2，新沂市下辖16个镇，城镇的分布如图1所示。三个重点城镇为新安镇、北沟镇和唐店镇，地处新沂市中北部，镇区城镇用地面积为31.06km2，占全市城镇用地的56%，这里也是新沂市的中心城区。

图1 新沂市现状城镇用地分布

从地形坡度、植被覆盖、邻域、区位条件、建设密度和交通可达性六个方面建立评价指标体系（表1）。参考相关标准及前人的研究经验，本文将阻力系数分为5级，值越大，阻力越大。各评价因子权重的确定采用层次分析法，判断矩阵的一致性比例为0.097〈0.1，通过一致性检验，最终整理出的因子权重及分值。

运用自然断点法[26-27]，对阻力值划分为5类，并对各个区间进行面积统计，结合2014年土地利用类型图分析各个区间适宜的或者主要地类，见表2。

最小阻力值区间1.29～2.49所占的区域面积为58.97km2，与现状城镇用地面积相差3.02km2，分布上基本相符，表明在已有的基础上进行城镇建设所遇到的阻力较小；阻力值区间2的面积为240.59km2，主要的地类为村庄，并且在地理位置上，大多与城镇用地距离较近，或在交通线路附近，体现了距离城镇、道路近，城镇扩张阻力小，更易于村庄向城镇发展的趋势；阻力值区间3与阻力值区间2在分布上具有一定的相关性，区间3在区间2的附近，并且大致沿着交通道路分布，体现了道路对城市发展的影响；阻力值区间4的阻力较大，主要的地类是耕地；阻力值区间5的阻力最大，除了绝对阻力（河湖水库等）外还包括其附近一定的区域。

表2 阻力值情况

在城镇扩张综合阻力评价结果的基础上，运用ArcGIS中的Cost-distance模块，以现状城镇用地为扩张“源”，计算每个单元到成本阻力面上最近“源”的最小累积成本距离，生成新沂城镇扩张最小累积阻力面。

运用ArcGIS中的Reclassify功能生成城镇扩张最小累积阻力面的分级图，计算各级的面积。参考《新沂市土地利用总体规划（2006—2020）》中对建设用地的空间管制要求，其中的限制建设区是区内土地主导用途为农业生产空间，开展土地整治和基本农田建设的主要区域，区内禁止城镇村建设，严格控制线型基础设施和独立建设项目用地；禁止建设区是区内土地的主导用途为生态与环境保护空间，严格禁止与主导功能不符的各项建设。从对限制建设区和禁止建设区的管制规则可以看出，限制建设区和禁止建设区这两个区域是城镇用地扩张的高阻力区间，这两个区规划到2020年的区域面积分别为1346.14km2、27.55km2，合计1373.69km2。对最小累计阻力分级面积进行逆向运算，即从高阻力值到低阻力值对各级面积进行累积计算，最小累积阻力值从1000～20875区间的面积为1431.24km2，与限制建设区和禁止建设区面积相似，提取最小累积阻力值在1000以内的边界即为城镇扩张的刚性增长边界。

2.2 新沂市城镇扩张模拟结果

通过IDRISI中时序变化模块中的Markov计算得出概率转换规则，以2009年地类作为初始状态，2012年地类作为结束状态，迭代参数为3，比例误差设置为0.15，运行后得出土地利用类型转移概率和面积矩阵。利用Markov-CA模块进行2014年的土地利用模拟预测，为了验证模型的精度，将模拟的和实际的2014年土地利用类型分布图进行精度和准确性评价，计算Kappa指数为0.95，说明模型模拟的效果很好。基于以上转换规则，以2014年经过格式处理的地类分布图作为起始数据，模拟出2020年的土地利用类型分布图。2020年新沂市模拟的城镇用地面积为69.96km2，比2014年实际城镇用地面积增加了14.01km2，提取其中城镇用地的分布，其边界即为城镇弹性增长边界。

2.3 新沂市城市增长边界的综合划定

通过对新沂市绝对阻力分析，提取城镇扩张的刚性边界，边界范围内即为城镇禁止建设区域；在相对阻力MCR模型模拟过程，提取最小累积阻力值在0～1000的边界，即为新沂市刚性增长边界；利用IDRISI软件中的Markov-CA模块，模拟2020年的土地利用类型分布，提取其中城镇用地的分布边界，则为新沂市2020年增长边界。中心城区的边界划定体现出新沂市城市总体规划所确定的发展方向即“南拓、东优、西延、北控”的总体思路；刚性边界以外的区域主要作为农业生产空间和生态与环境保护空间，总面积为1431.24km2。新沂市城镇增长边界综合划定如图2所示。

图2 新沂市城镇增长边界

3 结论

本研究将MCR模型和Markov-CA相结合划定UGB，既从城镇生态环境安全角度去划定城

镇不可建设的区域以及明确城镇用地适宜性的程度，同时考虑了城镇用地本身内在的发展规律。此外，本研究对于新沂市未来发展具有一定的政策启示作用。新沂未来要加强寻求资源

开发与环境保护的动态平衡点，合理地配置资源要素。具体而言，依据新沂市社会经济发展水平、资源环境承载力的空间分布及其演化状况，可以考虑在新沂市设立以新安镇为中心的创新提升区，注重拓展其他乡镇为新的经济中心区、旅游观光区文化商贸区等亚区，力求在推进特色产业及中心镇建设的同时，实现县域空间均衡协调可持续发展[28]。

当然，本文探索的方法仍然有待完善与深入研究。例如，本研究重点关注城镇扩张的阻力作用，而对于城镇用地的内在转换规则考虑不够。此外，由于数据收集与技术等的局限，本文在城镇用地扩张综合阻力评价的外部因子选取时，未考虑社会、经济等因素，在进一步研究中应尽可能完善，不断提高该方法的科学性。同时，今后有必要将该模型与其他模型进行对比分析，综合多种模型的优势划定UGB。