CN112131731A - 一种基于空间特征向量滤波的城市生长元胞模拟方法 - Google Patents

Abstract

本发明涉及一种基于空间特征向量滤波的城市生长元胞模拟方法，包括：利用马氏距离监督分类方法对卫星遥感图像进行监督分类，获取城市土地利用图和驱动因子，对其进行预处理；获取城市土地利用变化和驱动因子的样本数据；利用空间特征向量变化系数模型M‑SVC对采样数据进行训练，建立城市土地利用的元胞自动机CA转换规则；检验M‑SVC对转换规则残差聚集性消除的效果；利用CA转换规则建立M‑SVC‑CA模型，并计算产生城市土地转换概率图；使用M‑SVC‑CA模型模拟预测城市空间格局和扩张过程；对模拟结果进行评价；利用终态和变化两类指标对模拟结果进行综合评价；输出并保存模拟结果。与现有技术相比，本发明具有能够消除模型残差自相关效应、模拟精度高等优点。

Description

一种基于空间特征向量滤波的城市生长元胞模拟方法

技术领域

本发明涉及城市生长模拟技术领域，尤其是涉及一种基于空间特征向量滤波的城市生长元胞模拟方法。

背景技术

元胞自动机(CA)模型由空间、离散有限状态元、邻域、转换规则和时间组成。CA模型使用局部交互来模拟复杂***的演变，这与城市土地利用的复杂性高度契合。在过去的二十年中，大多数工作都采用经典统计学方法、人工神经网络和元启发式算法获取转换规则。

空间驱动因子之间不仅存在多重共线性、也同时存在空间自相关，导致这些方法没有充分量化空间驱动因子对城市土地利用变化的影响。模糊逻辑方法的参数具有主观性、而人工智能方法得到的参数往往缺乏物理意义。进化计算和启发式算法，如粒子群优化和遗传算法，不仅计算权重来量化空间驱动因子对土地利用转移概率的贡献，而且还降低了空间驱动因子之间多重共线性造成的不利影响。但是，启发式方法通常只是提供在合理计算成本下的局部最优解，搜索范围的变化可能导致不同的结果和计算成本急剧增加。城市土地利用过程中同时存在空间异质性、空间不稳定性和空间依赖，复杂的地理过程给模型构建带来了巨大挑战。常见的逻辑回归模型(logit)只能提供全局的最佳拟合关系，因此不能解释复杂地理过程中的多重共线性、自相关效应、空间不稳定性和空间异质性。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种能够消除模型残差自相关效应、模拟精度高的基于空间特征向量滤波的城市生长元胞模拟方法。

本发明的目的可以通过以下技术方案来实现：

一种基于空间特征向量滤波的城市生长元胞模拟方法，包括：

步骤1：利用马氏距离监督分类方法对卫星遥感图像进行监督分类，获取城市土地利用图和驱动因子，对其进行预处理；

步骤2：获取城市土地利用变化和驱动因子的样本数据；

步骤3：利用空间特征向量变化系数模型M-SVC对采样数据进行训练，建立城市土地利用的元胞自动机CA转换规则；

步骤4：检验M-SVC对转换规则残差聚集性消除的效果；

步骤5：利用CA转换规则建立M-SVC-CA模型，并计算产生城市土地转换概率图；

步骤6：使用M-SVC-CA模型模拟预测城市空间格局和扩张过程；

步骤7：利用总体精度OA和图形优度FOM对模拟结果进行综合评价；

步骤8：利用终态和变化两类指标对模拟结果进行综合评价；

步骤9：输出并保存模拟结果。

优选地，所述的步骤1具体为：

步骤1-1：对所有遥感影像、矢量地图数据、栅格驱动因子数据进行几何校正和投影转换，确保所有数据的空间参考一致；

步骤1-2：利用马氏距离监督分类方法对卫星遥感影像分类，获取城市土地利用图；使用欧式距离方法对矢量数据进行计算，得到栅格格式的驱动因子图，对所有驱动因子进行归一化处理；

步骤1-3：运用方差膨胀因子VIF法筛选因子，去除因子间的多重共线性；

步骤1-4：对城市土地利用图和驱动因子图进行重采样，使得所有数据的空间分辨率一致，如统一为30m分辨率。

优选地，所述的步骤2具体为：

对预处理之后的城市土地利用和驱动因子图进行***采样，为M-SVC-CA模型训练提供样本数据。

