CN110991497B - 一种基于bsvc方法的城市土地利用变化模拟元胞自动机方法 - Google Patents

Abstract

本发明涉及一种基于BSVC方法的城市土地利用变化模拟元胞自动机方法。包括：1)对遥感影像进行监督分类，获得城市土地利用分类图；2)获取城市土地利用变化驱动因子数据，进行预处理；通过随机分层抽样方法获取土地利用图和驱动因子的有效样本点；3)利用贝叶斯空间变系数(BSVC)方法，对有效样本点进行训练获取元胞自动机的转化规则；4)根据BSVC已经建立的CA转化规则，获取城市土地利用转化概率；5)建立基于BSVC的CA模型(CA_BSVC)；6)利用CA_BSVC模型对城市土地利用进行模拟应用和验证分析，通过总体精度(OA)和图形优度(FOM)实现评估；7)输出并保存模拟结果。与现有技术相比，本发明具有更高的模拟精度和更好的城市土地利用变化模拟能力。

Description

一种基于BSVC方法的城市土地利用变化模拟元胞自动机方法

技术领域

本发明涉及一种城市土地利用变化模拟元胞自动机方法，尤其是涉及一种基于BSVC方法的城市土地利用变化模拟元胞自动机方法。

背景技术

城市增长反映了自然地表向城市人工区域的变化，对人类社会和生态环境产生了重大影响。在城市增长模拟中，一些模型考虑了空间的非平稳性，但现有模型方法不能充分反映城市动态的空间异质性特征。

地理元胞自动机(CA)是自底向上、自组织的城市模拟模型，用于城市土地利用历史模式重建和未来场景预测。为了定义基于驱动因子的转化规则，CA模型通常使用反映城市动态的生物物理、社会经济和基础设施的驱动因子因素。基于驱动因子的转化规则通常使用逻辑回归(LR)来构建，它假设每个因素对城市空间增长的影响是相同的。这意味着逻辑回归推导的参数具有空间不变性，不能反映城市增长的空间非平稳性。为了解决空间非平稳性问题，一些学者发展了一些空间方法来构造CA模型。例如，基于分区的CA建模将空间划分为几个子区域，然后使用LR来处理城市增长的空间差异。采用K-means和KNN-cluster根据进化速率进行区域划分，其性能优于传统CA模型。通过在每个子区域采用唯一的转化规则，分区方法可以显式和隐式地反映空间非平稳。这些方法虽然可以表示空间的非平稳性，但由于使用了空间不变的参数，不能充分表达城市增长的空间非平稳性。

空间变化系数模型(SVCM)通过允许系数随空间和位置的变化而变化，可以生成具有空间非平稳性的转化规则，从而可以明确城市增长的空间非平稳性。在SVCM中，地理加权回归(GWR)是一种典型的方法，在没有预先定义区域分类的情况下，将空间非平稳性纳入参数化。GWR是一种局部统计模型，它使用样本的距离加权子集为每个点生成回归系数，这些系数对样本的位置很敏感。近期技术，GWR模型与CA模型的结合能够较传统CA模型更准确地模拟城市扩张。虽然该方法可以考虑空间的非平稳性，但仍存在一些可能导致建模困难的特点。首先，由于GWR回归依赖于周围样本的属性，数据集异常值可能通过空间飞地效应影响回归系数。其次，非恒定方差在空间中广泛存在，可能无法满足正态误差的假设。相比之下，贝叶斯空间变化系数(BSVC)方法可以检测非恒定方差，检测异常值的影响，从而缓解系数估计的偏差。因此，BSVC方法应更适合与CA模型结合进行城市增长模拟。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种基于BSVC方法的城市土地利用变化模拟元胞自动机方法。

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

一种基于BSVC方法的城市土地利用变化模拟元胞自动机方法，该方法包括以下步骤：

步骤1：对遥感影像进行监督分类获得初始和结束年份的土地利用分类图；

