前言

探索元胞自动机用于城市规划，是由于前不久在CSDN上看到相关案例后大开眼界，兴趣使然，想对家乡做一个城市发展预测，遂在巨人的肩膀上做一些探索与更正。文章末尾有这些案例的链接，感谢并致敬这些先行者。

背景

概念

元胞自动机(cellular automata，CA) 是一种时间、空间、状态都离散，空间相互作用和时间因果关系为局部的网格动力学模型，具有模拟复杂系统时空演化过程的能力。

不同于一般的动力学模型，元胞自动机不是由严格定义的物理方程或函数确定，而是用一系列模型构造的规则构成。凡是满足这些规则的模型都可以算作是元胞自动机模型。因此，元胞自动机是一类模型的总称，或者说是一个方法框架。其特点是时间、空间、状态都离散，每个变量只取有限多个状态，且其状态改变的规则在时间和空间上都是局部的。

城市模型

1950s~1960s

静态、均衡、宏观

1980s+

动态、微观

环境

Windows7 x64

MATLAB R2016a

特点

i)简单、直观

ii)离散

iii)灵活、开发

iv)易与GIS、遥感数据处理系统结合

用法

输入

图像

f(x)

CA for city planning

输出

城市规模仿真模拟结果

参数

针对程序中的可修改参数给出释义。

1)扩散因子

代码：

diffuseFactor=0.2; % 扩散因子

影响城市化过程中，元胞的扩散速率

输入

扩散因子数值

输出

定义扩散因子

Tip

输入应为正数，接受＞1的数值

2)繁殖因子

代码：

proliferateFactor=0.2; % 繁殖因子

影响城市化过程中，元胞的繁殖速率

输入

繁殖因子数值

输出

定义繁殖因子

Tip

输入应为正数，接受＞1的数值

3)传播因子

代码：

propagateFactor=0.2; % 传播因子

影响城市化过程中，元胞的传播速率

输入

传播因子数值

输出

定义传播因子

Tip

输入应为正数，接受＞1的数值

4)城市化因子

代码：

Urbanization=0.0004; % 城市化因子

影响城市化速率

输入

城市化因子数值

输出

定义城市化因子

Tip

输入应为正数，接受＞1的数值，但数值过大会导致求解失真，下图为城市化因子为2，迭代次数为3的结果，是不正确的结果。

城市化因子为2

代码

使用的图片

i)莆田市地图

莆田市地图

ii)经过二值化处理的莆田市地图

莆田市地图二值化结果

1)CityCA.m

clear;close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Author: Myoontyee.Chen

%% Data：20181227

%% License：BSD 3.0

%% CA - City Planning

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 图像预处理

I=imread('city.png');

level=graythresh(I); % 图像灰度处理

bw=im2bw(I,level); % 图像二值化处理

I=bw;

cells=double(I);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% button定义

plotbutton=uicontrol('style','pushbutton',...

'string','Run',...

'fontsize',12,...

'position',[100,400,50,20],...

'callback','runModel=1;');

erasebutton=uicontrol('style','pushbutton',...

'string','Stop',...

'fontsize',12,...

'position',[300,400,50,20],...

'callback','stopModel=1;');