优选地，所述的步骤3具体为：

步骤3-1：根据两个时期的城市土地利用图采样数据，得到土地利用的变化信息，生成土地变化信息因变量y：发生城市土地利用变化的点y记作1，未发生变化的y记作0；

步骤3-2：利用R语言从样本点的空间权阵中提取并筛选空间特征向量及对应的空间特征值；

步骤3-3：利用R语言，结合提取到的空间特征向量及特征值建立M-SVC模型，对样本点数据进行训练，获取基于M-SVC的城市土地利用转换规则。

更加优选地，所述的步骤3-2具体为：

构建样本点的空间权重矩阵C：

其中，C_(i,j)为C中的第(i,j)项；d_(i,j)为样本点i与j之间的欧氏距离；r为覆盖样本点的最小生成树中最长的距离；

构建空间权重矩阵C的变形矩阵MCM并进行分解：

MCM＝EΛE^T

其中，M＝(I-11^T/n)，M为以某一变量为中心的投影矩阵，I和1分别为n*n代为矩阵和n*1单位向量，n为样本容量；E为提取得到的特征向量组成的矩阵，E＝(E₁,E₂,E₃,…,E_n)；Λ为以特征值λ＝(λ₁,λ₂,λ₃,…,λ_n)为对角线的n*n对角矩阵。

优选地，所述的步骤4具体为：

步骤4-1：使用赤池信息准则AIC评估模型的拟合性能；

步骤4-2：使用残差的统计学指数检测引入空间特征滤波后土地利用转换规则的残差空间分布情况。

优选地，所述的步骤5具体为：

利用优化的CA转换规则构建基于M-SVC的CA模型，并计算土地利用转换概率，具体方法如下：

任意元胞i在时刻t+1的状态由其在时刻t的状态和元胞邻域决定，则元胞状态的转换可表达为：

其中，f()是总体转换规则；Div_i代表驱动因子的影响；Nei_i代表邻域的影响；Res()代表限制因素；Count是元胞数量；

对于邻域影响的评估，通常采用Moorem×m或者圆形作为邻域，邻域的影响可表达为：

其中，Cot表示m×m领域范围内城市元胞的数量；

考虑了限制因素、邻域缩放、概率缩放后的总体转换概率可表示：

其中，TIP的取值范围为(0,0.1)，表示对转换概率的缩放；LAP的取值范围为(0.5,1.0)，表示对邻域影响的缩放；Res()的取值为0或1，Res()＝0表示不能转化为城市元胞，Res()＝1表示可以转化为城市元胞；

使用M-SVC构建的驱动因子对土地利用的影响P_Div可表达为：

其中，y为因变量；x_k和β_k分别为自变量与对应系数向量；E为选取空间特征向量组成的矩阵；ε为随机噪声；

自变量x_k包含空间恒定部分与空间变化部分，对应回归系数向量β_k也包括空间恒定部分β_k,01及随空间变化部分Eγ_k；γ_k与γ均为满足正态分布的随机系数，其对应的矩阵Λ(α)与Λ(α_k)分别为特征值乘以比例常数α或α_k后组成的对称矩阵；

测算得到元胞全局转换概率P_overall后，与设定的阈值P_thd进行比较，以确定元胞i在下一时刻的状态是否发生转变：

其中，未知的阈值P_thd采用模拟逼近的方法进行探测，即先拟定一个阈值对目标土地进行判断，模拟出该阈值下的土地利用预测状况，若预测的城市元胞数大于目标年份的实际城市元胞数，则提高阈值；若预测的城市元胞数小于实际城市元胞数，则降低阈值，再进行判定模拟，直至预测元胞数与实际元胞数误差在1％以内，则可以确定该阈值及对应的模拟结果符合要求。

优选地，所述的步骤6具体为：

步骤6-1：结合初始土地利用图、土地转换概率、邻域影响和限制因素，在UrbanCA模拟框架下，产生M-SVC-CA模拟；

步骤6-2：将距离待模拟时间点最近的真实土地利用作为初始年份，利用M-SVC-CA迭代N次，N＝终止年份-初始年份，模拟和预测土地利用格局。

优选地，所述的步骤7具体为：

使用总体精度OA、Kappa模拟系数、图形优度FOM评估M-SVC-CA的模拟精度；

所述的总体精度OA通过比较终止年份的真实和模拟城市土地利用格局获取；

