CN102841963B - 元胞自动机模型尺度敏感性的探测方法 - Google Patents

Abstract

本发明提供的是一种利用正交试验设计来探测土地利用变化元胞自动机模拟尺度敏感性的方法，具体是：首先利用遥感影像来解译提取土地利用信息并进行整合，然后确定元胞自动机模型中需要考查的尺度因素及水平，选取合适的测试尺度敏感性的正交表；在此基础上，模拟研究区域的土地利用变化情况，进而基于精度评判指标Kappa系数，分析出元胞自动机模型中的最优水平组合，并对各尺度因素之间尚存在的相关性进行检验，从而全面客观地分析土地利用变化元胞自动机模拟过程中存在的敏感性。本发明填补了现有方法在土地利用变化元胞自动机模拟过程尺度敏感性分析与表达上的空缺，对确保模拟结果的正确性及高效性具有非常重要的指导意义。

Description

元胞自动机模型尺度敏感性的探测方法

技术领域

本发明涉及探测方法技术领域，特别是一种利用正交试验设计探测土地利用变化元胞自动机模拟尺度敏感性的方法。

背景技术

土地利用变化是近年来全球变化研究的重要领域和热点问题。由于土地利用变化过程错综复杂，因而以简化和抽象为特征的各种模型对于理解和预测土地利用变化过程具有不可替代的作用。元胞自动机（Cellular Automata，简称CA）模型具有“自下而上”的研究思路、强大的复杂计算功能以及固有的并行计算能力，这使得其在模拟土地利用变化方面具有自然性、合理性和可行性，从而成为土地利用变化模拟表达的重要方法和主流工具。

上世纪70年代，Tobler首次正式将CA作为地理模型对当时美国五大湖边底特律地区城市的迅速扩张进行模拟研究。80年代以来，伴随着CA模型理论的深入，美国圣巴巴拉加州大学地理系的Couclelis对CA在地理学中的应用潜力从理论上进行了充分的阐述，奠定了CA模型在地理学中应用的理论框架。90年代以后，越来越多的CA模型应用于城市增长的模拟研究中。国内开展土地利用CA模拟研究工作相对于国外要晚一些，但在90年代后也逐渐活跃起来，比较有影响的是周成虎、黎夏等人。周成虎等人在Batty的CA模型基础上，构建了面向对象的随机异构模型，通过引入GIS空间数据库、遥感土地分类等数据来模拟和预测城市土地发展模型，并于1999年出版了《地理元胞自动机》著作。黎夏等人结合主成分分析和神经网络等技术成功地对东莞市的城市用地扩展进行了模拟。

综合来看，上述这些研究主要集中在获取CA模型的转换规则上，而这些丰富的转换规则基本能够满足从多个角度实现对城市土地利用空间格局全面综合模拟的需要。但是，随着城市土地利用变化过程存在着对尺度依赖性逐渐被认识，利用CA来研究土地利用变化的瓶颈问题，已经从转换规则的构建向尺度敏感性问题研究而过渡。

根据文献，当前对于CA模型尺度敏感性问题主要集中于单一尺度的敏感性研究。Bethng-wolff和Wu用多尺度试验方法对菲尼克斯城镇发展模型进行了校准与确定，文中表明，输入数据的空间分辨率越高，CA模型的精度越高。CA模型其他尺度的敏感性研究相对较少。Liu和Andersson探究了时间间隔对基于CA的城市发展模型的影响。Verda Kocabas和Suzana Dragicevic对CA模型的邻域结构的敏感性进行了分析。

国内黎夏等人对影响该模型模拟精度的多种因素进行了较为全面的理论概述。尹长林等人以长沙市城市增长模型为例，对CA模拟的尺度问题进行了研究，认为CA城市增长模型只有在一定的尺度范围内才具有较高的模拟精度，并且模型对尺度具有一定的敏感性。

此外，上述研究还存在以下几个问题：（1）元胞自动机模型中尺度间关系和尺度组合对各个尺度敏感性的影响以及对模拟精度影响的研究普遍缺乏；（2）对于元胞自动机模拟的尺度选择，缺乏组合优化选取方法的研究；（3）对元胞自动机模型中不同尺度因素之间尚存在的相关性缺乏检验分析。

发明内容

本发明所要解决的技术问题是：提供一种元胞自动机模型尺度敏感性的探测方法，该方法利用正交试验设计的思路，可全面、客观探讨和解决元胞自动机模型应用于土地利用变化模拟中存在的尺度敏感性问题。

本发明解决其技术问题采用以下的技术方案：

本发明提供的元胞自动机模型尺度敏感性的探测方法，是一种利用正交试验设计来探测土地利用变化元胞自动机模拟尺度敏感性的方法，具体是：首先利用遥感影像来解译提取土地利用信息并进行整合，再确定元胞自动机模型中需要考查的尺度因素及水平，据此选取测试尺度敏感性的正交表，根据正交试验方案模拟研究区域的土地利用变化情况，然后基于精度评判指标Kappa系数，分析出元胞自动机模型中的最优水平组合，并对各尺度因素之间尚存在的相关性进行检验，从而实现全面客观地分析土地利用变化元胞自动机模拟过程中存在的敏感性问题。