步骤2：获取城市土地利用变化驱动因子数据，经预处理后获取土地利用分类图和驱动因子数据的有效样本点；

步骤3：利用贝叶斯空间变系数方法，对有效样本点进行训练获取元胞自动机的转化规则；

步骤4：利用根据BSVC训练建立的元胞自动机的转化规则，获取城市土地利用转化概率；

步骤5：综合转化概率、元胞领域、随机因素和限制因素建立基于BSVC的CA模型，即CA_BSVC模型；

步骤6：利用CA_BSVC模型对城市土地利用进行模拟应用和验证分析，得到CA_BSVC模型模拟结果；

步骤7：对CA_BSVC模型及其模拟结果，分别从规则拟合和模拟结果两方面进行精度评定，并输出保存模拟结果。

进一步地，所述的步骤1包括以下分步骤：

步骤11：获取建模所需的卫星遥感影像和矢量地图数据，对其进行空间参考统一和几何校正；

步骤12：利用两期卫星遥感图像，基于光谱角监督分类方法获取初始年份和结束年份的城市的土地利用分类图。

进一步地，所述的步骤2包括以下分步骤：

步骤21：选取影响城市土地利用变化的驱动空间因子数据，所述驱动空间因子数据包括到高速公路、铁路、地铁和一级公路及高程、人口数据；

步骤22：通过遥感图像数据、行政区划图、道路交通图在ArcGIS中利用欧几里德距离获取到高速公路、铁路、地铁和一级公路距离的变量；

步骤23：利用随机分层抽样方法对土地利用分类图和因子图层进行采样，获取用于为CA规则转化提供可靠样本数据点的土地利用分类图和驱动因子数据的有效样本点。

进一步地，所述的步骤4具体包括：利用根据BSVC训练建立的元胞自动机的转化规则在设定空间分辨率下获取空间变量影响下土地的转化概率分布，获取土地的转化概率的方法包括：

假定s表示元胞状态是否转变，从时间t到t+1元胞状态从非城市转变为城市，则y记为1；从时间t到t+1元胞状态未发生改变，则s记为0；

利用获取的空间变量数据进行土地的转化概率的测算。

进一步地，所述的步骤5中CA_BSVC模型的核心问题为确定下一步是否将一个单元从一种状态转换到另一种状态，对应描述公式为：

式中，

表示t+1时刻的元胞状态，

表示t时刻的元胞状态，f表示计算总转移概率P_g的综合转移规则，P_d表示基于驱动因子的土地利用转化概率，

表示领域的影响，Con表示空间抑制函数。

进一步地，所述的步骤5中CA_BSVC模型的基于驱动因子的土地利用转化概率和总转换概率，其计算公式为：

式中，a_i表示第i个驱动因子的权重，x_i表示第i个驱动因子，ε表示拟合残差，TIP表示时间增量参数，LAP表示局部调整参数。

进一步地，所述的步骤5中CA_BSVC模型的领域采用Moore领域，其描述公式为：

式中，

表示m×m领域范围内的城市细胞的总数，(j≠i)表示中心细胞i不参与计算。

进一步地，所述的第i个驱动因子的权重，其计算公式为：

式中，W_ij表示行标准化矩阵，I_k表示k行k列的标准矩阵，u_i表示随机扰动，β_i表示位置i处的模型系数矩阵，σ²表示方差，δ²表示控制平滑影响的比例因子，W_i表示位置i中的(n×n)空间矩阵，X表示依赖变量矩阵(n×k)。

进一步地，所述的步骤6具体包括：利用UrbanCA软件和MTALAB语言实现CA_BSVC模型的模拟与预测流程，选用某年的土地利用格局为初始状态利用CA模型运行M次，其中M表示初始与结束的年份差，得到土地利用变化的模拟及预测结果。

进一步地，所述的步骤7包括以下分步骤：

步骤71：通过与遥感分类的土地利用格局进行比较，对CA_BSVC模型模拟结果进行精度计算并评定，所述精度计算的指标包括：图形优度FOM和总体精度OA；

