CN104361255A - 一种改进元胞自动机城市扩张模拟方法 - Google Patents

Abstract

本发明公开了一种改进元胞自动机城市扩张模拟方法，首先构建指标体系对城镇化发展水平以及土地利用潜力进行评价，得到一幅土地发展潜力图。基于发展潜力图，对非城市像元的潜力分值进行从大到小排序，并用数组记录下排序后的分值以及对应的像元位置。然后遍历数组的每一项，直到找到满足邻域影响的最大分值像元并将该像元转换为城市用地类型，与此同时城市元胞数量增加1，然后更新该元胞所有邻域的影响，从数组中删除该项并跳出循环，判断元胞增长数量是否满足要求，满足则得到模拟结果，不满足则重复上述循环，直到城市用地增长数量满足为止。本发明解决了传统元胞自动机因转换阈值的不确定性而导致模拟结果不确定的难题。

Description

一种改进元胞自动机城市扩张模拟方法

技术领域

本发明属于城市规划领域，具体涉及一种改进元胞自动机城市扩张模拟方法。

背景技术

城市扩张元胞自动机模型是由一系列元胞来代表不同土地利用类型，并通过邻域状态和转换规则来控制非城市用地向城市用地类型的转换，以此来模拟城市扩张的离散时空动态模型。利用CA模拟城市空间结构的演变过程是由转换规则来控制的。转换规则是元胞自动机模型的核心，它决定了元胞自动机的动态转化过程。以logistic-CA模拟城市扩张模拟为例，最终的转换规则由影响城市空间结构演变的诸多自然环境和社会经济因素决定(地形、高程、离市中心的距离、离道路的距离等空间变量以及邻近范围的影响等)，这些因素通过某种方式构成一种综合“效用值”来影响非城市用地向城市用地转变的概率。在以往的城市CA模型中，如何设置一个合理的转换阈值来决定最终的模拟结果往往是城市扩张模拟中的难点。阈值设置不同会导致最终的模拟结果也大不相同。阈值设置过大可能导致CA模拟的城市元胞增长数量永远也达不到预设的增长数量；阈值设置过小则会使得元胞迭代的次数过少导致精度大大降低。因此迫切需要一种方法来解决这种因为阈值不同而导致最终结果也不确定的难题。

发明内容

为了解决上述的技术问题，本发明提供了一种改进元胞自动机城市扩张模拟方法。

本发明所采用的技术方案是：一种改进元胞自动机城市扩张模拟方法，其特征在于，包括以下步骤：

步骤1：收集研究区域数据，包括数据集A和数据集B；

所述的数据集A包括Y1年土地利用分类图、Y2年土地利用分类图和Y1年统计年鉴数据；

所述的数据集B包括Y1年DEM、TM校正后影像、地类图斑图层数据、重要基础设施数据；

其中Y1年统计年鉴数据包括城市人口、总人口、城市建设用地面积、辖区总面积、GDP、固定资产投资额、地方财政收入、第二产业产值、第三产业产值、城市居民人均可支配收入、城市居民人均住房建筑面积、城市职工社会保障覆盖率；重要基础设施包括城市道路、城市中心、商场、学校、医院；

步骤2：利用数据集A对Y1年城镇化发展水平进行评价；

步骤3：利用数据集B对Y1年非城市用地发展为城市用地的潜力进行评价；

步骤4：将步骤2和步骤3得到的结果进行“影响叠加”，得到最终的土地发展潜力面；

步骤5：将步骤4得到的土地发展潜力面上的所有非城市元胞的潜力值按照从大到小排序，用数组记录下潜力值以及对应的元胞位置；

步骤6：从头遍历步骤5得到的数组，判断遍历的每一项元胞的邻域个数；

如果邻域个数满足增长需要的最小邻域个数，则将该元胞状态由非城市变成城市，并修改该元胞所有邻域的影响，同时数组中删除该项，转步骤7；