所述的图形优度FOM通过比较真实和模拟的初始-终止年份城市土地利用获取，将比较结果分解为初始城市Initial urban、实际和模拟均为非城市CR、实际和模拟均为城市Hit、实际为非城市模拟为城市False、实际为城市模拟为非城市Miss五类；

所述的Kappa模拟系数通过将模拟精度分解为数量一致性的K_transition和位置一致性K_transloc获取。

优选地，所述的步骤9具体为：在既有的地理信息软件环境下输出并保存模拟结果。

与现有技术相比，本发明具有以下优点：

一、消除模型空间自相关效应：本发明中的城市生长元胞模拟方法提出了一种基于ESF和空间非平稳M-SVC方法的地理CA模型M-SVC-CA，M-SVC-CA结合空间特征向量滤波，将城市扩张的变量划分为空间相关部分与空间无关部分，使变量独立于相邻区域的影响而进行研究；两类变量的分离保证了线性回归模型中的残差独立性，可以有效表达城市土地利用的空间自相关，有效地消除由于驱动因子的选取产生的空间自相关、城市扩张现象本身的空间自相关和空间异质性导致的模型残差自相关效应。

二、贴合地理要素表达的空间差异性：选用空间非平稳性模型M-SVC，将地理坐标作为自变量之一，在M-SVC-CA模型中解算得到随空间变化的因子系数，通过因子系数的空间分布可以对各个因子的作用进行分析解释，并分区域进行剖析，进而对城市扩张的规则进行深一步挖掘。

三、模拟精度高：本发明中的城市生长元胞模拟方法在CA建模获取转换规则的过程中引入空间特征向量滤波变化系数M-SVC构建M-SVC-CA，用于提高城市土地利用模拟模型捕捉土地利用变化与驱动因子之间关系的能力，对城市扩张过程及扩张结果的预测能力均有显著的提高，有效提高了模型的模拟精度。

附图说明

图1为本发明中城市生长元胞模拟方法的流程示意图；

图2为本发明实施例中研究区域的示意图；

图3为本发明实施例中区域影响因子的示意图；

图4为本发明实施例中M-SVC-CA模型的驱动因子作用权重分布示意图；

图5为本发明实施例中M-SVC-CA、ESF-CA和LR-CA模型模拟的2014年城市土地利用格局示意图；

图6为本发明实施例中M-SVC-CA、ESF-CA和LR-CA模型模拟的2014年城市土地利用格局的比较评估示意图；

图7为本发明实施例中M-SVC-CA、ESF-CA和LR-CA模型模拟的2019年城市土地利用格局；

图8为本发明实施例中M-SVC-CA、ESF-CA和LR-CA模型模拟的2019年城市土地利用格局的比较评估示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都应属于本发明保护的范围。

一种基于空间特征向量滤波的城市生长元胞模拟方法，其流程如图1所示，包括：

步骤1：利用马氏距离监督分类方法对卫星遥感图像进行监督分类，获取城市土地利用图和驱动因子，对其进行预处理，具体为：

步骤1-1：对所有遥感影像、矢量地图数据、栅格驱动因子数据进行几何校正和投影转换，确保所有数据的空间参考一致；

步骤1-2：利用马氏距离监督分类方法对卫星遥感影像分类，获取城市土地利用图；使用欧式距离方法对矢量数据进行计算，得到栅格格式的驱动因子图，对所有驱动因子进行归一化处理；

步骤1-3：运用方差膨胀因子VIF法筛选因子，去除因子间的多重共线性；

步骤1-4：对城市土地利用图和驱动因子图进行重采样，使得所有数据的空间分辨率一致。

步骤2：获取城市土地利用变化和驱动因子的样本数据，具体为：对预处理之后的城市土地利用和驱动因子图进行***采样，为M-SVC-CA模型训练提供样本数据。

步骤3：利用空间特征向量变化系数模型M-SVC对采样数据进行训练，建立城市土地利用的元胞自动机CA转换规则，具体为：

步骤3-1：根据两个时期的城市土地利用图采样数据，得到土地利用的变化信息，生成土地变化信息因变量y：发生城市土地利用变化的点y记作1，未发生变化的y记作0；

步骤3-2：利用R语言从样本点的空间权阵中提取并筛选空间特征向量及对应的空间特征值，具体为：

构建样本点的空间权重矩阵C：