步骤72：将CA_CEO模型模拟结果与遥感分类结果进行叠加并评定，所叠加的结果包括：实际和模拟均为城市Hit、实际为非城市模拟为城市False、实际为城市模拟为非城市Miss、实际和模拟均为非城市CR及水体Water；

步骤73：在GIS软件中输出并保存模拟结果。

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

(1)本发明考虑贝叶斯空间变系数，构建了基于贝叶斯空间变系数(BSVC)的CA模型CA_BSVC，并应用于城市土地利用变化的模拟。通过将该模型与地理加权回归(GWR)的CA模型CA_GWR的模拟结果进行对比。CA_BSVC在总体精度和图形优度方面均优于CA_GWR。基于元胞自动机与贝叶斯空间变系数的城市土地利用变化模型模拟的方法能够较好的模拟和预测城市土地利用变化动态。

(2)本发明涉及一种基于BSVC方法的城市土地利用变化模拟元胞自动机方法，采用的影响因素主要为空间距离变量与自然要素及社会经济要素变量，且采用的贝叶斯空间变系数CA模型(CA_BSVC)模型的模拟效果优于地理加权回归CA模型(CA_GWR)，能够较好的完成对城市土地利用变化的模拟。

附图说明

图1为本发明的流程图；

图2为实施例案例研究区域图；

图3为城市土地利用驱动因子图；

图4为真实城市土地利用变化和CA_BSVC、CA_GWR两种模型的转化概率图；

图5为真实城市土地分布和CA_BSVC、CA_GWR两种模型模拟的城市土地分布图。

具体实施方式

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

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

如图1所示为本发明基于贝叶斯空间变系数(BSVC)方法的城市土地利用变化模拟元胞自动机方法的整体流程，包括以下步骤：

1)利用光谱角方法对遥感影像进行监督分类，获得初始和结束年份的城市土地利用分类图；

2)获取城市土地利用变化驱动因子数据，并进行预处理；通过随机分层抽样方法获取土地利用图和驱动因子的有效样本点；

3)利用贝叶斯空间变系数(BSVC)方法，对有效样本点进行训练获取元胞自动机(CA)的转化规则；

4)根据BSVC已经建立的CA转化规则，获取城市土地利用转化概率；

5)综合转化概率、元胞邻域、随机因素、限制因素建立基于BSVC的CA模型(CA_BSVC)；

6)利用CA_BSVC模型对城市土地利用进行模拟应用和验证分析，通过总体精度(OA)和图形优度(FOM)实现评估；

7)输出并保存模拟结果。

所述步骤1)具体为：

8)获取建模所需的卫星遥感影像和矢量地图数据，对其进行空间参考统一和几何校正；

9)利用两期卫星遥感图像，基于光谱角监督分类方法获取初始年份和结束年份的城市土地利用分类图；

所述步骤2)具体为：

10)首先选取影响城市土地利用变化的空间因子，包括到高速公路、铁路、地铁和一级公路及高程、人口数据；

11)通过遥感图像数据、行政区划图、道路交通图在ArcGIS中利用欧几里德距离获取到高速公路、铁路、地铁和一级公路距离的变量。

所述步骤3)具体为：

12)首先利用随机分层抽样方法对土地利用图、因子图层进行采样，为CA规则转化提供可靠的样本点数据；

所述步骤4)具体为；

13)利用贝叶斯空间变系数(BSVC)，对获取的有效采样点数据训练CA转化规则；

14)利用已建立的CA转化规则，在30m空间分辨率下获取空间变量影响下的转化概率分布。获取土地转变概率数据方法如下：

假定s表示元胞状态是否转变，从时间t到t+1元胞状态从非城市(Non-urban)转变为城市(Urban)，s记为1；从时间t到t+1元胞状态未发生改变，则s记为0。利用获取的空间变量数据进行土地转变概率的测算。

所述步骤5)具体为；

