元胞自动机实现人群疏散模型：MATLAB仿真全流程解析

前阵子有个学弟找我问毕业设计，题目是“基于元胞自动机的人口疏散模型MATLAB实现”。他最开始的理解特别乐观：把房间画成网格，人涂成几个格子，设定出口，然后一运行就能看到人流往门口涌，最后做两张曲线图收工。我听完就跟他说，画格子确实不难，但真正的功夫在三个地方：一是网格怎么离散化才合理，二是行人每一步“往哪走”的规则怎么定，三是仿真结束后怎么从数据里读出堵塞、瓶颈和疏散瓶颈。这篇文章我就把自己做这类建模时的完整思路、MATLAB代码骨架、参数实验和踩坑记录摊开讲，适合正在做课程设计、毕业设计，或者准备用MATLAB快速验证疏散方案的读者参考。

1. 为什么选元胞自动机：从宏观排队到微观自组织的思路转变

1.1 三类主流疏散建模方法对比

在做人群疏散建模之前，最好先想清楚一个最基本的问题：为什么大家普遍选择元胞自动机，而不是直接上流体方程或者社会力模型？我用一个简单对比来说明白。

我做疏散类的项目，很多时候还要跑不同密度、不同出口宽度、不同人群构成下的几十组对照实验。用社会力模型当然更精细，但每一步都要解微分方程组，扫描参数时计算成本很高。流体模型又抹掉了个体差异，恰恰疏散问题最核心的“个体抢位、拥挤堵塞”就看不见了。元胞自动机把空间和时间都离散化，用非常简单的局部规则就能复现很多宏观现象，是效率和解释力之间比较平衡的选择。

1.2 元胞自动机的三个组成要素

元胞自动机之所以好上手，是因为它只有三个要素需要你定义清楚：格子、状态、规则。

• 格子（Cell）：把二维疏散空间切成方形网格，每一个格子代表一个可以被行人占据的空间单元。
• 状态（State）：每个格子在同一时刻只能处于一种状态。简化版里通常有四种：空地、行人在此、墙/障碍物、出口。
• 规则（Rule）：行人的每一步决策只依赖自己和周围局部邻域内的格子状态，不依赖全局信息。

你可以把整个过程想象成在棋盘上模拟数百枚棋子移动：每个棋子只看得见自己周围两三步的情形，它根据“哪边离出口更近”“旁边有没有人挡路”“哪个方向已经被信息素标记过”这几个信息决定下一步。所有棋子同时决策、同时更新，久而久之，人群的宏观形态就自己长出来了。

这种“局部规则 + 并行更新”的思维方式跟MATLAB特别合拍。MATLAB里二维数组天生就是网格，矩阵操作能高效处理状态更新，imagesc和绘图函数又能快速把仿真过程变成动态画面，很适合做原型验证和教学演示。

2. 场景建模与静态场：把真实空间翻译成矩阵

2.1 网格尺寸与状态编码

建模型的第一步，不是写代码，而是确定“一个格子代表现实里多大一块地”。行人动力学研究里，一个成年人站立时的投影面积大约在0.3～0.5平方米，因此很多文献把网格边长取在0.4米左右。也就是说，一个格子差不多站一个人，密度高时人挨着人，密度低时格子空着，比较符合直觉。

网格尺寸一旦定下来，就要把真实场景翻译成MATLAB能处理的矩阵。我习惯用这几个数值表示状态：

• 0：空地，行人可以进入；
• -1或者inf标记：墙/障碍物，行人绝对不可穿越；
• 1：行人占据；
• 2：出口格子。

用数值标记而不是建立复杂的对象结构，原因很简单：MATLAB矩阵操作快，而且后面算静态场、检测邻居状态时，直接做逻辑比较就行了。比如判断某个格子能不能走，就要同时满足不在墙内、没有被行人占用、距离场值不是无穷大三个条件，拼成一个逻辑表达式就可以向量化处理。

2.2 静态场必须绕开障碍物

如果说网格是疏散模型的“舞台”，那么静态场就是行人的“指南针”。静态场通常用S(i,j)表示：从每个格子到最近出口的行走距离。出口处的 S 值设为0，离出口越远的格子 S 值越大。行人移动时总倾向于走到 S 值更小的格子，这样就产生了“向出口方向走”的吸引力。

这里有一个特别容易翻车的坑：很多人图省事，直接用欧氏距离算静态场。在一个没有任何隔断的开阔房间里，欧氏距离没问题；但一旦场景里有墙、有拐弯，欧氏距离会直接穿过墙体，把墙另一侧计算成“很近”，行人就会表现出穿墙而过的诡异行为。所以静态场必须考虑障碍物的绕行距离。