本发明提供的上述元胞自动机模拟尺度敏感性的方法，可以采用包括以下步骤的方法：

（1）利用遥感影像解译提取研究区域的土地利用类型并对其进行整合，整合后的数据为土地利用分类信息，并将其作为土地利用变化元胞自动机模拟的输入数据；

（2）确定元胞自动机模型中需要考查的尺度因素，针对这些因素采用适当的水平数及合适的正交表并进行表头设计，编制正交试验方案；

（3）根据正交试验方案，完成所有元胞自动机模拟研究区域的土地利用变化情况的试验，并利用Kappa系数作为模拟结果的精度评判标准；

（4）对正交试验设计结果进行分析，探测出元胞自动机模型中的最优水平组合，并对各尺度因素之间尚存在的相关性进行检验。

本发明可以采用以下方法获取所需研究区域的图斑的分类信息，其步骤包括：

（1）用Erdas软件对所研究区域的遥感影像进行解译，依据所需研究的土地类型进行监督分类，获取各类土地的分类信息；

（2）将分类后的遥感影像导入ArcGIS中，获取研究区域的遥感影像的ASCII码文件；再利用ArcToolbox/Conversion Tools/From Raster/Raster to ASCII，完成转换；然后打开转换后的记事本文件，将里面的头文件去掉，得到所述遥感影像的分类信息矩阵。

本发明可以采用以下方法编制正交试验方案，其步骤包括：

（1）确定尺度因素和水平：

元胞自动机模型具有尺度敏感性，影响尺度敏感性的因素包括元胞大小、邻域结构和邻域大小,根据所掌握的信息资料和相关知识来确定所述尺度敏感性因素的水平数；

（2）确定正交表：

基于所述尺度敏感性因素的水平数，以及考查两两尺度因素之间是否交互作用，来确定采用何种正交表；

（3）表头设计并编制正交试验方案：

根据所确定的正交表来确定试验次数，并参照正交表中的二列间交互作用列表来确定正交表中每一列需要放置的尺度因素，确定元胞大小放置第一列，邻域结构放置第二列，元胞大小和邻域结构的交互作用放置第三列，邻域大小放置第四列，元胞大小和邻域大小的交互作用放置第五列，邻域结构和邻域大小的交互作用放置第六列，第七列为空列，用来考查估计试验误差的大小；随后，把正交表中安排各因素的列中的每个水平数字换成该因素的实际水平值，便形成了正交试验方案，所述各因素的列不包含欲考查的交互作用列。

本发明可以采用二水平L₈(2⁷)正交表。

不能发明可以采用元胞自动机模型完成所有土地利用变化模拟试验，其步骤包括：

（1）按照所编制的正交试验方案中不同因素组合来完成元胞自动机模拟试验，具体是：根据正交表中每一行各个因素的实际水平值，采用其对应的元胞自动机模型代码，模拟得出各个因素水平下对应的土地利用变化情况；

（2）根据Kappa系数计算公式代码获得每一个模拟结果的精度评判指标Kappa系数值。

本发明可以根据Kappa系数展开试验结果分析，然后采用正交试验设计中的极差分析法完成结果讨论，其步骤包括：

（1）根据试验结果Kappa系数，计算极差分析结果的各个参数；

（2）总结上述分析计算结果，得出结论，从而寻求最优水平组合的元胞自动机模型，并检验各尺度因素之间是否存在相关性。

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

（1）本发明并不仅仅局限在单一尺度的敏感性分析，而是全面综合地讨论了元胞自动机模型对邻域结构、元胞大小和邻域大小三个尺度因素的敏感性。

例如，在本发明试验中，通过正交试验设计中极差分析法计算得到元胞大小、邻域结构和邻域大小三个因素列的极差值分别为：11.195，2.199，3.873。这些值可以可以说明元胞自动机模型对元胞大小、邻域结构和邻域大小三个尺度因素都具有敏感性，缺少任何一个因素的分析都不能全面研究模型的尺度敏感性。

（2）模型中三个尺度因素之间的关系和尺度组合对各个尺度敏感性的影响以及对模拟精度的影响在本发明的实施下可以更为全面、客观地表达出来。