15)基于步骤14)的数据建立基于BSVC和GWR的地理CA模型，具体如下：CA模型确定下一步是否将一个单元从一种状态转化到另一种状态，t+1时刻的元胞状态可以表示为：

式中，

表示t+1时刻的元胞状态，

表示t时刻的元胞状态，f表示计算总转移概率P_g的综合转移规则，P_d表示基于驱动因子的土地利用转化概率，

表示领域的影响，Con表示空间抑制函数。

利用时间增量参数(TIP)来抵消基于驱动因子的转化概率的衰减、局部调整参数(LAP)来抵消邻域影响的增加。因此，基于驱动因子的转化概率和总转移概率可由以下公式给出：

式中，a_i表示第i个驱动因子的权重，x_i表示第i个驱动因子，ε表示拟合残差，TIP表示时间增量参数，其取值范围为0.0～0.1，LAP表示局部调整参数，其取值范围为0.5～1。

对于邻域影响的评估，CA模型较多地采用正方形或圆形的规则邻域，如m×m的Moore邻域可以表达为：

式中，

表示m×m领域范围内的城市细胞的总数，(j≠i)表示中心细胞i不参与计算，本发明选取Moore5×5作为元胞邻域。

限制因素Con表示元胞受到某种限制而不能发展并转化为城市元胞，这种限制包括大面积水体、基本农田、生态保护区、公园、绿地等。Con可以表达为：

Con＝Bin(cell_i(t)～available)

式中，Con取值0或1，0表示该元胞不可以发展为城市元胞，1表示可以发展为城市元胞。

随机因素R用于模拟不确定因素导致的元胞状态转变，如某元胞在附近没有城市元胞的状态下通过随机因素提高发展概率，从非城市状态转化为城市状态。随机因素R表达为：

R＝1+(lnr)^a

式中，r表示0到1之间的随机数，a表示随机因素R的控制参数，取值为0和10之间的整数。

由空间变量决定的基于驱动因子的土地利用转化概率P_d是转换规则的核心部分，它代表了这些因素对土地利用的影响，并通过概率的方式影响到下一时刻的元胞状态。如果利用地理加权回归(GWR)去获取CA参数，第i个驱动因子的权重a_i内可通过下式计算：

W_iY＝W_iXβ_i+ε_i

式中，W_i表示位置i中的(n×n)空间矩阵；Y是解释变量向量(n×1)；X是依赖变量矩阵(n×k)；β_i表示位置i处的模型系数矩阵；ε_i表示随机误差。

如果利用贝叶斯空间变系数(BSVC)获取CA参数，第i个驱动因子的权重a_i可修改为：

式中，W_ij表示行标准化矩阵，I_k表示k行k列的标准矩阵，u_i表示随机扰动，β_i表示位置i处的模型系数矩阵，σ²表示方差，δ²表示控制平滑影响的比例因子，W_i表示位置i中的(n×n)空间矩阵，X表示依赖变量矩阵(n×k)。

根据上面公式计算得出土地利用的总转换概率P_g。在实际计算中，利用MTALAB语言完成GWR与BSVC的参数计算，并将计算结果与设定的阈值P_thd进行比较，判断元胞是否会在一个时刻发生转化。当元胞i的转换概率P_gi,t大于设定阈值P_thd时转化为城市类型，否则元胞保持状态不变，对应描述公式为：

所述步骤6)具体为；

16)利用UrbanCA软件和MTALAB语言实现CA_BSVC的模拟与预测流程，选用某年的土地利用格局为初始状态利用CA模型运行M次(初始与结束的年份差)得到土地利用变化的模拟及预测结果。

17)将基于两种模型模拟以及预测的土地利用变化结果输出。

所述步骤7)具体为；

18)对CA_BSVC模型模拟的土地利用结果分别从规则拟合精度及模拟结果两方面进行模拟精度的评定