如果邻域个数不满足要求，则移动到数组下一项，并回转执行步骤6中的判断；

步骤7：判断模拟增长的城市元胞数量是否已经达到总控制增长数量要求；

如果达到了数量要求，则得到最终的模拟结果，转步骤8；

如果没有达到数量要求，则回转执行所述的步骤6；

步骤8：对得到的模拟结果与真实结果对比进行精度评价，包括kappa评价、生产者精度评价、用户精度评价。

作为优选，步骤2的具体实现包括以下子步骤：

步骤2.1：对城镇化发展水平评价选取三个一级指标，包括城市化率、经济产业水平、居民生活水平；并对上述三个一级指标设置一级指标权重W1、W2、W3，且W1+W2+W3＝1；

步骤2.2：对步骤2.1中所述的三个一级指标下设置二级指标；其中城市化率下设置二级指标为土地城市化率、人口城市化率，对应权重为w11、w12，且w11+w12＝1；经济产业水平下设置二级指标为人均GDP、人均固定资产投资额、人均地方财政收入、第二产业占GDP比例、第三产业占GDP比例，对应权重为w21、w22、w23、w24、w25，且w21+w22+w23+w24+w25＝1；居民生活水平下设置二级指标为城市居民人均可支配收入、城市居民人均住房建筑面积、城市职工社会保障覆盖率，对应权重为w31、w32、w33，且w31+w32+w33＝1；

步骤2.3：对上述评价二级指标计算值数据标准化后进行加权求和，得到研究区每个下一级行政区的城市发展质量评价。

作为优选，步骤3的具体实现包括以下子步骤：

步骤3.1：计算研究区每个元胞坡度值(Slope_score)、归一化植被覆盖指数(NDVI_score)、与所述的重要基础设施的距离；

步骤3.2：步骤3.1中得到的每个元胞与各个重要基础设施的距离值带入到距离衰减函数中计算每个元胞的交通通达性(Traffic_score)；

步骤3.3：分别对得到的Slope_score、NDVI_score、Traffic_score分值得按照从小到大间隔划分为5个分值段，并对每个分值段赋予1、3、5、7、9分数；

步骤3.4：将所有图斑类型划分为5大类，每类赋予一个分值(TB_score)；

步骤3.5：利用专家打分法对坡度、NDVI、交通通达性、图斑进行权重打分，分别为w_slope、w_ndvi、w_traffic、w_tb；

步骤3.6：基于上述步骤的打分分值以及对应权重，加权求和计算每个元胞土地转变潜力评价值。

作为优选，步骤3.1中的所述的坡度值采用ArcGIS的分析工具计算获取，归一化植被覆盖指数利用TM影像的三四波段采用如下公式计算获取：

NDVI＝(Band4-Band3)/(Band4+Band3)

其中NDVI表示得到的植被覆盖指数，Band4为TM第四波段值，Band3为TM第3波段值。

作为优选，步骤3.2中所述的距离衰减函数为：

$U_a=\frac{1}{5}\×e^{-{\mathrm{\β}}_1{D}_{\mathrm{road}}}+\frac{1}{5}\×e^{-{\mathrm{\β}}_2{D}_{\mathrm{center}}}+\frac{1}{5}\×e^{-{\mathrm{\β}}_3{D}_{\mathrm{cmc}}}+\frac{1}{5}\×e^{-{\mathrm{\β}}_3{D}_{\mathrm{hos}}}+\frac{1}{5}\×e^{-{\mathrm{\β}}_3{D}_{\mathrm{school}}};$

其中U_a表示交通通达性；D_road、D_center、D_cmc、D_hos、D_school分别表示所求元胞到城市道路、城市中心、商业中心、医院、学校的距离；β₁、β₂、β₃、β₄、β₅为对应的空间衰减系数。

作为优选，步骤3.5中所述的w_slope、w_ndvi、w_traffic、w_tb为：

作为优选，步骤3.6中所述的的土地转变潜力评价计算公式如下所示：

