考虑到多种生物之间具有复杂捕食和共生关系，且不同区域的水环境之间也会有相互作用，常见的方程不容易简单描述，故采用元胞自动机进行模拟，以总藻类密度作为指标，分析虾、鱼、蟹、藻生长情况。
本文所涉及的建模思路、MATLAB源代码均参照于华中农大特等奖——《淡水养殖池塘水华发生及池水自净化研究》（提取码：7qog ）
元胞自动机（𝐶𝑒𝑙𝑙𝑢𝑙𝑎𝑟 𝐴𝑢𝑡𝑜𝑚𝑎𝑡𝑜𝑛，𝐶𝐴）是由𝑉𝑜𝑛 𝑁𝑒𝑢𝑚𝑎𝑛𝑛提出的一种在时间、空间、状态上都离散的动力系统，利用相邻元胞之间的互相作用关系，能够以极简单的规则模拟复杂的行为，很适合用于生态养殖池的水质分析。 一个完整的元胞自动机系统通常包括元胞维数(𝐷)、邻居（𝑁）、状态（𝑆）及演化规则（𝑓），记为𝐴 = (𝐷, 𝑁, 𝑆, 𝑓)。其中，元胞是构成元胞自动机的最小单位，在本文中设置为 S = ( I , I I , I I I , I V , V ) S=(I,II,III, IV, V) S=(I,II,III,IV,V)五种藻类密度状态，以总藻类密度作为衡量指标，具体划分如下， v v v表示该元胞的总藻类密度值。
考虑到三维的元胞自动机太复杂，为了简单起见，本文采用维度为D=2。考虑到水体之间极强的流动性和分子扩散性，采用𝑀𝑜𝑜𝑟𝑒型邻居会更合理。
Moore型邻居 具体参数设定
此次共有5个脚本文件，其中fish.m为脚本文件，其他4个都是函数文件，模拟的情况是对于无水华的情况进行模拟，下面分别做详细说明。
fish.m 代码如下（示例）：
clear; clc; rho=50;length=50; fish=[500 300 500];%300kg鱼、500kg虾、500kg蟹 lambda=0.8; T=10;%10周 alpha=[0.00001,0.00015,0.00035,0.00025,0.0002];%不同状态下的鱼儿们消耗系数 range=[15,50,150,500]; %==========画方块图 area=MyCreate(length,rho,lambda); DrawArea(area,range); n=length;%n=50 plot([(0:n)',(0:n)']+0.5,[0,n]+0.5,'k');%k代表黑色 plot([0,n]+0.5,[(0:n)',(0:n)']+0.5,'k'); axis image; set(gca,'xtick',[]); set(gca,'ytick',[]);%消除x,y轴 figure; %========= data=zeros(T,3);%10行3列的零 for i=1:T [area,fish]=MyChange(area,fish,range,alpha); DrawArea(area,range); drawnow; data(i,:)=fish; end plot([(0:n)',(0:n)']+0.5,[0,n]+0.5,'k'); plot([0,n]+0.5,[(0:n)',(0:n)']+0.5,'k'); axis image; set(gca,'xtick',[]); set(gca,'ytick',[]); figure; %========== plot(data,'LineWidth',2); grid on; legend('虾','鱼','蟹'); title('鱼虾蟹数量趋势图'); xlabel('周'); ylabel('公斤');DrawArea.m 代码如下（示例）：
function DrawArea(area,range) temp=area; n=length(temp); Area(:,:,1)=temp; Area(:,:,2)=temp; Area(:,:,3)=temp; %给每个方格块上色，水藻密度越大，颜色越绿 for i=1:n for j=1:n if temp(i,j)<range(1) Area(i,j,:)=[0 230 214]; elseif temp(i,j)<=range(2) Area(i,j,:)=[139 255 223]; elseif temp(i,j)<=range(3) Area(i,j,:)=[147 255 170]; elseif temp(i,j)<=range(4) Area(i,j,:)=[91 240 91]; else Area(i,j,:)=[4 176 58]; end end end Area=uint8(Area); p=imagesc(Area); hold on; endMyRound.m 代码如下（示例）：
function total=MyRound(area,i,j) n=length(area); Area=ones(n+2)*area(i,j); Area(2:n+1,2:n+1)=area; i=i+1; j=j+1; %根据Moore型邻居 赋予权重 权重之和为1 total=0.2*Area(i,j)+0.1*Area(i+1,j)+0.1*Area(i-1,j)+... 0.1*Area(i,j+1)+0.1*Area(i,j-1)+0.1*Area(i+1,j+1)+... 0.1*Area(i+1,j-1)+0.1*Area(i-1,j+1)+0.1*Area(i-1,j-1); endMyCreate.m 代码如下（示例）：
function area=MyCreate(length,rho,lambda) area=ones(length)*rho;%length=50 rho=50 area=area.*(1+(rand(length)-0.5)*lambda);%初始水藻的密度大都是50-150 属于无水华的情况 endMyChange.m 代码如下（示例）：
function [area,fish]=MyChange(area,fish,range,alpha) n=length(area); for i =1:n for j=1:n %周围9个方块共同impact if MyRound(area,i,j)<range(1)%15 area(i,j)=MyRound(area,j,j)*(1-rand*mean(fish)*alpha(1)); %(1-rand*mean(fish)*alpha(1))鱼虾蟹等影响系数 elseif MyRound(area,i,j)<range(2)%50 area(i,j)=MyRound(area,j,j)*(1-rand*mean(fish)*alpha(2)); elseif MyRound(area,i,j)<range(3)%150 area(i,j)=MyRound(area,j,j)*(1-rand*mean(fish)*alpha(3)); elseif MyRound(area,i,j)<range(4)%500 area(i,j)=MyRound(area,j,j)*(1-rand*mean(fish)*alpha(4)); else area(i,j)=MyRound(area,j,j)*(1-rand*mean(fish)*alpha(5)); end end end avg=mean(mean(area)); %成长系数设置的很迷惑 %蟹类主要以鱼虾等动物尸体为食物，所以其数目会受到鱼虾的影响 if avg<range(1) fish=fish*(1+(rand-0.5)*0.01); elseif avg<range(2) fish=fish.*[1.006 1.005 1.007]; elseif avg<range(3) fish=fish.*[1.013 1.011 1.013]; elseif avg<range(4) fish=fish.*[1.03 1.03 1.03]; else fish=fish.*[1.05 1.04 1.05]; end fish=fish+fish.*unifrnd(-0.12,0.1,1,3);%有随机性，更有说服力 end根据以上结果可以知道，经过10周的模拟，水质有了很大的改善，鱼虾蟹的数量有不同程度的下降，说明此时藻类的含量已经不足以支持鱼虾蟹的生长，需要人工投喂。
