基于城市扩展模拟与生态—规模约束的深圳市城市开发边界划定

赵 红 蕊，刘 欣 桐，姚 毅，郑 卓 凡

(清华大学土木工程系交通工程与地球空间信息研究所，北京 100084；清华大学3S中心，北京 100084)

0 引言

改革开放以来，中国城镇化进程实现了跨越式发展，但部分城市无序蔓延和过度扩张导致生态系统功能退化、土地资源紧张、城市空间格局混乱等问题，不利于城市可持续发展[1-3]。为此，我国目前大力倡导完善城市发展体系，科学规划城市发展。城市开发边界(Urban Growth Boundary,UGB)的划定不仅影响城市未来的空间发展方向，对城市发展起到控制和引导作用，更是一种可融入科学的空间规划方法和可持续发展理念的空间治理工具[4,5]。UGB划定方法可分为两类：1)反推法，即借鉴反规划理念，从环境生态安全角度，通过城市用地适宜性评价或生态安全格局构建等，明确城市禁止开发建设的区域，反向划定UGB[6,7]，能有效保护自然生态环境，但对城市内部扩张趋势考虑不足；2)动态模拟法，将城市视为不断变化的复杂系统，综合考虑客观条件下城市内部扩张趋势，基于地理仿真模型划定UGB[8,9]，可满足城市未来发展内在需求，但对生态保护关注不够。因此，将以上两种方法结合，同时考虑生态保护和城市扩展双重需求的综合法成为目前UGB划定的研究热点和发展趋势[10,11]。

元胞自动机(Cellular Automata,CA)是动态模拟法中的一种重要基础模型，采用自下而上的研究思路，具有强大的复杂计算能力和时空动态特征，非常适用于城市复杂系统行为的动态模拟[12,13]；同时，CA模型具有较好的开放性，研究者能根据需要加入其他方法或模型，对标准模型进行优化改进，利用CA研究城市土地利用动态演化过程已成为城市空间动态模拟的主流方向之一[14-16]。在城市扩展中，CA动态模拟的关键在于元胞转换规则的确定，基于城市实际情况，利用多源时空地理数据构建多城市扩展因子驱动与多元胞状态转换作用的城市扩展模拟模型是划定UGB的重要基础[17]。CA模型迭代终止条件的科学性易被研究者忽略，却对模拟结果产生较大影响[18]，基于对历史数据的学习确定迭代终止次数，可能导致城市规模盲目扩张等问题，与城市科学发展理念不符，因此，需要建立在一定规模约束下的城市扩展模型。

在可持续发展背景下，还需充分考虑城市开发与生态环境之间的关系，避免或减轻因城市扩张而引发的生态环境问题。生态约束可分为刚性约束和弹性约束两类：刚性约束是指相关法律法规严格管控的非城市建设区域，具有确定的区域边界，区域内不允许任何开发建设行为，如城市生态保护红线、重点保护区等[7]；弹性约束是指通过生境质量评价或生态敏感性评价等方法对生态环境进行综合评价，得到连续变化的城市建设适宜性结果[3,11]。在CA模拟城市扩展的基础上，考虑两种约束与CA模型结合的方式，构建普适性的UGB划定框架，可为推动城市可持续发展提供支持。鉴于此，本文以深圳市为研究区，兼顾生态保护和城市扩展，以CA模型为基础，探索一种在生态刚性约束、弹性约束以及城市规模约束下，符合可持续发展理念的UGB划定方法，以期为协调深圳市合理开发和生态保护，引导城市健康发展提供参考。

1 研究区概况与数据来源