其中，C_(i,j)为C中的第(i,j)项；d_(i,j)为样本点i与j之间的欧氏距离；r为覆盖样本点的最小生成树中最长的距离；

构建空间权重矩阵C的变形矩阵MCM并进行分解：

MCM＝EΛE^T

其中，M＝(I-11^T/n)，M为以某一变量为中心的投影矩阵，I和1分别为n*n代为矩阵和n*1单位向量，n为样本容量；E为提取得到的特征向量组成的矩阵，E＝(E₁,E₂,E₃,…,E_n)；Λ为以特征值λ＝(λ₁,λ₂,λ₃,…,λ_n)为对角线的n*n对角矩阵。

步骤3-3：利用R语言，结合提取到的空间特征向量及特征值建立M-SVC模型，对样本点数据进行训练，获取基于M-SVC的城市土地利用转换规则。

步骤4：检验M-SVC对转换规则残差聚集性消除的效果，具体为：

步骤4-1：使用赤池信息准则AIC评估模型的拟合性能；

步骤4-2：使用残差的统计学指数，Moran’s I、z分数、p值等，检测引入空间特征滤波后土地利用转换规则的残差空间分布情况。

步骤5：利用CA转换规则建立M-SVC-CA模型，并计算产生城市土地转换概率图，具体为：

利用优化的CA转换规则构建基于M-SVC的CA模型，并计算土地利用转换概率，具体方法如下：

任意元胞i在时刻t+1的状态由其在时刻t的状态和元胞邻域决定，则元胞状态的转换可表达为：

其中，f()是总体转换规则；Div_i代表驱动因子的影响；Nei_i代表邻域的影响；Res()代表限制因素；Count是元胞数量；

对于邻域影响的评估，通常采用Moorem×m或者圆形作为邻域，邻域的影响可表达为：

其中，Cot表示m×m领域范围内城市元胞的数量；

考虑了限制因素、邻域缩放、概率缩放后的总体转换概率可表示：

其中，TIP的取值范围为(0,0.1)，表示对转换概率的缩放；LAP的取值范围为(0.5,1.0)，表示对邻域影响的缩放；

限制因素Res()表示元胞受到某种限制而不能发展并转化为城市元胞，这种限制包括大面积水体、基本农田、生态保护区、以及公园和绿地等，Res()的取值为0或1，Res()＝0表示不能转化为城市元胞，Res()＝1表示可以转化为城市元胞；

使用M-SVC构建的驱动因子对土地利用的影响P_Div可表达为：

其中，y为因变量；x_k和β_k分别为自变量与对应系数向量；E为选取空间特征向量组成的矩阵；ε为随机噪声；

自变量x_k包含空间恒定部分与空间变化部分，对应回归系数向量β_k也包括空间恒定部分β_k,01及随空间变化部分Eγ_k；γ_k与γ均为满足正态分布的随机系数，其对应的矩阵Λ(α)与Λ(α_k)分别为特征值乘以比例常数α或α_k后组成的对称矩阵。故系数向量中的空间变化部分可以对自变量作用的空间差异性进行描述，从而解决模型的不稳定问题。

测算得到元胞全局转换概率P_overall后，与设定的阈值P_thd进行比较，以确定元胞i在下一时刻的状态是否发生转变：

其中，未知的阈值P_thd采用模拟逼近的方法进行探测，即先拟定一个阈值对目标土地进行判断，模拟出该阈值下的土地利用预测状况，若预测的城市元胞数大于目标年份的实际城市元胞数，则提高阈值；若预测的城市元胞数小于实际城市元胞数，则降低阈值，再进行判定模拟，直至预测元胞数与实际元胞数误差在1％以内，则可以确定该阈值及对应的模拟结果符合要求。

步骤6：使用M-SVC-CA模型模拟预测城市空间格局和扩张过程，具体为：

步骤6-1：结合初始土地利用图、土地转换概率、邻域影响和限制因素，在UrbanCA模拟框架下，产生M-SVC-CA模拟；

步骤6-2：将距离待模拟时间点最近的真实土地利用作为初始年份，利用M-SVC-CA迭代N次，N＝终止年份-初始年份，模拟和预测土地利用格局。

步骤7：利用总体精度OA和图形优度FOM对模拟结果进行综合评价，具体为：

使用总体精度OA、Kappa模拟系数、图形优度FOM评估M-SVC-CA的模拟精度；

