从零用MATLAB实现元胞自动机森林火灾模拟：规则、可视化与参数扫描-北京尧图网络

简介基于元胞自动机的森林火灾模拟MATLAB代码包面向本科、硕士阶段需要理解元胞自动机建模与仿真应用的教研学习者可用于课程设计或算法演示。资源内含2个文件包含一个实现森林火灾传播规则的MATLAB主程序以及对应的效果示意图整体压缩包仅6KB轻量易用便于快速运行与修改参数观察演化过程。代码在MATLAB 2019a环境下开发包含完整的元胞规则实现可通过调整初始化条件观察火势动态蔓延过程并输出可视化图形。目前已有173人学习使用对初学者理解火灾蔓延的随机性、相邻元胞状态更新等核心机制具有直观帮助。通过该代码可掌握元胞自动机的基本建模方法并结合可视化结果分析不同初始条件下的模拟效果适合作为元胞自动机专题的入门参考。1. 把“烧到哪、何时烧、烧多大”变成格子游戏的模型防灾部门拿到一批火点坐标第一反应是想知道再烧 12 小时会到哪。真做山火传播的物理模型要风速场、燃料湿度、坡度坡向、雷击点密度一起上复杂度跟写气候模式差不多。元胞自动机把复杂传播过程压缩成一句工程师能听懂的话每个小格子下一时刻烧不烧只取决于它自己和周围邻居的状态。基于元胞自动机模拟森林火灾核心就是“树、火、空”三种状态加一条带概率的转移规则却能推出蔓延速度、火线形态、过火面积这些宏观指标。这份 MATLAB 代码从零搭建完整模拟器规则怎么设、参数怎么调、结果怎么统计一步到位。适合做算法课设、复杂系统入门也想快速验证“砍出隔离带有没有用”的人。2. 元胞自动机与森林火灾模型的建模思路2.1 为什么森林火灾这道题会落到元胞自动机头上真实山火的物理过程包含热辐射、热对流、飞火、水分蒸发直接求偏微分方程组非常痛苦。而元胞自动机的核心假设是一个单元格的状态演化只依赖有限邻域内的局部状态这正是火蔓延的主要物理事实——一片林子的点燃能力在时间尺度上主要影响邻近格点。森林火灾在元胞自动机里天然匹配因为火势传播本质上是一个离散化、随机化、并行的局部过程。每个单元格不需要知道自己距离火场中心多远只需要看周围有没有火星溅过来这和森林火线的实际行为是一致的。和微分方程模型相比元胞自动机的三个优势很明显不需要推导火焰传播偏微分方程状态规则直接写赋值语句天然支持并行性和局部异质性坡度、植被类型可以挂到每个格点属性里概率参数化容易对齐观测数据比如把 p 理解成“邻居有火时这棵树被点燃的概率”代价是物理保真度有限不能模拟火羽流的辐射换热细节。这类模型的价值在“相对趋势”和“临界行为”不在精确预报。教学和前期预研用它验证蔓延趋势工程评估再上精细模型是常见分工。2.2 森林火灾元胞自动机的三件套格点、邻域与状态任何元胞自动机都要先定义三样东西格点排布、邻域形状、状态值。森林火灾模型使用正方形网格每个格点就是一个最小面积的林地单元。状态我用整数编码状态编码含义颜色映射0空地 / 烧过之后的灰地灰白色1健康树木绿色2正在燃烧红色邻域有两种常用选择。Von Neumann 邻域看上下左右四个方向模拟“火只能沿着共用边蔓延”。Moore 邻域看周围 3×3 的全部 8 个格子模拟飞火和热辐射往斜向扩散。森林火灾中斜向蔓延不能忽略所以我通常用 Moore 邻域。边界处理上固定边界直接把越界邻居视为空地也就是模拟区域边上的格子没有树木资源。别用周期边界——火焰绕场一圈会出现从对面重新点燃的非物理效果。2.3 树、火、空的转移规则与参数对照表森林火灾元胞自动机的经典规则表长这样当前状态下一时刻状态触发条件火2空地0燃烧只持续一个时间步熄灭后变成灰地空地0树1以 g 概率重新长树树1火2邻居有火且随机数小于 p树1火2无邻居火但随机数小于 f雷击等独立火源树1树1否则保持存活注意 f 和 p 的实验区差异极大。f 是全局背景着火率设大了整个森林会在 100 步内到处起火模拟实验完全失真。一般的做法是 f 设在 0.0001 到 0.01 之间让少数几个雷击点在森林里自然发展。p 是火线推进概率p 小于 0.2 时火势经常断层熄灭p 大于 0.5 时火线会被视为强传播表上这三个参数就是你要第一轮做敏感性扫描的对象。% 单格状态转移的伪代码示意 switch old(i, j) case 2 new(i, j) 0; % 火停在当前时间步熄灭 case 0 if rand g, new(i, j) 1; end case 1 if hasFireNeighbor(old, i, j) rand p new(i, j) 2; elseif rand f new(i, j) 2; end end上面这段用于解释规则逻辑真正跑批量实验的函数会把hasFireNeighbor展开成数组运算避免逐格调用子函数。2.4 同步更新最容易踩坑的一环元胞自动机要求所有格点基于“上一时刻”的状态同步更新。如果边遍历边更新火会像推土机一样按循环方向单方向流动而且同一棵树在同一时间步内可能被反复判断多次相当于把火速放大了到原来的好几倍。正确写法一定是“读旧写新”演化前复制一份 old每个格点的判断从 old 上取状态计算结果写入 new结束后用 new 替换当前网格。这个看起来很小的一个赋值动作决定模拟结果的合理性。old grid; new grid; % 先复制再基于 old 修改 new % ... 逐格推演 ... grid new;3. 用 MATLAB 写森林火灾元胞自动机最小代码与可视化3.1 运行环境与数据结构选择这段代码不依赖任何工具箱MATLAB R2016a 及以后版本都能跑。不需要给优化工具箱或深度学习工具箱配置任何内容安装好基础 MATLAB 环境即可。我用rand生成随机数用imagesc做热区图显示核心数据结构就是一个整数矩阵。grid zeros(N, N); % 0空地, 1树, 2火选择始终矩阵有两点理由索引快MATLAB 在连续内存上的整数矩阵切片性能远好于 cell 数组显示直接imagesc接受整数矩阵颜色映射的序号正好对应状态码不需要转成 RGB 三维数组。3.2 最小可运行的 MATLAB 主程序把下面代码保存为forest_fire_ca.m直接运行。代码模拟第 2 章规则表使用 Moore 邻域和同步更新。% forest_fire_ca.m 元胞自动机森林火灾模拟 clear; clc; close all; %% 参数配置 N 150; % 网格边长 p 0.3; % 被邻居点燃概率 f 0.0005; % 每时间步雷击自然着火概率 g 0.03; % 空地长树概率 rho0 0.6; % 初始树木占比 Tmax 400; % 最大迭代步数 %% 初始化随机种树选择一个树格点火 grid zeros(N, N); grid(rand(N, N) rho0) 1; [rows, cols] find(grid 1); if isempty(rows) error(没有树木请增大 rho0); end startIdx randi(length(rows)); grid(rows(startIdx), cols(startIdx)) 2; %% 显示初始化 h imagesc(grid); axis equal tight; colormap([0.9 0.9 0.9; 0 0.5 0; 1 0 0]); % 空地, 树, 火 title(元胞自动机森林火灾模拟); %% 同步演化主循环 for t 1:Tmax old grid; new grid; % 同步更新读 old写 new for i 1:N for j 1:N if old(i, j) 2 new(i, j) 0; elseif old(i, j) 0 if rand g new(i, j) 1; end else % old(i,j) 1检查 Moore 邻域有无火 i2 max(1, i-1):min(N, i1); j2 max(1, j-1):min(N, j1); hasFire any(old(i2, j2) 2, all); if hasFire rand p new(i, j) 2; elseif rand f new(i, j) 2; end end end end grid new; set(h, CData, grid); title(sprintf(步数%d 过火面积%d, t, sum(grid(:) 0, all))); drawnow limitrate; pause(0.02); end代码里几个要点拆开说明。grid(rand(N,N) rho0) 1先按密度种树火源随机选一棵树的格子保证起始点有意义而不是点在一空地。hasFire any(old(i2, j2) 2, all)检查 3×3 邻域内是否存在火。max/min截断让边界格点自动退化为短邻域免去单独的边界条件分支。主循环先old grid; new grid;中间的判断只读取old写入new。这是 2.4 节强调的同步更新删除new中间变量直接改grid火线会沿遍历顺序传播出明显伪影排查这类问题时要先检查这里。用sum(grid(:)0)统计过火面积是因为烧完的格点会变成 0而本来就是空地的格点也会被算入初始空地面积需要扣掉更严谨的指标可以单独维护一个 burned 布尔矩阵。3.3 参数表这几个参数先按参考值跑参数参考范围初始建议物理含义p0.1 ~ 0.60.3热辐射引燃邻居的概率控制火线推进速度f1e-5 ~ 1e-20.0005雷击/未知火源独立点燃单格的概率g0.001 ~ 0.20.03火烧迹地植被恢复速度rho00.3 ~ 0.90.6初始植被覆盖率N50 ~ 500150网格边长决定模拟空间分辨率f 和 g 不在一个量级这是正常的。背景着火率太高会让火场从几十个点同时冒出来f 先保持小量级重点先调 p 和 rho0。3.4 录制动画把模拟结果存成 AVI 或 GIF上面代码用drawnow limitrate做实时预览量级到 300×300 时刷新还能跟上。想要保存结果可以在循环内加录制帧v VideoWriter(fire_sim.avi); open(v); % 主循环开始时插入 frame getframe(gcf); writeVideo(v, frame); % 主循环结束收尾 close(v);GIF 同理用rgb2ind压缩色板[imind, cm] rgb2ind(frame2im(getframe(gcf)), 256); if t 1 imwrite(imind, cm, fire_sim.gif, gif, Loopcount, inf); else imwrite(imind, cm, fire_sim.gif, gif, WriteMode, append); end帧率取决于 pause 设置GIF 文件大小可以控制在 2 到 5 MB。4. 深入模拟风场、双参数扫描与概率标定4.1 加风不对称的 Moore 邻域权重基础模型在无风场景下火线呈近似圆形向外扩散。实际山火常有主风向下风侧火线延烧显著加快火场长轴指向风方向。实现上不改状态转移逻辑只需要在邻域检查时加上风向权重。% 权重矩阵 W对应 Moore 邻域 3x3 W ones(3, 3); W(2, 3) 1.8; % 东风右侧邻居权重加大 W(2, 1) 0.5; % 上风邻居削弱 windDir east; windStrength 1.8;判断格点的燃烧概率时把原来的布尔型hasFire改成加权火场强度subGrid old(i2, j2) 2; % 邻域火标记 subWeight W((2 - (i - i2(1))):(2 (i2(end) - i)), ... (2 - (j - j2(1))):(2 (j2(end) - j))); fireIntensity sum(sum(subGrid .* subWeight)); if fireIntensity 0 rand min(1, p * fireIntensity) new(i, j) 2; elseif rand f new(i, j) 2; end要注意边界格点取不到完整 3×3 邻域这里用fireIntensity两维索引裁剪出一个和subGrid等大的矩阵让边界格点也能安全计算加权和。风向参数没有统一标准一般让下风侧权重大约 1.5 到 2.5上风侧 0.3 到 0.7。权重越大火线沿风拉伸越明显。实测 p 较高情况下windStrength2.0会让火场长轴与短轴比达到 1.8 左右和实际卫星图上看到的火烧迹地轮廓比较接近。4.2 把单次模拟变成批量实验p 与 rho0 双参数扫描单次模拟只能看个案研究“哪些参数配置下火会失控”需要跑三维统计。把主循环包进函数输出过火比例。function [burnedRatio] runFire(p, rho0, f, g, N, Tmax) grid zeros(N, N); grid(rand(N, N) rho0) 1; [~, cols] find(grid 1); if isempty(cols) burnedRatio 0; return; end startIdx randi(numel(cols)); [r, c] find(grid 1); grid(r(startIdx), c(startIdx)) 2; initialTrees sum(grid(:) 1) 1; for t 1:Tmax old grid; new grid; % 省略逐格判断与第 3 章代码相同 % ... grid new; if sum(grid(:) 2, all) 0 break; end end burnedRatio 1 - sum(grid(:) 1, all) / initialTrees; end画出热力图pList 0.1:0.05:0.6; rhoList 0.3:0.05:0.8; B zeros(length(pList), length(rhoList)); for i 1:length(pList) for j 1:length(rhoList) B(i, j) runFire(pList(i), rhoList(j), 0.0005, 0.03, 120, 300); end end imagesc(rhoList, pList, B); xlabel(初始树木占比 rho0); ylabel(引燃概率 p); colorbar;建议先在 10 组参数上跑通再扩大网格。双循环界面必须全套跑完中间断掉前面结果全丢可以边算边存为一个矩阵文件随时保存。双参数扫描在模型分析和答辩里是最常被评价的部分。单次动画只说明“看起来像火灾”参数平面则显示出临界区域rho0 低时火自熄灭rho0 超过某个阈值通常 0.5 到 0.6后偶发火场才会持续扩散到大面积这对应了真实森林中的连通渗流效应。4.3 让 p 对应真实数据用小样方实验完成概率标定模拟参数 p 不是随便拍的常见做法是将 p 对齐到实际的蔓延速度。测得有效风速下火线在 5 分钟内前进了 200 米模拟中一个格点边长若设为 10 米每时间步 1 分钟那么每步火线需要推进 4 个格点。p 和推进速度之间是非线性关系因为邻域内有火的格点有多个从概率角度看多方向引燃会叠加。一般做法是先在 50×50 小网格上用二分法搜索 p让模拟火线的前锋推进速度匹配真实速度然后再把 p 拿到完整模拟区域使用。% 二分标定 p目标是每 20 步火线推进 10 格 target_gps 10 / 20; % 每步推进格点数 p_low 0.05; p_high 0.8; for iter 1:12 p_mid (p_low p_high) / 2; gps runFrontSpeed(p_mid); % 返回每步平均前锋推进格数 if gps target_gps p_high p_mid; else p_low p_mid; end end实际工程中常用这个标定好的 p 值先模拟无风场景再加风做相对比较。想完全复现某场山火的所有蔓延细节元胞自动机不够用要做更细的空间异质性参数比如把每块格子的植被含水率写成独立的燃料湿度矩阵。5. 让模型可信的几个进阶操作5.1 固定边界之外的改进第 3 章代码直接把边界外视为空地这被称为吸收边界。适合障碍物边缘明确的小区域在大尺度森林场景中边界外围仍有树木火在边界处被强行熄灭会让过火面积偏低。改进有两条路线。格点密度较小时让边界外格子继承火场状态相当于外推虚拟邻域更严谨的做法是加缓冲区在模拟区域外多画 N/10 圈格子这些格子的参数与边界内部一致统计时只统计内部区域。缓冲区格点也会被点燃避免边界对火场形态的收缩作用。还有一个常被忽略的点初始火源的位置选择。放在森林中心会得到近似圆形火场放在角落则因为边界截断整体蔓延面积明显变小。做参数对比实验必须固定初始火点坐标或者设置多种种子点取均值否则几个实验的火源位置差异会淹没参数变化带来的影响。5.2 构建小型蒙特卡洛实验单次模拟受随机数影响很大特别是 f 引发的首次点燃位置和时间。正规实验应运行多次实验取统计指标。MATLAB 里的便捷方式是用parfor替换for。Nmc 100; burned zeros(1, Nmc); parfor k 1:Nmc burned(k) runFire(0.3, 0.6, 0.0005, 0.03, 150, 300); end fprintf(平均过火比例: %.3f ± %.3f\n, mean(burned), std(burned));在 runFire 内不要调用 figure 或 imshow这样并行效率不受图形锁限制。随机数发生器在每个 worker 上默认独立不用手工处理种子。统计分析时过火面积分布通常不是正态的强风或高概率 p 下会出现双峰形态要么火刚开始就熄灭小面积要么突破临界点后烧掉大部分区域。这种情况只看均值的意义有限建议直接画直方图。拟合双峰高斯分布可以算出临界概率告诉你说这个参数附近不能外推。5.3 快速检查规则实现是否出错的三个验证实验每写一个模型都要有“对拍”方法元胞自动机没有解析解可对照但边界条件验证很可靠。把 f 设为 0、p 设为 1其他不变。正常实现应该看到火从起点沿连通林地向所有方向推进不可燃的空地树挡住所有传播森林内部变成无法逾越的隔离带。任何一个金地被烧掉基本可判断旧矩阵污染了更新逻辑或邻域截断写错方向。把 g 设为 0、f 保持极小值。模拟结束后所有被烧过的地方长时间保持灰地不会长树。如果出现树重新出现的问题说明空地分支的随机生长条件写错了。最后一个验证火源熄灭后统计过火面积应保持稳定。正确同步更新中面积为火后不会下降因为灰地不会重新长树若你调大了 f灰地也可能被雷击重新点燃这个场景的过火面积就不可统计了所以做保存实验时 f 应保持小量级。这三项检查通过之后再重新跑一批 p0.4、rho00.7、东风权重 1.8 的模拟观察火线长轴方向是否沿风偏转和蔓延速度是否随步数稳定用于进一步校核你的邻域权重矩阵切边索引。本文还有配套的精品资源点击获取

