【元胞自动机】基于元胞自动机模拟森林大火matlab源码

一.元胞自动机定义和适用范围

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

元胞自动机——应用于森林火灾和传染病场景

最近接触了元胞自动机模型，做了一些资料搜查，并进行学习，推荐这篇文章澳洲变燠洲，考拉成烤拉！澳大利亚山火为什么难以控制？\ 以下对所学进行记录。

森林火灾元胞自动机原理

在元胞自动机模型中，空间被离散成网格，每一个网格被称为元胞。森林火灾元胞有三种状态：树，火(正在燃烧的树)和空(空地)状态。元胞下一时刻状态的更新规则如下：\ 树变火：一棵树，其上下左右若有一个状态为火，下一刻就会变成火。或者一棵树遇上闪电，下一刻就会变成火。由于遇上闪电着火的概率Plight很小。

火变空：火在下一时刻会变成空。

空变树：空地下一时刻会以一个很小的概率Pgrowth长出新树。

改进模型会考虑树的对角位置有没有着火。或者会考虑风向(比如吹西风(火从东吹向西)，火的西边着火的机率会变大(顺风)，火的东边着火的几率变小(逆风))，这里盗个图：\ 图a是基础元胞自动机，图b是考虑对角的元胞自动机，图c是吹西风的元胞自动机。\ 本人这里对基础元胞自动机，考虑对角情况的元胞自动机，考虑吹西风的元胞自动机三种模型进行仿真，仿真结果如下：\ \ \

```

P=[];

for m=1:10

clear D T firetimelightning firetimeitself aspect tdata Index;

%% Orginal

%the first 3 dimension is RGB，R is the fire，G is the tree.

%Black is the meaning of no tree.

global n D T Y firetimelightning firetimeitself firetimedemend pull aspect count1 tdata Pulltimes

n=500; % the length of the forest matrix

D=zeros(n);

T=zeros(n);

Y=zeros(n,n,3); % draw the picture by matrix Y(RGB)

firetimelightning=0;

firetimeitself=0;

firetimedemend=0;

times=1;

pull=1/times; % the rate of pull the fire out

aspect=ceil(rand(1)*4); % 1 is right,2 is up,3 is left and 4 is down.

count_1=0;

tdata=[]; % the day by each fire happened.

Pull_times=0;

f1=1/1000; % f1 is the probability in ceil when it being struck by lightning.

f2=1/500; % f2 is the probability in ceil when it being fired itself.

Z=Terrain();

[scale_b,S]=Forest(Z);

Tem=S;

Yi=imshow(Y);

set(gcf,'DoubleBuffer','on');

% set up the double cache to prevent the flash in palying animation constantly

t=0;

tp=title(['T = ',num2str(t)]);

%ap=title(['aspect = ',num2str(aspect)]);

%while 1

% t=t+1;

for t=1:2400

%% Fire in the early time

if rem(t,50)==0&&t<1000

Fire_Demand(S);

end

%% Lightning

if rand

OriginFireLightning(S);

set(Yi,'CData',Y);

end

%% Fire itself

if t>10

OriginFireCritical(S,t,f2);

set(Yi,'CData',Y);

end

%% FireRule1

[a,b]=FireRule1(S);

%% FireRule2

S=FireRule2(S,a,b);

set(Yi,'CData',Y);

set(tp,'string',['T = ',num2str(fix(t/(24))),'D = ',num2str(t),'h '])

%set(ap,'string',['aspect = ',num2str(aspect)])

pause(2e-6)

end

scale_n=sum(sum(Y(:,:,2)));

firecount=firetimeitself+firetime_lightning;

Lostarea=scaleb-scale_n;

Lostrateall=Lostarea/scaleb;

Lostrateaverage=Lostrateall/(fire_count);

Tem=Tem-S;

Tem=Tem.*Tem;

%Lost_Value=sum(sum(Tem));

%N=[pull Lostrateall Lostrateaverage fire_count]

P=[P;firetimedemend Pulltimes Lostrateall firetimeitself firetime_lightning]

%figure (2)

%plot(tdata(:,1),tdata(:,2));

end ```