所述的总体精度OA通过比较终止年份的真实和模拟城市土地利用格局获取；

所述的图形优度FOM通过比较真实和模拟的初始-终止年份城市土地利用获取，将比较结果分解为初始城市Initial urban、实际和模拟均为非城市CR、实际和模拟均为城市Hit、实际为非城市模拟为城市False、实际为城市模拟为非城市Miss五类；

所述的Kappa模拟系数通过将模拟精度分解为数量一致性的K_transition和位置一致性K_transloc获取。

步骤8：利用终态和变化两类指标对模拟结果进行综合评价；

步骤9：在既有的地理信息软件环境下输出并保存模拟结果。

M-SVC是一种基于空间特征向量滤波的空间统计学方法。空间特征向量滤波(ESF)，又称Moran特征向量映射，是一种在空间自相关的情况下估计回归系数的方法。该滤波将原始变量分为空间相关部分和非空间相关部分，使变量独立于相邻区域的影响进行研究，能够有效消除空间相关性的影响。同时，线性回归模型的空间非平稳性模型，提供了自变量系数随空间变化情况下的回归分析，可以对驱动因子作用的空间差异性进行捕捉。将ESF引入CA转换规则的获取可以有效消除模型残差中的空间自相关，同时考虑了空间非平稳性，可以提高CA建模的精度。

下面结合附图和具体实施例对本发明进行详细说明。

(1)本实施例中的研究区域如图2所示，从USGS EarthExplorer获取原始遥感影像，包括于2009年4月4日拍摄的两幅Landsat-7 TM图像(30m)，于2014年3月16日拍摄的两幅Landsat-8 OLI图像(30m)，以及于2019年4月15日拍摄的两幅Landsat-8 OLI图像(30m)。所得图像在Envi 5.3环境中经过辐射校正、图像镶嵌处理并按苏州市轮廓裁剪，分别得到2009、2014、2019三个年份的苏州市遥感影像，其中2009年图像受卫星设备影响存在黑色条带，故预先进行了去条带处理。利用马氏距离方法进行采样后的监督分类，将区域用地分为城市用地、非城市用地及水体部分，并在ArcGIS中对分类结果进行目视手动改正；为方便后续城市变化分析，对三幅地图进行水体部分恒定，城市用地不可逆向转换为非城市用地的处理，得到2009、2014、2019年的苏州市城市土地利用格局

(2)获取城市位置和交通路网数据，基于欧氏距离得到市中心、到城镇中心、到地铁线路、到火车站、到高速公路、到主要道路、到日常设施、到旅游景点等的距离，将这些距离通过经过归一化处理后，作为城市土地利用变化的影响因素。获取地形(DEM)、人口和经济(GDP)数据，作为土地利用变化的地形和社会经济影响因素，如图3所示。