深圳市是珠江三角洲都市圈核心组成城市(图1)，地处亚热带区域，受东南季风影响，全年高温多雨；地形以平原和低起伏度丘陵为主，梧桐山等少部分山地位于深圳的东部和北部，核心城市建成区位于西南部，南山、罗湖、福田三区为城市中心区域。2020年深圳市地区生产总值仅次于上海和北京，是全国性的经济中心之一[19]。第七次全国人口普查公报显示，深圳市常住人口为1 756.01万人，相比第六次全国人口普查增加713.61万人。深圳市使用经济激励机制促进生态环境保护，于2005年出台了“深圳市生态保护红线”，推动建立多层次、多元素的立体生态安全监测体系。

图1 深圳市位置及2015年土地利用现状Fig.1 Location and land use status of Shenzhen in 2015

研究数据及其来源如表1所示。需要说明的是，人工不透水区域与城市建成区的定义虽然不完全相同，但人工不透水区域是城市建成区的重要标志，已广泛用于城市化进程分析[20]，因此，本文采用1985-2018年全球人工不透水区域(GAIA)[21]数据中深圳市2015年和2018年人工不透水区域分别表征当年实际城市建成区。所有数据统一为WGS_1984_UTM_Zone_49N投影坐标系，栅格数据的空间分辨率均重采样为30 m×30 m。

表1 研究数据Table 1 Description of research data

2 研究方法

2.1 生态环境约束

目前深圳市已完成生态控制红线的划定，本文以深圳市官方发布的2013年深圳市生态保护红线[25]作为深圳市UGB划定的刚性约束条件。生境质量体现了区域环境为其内部生物生存和生产提供的适宜条件，可代表区域生物多样性和生态系统服务水平[26]，与城市可持续发展息息相关，本文采用生境质量评价结果作为UGB划定的弹性约束。InVEST模型[27]的Habitat Quality模块因所需数据便于获取、参数设置简洁灵活等优势[28]，广泛应用于生境质量评价，计算公式如下：

(1)

(2)

式中：Dxj为生境类型j中栅格x的生境退化度；r为不同类型的威胁源，由于城市建设、工业污染等人为干扰对深圳市生境质量的改变起决定性作用[29]，参考文献[30-32]并结合研究区实际情况，选取人口、GDP、城市建成区、城市路网作为生态用地的威胁源；R为威胁源总数；y、Yr分别为威胁源r所占的栅格和栅格总数；Wr为威胁源r的权重；ry为栅格y的威胁值；irxy为栅格y中威胁源r对栅格x的影响；βx为栅格x的可达性；Sjr为生境类型j对威胁源r的敏感程度；Qxj为生境类型j中栅格x的生境质量指数；Hj为生境类型j的生境适宜性；k为半饱和常数，一般取生境退化度最大值的一半；z为归一化常量，一般取值为2.5。

由于建设用地扩展强度与生境质量变化空间分布呈明显的空间负相关关系，城市建设用地迅速扩张会加剧生境质量退化趋势[33,34],因此，在UGB划定过程中，高生境质量的生态用地应得到充分保护，新增城市用地应优先选择中低生境质量地区。

2.2 基于CA的城市扩展模型

2.2.1 城市扩展模型构建 CA是一种离散的动力学模型，其基本计算单元为可以改变状态的元胞，所有元胞组成的域为作用空间，通过定义多种影响元胞状态的规则改变元胞状态，其数学表达式为：

CA=(Ld,S,N,f)

(3)

式中：L为由元胞组成的格网空间，基本单位为一个元胞，本文中代表栅格单元；d为L的维数，取d=2；S为各种元胞状态，不同土地利用类型定义为不同的元胞状态；N为元胞的邻域集合；f为转换规则的映射函数，是CA模型的核心，用以控制元胞状态在城市和非城市间变化。元胞的转化概率Pvalue决定元胞状态的改变，3种元胞状态变化的作用方式分别对应3种概率值：城市扩展驱动因子对元胞状态改变的促进作用Plogistic、元胞本身不同用地类型的非城市用地向城市用地的转换影响Pweight、邻域内城市元胞对非城市元胞转换为城市元胞的促进作用Pneighbor。

Pvalue=Plogistic×Pweight×Pneighbor

(4)

