1、中山大学遥感与地理信息工程系2009.07.16劳春华trycourlchqqftp:/202.116.70.210gisgis一、CA概念CACA英文全称是英文全称是Celluar AutomataCelluar Automata，中文译名为元胞自动机，又有，中文译名为元胞自动机，又有人称之为细胞自动机。人称之为细胞自动机。CACA是一种时间、空间、状态都是一种时间、空间、状态都离散离散，( (空间上的空间上的) )相互作用和相互作用和( (时时间上的间上的) )因果关系皆因果关系皆局部局部的的格网格网动力学模型。具有模拟复杂系统动力学模型。具有模拟复杂系统时空演化过程的能力。时空演化过程的
2、能力。19481948年，数学家年，数学家Von NeumannVon Neumann首次提出元胞自动机首次提出元胞自动机(CA)(CA)的概念。的概念。二、CA组成t时刻状态t+1时刻状态转换规则CACA由由 “元胞元胞” 、“邻域邻域”和和“转换规则转换规则” 三部分组成，元胞具有三部分组成，元胞具有“状态状态”属性属性例如例如12碰上奇数+1碰上偶数+356碰上奇数+1碰上偶数+3碰上奇数+1碰上偶数+3元胞状态由元胞状态由1 1经过三次转换迭代变成经过三次转换迭代变成6 6。如果任由元胞演变下去，将会。如果任由元胞演变下去，将会产生一个复杂的无穷数列。产生一个复杂的无穷数列。三、CA分
3、类元胞自动机的构建没有固定的数学公式，构成方式繁杂，变种元胞自动机的构建没有固定的数学公式，构成方式繁杂，变种很多，行为复杂，故其分类难度也较大很多，行为复杂，故其分类难度也较大 。基于不同的出发点，元胞自动机可有多种分类基于不同的出发点，元胞自动机可有多种分类 。其中，最具影。其中，最具影响力的当属响力的当属S. WolframS. Wolfram在在8080年代初做的年代初做的基于动力学行为的元胞自基于动力学行为的元胞自动机分类动机分类，而，而基于维数的元胞自动机分类基于维数的元胞自动机分类也是最简单和最常用也是最简单和最常用的划分。的划分。 三、CA分类基于动力学行为的元胞自动机基于动力
4、学行为的元胞自动机(1)(1)平稳型平稳型: :自任何初始状态开始自任何初始状态开始, ,经过一定时间运行后，元胞空间趋于一个空间平稳的经过一定时间运行后，元胞空间趋于一个空间平稳的构形，这里空间平稳即指每一个元胞处于固定状态。不随时间变化而变化。构形，这里空间平稳即指每一个元胞处于固定状态。不随时间变化而变化。(2)(2)周期型周期型: :经过一定时间运行后，元胞空间趋于一系列简单的固定结构经过一定时间运行后，元胞空间趋于一系列简单的固定结构(Stable (Stable Paterns)Paterns)或周期结构或周期结构(Perlodical Patterns)(Perlodical P
5、atterns)。由于这些结构可看作是一种滤波器。由于这些结构可看作是一种滤波器(Filter)(Filter)，故可应用到图像处理的研究中。，故可应用到图像处理的研究中。(3)(3)混沌型混沌型: :自任何初始状态开始，经过一定时间运行后，元胞自动机表现出混沌的非周自任何初始状态开始，经过一定时间运行后，元胞自动机表现出混沌的非周期行为，所生成的结构的统汁特征不再变止，通常表现为分形分维特征。期行为，所生成的结构的统汁特征不再变止，通常表现为分形分维特征。(4)(4)复杂型复杂型: :出现复杂的局部结构，或者说是局部的混沌，其中有些会不断地传播。从另出现复杂的局部结构，或者说是局部的混沌，其
6、中有些会不断地传播。从另一角度，元胞自动机可视为动力系统，因而可将初试点、轨道、不动点、周期轨和终极一角度，元胞自动机可视为动力系统，因而可将初试点、轨道、不动点、周期轨和终极轨等一系列概念用到元胞自动机的研究中轨等一系列概念用到元胞自动机的研究中 三、CA分类基于维数的元胞自动机基于维数的元胞自动机一维元胞自动机一维元胞自动机二维元胞自动机二维元胞自动机三维元胞自动机三维元胞自动机高维元胞自动机高维元胞自动机四、CA应用CACA应用应用社会学社会学 生物学生物学 生态学生态学 数学数学 物理学物理学 化学化学 地理学地理学 研究经济危机的形成与爆发过程研究经济危机的形成与爆发过程 等等肿瘤细
7、胞的增长机理和过程模拟肿瘤细胞的增长机理和过程模拟 等等生物群落的扩散模拟生物群落的扩散模拟 等等研究数论和并行计算研究数论和并行计算 等等用于磁场、电场等场的模拟，以及热扩用于磁场、电场等场的模拟，以及热扩散、热传导和机械波的模拟散、热传导和机械波的模拟 等等海上石油泄露后的油污扩散、工厂周围海上石油泄露后的油污扩散、工厂周围废水、废气的扩散等过程的模拟废水、废气的扩散等过程的模拟 四、CA应用地理学上的应用CACA应用应用土地利用变化土地利用变化城市扩展城市扩展人口迁移人口迁移 火灾蔓延火灾蔓延 沙漠化沙漠化 洪水掩没洪水掩没 交通控制交通控制 五、生命游戏模型最经典的CA模型Martin
8、 C(1970,1971)Martin C(1970,1971)将生命游戏规则引入到数字游戏中。该游将生命游戏规则引入到数字游戏中。该游戏通过分布在二维空间网格上的细胞来发挥作用。每个细胞只戏通过分布在二维空间网格上的细胞来发挥作用。每个细胞只以一种状态存在（以一种状态存在（0 0或或1 1）, ,并且在下个时刻的状态由当前状态以并且在下个时刻的状态由当前状态以及与它最近的及与它最近的8 8个邻居的状态共同决定。个邻居的状态共同决定。五、生命游戏模型最经典的CA模型定义了如下定义了如下3 3种转换规则：种转换规则：生存规则，周围有生存规则，周围有2 2个或者个或者3 3个活着的邻居细胞，该活着
9、的细胞个活着的邻居细胞，该活着的细胞将在下一时刻继续生存；将在下一时刻继续生存；死亡规划，周围活着的细胞有死亡规划，周围活着的细胞有3 3个以上，或者少于个以上，或者少于2 2个，该活着个，该活着的细胞将在下一时刻死亡；的细胞将在下一时刻死亡；繁殖规则，周围存活邻居数达到繁殖规则，周围存活邻居数达到3 3个，该死亡细胞在下一时刻被个，该死亡细胞在下一时刻被激活过来激活过来五、生命游戏模型最经典的CA模型从数学模型的角度看，该模型将平面划分成方格棋盘，每个方格从数学模型的角度看，该模型将平面划分成方格棋盘，每个方格代表一个元胞。元胞状态：代表一个元胞。元胞状态：0 0死亡，死亡，1 1活着；领域
10、半径：活着；领域半径：MooreMoore型；演化规则型；演化规则五、生命游戏模型最经典的CA模型演示五、生命游戏模型最经典的CA模型五、基于空间数据挖掘的CA模型遥感影像：T1遥感影像：T2空间数据挖掘算法CA转换规则T时刻状态（T+1）时刻状态逻辑回归CA神经网络CA决策树CA蚁群CA支持向量机CA五、基于空间数据挖掘的CA模型逻辑回归五、基于逻辑回归的CA模型nnxbxbxbaz2211zgep11 逻辑回归不同于线性回归，它研究的是一个事件发生的概率，与其他因素之间的关系。根据随机试验的结果，通过最大似然法对回归参数进行估计。五、基于逻辑回归的CA模型Logistic CA主要由三大部
11、分组成，分别是全局性开发概率和局主要由三大部分组成，分别是全局性开发概率和局部作用的邻域影响以及随机项。这三部分相乘，得出最终转换部作用的邻域影响以及随机项。这三部分相乘，得出最终转换概率。当转换概率大于给定阈值，发生由非城市用地到城市用概率。当转换概率大于给定阈值，发生由非城市用地到城市用地的转变，否则不发生转变。地的转变，否则不发生转变。nnxbxbxbaz2211zgep11)ln(1RARAppgthresholdthresholdtppdUnDevelopeppDevelopedS,1五、基于逻辑回归的CA模型准备数据操作流程处理数据编写代码模拟输出五、基于逻辑回归的CA模型准备数据
12、数据准备数据准备20192019年东莞市土地利用分类数据年东莞市土地利用分类数据(2019.img)(2019.img)东莞市市中心点数据东莞市市中心点数据(Prop.shp)(Prop.shp)东莞市镇中心点数据东莞市镇中心点数据(Town.shp)(Town.shp)东莞市铁路线数据东莞市铁路线数据(Rail.shp)(Rail.shp)东莞市高速公路数据东莞市高速公路数据(Express.shp)(Express.shp) 东莞市一般公路数据东莞市一般公路数据(Road.shp)(Road.shp)以东莞市以东莞市2019年到年到2019年为例年为例20192019年东莞市土地利用分类数
13、据年东莞市土地利用分类数据(2019.img)(2019.img)五、基于逻辑回归的CA模型数据处理2019.img2019.img2019.img2019.imgTown.shpTown.shpRail.shpRail.shpExpress.shpExpress.shp Road.shpRoad.shpUrban2019.imgUrban2019.imgUrban2019.imgUrban2019.imgDisTown.imgDisTown.imgDisRail.imgDisRail.imgDisExpress.imgDisExpress.imgDisRoad.imgDisRoad.imgU
14、rbanChange.imgUrbanChange.imgProp.shpProp.shpDisProp.imgDisProp.imgUrban2019.txtUrban2019.txtUrban2019.txtUrban2019.txtdianData.shpdianData.shp五、基于逻辑回归的CA模型数据处理UrbanChange.imgUrbanChange.imgdianData.shpdianData.shpTown.shpTown.shpRail.shpRail.shpExpress.shpExpress.shp Road.shpRoad.shpDisTown.imgDisT
15、own.imgDisRail.imgDisRail.imgDisExpress.imgDisExpress.imgDisRoad.imgDisRoad.imgProp.shpProp.shpDisProp.imgDisProp.imgdianValue.dbfdianValue.dbfDisTown.imgDisTown.imgDisRail.imgDisRail.imgDisExpress.imgDisExpress.imgDisRoad.imgDisRoad.imgDisProp.imgDisProp.imgdianValue.dbfdianValue.dbf五、基于逻辑回归的CA模型数据
16、处理Zfile.imgZfile.imgPgFile.imgPgFile.img五、数据处理-获取UrbanChange.img加载2019年和2019年遥感分类图2019年遥感分类图2019年遥感分类图通过栅格运算，计算出2019年和2019年城市和非城市遥感分类图2019年和2019年城市和非城市遥感分类图如右图所示从下图可以看出，影像分辨率太高，行列数太多，可进行重采样，适当调低分辨率左图是重采样对话框，我们把分辨率调成85.5米可以看出，分辨率已经调成了85.5米打开打开2019年和年和2019年属性表，年属性表，发现取值只有发现取值只有0和和1,我们把这两我们把这两年数据进行合成年数
17、据进行合成合成后的数据，如下对合成后的数据进一步处理，得到对合成后的数据进一步处理，得到2019年和年和2019年城市变化遥感图，年城市变化遥感图，1为新增的，为新增的，0为不变为不变的，的，2为为01年是城市的，年是城市的，05年还是城市年还是城市下图是进一步处理好的数据下图是进一步处理好的数据导出导出01年到年到05年城市变化遥感数据，取名为年城市变化遥感数据，取名为UrbanChange.img打开打开erdas9.2,对对UrbanChange.img进行采点，首先把进行采点，首先把Urbanchange.img的的Layer Type改成改成thematic打开打开UrbanChan
18、ge.img，我们可以看到它本来的，我们可以看到它本来的Layer Type是是Continuous把把UrbanChange.img的的Layer Type改成改成ThematicClassifier-Accuracy Assessment，打开右下图窗口，打开右下图窗口打开打开UrbanChange.img文件，文件，Edit-Create/Add Random Points，打开生成随机点窗口打开生成随机点窗口点击点击Select Classes，打开属性编辑窗口，选打开属性编辑窗口，选择择1,设置采样点和搜索设置采样点和搜索数，这里采数，这里采5000个点，个点，生成的随机点如右图生成
19、的随机点如右图所示所示把采到的点输出为把采到的点输出为dat数据，这里命名为数据，这里命名为diandata.dat利用同样的方法，对利用同样的方法，对0值进行采样，这里采值进行采样，这里采20000个点个点输出为输出为diandata2.dat在我的电脑中看在我的电脑中看到点数据文件如到点数据文件如下下在在excel中打开中打开把把diandata2.dat中的数据合中的数据合到到diandata.dat中来中来在第一行中插入一行，输入在第一行中插入一行，输入x,y作为标题名作为标题名保存成保存成csv格式，用记事本打格式，用记事本打开，如右图所示开，如右图所示在我的电脑中，直接把在我的电脑
20、中，直接把diandata.csv改成改成diandata.txt，使用，使用arcMap加载该点数据加载该点数据arcMap-tools-Add XY Data，打开窗口如右图所示，打开窗口如右图所示打开打开diandata.txt,如右图所示，如右图所示，这时点数据没有投影，点这时点数据没有投影，点edit按钮，为点数据加投影按钮，为点数据加投影点点Import按钮，选择按钮，选择UrbanChange.img，把其投影导进来，把其投影导进来导进投影如上图所示导进投影如上图所示确认后，确认后，arcMap根据点坐标生成矢量根据点坐标生成矢量点数据，如上图点数据，如上图点数据，放大图点数据，
21、放大图把点数据导出保存为把点数据导出保存为diandata.shp生成矢量点数据后，把市中心、镇中心、生成矢量点数据后，把市中心、镇中心、铁路、高速公路、一般公路的矢量数据铁路、高速公路、一般公路的矢量数据加进来，准备生成空间距离栅格数据加进来，准备生成空间距离栅格数据下图是加进来的数据，用于生成空间距下图是加进来的数据，用于生成空间距离变量栅格数据离变量栅格数据设定栅格运算的范围设定栅格运算的范围为为UrbanChange.img的范围，的范围，cell大小为大小为UrbanChange.img的的大小大小开始计算离市中心距开始计算离市中心距离，生成栅格数据离，生成栅格数据生成的离市中心距离
22、生成的离市中心距离栅格数据如上栅格数据如上生成离镇中心空间距生成离镇中心空间距离栅格数据离栅格数据生成的离镇中心距离生成的离镇中心距离栅格数据如上栅格数据如上生成离铁路空间距离生成离铁路空间距离栅格数据栅格数据生成的离铁路距离栅生成的离铁路距离栅格数据如上格数据如上生成的离高速公路空生成的离高速公路空间距离栅格数据如上间距离栅格数据如上生成的离一般公路空生成的离一般公路空间距离栅格数据如上间距离栅格数据如上下图是生成的栅下图是生成的栅格数据格数据为了消除量纲影为了消除量纲影响，可对空间距响，可对空间距离栅格数据进行离栅格数据进行归一化处理归一化处理归一化离市中心归一化离市中心距离栅格数据距离栅
23、格数据导出已经完成归一化的导出已经完成归一化的数据，存为数据，存为DisProp_gyh.img用同样的方法，归一化用同样的方法，归一化其它空间距离变量栅格其它空间距离变量栅格数据，如左图所示数据，如左图所示Spatial Analyst ToolExtractionSample，对已经归一，对已经归一化的栅格数据和化的栅格数据和UrbanChange.img进行采进行采样，结果存为样，结果存为DianValue.dbf从我的电脑上看采样好的数据从我的电脑上看采样好的数据在在spss中打开采样好的数据，其中，只有中打开采样好的数据，其中，只有列列z_z_z2c1到到6是有用的是有用的按顺序把列
24、名改好，顺序为采样的时候，按顺序把列名改好，顺序为采样的时候，添加数据的顺序添加数据的顺序返回数据视图，发现有些点出现误差，返回数据视图，发现有些点出现误差，UrbanChange的值为的值为2,应该去掉，应该去掉，data-Select Cases选择选择UrbanChange的值不等的值不等2的行的行选择删掉未选中的数据选择删掉未选中的数据返回数据视图中，这时，数据返回数据视图中，这时，数据已经是可用的了已经是可用的了Analyze-Regression-Binary Logistic.，进行二项，进行二项逻辑回归分析逻辑回归分析回归出来的系数的误差如下表所示，回归出来的系数的误差如下表所
25、示，在在ArcMap中进行栅中进行栅格运算，格运算，算出算出Z值值算出的算出的Z值如上图所示值如上图所示导出成导出成Zfile.img文件文件再进一步算出再进一步算出Pg值值Pg值数据如上图所示值数据如上图所示导出为导出为PgFile.img文件文件将将PgFile.img转换成转换成PgFile.txt右图是右图是2019年和年和2019年城年城市和非城市分类图市和非城市分类图进行模拟的时候，可以把进行模拟的时候，可以把水体加进来，取值为水体加进来，取值为2。新的栅格图的取值为：新的栅格图的取值为：0,非城市；非城市；1城市；城市；2水体水体把合成的栅格图导出为把合成的栅格图导出为Urban
26、2019.img和和Urban2019.img再把再把Urban2019.img和和Urban2019.img转换成转换成Urban2019.txt和和Urban2019.txt，作为模拟，作为模拟的输入数据的输入数据从从2019.img中提取中提取开发适宜性数据，这开发适宜性数据，这里提取水体出来，取里提取水体出来，取值为值为0,其它为其它为从从2019.img中提取开发适中提取开发适宜性数据，这里提取水体宜性数据，这里提取水体出来，取值为出来，取值为0,其它为其它为提取出来的土地适宜性数据，如上图所示，这提取出来的土地适宜性数据，如上图所示，这里也可以把保护区的数据加进来里也可以把保护区的
27、数据加进来把土地适宜性文件把土地适宜性文件导出来，取名为导出来，取名为LandSuitable.img把把LandSuitable.img转换成转换成LandSuitable.txt，作为模拟时，作为模拟时的输入数据的输入数据五、基于逻辑回归的CA模型编写代码输入输入Urban2019.txtUrban2019.txtUrban2019.txtUrban2019.txtPgFile.txt PgFile.txt LandSuitable.txt LandSuitable.txt UrbanSimulate2019.txt UrbanSimulate2019.txt CACA迭代迭代 输出输出运
28、算运算五、基于逻辑回归的CA模型核心代码变量Public data(,) As Int32 2019urban.txt 数据，以列行存储数据，以列行存储Public dataFinal(,) As Int32 2019urban.txt 数据，以列行存储数据，以列行存储Public tempData(,) As Int32 临时数据临时数据Public PgData(,) As Double PgFile.txt 数据数据，以列行存储，以列行存储Public suitableData(,) As Double LandSuitable.txt 数据数据，以列行存储，以列行存储Public no
29、DataValue As Int32 无值数据无值数据Public xCor() As Int32 变化变化元胞的列坐标元胞的列坐标Public yCor() As Int32 变化变化元胞的行坐标元胞的行坐标Public upData() As Int32 变化变化元胞的数据值元胞的数据值 Public chgNumber As Int32 变化的点变化的点Public rdm As Random 产生随机数类产生随机数类Public realUrbanNumber As Int32 实际城市数目实际城市数目Public simUrbanNumber As Int32 模拟城市数目模拟城市数
30、目五、基于逻辑回归的CA模型核心伪代码for 每一行每一行 for 每一列每一列 if data(列列,行行)=NoData or data(列列,行行)=1 or data(列列,行行)=2 then tempdata(列列,行行)= data(列列,行行) else 计算领域影响计算领域影响con 计算随机因子影响计算随机因子影响rdmdata 读取土地适宜性因子读取土地适宜性因子suitabledata 读取读取PgFile.txt中的开发概率中的开发概率Pg 计算总开发概率计算总开发概率P=con*rdmdata*suitabledata*Pg if PPthreshold then
31、tempdata(列列,行行)=1 else tempdata(列列,行行)=data(列列,行行) end ifend forend for五、基于逻辑回归的CA模型核心代码 Public Sub stickOne() Ks += 1 Ks为迭代次数为迭代次数 Dim i As Int32, j As Int32 For j = 0 To rows - 1 For i = 0 To cols - 1如果该元胞值处于无数据状态或者已经是城市或者是水体，则值不变如果该元胞值处于无数据状态或者已经是城市或者是水体，则值不变 If data(i, j) = noDataValue Or data(i
32、, j) = 1 Or data(i, j) = 2 Then tempData(i, j) = data(i, j)否则，计算该元胞城市开发概率否则，计算该元胞城市开发概率 Else -第一步，第一步，计算领域影响计算领域影响- Dim con As Double = 0 Dim tempI As Int32, tempJ As Int32 tempI = i - 1 If tempI = 0 Then tempJ = j - 1 If tempJ = 0 Then If data(tempI, tempJ) = 1 Then con += 1 End If tempJ = j If dat
33、a(tempI, tempJ) = 1 Then con += 1 tempJ = j + 1 If tempJ = 0 Then If data(tempI, tempJ) = 1 Then con += 1 End If tempJ = j + 1 If tempJ = rows - 1 Then If data(tempI, tempJ) = 1 Then con += 1 End If tempI = i + 1 If tempI = 0 Then If data(tempI, tempJ) = 1 Then con += 1 End If tempJ = j If data(temp
34、I, tempJ) = 1 Then con += 1 tempJ = j + 1 If tempJ = 1 Then rungDa = rungDa - 0.00001 rdmData = Pow(-Log(rungDa), Rfa) + 1 -读取城市发展适宜因子读取城市发展适宜因子- Dim suitable As Double 城市发展适宜因子城市发展适宜因子 suitable = suitableData(i, j) -读取空间变量发展概率读取空间变量发展概率Pg- Dim Pg As Double 空间变量发展概率空间变量发展概率 Pg = PgData(i, j) 计算城市开发概
35、率计算城市开发概率P Dim p As Double 总发展概率总发展概率 p = Pg * con * suitable * rdmData 五、基于逻辑回归的CA模型核心代码 根据城市发展概率判断该元胞是否开发为城市根据城市发展概率判断该元胞是否开发为城市 If p Pthreshold And rungDa = 1.0 / Ks Then tempData(i, j) = 1 Else tempData(i, j) = data(i, j) End If End If NextNext将临时数据将临时数据tempdata赋回给赋回给data，进行下一次迭代，进行下一次迭代 For i = 0 To cols - 1 For j = 0 To rows - 1 data(i, j) = tempData(i, j) Next Next end