最简单的绕障碍物计算方法是广度优先搜索（BFS），从所有出口格子出发，逐层向周围扩展。BFS在离散网格上的语义很直接：步数就是距离，墙不能进入。下面是我常用的函数，输入墙掩膜和出口掩膜，输出每个非墙格子的静态场值：

function S = computeStaticField(M, N, wallMask, exitMask) % 计算带障碍物的静态场 % wallMask: 1表示墙, 0表示可通行 % exitMask: 1表示出口, 0表示非出口 INF_VAL = inf; S = zeros(M, N); S(:) = INF_VAL; queue = zeros(M * N, 2); head = 1; tail = 1; % 多出口同时作为BFS起点 [ex, ey] = find(exitMask); for k = 1:length(ex) S(ex(k), ey(k)) = 0; queue(tail, :) = [ex(k), ey(k)]; tail = tail + 1; end dirs = [1 0; -1 0; 0 1; 0 -1; 1 1; 1 -1; -1 1; -1 -1]; while head < tail pos = queue(head, :); head = head + 1; for d = 1:size(dirs, 1) ni = pos(1) + dirs(d, 1); nj = pos(2) + dirs(d, 2); if ni >= 1 && ni <= M && nj >= 1 && nj <= N ... && ~wallMask(ni, nj) && isinf(S(ni, nj)) S(ni, nj) = S(pos(1), pos(2)) + 1; queue(tail, :) = [ni, nj]; tail = tail + 1; end end end end

注意这里采用的是8邻域扩展，也就是人可以向上下左右和四个对角线方向移动。实际场景中如果只允许4方向行走，把dirs改成上下左右四个方向就行。BFS在多出口场景下特别方便，多个出口同时出发，每个格子记录的是到其中任意一个出口的最近步数。

2.3 一个小房间场景的初始化

写一个最简单的算例：60×60的网格，左墙中央和右墙中央各开一个出口，屋子里随机分布行人。初始化代码大致如下：

M = 60; N = 60; wallMask = zeros(M, N); exitMask = zeros(M, N); % 四周做墙 wallMask(1, :) = 1; wallMask(M, :) = 1; wallMask(:, 1) = 1; wallMask(:, N) = 1; % 左右出口：墙中间各留3个格子 exitMask(30, 1:3) = 1; exitMask(30, N-2:N) = 1;

这里有一件事必须提前说：出口宽度在网格里就是几个格子。3个格子大约对应现实1.2米宽，这比单人通过的窄门要宽。疏散模型最后算出来的“疏散时间”对出口宽度极其敏感，所以出口设置要贴近真实场景，而不是随便填数字。

接下来给区域内随机放行人。放行人的时候一定要避免重叠，我通常会先取出所有可用的空地索引，再从中随机抽取指定数量的位置：

rho = 0.4; % 行人密度 ped = zeros(M, N); canWalk = ~wallMask & ~exitMask & isfinite(computeStaticField(M, N, wallMask, exitMask)); freeIdx = find(canWalk); numPed = round(rho * length(freeIdx)); selIdx = freeIdx(randperm(length(freeIdx), numPed)); ped(selIdx) = 1;

因为静态场里墙的值为inf，isfinite能排除墙和不可达区域。行人密度不要设置得过高，否则初始化阶段就可能大量重叠出错。

3. 行人怎么走：移动规则与冲突消解的细节

3.1 从静态场到转移概率

元胞自动机疏散模型最核心的部分，是决定行人每一步移动到哪个格子。最经典的做法是场域模型：行人处在当前格子 i 时，会对周围可通行的邻居格子 j 计算一个“吸引力权重”，然后按照权重随机选择。权重公式之一如下：

[ P_{ij} = \frac{\exp\left(-\beta (S_j - S_i)\right)}{\sum_{k \in neighbors} \exp\left(-\beta (S_k - S_i)\right)} ]

其中 (S_i)、(S_j) 是静态场值，(\beta) 是灵敏度系数。这个公式的意思很直观：如果邻居格子比当前格子更靠近出口，即 (S_j < S_i)，那么指数内部是正的，权重就会大于1，被选中的概率更高。(\beta) 越大，行人越“目标明确”，几乎只会朝出口方向走；(\beta=0) 时，行人完全随机游走，像失去方向感的人一样。

候选邻居集合要满足几个条件：在网格范围内、不是墙、没有其他行人占据、静态场值有限。这里多强调一句：我是一个个条件写逻辑判断的，不要图省事忽略“行人占据”这个条件，否则两个行人在同一时间步会重到同一个格子里。