(1)城市扩展驱动因子影响模型。城市扩展是多种驱动因素共同作用的结果，CA模型需要建立各驱动因子与非城市元胞向城市元胞转换间的关系。城市的地理条件决定了城市扩展的空间格局，是城市扩展的基础因素；城市政策直接影响短期内城市的发展方向，是城市扩展的直接因素；城市交通、资源分布、基础设施建设等对城市的发展具有引力作用[35,36]；自然灾害方面，根据深圳市地质灾害防治规划，对于地质灾害易发区主要以重点监控与防治为主，而不是禁止城市建设，故暂不考虑地质灾害对城市扩展的影响。综合上述分析及前人研究[3,14，15]，选取高程、坡度、起伏度、距水域距离、距铁路距离、距公路距离、科研教育点密度、大型购物广场密度、人口密度、距城市副中心距离和距城市主中心距离共11个驱动因子(图2)。其中，坡度和起伏度由DEM数据处理得到，距水域、铁路、公路的距离为各栅格距每个驱动因子的最近欧氏距离，距城市主中心、副中心距离为各栅格距城市主中心和副中心的最近欧氏距离。Logistic回归可以确定自变量与因变量间相关性，将其用于城市扩展模拟领域，以确定城市扩展驱动因子与城市扩展间的关联性。根据Logistic回归，构建各驱动因子对元胞由非城市状态向城市状态转变的概率计算公式(式(5))，进而计算11种城市扩展驱动因子的权重参数(表2)。

表2 城市扩展驱动因子权重Table 2 Weights of urban expansion driving factors

图2 城市扩展驱动因子Fig.2 Driving factors of urban expansion

(5)

式中：A0为常数，设为0；Xi、Ai(i=1,2,3,…)分别为各城市扩展驱动因子值及其影响权重。

(2)不同类型元胞转化模型。不同类型非城市用地向城市用地转变的强度不同，统计相同时间间隔内各类型非城市元胞向城市元胞转换数量占非城市元胞向城市元胞转换总数的比例，计算公式为：

(6)

式中：U→C为不同类型非城市用地向城市用地转换。

(3)邻域作用模型。相邻元胞对元胞状态的变化存在一定的制约或驱动作用，考虑包含24个邻域元胞的7×7型冯诺伊曼邻域，邻域内城市元胞对非城市中心元胞状态变化的影响可定量表达为：

Pneighbor=Sum(Ci,j=city)/24

(7)

式中：Sum(Ci,j=city)为24邻域内城市元胞总数。

2.2.2 随机因素约束模型 由于影响城市扩展的因素众多且复杂，无法用严格的物理或数学模型表达，且无法忽略政治、经济等因素对城市扩展的影响，因此，需要引入随机因素的影响Pthred(随机阈值)以反映城市扩展过程中的不确定性，计算公式为：

Pthred=α(1-δ2)

(8)

式中：δ为0到1之间的随机数；α为控制随机因素对于城市扩展模拟影响大小的参数，本文设定为2。

2.3 UGB划定

以CA模型为基础，将城市规模作为规模约束，并分别以生态保护红线和基于InVEST模型得到的生境质量评价结果作为生态刚性约束和弹性约束，划定UGB。具体流程如下：1)以生态保护红线作为刚性约束，其内部区域在CA模型进行城市扩展模拟过程中不允许发展为城市；2)将生境质量评价结果按照生境质量指数由低到高划分为5个等级，对于第4级和第5级的非城市元胞，不予扩展为城市；3)对上述步骤约束外的元胞按照式(4)计算扩展概率Pvalue，与随机阈值Pthred进行比较，当非城市元胞的Pvalue大于Pthred时，则该元胞转换为城市元胞；4)统计每次扩展后的城市元胞数是否达到规模约束值，若达到，则结束扩展，若此时扩展模拟得到的城市元胞数超过规模约束值，则优先排除最后一次扩展概率小的元胞，使城市元胞数与规模约束值相等，并输出扩展模拟结果；若未达到规模约束值，则重复步骤3)进行迭代；5)UGB为集中连片城市建设用地的连续轮廓，应尽可能多包含大面积连片的城市建设用地[37]，基于此，对步骤4)得到的城市扩展模拟结果进行面积阈值划分、边缘检测等后处理流程，提取UGB。