$P_{\mathrm{score}\_i}={\mathrm{\Σ}}_{j=1}^n{x}_{\mathrm{ij}}*w_{\mathrm{ij}}$

其中P_{scpre_i}为第i个元胞的土地转变潜力的评价，x_ij为第i个元胞第j项指标赋予的分值，w_ij为第i个元胞第j项指标权重,n为指标项数。

作为优选，步骤4中所述的土地发展潜力面，其土地发展潜力面的每个元胞最终值的由下面公式获取：

R_{score_i}＝P_{score_i}*δ

其中R_{score_i}为每个非城市元胞最终转变为城市元胞的潜力值，P_{score_i}为对每个非城市元胞土地转变潜力的评价值，由步骤3中的步骤3.6获取，δ为该元胞所处的区域的城市发展质量评价值，由步骤2中的步骤2.3获取。

作为优选，步骤6中所述的修改该元胞所有邻域的影响，假设一次迭代过程后，元胞集合状态由a变成b，那么将状态b中元胞M邻域内的所有邻居的邻域影响值加1，则下一次遍历到b中除元胞M外的其他元胞的任一个时候判断它新的邻域个数都以这个更新后的邻域影响值数计算。

作为优选，步骤7中所述的总控制增长数量由马尔科夫链预测方法得到。

本发明与传统元胞自动机模拟城市扩张相比迭代方式存在不同，传统元胞自动机一次迭代过程为逐行(列)遍历，将满足转变规则(阈值)的元胞变成城市元胞，遍历完所有行(列)后则一次迭代完成，一般而言一次迭代产生远大于1个的新增城市元胞。而本发明的每次迭代过程则是在整个潜力面上选择具有满足邻域要求且潜力分值最大的元胞转变成城市元胞，并调整该元胞转变状态后对它邻域的影响，一次迭代只产生一个新增城市元胞。

本发明的优点在于解决了传统元胞自动机因转换阈值的不确定性而导致模拟结果不确定的难题。本发明所改进的元胞自动机模拟结果唯一，且通过验证该方法与传统方法相比在模拟精度上也有较大的提高。

附图说明

图1：本发明的方法流程图；

图2：本发明实施例的2003年原始分类图；

图3：本发明实施例的2013年原始分类图；

图4：本发明实施例的武汉主城区七个辖区城市化发展质量评价图；

图5：本发明实施例的最终土地转变潜力面；

图6-1：本发明实施例的利用改进CA模拟的2013年结果图；

图6-2：本发明实施例的利用传统CA模拟的2013年结果图；

图7-1：本发明实施例的利用传统CA在模拟2013年结果的时候每次迭代过程kappa图；

图7-2：本发明实施例的利用改进CA在模拟2013年结果的时候每次迭代过程kappa图。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

请见图1，本发明所采用的技术方案是：一种改进元胞自动机城市扩张模拟方法，包括以下步骤：

步骤1：以武汉主城区(包括江汉区、青山区、汉阳区、武昌区、洪山区、硚口区、江岸区)为研究区，收集数据，包括数据集1：2003年、2013年土地利用分类图(见图2,图3)、2003年统计年鉴数据(包括城市人口、总人口、城市建设用地面积、辖区总面积、GDP、固定资产投资额、地方财政收入、第二产业产值、第三产业产值、城市居民人均可支配收入、城市居民人均住房建筑面积、城市职工社会保障覆盖率)；数据集2：2003年DEM、TM校正后影像、地类图斑图层数据、重要基础设施数据(城市道路、城市中心、商场、学校，医院)。

