CN106021751A - 基于ca和sar的海岸带土地利用变化模拟方法 - Google Patents

Abstract

本发明涉及一种基于CA和SAR的海岸带土地利用变化模拟方法，包括以下步骤：1)对土地利用初始年份和结束年份进行遥感分类，基于矢量地图数据和卫星遥感影像获得影响土地利用转化的各空间变量因子数据；2)获取土地利用初始年份状态值、结束年份状态值以及各因子样本数据；3)对有效样本点建立CA规则转换，产生空间变量影响下土地的转化概率图；4)建立基于SAR的地理SARCA模型；5)利用已建立的SARCA模型进行土地利用变化的模拟和预测，并输出土地利用变化的模拟结果；6)对SARCA分别从规则拟合精度及模拟结果两方面进行模拟精度的评定；7)输出并保存模拟结果。与现有技术相比，本发明具有模拟效果好等优点。

Description

基于CA和SAR的海岸带土地利用变化模拟方法

技术领域

本发明涉及一种海岸带土地利用变化模拟方法，尤其是涉及一种基于CA(元胞自动机)与SAR(空间自回归)的海岸带土地利用变化模拟方法。

背景技术

对土地利用变化模拟的方法一般如下：

元胞自动机(CA)是一种能够“自下而上”自动演化的模型，常用于模拟地理复杂系统及其动态演变过程，在城市扩展和土地利用变化领域尤其受到重视。国内外文献从转换规则、元胞邻域、尺度影响、模型评价、案例应用、精度评定等多角度进行研究。元胞下一时刻的状态由转换规则决定的，因此转换规则是CA建模的核心和难点、地理系统模拟科学问题之所在，该规则直接决定了CA模型的特征并影响其模拟性能。

目前有很多方法可以用来获取转换规则，简单地归纳为两类：显性规则和隐性规则。显性规则指的是CA参数可以类比为空间变量的权重或规则具有明确的物理意义，通过这些权重能够清晰了解每一个变量对土地利用变化的贡献，Logistic回归是其典型代表；隐性规则指的是构建CA规则的参数没有明确的物理意义，神经网络是其典型代表。除了Logistic回归，显性规则还包括主成分分析、偏最小二乘回归、核方法、粗集、模糊集和案例推理等，也包括全局最优化方法，如遗传算法、粒子群算法和模拟退火算法等。除了神经网络，隐性规则还包括蚁群智能、免疫智能、支持向量机、系统动力学、CA-Markov等。国内外文献中涉及的CA模型远不止上述所列，这些模型各有其优势与适用性，成功地应用在国内外不同区域，实质性地推动了CA理论与应用的发展。

空间实体之间存在空间自相关特性，土地利用与城市格局也不例外。地理学第一定律表明，相邻的地理实体有较高的相似性，距离较远的实体相似性较小。这意味着土地利用在空间上可能存在聚类簇，在这些聚类簇内土地利用与空间变量之间的关系是高度相似的。因此，构建CA模型时如果考虑空间自相关因素，所建立的模型则能很好地体现相关性及距离因素的影响。空间自回归(Spatial AutoRegression，SAR)是一种考虑空间自相关特性的统计学方法，在土地利用变化分析中已经得到应用，将其用于CA规则的获取能够使建立的模型融合空间自相关特性。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种模拟效果好的基于CA和SAR的海岸带土地利用变化模拟方法，采用的影响因素主要为空间距离变量，且采用的空间自回归的地理CA模型SARCA的模拟效果优于传统的LogCA模型，且能较好的完成对海岸带土地利用变化的模拟。

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

一种基于CA和SAR的海岸带土地利用变化模拟方法，包括以下步骤：

1)对土地利用初始年份和结束年份进行遥感分类，基于矢量地图数据和卫星遥感影像获得影响土地利用转化的各空间变量因子数据；

2)通过采样获取研究区域有效样本点，获取土地利用初始年份状态值、结束年份状态值以及各因子样本数据；

3)根据步骤1)和2)的数据对有效样本点建立CA规则转换，产生空间变量影响下土地的转化概率图；

4)利用已建立的转换规则，建立基于SAR的地理SARCA模型；

5)利用已建立的SARCA模型进行土地利用变化的模拟和预测，并输出土地利用变化的模拟结果；

6)对SARCA分别从规则拟合精度及模拟结果两方面进行模拟精度的评定；

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

所述的步骤1)具体为：

11)利用两期遥感图像，基于马氏距离监督分类方法获取初始年份和结束年份的遥感分类土地利用图；

12)选取影响土地利用变化的空间变量因子，包括到市中心、区中心、主干道、农用地、中环线、外环线的距离；

13)通过遥感图像数据、行政区划图、道路交通图在ArcGIS中利用欧几里德距离获取到市中心、区中心、主干道、农用地、中环线、外环线距离的变量。