例如，上述三个尺度因素的极差值从大到小排序为：元胞大小>邻域大小>邻域结构。而由正交试验设计的特点来看，极差值越大，代表这个因素对试验指标的影响越大。因此可得出结论：这三个尺度因素对元胞自动机模型的敏感性从大到小的顺序为：元胞大小>邻域大小>邻域结构。在正交试验设计结果分析中，因素同水平指标值之和Ⅰ、Ⅱ大小顺序可以确定因素各个水平对因素的影响程度，最大的同水平指标值之和表明其对应的水平是该因素的最优水平。计算得到：元胞大小尺度因素两个水平对应的指标值之和分别为14.949和3.754；邻域结构尺度因素两个水平对应的指标值之和分别是10.451和8.252；邻域大小尺度因素两个水平对应的指标值之和分别是11.288和7.415。因此，可以确定元胞大小尺度因素的最优水平是第1个水平，即30m大小的正方形元胞；邻域结构尺度因素的最优水平是第1个水平，即图2所示的邻域结构；邻域大小尺度因素的最优水平是第1个水平，即5*5大小。

（3）正交试验设计是用于多因素试验的一种方法，它是从全面试验中挑选出部分有代表的点进行试验，将其运用在土地利用变化元胞自动机模拟中用以探测及选取最优水平组合尺度因素具有很高的效率。

例如本发明中所涉及的试验借助正交试验设计只需要完成8次（本发明中正交试验方案的试验号数：8次）模拟就可以找出三个尺度因素的显著性以及各个尺度因素的最优水平，同时还可以确定三尺度因素两两之间是否存在相关性。将正交试验设计用于土地利用变化元胞自动机模拟中，得出的结论非常具有信服力。

（4）实现了元胞自动机模型中不同尺度因素之间尚存在的相关性检验研究，填补了现有研究在全面分析及表达土地利用变化元胞自动机模拟过程尺度敏感性上的空缺，有助于更深入挖掘土地利用变化元胞自动机模拟的尺度敏感性问题。

例如，本发明通过试验分析，计算得到元胞大小×邻域结构交互列因素的极差值为0.045；元胞大小×邻域大小交互列因素的极差值为0.031；邻域结构×邻域大小交互列因素的极差值为0.223。这三个数字都非常小，可以充分证明元胞大小、邻域结构和邻域大小三个尺度因素是独立存在的，两两之间不存在交互性。本发明用定量的方法证明了元胞自动机模型中元胞大小、邻域结构和邻域大小之间不存在相关性。而这些相关性的存在与否在其他文献中都是没有被研究过的，因此这项相关性检验研究是本发明中最大的创新点。

附图说明

图1是本发明的技术路线图。

图2是元胞自动机模型中第一种典型的邻域结构示意图：以5*5邻域大小形式展示。

图3是元胞自动机模型中第二种典型的邻域结构示意图：以5*5邻域大小形式展示。

图4是武汉市1987年的土地利用分类结果图。

图5是武汉市1996年的土地利用分类结果图。

图6是武汉市2005年的土地利用分类结果图。

具体实施方式

下面结合实施例及附图对本发明作进一步说明，但并不局限于下面所述内容。

本发明针对现有研究在土地利用变化元胞自动机模拟过程中尺度敏感性分析与表达上的缺陷和不足，提出一种按照技术路线图1所示的研究方法，即利用正交试验设计探测土地利用变化元胞自动机模拟尺度敏感性的方法，该方法是：首先利用遥感影像来解译提取土地利用信息并进行整合，然后确定元胞自动机模型中需要考查的尺度因素及水平，选取合适的测试尺度敏感性的正交表；在此基础上，模拟研究区域的土地利用变化情况，进而基于精度评判指标Kappa系数，分析出元胞自动机模型中的最优水平组合，并对各尺度因素之间尚存在的相关性进行检验，从而全面客观地分析土地利用变化元胞自动机模拟过程中存在的敏感性。

本发明的技术核心内容是将正交试验设计整套思路运用于土地利用变化元胞自动机模拟中，探测出元胞自动机模型的尺度敏感性。因此，根据元胞自动机模型自身的性质及其在土地利用变化模拟中的应用特点合理地实施正交试验设计是非常重要的。

本发明提供的方法采用包括以下步骤的方法：

一. 数据准备并整合基础数据

利用遥感影像来解译提取土地利用基础数据并进行整合，整合后的数据为所需研究区域的图斑的分类信息。其步骤包括：

1. 用Erdas软件对所研究区域的遥感影像进行解译，依据所需研究的土地类型进行监督分类，获取各类土地的分类信息。

在利用遥感影像进行土地利用变化元胞自动机模拟时，其实质是根据遥感影像的分类信息获取若干年后研究区域的土地利用类型变化情况，因此在模拟土地利用变化之前，需要为每一期的遥感影像进行分类，附上不同的值。

2. 将分类后的遥感影像导入ArcGIS中，获取研究区域影像的ASCII码文件。利用ArcToolbox/Conversion Tools/From Raster/Raster to ASCII，完成转换。打开转换后的记事本文件，将里面的头文件去掉，得到研究区域遥感影像的分类信息矩阵。

得到的分类信息矩阵大小与利用Erdas解译提取得到的分类信息遥感影像大小相同，矩阵中值等于不同类型用地对应的代码值。

二. 编制正交试验方案