3 结果分析

3.1 生境质量评价

应用InVEST模型对2015年深圳市生境质量进行评估，并将生境质量划分为0.0～0.2、0.2～0.4、0.4～0.6、0.6～0.8、0.8～1.0共5个等级(图3)，等级越高，表明生境质量越好。由图3可知，生境质量等级高的地区多处于开阔水域、森林、湿地公园、风景名胜区以及周边，这些区域与生态控制红线区域重合度较高，应予以保护，生境质量等级低的地区多处于城市集中建成区，如罗湖区、南山区、福田区；生境质量等级为中低的地区主要位于城市集中建设区的周边区域，如宝安区、龙华区等，为未来城市用地增长的潜在区域。

图3 2015年深圳市生境质量评价结果Fig.3 Evaluation result of habitat quality of Shenzhen in 2015

3.2 2018年城市扩展模拟与精度检验

为区分不同土地利用类型向城市元胞转换的差异，将深圳市土地利用分类数据设置为耕地、林地、草地、水域、城市建成区及其他6种元胞状态。以2015年深圳市土地利用数据(图1)和建成区空间分布图分别定义元胞地类的初始状态和城市与非城市元胞的初始状态。以深圳市生态保护红线[25]作为生态刚性约束条件，以2015年深圳市生境质量评价结果(图3)作为生态弹性约束条件，进行城市扩展模拟。

首先，为验证模型精度，以2018年为模拟目标年份，以2018年深圳市辖区建成区土地面积928 km2[38]作为CA模型迭代终止条件，得到2018年深圳市建成区模拟结果(图4)，将人工不透水面产品(GAIA)[21]中2018年深圳市人工不透水区域视为实际城市建成区与扩展模拟结果进行对比(图5)。由于生态保护红线为刚性约束，新增城市用地均位于生态保护红线外，实际为城市用地而模拟为非城市建设用地(图5中蓝色区域)中43.56%位于生态保护红线内，与生态保护政策不符。经统计(表3)，2018年建设用地模拟的元胞数量为1 066 831个，实际建设用地的元胞数量为1 066 762个，建设用地模拟正确率为97.20%，模拟总精度为97.42%，Kappa系数为0.948，FoM值为0.13。

图4 2018年深圳城市扩展模拟结果Fig.4 Urban expansion simulation result of Shenzhen in 2018

图5 2018年深圳城市扩展模拟结果与实际对比Fig.5 Comparison between simulated urban expansion and actual built-up area of Shenzhen in 2018

表3 2018年深圳城市扩展模拟精度验证Table 3 Accuracy validation of urban expansion simulation result of Shenzhen in 2018

3.3 2025年城市扩展模拟与UGB划定

以2015年为基期，模拟2025年深圳市城市扩展情况。《深圳市国土空间总体规划(2020-2035年)》(草案)[39]中指出,2035年深圳市建设用地规模目标为1 105 km2，结合2015-2019年深圳市实际建成区土地面积[38]进行多项式拟合，估计2025年深圳市建设用地规模目标为1 018.04 km2，以此作为CA模型迭代的终止条件，得到2025年深圳城市扩展模拟结果(图6)。由新增城市像元在深圳市各行政区分布占比(表4)可知，新增城市像元主要集中在宝安区和龙岗区，是未来主要城市扩展区域，罗湖区、南山区、福田区3个中心区现有城市发展水平已经较高，因此新增城市像元占比相对较少，大鹏新区和盐田区主要为生态保护用地，受生态约束，新增城市像元较少。该模拟结果与深圳城市规划文件[40]中指出的未来新增建设用地的主要区域相符。