Iterations=uicontrol('style','text',...

'string','1',...

'fontsize',12,...

'position',[20,400,50,20]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Initialization参数初始化

cells(33,44)=2;

cells(88,31)=2;

cells(33,80)=2;

%% 可修改参数

diffuseFactor=0.2; % 扩散因子

proliferateFactor=0.2; % 繁殖因子

propagateFactor=0.2; % 传播因子

Urbanization=0.0004; % 城市化因子

%% 建议默认参数

sch=[];skz=[];t=1;

[x,y]=size(cells);

runModel=0;

stopModel=0;

stop=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 法则定义

while (stop==0)

%% 按下Run

if(runModel==1)

for i=2:x-1

for j=2:y-1

if(cells(i,j)~=1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 自然增长

if(cells(i,j)==0)

if(rand

cells(i,j)=2;

end

if(aroundCenter(i,j,cells))

if(rand

cells(i,j)=2;

end

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 边界增长

if(cells(i,j)==2 && rand

if(existCase(i,j,cells))

m=1+3*rand;

switch m

case 1

ii=i-1;jj=j;

case 2

ii=i;jj=j-1;

case 3

ii=i;jj=j+1;

otherwise 4

ii=i+1;jj=j;

end

if(canCity(ii,jj,cells))

cells(ii,jj)=2;

end

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 新扩展中心型

if(cells(i,j)==2 && existCase(i,j,cells))

if(rand

if(canCity(i,j,cells))

cells(i,j)=3;

end

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end

end

end

ch=0;kzch=0;

for i=1:x

for j=1:y

if(cells(i,j)==2) ch=ch+1;end

if(cells(i,j)==3) kzch=kzch+1;end

end

end

sch(t)=ch;skz(t)=kzch;

t=t+1;

[A,B]=size(cells);

Area(1:A,1:B,1)=zeros(A,B);

Area(1:A,1:B,2)=zeros(A,B);

Area(1:A,1:B,3)=zeros(A,B);

for i=1:A

for j=1:B

if cells(i,j)==1

Area(i,j,:)=[1,1,1]; % 黑色

elseif cells(i,j)==0

Area(i,j,:)= [255, 255, 255]; % 白色

elseif cells(i,j)==3

Area(i,j,:)= [255,0,0]; % 红色

elseif cells(i,j)==2

Area(i,j,:)= [255,177,0]; % 橙色

end

end

end

pause(0.0005);

Area=uint8(Area);

imagesc(Area);

axis equal;

axis tight;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 迭代次数记录

IterationNum=1+str2num(get(Iterations,'string'));

set(Iterations,'string',num2str(IterationNum));

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 按下Stop

if stopModel==1

runModel=0;

stopModel=0;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 暂停并刷新图像

drawnow

end

2)existCase.m

function result= existCase(i,j,cells)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Author: Myoontyee.Chen

%% Data：20181227

%% License：BSD 3.0

%% existCase - Correspondence rule

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

a=0;

if(cells(i-1,j)==2) a=a+1; end

if(cells(i,j-1)==2) a=a+1;end

if(cells(i,j+1)==2) a=a+1;end

if(cells(i+1,j)==2) a=a+1;end

if(a>=2)

result=1;

else

result=0;

end

end

3)canCity.m

function result=canCity(i,j,cells)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Author: Myoontyee.Chen

%% Data：20181227

%% License：BSD 3.0

%% canCity - Correspondence rule

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

s=0

if(cells(i-1,j-1)==1) s=s+1;end

if(cells(i-1,j)==1) s=s+1;end

if(cells(i-1,j+1)==1) s=s+1;end

if(cells(i,j-1)==1) s=s+1;end

if(cells(i,j+1)==1) s=s+1;end

if(cells(i+1,j-1)==1) s=s+1;end

if(cells(i+1,j)==1) s=s+1;end

if(cells(i+1,j+1)==1) s=s+1;end

if(s>=4)

result=0;

else

result=1;

end

end

4)aroundCenter.m

function a=aroundCenter(i,j,cells)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Author: Myoontyee.Chen

%% Data：20181227

%% License：BSD 3.0

%% aroundCenter - Correspondence rule

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

a=0;

if(cells(i-1,j)==3) a=1; end

if(cells(i,j-1)==3) a=1;end

if(cells(i,j+1)==3) a=1;end

if(cells(i+1,j)==3) a=1;end

end

输出

城市发展模拟过程如下图：

迭代次数2

迭代次数671

其他

结果与真实情况相差较大，内部算法有待改善。

若你来自兰州，也可以用文件夹里的兰州地图跑一遍程序，是有意思的小例子。

参考