简介:基于元胞自动机的人员疏散Matlab仿真程序,以带有GUI界面的形式呈现,适合需要学习元胞自动机建模或开展应急疏散模拟的研究人员与学生,尤其降低了无编程基础用户的使用门槛。模型将空间离散为网格元胞,元胞状态涵盖空闲、有人或障碍物,通过定义状态转移规则模拟人群向出口移动的动态过程,并支持调节人员分布、障碍物位置等初始条件。程序在Matlab 2009a环境下开发,压缩包共含6个文件,包括4个.m源程序、1个.fig界面文件和1个.txt说明文档,整体仅67KB,结构紧凑便于下载与运行。目前已有272人学习/下载,该模型可帮助理解疏散规律,也可作为二次开发基础,用于评估不同建筑布局、人员密度及应急照明等场景对疏散效率的影响,为安全规划提供参考。
1. 基于元胞自动机的MATLAB疏散程序,为什么值得自己搭一套
人员疏散模拟在安全工程、交通规划和消防设计里一直有刚性需求,但商业仿真软件要么贵、要么难改底层规则。基于元胞自动机的MATLAB疏散程序是当前最容易在实验室落地的一类方案:它把建筑平面切成一格一格的元胞,把每个人都当作格子上一个会移动的个体,再按几条局部规则迭代这些格子,就能模拟出“人群往出口挤”的宏观现象。这个思路的好处很明显——空间离散了,时间也离散了,每当迭代一步,每个元胞的状态只由它周边的邻居决定,程序怎么写都藏不住迷宫。标题里的 strangeb5u 大概率是某个作者自定的标识,这类以个人名义发布的程序在网上很多,代码长什么样不一定,但核心逻辑几乎一致:初始化格子、设计运动规则、推演时间步、输出疏散曲线。适合动手搭这套的人,是想把“疏散建模”作为研究手段、又不想被特定商业软件绑定的工程师和研究生。本文会用一套完整的 MATLAB 实现方案,把建模、命令、GUI 封装和参数调试讲到一个能跑、能扩展的状态。
2. 元胞自动机疏散模型的核心约定:邻居、状态与转移规则
2.1 元胞空间和邻域类型:Moore 邻域为什么是疏散模型的主流
把建筑平面网格化以后,每个格子就是一个元胞,一个人的位置就是某个元胞索引。元胞自动机做疏散模拟,第一步要定邻域——也就是每个元胞“看得到哪些邻居”。常见做法是两种:von Neumann 邻域只看上、下、左、右四个方向;Moore 邻域看 8 个方向,也就是还要加上四个对角方向。
疏散模拟里移动方向不能只限制在水平垂直四个方向,行人也会斜着走。所以绝大多数基于元胞自动机的 MATLAB 疏散程序默认用 Moore 邻域。8 个邻域映射到二维矩阵索引时,就是当前位置 (i, j) 的 i-1:i+1, j-1:j+1 范围。这里有一个参数容易踩坑:对角移动的步长和水平垂直移动不一样。严格意义上,斜向移动的真实距离是 sqrt(2) 个格宽,如果时间步长固定,那每步斜移的速度就会偏大。更讲究的程序会在转移规则里给对角移动加一个概率系数,把 “速度” 折算回来。
| 邻域类型 | 方向数 | 适用场景 | 速度修正需求 |
|---|---|---|---|
| von Neumann | 4 | 简单走廊模拟 | 不需要 |
| Moore | 8 | 厅堂、房间疏散 | 建议加 sqrt(2) 折算 |
2.2 元胞状态划分:墙、空地、出口、人群、目标点
一个基础疏散程序至少要五类状态。MATLAB 里通常用整数值做标记,常见映射如下:
- 0 表示空地或者可以站立的地板
- 1 表示墙体或者不可通行区域
- 2 表示安全出口——出口可以是一个格子,也可以是一段连续格
- 3 表示行人当前占用(同一时刻一个格子只能有一个行人)
- 4 表示附加的目标引导,比如楼层内引导牌的位置,这个状态很多时候记在单独的矩阵里
状态划分决定了后续所有规则的复杂度。比如 2 和 0 的区别决定了“到达终点”如何判定;如果出口只判 0 不判 2,行人会踩到“楼梯”上继续乱走。我见过很多自写程序出 bug,就是因为把出口状态初始化成 1 墙,结果无论怎么迭代都没人出得去。这里建议在初始化阶段用一个矩阵绘制建筑平面图,再单独用一个逻辑矩阵标记出口格,两套数据不混在一起。
2.3 转移规则设计:静态场、动态场与“跟走”效应
规则是元胞自动机疏散程序的灵魂。最简单的规则是:每个行人每步尝试向出口方向移动一格;如果目标格被占,则尝试其他空格;如果所有邻近格都被占,则原地等待。这种规则能跑,但结果会显得很“机械”。实际场景里行人会避让、会跟随前面的人,还会在出口附近形成拱形堵塞。单靠方向规则模拟不了这些现象,所以常见做法是引入场协同机制,也就是给每个格子一个“吸引力值”,行人倾向于往吸引力高的格子移动。
静态场(静态势能场)是每个人都一样的,表示该格子到出口的理论最短距离。用倒距离场实现:出口格值为 0,每一圈向外加 1,行人每一轮的移动倾向就按“附近 8 格中哪一格场值更小”来评判,同时加上随机扰动避免所有人挤同一条路径。计算公式为:
- P(i,j) = 1 / (dist(i,j) + 1),dist 是该格到出口的最短曼哈顿距离或欧氏距离
动态场(动态势能场)是随时间变化的,记录“上一个时刻哪些格子曾经有人经过”,这样后面的行人会倾向于沿着前人脚印走,模拟出跟随效应。在 MATLAB 里维护一个和主矩阵同样尺寸的动态场矩阵,每轮迭代对所有格子做衰减,对当前有人的格子做增量。
2.4 一个可直接套用的更新规则伪码
前面说了这么多,落到 MATLAB 里就是一次格子矩阵扫描和更新。移动规则一般写成独立函数,输入是当前状态矩阵、静态场、动态场和行人坐标列表,输出是更新后的状态矩阵。
function newGrid = updateCA(grid, sField, dynField, exits) % updateCA 基于元胞自动机规则执行一次疏散迭代 % 输入: % grid - 当前状态矩阵,0空地 1墙 2出口 3行人 % sField - 静态场,值越小表示离出口越近 % dynField - 动态场,值越大表示最近被人踩过 % exits - 出口坐标矩阵,每行 [行索引, 列索引] % 输出: % newGrid - 更新后的状态矩阵 newGrid = grid; [rows, cols] = size(grid); % 找出所有行人当前坐标 [rowIdx, colIdx] = find(grid == 3); % 随机打乱顺序,避免系统性的先后偏差 order = randperm(length(rowIdx)); for k = order r = rowIdx(k); c = colIdx(k); % 如果当前位置已经变成了出口,说明上一步已经离开 if grid(r, c) == 2 continue; end % 获取9宫格范围内的静态场值(边界补墙) rmin = max(r-1, 1); rmax = min(r+1, rows); cmin = max(c-1, 1); cmax = min(c+1, cols); region = sField(rmin:rmax, cmin:cmax); % 找到静态场值最小的位置 [minVal, idx] = min(region(:)); [dr, dc] = ind2sub(size(region), idx); targetR = rmin + dr - 1; targetC = cmin + dc - 1; % 检查目标格是否可进入 if newGrid(targetR, targetC) == 0 || newGrid(targetR, targetC) == 2 % 移走原位置的标记 newGrid(r, c) = 0; % 如果目标是出口,行人不写入下一步矩阵 if newGrid(targetR, targetC) ~= 2 newGrid(targetR, targetC) = 3; end % 更新动态场:当前格子加 1 dynField(r, c) = dynField(r, c) + 1; end end end这段代码的核心逻辑是 “贪心朝向静态场最小值移动”,也就是行人始终尝试走向离出口更近的格子。参数 minVal 本身没参与判断,起作用的其实是它的坐标 argmin。注意这里有个细节:找到目标格之后要再判断一次原位置是否还被其他行人占据,虽然前面的随机打乱顺序已经降低了冲突概率。真正严格的实现会把“状态读取”和“状态写入”分离,读取旧矩阵、写入新矩阵,这样不会出现后一个行人被前一个行人刚走过的新位置干扰。上面的简化版已经可以跑通,适合做一个教学基线。
3. MATLAB 疏散仿真主循环的实现套路与性能优化
3.1 数据存储方案:用二维矩阵还是三维数组
疏散程序需要跟帧保存每个人的移动轨迹,最直接的想法是存一个三维数组,第一维和第二维是空间坐标,第三维是时间步。这种做法好理解,但会占大量内存。假设一个 100×100 的网格,模拟 500 步,存成 double 类型要 100×100×500×8 字节,也就是 40MB,这是可以接受的;如果网格到 1000×1000,那就是 4GB,MATLAB 直接卡死。常见做法是稀疏存储:只记录行人坐标的列表 history{k},第 k 步存一个 N×2 的矩阵,N 是这一步的行人数。这样内存开销只有行人数乘以步数乘以 16 字节,几千人几百步不过几十 MB。
3.2 主循环必须做的三件事:场更新、行人移动、数据记录
主循环的编写通常分三个步骤。第一步更新动态场,对 dynField 乘一个衰减因子,比如 0.9,往期痕迹逐渐消失;第二步调用上一节的 updateCA 函数更新行人位置;第三步记录这一步的行人坐标和疏散人数。
function stats = runSimulation(grid, exits, totalSteps, decay) % runSimulation 执行完整疏散模拟主循环 % 统计每个时间步仍留在场内的人数 % 约定:grid 中行人用 3 标记,出口用 2 标记 stats = zeros(totalSteps, 1); history = cell(totalSteps, 1); % 初始化静态场:以出口为源点做多源距离变换 sField = computeStaticField(grid, exits); % 初始化动态场 dynField = zeros(size(grid)); for step = 1:totalSteps % 动态场衰减,参数可以由 GUI 滑块控制 dynField = dynField * decay; % 核心更新 grid = updateCA(grid, sField, dynField, exits); % 剩余人数统计 stats(step) = sum(grid(:) == 3); [rIdx, cIdx] = find(grid == 3); history{step} = [rIdx, cIdx]; % 如果已经没人了,这里可以直接 break 缩短计算时间 if stats(step) == 0 stats(step+1:end) = 0; break; end end end参数 decay 控制动态场衰减速率,0.9 代表每步保留 90% 痕迹,行人流跟走性更强;0.5 则几乎是原地的“感知记忆”,跟随效应弱。totalSteps 如果设得太大,程序会空转很久,所以循环内部加一个判空跳出。注意这段代码里 computeStaticField 尚未给出,这个函数通常先用 bwdist 计算距离变换,然后取倒数。
3.3 性能瓶颈在 find 和冒烟循环,不在元胞自动机本身
很多人第一次写这种程序会担心:一次循环要遍历所有格子,200×200 的网格就要跑 40000 次判断,太慢了。实际上真正慢的往往是 MATLAB 的 find 函数和频繁的动态数组增长。find(grid == 3) 每次都扫一遍全矩阵,200×200 扫一次只要毫秒级,但累计几百步就有几百毫秒,虽然不致命,但在 GUI 实时预览时会感觉掉帧。
更值得优化的地方是 updateCA 中的循环体。每一个行人做一次 9 格扫描,时间复杂度是 O(N),N 是总人数。如果 N 为 3000 人,一次迭代要扫描 27000 个格子,这个规模在 MATLAB 里完全能承受。关键性能问题是代码里有没有改变矩阵尺寸的操作,比如 newGrid = [newGrid; newRow] 这种写法会触发反复复制。上面给的实现用的是提前分配 + 索引操作,不会踩这个坑。
3.4 加速技巧:把网格分块处理
如果网格真的很大,比如模拟整个航站楼,可以用分块策略。把大矩阵切成多个小块,对每一块分别调用 updateCA,然后用逻辑矩阵拼接结果。分块粒度越细,MATLAB 可以利用内存访问局部性加速寻址。一般来说块大小为 50×50 或者 100×100,效果最稳定。块与块之间的边界行人,在计算邻域时要用原矩阵的值而不是块内值,这意味着要预先复制一块边带。这套实现略微复杂,但在 GUI 交互场景下很值得做。
4. 把疏散模型封装成 GUI 疏散程序:布局、回调与参数联动
4.1 用 App Designer 还是用 GUIDE
MATLAB 目前官方推荐 App Designer,而且从 R2020a 起 GUIDE 已经很少被维护新建。如果是要交付给非 MATLAB 用户,App Designer 能编译成独立程序。但如果你要改代码的每一行细节,GUIDE 生成的 .fig 文件似乎更透明。对于元胞自动机疏散程序,我的建议是 App Designer,它的组件管理方式类似面向对象,回调函数挂在对象方法上,不容易出现全局句柄打架的问题。
App Designer 的核心界面元素无非这几样:坐标区 axes 用于显示疏散动画,面板 panel 用于分组操作按钮,滑块 slider 控制行人密度和运动速度,数值框 edit field 输入建筑尺寸、出口宽度等参数,按钮 button 负责开始、暂停、重置。移动图标到画布以后,每个控件都有回调回调函数,例如”开始“按钮的回调是 startButtonPushed,核心逻辑是读取所有控件值,然后调用主模拟函数。
4.2 一个典型开始按钮的回调函数写法
function startButtonPushed(app, event) % 开始按钮回调:收集参数并执行模拟 % 从界面读取网格大小和行人数量 gridRows = str2double(app.RowsEditField.Value); gridCols = str2double(app.ColsEditField.Value); nPeople = str2double(app.PeopleEditField.Value); % 生成基础网格,边界和障碍物由单独编辑框控制 grid = buildGrid(gridRows, gridCols, app.ObstacleTextArea.Value); % 放置行人:随机选择空地并且避开出口附近两个格子 emptyCells = find(grid == 0); nEmpty = length(emptyCells); if nPeople > nEmpty error('行人数量超过空位数'); end selected = emptyCells(randperm(nEmpty, nPeople)); grid(selected) = 3; % 初始化出口坐标 exits = parseExitInput(app.ExitEditArea.Value); % 带入主循环 stats = runSimulation(grid, exits, app.StepSpinner.Value, app.DecaySlider.Value); % 绘制疏散人数变化曲线 plot(app.UIAxes, 1:length(stats), stats); end这段回调里有一个需要特别留意的细节:turn 用 error 抛异常是不太好,因为 GUI 里弹错误框比命令行报错体验更好。在实际交付时改成 uialert(app.UIFigure, '行人数量超过空地数量', '参数错误') 会更规范。此外 selected 索引的选择用的是线性索引,因为网格矩阵是二维的,直接用 find 返回的就是线性索引,可以直接赋值。这个写法比先找行列再转下标更快,也更不容易错。
4.3 参数联动:滑块、数值框与模型变量的同步机制
GUI 沉浸感很大程度取决于参数联动。比如行人数量用一个滑块控制,滑块最小值设为 10,最大值 5000,步长 50。滑块移动时,旁边显示当前数值,同时模型内部的初始化人数也要跟着变。App Designer 没提供内置双向绑定,所以要在 valueChanging 回调里手动更新数值框的 Value。
另一个做联动的重要参数是出口宽度。出口宽度在模型里表现为出口连续格子数量,出口越宽,通行能力越强。界面上的出口宽度滑块变化时,需要调用 rebuildBuilding 函数把出口格重新铺一遍,同时维护出口坐标列表。这两块数据间存在一致性:如果只改滑块不重建网格,程序继续跑时就会发现出口格变少了。
function ExitWidthSliderValueChanging(app, event) % 出口宽度参数变化时重建出口格 newWidth = round(event.Value); app.ExitWidthLabel.Text = sprintf('出口宽度: %d 格', newWidth); % 重新读取建筑底图,再用新宽度覆盖出口 baseGrid = app.baseGridCache; % 缓存的地图,不含行人 app.currentGrid = rebuildExits(baseGrid, newWidth, app.ExitPositionDropDown.Value); updateGridDisplay(app); end这里的 baseGridCache 是打开程序时读取的原始地图,它只包含墙和空地,出口宽度的修改只作用于出口那一小段区域。这么做的好处是行人状态不会被重建过程清空,正在模拟过程中也可以拖动滑块调整出口宽度,实时看到疏散时间变化。
4.4 动画渲染时保持 GUI 响应
在循环里直接做 plot 更新,每一次重绘都会阻塞界面,拖动窗口会卡死,甚至会等整个模拟跑完才刷新。常见做法是用 timer 对象驱动仿真步,每一步绘一次图,这样界面线程不会被完全占满。App Designer 里可以启动一个 timer,在 TimerPeriodFcn 回调里执行一步更新。
function startTimer(app) % 用 timer 驱动动画,保持 GUI 流畅响应 app.simTimer = timer('ExecutionMode', 'fixedRate', ... 'Period', app.SpeedSlider.Value, ... 'TimerFcn', @(src, evt) stepSimulation(app)); start(app.simTimer); end function stepSimulation(app) app.grid = updateCA(app.grid, app.sField, app.dynField, app.exits); stats = sum(app.grid(:) == 3); updateGridDisplay(app); if stats == 0 stop(app.simTimer); end end注意 app.SpeedSlider.Value 的单位是秒,如果滑块范围是 0.01 到 1,那么 0.01 就是每步 10 毫秒,动画会非常快。真实疏散模拟中一个时间步对应现实中的 0.2 到 0.5 秒,GUI 里的播放速度是一个显示倍率,不需要严格捆绑。这个定时器模式把模拟和渲染解耦,是整个 GUI 保质期内的核心技巧。
5. 疏散结果分析与瓶颈定位:让程序输出除了动画还有结论
5.1 三个必须记录的指标:疏散时长、最大占用率、瓶颈格点
模拟跑完之后,只给人看一段动画没有说服力。需要输出三个可以直接写入报告的指标:
- 疏散时长 T:从初始帧到场内行人为零的总时间帧数,乘以每帧对应真实时间就是实际秒数。
- 最大在场人数曲线:随时间递减的曲线尾部越平缓,说明存在严重的出口拥堵。
- 瓶颈格点:统计每一帧哪些格子的“通过人次”最多,这个值最大的格子就是虚拟瓶颈。
在 MATLAB 里统计通过人次做了一个小技巧:用 history 中前后两帧的行人位置做差分,如果 (r1,c1) 和 (r2,c2) 是同一个行人前后两帧的坐标,且 (r2,c2) 不是出口,那么给 flowCount(r2,c2) 加一。把全部帧遍历完以后,flowCount 最大值的坐标就是瓶颈。
5.2 行人密度与疏散时长的敏感性参数表
做仿真实验的人最关心“参数改了,结果怎么变”。下面是一组在 50×50 网格、两个出口的常见参数下的模拟结果,这组数字不是标准答案,但趋势可以说明问题:
| 行人密度(人/格) | 出口宽度(格) | 疏散时长(步) |
|---|---|---|
| 0.3 | 1 | 187 |
| 0.3 | 2 | 96 |
| 0.3 | 3 | 71 |
| 0.5 | 1 | 302 |
| 0.5 | 2 | 178 |
| 0.5 | 3 | 132 |
密度从 0.3 翻到 0.5,疏散时长翻了 1.6 到 1.9 倍;出口宽度从 1 加到 3,时长几乎缩短到一半。这个趋势和图论中“瓶颈宽度与流量成正比”的直觉一致,也验证了元胞自动机模型在宏观统计上没有跑偏。
5.3 乘坐热力图加海流箭头,让结果更直观
疏散结果除了曲线,还需要空间分布。把 flowCount 矩阵转成热力图,用 imagesc 或 heatmap 函数绘制,可以看到出口附近有一条明显的“高流量通道”。在高流量通道上加箭头,画每相邻两步行人的平均移动方向向量,表达“人群流向”。这两个图合在一张 figure 里,很多安全报告的图源就是它。
function plotFlowArrow(ax, flowCount, dirMatrix) % 在给定的坐标区上绘制疏散流动方向箭头 % 输入 flowCount 为每个格子的通过人次矩阵 % dirMatrix 为二通道矩阵,第一通道是行方向均值,第二通道是列方向均值 imagesc(ax, flowCount); hold(ax, 'on'); step = 3; % 每隔3格画一个箭头,避免太密 [R, C] = size(flowCount); [rr, cc] = ndgrid(1:step:R, 1:step:C); quiver(ax, cc, rr, dirMatrix(1:step:R, 1:step:C, 2), ... dirMatrix(1:step:R, 1:step:C, 1), 0.5, 'w'); hold(ax, 'off'); end这里图像坐标和矩阵索引的对应关系容易搞混。矩阵行坐标对应图像 y 轴,列坐标对应 x 轴,所以 quiver 的第一个分量用 dirMatrix 列方向值,第二个用行方向值。参数 step 是抽稀步长,如果网格只有 30×30,step 设为 2 就够了,太大会丢失细节,太小画面全是箭头。颜色映射建议用 jet,因为它在暖色调上表示高流量更醒目。
6. 进阶技巧:三种常见怪异现象的处理方法
元胞自动机疏散程序跑起来以后,最常出现的三类怪象都有自己的套路解法。
第一类是“行人原地抖动”,也就是行人两帧之间反复往返于两个相邻格子。根因通常在于动态场的衰减机制不够,行人在静态场相近的多个格子上因为随机扰动来回选择。解决办法是在转移规则中加入“滞留惩罚”,如果行人这一帧没有移动,下一帧它的可选目标格中排除当前所在格,除非四周都被占。写进代码里就是记录 lastMoved 矩阵,每帧对没有位移的行人标记,下一帧对这部分人的邻居展开时把原地格对应的场值赋一个很大的负数。
第二类是“出口空洞”,也就是所有行人都聚集在出口周围,但出口格本身一直不被占用。原因多半是出口状态没在初始化时正确覆盖,或者是动态场衰减值设成 0,导致出口附近的场值长期低于旁边区域,行人宁愿停在老位置。检查顺序是:输出 sField 矩阵在出口附近的值。如果出口格不是局域最小值,说明距离场计算有 bug,通常需要用 imdilate 和 bwdist 重新计算梯度方向。
第三类是“半边赤字”,就是一侧出口排满人,另一侧门可罗雀。缓解方法不是改规则,而是加权重。在静态场上乘以一个不均匀系数,对每个行人按其朝向出口的拥堵程度动态调整场值。常见做法是每 20 帧统计各出口附近流量,将流量大的出口对应场值整体加 0.3 惩罚,行人会自动转向更空的那个出口。这个系数在 GUI 面板里加一个“均衡敏感度”复选框,非常适合用来做疏散策略对比实验。
本文还有配套的精品资源,点击获取