3.2 并行更新与冲突消解

疏散模拟必须采用并行更新，而不是串行更新。所谓并行更新，是指所有行人在同一时间步先各自想好自己的目标格子，然后统一执行移动。如果按顺序一个一个人更新，前面的人先动了，后面的人看到的世界已经变了，结果就会产生“我先抢到、你被迫等”的人为偏差。

并行更新自然会引出一个问题：两个行人都想进同一个空格子怎么办？这就叫冲突。处理冲突最简单的方法：在每一轮移动时，把行人的处理顺序随机打乱，先处理的人获得该格子，后处理的人发现目标被占就只能留在原地。MATLAB里用randperm打乱索引即可。

更严格的做法是先把所有目标位置相同的行人收集到一个集合里，再让集合内部按概率随机选取一个赢家。我的经验是，如果只是课程设计或方案验证，随机顺序法已经能体现出堵塞涌现；但如果你想发论文，最好做严格版本。

3.3 速度差异与恐慌参数

现实里不是所有人都走得一样快。老人、儿童、行动不便者速度明显慢于普通成年人。在元胞自动机里，表示速度差异有个常用技巧：给每个行人一个移动概率，只有随机数小于这个概率时才尝试移动。普通人取0.9，慢速行人取0.5，这样同样一个时间步内，慢速行人移动次数更少。

恐慌程度怎么体现？可以通过设置不同的 (\beta) 值来模拟。(\beta) 低的人表现得慌乱、乱跑，(\beta) 高的人目标明确。还可以给模型增加一个“从众项”：如果某个邻居格子上有信息素，行人会更倾向于跟着走，这部分属于动态场（Dynamic Floor Field）的范畴。动态场的思路是行人走过的地方留下一定“痕迹”，痕迹随时间挥发，形成对后续行人的吸引。它会自然产生出口前“先到的人走掉、后来的人跟着旧路”的现象。我建议初学者先把静态场模型跑通，再去扩展动态场，否则两个场叠加在一起，出了问题很难定位。

4. 主仿真循环的MATLAB骨架

4.1 参数区与行人初始化

下面是完整的仿真主循环骨架，包含参数设置、初始化和每个时间步的更新。你可以直接复制后根据自己场景改参数。

% ============ 参数设置 ============ M = 60; % 网格行数 N = 60; % 网格列数 beta = 3; % 场域灵敏度 p_move = 0.9; % 正常行人移动概率 p_slow = 0.5; % 慢速行人移动概率 slowRatio = 0.2; % 慢速行人比例 simTime = 500; % 最大仿真步数 visStep = 5; % 每隔几帧绘制一次 rng(42); % 保证可复现 % 墙、出口、静态场 wallMask = zeros(M, N); exitMask = zeros(M, N); wallMask(1, :) = 1; wallMask(M, :) = 1; wallMask(:, 1) = 1; wallMask(:, N) = 1; exitMask(30, 1:3) = 1; exitMask(30, N-2:N) = 1; S = computeStaticField(M, N, wallMask, exitMask); % 行人状态 ped = zeros(M, N); canWalk = ~wallMask & ~exitMask & isfinite(S); freeIdx = find(canWalk); rho = 0.4; numPed = round(rho * length(freeIdx)); selIdx = freeIdx(randperm(length(freeIdx), numPed)); ped(selIdx) = 1;

这里有几个细节值得解释。rng(42)是为了固定随机种子，否则每次运行结果都不一样，后续很难对比参数。行人初始位置不仅避开墙，还要避开出口格子，避免一开始就站在出口上导致疏散时间被低估。慢速行人的标记可以用一个同尺寸的矩阵来存，但要注意它只对“当前时刻是行人”的格子有效：

slowFlag = zeros(M, N); slowFlag(selIdx(rand(length(selIdx), 1) < slowRatio)) = 1;

随机的慢速标记只给行人的初始位置赋值，后面行人移动时，得把这个标记同步搬到新位置。处理方式是在每次移动时，把slowFlag里旧位置的标记一起搬到新位置，再清空旧位置。

4.2 主循环与移动逻辑

