前阵子有个学弟找我问毕业设计题目是“基于元胞自动机的人口疏散模型MATLAB实现”。他最开始的理解特别乐观把房间画成网格人涂成几个格子设定出口然后一运行就能看到人流往门口涌最后做两张曲线图收工。我听完就跟他说画格子确实不难但真正的功夫在三个地方一是网格怎么离散化才合理二是行人每一步“往哪走”的规则怎么定三是仿真结束后怎么从数据里读出堵塞、瓶颈和疏散瓶颈。这篇文章我就把自己做这类建模时的完整思路、MATLAB代码骨架、参数实验和踩坑记录摊开讲适合正在做课程设计、毕业设计或者准备用MATLAB快速验证疏散方案的读者参考。1. 为什么选元胞自动机从宏观排队到微观自组织的思路转变1.1 三类主流疏散建模方法对比在做人群疏散建模之前最好先想清楚一个最基本的问题为什么大家普遍选择元胞自动机而不是直接上流体方程或者社会力模型我用一个简单对比来说明白。方法空间形式典型特征优势劣势连续流体模型连续把人流当作流体场强调密度、速度、流量数学形式成熟适合分析走廊、大门等大尺度通道个体行为被平均化看不出排队、冲突、绕行社会力模型连续每个行人受自驱力、排斥力、墙边界力作用轨迹精细能模拟拥挤、恐慌、身体接触计算量大参数多调参是个体力活元胞自动机离散空间被均匀网格切分行人状态离散按局部规则更新实现简单、速度快能涌现出拱形堵塞、通道振荡等现象精度有限个体轨迹不够平滑我做疏散类的项目很多时候还要跑不同密度、不同出口宽度、不同人群构成下的几十组对照实验。用社会力模型当然更精细但每一步都要解微分方程组扫描参数时计算成本很高。流体模型又抹掉了个体差异恰恰疏散问题最核心的“个体抢位、拥挤堵塞”就看不见了。元胞自动机把空间和时间都离散化用非常简单的局部规则就能复现很多宏观现象是效率和解释力之间比较平衡的选择。1.2 元胞自动机的三个组成要素元胞自动机之所以好上手是因为它只有三个要素需要你定义清楚格子、状态、规则。格子Cell把二维疏散空间切成方形网格每一个格子代表一个可以被行人占据的空间单元。状态State每个格子在同一时刻只能处于一种状态。简化版里通常有四种空地、行人在此、墙/障碍物、出口。规则Rule行人的每一步决策只依赖自己和周围局部邻域内的格子状态不依赖全局信息。你可以把整个过程想象成在棋盘上模拟数百枚棋子移动每个棋子只看得见自己周围两三步的情形它根据“哪边离出口更近”“旁边有没有人挡路”“哪个方向已经被信息素标记过”这几个信息决定下一步。所有棋子同时决策、同时更新久而久之人群的宏观形态就自己长出来了。这种“局部规则 并行更新”的思维方式跟MATLAB特别合拍。MATLAB里二维数组天生就是网格矩阵操作能高效处理状态更新imagesc和绘图函数又能快速把仿真过程变成动态画面很适合做原型验证和教学演示。2. 场景建模与静态场把真实空间翻译成矩阵2.1 网格尺寸与状态编码建模型的第一步不是写代码而是确定“一个格子代表现实里多大一块地”。行人动力学研究里一个成年人站立时的投影面积大约在0.30.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;因为静态场里墙的值为infisfinite能排除墙和不可达区域。行人密度不要设置得过高否则初始化阶段就可能大量重叠出错。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) 越大行人越“目标明确”几乎只会朝出口方向走(\beta0) 时行人完全随机游走像失去方向感的人一样。候选邻居集合要满足几个条件在网格范围内、不是墙、没有其他行人占据、静态场值有限。这里多强调一句我是一个个条件写逻辑判断的不要图省事忽略“行人占据”这个条件否则两个行人在同一时间步会重到同一个格子里。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秒。这个换算关系要在报告里写清楚否则结果没有实际参考价值。第五个坑是只跑一次就下结论。元胞自动机里的随机性很强初始位置、冲突消解顺序都会影响结果。参数对比时要多跑几组取均值不然很可能会被单次实验的偶然性带偏。我在实际做这类项目时最大的感受是不要把精力全花在花哨的动态画面上模型规则本身的合理性才是灵魂。先跑通一个最简单场景再逐步加出口、加障碍、加慢速行人每加一个因素就做一轮对比实验这样出来的结果既有层次也经得起问。如果你按照上面的框架把代码跑通再往里面加入自己的场景约束就能很快得到一个可复现、可解释的疏散分析工具。