简介「元胞自动机—森林火灾模型MATLAB代码.pdf」面向学习复杂系统模拟与元胞自动机建模的高校学生及科研入门者用一份可直接运行的MATLAB脚本演示森林火灾的起火、蔓延与自然恢复过程。资源包仅含1个PDF文件压缩后约64KB内容为带注释的完整代码清单。代码以二维矩阵表示森林状态0为空地、1为绿树、2为燃烧树木按四条规则迭代演化燃烧树转为空地、绿树被相邻火点引燃、空地以概率p0.3长出新树、绿树以概率f6e-5因闪电自燃。邻居统计借助矩阵平移叠加实现while循环配合图像刷新可实时观察火势扩散。读者可借此掌握元胞自动机元胞、状态、邻域、规则四要素的落地写法学习向量化邻域求和、随机概率判定与动态可视化技巧并可通过调节生长与闪电概率观察火灾规模差异。目前已有1072人学习。1. 元胞自动机森林火灾模型在 MATLAB 里到底在算什么森林火灾模型最反直觉的一点是它不预测某一场火从哪里烧到哪里而是把“树生长—闪电点燃—火烧连片—空地重新长树”压成几条概率规则看系统自己会不会走到一个临界状态。元胞自动机把林地切成方格每个格子只有空地、树、着火三种状态下一时刻的状态只由当前时刻自己和邻居决定。MATLAB 适合做这件事因为矩阵就是网格conv2、逻辑索引和imagesc能把规则和可视化串成几十行代码。这个标题对应的内容适合做复杂系统、生态建模、概率仿真的人也适合拿 MATLAB 做教学示例的读者。真正要盯住的是三个参数空地长树的概率、闪电点燃的概率、邻域半径它们决定了火是小打小闹还是烧穿整片网格。2. 森林火灾模型的元胞规则与状态编码2.1 三种状态与四个演化规则最常用的状态编码是整数0 表示空地1 表示树2 表示着火。同步更新时每个格子按当前时刻的邻居状态计算下一时刻状态不能一边算一边改。四个规则通常这样写着火的格子下一时刻变成空地树如果邻居里有着火下一时刻变成着火空地以概率 p 长出树树以概率 f 被闪电击中并着火。这里把 p 叫生长概率把 f 叫闪电概率f 一般远小于 p因为闪电点燃在真实林火里是小概率事件但正是这个小概率让系统不会永久停在全是树的状态。同步更新是森林火灾模型里最容易被忽略的细节。如果写成顺序更新先遍历到的树被点燃后同一轮里它右边的树又会被它点燃火会在一轮内“走”很远导致蔓延速度被严重高估。同步更新要求先把所有格子的新状态算完再整体替换。MATLAB 里可以用逻辑索引一次性完成不必写双重循环逐个格子判断。2.2 邻域传播用 conv2 还是双重循环邻域定义直接决定火的蔓延形状。4 邻域只考虑上下左右火边界更方正8 邻域把对角也算进去火边界更圆。两者在临界行为上会有差异做参数扫描时要固定一种不能中途换。计算邻居里有多少个着火格子常见做法有两种双重循环和conv2。双重循环直观但 MATLAB 里循环慢N200、迭代几千步时差距很明显。conv2把邻域核当成卷积核一次矩阵运算就能得到每个格子的着火邻居数代码短、速度快。邻域类型核矩阵蔓延形状适用场景4 邻域[0 1 0; 1 0 1; 0 1 0]方正边界沿轴向扩展规则简单、强调轴向传播8 邻域[1 1 1; 1 0 1; 1 1 1]圆润对角也能传火更接近自然火蔓延边界零填充卷积默认边缘邻居少火不易出界固定边界边界周期需手动补边左右上下相连无限大林地近似% 用 conv2 统计每个格子的 8 邻域着火数 kernel [1 1 1; 1 0 1; 1 1 1]; neighbor_fire conv2(double(grid 2), kernel, same); % neighbor_fire(i,j) 就是 (i,j) 周围 8 个格子中着火格的数量conv2的第三个参数same表示输出尺寸和输入一致方便直接和grid做逻辑运算。double(grid 2)把着火格子变成 1其余变成 0。核矩阵中心为 0因为格子自己不能算自己的邻居。若用 4 邻域把核换成[0 1 0; 1 0 1; 0 1 0]即可。边界默认按零填充相当于林地外面没有火边缘格子的邻居数会少一些。如果想让火从另一边绕回来可以在卷积前用padarray补边算完再裁掉但多数森林火灾模型用固定边界就够了。2.3 概率参数 p、f 与边界条件怎么定参数没有唯一标准值但有一组常用的起步范围。网格 N 取 100 到 300初始树密度取 0.5 到 0.7p 取 0.001 到 0.05f 取 1e-6 到 1e-4。p 越大空地恢复成树越快系统树密度越高f 越大闪电点燃越频繁火事件更多但每次火烧面积可能更小。做临界现象观察时常把 f 固定得很小然后扫描 p看火事件大小分布是否出现长尾。参数含义典型范围调大后的效果N网格边长100~300计算量按 N² 增长火事件更充分p空地长树概率0.001~0.05树密度上升火更容易连片f闪电点燃概率1e-6~1e-4点火频率上升稳态树密度下降rho0初始树密度0.5~0.7影响达到稳态前的暂态长度邻域4 或 8固定选一8 邻域火蔓延更快边界条件要和统计量一起考虑。固定边界下火靠近边缘时邻居少可能提前熄灭周期边界下火可以从一侧烧到另一侧适合研究无限大系统的统计性质。如果只是做教学演示固定边界加imagesc已经足够。若要做树密度、火事件大小的定量统计建议用周期边界并且每次实验换随机种子跑多组取平均。闪电概率很小时系统可能几百步都不着火这不是代码错了而是 f 太小、等待时间太长可以先把 f 临时调大观察规则是否正确再调回小值做正式实验。3. 用 MATLAB 搭一个可运行的森林火灾模型最小代码3.1 网格初始化与状态矩阵先确定网格尺寸和初始状态。rand(N) rho0生成逻辑矩阵再转成 double得到 0 和 1其中 1 表示树0 表示空地。初始时随机选一个格子设成 2表示第一把火。这样系统不会一开始就静止能立刻看到火蔓延。初始化时不要把所有树都设成 2否则第一步全图着火看不到传播过程。rng(1)固定随机种子方便复现同一组结果做参数扫描时则要换种子或跑多次平均。N 150; % 网格边长 rho0 0.6; % 初始树密度 p 0.01; % 空地长树概率 f 1e-5; % 闪电点燃概率 rng(7); % 固定随机种子方便复现 grid double(rand(N) rho0); % 0 空地1 树 grid(randi(N*N)) 2; % 随机点燃一棵树这段初始化里rand(N) rho0产生约 rho0 比例的树。randi(N*N)返回 1 到 N² 的整数用来选一个线性索引位置。若想从边界点燃可以把索引改成sub2ind([N N], 1, randi(N))。初始树密度不建议取 1因为全树状态下第一把火会烧掉几乎整个网格统计上不好区分是模型临界还是初始条件太极端。3.2 单步更新函数怎么写把单步更新写成独立函数主循环只负责调用和画图。下面这个函数输入当前网格、p 和 f输出下一时刻网格和本步燃烧格数。燃烧格数用来统计火事件大小后面参数扫描会用到。函数内部先算着火邻居数再依次处理着火变空地、树被点燃、空地长树、闪电点燃。顺序不能乱如果先长树再判断点燃新长出来的树在同一轮里也可能被闪电击中概率含义会变。function [grid, burned] step_fire(grid, p, f) N size(grid, 1); tree (grid 1); burning (grid 2); kernel [1 1 1; 1 0 1; 1 1 1]; neighbor_fire conv2(double(burning), kernel, same); new_grid grid; new_grid(burning) 0; % 着火格变空地 new_grid(tree neighbor_fire 0) 2; % 邻居有火树被点燃 empty (new_grid 0); new_grid(empty (rand(N) p)) 1; % 空地长树 tree_now (new_grid 1); new_grid(tree_now (rand(N) f)) 2; % 闪电点燃 burned sum(burning(:)); % 本步燃烧格数 grid new_grid; end参数说明grid是 N×N 整数矩阵取值 0、1、2p和f是标量概率burned是本步从树变成火的格子数量也就是当前步火的大小。逻辑说明neighbor_fire 0表示至少有一个着火邻居tree neighbor_fire 0就是“树且邻域有火”的格子。rand(N) p和rand(N) f分别给每个空格和每棵树独立抽一次概率保证同步更新。注意new_grid在长树之后又用于闪电判断所以闪电只能点燃本轮之前已经存在的树和本轮新长出的树若不想让新树被闪电点燃可以把闪电判断挪到长树之前。3.3 主循环与 imagesc 可视化主循环负责迭代、调用step_fire、刷新图像。imagesc把整数矩阵映射成颜色colormap定义三行颜色分别对应 0、1、2。caxis固定颜色范围否则 MATLAB 会根据当前矩阵最大值自动缩放火少的时候颜色会跳。drawnow limitrate比每步drawnow快适合长时间动画。若要做 matlab画图 导出可以在循环里抓帧写进VideoWriter。figure; colormap([0.92 0.92 0.92; 0.10 0.55 0.15; 1.00 0.15 0.10]); % 空地/树/火 caxis([0 2]); axis equal tight; axis off; T 3000; for t 1:T [grid, burned] step_fire(grid, p, f); imagesc(grid); caxis([0 2]); title(sprintf(t%d, burned%d, t, burned)); drawnow limitrate; endcolormap的三行分别对应 0、1、2顺序不能反。caxis([0 2])把颜色映射固定在 0 到 2保证空地、树、火颜色稳定。title里显示当前步和燃烧格数方便观察火事件。若火事件很少burned大多数时候是 0可以只在burned 0时打印或保存帧避免生成大量无意义图像。T取 3000 到 10000取决于 f 的大小f 越小需要越长的等待时间才能看到闪电点燃。4. 跑参数扫描从树密度、燃烧面积找临界点4.1 统计量设计单次动画只能看个热闹要判断临界行为得设计统计量。最常用的三个稳态树密度rho_tree即非暂态阶段树格数占 N² 的比例平均火事件大小mean_fire即每次burned 0时燃烧格数的均值火事件频率event_rate即着火步数占总步数的比例。树密度反映系统积累了多少燃料平均火事件大小反映燃料连通程度火事件频率反映点火概率和燃料恢复速度的平衡。把这三个量放在同一张表里比只看动画有信息量得多。统计时要舍弃前一段暂态。初始树密度和稳态树密度可能差很多前几百步的数据会污染均值。常见做法是前 20% 步数不统计只统计后 80%。如果 f 非常小火事件本身就很稀疏需要把 T 拉长到几万步或者把 N 调大否则平均火事件大小会很不稳定。4.2 参数扫描脚本下面脚本扫描 p固定 f、N 和 T每个 p 跑一次。tree_sum累加每步树格数fire_area记录每次火事件的燃烧格数fire_events记录火事件次数。最后把结果整理成表格。若要更稳的统计可以把每个 p 重复 5 到 10 次换随机种子后取平均MATLAB 的parfor可以并行加速前提是安装了 Parallel Computing Toolbox没有的话用普通for也能跑只是慢一些。p_list 0.005:0.005:0.05; f 1e-5; N 120; T 8000; burn_in round(0.2 * T); results zeros(numel(p_list), 4); for k 1:numel(p_list) p p_list(k); rng(100 k); grid double(rand(N) 0.5); grid(randi(N*N)) 2; tree_sum 0; fire_area []; fire_events 0; for t 1:T [grid, burned] step_fire(grid, p, f); if t burn_in tree_sum tree_sum sum(grid(:) 1); if burned 0 fire_area(end 1) burned; fire_events fire_events 1; end end end rho_tree tree_sum / (T - burn_in) / N / N; mean_fire mean(fire_area); event_rate fire_events / (T - burn_in); results(k, :) [p, rho_tree, mean_fire, event_rate]; end Tbl array2table(results, VariableNames, ... {p, rho_tree, mean_fire, event_rate}); disp(Tbl);p_list是扫描的生长概率。burn_in控制暂态步数这里取总步数的 20%。tree_sum只累加暂态之后的树格数fire_area只记录burned 0的步避免大量 0 值拉低均值。event_rate是火事件步数除以统计步数。array2table把矩阵转成带列名的表格方便直接看和导出。rng(100 k)让每个 p 的随机种子不同避免同一随机序列影响不同参数。4.3 结果可视化与相变观察把表格画成图能直观看到趋势。用subplot把树密度、平均火事件大小、火事件频率画在三张子图上横轴都是 p。树密度随 p 增大而上升平均火事件大小在某个 p 附近开始快速增大说明燃料连成了大片火事件频率可能先升后降因为树多了火容易烧但烧完空地恢复也需要时间。这个快速增大的位置就是临界区的粗略信号。注意不要把它当成精确相变点有限尺寸和随机性会让曲线平滑。figure; subplot(3,1,1); plot(Tbl.p, Tbl.rho_tree, -o, LineWidth, 1.2); ylabel(稳态树密度); grid on; subplot(3,1,2); plot(Tbl.p, Tbl.mean_fire, -s, LineWidth, 1.2); ylabel(平均火事件大小); grid on; subplot(3,1,3); plot(Tbl.p, Tbl.event_rate, -^, LineWidth, 1.2); xlabel(生长概率 p); ylabel(火事件频率); grid on;-o、-s、-^分别用圆圈、方块、三角标记数据点。LineWidth加粗线条grid on打开网格。三张图共享横轴 p方便对齐观察。若平均火事件大小出现长尾可以改用对数纵轴set(gca, YScale, log)。若要做更严格的临界分析需要统计火事件大小分布并看它是否服从幂律那就要把每次火事件的burned全部保存下来而不是只存均值。保存时用fire_area的完整数组每个 p 存一个元胞或 CSV后续再拟合。观察量随 p 增大的趋势可能含义稳态树密度上升空地恢复快燃料积累多平均火事件大小临界区附近快速上升树连通性增强火容易连片火事件频率先升后降或趋于平稳点火概率与燃料恢复的平衡最大火事件偶尔出现接近 N²系统接近自组织临界5. 加速、验证与几个容易踩的坑5.1 向量化之外的加速技巧conv2已经把邻域计算向量化了真正的瓶颈常在rand(N)和逻辑索引。每步生成两个 N×N 随机矩阵T10000、N300 时开销不小。可以只在需要的位置生成随机数先找出空格索引再对索引向量抽rand(numel(empty),1) p树的位置同理。这样随机数数量从 N² 降到实际空格数或树数稀疏林地时提升明显。另一个技巧是把imagesc刷新频率降低比如每 10 步画一次计算和绘图分开统计结果不受影响。5.2 用守恒量验证更新逻辑森林火灾模型没有严格守恒量但可以检查树格数变化树增加数等于空地长树数树减少数等于本步被点燃的树数。把sum(new_grid1) - sum(grid1)和sum(empty grow) - burned对比若不等说明更新顺序有重叠比如先长树再判断点燃导致同一格被重复处理。另一个验证是关闭闪电f0系统最终会接近全树火事件消失再把 p0、f0系统只剩初始火蔓延烧完即停。这两个极端能快速暴露规则写错。5.3 三个高频坑第一个坑是colormap和caxis不匹配导致火显示成绿色或树显示成红色。三行 colormap 必须对应 0、1、2并且每帧都设caxis([0 2])否则 MATLAB 自动缩放会让颜色跳。第二个坑是conv2边界零填充让边缘火提前熄灭若做周期边界记得手动补边再卷积。第三个坑是 f 太小导致长时间看不到火误以为代码卡死先把 f 临时调到 1e-3 验证规则再调回 1e-5 跑正式实验。导出动画时用VideoWriter写 MP4帧率取 10 到 20配合drawnow limitrate既能看到火蔓延又不会拖慢主循环。本文还有配套的精品资源点击获取
