技术领域
[0001] 本
发明涉及一种城市土地利用模拟方法，尤其是涉及一种基于空间误差的城市土地利用模拟元胞自动机方法。
背景技术
[0002] 元胞自动机(CA)模型由空间、离散有限状态元、邻域和转换规则组成。CA模型使用局部交互来模拟复杂系统的演变，这与城市土地利用的复杂性高度契合。在过去的二十年中，大多数工作都采用经典统计学方法、
人工神经网络和元启发式
算法获取转换规则。空间驱动因子之间不仅存在多重共线性、也同时存在空间自相关，导致这些方法没有充分量化空间驱动因子对城市土地利用变化的影响。
模糊逻辑方法的参数具有主观性、而
人工智能方法得到的参数往往缺乏物理意义。进化计算和启发式算法(如
粒子群优化和
遗传算法)，不仅计算权重来量化空间驱动因子对土地利用转移概率的贡献，而且还降低了空间驱动因子之间多重共线性造成的不利影响。但是，启发式方法通常只是提供在合理计算成本下的局部最优解，搜索范围的变化可能导致不同的结果和计算成本急剧增加。
[0003] 城市土地利用过程中同时存在空间异质性和空间依赖，复杂的地理过程给模型构建带来了巨大挑战。逻辑回归模型(logit)只能提供全局的最佳拟合关系，因此不能解释复杂地理过程中的多重共线性、自相关效应和空间异质性。SEM是一种最基本的空间
自回归模型，它可用于探索影响土地利用变化的因素、解释土地利用变化过程，建模残差的分布会产生空间相关性，而将空间自相关项引入CA转换规则的获取可以有效消除模型残差中的空间自相关、提高CA建模的
精度。
发明内容
[0004] 本发明的目的就是为了克服上述
现有技术存在的
缺陷而提供一种基于空间误差的城市土地利用模拟元胞自动机方法。
[0005] 本发明的目的可以通过以下技术方案来实现：
[0006] 一种基于空间误差的城市土地利用模拟元胞自动机方法，该方法包括以下步骤：
[0007] 步骤1：对卫星遥感图像进行监督分类，获取城市土地利用图和驱动因子，并进行预处理；
[0008] 步骤2：对经过预处理后的城市土地利用图和驱动因子进行系统抽样，获取城市土地利用变化和驱动因子的样本点数据；
[0009] 步骤3：利用空间误差模型对
采样的样本点数据进行训练，建立城市土地利用的元胞自动机转换规则；
[0010] 步骤4：采用对数似然值、赤池信息准则评估转换规则的拟合性能，检验SEM对转换规则残差聚集性消除的效果，得到优化的CA转换规则；
[0011] 步骤5：利用优化的CA转换规则建立CASEM新模型，并计算产生城市土地转换概率图；
[0012] 步骤6：基于遥感分类的历史土地利用数据，使用
马尔科夫链模型预测未来每个时间点城市元胞数量；
[0013] 步骤7：结合初始土地利用图、土地转换概率、邻域影响和限制因素，使用CASEM模拟预测城市空间格局和扩张过程；
[0014] 步骤8：针对使用CASEM模型模拟预测的模拟结果利用终态和变化两类指标进行综合评价，并输出保存模拟结果。
[0015] 进一步地，所述的步骤1包括以下分步骤：
[0016] 步骤11：对所有遥感影像、矢量地图数据、栅格驱动因子数据进行几何校正和投影转换，确保所有数据的空间参考一致；
[0017] 步骤12：利用
支持向量机监督分类方法对卫星遥感影像分类，获取城市土地利用图；使用欧式距离方法对矢量数据进行计算，得到栅格格式的驱动因子图，对所有驱动因子进行归一化处理；
[0018] 步骤13：对城市土地利用图和驱动因子图进行重采样，使得所有数据的空间
分辨率一致。
[0019] 进一步地，所述的步骤3包括以下分步骤：
[0020] 步骤31：根据两个时期的城市土地利用图采样的样本点数据，得到土地利用的变化信息，生成土地变化信息因变量y：发生城市土地利用变化的点y记作1，未发生变化的y记作0；
[0021] 步骤32：利用R语言建立SEM模型，对样本点数据进行训练，建立城市土地利用的元胞自动机转换规则。
[0022] 进一步地，所述的步骤4包括以下分步骤：
[0023] 步骤41：使用Moran’s I检测引入空间回归模型后土地利用转换规则的残差中是否存还在空间自相关；
[0024] 步骤42：采用对数似然值、赤池信息准则评估转换规则的拟合性能，检验SEM对转换规则残差聚集性消除的效果，得到优化的CA转换规则。
[0025] 进一步地，所述的步骤5中CASEM新模型元胞状态的转换，其描述公式为：
[0026]
[0027] 式中， 表示t+1时刻的元胞状态， 表示t时刻的元胞状态，f表示总体转换规则，Neii表示领域的影响，Res()表示限制因素，Divi表示驱动因子的影响，Count表示元胞数量。
[0028] 进一步地，所述的步骤5中CASEM新模型的总体转换概率，其描述公式为：
[0029]
[0030] 式中，Poverall表示总体转换概率，PDiv表示驱动因子对土地利用的影响，PNei表示领域对土地利用的影响，TIP表示对转换概率的缩放，LAP表示对领域影响的缩放。
[0031] 进一步地，所述的步骤5中的CASEM新模型中的驱动因子对土地利用的影响，其描述公式为：
[0032]
[0033] 式中，β0表示模型的截距，Varj表示驱动因子j，βj表示驱动因子j的权重，I表示单位矩阵，λ表示SEM模型参数，W表示空间邻接矩阵，ε表示模型的残差。
[0034] 进一步地，所述的步骤6包括以下分步骤：
[0035] 步骤61：根据起始年份和终止年份的城市土地利用图，利用马尔科夫链模型得到历史土地利用类型之间的转换概率；
[0036] 步骤62：根据起始年份各类型的土地利用数量和土地转换概率，利用马尔科夫链模型计算未来时间点的各类型土地利用数量。
[0037] 进一步地，所述的步骤7包括以下分步骤：
[0038] 步骤71：结合初始土地利用图、土地转换概率、邻域影响和限制因素，在UrbanCA模拟
框架下，产生CASEM模拟；
[0039] 步骤72：将距离待模拟时间点最近的真实土地利用作为初始年份，利用CASEM
迭代N次，其中N表示初始与结束的年份差，模拟和预测土地利用格局和扩张过程。
[0040] 进一步地，所述的步骤8包括以下分步骤：
[0041] 步骤81：使用总体精度OA、图形优度FOM、Kappa模拟系数评估使用CASEM模型模拟预测的模拟结果；
[0042] 步骤82：将使用CASEM模型模拟预测的模拟结果与真实的初始-终止年份城市土地利用比较并评定，所比较并评定的结果包括：初始城市Initial urban、实际为非城市模拟为城市False、实际和模拟均为城市Hit、实际和模拟均为非城市CR、实际为城市模拟为非城市Miss；
[0043] 步骤83：在地理信息
软件环境下输出并保存模拟结果。
[0044] 与现有技术相比，本发明具有以下优点：
[0045] (1)本发明一种基于空间误差的城市土地利用模拟元胞自动机方法，该方法中提出了一种基于SEM的地理CA模型(CASEM)，CASEM结合空间邻接矩阵可以有效表达城市土地利用的空间自相关和空间异质性，有效地消除由于驱动因子的选取产生的空间自相关、城市扩张现象本身的空间自相关和空间异质性导致的模型残差自相关效应。本发明在CA建模获取转换规则的过程中引入空间计量经济学模型SEM构建CASEM，用于提高城市土地利用模拟模型捕捉土地利用变化与驱动因子之间关系的能
力。与传统CA模型相比，由SEM获取的转换规则的残差是随机分布的，证明CASEM可以消除模型残差中的自相关。基于空间误差的城市土地利用模拟元胞自动机方法，能够较好的模拟和预测城市土地利用变化。
[0046] (2)本发明一种基于空间误差的城市土地利用模拟元胞自动机方法，以社会经济、地形和空间临近等作为驱动因子，本发明提出的CASEM相比于传统不消除残差聚集的CA模型，在城市土地利用变化和城市土地利用模拟预测中，能够获取更准确的转换规则、模拟结果精度也更高。
附图说明
[0048] 图2为发明案例研究区域图；
[0049] 图3为本发明的案例区域影响因子图；
[0050] 图4为CAlogit和CASEM的模型残差分布图；
[0051] 图5为CAlogit模型模拟的2000-2010城市土地利用格局；
[0052] 图6为CASEM模型模拟的2000-2010城市土地利用格局；
[0053] 图7为CAlogit模型模拟的2011-2017城市土地利用格局；
[0054] 图8为CASEM模型模拟的2011-2017城市土地利用格局。
具体实施方式
[0055] 下面将结合本发明
实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都应属于本发明保护的范围。
[0056] 本发明可以通过以下技术方案来实现：
[0057] 如图1所示为本发明一种基于空间误差的城市土地利用模拟元胞自动机方法的整体流程，包括以下步骤：
[0058] 1)利用支持向量机分类器对卫星遥感图像进行监督分类，获取城市土地利用图和驱动因子，对其进行预处理；
[0059] 2)对城市土地利用图和驱动因子进行系统抽样，获取城市土地利用变化和驱动因子的样本数据；
[0060] 3)利用空间误差模型(SEM)对采样数据进行训练，建立城市土地利用的元胞自动机(CA)转换规则；
[0061] 4)采用对数似然值、赤池信息准则(AIC)评估转换规则的拟合性能，检验SEM对转换规则残差聚集性消除的效果；
[0062] 5)利用构建的CA转换规则建立CASEM新模型，并计算产生城市土地转换概率图；
[0063] 6)基于遥感分类的历史土地利用数据，使用马尔科夫链模型预测未来每个时间点城市元胞数量；
[0064] 7)结合初始土地利用图、土地转换概率、邻域影响和限制因素，使用CASEM模拟预测城市空间格局和扩张过程；
[0065] 8)利用终态和变化两类指标对模拟结果进行综合评价；
[0066] 9)输出并保存模拟结果。
[0067] 所述的步骤1)具体为：
[0068] 9)对所有遥感影像、矢量地图数据、栅格驱动因子数据进行几何校正和投影转换，确保所有数据的空间参考一致；
[0069] 10)利用支持向量机监督分类方法对卫星遥感影像分类，获取城市土地利用图；使用欧式距离方法对矢量数据进行计算，得到栅格格式的驱动因子图，对所有驱动因子进行归一化处理；
[0070] 11)对城市土地利用图和驱动因子图进行重采样，使得所有数据的空间分辨率一致，如果统一为30m分辨率；
[0071] 所述的步骤2)具体为：
[0072] 12)对预处理之后的城市土地利用和驱动因子图进行系统抽样，为CASEM模型训练提供样本数据；
[0073] 所述的步骤3)具体为；
[0074] 13)根据两个时期的城市土地利用图采样数据，得到土地利用的变化信息，生成土地变化信息因变量y：发生城市土地利用变化的点y记作1，未发生变化的y记作0；
[0075] 14)利用R语言建立SEM，对样本点数据进行训练(本发明中其特点在于将因变量y也作为一个自变量，并为该变量赋权)，获取基于SEM的城市土地利用转换规则；
[0076] 所述的步骤4)具体为；
[0077] 15)使用Moran’s I检测引入空间回归模型后土地利用转换规则的残差中是否存还在空间自相关；
[0078] 16)使用对数似然值、赤池信息准则(AIC)评估模型的拟合性能；
[0079] 所述的步骤5)具体为；
[0080] 17)利用优化的CA转换规则构建基于SEM的CA模型，并计算土地利用转换概率。方法如下：
[0081] 任意元胞i在时刻t+1的状态由其在时刻t的状态和元胞邻域决定，则元胞状态的转换可表达为：
[0082]
[0083] 式中， 表示t+1时刻的元胞状态， 表示t时刻的元胞状态，f表示总体转换规则，Neii表示领域的影响，Res()表示限制因素，Divi表示驱动因子的影响，Count表示元胞数量。
[0084] 对于邻域影响的评估，通常采用Moore m×m或者圆形作为邻域。邻域的影响可表达为：
[0085]
[0086] 式中，中心元胞i不参与计算，Cot表示m×m邻域范围内城市元胞的数量，本发明采用Moore 5×5作为元胞邻域。考虑了限制因素、邻域缩放、概率缩放后的总体转换概率可表示：
[0087]
[0088] 式中，Poverall表示总体转换概率，PDiv表示驱动因子对土地利用的影响，PNei表示领域对土地利用的影响，TIP表示对转换概率的缩放，TIP的取值范围为(0.0,0.1)，LAP表示对领域影响的缩放，LAP的取值范围为(0.5,1.0)。
[0089] 本发明中TIP设置为0.01，LAP设置为0.8。限制因素Res()表示元胞受到某种限制而不能发展并转化为城市元胞，这种限制包括大面积
水体、基本农田、生态保护区、以及公园和绿地等。Res()的取值范围为0和1，0表示不能转化为城市元胞，1表示可以转化为城市元胞。使用SEM构建的驱动因子对土地利用的影响可表达为：
[0090]
[0091] 式中，β0表示模型的截距，Varj表示驱动因子j，βj表示驱动因子j的权重，I表示单位矩阵，λ表示SEM模型参数，W表示空间邻接矩阵，ε表示模型的残差。
[0092] 本发明基于土地利用变化的一阶Rook关系构建空间邻接矩阵。测算得到元胞总体转换概率后，与设定的
阈值Pthd(其值在0～1之间)进行比较，以确定元胞i在下一时刻的状态是否发生转变：
[0093]
[0094] 所述的步骤6)具体为；
[0095] 18)根据起始年份和终止年份的城市土地利用图，利用马尔科夫链模型得到历史土地利用类型之间的转换概率；
[0096] 19)根据起始年份各类型的土地利用数量和土地转换概率，利用马尔科夫链模型计算未来时间点的各类型土地利用数量；
[0097] 所述的步骤7)具体为；
[0098] 20)结合初始土地利用图、土地转换概率、邻域影响和限制因素，在UrbanCA模拟框架下，产生CASEM模拟；
[0099] 21)将距离待模拟时间点最近的真实土地利用作为初始年份，利用CASEM迭代N次(N＝终止年份-初始年份)模拟和预测土地利用格局；
[0100] 所述的步骤8)具体为；
[0101] 22)使用总体精度(OA)、图形优度(FOM)、Kappa模拟系数评估CASEM的模拟精度；总体精度通过比较真实和模拟终止年份的城市土地利用格局，将比较结果分解为差异率和一致率；FOM通过比较真实的初始-终止年份城市土地利用，将比较结果分解为初始城市(Initial urban)、实际和模拟均为非城市(CR)、实际和模拟均为城市(Hit)、实际为非城市模拟为城市(False)、实际为城市模拟为非城市(Miss)等五个大类；Kappa模拟系数将模拟精度分解数量一致性的Ktransition和
位置一致性Ktransloc。
[0102] 所述的步骤9)具体为：
[0103] 23)在地理信息软件环境下输出并保存模拟结果。
[0104] 本发明的实际具体实施例如下：
[0105] 1)本案例的位置如图2所示。通过GlobeLand30提取郑州市2000、2010年的土地利用数据，经重分类得到郑州市2000、2010年城市土地利用格局。支持向量机方法对Landsat 8/OLI于2017年12月24日获取的郑州市遥感影像进行解译，得到2017年城市土地利用格局；
[0106] 3)获取城市位置和交通路网数据，基于欧氏距离得到到市中心、到区县中心、到
铁路、到高速路的距离，将这些距离通过经过归一化处理后，作为城市土地利用变化的影响因素。获取数字高程模型(DEM)和经济(GDP)数据，作为土地利用变化的地形和社会经济影响因素，如图3所示；
[0107] 4)对驱动因子以和土地利用数据进行几何校正和重采样，保证空间参考和空间分辨率一致；
[0108] 5)利用系统抽样方法获取郑州市2000-2010年间，城市土地利用变化信息和驱动因子的样本数据；
[0109] 3)利用R语言实现logit和SEM算法，采用已获取的样本数据校准logit和SEM算法，基于模型残差最小化和随机分布的原则，进行转换规则优化。表1为转换规则参数，表2为模型拟合性能和残差统计分布；
[0110] 4)如图4所示，SEM可以消除模型残差中的自相关，SEM建立的模型对数似然值更大并且AIC更小，表明由SEM获取的转换规则拟合能力更优；
[0111] 5)利用logit建立CAlogit、利用SEM建立CASEM并生成转换概率图，将2000年城市土地利用作为初始地图，设置迭代次数为10，设置邻域为Moore 5×5，图5和图6分别显示校准阶段模拟两个模型模拟2000-2010年城市土地利用格局；
[0112] 6)利用与步骤5相同的转换概率图，将2010年土地利用作为初始地图，设置迭代次数为7，设置邻域为Moore 5×5，得到图7和图8分别显示的验证阶段两个模型模拟的2010-2017年城市土地利用格局；
[0113] 7)表3显示了使用OA、FOM、Kappa评估校准阶段和验证阶段模型的模拟精度；很显然，在校准和验证阶段，CASEM的精度高于CAlogit模型、而错误率低于CAlogit模型，有力证明了本发明的有效性。
[0114] 表1由logit和SEM获取的CA参数
[0115]模型 截距 Durban Ddistr Drailway Dhighway GDP DEM
CAlogit 0.052 0.046 -0.099 -0.048 -0.061 0.287 -0.076
CASEM 0.045 0.052 -0.096 -0.044 -0.063 0.299 -0.076
[0116] 表2 CAlogit和CASEM转换规则的拟合效果
[0117]
[0118] 表3 CAlogit和CASEM模拟土地利用的精度
[0119]
[0120] 以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉
本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到各种等效的
修改或替换，这些修改或替换都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以
权利要求的保护范围为准。