所述的步骤2)具体为：

14)设定采样分辨率，对土地利用图、因子图层进行采样，通过ArcGIS生成Fishnet采样点，获取土地利用初始年份状态值、结束年份状态值以及各因子样本数据，为CA规则转换提供可靠数据。

所述的步骤3)具体为；

15)利用空间自回归方法SAR，对获取的有效采样点数据以及空间变量因子数据训练CA转换规则；

16)利用已建立的CA转换规则，在30m空间分辨率下获取空间变量影响下的土地转化概率P_di分布。

所述的获取土地转化概率数据方法如下：

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

所述的步骤4)具体为；

17)基于步骤16)的数据建立基于SAR的地理CA模型：

任意元胞i在时刻t+1的状态由其在时刻t的状态和元胞邻域决定，则非城市Non-urban元胞转化为城市Urban元胞的全局概率表达为P_gi，t，可表达为：

P_gi，t＝P_di×P_Ni，t×Con×R (1)

式中，P_di是由空间变量决定的土地利用转化概率，P_Ni，t是半径为N的邻域影响，Con是限制因素，R是随机因素；

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

$P_{Ni,t}=\frac{{\mathrm{\Σ}}_j^{m\×m}\left(S_{i,t=Urban}\right)(j\≠i)}{m\×m-1}---\left(2\right)$

式中，中心元胞i不参与计算，表示m×m邻域范围内城市元胞的总数，选取Moore 5×5作为元胞邻域；

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

Con＝Bin(cell_i(t)～available) (3)

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

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

R＝1+(lnr)^a (4)

式中，r是0到1之间的随机数，α是随机因素R的控制参数，取值是0和10之间的整数；

由空间变量决定的转化概率P_dt是转换规则的核心部分，代表了这些因素对土地利用的影响，并通过概率的方式影响到下一时刻的元胞状态，概率P_dt表达为：

$P_{di}=\frac{1}{1+\exp (-z_i)}---\left(5\right)$

式中，z_i是空间变量对土地利用转变的影响，z_i通过logistic回归获取，z_i表达为：

$z_i=a_0+{\mathrm{\Σ}}_{j=1}^k{a}_j{x}_j---\left(6\right)$

式中，a₀是常数，a_j是CA转换规则的参数，亦即空间变量x_j的权重，其中j＝1，2，...，k；

如果利用空间自回归(SAR)去获取CA参数，则z_i表达为：

$\left\{\begin{array}{c}{z}_i={\mathrm{\ρWz}}_i+X\mathrm{\β}+\mathrm{\ϵ}\\ \mathrm{\ϵ}~N(0,{\mathrm{\σ}}^2{I}_n)\end{array}\right.---\left(7\right)$

式中，W是z_i的标准化空间权重矩阵，ρ是空间滞后Wz_i的参数，X＝(x_i，..，x_k)是向量形式的空间变量，β是空间变量X的参数向量，即空间变量的权重向量，其意义等同于式(6)中的a_j，σ²和I_n分别是残差项ε的方差和空间自相关系数；基于空间邻接关系定义空间权重W，具体地采用1阶Q邻接(Queen Contiguity)权重矩阵；

通过式(7)得到式(1)等式右侧的4项，即可测算元胞i的全局转化概率P_gi，t；在实际计算中，logistic回归用R语言(RGui)实现，SAR则用GeoDa实现。

计算得到元胞全局转换概率P_gi，i之后，与设定的阈值P_thd(其值在0～1之间)进行比较，以确定元胞i在下一时刻的状态是否发生转变，当元胞i的转换概率P_gi，t大于设定阈值P_thd时转化为城市类型，否则元胞保持状态不变：

$S_{i,t+1}=\left\{\begin{array}{cc}Urban& if{P}_{gi,t}\&GreaterEqual;P_{thd}\\ Nochange& if{P}_{gi,t}<P_{thd}\end{array}\right.---\left(8\right)$

所述的步骤5)具体为；

18)采用GeoDa实现SAR，选用某年的土地利用格局为初始状态利用SARCA模型运行M次(初始与结束的年份差)得到土地利用变化的模拟及预测结果；

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

所述的步骤6)具体为；

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

首先通过与遥感分类的土地利用格局进行比较，对模拟结果进行精度计算，主要指标有总体精度(Overall Accuracy)、Kappa系数和景观指数(Landscape Metrics)等。将总体精度分解为城市(Hit)和非城市(Correct Rejection)两类、将错误分解为忽略性(Miss或Omission)和替代性(False或Commission)两类，利用Pontius方法将错误同时分解为数量错误(Quantity Error)和位置错误(Allocation Error)。忽略性错误指的是实际为城市但模拟为非城市，即CA模型未能捕捉到的城市元胞；替代性错误指的是实际为非城市但模拟为城市，即CA模型错误地增加的城市元胞；

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

