耦合时空高斯模型和斑块元胞自动机的城镇扩张模拟方法

1.本发明属于大数据挖掘与应用技术领域，涉及一种城镇扩张模拟方法，具体涉及一种耦合时空高斯模型和斑块元胞自动机的城镇扩张模拟方法。


背景技术：

2.研究城镇扩张的历史特征、未来趋势和动力机制对于认识人类活动对地球系统的影响，促进全球可持续发展具有重要意义。元胞自动机(cellular automata,ca)模型可模拟和分析城镇用地扩张的动态过程和时空模式，可深化对城市系统演化规律的认知，帮助预判城镇空间未来演化趋势，最终为制定城镇化管理政策提供决策支撑。
3.当前多数基于ca的城市扩张模拟模型(也称为城市ca模型)主要包括三大模块：土地需求评估模块、适宜性评价模块和空间布局模块。需求评估模块用于评估城镇扩张的总需求或年度需求；适宜性评价模块则基于自然地理和社会经济因素应用数理模型、深度学习模型或其他模型方法评价研究区每个空间单元被开发为城市用地的适宜性(即城市建设适宜性)，适宜性一般通过概率的形式表达；空间布局模块基于ca模型的元胞状态、转换规则、邻域定义、空间单元、时空间约束、模型随机性等组件在城市建设适宜性评价的基础上将一定量的城市土地分配于空间之中。
4.但城市ca模型中的三个模块通常是松散耦合的。多数情况下，城市发展的年度土地需求是利用外部模型独立评估的，并直接输入空间布局模块，然后空间布局模块依据城市建设适宜性将年度总需求分配于空间之上。这一过程并没有考虑城市内部不同区域城市土地年度需求及其变化的时空异质性。然而，在现实世界中，城市空间扩张不仅受到底层微观尺度上城市土地利用空间布局的影响，还受到宏观尺度上城市土地需求时空分布的影响。
5.城市ca模型中反映城市土地需求随时间变化的途径通常是基于时间序列分析模型评估城市土地年度需求，然后通过微观层面的空间布局模块按照城市建设适宜性将年度需求落于空间之上。这类模型的代表包括markov-ca、sd-ca等经典模型。一般来说，将年度城市土地需求直接输入空间布局模块的耦合方式考虑了城市土地需求的时间差异性，但忽略了城市土地需求的空间异质性。
6.城市土地需求的空间异质性是一个重要但经常被城市ca模型忽视的方面，它对准确模拟城市扩张过程和模式有重要影响。城市土地需求的空间异质性表征城市内部不同区域城市化建设强度的空间差异性。研究表明靠近城市中心的地区虽然基础设施十分完善，但因可开发空间有限，开发成本高，其年度新增城市土地需求往往较少；而城市边缘地区(城乡结合部)因有大量可开发空间，且基础设施相对完善，往往是城市建设的热点区域，新增城市土地需求较大；而靠近城市外围的区域(包括广大的农村地区)因为基础设施不完善，新增城市用地需求也往往较少。整体来看，年度新增城市土地需求的时空分布呈现出明显的波纹效应，随时间推移逐步由城市中心向外围推移。
7.在模拟城市扩张时，城市土地需求的空间异质性不应被忽视。在利用ca模型模拟
大尺度区域(如全球、全国或区域尺度)城市空间扩张时，城市土地需求的空间异质性通常通过评估每个下辖行政单元的城市土地需求来表达。而在模拟单个城市(或都市区)的空间扩张时，现有ca模型一般通过如下方式考虑城市土地需求的空间异质性：首先，将研究区域分解为更小的行政单位(如县、镇或村)，或划分为规则(或不规则)的小区域；然后，使用外部模型评估各个行政单位或小区域的城市土地需求(可以是总需求也可以是年度需求)；最后以小区域为基本单元基于空间布局模块依据城市建设适宜性将土地需求分配于空间之上。由于各个小区域内的城市土地需求是独立估计的，因此这一类处理方式忽略了城市土地需求时空分布在整个研究区内表现出的空间联系内在规律与宏观模式特征。


技术实现要素：