从零用MATLAB实现元胞自动机森林火灾模拟：规则、可视化与参数扫描

相关推荐

合肥美的太阳能故障报修电话|加热慢上门排查|欧米到家服务热线

pandoc YAML 元数据块标量类型解析：数字与布尔值如何映射到 Meta 类型

ComfyUI-Workflows-ZHO：直接加载现成的 ComfyUI 工作流模板

更多资讯

Hudi 并发控制：深入理解乐观并发与多作业写入场景

Hudi 与 Hive/Spark SQL 集成：构建高效湖仓一体解决方案

Agentic Awesome Skills 之 NestJS 后端模式：用 Screaming Architecture 生成工程级 Backend Rules

用DOCXReadWrite在Delphi中轻松读写Word文档

eBPF实战：打造进程级CPU功耗监控与分析方案

Trae CN 连上 TaoToken，项目规则才真正生效

今日更新

深入解析Transformer多头注意力机制与工程优化

OpenClaw 的 Skills 跑学习任务，模型通道改到 TaoToken 通道行不行？

ChatGPT报错Oops, an error occurred! 全链路排查指南

本周精选

深入解析Transformer多头注意力机制与工程优化

OpenClaw 的 Skills 跑学习任务，模型通道改到 TaoToken 通道行不行？

ChatGPT报错Oops, an error occurred! 全链路排查指南

本月热门

基于YOLO的焊缝缺陷检测：从数据标注到工业部署全流程实践

分布式服务的决策记录

创业产品的首版范围

读完文章，想聊聊您的网站？