元胞自动机是一种时间和空间都离散的动力系统，散布在格网里的每一个元胞都取有限且已知的离散状态（在本发明中相当于遥感影像分类信息矩阵里面的具体且有限的数值），通过遵循相同的转换规则作同步更新，格网里的大量元胞就通过这种简单的相互作用而构成动态系统的演化。元胞自动机模型的组成元素是元胞、元胞空间、邻域和转换规则。元胞相当于本研究中的一个个像元；元胞空间相当于本研究中的遥感影像覆盖的范围；邻域则包括邻域结构和邻域大小，目前国际上有两种通用的邻域结构，分别为冯诺依曼邻域结构（见图2所示）和摩尔邻域结构（见图3所示）；转换规则则是根据元胞当前状态及其邻域状态来确定其下一时刻状态的动力函数。通过元胞自动机模型自身的这些特点，我们可以确定其具有尺度敏感性，不同尺度的元胞、元胞空间以及邻域对模型的准确应用都会带来一定程度的影响。因此，必须全面、客观地研究模型的尺度敏感性。首先需要确定模型要考查的尺度因素，然后针对这些因素选取适当的水平，在此基础上选用合适的正交表并进行表头设计，最后编制正交试验方案，参见表1。其步骤包括：

1. 确定尺度因素和水平：

元胞自动机模型具有尺度敏感性，在土地利用变化模拟时，受到元胞大小、邻域结构和邻域大小的影响。因此，本发明中需要考查的尺度因素就是元胞大小、邻域结构和邻域大小。根据所掌握的信息资料和相关知识，若要考虑交互作用，最好选择两水平表，实际工作中一般也不在三水平以上的正交表中安排交互作用。一般情况下，混合正交表也不能考查交互作用。因此，本发明将元胞大小、邻域结构和邻域大小三个尺度因素的水平数均定为两个。邻域结构因素的两个水平分别是冯诺依曼邻域结构和摩尔邻域结构；元胞大小因素的两个水平分别是30m和60m；邻域大小因素的两个水平分别是5*5和9*9。

2. 选用正交表：

本发明中的三个尺度因素的水平数相等，均为两个。因此，需要选用二水平正交表。同时，为了考查三个尺度因素两两之间是否交互作用，则需要至少包含6列的正交表。查阅最新修订的常用正交表，发现L₈(2⁷)正交表可以满足我们的需求。

3. 表头设计并编制正交试验方案：

根据上面选定的L₈(2⁷)正交表，可以确定本发明试验次数为八次。参照L₈(2⁷)二列间交互作用列表来确定正交表中每一列需要放置的尺度因素，确定元胞大小放置第一列，用A来表示；邻域结构放置第二列，用B来表示；元胞大小和邻域结构的交互作用放置第三列，用A×B来表示；邻域大小放置第四列，用C来表示；元胞大小和邻域大小的交互作用放置第五列，用A×C来表示；邻域结构和邻域大小的交互作用放置第六列，用B×C来表示；第七列为空列，用来考查估计试验误差的大小。随后，把正交表中安排各因素的列（不包含欲考查的交互作用列）中的每个水平数字换成该因素的实际水平值，便形成了正交试验方案。

三. 土地利用变化模拟

首先利用上述编制的方案展开试验，分别利用元胞自动机模型模拟得出各个不同因素及水平下研究区域的土地利用变化情况，并计算出模拟结果对应的Kappa系数值。其步骤包括:

1. 编写matlab代码：

基于元胞自动机模型来模拟土地利用变化情况，编写matlab代码。

所述元胞自动机模型实现土地利用变化模拟主要包括两大部分：一是模拟土地利用变化；二是计算模拟结果的精度。在模拟土地利用变化中，发挥核心作用的是元胞自动机模型中的转换规则。转换规则即为根据元胞当前状态及其邻域状态来确定其下一时刻状态的动力函数，而元胞大小、邻域结构和邻域大小这三种尺度因素对在元胞自动机模型中实现转换规则具有较大的影响，也就进一步证明这些尺度因素对模型的敏感性是必须要研究的问题。为尽量减少在获取转换规则的过程中，因为需要计算城市发展适宜性而必须完成的一系列数据转换和叠置运算所导致的误差传递和不确定性，本发明特采用比较客观、人为干扰因素较少的基于统计学的转换规则来实现土地利用变化模拟，主要是采用包括以下步骤的计算方法得到土地利用变化模拟情况：

（1）对元胞的组合状态进行编码：

从遥感影像分类信息矩阵的第一行和第一列开始，分别扫描t1和t2 两时刻整个元胞空间，检测出t1时刻任一中心元胞Cij在其邻域范围内的组合状态码u，同时检测出对应位置处中心元胞Cij在t2时刻的状态码值v及其组合状态码。上述状态码的取值范围为遥感影像分类信息矩阵里面的具体数值，t1和t2分别对应的是两期遥感影像所在的时间，t1是稍早一期的时间。

（2）建立状态转移频数表：