步骤2：利用数据集1对2003年武汉市主城区的七个下辖区的城镇化发展水平进行评价。评价步骤分为单指标计算、数据标准化、权重确定、综合评价四个步骤。需要计算的单指标为人口城市化率、土地城市化率、人均GDP、人均固定资产投资额、人均地方财政收入、第二产业占GDP比例、第三产业GDP比例；数据标准化为消除不同量纲影响；权重确定用专家打分法，其中一级指标为城市化率、经济产业水平、居民生活水平分别设置权重为0.25、0.60、0.15，城市化率下的二级指标土地城市化率、人口城市化率，分别设置权重为0.60、0.40；经济产业水平下二级指标人均GDP、人均固定资产投资额、人均地方财政收入、第二产业占GDP比例、第三产业占GDP比例，分别设置权重为0.25,0.25,0.10,0.25,0.15；居民生活水平下二级指标为城市居民人均可支配收入、城市居民人均住房建筑面积、城市职工社会保障覆盖率，分别设置权重为0.50,0.25,0.25。综合评价利用指标加权求和的方式，最终得到武汉主城区七个下辖区的城市化发展质量水平评价图，见图4所示。其中武昌区、汉阳区、江岸区、江汉区、洪山区、青山区、硚口区对应的城市发展评价值分别为：0.4252、0.5937、0.4666、0.661、0.2623、0.4126、0.4517。

步骤3：利用数据集2对研究区内所有非城市用地发展为城市用地的潜力进行评价。选取地形坡度、植被覆盖度、图斑类型、交通通达性作为土地转变潜力评价指标。地形坡度利用ArcGIS分析工具处理DEM获得；植被覆盖度利用TM影像的的第三和四波段按照公式就算获取；图斑类型由遥感影像解译获取；交通通达性由ArcGIS的Near工具算的的待求点到城市道路、城市中心、医院、学校、商业中心的距离代入到距离衰减函数计算获得。设置上述指标的权重分别为0.10、0.15、0.40、0.35。将得到的指标分值按照如下表所示进行附分值：

则每个非城市元胞的土地转变潜力评价可由下述公式计算所得：

$P_{\mathrm{score}\_i}={\mathrm{\Σ}}_{j=1}^n{x}_{\mathrm{ij}}*w_{\mathrm{ij}}$

其中P_{score_i}为第i个元胞的土地转变潜力的评价，x_ij为第i个元胞第j项指标的赋予的分值，w_ij为第i个元胞第j项指标权重,n为指标项数。

步骤4：将步骤2和步骤3得到的结果进行“影响叠加”，按照如下公式计算得到每个非城市元胞最终的土地转变潜力分值从而得到最终的土地发展潜力面，见图5。

R_{score_i}＝P_{score_i}*δ

其中R_{score_i}为每个非城市元胞最终转变为城市元胞的潜力值，P_{score_i}为对每个非城市元胞土地转变潜力的评价值。δ为该元胞所处的区域的城市发展质量评价值。

步骤5：将步骤4得到的土地发展潜力面上的所有非城市元胞的潜力值按照从大到小排序，用动态数组List记录下潜力值以及对应的元胞位置(以元胞所处的行列号表示)。

步骤6：从头遍历步骤5得到的数组List，判断遍历的每一项元胞的邻域个数。判断邻域个数满足增长需要的最小邻域个数，则将该元胞状态由非城市变成城市(即将该元胞值设为城市元胞的像元值)，并修改该元胞所有邻域的影响(该元胞所有邻居的邻域影响值加1，下一次遍历到这些邻居中的任一个时候判断它新的邻域个数都以这个更新后的邻域影响值数计算)，同时List中数组中删除该项。转步骤7；

如果邻域个数不满足要求，则移动到数组下一项，重复步骤6判断条件。

步骤7：跳出步骤6循环，判断模拟增长的城市元胞数量是否已经达到总控制增长数量(马尔科夫预测结果为增长数量为173950个)要求；

如果达到了数量要求，则得到最终的模拟结果(见图6-1)，同时为了与传统CA方法对对比，我们也用传统CA模拟了城市扩张结果(见图6-2)。转步骤8；

如果没有达到数量要求，则重复步骤6-7。

