简介：「元胞自动机—森林火灾模型MATLAB代码.pdf」面向学习复杂系统模拟与元胞自动机建模的高校学生及科研入门者，用一份可直接运行的MATLAB脚本演示森林火灾的起火、蔓延与自然恢复过程。资源包仅含1个PDF文件，压缩后约64KB，内容为带注释的完整代码清单。代码以二维矩阵表示森林状态，0为空地、1为绿树、2为燃烧树木，按四条规则迭代演化：燃烧树转为空地、绿树被相邻火点引燃、空地以概率p=0.3长出新树、绿树以概率f=6e-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 里循环慢，N=200、迭代几千步时差距很明显。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.7；p 取 0.001 到 0.05；f 取 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、2；p和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 随机矩阵，T=10000、N=300 时开销不小。可以只在需要的位置生成随机数：先找出空格索引，再对索引向量抽rand(numel(empty),1) < p，树的位置同理。这样随机数数量从 N² 降到实际空格数或树数，稀疏林地时提升明显。另一个技巧是把imagesc刷新频率降低，比如每 10 步画一次，计算和绘图分开，统计结果不受影响。
5.2 用守恒量验证更新逻辑
森林火灾模型没有严格守恒量，但可以检查树格数变化：树增加数等于空地长树数，树减少数等于本步被点燃的树数。把sum(new_grid==1) - sum(grid==1)和sum(empty & grow) - burned对比，若不等，说明更新顺序有重叠，比如先长树再判断点燃导致同一格被重复处理。另一个验证是关闭闪电f=0，系统最终会接近全树，火事件消失；再把 p=0、f=0，系统只剩初始火蔓延，烧完即停。这两个极端能快速暴露规则写错。
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，既能看到火蔓延，又不会拖慢主循环。
本文还有配套的精品资源，点击获取