history = zeros(1, simTime); nb = [0 0; 1 0; -1 0; 0 1; 0 -1; 1 1; 1 -1; -1 1; -1 -1]; for t = 1:simTime % 移除已经到达出口的行人 arrived = exitMask & (ped == 1); evacuated = sum(arrived(:)); ped(arrived) = 0; history(t) = sum(ped(:) > 0); if history(t) == 0 && t > 2 break; end % 找到当前所有行人 [pr, pc] = find(ped == 1); n = length(pr); goalr = zeros(n, 1); goalc = zeros(n, 1); for k = 1:n i = pr(k); j = pc(k); candR = []; candC = []; candP = []; for d = 1:9 ni = i + nb(d, 1); nj = j + nb(d, 2); if ni >= 1 && ni <= M && nj >= 1 && nj <= N ... && ~wallMask(ni, nj) && ped(ni, nj) == 0 ... && isfinite(S(ni, nj)) candR(end + 1) = ni; candC(end + 1) = nj; candP(end + 1) = exp(-beta * (S(ni, nj) - S(i, j))); end end if isempty(candR) continue; end % 按概率采样 candP = candP / sum(candP); cdf = cumsum(candP); sel = find(rand <= cdf, 1, 'first'); goalr(k) = candR(sel); goalc(k) = candC(sel); end % 并行更新与冲突消解 occupied = ped > 0; order = randperm(n); for k = order if goalr(k) == 0 continue; end gi = goalr(k); gj = goalc(k); if ~occupied(gi, gj) ped(pr(k), pc(k)) = 0; ped(gi, gj) = 1; occupied(gi, gj) = true; % 同步慢速标记 if slowFlag(pr(k), pc(k)) > 0 slowFlag(gi, gj) = 1; slowFlag(pr(k), pc(k)) = 0; end end end % 可视化 if mod(t, visStep) == 0 imagesc(ped + 2 * wallMask + 3 * exitMask); colormap([1 1 1; 0 0 0; 1 0 0; 0 1 0]); % 空地白, 行人黑, 墙红, 出口绿 axis equal tight; title(sprintf('t = %d, 剩余 %d 人', t, sum(ped(:) > 0))); drawnow; end end

主循环里nb的第一项是[0 0]，代表“停留在原地”。这样即使周围没有合适的目标，候选集合里也至少有一个选项，不会出现概率和为零的边界情况。如果你希望行人必须移动，就把这一项去掉。

还有一点容易忽略：出口格子的静态场值为0，所以行人会先走到出口格子上，然后在下一轮开始时被移除。因此疏散曲线里所有逃出人数都有约一个时间步的滞后，这种滞后在宏观统计上基本可以忽略。

4.3 可视化的一些心得

imagesc是最快的可视化方式，但图像比较粗糙。如果想做更直观的动态图，可以在imagesc的基础上叠加plot画出墙体，或者改用scatter画行人点，并按墙体颜色填充背景。做演示时，我会在每帧把剩余人数写到标题里，这样仿真过程本身就可以当一张“疏散曲线实时图”。

另外一个建议：不要在循环里每帧都重算整个imagesc，否则高密度场景会卡到没法看。设置每隔5步或10步绘制一次，速度就能接受。如果想让输出更顺滑，可以配合VideoWriter把每一帧写入视频文件，最终生成一段“疏散过程.mp4”，这是答辩时最能直观展示建模效果的材料。

5. 参数调优与疏散瓶颈：实测中看到的那些现象

5.1 密度升高后疏散时间的变化

把仿真跑起来之后，第一个应该做的事情不是去调参数，而是做几组密度扫描。我建议固定出口宽度、固定 (\beta)，只改变初始行人密度，从0.1一直试到0.8。你会看到疏散时间并不是从低到高线性增加的，而是在某个临界密度附近开始快速上升。

原因是高密度下出口前会出现“拱形结构”：大量行人同时涌向一个出口，大家都想挤进去，结果谁也走不快，反而在出口外形成半圆形堵塞。元胞自动机虽然把人和空间都离散化了，但只要冲突消解规则存在，这个现象就会自然涌现出来。能涌现，才说明模型具备解释力，而不是简单叠加。

如果画出密度-疏散时间曲线，你会发现密度增高到某个阈值时曲线会明显变陡。这个阈值对应的就是现实场景里“出口通行能力达到上限”的时刻。做建筑设计疏散分析时，这个阈值就很有参考价值。

5.2 出口宽度是真正的瓶颈

出口宽度的灵敏度往往比我预想的还要高。把出口从1个格子扩大到3个格子，疏散时间可能会缩短一半以上；从3个格子扩大到5个格子，改善幅度就明显减小。这说明疏散瓶颈并不完全取决于出口面积，而跟行人的排队机制有关。

我在仿真里实测过，当出口只有一个格子时，行人目标选择几乎变成“单行道”，出口前方会形成明显的排队；当出口宽度增加到3格以上后，行人就能并行通过，堵塞消失。做设计时，这个现象提醒我们：单纯加宽出口在一定范围内效果显著，但超过阈值后边际收益会下降。