步骤8：对得到的模拟结果(2013年模拟结果)与真实结果(2013年真实分类结果)对比进行精度评价，并与传统元胞自动机结果做对比，见图7-1和图7-2。由kappa检验结果可知改进型元胞自动机的kappa达到0.8129919，与之对比的传统的元胞自动机的kappa只有0.736517，因此模拟精度上提高将近8个百分点。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

Claims (10)

1.一种改进元胞自动机城市扩张模拟方法，其特征在于，包括以下步骤：

步骤1：收集研究区域数据，包括数据集A和数据集B；

所述的数据集A包括Y1年土地利用分类图、Y2年土地利用分类图和Y1年统计年鉴数据；

所述的数据集B包括Y1年DEM、TM校正后影像、地类图斑图层数据、重要基础设施数据；

其中Y1年统计年鉴数据包括城市人口、总人口、城市建设用地面积、辖区总面积、GDP、固定资产投资额、地方财政收入、第二产业产值、第三产业产值、城市居民人均可支配收入、城市居民人均住房建筑面积、城市职工社会保障覆盖率；重要基础设施包括城市道路、城市中心、商场、学校、医院；

步骤2：利用数据集A对Y1年城镇化发展水平进行评价；

步骤3：利用数据集B对Y1年非城市用地发展为城市用地的潜力进行评价；

步骤4：将步骤2和步骤3得到的结果进行“影响叠加”，得到最终的土地发展潜力面；

步骤5：将步骤4得到的土地发展潜力面上的所有非城市元胞的潜力值按照从大到小排序，用数组记录下潜力值以及对应的元胞位置；

步骤6：从头遍历步骤5得到的数组，判断遍历的每一项元胞的邻域个数；

如果邻域个数满足增长需要的最小邻域个数，则将该元胞状态由非城市变成城市，并修改该元胞所有邻域的影响，同时数组中删除该项，转步骤7；

如果邻域个数不满足要求，则移动到数组下一项，并回转执行步骤6中的判断；

步骤7：判断模拟增长的城市元胞数量是否已经达到总控制增长数量要求；

如果达到了数量要求，则得到最终的模拟结果，转步骤8；

如果没有达到数量要求，则回转执行所述的步骤6；

步骤8：对得到的模拟结果与真实结果对比进行精度评价，包括kappa评价、生产者精度评价、用户精度评价。

2.根据权利要求1所述的改进元胞自动机城市扩张模拟方法，其特征在于：步骤2的具体实现包括以下子步骤：

步骤2.1：对城镇化发展水平评价选取三个一级指标，包括城市化率、经济产业水平、居民生活水平；并对上述三个一级指标设置一级指标权重W1、W2、W3，且W1+W2+W3＝1；

步骤2.2：对步骤2.1中所述的三个一级指标下设置二级指标；其中城市化率下设置二级指标为土地城市化率、人口城市化率，对应权重为w11、w12，且w11+w12＝1；经济产业水平下设置二级指标为人均GDP、人均固定资产投资额、人均地方财政收入、第二产业占GDP比例、第三产业占GDP比例，对应权重为w21、w22、w23、w24、w25，且w21+w22+w23+w24+w25＝1；居民生活水平下设置二级指标为城市居民人均可支配收入、城市居民人均住房建筑面积、城市职工社会保障覆盖率，对应权重为w31、w32、w33，且w31+w32+w33＝1；

步骤2.3：对上述评价二级指标计算值数据标准化后进行加权求和，得到研究区每个下一级行政区的城市发展质量评价。

3.根据权利要求2所述的改进元胞自动机城市扩张模拟方法，其特征在于：步骤3的具体实现包括以下子步骤：

步骤3.1：计算研究区每个元胞坡度值(Slope_score)、归一化植被覆盖指数(NDVI_score)、与所述的重要基础设施的距离；

步骤3.2：步骤3.1中得到的每个元胞与各个重要基础设施的距离值带入到距离衰减函数中计算每个元胞的交通通达性(Traffic_score)；