将两时刻在中心元胞Cij邻域范围内对应的组合状态码值一一建立对应关系，统计出状态转移频数s(v,m)。其中v为t1时刻中心元胞Cij邻域范围内与t2时刻中心元胞Cij相同的的状态码值，其值唯一；m为对应位置处在t2 时刻的状态码值，其值不唯一，变化范围是整个遥感影像分类信息矩阵里面的数值。将状态转移频数进行归一化处理，得到状态转移概率p(v,m)，p(v,m)代表的是t2 时刻的中心元胞Cij在欲模拟的下一时刻t3时状态改变为m时的概率。

（3）确定随机预测模型：

将上述获取的多个状态转移概率转换为概率向量P=[p1, p2, …, pk]，k代表的是m值的个数，并且p1, p2,…, pk按照从小到大的顺序排列。对概率向量的元素做累加计算，可以得到概率累加向量PS=[ p1S, p2S, …, pkS]，式中pfS=                                                ，f=1，2，…，k。这样PS中的k个元素就形成了区间（0,1）的一种分割，k个子区间(0,p1S], (p1S,p2S], …, (p(k-1)S, pkS]分别对应结果状态1，2，…k。最终通过系统产生均匀随机数来实现随机预测，即对t2时刻下一个特定的中心元胞Cij的的条件状态组合码值，就产生一个在（0,1）上均匀分布的随机数r，如果r∈(p(k-1)S, pkS]，则判定t3时刻该元胞的状态为k。资料表明，该方法可以充分表达元胞状态变化的多样性和土地利用空间结构的复杂性。

试验结果模拟完成后，就需要对模拟结果进行精度评价，主要是采用包括以下步骤的计算方法得到精度评定系数：将Kappa系数作为精度评判标准，求出模拟结果的Kappa系数值。Kappa系数可用来判断两幅遥感影像的一致性，其公式为：Kappa=(P₀-Pc)/(1-Pc)，式中：P₀是两期遥感影像上土地利用类型的观测一致率；Pc是两期遥感影像上土地利用类型的期望一致率。

在对遥感影像进行处理后，我们得到的是分类信息矩阵。而matlab的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用matlab来解算此类矩阵问题问题要比用C，FORTRAN等语言完成相同的事情简捷得多。

根据以上计算方法，编写matlab代码，主要涉及到五个函数，分别是MooreNeighborhood函数、Neumann函数、topo函数、result函数和Kappa函数。其中MooreNeighborhood和Neumann函数主要是实现邻域提取功能，分别提取摩尔邻域结构和冯诺依曼邻域结构；topo函数主要是实现遥感影像分类信息矩阵的拓展功能，以正确表达元胞空间中四个边缘上的所有元胞的状态变化情况；result函数则是采用上述元胞自动机模型转换规则函数实现土地利用变化模拟功能；Kappa函数则是实现模拟结果的精度评价功能。

matlab代码具体如下，其中，Neumann功能函数邻域结构提取代码是以图2所示的冯诺依曼邻域结构为例，MooreNeighborhood功能函数邻域结构提取代码是以图3所示的摩尔邻域结构为例。

1）Neumann功能函数：

function Pm05=Neumann(num,mt87,mt96)

%根据1987年和1996年的分类信息矩阵，利用冯诺依曼邻域结构模拟2005年的土地利用变化情况

Pm05=zeros(size(mt96));

[mt87Topo,mt96Topo]=topo(num,mt87,mt96);

[r,c]=size(mt96);

for i=1:r

    for j=1:c

        if mt96(i,j)==0

            Pm05(i,j)=0;

        else

%得到以i,j为中心的冯诺依曼邻域

        mtTemp87=mt87Topo(i:i+num-1,j:j+num-1);          mtTemp96=mt96Topo(i:i+num-1,j:j+num-1);

        e=(num-1)/2;

        k=(num+1)/2;

        for x=1:e

            for y=1:k-x

                mtTemp87(x,y)=0;

                mtTemp87(x,num+1-y)=0;

                mtTemp87(num+1-x,y)=0;

                mtTemp87(num+1-x,num+1-y)=0;

            end

        end

             Pm05=result(mt96(i,j),num,mtTemp87,mtTemp96,i,j,Pm05);

        end

    end

end

dlmwrite(['Neumann',num2str(num),'.txt'],Pm05,' ');

end

2）MooreNeighborhood功能函数：

function Pm05=MooreNeighborhood(num,mt87,mt96)

%根据1987年和1996年的分类信息矩阵，模拟2005年的土地利用变化情况

Pm05=zeros(size(mt96));

[r,c]=size(mt96);

[mt87Topo,mt96Topo]=topo(num,mt87,mt96);

for i=1:r

    for j=1:c

        if mt96(i,j)==0

            Pm05(i,j)=0;

        else

%得到以i,j为中心的摩尔邻域

     mtTemp87=mt87Topo(i:i+num-1,j:j+num-1);  

     mtTemp96=mt96Topo(i:i+num-1,j:j+num-1);

     Pm05=result(mt96(i,j),num,mtTemp87,mtTemp96,i,j,Pm05);

        end

    end