8.针对现有城镇扩张模拟技术存在的缺陷，本发明提出了一种耦合时空高斯模型和斑块元胞自动机的城镇扩张模拟方法，提出一种时空高斯函数表达城市土地需求在城市内部的时空异质性分布规律，即城市土地需求在时间和空间维度上表现出的从城市中心向城市边缘圈层波动的宏观模式特征。
9.本发明的方法所采用的技术方案是：一种耦合时空高斯模型和斑块元胞自动机的城镇扩张模拟方法，包括以下步骤：
10.步骤1：从城市中心点向外围区域创建等间距同心环，将研究区划分为n个等间距同心环带；
11.步骤2：根据年度t对应的新增城市土地密度高斯函数计算出研究区年度t内各环带i中新增城市土地预期面积a
exp,i,t
及年度t新增城市土地预期总面积a
exp,texp,t
12.步骤3：根据年度t对应的新增城市斑块频率高斯函数抽取一个同心环带i用于执行土地开发行为，即生成城市土地斑块(新增城市斑块频率大小决定该环带被抽取的概率大小)；
13.步骤4：检查年度t内第k次土地开发行为之前同心环带i内已布局的新增城市土地总面积a
sim,i,t,k
，与年度t内环带i内预期模拟新增城市土地年总面积a
exp,i,t
进行对比，按如下规则确定年度t内第k次城市土地开发行为需要在环带i内新生成的城市土地斑块总面积a
i,t,k
；
14.若a
sim,i,t,k
≥a
exp,i,t
，则a
i,t,k
取值a1(50》a1》0)；否则继续判定，若a
exp,i,t-a
sim,i,t,k
≥a2，则a
i,t,k
取值a2，否则a
i,t,k
取值a
exp,i,t-a
sim,i,t,k
；其中，a
i,t,k
表示年度t内的第k次城市土地开发行为需要在环带i内新生成的城市土地斑块总面积，a1和a2是为常数；其中，a2为a
i,t,k
的默认值，且0《a1《a2；
15.步骤5：在确定a
i,t,k
后，调用斑块生成引擎生成一定数量的城市土地斑块，直到在同心环带i内新生成的斑块总面积达到a
i,t,k
；
16.步骤6：检查年度t内所有环带中已布局的新增城市土地总面积a
sim,t
，如果a
sim,t
大于等于a
exp,t
，则进入t 1年度；如果a
sim,t
小于a
exp,t
，则重复执行步骤3-步骤5，直到a
sim,t
大于等于a
exp,t
；
17.步骤7：重复执行步骤3-步骤6，直到运行到预设模拟的最后一个年份t。
18.相对于现有技术，本发明的有益效果是：
19.(1)在宏观上控制新增城市土地需求的时空异质性分布规律，即使用高斯模型表达了城市扩张模拟中城市土地需求从城市中心到城市外围呈现的波形分布规律；
20.(2)在微观尺度上根据城市建设适宜性，使用基于斑块的城市分配引擎生成城市斑块，并可以通过调整参数控制有机生长和自然生长斑块的面积比例及斑块形状，模拟具有不同空间形态的城市景观。
附图说明
21.图1为本发明实施例的流程图；
22.图2为本发明实施例的研究期新增城市土地密度/斑块频率拟合高斯函数族示意图；
23.图3为本发明实施例的斑块生成引擎运行流程图；
24.图4为本发明实施例的有机生长/自然生长斑块示意图；
25.图5为本发明实施例的有机生长/自然生长斑块面积对数正态分布示意图；
26.图6为本发明实施例的有机生长/自然生长初始单元候选单元集示意图；
27.图7为本发明实施例的有机生长/自然生长初始单元及其邻域示意图。
具体实施方式
28.为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。
29.请见图1，本发明提供的一种耦合时空高斯模型和斑块元胞自动机的城镇扩张模拟方法，包括以下步骤：
30.步骤1：从城市中心点向外围区域创建等间距同心环，将研究区划分为n个等间距同心环带；
31.步骤2：根据年度t对应的新增城市土地密度高斯函数计算出研究区年度t内各环带i中新增城市土地预期面积a
exp,i,t
及年度t新增城市土地预期总面积a
exp,texp,t
32.本实施例中，在一个模拟周期内，各同心环带中新增城市土地密度和新增城市斑块的频率分布的观测值与同心环带到城市中心的距离经拟合满足下式表示的时空高斯函数分布。请见图2，新增城市土地密度高斯函数为：
[0033][0034]
式中：den(i,d,t)表示年度t内环带i中新增城市土地密度，d为同心环带i到城市中心距离，a
t
、b
t
、c
t
为年度t所对应高斯函数参数；
[0035]
年度t内同心环带i中新增城市土地密度等于环带i内新增城市土地面积与该环带内可开发土地面积之比，其中，环带i内可开发土地面积为该环内土地总面积扣除水体面积。
[0036]
步骤3：根据年度t对应的新增城市斑块频率高斯函数抽取一个同心环带i用于执
行土地开发行为，即生成城市土地斑块(新增城市斑块频率大小决定该环带被抽取的概率大小)；
[0037]
本实施例中，请见图2，新增城市斑块频率高斯函数为：
[0038][0039]
式中，fre(i,d,t)表示年度t内环带i中新增城市斑块频率，d为同心环带i到城市中心距离，a
t
、b
t
、c
t
为年度t所对应高斯函数参数；
[0040]
年度t内同心环带i中新增城市斑块密度等于环带i中斑块质心的数量与年度t内研究区域新增城市斑块总数之比。
[0041]
步骤4：检查年度t内第k次土地开发行为之前同心环带i内已布局的新增城市土地总面积a
sim,i,t,k
，与年度t内环带i内预期模拟新增城市土地年总面积a
exp,i,t
进行对比，按如下规则确定年度t内第k次城市土地开发行为需要在环带i内新生成的城市土地斑块总面积a
i,t,k
；
[0042]
若a
sim,i,t,k
≥a
exp,i,t
，则a
i,t,k
取值a1(50》a1》0)；否则继续判定，若a
exp,i,t-a
sim,i,t,k
≥a2，则a
i,t,k
取值a2，否则a
i,t,k
取值a
exp,i,t-a
sim,i,t,k
；其中，a
i,t,k
表示年度t内的第k次城市土地开发行为需要在环带i内新生成的城市土地斑块总面积，a1和a2是为常数；其中，a2为a
i,t,k
的默认值，且0《a1《a2；
[0043]
步骤5：在确定a
i,t,k
后，调用斑块生成引擎生成一定数量的城市土地斑块，直到在同心环带i内新生成的斑块总面积达到a
i,t,k
；
[0044]
请见图3，本实施例中，步骤5的具体实现包括以下子步骤：
[0045]
步骤5.1：自定义总面积为a
itk
的斑块中有机生长斑块和自然生长斑块的面积比例；其中，有机生长斑块表示新生成的斑块与已开发的城市土地斑块具有空间连接关系，自然生长斑块表示新生成的斑块与已开发的城市土地斑块不具有空间连接关系；
[0046]
本实施例中，斑块生成引擎可生成两类城市土地斑块：有机生长斑块和自然生长斑块。其中，有机生长斑块表示新生成的斑块与已开发的城市土地斑块具有空间连接关系，自然生长斑块表示新生成的斑块与已开发的城市土地斑块不具有空间连接关系，如图4。在准备生成一个城市土地斑块之前，斑块生成引擎通过随机方式确定该斑块所属类型。用户可自定义总面积为a
itk
的斑块中有机和自然生长斑块的面积比例，该比例可控制城市土地布局的紧凑性和连续性。
[0047]
步骤5.2：根据有机生长斑块和自然生长斑块面积各自对应的对数正态分布，随机产生一个数值作为该斑块的预期面积大小；其中，请见图5，本实施例假设城市土地斑块面积服从均值为μ方差为ε的对数正态分布；
[0048]
步骤5.3：将同心环带i全部栅格单元按城市建设适宜性概率从高到低排序，截取适宜性概率最高的前a
itk
*β个栅格单元作为斑块生长初始栅格单元的候选单元集；其中，β为用户自定义参数；
[0049]
其中，请见图6，若预备生产的栅格单元为有机生长斑块，则候选单元集与已开发的城市土地斑块具有空间连接关系；若预备生产的栅格单元为自然生长斑块，则候选单元集与已开发的城市土地斑块不具有空间连接关系；
[0050]
本实施例中，基于自然地理因素和社会经济因素评估研究区每个栅格单元的城市
建设适宜性概率。常用的自然地理因素主要是地形地貌，如高程、坡度等。常用的社会经济因素主要是区位交通(包括：到市、县、镇中心的距离；到国道、省道、高速公路、港口等交通设施的距离；到学校、医院等公共服务设施的距离等)和社会经济条件(人口和gdp等)，也有社会政策制度及规划，比如“三区三线”。常用评估方法主要有统计学习模型(支持向量机、最小二乘支持向量机)、集成学习模型(如随机森林、深度森林)、神经网络模型(如多层感知器、人工神经网络、深度学习、卷积神经网络)等。
[0051]
步骤5.4：随机从候选单元集中选择一个栅格单元，将其城市建设适宜性概率与一个0～1之间的服从均匀分布的随机数进行对比，若该栅格单元的城市建设适宜性概率大于随机数，则确定该栅格单元为新增城市斑块的初始栅格单元(例如随机数为0.75，则适宜性概率大于0.75的候选单元成为初始栅格单元，如图7)；否则继续测试候选单元集中的下一个栅格单元，直到确定一个初始栅格单元为止；
[0052]
步骤5.5：将初始栅格单元3
×
3邻域范围内的可开发栅格单元放入一个邻域单元集中，作为拟生产斑块从初始栅格单元开始扩展生长的备选单元集；使用和步骤5.4中相同的随机生存测试过程从邻域单元集中随机选择一个栅格单元作为拟生成斑块的一部分，并将该单元的3
×
3邻域范围内的可开发栅格单元放入邻域单元集中；若新加入的邻域单元首次被检测到已位于邻域单元集中，则将这些重复邻域单元的城市建设适宜性概率乘以一个[0,2]之间的小数：该小数属于[0,1]时，将使得拟生成斑块趋向于呈现条带状；该小数属于[1,2]时，将使得拟生成斑块趋向于紧凑的圆形；迭代执行上述过程，直到扩展生成的斑块达到预期面积大小；
[0053]
步骤5.6：重复执行步骤5.1-步骤5.5，直到生成的斑块总面积达到预期总面积a
itk
。
[0054]
步骤6：检查年度t内所有环带中已布局的新增城市土地总面积a
sim,t
，如果a
sim,t
大于等于a
exp,t
，则进入t 1年度；如果a
sim,t
小于a
exp,t
，则重复执行步骤3-步骤5，直到a
sim,t
大于等于a
exp,t
；
[0055]
步骤7：重复执行步骤3-步骤6，直到运行到预设模拟的最后一个年份t。
[0056]
本发明提出的耦合时空高斯模型和斑块元胞自动机的城镇扩张模拟方法，在对大都市区(武汉，1991-2017年)的研究试验中，实施后效果如下。
[0057]
1.城市建设适宜性评价；
[0058]
选取自然地理因素(高程、坡度)、社会经济因素(核密度；到城市中心、城市副中心、镇中心的距离；到机场、火车站、高速路入口、主要道路、河流的距离)以及规划约束因素(永久基本农田)，使用深度森林(deep forest)模型进行城市建设适宜性评估，模型训练精度为0.9344。
[0059]
2.时空高斯模型精度；
[0060]
时空高斯模型的拟合优度(r2)值均在0.8以上，能够较好地表达新增城市土地的时空分布规律。
[0061]
3.本发明提出模型的精度；
[0062]
将本发明提出的模型与参考模型(没有耦合时空高斯模型的相同斑块元胞自动机模型)进行对比实验，使用模式级相似性的模糊kappa(k-fuzzy)指标和单元级一致性的品质因数(fom)指标来评估模型精度。
[0063]
(1)研究期历年模型整体精度评价
[0064]
参考模型的fom指数值略高于本发明提出的模型的值，而本发明提出的模型的k-fuzzy指数值大于参考模型的值，这表明所提出的模型在单元级一致性方面产生了相对较差的模拟结果，而在捕捉城市发展的一般空间格局方面优于参考模型。
[0065]
(2)研究期历年各同心环带内精度评价
[0066]
计算了两个模型历年在每个同心环带内模拟城市景观的fom指数值，结果表明，本发明提出的模型在城市土地快速增长的内城和郊区通常比参考模型提供更高的单元级一致性，而参考模型在中心城区和城市腹地具有更好的模拟精度。这表明本发明提出的模型在模拟跨越式城市扩张方面优于参考模型，尤其在模拟由政策驱动的新城市中心的出现方面具有较大潜力。
[0067]
应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。