首先通过与遥感分类的土地利用格局进行比较，对模拟结果进行精度计算，主要指标有图形优度(FOM)和总体精度(OA)。将总体精度分解为城市(Hit)和非城市(CR)两类、将错误分解为忽略性(Misses)和替代性(False)两类，忽略性错误指的是实际为城市但模拟为非城市，即CA模型未能捕捉到的城市元胞；替代性错误指的是实际为非城市但模拟为城市，即CA模型错误地增加的城市元胞。

将模拟结果与遥感分类结果进行叠加，叠加结果包含5类：实际和模拟均为城市(Hit)、实际为非城市模拟为城市(False)、实际为城市模拟为非城市(Miss)、实际和模拟均为非城市(CR)、以及水体(Water)。根据目视判别显示，对比CA_BSVC模型的模拟结果与实际分类结果的差异。

所述步骤8)具体为：

19)在GIS软件中输出并保存模拟结果。

本发明的实际具体实施例如下：

以2008-2018年间合肥城市土地利用为案例，本案例区域位置如图2所示。为验证CA_BSVC模型在土地利用变化模拟中的有效性，案例中将基于地理加权回归的CA模型(CA_GWR)作为比较对象，模拟了同期城市土地利用变化过程，结果表明CA_BSVC的模拟效果优于CA_GWR模型。基于贝叶斯空间变系数的城市土地利用变化模拟元胞自动机方法，包括以下步骤：

1)首先选择合肥市2008、2018年的遥感图像数据，以及行政区划图、道路交通图，作为训练CA规则转化以及获取土地转变概率的基础数据；

2)根据遥感图像数据利用随机分层抽样方法对各空间变量的值、土地利用初始年份和结束年份状态值进行采样处理；

3)利用各年遥感图像以及行政区划图层、道路交通图层，然后利用欧几里德距离计算出到高速公路、铁路、地铁和一级公路的距离及高程、人口数据(表1)；

4)在利用利用光谱角方法对合肥遥感图像进行监督分类，从而对土地利用格局进行解译；

5)利用随机分层抽样方法得到的有效样本点以及各空间变量值，在实际计算中，利用MATLAB语言进行贝叶斯空间变系数(BSVC)和地理加权回归(GWR)的实现。表1显示了两种模型的拟合优度，表明CA_BSVC模型更好；表2显示了两种模型下的转化规则参数，图3显示了两种模型下的城市土地转变；

6)利用获取的土地转变概率和CA转化规则(图4)，建立基于BSVC和GWR的地理CA模型CA_BSVC和CA_GWR；

7)以2008年状态为初始值，分别利用CA_BSVC和CA_GWR模型运行10次进行，从而预测2018年土地利用变化(图5)；

8)将模拟预测后的结果与真实城市增长进行对比，并分析总体精度(OA)和图形优度(FOM)的变化；表2显示，相比基于GWR的CA模型，本发明的总体精度提高0.5％、而关注变化的图形优度提高23.5％；

9)输出并保存可视化的结果。

表1 CA_GWR和CA_BSVC模型参数综合统计

表2 CA_GWR和CA_BSVC的精度结果分析

注意：变化＝(BSVC-GWR)/GWR×100

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

Claims (4)

1.一种基于BSVC方法的城市土地利用变化模拟元胞自动机方法，其特征在于，该方法包括以下步骤：

步骤1：对遥感影像进行监督分类获得初始和结束年份的土地利用分类图；

步骤2：获取城市土地利用变化驱动因子数据，经预处理后获取土地利用分类图和驱动因子数据的有效样本点；

步骤3：利用贝叶斯空间变系数方法，对有效样本点进行训练获取元胞自动机的转化规则；

步骤4：利用根据BSVC训练建立的元胞自动机的转化规则，获取城市土地利用转化概率；

步骤5：综合转化概率、元胞领域、随机因素和限制因素建立基于BSVC的CA模型，即CA_BSVC模型；

步骤6：利用CA_BSVC模型对城市土地利用进行模拟应用和验证分析，得到CA_BSVC模型模拟结果；