end

dlmwrite(['Moore',num2str(num),'.txt'],Pm05,' ');

end

3）topo功能函数：

%将矩阵进行拓展，将外拓展的值赋予最外层的值，实现topo()函数功能

function [m87Topo,m96Topo]=topo(num,mt87,mt96)

[r,c]=size(mt87);

m1=int8(zeros((num-1)/2,c));

m2=int8(zeros(r+num-1,(num-1)/2));

m87Topo=cat(1,m1,mt87,m1);

m87Topo=cat(2,m2,m87Topo,m2);

m96Topo=[m1;mt96;m1];

m96Topo=[m2,m96Topo,m2];

end

4）result功能函数：

%第一个参数是96年的矩阵，第三四个参数分别是拓展后得到的矩阵，

% Pm05为最终结果

function Pm05 = result(mt96z,num,mtTemp87,mtTemp96,i,j,Pm05)

     b=zeros(6,6);

     for t1=1:num

         for t2=1:num

              u=mtTemp87(t1,t2);

              v=mtTemp96(t1,t2);

              if u~=0&&v~=0

              b(u,v)=b(u,v)+1;

              end

         end

     end

     c=b(mt96z,:);

     total=0;

     p=zeros(1,6);

     if c==0

         Pm05(i,j)=mt96z;

     else

         r=rand();

         for x=1:6

             total=total+c(x);

         end

         p(1)=c(1)/total;

         if r<p(1)

             Pm05(i,j)=1;

         else

         for x=2:6

             p(x)=c(x)/total+p(x-1);

             if r<p(x)

                 Pm05(i,j)=x;

                 break;

             end

         end

         end

     end

end

5）Kappa功能函数：

%mtPred是通过计算预测得到的矩阵，mtReal是真实矩阵，classes是土地利用类型数

function K = Kappa(mtPred,mtReal,classes)

%矩阵像元总数

s=0;

Preclass=zeros(1,classes);         

Realclass=zeros(1,classes);

[r,c]=size(mtPred);

N=size(find(mtPred~=0),1);

for i=1:r

    for j=1:c

        if mtPred(i,j)~=0

        %对应像元值相等的像元总数

        if mtPred(i,j)==mtReal(i,j)

            s=s+1;

        end

        %预测矩阵和真实矩阵每一类像元总数

       for m=1:classes

           if mtPred(i,j)==m

               Preclass(m)=Preclass(m)+1;

           end

           if mtReal(i,j)==m

               Realclass(m)=Realclass(m)+1;

           end

       end

        end

    end

end

P=Preclass*Realclass;

K=(N*s-P)/(N^2-P);

end

2. 调用函数完成土地利用变化元胞自动机模拟：

将步骤1得到的1987年和1996年的遥感影像分类信息矩阵导入matlab软件之中，调用函数求取2005年该研究区域的土地利用变化情况并进行模拟精度评价。调用函数的代码只有两个，具体如下：

MooreNeighborhood(num,mt87,mt96)；

%num,mt87,mt96在调用函数时均为具体的值，其中：

%num为邻域大小的值，在本发明中主要是5和9两个值；

%mt87是研究区域1987年的遥感影像分类信息矩阵；

%mt96是研究区域1996年的遥感影像分类信息矩阵；

K = Kappa(mtPred,mtReal,classes)；

%mtPred,mtReal,classes在调用函数时均为具体的值，其中：

%mtPred为上述模拟得到的结果矩阵；

%mtReal是研究区域2005年的遥感影像分类信息矩阵；

%classes是研究区域遥感影像的分类数目，在本发明中该值为6。

四、正交试验设计结果分析

在正交试验设计结果分析过程中，应将步骤二中的编制正交试验方案和步骤三中的土地利用变化模拟的模型方法共同作用后最终得到的模拟结果精度数据Kappa系数值，完整记录于试验结果分析表中，进而运用极差分析法对试验结果展开讨论：分清元胞自动机模型中各尺度因素及其交互作用的主次顺序，判断各尺度因素对试验指标影响的显著程度，找出试验因素的优水平和试验范围内的最优组合，分析尺度因素与试验指标之间的关系，了解我们欲探讨的元胞大小、邻域结构和邻域大小两两之间是否存在交互作用。具体包括以下步骤：

1. 将步骤三中最终得到的Kappa系数数据值完整记录于试验结果分析计算表中；

2. 对上述试验结果分析表中的数据，利用极差分析法展开结果讨论，并将结果记录于试验结果分析表中。

极差分析中需要计算的物理量有各列因素的同水平指标值和，即Ⅰ和Ⅱ。在此基础上，求出每一个因素的极差，即为该因素的同水平指标值和的最大值与最小值的差值。最后，根据试验指标Kappa系数值计算得出T值，即8次试验中所得到的总的指标值。