与现有技术相比，本发明能够利用空间变量因子数据以及CA规则转换建立SARCA模型对土地利用变化进行模拟结果预测，同时为检测SARCA模型的性能，利用Logistic回归建立CA模型(LogCA)并模拟同时期土地利用变化，在规则拟合精度及模拟结果两方面与SARCA模型进行比较。CA规则的拟合结果表明，SAR在AIC、残差的描述性统计量和空间自相关指标等方面均优于Logistic回归。因此SARCA模型在模拟以及预测土地利用变化结果方面具有较高的合理性。

附图说明

图1为本发明的流程图；

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

图3为空间变量决定的土地利用转化概率P_di示意图；

图4为LogCA和SARCA两种模型下2015年的预测结果示意图。

具体实施方式

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

以2005-2015年间上海城市土地利用为案例区域，本案例区域位置如图2所示。为验证SARCA模型在土地利用变化模拟上的有效性，案例中将基于Logistic回归的CA模型(LogCA)作为比较对象，也模拟了上海市同期城市土地利用变化过程，结果表明SARCA的模拟效果优于传统LogCA模型。基于元胞自动机与空间自回归的海岸带土地利用变化模拟模型方法，包括以下步骤：

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

2)根据遥感图像数据利用系统抽样方法对每1km²网络抽取一个点，获取各空间变量的值、土地利用初始年份和结束年份状态值，上海全域共获有6701个有效样本点；

3)利用各年遥感图像以及行政区划图层、道路交通图层、然后利用空间分析工具中的欧几里德距离(Euclidean Distance)计算出各区域到市中心、到区中心、到主干道路、到农用地、到中环线以及到外环线的距离数据；

4)在ENVI软件中利用监督分类中的最大似然分类法，对上海遥感图像进行监督分类从而对土地利用格局进行解译；

5)利用系统抽样方法得到的有效样本点以及各空间变量值，在实际计算中，用R语言进行Logistic回归的实现，用GeoDa来实现SAR分别获取土地转变概率。表1显示了两种模型的拟合优度，表明SARCA模型更好，表2显示了两种模型下的转换规则参数，图3显示了两种模型下的土地转变；

6)利用获取的土地转变概率以及CA转换规则建立基于SAR的地理CA模型和LogCA模型；

7)以2005年状态为初始值分别利用SARCA模型以及LogCA模型运行10次进行从而预测2015年土地利用变化；

8)将模拟预测后的土地利用变化值在ArcGIS中加载，得到2015年图4的土地利用变化图，同时模拟结果的精度显示SARCA显著优于传统LogCA模型；

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

表1 CA转换规则获取的拟合效果

表2 用于测算土地转变概率的空间变量

表3 2015年两种模型的预测精度

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

Claims (8)

1.一种基于CA和SAR的海岸带土地利用变化模拟方法，其特征在于，包括以下步骤：

1)对土地利用初始年份和结束年份进行遥感分类，基于矢量地图数据和卫星遥感影像获得影响土地利用转化的各空间变量因子数据；

2)通过采样获取研究区域有效样本点，获取土地利用初始年份状态值、结束年份状态值以及各因子样本数据；

3)根据步骤1)和2)的数据对有效样本点建立CA规则转换，产生空间变量影响下土地的转化概率图；

4)利用已建立的转换规则，建立基于SAR的地理SARCA模型；

5)利用已建立的SARCA模型进行土地利用变化的模拟和预测，并输出土地利用变化的模拟结果；

6)对SARCA分别从规则拟合精度及模拟结果两方面进行模拟精度的评定；

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

2.根据权利要求1所述的一种基于CA和SAR的海岸带土地利用变化模拟方法，其特征在于，所述的步骤1)具体为：

11)利用两期遥感图像，基于马氏距离监督分类方法获取初始年份和结束年份的遥感分类土地利用图；

12)选取影响土地利用变化的空间变量因子，包括到市中心、区中心、主干道、农用地、中环线、外环线的距离；

13)通过遥感图像数据、行政区划图、道路交通图在ArcGIS中利用欧几里德距离获取到市中心、区中心、主干道、农用地、中环线、外环线距离的变量。

3.根据权利要求2所述的一种基于CA和SAR的海岸带土地利用变化模拟方法，其特征在于，所述的步骤2)具体为：

14)设定采样分辨率，对土地利用图、因子图层进行采样，获取土地利用初始年份状态值、结束年份状态值以及各因子样本数据，为CA规则转换提供可靠数据。

4.根据权利要求3所述的一种基于CA和SAR的海岸带土地利用变化模拟方法，其特征在于，所述的步骤3)具体为；

15)利用空间自回归方法SAR，对获取的有效采样点数据以及空间变量因子数据训练CA转换规则；

16)利用已建立的CA转换规则，在30m空间分辨率下获取空间变量影响下的土地转化概率P_di分布。

