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

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

森林火灾元胞自动机原理

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

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

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

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

代码如下：

m = 300;
n = 300; % 表示森林的矩阵行列 m x nPlight = 5e-6;    % 闪电概率
Pgrowth = 1e-2;   % 生长概率% 邻居方位 d 和点燃概率 p
d1 = {[1,0], [0,1], [-1,0], [0,-1]};%最下一行往最上，最右一列往最左，最前一列往最下，最左一列往最右，circshift(S,d{j})
p1 = [    1,     1,      1,      1];% % 改进元胞自动机
d2 = {[1,0], [0,1], [-1,0], [0,-1], [1,1], [-1,1], [-1,-1], [1,-1]};%考虑对角
p2 = [ones(1,4), ones(1,4)*(sqrt(1/2)-1/2)];% % 考虑风的情况,西风,因此西边方向概率为0.3,东方向概率为1,增加对角（0.3），东二格子（0.8）
d3 = {[1,0], [0,1], [-1,0], [0,-1], [1,1], [-1,1], [-1,-1], [1,-1], [0,-2]};
p3 = [ 0.80,  0.30,   0.80,   1.00,  0.12,   0.12,    0.30,   0.30,    0.8];% 空=0, 火=1, 树=2
E = 0;
F = 1;
T = 2;
S = T*(ones(m,n));%全树
S(m/2,n/2)=1;imh = image(cat(3, S==F, S==T, zeros(m,n)));
axis image;
%plot box fits tightly around the data
for t = 1:3000
% 计算邻居中能传播着火的个数sum = zeros(size(S));note='  吹西风';p=p3;d=d3;for j = 1:length(d)sum = sum + p(j) * (circshift(S,d{j})==F);endisE = (S==E); isF = (S==F); isT = (S==T);     % 找出三种不同的状态ignite = rand(m,n)<sum | (rand(m,n)<Plight);  % 着火条件% 规则 1: 着火FRule1 = F*(isT & ignite);% 规则 2: 烧尽ERule2 = F*isF - F*isF;%这里的结果是0矩阵，除去新着火，新生及原来没烧的树木，剩下的就是烧尽的空（0）% 规则 3: 新生及原来没烧的树木TRule3 = T*(isE & rand(m,n)<Pgrowth)+T*(isT & ~ignite) ;S = Rule1 + Rule2 + Rule3;Fnum=length(find(S==F));Tnum=length(find(S==T));Enum=length(find(S==0));%返回true的索引，一列一列来figure(1);set(imh, 'cdata', cat(3, isF, isT, zeros(m,n)) )%红绿蓝str1=['树:',num2str(Tnum),' 火:',num2str(Fnum),' 空:',num2str(Enum)];str2=['时间：',num2str(t),'次',note];title({str1,str2});%title{[],[]}换行显示drawnow%Use this command if you modify graphics objects and want to see the updates on the screen immediately.pause(0.1);end

传染病元胞自动机原理

传染病元胞自动机中包含四种元胞状态，分别为易感、潜伏、感染和移除（死亡或治愈）。元胞自动机的规则如下：

如果易感元胞的东、南、西、北四个邻居中有潜伏元胞，则易感元胞下一时刻以概率P0转变为潜伏元胞。

潜伏元胞经过T0个时间步后转变为感染元胞。

感染元胞经过C个时间步后转变为移除元胞。

估计T0=7;C=5;
gamma=1/C=1/5,gamma是治愈率,C是治愈周期。
beta=kb,k是每个病人每天接触的人，b是传染概率,beta是每个病人每天有效接触的人，P0=beta/(gamma4*T0)，参数估计内容更详细请见武汉加油——传染病模型拟合以及大揭秘！武汉肺炎传播四大数学模型（附代码）

仿真结果如下：
代码如下：

%%%%%%%%%%%%%%%%参数估计###########33
d=0.03;%死亡率
global C;%治愈时间
C=5;
global N gamma
gamma=1/C;
N=13000000;t0=0;
t1=35;
t2=41;
I0=0;
I1=log(1182);
I2=log(2758);global k;
k=1;
global k1;
k1=1.5;
t=[t0,t1,t2];
I=[I0,I1,I2];
p=polyfit(t,I,1);
b=(p(1)+gamma)/k;
global beta beta2;
beta=b*k;
beta2=b*k1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m = 300;
n = 300; % 表示森林的矩阵行列 m x n%Plight = 5e-6;    % 闪电概率
%Pgrowth = 1e-2;   % 生长概率% 邻居方位 d 和点燃概率 p0
global T0;
T0=7;%潜伏天数,R0=beta2/gamma,周围四个有患病危险
p0=(beta2/gamma)/(4*T0);
d1 = {[1,0], [0,1], [-1,0], [0,-1]};%最下一行往最上，最右一列往最左，最前一列往最下，最左一列往最右，circshift(S,d{j})
p1 = [p0,p0,p0,p0];% % 改进元胞自动机
d2 = {[1,0], [0,1], [-1,0], [0,-1], [1,1], [-1,1], [-1,-1], [1,-1]};%考虑对角
p2 = [ones(1,4), ones(1,4)*(sqrt(1/2)-1/2)];% 易感=0, 潜伏=1, 患病=2，移除=3
S = 0;
E = 1;
I = 2;
R = 3;
people = zeros(m,n);
people(ceil(m/2),ceil(n/2))=1;imh = image(cat(3, people==I, people==R, people==E));
axis image;
%plot box fits tightly around the data
E_time=(people==E);
I_time=(people==I);
for t = 1:3000
% 计算邻居中能传播着火的个数sum = zeros(size(people));%note='  传染';p=p1;d=d1;for j = 1:length(d)sum = sum + p(j) * (circshift(people,d{j})==E);end%s上一步isS = (people==S); isE = (people==E);isI = (people==I);     % 找出三种不同的状态isR = (people==R);if(E*isE+S*isS+I*isI+R*isR~=people)disp('error');end%本步incubate = rand(m,n)<sum;E_time=E_time+isE+incubate;    I_time=I_time+isI+(E_time>=T0);E_time=E_time-T0*(E_time>=T0);remove=(I_time>=C);%移除条件;    I_time=I_time-C*(I_time>=C) ; % 变潜伏    Rule1 = E*(E_time~=0);    % 规则 2:变感染Rule2 =I*(I_time~=0);% 规则 3: 变移除Rule3 =R*(remove)+R*(isR);  %isR是上一步的 people = Rule1 + Rule2 + Rule3;Snum=length(find(people==S));Enum=length(find(people==E));Inum=length(find(people==I));%返回true的索引，一列一列来Rnum=length(find(people==R));figure(1);set(imh, 'cdata', cat(3, isI, isR, isE) )%红绿蓝str1=['健康者(黑):',num2str(Snum),' 潜伏者（蓝）:',num2str(Enum),' 感染者（红）:',num2str(Inum),' 移除者（绿）:',num2str(Rnum)];str2=['时间：',num2str(t),'次'];title({str1,str2});%title{[],[]}换行显示drawnow%Use this command if you modify graphics objects and want to see the updates on the screen immediately.pause(0.1);end