3. 分析上述试验分析表中的数据，对各尺度因素及其交互作用的主次顺序、各尺度因素对试验指标影响的显著程度、试验因素的优水平和试验范围内的最优组合、尺度因素与试验指标之间的关系以及交互作用情况进行分析，从而达到本发明的目的。

本发明提供的上述方法可以为全面、客观地探讨元胞自动机模型应用于土地利用变化模拟中存在的尺度敏感性问题提供理论与实际研究参考依据。下面以具体实例说明：

（1）以1987、1996和2005年武汉市Landsat TM遥感影像数据为基础数据，在Erdas软件中对三幅影像进行前期的几何校正和波段合成后，就着重利用监督分类的方法解译提取土地利用基础数据。结合武汉市实际情况，将研究区域的土地利用类型主要分为水域、建筑用地、林地、农业用地、草地和未利用地六大类，各用地类型对应的代码值依次为：

1代表水域

2代表建筑用地

3代表林地

4代表农业用地

5代表草地

6代表未利用地

解译后得到的土地利用分类影像图如图6所示。

（2）本发明中所使用的数据分辨率是30m，为考虑元胞大小对模型模拟结果的影响，需要对数据进行重采样，得到元胞大小为60m的数据。主要是将这三期已分类的遥感影像数据导入Erdas中，通过interpreter—>utilities—>degrade模块达到目的，得到三期重采样分类信息矩阵。

（3）将监督分类后的遥感影像（图4-图6所示）导入ArcGIS中，分别获取这三幅分类影像的ASCII码文件。分别用记事本打开转换后的ASCII文件，将里面的头文件去掉，得到研究区域三期遥感影像分类信息矩阵，分类信息矩阵里面对应的数值就是上述步骤中各用地类型所对应的数字。

（4）参照表1所示的土地利用变化元胞自动机模拟敏感性探测试验8组试验方案，选择对应的元胞自动机模型尺度因素，准备进行模拟试验。这里需要注意的是，每一行试验的顺序都可以是随机的。

（5）根据试验方案表1，完成8组土地利用变化模拟试验。将研究区域三期遥感影像的分类信息矩阵以及经过重采样操作后得到的60m分辨率分类信息矩阵分别导入matlab软件，基于上述元胞自动机模型模拟土地利用变化的算法，选择不同的邻域结构、元胞大小和邻域大小对应的模拟代码，分别模拟出8个试验号对应的研究区域在2005年的土地利用变化情况。这里，邻域结构分别选择图2-图3所示的两种典型的邻域结构；元胞大小分别选择30m和60m；邻域大小所代表的数字分别选择5和9。

（6）将模拟产生的8组文本数据和武汉市2005年实际的土地利用分类信息矩阵一同导入matlab中，依据Kappa函数分别计算模拟结果的精度值，得到的结果参见表1中的Kappa系数值。

（7）利用极差分析法对试验结果展开讨论，分析计算结果数据见表2，为方便后续正交设计分析参数的计算，表2中Kappa系数的值均减去50，此操作对结果分析不会造成影响。

对表2所示的土地利用变化元胞自动机模拟敏感性探测试验结果展开分析如下：

正交试验设计用于土地利用变化元胞自动机模拟中，能够全面、客观地探测出模型中存在的尺度敏感性程度。在给定试验因素及水平下，利用极差分析方法，本发明可以定量评价模型中各尺度因素及其交互作用的主次顺序、各尺度因素对试验指标影响的显著程度、试验因素的优水平、试验范围内的最优组合以及尺度因素与试验指标之间的关系。极差分析法可以圆满迅速的求得试验的优化成果——主次因素、优水平、优搭配及最优组合，极差值越大表示因素对试验指标的影响越大，极差值越小说明该因素对试验指标影响不显著；在计算极差值过程中产生的因素同水平指标值之和中的最大值对应的水平为改因素的最优水平。

探测武汉市土地利用变化元胞自动机模拟尺度敏感性，主要分为四个步骤：第一，对原始数据分类并确定分类信息矩阵；第二，编制正交试验方案；第三，模拟土地利用变化情况并计算模拟结果的精度评判指数Kappa系数；第四，计算各因素的极差值及同水平指标值之和，根据计算结果对试验结果展开分析讨论。

从表2中可以看出，交互作用列A×B、A×C及B×C的极差值R分别是0.045、0.031和0.023。这三个值都非常小，说明因素A和B、A和C以及B和C之间不存在交互作用，这三个因素是独立存在的，两两之间不存在相关性。这一点从空列即正交表的第7列也可以看出来，该列上的极差值R为0.073，这个值非常小，也验证了无交互作用存在的假设，可作为试验误差来处理。