5.根据权利要求4所述的一种基于CA和SAR的海岸带土地利用变化模拟方法，其特征在于，所述的获取土地转化概率数据方法如下：

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

6.根据权利要求4所述的一种基于CA和SAR的海岸带土地利用变化模拟方法，其特征在于，所述的步骤4)具体为；

17)基于步骤16)的数据建立基于SAR的地理CA模型：

任意元胞i在时刻t+1的状态由其在时刻t的状态和元胞邻域决定，则非城市Non-urban元胞转化为城市Urban元胞的全局概率表达为P_gi，t，可表达为：

P_gi，t＝P_di×P_Ni，t×Con×R (1)

式中，P_di是由空间变量决定的土地利用转化概率，P_Ni，t是半径为N的邻域影响，Con是限制因素，R是随机因素；

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

$P_{Ni,t}=\frac{{\mathrm{\&Sigma;}}_j^{m\&times;m}\left(S_{i,t=Urban}\right)(j\&NotEqual;i)}{m\&times;m-1}---\left(2\right)$

式中，中心元胞i不参与计算，表示m×m邻域范围内城市元胞的总数，选取Moore 5×5作为元胞邻域；

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

Con＝Bin(cell_i(t)～available) (3)

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

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

R＝1+(lnr)^a (4)

式中，r是0到1之间的随机数，α是随机因素R的控制参数，取值是0和10之间的整数；

由空间变量决定的转化概率P_dt是转换规则的核心部分，代表了这些因素对土地利用的影响，并通过概率的方式影响到下一时刻的元胞状态，概率P_dt表达为：

$P_{di}=\frac{1}{1+\exp (-z_i)}---\left(5\right)$

式中，z_i是空间变量对土地利用转变的影响，z_i通过logistic回归获取，z_i表达为：

$z_i=a_0+{\mathrm{\&Sigma;}}_{j=1}^k{a}_j{x}_j---\left(6\right)$

式中，a₀是常数，a_j是CA转换规则的参数，亦即空间变量x_j的权重，其中j＝1，2，...，k；

如果利用空间自回归(SAR)去获取CA参数，则z_i表达为：

$\left\{\begin{array}{c}{z}_i={\mathrm{\&rho;Wz}}_i+X\mathrm{\&beta;}+\mathrm{\&epsiv;}\\ \mathrm{\&epsiv;}~N(0,{\mathrm{\&sigma;}}^2{I}_n)\end{array}\right.---\left(7\right)$

式中，W是z_i的标准化空间权重矩阵，ρ是空间滞后Wz_i的参数，X＝(x_i，...，x_k)是向量形式的空间变量，β是空间变量X的参数向量，即空间变量的权重向量，其意义等同于式(6)中的a_j，σ²和I_n分别是残差项ε的方差和空间自相关系数；基于空间邻接关系定义空间权重W；

通过式(7)得到式(1)等式右侧的4项，即可测算元胞i的全局转化概率P_gi，t；

计算得到元胞全局转换概率P_gi，t之后，与设定的阈值P_thd进行比较，以确定元胞i在下一时刻的状态是否发生转变，当元胞i的转换概率P_gi，t大于设定阈值P_thd时转化为城市类型，否则元胞保持状态不变：

$S_{i,t+1}=\left\{\begin{array}{ccc}Urban& if& P_{gi,t}\&GreaterEqual;P_{thd}\\ Nochange& if& P_{gi,t}<P_{thd}\end{array}\right.---\left(8\right).$

7.根据权利要求6所述的一种基于CA和SAR的海岸带土地利用变化模拟方法，其特征在于，所述的步骤5)具体为；

18)采用GeoDa实现SAR，选用某年的土地利用格局为初始状态利用SARCA模型运行M次得到土地利用变化的模拟及预测结果；

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

8.根据权利要求7所述的一种基于CA和SAR的海岸带土地利用变化模拟方法，其特征在于，所述的步骤6)具体为；

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

首先通过与遥感分类的土地利用格局进行比较，对模拟结果进行精度计算，指标包括总体精度Overall Accuracy、Kappa系数和景观指数Landscape Metrics；

将总体精度分解为城市Hit和非城市Correct Rejection两类、将错误分解为忽略性和替代性两类，利用Pontius方法将错误同时分解为数量错误和位置错误；其中忽略性错误指的是实际为城市但模拟为非城市，即CA模型未能捕捉到的城市元胞；替代性错误指的是实际为非城市但模拟为城市，即CA模型错误地增加的城市元胞；

将模拟结果与遥感分类结果进行叠加，叠加结果包含5类：实际和模拟均为城市、实际为非城市模拟为城市、实际为城市模拟为非城市、实际和模拟均为非城市、以及水体；根据目视判别显示，SARCA模型的模拟结果与实际分类结果的差异。