简介这份PDF文档围绕元胞自动机Cellular Automata在森林火灾传播模拟中的应用展开面向正在学习复杂系统建模、MATLAB仿真或离散动力学的大学生与研究人员也适合作为课程设计、实验报告与课堂演示的参考资料。包内仅含1个PDF文件压缩包约64KB体量轻便便于离线阅读与随时代码对照。文档以二维矩阵S表示森林状态0空地、1绿树、2燃烧系统讲解燃烧、蔓延、生长与闪电起火四条演化规则并给出p0.3、f6e-5等关键参数的设置依据与实现思路还介绍了燃烧树与绿树着色矩阵R、G以及彩色显示矩阵C的图形化处理方法帮助读者把抽象的元胞邻居叠加运算转化为直观的动态画面。借助这份材料读者可以快速理解火势扩散与森林自然恢复的耦合机制掌握用MATLAB循环迭代模拟时空演化的编码套路并进一步探讨生长概率、闪电概率对火灾行为的影响。目前已有1072人学习下载是入门元胞自动机仿真的一份实用参考。1. 从一次蔓延实验说起元胞自动机森林火灾模型在MATLAB里到底算什么一片林子初始时树木随机撒在网格上你每按一次时间步空地随机长出树树以很小的概率被闪电点着火再沿着上、下、左、右四个邻格烧过去。几十步之后屏幕上的绿色会停在某个密度不再下降偶尔一次大火把整片森林烧穿紧接着又慢慢长回来。这件事反直觉的地方在于没有任何参数在变化系统自己就卡在了「随时可能发生雪崩」的边缘。元胞自动机森林火灾模型做的就是把这种自组织临界现象搬到离散网格上用最少的规则复现最像真实火灾的统计规律。适合的人群也很明确——学复杂系统、做元胞自动机入门、需要一套能跑能改的 MATLAB 代码的人。标题里的三个词元胞自动机是方法森林火灾模型是对象MATLAB 代码是落地形式接下来就按这三层拆开讲透。2. 森林火灾模型的元胞自动机规则设计与状态编码2.1 三种状态的元胞定义与邻域选择森林火灾模型的状态空间很小典型做法只用三个整数表示一个格子当前是什么0 代表空地1 代表活着的树2 代表正在燃烧。所有格子在同一时刻更新这就是元胞自动机里说的同步更新。邻域的选择直接决定火势传播的形状常见两种冯诺依曼邻域是上下左右四格摩尔邻域是含对角线在内的八格。四邻域烧出来的火边界比较规整八邻域扩散更快、更接近「风吹火」的观感。写代码时我一般先用四邻域规则好推导、统计量也好解释需要时再换成八邻域。状态编码含义下一步的转移条件0空地以概率 f 长成树(1)否则保持空地1树邻域有火(2) 或以概率 p 被闪电击中则变为火(2)否则保持2火恒定在下一步变为空地(0)即烧完熄灭注意状态 2 只存在一个时间步这是模型不积累「烧焦」状态、维持可解析性的关键别顺手写成「火烧若干步才熄灭」。2.2 着火概率 p 与树木生长概率 f 的角色模型里有两个核心参数一个管「毁」一个管「生」。着火概率 p是每棵树在每个时间步被闪电单独点燃的概率取值通常很小量级在 1e-4 到 1e-1 之间生长概率 f是每块空地每步长出树苗的概率量级同样在 1e-3 到 1e-1。两个参数的比例关系决定了系统最终停在哪p 相对 f 越大火越频繁稳态森林密度越低p 相对 f 越小树越容易长满攒到一定程度后一次大火烧光反而出现典型的临界雪崩。很多人第一次做会直觉认为 p 越大烧得越干净但如果你把 p 从 0 慢慢往上加会发现森林密度并不是线性掉下去而是在某个区间突然塌方这正是自组织临界在参数扫描里的影子。还有一类变体把「闪电引燃概率」和「火势蔓延概率」分开树被点着用 p邻格火传给树用另一个概率 q。绝大多数 MATLAB 实现里 q 直接取 1让火一步传一格简化掉一个自由度。你要复现文献里的曲线第一步就得确认对方用的是哪一套规则别拿单参数模型去对双参数模型的结果。2.3 边界条件开放边界与环形边界边界看起来是小细节却会显著影响大火的规模分布。开放边界下网格四周是「悬崖」火到边上就没了传播方向边界附近的树等同于被保护会让统计出的火烧面积偏小。环形边界则是把上下、左右卷起来格子坐标对 L 取模整片森林没有一个「边缘」这正是统计临界指数时想要的无限大平面近似。MATLAB 里实现环形边界最省事的写法不是 if 判断而是用 mod 或索引拼接% 环形边界下的四邻域火源L 为网格边长 idx (x) mod(x-1, L) 1; % 把坐标映射回 1..L spread (G(idx(rr-1), cc) 2) | (G(idx(rr1), cc) 2) | ... (G(rr, idx(cc-1)) 2) | (G(rr, idx(cc1)) 2);逻辑说明mod把第 0 行映射到第 L 行、第 L1 行映射到第 1 行等价于把平面卷成圆环rr、cc是网格坐标矩阵。参数说明L必须和网格实际尺寸一致改网格尺寸时别漏改这里的 L否则边界会错位一格统计出来的密度会带上周期性的假象。3. MATLAB代码实现从网格初始化到同步更新3.1 参数设定与网格初始化先定几个一眼能改的参数再生成初始森林。初始条件用「随机撒一半树」比「全空」更快进入稳态因为全空要等很久才有结构。固定随机种子是个好习惯否则每次跑出来的曲线都不一样排查问题时抓不住对照。function G initForest(L, treeFrac, seed) % 初始化森林网格 % L 网格边长 % treeFrac 初始树木占比0~1 % seed 随机种子固定后结果可复现 rng(seed); G zeros(L, L); G(rand(L) treeFrac) 1; % 1 表示树0 表示空地 end逻辑说明zeros建全空网格rand(L) treeFrac生成逻辑矩阵把满足条件的格子置 1。参数说明treeFrac取 0.4 到 0.6 之间收敛较快取 0 就是从光地开始长适合观察生长阶段。3.2 同步更新的向量化写法森林火灾模型的更新必须整批算、最后一次性赋值这就是同步更新的精髓。逐格 for 循环不仅慢还会因为「先改的格子影响后改的格子」引入顺序依赖破坏模型的统计性质。MATLAB 里我用掩码加逻辑索引来实现function G stepForest(G, p, f) % 单步同步更新 % p 每棵树被闪电点燃的概率 % f 每块空地长树的概率 L size(G, 1); burning (G 2); % 当前燃烧的格子 % 四邻域火势传播受火来源是原网格的燃烧状态 spread false(L); spread(2:end,:) spread(2:end,:) | burning(1:end-1,:); spread(1:end-1,:) spread(1:end-1,:) | burning(2:end,:); spread(:,2:end) spread(:,2:end) | burning(:,1:end-1); spread(:,1:end-1) spread(:,1:end-1) | burning(:,2:end); light (G 1) (rand(L) p); % 闪电随机点燃 grow (G 0) (rand(L) f); % 空地长树 newG G; % 先复制再统一改 newG(burning) 0; % 火一步后熄灭 newG(G 1 spread) 2; % 邻火引燃 newG(light) 2; % 闪电引燃 newG(grow) 1; % 空地长树 G newG; end逻辑说明spread全部基于旧的burning计算保证一步内火只从旧火源传一格三条赋值语句分别处理熄灭、点燃、生长全都作用在G的副本newG上最后整体覆盖杜绝顺序依赖。参数说明rand(L)每调用一次消耗一段随机数别把light和grow写成同一个随机矩阵否则点燃和生长会被人为绑定p越大整体越干f越大恢复越快。一个绕不开的实现细节这里用的是开放边界spread(2:end,:)这类写法会让第一行和最后一行不会从「对侧」接收火源。想要环形边界就按 2.3 节的mod索引重写传播部分。3.3 可视化与动画录制跑起来看不见等于没跑。用image加自定义色图最省事空地浅灰、树深绿、火红色三色分明。参数建议值作用L100 ~ 200网格边长200 以上单步明显变慢p1e-3 ~ 1e-1闪电引燃概率f1e-3 ~ 1e-1空地生长概率steps300 ~ 1000观察稳态需要的步数frameStep1 ~ 5每隔几步渲染一帧cmap [0.88 0.88 0.88; % 0 空地 0.13 0.55 0.13; % 1 树 0.85 0.15 0.10]; % 2 火 figure; for t 1:steps G stepForest(G, p, f); if mod(t, frameStep) 0 image(G 1); % 1 是因为 colormap 索引从 1 开始 colormap(cmap); axis equal off; title(sprintf(step %d, t)); drawnow; end end逻辑说明G 1把 0/1/2 映射成 1/2/3 去索引色图行drawnow触发实时刷新。参数说明frameStep越大动画越流畅但细节越少录制 GIF 时把它设成 1 或 2用getframe逐帧存下来更稳。4. 参数扫描与临界行为的MATLAB实验4.1 着火概率p对稳态森林密度的影响单跑一次只能看到一条时间曲线想做参数扫描就得把「跑一步」封装成「跑全程并返回统计量」。我关心的量有两个稳态森林密度最后若干步树占格子的比例和最大单次火烧面积。前者描述整体稀疏程度后者描述雪崩规模。忽略前 20% 的步数作为预热只统计后半段的平均能避开初始条件带来的偏差。4.2 批量参数扫描脚本把主循环包成函数外面套一层 p 的循环就能扫出一条密度-p 曲线。function [dens, maxBurn] runForestFire(L, p, f, steps, warmup) % 返回稳态森林密度与最大火烧面积 % warmup 预热步数不参与统计 G initForest(L, 0.5, 42); treeCount zeros(steps, 1); maxBurn 0; for t 1:steps nBurnBefore sum(G(:) 2); G stepForest(G, p, f); treeCount(t) sum(G(:) 1) / L^2; if t warmup maxBurn max(maxBurn, nBurnBefore); end end dens mean(treeCount(warmup1:end)); endpList 0.005:0.005:0.2; dens zeros(size(pList)); for k 1:numel(pList) dens(k) runForestFire(150, pList(k), 0.05, 600, 200); fprintf(p%.3f density%.4f\n, pList(k), dens(k)); end plot(pList, dens, -o); xlabel(p); ylabel(steady density);逻辑说明runForestFire内部固定种子、固定初始条件只让 p 变化保证变量隔离warmup之后才累计最大火烧面积。参数说明pList步长取 0.005 够看清拐点f固定成 0.05 是为了让 p/f 比例成为唯一自变量steps和warmup要一起放大否则小 p 时系统还没进入稳态就开始统计曲线会整体偏低。实际跑出来密度会随 p 增大而单调下降但在中小 p 区间下降很缓说明系统在用「攒树—大火—清空—再攒」的方式自我维持把 p 再往上推森林被频繁的小火反复打断几乎长不起来。这条曲线的形状就是判断代码规则写没写对的第一个信号。4.3 常见坑随机数种子、更新顺序、边界三个坑几乎每次都会绊住人。一个是更新顺序改成逐格 for 循环后结果变样因为火在同一个时间步内沿着循环方向连烧了好几格实际把时间步拉长了。第二个是随机数种子忘了rng或者每条曲线用不同种子对比出来的差异可能全是噪声。第三个是边界开放边界会让靠近边缘的火提前熄灭统计火烧面积偏小如果论文里要求临界指数必须换成环形边界。还有内存坑L 从 100 提到 500网格从 1 万格涨到 25 万格rand(L)每步都新建矩阵递归调用几百步下来会明显变慢此时该把rand(L)拆成按需生成或者用randi一次性给足随机数。5. 进阶用守恒量与分布函数验证森林火灾模型的结果跑出图像并不等于模型对了得拿可观测量去卡。常用的验证手段是看火烧面积的分布函数临界状态下单次火烧掉 s 个格子的概率 P(s) 近似服从幂律 P(s) ~ s^(-τ)在双对数坐标下画出来是一段直线。用 MATLAB 统计的方法很简单在runForestFire里每步记录「上一步烧掉的格子数」跑完之后做直方图取对数拟合斜率。sizes sizes(sizes 0); % 去掉没着火的步 edges logspace(0, log10(max(sizes)), 30); counts histcounts(sizes, edges); centers sqrt(edges(1:end-1) .* edges(2:end)); pS counts / sum(counts) ./ diff(edges); loglog(centers, pS, o); pfit polyfit(log(centers(centers5)), log(pS(centers5)), 1); fprintf(拟合指数 tau %.3f\n, -pfit(1));逻辑说明histcounts在对数分箱里统计频数centers取几何中心避免对数坐标下偏左polyfit只对中段剔除小面积噪声和大面积截断做线性拟合斜率取负就是幂指数。参数说明分箱数 30 是经验值太少拟合不稳太多末端箱子只有一两个样本拟合区间下限建议取 5 到 10上限贴近max(sizes)之前因为有限网格在大火烧穿边界处必然偏离幂律。想再往前走一步可以固定 f、只调 p观察拟合出的 τ 是否稳定稳定说明系统处在临界区随 p 乱跳说明还没到位或者网格太小。另一个省事但有效的检查是守恒量监测总面积固定为 L²空地、树、火三者占比之和恒为 1把三条曲线画在一张图上看是否始终相加为 1能立刻发现状态转移写漏了分支。把这几步做完你手里的森林火灾模型就不只是会动的彩色方块而是能对得上量纲、扛得住参数扰动的可复现实验平台。本文还有配套的精品资源点击获取