图6 2025年深圳城市扩展模拟结果Fig.6 Urban expansion simulation result of Shenzhen in 2025

表4 2025年深圳城市扩展模拟新增城市像元各区占比Table 4 Distribution of new pixels in simulated urban expansion area for each district in Shenzhen in 2025

为检验生境质量在城市扩展模拟中的约束效果，在未考虑生境质量约束情况下模拟2025年城市扩展结果，统计生境质量约束前后新增城市像元覆盖生境质量等级情况(表5)。由表5可知，生境质量约束前，新增城市像元较均匀分布在生境质量1-4等级，而考虑生境质量约束后，生境质量4级和5级处无新增城市像元，新增像元主要分布在生境质量前3个等级，使高生境质量用地得到较好保护。以2025年深圳城市扩展模拟结果(图6)为基础，按照集中连片的原则，划定2025年深圳市UGB(图7)。

表5 生境质量约束前后新增城市像元覆盖生境质量等级统计Table 5 Statistics of habitat quality grades covered by new urban pixels with and without habitat quality constraint

图7 2025年深圳市UGB划定Fig.7 Delimitation of urban growth boundary in Shenzhen in 2025

4 讨论与总结

针对城市盲目扩张与无序蔓延等问题，本研究基于CA模型，通过设定元胞状态和转换规则，同时将生态约束分为刚性约束和弹性约束两类，构建多城市扩展因子驱动与多元胞状态转换作用下的深圳城市扩展模拟模型；以2015年为基期，分别使用生态保护红线和生境质量评价结果作为生态刚性约束和生态弹性约束，模拟2018年深圳城市建设区域的结果显示，建设用地的正确率为97.20%，模拟总精度为97.42%，Kappa系数为0.948，FoM值为0.13，表明预测结果达到极佳水平[41]。另外，以深圳市人工不透水区域为精度检验对比区域与实际城市建设区域不完全相符，会对精度检验结果产生一定影响。此外，由于城市实际建设过程中会受政治、经济等多因素影响，可能会对生态因素的重要性考虑不足，而本文在城市扩展模拟过程中，考虑生态刚性约束和弹性约束条件，力求模拟出兼顾城市扩展和生态保护的、科学的城市扩展情形，其结果会与实际城市建设区域存在一定的差异。以2015年为基期，模拟2025年城市扩展的结果显示，新扩展区域主要集中在宝安区和龙岗区，作为未来主要城市发展区域，大鹏新区和盐田区主要为生态保护用地，在扩展过程中要重点考虑生态约束，该模拟结果与深圳市城市规划文件[40]中指出的未来新增建设用地的主要区域一致。考虑生境质量约束模式可有效改变城市扩展模拟结果，使新增城市像元避开生境质量较好的区域，避免或减少城市建设对生态环境的破坏。在城市扩展模拟的基础上，按照集中连片的原则划定2025年深圳市UGB，可为引导深圳市健康发展提供参考。

本研究提出的基于城市扩展模拟与生态—规模约束的UGB划定方法，兼顾城市扩展和生态资源保护，适用于大部分城市的开发边界划定。针对本文中存在的不足，后续研究将对以下问题进一步探索：1)虽然人工不透水区域可作为城市建成区的重要标志，但二者在区域分布上存在一定差距，以2015年不透水区域作为建成区初始状态，其部分区域位于生态保护红线内，导致虽然新增城市元胞均落在保护红线外，但最终城市扩展模拟结果包含一部分生态保护红线内部区域；同时，以2018年不透水区域作为城市扩展模拟精度检验实际参照，对检验结果会产生一定影响，今后需结合夜间灯光等数据提取更高精度的实际城市建成区。2)目前仅以生境质量指数作为生态弹性约束条件，下一步可结合更多生态环境影响因素以及生态系统服务和景观格局指数等多种生态评价方法，确定深圳市生态环境综合评价方案。