(3)对得到的11项驱动因子进行筛选，采取方差膨胀因子VIF法，使用ArcGIS中的OLS分析工具，设定自变量为驱动因子，因变量为苏州土地变化观测值(若由非城市用地转为城市用地则值取1，若非城市用地不变化则为，计算各驱动因子的VIF值。设定阈值为7.5，每一次计算后剔除VIF最大的变量并重新计算VIF值，再进行剔除，直至所有变量VIF均小于7.5，则可以认为因子间基本不存在明显的多重共线性。此外，由于苏州市地形较平坦，坡度因子对城市变化影响较小，故筛除该变量。驱动因子作用权重分布图如图4所示，初始因子及筛选后因子的VIF值如表1所示。

表1驱动因子VIF值

(4)对驱动因子以和土地利用数据进行几何校正和重采样，保证空间参考和空间分辨率一致。

(5)利用***抽样方法获取苏州市2009-2014年间，城市土地利用变化信息和驱动因子的样本数据。

(6)利用R语言实现M-SVC、ESF和LR算法，采用已获取的样本数据校准三类算法，基于模型残差最小化和随机分布的原则，进行转换规则优化。表2为转换规则参数，表3为模型拟合性能和残差统计分布。

表2由三类模型确定的CA转换参数

表3三类模型转换规则的拟合效果

指数	M-SVC	ESF	LR
				Moran’s I	0.005511	0.006266	0.081179
Z	0.529466	0.599084	7.510485
				P	0.596482	0.549117	0
AIC	2309.410	2344.240	3115.478523
				Adjusted R<sup>2</sup>	0.084014	0.082126	0.027917

(7)建立M-SVC-CA、ESF-CA和LR-CA模型，生成转换概率图，将2009年城市土地利用作为初始地图，设置迭代次数为11，设置邻域为Moore5×5，探测确定2014年所需阈值，图5和图6分别显示三类模型模拟2014年城市土地利用格局的结果和评估。

(8)利用与步骤5相同的转换概率图，仍以2009年土地利用作为初始地图，设置迭代次数为11，设置邻域为Moore5×5，但减小阈值后重新探测并确定，得到图7和图8分别显示的三类模型模拟的2019年城市土地利用格局。

(9)表4和表5显示了2014年和2019年使用OA、FOM、Kappa评估校准阶段和验证阶段模型的模拟精度；很显然，M-SVC-CA对变化过程及变化结果的模拟精度均要高于ESF-CA和LR-CA模型，有力证明了本发明的有效性。

表4 2014年三类模型模拟土地利用的精度

表5 2019年三类模型模拟土地利用的精度

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到各种等效的修改或替换，这些修改或替换都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。

Claims (10)

1.一种基于空间特征向量滤波的城市生长元胞模拟方法，其特征在于，包括：

步骤1：利用马氏距离监督分类方法对卫星遥感图像进行监督分类，获取城市土地利用图和驱动因子，对其进行预处理；

步骤2：获取城市土地利用变化和驱动因子的样本数据；

步骤3：利用空间特征向量变化系数模型M-SVC对采样数据进行训练，建立城市土地利用的元胞自动机CA转换规则；

步骤4：检验M-SVC对转换规则残差聚集性消除的效果；

步骤5：利用CA转换规则建立M-SVC-CA模型，并计算产生城市土地转换概率图；

步骤6：使用M-SVC-CA模型模拟预测城市空间格局和扩张过程；

步骤7：对模拟结果进行评价；

步骤8：利用终态OA和变化FOM两类指标对模拟结果进行综合评价；

步骤9：输出并保存模拟结果。

2.根据权利要求1所述的一种基于空间特征向量滤波的城市生长元胞模拟方法，其特征在于，所述的步骤1具体为：

步骤1-1：对所有遥感影像、矢量地图数据、栅格驱动因子数据进行几何校正和投影转换，确保所有数据的空间参考一致；

步骤1-2：利用马氏距离监督分类方法对卫星遥感影像分类，获取城市土地利用图；使用欧式距离方法对矢量数据进行计算，得到栅格格式的驱动因子图，对所有驱动因子进行归一化处理；

步骤1-3：运用方差膨胀因子VIF法筛选因子，去除因子间的多重共线性；

步骤1-4：对城市土地利用图和驱动因子图进行重采样，使得所有数据的空间分辨率一致。

3.根据权利要求1所述的一种基于空间特征向量滤波的城市生长元胞模拟方法，其特征在于，所述的步骤2具体为：

对预处理之后的城市土地利用和驱动因子图进行***采样，为M-SVC-CA模型训练提供样本数据。

4.根据权利要求1所述的一种基于空间特征向量滤波的城市生长元胞模拟方法，其特征在于，所述的步骤3具体为：

步骤3-1：根据两个时期的城市土地利用图采样数据，得到土地利用的变化信息，生成土地变化信息因变量y：发生城市土地利用变化的点y记作1，未发生变化的y记作0；

步骤3-2：利用R语言从样本点的空间权阵中提取并筛选空间特征向量及对应的空间特征值；

步骤3-3：利用R语言，结合提取到的空间特征向量及特征值建立M-SVC模型，对样本点数据进行训练，获取基于M-SVC的城市土地利用转换规则。

5.根据权利要求4所述的一种基于空间特征向量滤波的城市生长元胞模拟方法，其特征在于，所述的步骤3-2具体为：

构建样本点的空间权重矩阵C：

其中，C_(i,j)为C中的第(i,j)项；d_(i,j)为样本点i与j之间的欧氏距离；r为覆盖样本点的最小生成树中最长的距离；

构建空间权重矩阵C的变形矩阵MCM并进行分解：

MCM＝EΛE^T

其中，M＝(I-11^T/n)，M为以某一变量为中心的投影矩阵，I和1分别为n*n代为矩阵和n*1单位向量，n为样本容量；E为提取得到的特征向量组成的矩阵，E＝(E₁,E₂,E₃,…,E_n)；Λ为以特征值λ＝(λ₁,λ₂,λ₃,…,λ_n)为对角线的n*n对角矩阵。

6.根据权利要求1所述的一种基于空间特征向量滤波的城市生长元胞模拟方法，其特征在于，所述的步骤4具体为：

步骤4-1：使用赤池信息准则AIC评估模型的拟合性能；

步骤4-2：使用残差的统计学指数检测引入空间特征滤波后土地利用转换规则的残差空间分布情况。

7.根据权利要求1所述的一种基于空间特征向量滤波的城市生长元胞模拟方法，其特征在于，所述的步骤5具体为：

利用优化的CA转换规则构建基于M-SVC的CA模型，并计算土地利用转换概率，具体方法如下：

任意元胞i在时刻t+1的状态由其在时刻t的状态和元胞邻域决定，则元胞状态的转换可表达为：

其中，f()是总体转换规则；Div_i代表驱动因子的影响；Nei_i代表邻域的影响；Res()代表限制因素；Count是元胞数量；

对于邻域影响的评估，通常采用Moorem×m或者圆形作为邻域，邻域的影响可表达为：

其中，Cot表示m×m领域范围内城市元胞的数量；

考虑了限制因素、邻域缩放、概率缩放后的总体转换概率可表示：

其中，TIP的取值范围为(0,0.1)，表示对转换概率的缩放；LAP的取值范围为(0.5,1.0)，表示对邻域影响的缩放；Res()的取值为0或1，Res()＝0表示不能转化为城市元胞，Res()＝1表示可以转化为城市元胞；

使用M-SVC构建的驱动因子对土地利用的影响P_Div可表达为：

γ～N(0,σ_λ ²Λ(α)),ε～N(0,σ²I)

β_k＝β_k,01+Eγ_k,

其中，y为因变量；x_k和β_k分别为自变量与对应系数向量；E为选取空间特征向量组成的矩阵；ε为随机噪声；

自变量x_k包含空间恒定部分与空间变化部分，对应回归系数向量β_k也包括空间恒定部分β_k,01及随空间变化部分Eγ_k；γ_k与γ均为满足正态分布的随机系数，其对应的矩阵Λ(α)与Λ(α_k)分别为特征值乘以比例常数α或α_k后组成的对称矩阵；

测算得到元胞全局转换概率P_overall后，与设定的阈值P_thd进行比较，以确定元胞i在下一时刻的状态是否发生转变：

其中，未知的阈值P_thd采用模拟逼近的方法进行探测，即先拟定一个阈值对目标土地进行判断，模拟出该阈值下的土地利用预测状况，若预测的城市元胞数大于目标年份的实际城市元胞数，则提高阈值；若预测的城市元胞数小于实际城市元胞数，则降低阈值，再进行判定模拟，直至预测元胞数与实际元胞数误差在1％以内，则可以确定该阈值及对应的模拟结果符合要求。

8.根据权利要求1所述的一种基于空间特征向量滤波的城市生长元胞模拟方法，其特征在于，所述的步骤6具体为：

步骤6-1：结合初始土地利用图、土地转换概率、邻域影响和限制因素，在UrbanCA模拟框架下，产生M-SVC-CA模拟；

步骤6-2：将距离待模拟时间点最近的真实土地利用作为初始年份，利用M-SVC-CA迭代N次，N＝终止年份-初始年份，模拟和预测土地利用格局。

9.根据权利要求1所述的一种基于空间特征向量滤波的城市生长元胞模拟方法，其特征在于，所述的步骤7具体为：

使用总体精度OA、Kappa模拟系数、图形优度FOM评估M-SVC-CA的模拟精度；

所述的总体精度OA通过比较终止年份的真实和模拟城市土地利用格局获取；

所述的图形优度FOM通过比较真实和模拟的初始-终止年份城市土地利用获取，将比较结果分解为初始城市Initial urban、实际和模拟均为非城市CR、实际和模拟均为城市Hit、实际为非城市模拟为城市False、实际为城市模拟为非城市Miss五类；

所述的Kappa模拟系数通过将模拟精度分解为数量一致性的K_transition和位置一致性K_transloc获取。

10.根据权利要求1所述的一种基于空间特征向量滤波的城市生长元胞模拟方法，其特征在于，所述的步骤9具体为：在既有的地理信息软件环境下输出并保存模拟结果。

Images