元胞自动机(Cellular Automata)与城市规划及其MATLAB实现——莆田市城市发展预测
github:https://github.com/Myoontyee/CA-city-planning-for-Putian-by-MATLAB
前言
探索元胞自动机用于城市规划，是由于前不久在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的结果，是不正确的结果。
代码
使用的图片
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<Urbanization)
cells(i,j)=2;
end
if(aroundCenter(i,j,cells))
if(rand<propagateFactor)
cells(i,j)=2;
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 边界增长
if(cells(i,j)==2 && rand<diffuseFactor)
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<proliferateFactor)
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
输出
城市发展模拟过程如下图：
其他
结果与真实情况相差较大，内部算法有待改善。
若你来自兰州，也可以用文件夹里的兰州地图跑一遍程序，是有意思的小例子。
参考
[1]基于元胞自动机的城市规划
[2]【matlab】基于元胞自动机的城市规划（程序修正）
[3]https://wenku.baidu.com/view/ddf76d8c50e2524de4187e45.html