5.3 灵敏度系数与人员构成的影响

(\beta) 参数对疏散时间的影响有些反直觉。在低密度场景中，(\beta) 从1增加到8，疏散速度提升并不大，因为行人本来就很容易找到出口；在高密度场景中，(\beta) 太高反而会让大家认为“所有方向都必须朝着出口”，在出口前形成严重挤压，堵塞加剧。

慢速行人的比例也值得测一测。我试过把20%的行人设为慢速，初始离散位置随机分布，最后疏散时间比全正常速度时多出15%到25%。如果慢速人恰好分布在出口附近，影响会更明显。这说明模型里的人员构成对结果有很大影响，做报告时应该把慢速人员比例作为不确定度的一部分展示出来，而不是只给一条确定曲线。

6. 结果怎么看：疏散曲线、热力图以及避坑清单

6.1 疏散曲线和平均疏散时间

仿真结束后，最直接的输出是“剩余人数随时间变化”的疏散曲线。前面的代码里我把每步剩余人数存到了history向量里，画出来即可：

figure; plot(1:simTime, history, 'LineWidth', 1.5); xlabel('时间步'); ylabel('剩余人数'); title('疏散曲线'); grid on;

观察疏散曲线时，通常会出现两个阶段：前段曲线下降较快，因为人群分散，大家都能顺利走到出口；后段曲线变得平缓，因为出口已经满了，排队成了主导因素。这两个阶段的分界点，对应的就是出口开始形成稳定队列的时刻。

单纯看一条疏散曲线还不能下结论，因为随机种子不同，结果会波动。我建议对同一组参数跑20次，取平均疏散时间和标准差。平均疏散时间比单次结果可靠得多，标准差也能反映人流组织的稳定性。比如某些参数下20次结果的标准差很大，说明系统对初始位置敏感，这时方案就不够稳健。

6.2 用热力图找“堵点”

想要知道疏散过程中哪里拥堵最严重，有一个非常直观的技巧：额外维护一个累积计数器visitCount，每轮更新时，对当前所有行人所在的格子加1。仿真结束后把visitCount画成热力图，就能看到哪些区域被频繁踩踏。

出口前方通常是热力图最亮的区域，因为所有行人都会汇聚过来；如果场景里有走廊拐角，拐角内侧也会非常亮。热力图的价值在于快速定位设计缺陷，比如某个角落明明有大量行人经过，却没有足够的疏散通道，那说明通道设置可能有问题。

imagesc(visitCount); colorbar; colormap(flipud(hot)); axis equal tight;

6.3 仿真前后的几个常见错误

最后列几个我踩过的坑，每一条都对应真实调试时间。

第一个坑是初始行人重叠。随机生成行人位置时如果没有做好去重，会导致多个行人在同一格子里，后续移动逻辑会出各种诡异bug。用randperm从可用空地索引中一次性抽取，是最简单的防重叠方案。

第二个坑是墙边行人被“困住”。当行人站在紧贴墙壁的位置时，它的候选邻居里可能有一部分落在墙内，剩余可行方向又全被占用，最后它哪儿也去不了。这不算模型bug，但会让仿真图像看起来像有人“卡死”在墙边。解决方法是在候选集合里始终保留“原地等待”这一项，也就是前面代码里nb包含[0 0]的原因。

第三个坑是静态场穿墙。场景复杂时，用欧氏距离算出来的静态场会让行人直接穿过隔墙。一定要用绕障碍物距离，哪怕场景里只有一道矮墙，也要在初始化阶段就处理好。

第四个坑是时间步与真实时间的对应关系。元胞自动机本身只输出时间步，不是秒数。要跟真实时间对应，通常需要根据行人平均速度反推一个时间步对应的秒数。比如设定0.4米一格、正常步行速度1.2米每秒，那么一个时间步大约对应0.33秒。这个换算关系要在报告里写清楚，否则结果没有实际参考价值。

第五个坑是只跑一次就下结论。元胞自动机里的随机性很强，初始位置、冲突消解顺序都会影响结果。参数对比时要多跑几组取均值，不然很可能会被单次实验的偶然性带偏。

我在实际做这类项目时，最大的感受是不要把精力全花在花哨的动态画面上，模型规则本身的合理性才是灵魂。先跑通一个最简单场景，再逐步加出口、加障碍、加慢速行人，每加一个因素就做一轮对比实验，这样出来的结果既有层次，也经得起问。如果你按照上面的框架把代码跑通，再往里面加入自己的场景约束，就能很快得到一个可复现、可解释的疏散分析工具。