三个因素列A、B和C的极差值R都比较大，分别为11.195、2.199和3.873，按照从大到小排序为A、C、B，说明元胞自动机模型对A、B和C代表的尺度因素都存在敏感性。此外，对A因素而言，Ⅰ的值为14.949，Ⅱ 的值为3.754，Ⅰ比Ⅱ值大，说明A的最优水平为Ⅰ对应的水平；对B因素而言，Ⅰ的值为10.451，Ⅱ 的值为8.252，说明其最优水平为Ⅰ对应的水平；C因素的Ⅰ值为11.288，Ⅱ值为7.415，反映出其最优水平也为Ⅰ对应的水平。因此可确定本发明中土地利用变化元胞自动机模拟敏感性探测正交试验的因素主次顺序为：元胞大小、邻域大小、邻域结构；最优水平组合为A₁B₁C₁ ，即邻域结构为图2中所示的冯诺依曼邻域结构，元胞大小为30m，邻域大小为5*5。

附表

表1

表2

Claims (5)

1. 一种元胞自动机模型尺度敏感性的探测方法，其特征是一种利用正交试验设计来探测土地利用变化元胞自动机模拟尺度敏感性的方法，具体是：首先利用遥感影像来解译提取土地利用信息并进行整合，再确定元胞自动机模型中需要考查的尺度因素及水平，据此选取测试尺度敏感性的正交表，根据正交试验方案模拟研究区域的土地利用变化情况，然后基于精度评判指标Kappa系数，分析出元胞自动机模型中的最优水平组合，并对各尺度因素之间尚存在的相关性进行检验，从而实现全面客观地分析土地利用变化元胞自动机模拟过程中存在的敏感性问题；

该探测方法采用包括以下步骤的方法：

（1）利用遥感影像解译提取研究区域的土地利用类型并对其进行整合，整合后的数据为土地利用分类信息，并将其作为土地利用变化元胞自动机模拟的输入数据；

（2）确定元胞自动机模型中需要考查的尺度因素，针对这些因素采用适当的水平数及合适的正交表并进行表头设计，编制正交试验方案；

（3）根据正交试验方案，完成所有元胞自动机模拟研究区域的土地利用变化情况的试验，并利用Kappa系数作为模拟结果的精度评判标准；

（4）对正交试验设计结果进行分析，探测出元胞自动机模型中的最优水平组合，并对各尺度因素之间尚存在的相关性进行检验；

采用以下方法编制正交试验方案，其步骤包括：

（1）确定尺度因素和水平：

元胞自动机模型具有尺度敏感性，影响尺度敏感性的因素包括元胞大小、邻域结构和邻域大小,根据所掌握的信息资料和相关知识来确定所述尺度敏感性因素的水平数；

（2）确定正交表：

基于所述尺度敏感性因素的水平数，以及考查两两尺度因素之间是否交互作用，来确定采用何种正交表；

（3）表头设计并编制正交试验方案：

根据所确定的正交表来确定试验次数，并参照正交表中的二列间交互作用列表来确定正交表中每一列需要放置的尺度因素，确定元胞大小放置第一列，邻域结构放置第二列，元胞大小和邻域结构的交互作用放置第三列，邻域大小放置第四列，元胞大小和邻域大小的交互作用放置第五列，邻域结构和邻域大小的交互作用放置第六列，第七列为空列，用来考查估计试验误差的大小；随后，把正交表中安排各因素的列中的每个水平数字换成该因素的实际水平值，便形成了正交试验方案，所述各因素的列不包含欲考查的交互作用列。

2.根据权利要求1所述的元胞自动机模型尺度敏感性的探测方法，其特征是采用以下方法获取所需研究区域的图斑的分类信息，其步骤包括：

（1）用Erdas软件对所研究区域的遥感影像进行解译，依据所需研究的土地类型进行监督分类，获取各类土地的分类信息；

（2）将分类后的遥感影像导入ArcGIS中，获取研究区域的遥感影像的ASCII码文件；再利用ArcToolbox/Conversion Tools/From Raster/Raster to ASCII，完成转换；然后打开转换后的记事本文件，将里面的头文件去掉，得到所述遥感影像的分类信息矩阵。

3.根据权利要求1所述的元胞自动机模型尺度敏感性的探测方法，其特征是采用二水平L₈(2⁷)正交表。

4.根据权利要求1所述的元胞自动机模型尺度敏感性的探测方法，其特征是采用元胞自动机模型完成所有土地利用变化模拟试验，其步骤包括：

（1）按照所编制的正交试验方案中不同因素组合来完成元胞自动机模拟试验，具体是：根据正交表中每一行各个因素的实际水平值，采用其对应的元胞自动机模型代码，模拟得出各个因素水平下对应的土地利用变化情况；

（2）根据Kappa系数计算公式代码获得每一个模拟结果的精度评判指标Kappa系数值。

5.根据权利要求4所述的元胞自动机模型尺度敏感性的探测方法，其特征是根据Kappa系数展开试验结果分析，然后采用正交试验设计中的极差分析法完成结果讨论，其步骤包括：

（1）根据试验结果Kappa系数，计算极差分析结果的各个参数；

（2）总结上述分析计算结果，得出结论，从而寻求最优水平组合的元胞自动机模型，并检验各尺度因素之间是否存在相关性。