步骤3.3：分别对得到的Slope_score、NDVI_score、Traffic_score分值得按照从小到大间隔划分为5个分值段，并对每个分值段赋予1、3、5、7、9分数；

步骤3.4：将所有图斑类型划分为5大类，每类赋予一个分值(TB_score)；

步骤3.5：利用专家打分法对坡度、NDVI、交通通达性、图斑进行权重打分，分别为w_slope、w_ndvi、w_traffic、w_tb；

步骤3.6：基于上述步骤的打分分值以及对应权重，加权求和计算每个元胞土地转变潜力评价值。

4.根据权利要求3所述的改进元胞自动机城市扩张模拟方法，其特征在于：步骤3.1中的所述的坡度值采用ArcGIS的分析工具计算获取，归一化植被覆盖指数利用TM影像的三四波段采用如下公式计算获取：

NDVI＝(Band4-Band3)/(Band4+Band3)

其中NDVI表示得到的植被覆盖指数，Band4为TM第四波段值，Band3为TM第3波段值。

5.根据权利要求3所述的改进元胞自动机城市扩张模拟方法，其特征在于：步骤3.2中所述的距离衰减函数为：

$U_a=\frac{1}{5}\×e^{-{\mathrm{\β}}_1{D}_{\mathrm{road}}}+\frac{1}{5}\×\frac{1}{5}{e}^{-{\mathrm{\β}}_2{D}_{\mathrm{center}}}+\frac{1}{5}\×e^{-{\mathrm{\β}}_3{D}_{\mathrm{cmc}}}+\frac{1}{5}\×e^{-{\mathrm{\β}}_3{D}_{\mathrm{hos}}}+\frac{1}{5}\×e^{-{\mathrm{\β}}_3{D}_{\mathrm{school}}};$

其中U_a表示交通通达性；D_road、D_center、D_cmc、D_hos、D_school分别表示所求元胞到城市道路、城市中心、商业中心、医院、学校的距离；β₁、β₂、β₃、β₄、β₅为对应的空间衰减系数。

6.根据权利要求3所述的改进元胞自动机城市扩张模拟方法，其特征在于：步骤3.5中所述的w_slope、w_ndvi、w_traffic、w_tb为：

7.根据权利要求3所述的改进元胞自动机城市扩张模拟方法，其特征在于：步骤3.6中所述的的土地转变潜力评价计算公式如下所示：

$P_{\mathrm{score}\_i}={\mathrm{\Σ}}_{j=1}^n{x}_{\mathrm{ij}}*w_{\mathrm{ij}}$

其中P_{score_i}为第i个元胞的土地转变潜力的评价，x_ij为第i个元胞第j项指标赋予的分值，w_ij为第i个元胞第j项指标权重,n为指标项数。

8.根据权利要求3所述的改进元胞自动机城市扩张模拟方法，其特征在于：步骤4中所述的土地发展潜力面，其土地发展潜力面的每个元胞最终值的由下面公式获取：

R_{score_i}＝P_{score_i}*δ

其中R_{score_i}为每个非城市元胞最终转变为城市元胞的潜力值，P_{score_i}为对每个非城市元胞土地转变潜力的评价值，由步骤3中的步骤3.6获取，δ为该元胞所处的区域的城市发展质量评价值，由步骤2中的步骤2.3获取。

9.根据权利要求1所述的改进元胞自动机城市扩张模拟方法，其特征在于：步骤6中所述的修改该元胞所有邻域的影响，假设一次迭代过程后，元胞集合状态由a变成b，那么将状态b中元胞M邻域内的所有邻居的邻域影响值加1，则下一次遍历到b中除元胞M外的其他元胞的任一个时候判断它新的邻域个数都以这个更新后的邻域影响值数计算。

10.根据权利要求1所述的改进元胞自动机城市扩张模拟方法，其特征在于：步骤7中所述的总控制增长数量由马尔科夫链预测方法得到。