步骤7：对CA_BSVC模型及其模拟结果，分别从规则拟合和模拟结果两方面进行精度评定，并输出保存模拟结果；

所述的步骤2包括以下分步骤：

步骤21：选取影响城市土地利用变化的驱动空间因子数据，所述驱动空间因子数据包括到高速公路、铁路、地铁和一级公路及高程、人口数据；

步骤22：通过遥感图像数据、行政区划图、道路交通图在ArcGIS中利用欧几里德距离获取到高速公路、铁路、地铁和一级公路距离的变量；

步骤23：利用随机分层抽样方法对土地利用分类图和因子图层进行采样，获取用于为CA规则转化提供可靠样本数据点的土地利用分类图和驱动因子数据的有效样本点；

所述的步骤4具体包括：利用根据BSVC训练建立的元胞自动机的转化规则在设定空间分辨率下获取空间变量影响下土地的转化概率分布，获取土地的转化概率的方法包括：

假定s表示元胞状态是否转变，从时间t到t+1元胞状态从非城市转变为城市，则y记为1；从时间t到t+1元胞状态未发生改变，则s记为0；

利用获取的空间变量数据进行土地的转化概率的测算；

所述的步骤5中CA_BSVC模型的核心问题为确定下一步是否将一个单元从一种状态转换到另一种状态，对应描述公式为：

式中，

表示t+1时刻的元胞状态，

表示t时刻的元胞状态，f表示计算总转移概率P_g的综合转移规则，P_d表示基于驱动因子的土地利用转化概率，Ω_it表示领域的影响，Con表示空间抑制函数；

所述的步骤5中CA_BSVC模型的基于驱动因子的土地利用转化概率和总转换概率，其计算公式为：

式中，a_i表示第i个驱动因子的权重，x_i表示第i个驱动因子，ε表示拟合残差，TIP表示时间增量参数，LAP表示局部调整参数；

所述的步骤5中CA_BSVC模型的领域采用Moore领域，其描述公式为：

式中，

表示m×m领域范围内的城市细胞的总数，(j≠i)表示中心细胞i不参与计算；

所述的第i个驱动因子的权重，其计算公式为：

式中，W_ij表示行标准化矩阵，I_k表示k行k列的标准矩阵，u_i表示随机扰动，β_i表示位置i处的模型系数矩阵，σ²表示方差，δ²表示控制平滑影响的比例因子，W_i表示位置i中的(n×n)空间矩阵，X表示依赖变量矩阵(n×k)。

2.根据权利要求1所述的一种基于BSVC方法的城市土地利用变化模拟元胞自动机方法，其特征在于，所述的步骤1包括以下分步骤：

步骤11：获取建模所需的卫星遥感影像和矢量地图数据，对其进行空间参考统一和几何校正；

步骤12：利用两期卫星遥感图像，基于光谱角监督分类方法获取初始年份和结束年份的城市的土地利用分类图。

3.根据权利要求1所述的一种基于BSVC方法的城市土地利用变化模拟元胞自动机方法，其特征在于，所述的步骤6具体包括：利用UrbanCA软件和MTALAB语言实现CA_BSVC模型的模拟与预测流程，选用某年的土地利用格局为初始状态利用CA模型运行M次，其中M表示初始与结束的年份差，得到土地利用变化的模拟及预测结果。

4.根据权利要求1所述的一种基于BSVC方法的城市土地利用变化模拟元胞自动机方法，其特征在于，所述的步骤7包括以下分步骤：

步骤71：通过与遥感分类的土地利用格局进行比较，对CA_BSVC模型模拟结果进行精度计算并评定，所述精度计算的指标包括：图形优度FOM和总体精度OA；

步骤72：将CA_CEO模型模拟结果与遥感分类结果进行叠加并评定，所叠加的结果包括：实际和模拟均为城市Hit、实际为非城市模拟为城市False、实际为城市模拟为非城市Miss、实际和模拟均为非城市CR及水体Water；

步骤73：在GIS软件中输出并保存模拟结果。

Images