简介基于Python的元胞自动机病毒传染模拟SEIR模型可视化项目面向对复杂系统建模、流行病仿真及数据可视化感兴趣的Python学习者。项目以元胞自动机为空间演化框架结合SEIR四状态易感、暴露、感染、康复转移规则实现病毒在二维网格上的动态传播模拟并通过Matplotlib输出时间序列图和热力图帮助理解传染过程及干预策略的影响。压缩包共3个文件封装为一个zip压缩包、一个SEIR.py源码脚本和一份README.md说明文档整体仅8KB结构紧凑。源码脚本涵盖数据处理、状态更新与绘图逻辑文档则对模型参数和运行方式作简要说明便于直接运行与二次修改。目前已有142人学习下载适合正在接触元胞自动机或SEIR模型的读者作为入门范例。通过调参可观察不同感染率、恢复率下的传播差异既能巩固Numpy矩阵操作与随机模拟技巧也能为后续探讨社交距离、免疫策略等实际问题提供可视化参照。1. 用元胞自动机模拟病毒传染SEIR 模型不该只是一条光滑曲线拿一份常规的 SEIR 流行病学报告你看到的是“易感、潜伏、感染、康复”四条随时间变化的曲线光滑、确定、像教科书一样干净。但真实世界的传染从来不是均匀混合的小区里哪栋楼先爆、学校食堂为什么成为超级传播点、一个路口封住之后疫情走向如何曲线都给不出答案。标题里的这套方案本质是把 SEIR 模型的每个“人”按到二维网格的格子中让传染只沿着邻居关系扩散再用 python 把网格渲染成动态可视化。它适合正在做课程设计、论文仿真或者想完成一次完整 python 数据分析与可视化实践的人。跑通不难难在参数口径和可视化直觉这篇就是照着这个顺序讲。2. SEIR 模型与元胞自动机把四条状态转换拆到网格上2.1 SEIR 四状态为什么中间要夹一个“暴露 E”SEIR 和更常见的 SIR 模型最核心的区别就是多了 ExposedE暴露态。在传统 SIR 里易感者 S 接触感染者 I 后直接变成 I相当于潜伏期为零SEIR 则要求 S 先进入 E 状态在 E 状态里“潜伏但不具备传染性”或传染性很弱经过一段时间后才真正成为 I。这个 E 状态不是形式主义它直接影响曲线形态。潜伏期的存在会让感染峰值来得更晚、更低、更平缓。做模拟的时候如果你发现自己的输出曲线比真实疫情数据“尖”很多、峰值时间靠前十有八九是 E 状态的转换概率设得不对或者直接用 SIR 模型替代了 SEIR。对防控仿真来说E 状态还代表“带毒但未发病的窗口期”隔离措施的模拟、核酸检测策略的讨论都得依赖这个中间态。在代码实现里四个状态就是四个整数。我一般习惯这样约定S 0E 1I 2R 3。这里有个小技巧——用 0-3 的连续整数而不是任意常量后面写 numpy 数组统计、画颜色映射时会省很多事。# 状态编码常量全项目统一引用 S, E, I, R 0, 1, 2, 3你可能觉得一行代码没什么信息量但它是一个项目中所有数组操作的基础。后面做 np.sum(grid I) 统计感染人数、用 ListedColormap 给四种状态配颜色都依赖这个统一编码。如果在这里用了 1、5、9、20 这种离散大数逻辑上没错但代码会难读得多。2.2 元胞自动机承载 SEIR状态、邻居和时间步元胞自动机CA的三个基本要素元胞网格、状态集合、邻居规则与时间步。把这个框架和 SEIR 对照网格上的每个格子是一个“人”或一个“人群单元”每个格子持有 S/E/I/R 四种状态之一时间步推进时格子根据邻居的状态按概率更新自己。这里最影响模拟效果的是邻居规则的选择。常用的有两类Von Neumann 邻域只取上下左右四个邻居Moore 邻域取周围 3×3 范围内除自身外的 8 个格子。Moore 邻域模拟的是“广场式混合”人在局部空间里会和四面八方的人接触Von Neumann 邻域更接近“廊道式传播”比如一条走廊两侧的房间、一条街道上的门面房。做一个直觉判断同样的感染概率 β 和康复概率 γ 下用 Moore 邻域模拟出的疫情扩散速度会显著快于 Von Neumann 邻域因为每个 S 状态格子面临的潜在传染源数量差了一倍。如果你模拟的是学校教室、商场这样的人员密集场所应该选 Moore如果你模拟的是宿舍楼道这种线性空间Von Neumann 更贴近。选错邻居规则调参怎么调都救不回来。2.3 为什么不用欧拉法解微分方程均匀混合假设的代价有人会问直接用 scipy.integrate 解 SEIR 微分方程不是更简单吗确实更简单但那是建立在“均匀混合”假设上的——每个易感者与每个感染者的接触概率完全相等。这在宏观统计上能输出一条合理的曲线却完全无法回答“空间上疫情是怎么蔓延开的”。CA 的优势在于空间异质性可以设置网格不同区域的人口密度不同可以让某个区域的格子状态被强制隔离可以观察感染热点如何在网格上移动。标题既然强调“可视化”那必然是希望看到空间上的爆发过程而不只是几个人口比例曲线。这也是我在做这类模拟时的选型逻辑要趋势、要参数敏感性微分方程够了要看空间、看干预、看局部爆发上元胞自动机。3. 网格初始化与疫情地图从空网格到第一帧感染点3.1 用 numpy 初始化网格状态分布与随机种子初始化是整个模拟的地基。常见做法是先生成一个全 S 的二维网格随机撒少量 I再按比例预置一部分 R 来代表疫苗接种或既往感染获得的免疫力。import numpy as np S, E, I, R 0, 1, 2, 3 GRID_SIZE 200 # 200x200 网格共 40000 个元胞 INIT_INFECTED 5 # 初始感染者数量 INIT_IMMUNE_RATE 0.03 # 初始免疫比例模拟疫苗覆盖率 SEED 20240101 # 调试时固定随机种子保证结果可复现 rng np.random.default_rng(SEED) # 1. 全网格初始化为易感状态 S grid np.full((GRID_SIZE, GRID_SIZE), S, dtypenp.int8) # 2. 先按比例随机设置免疫人群 R immune_mask rng.random((GRID_SIZE, GRID_SIZE)) INIT_IMMUNE_RATE grid[immune_mask] R # 3. 在剩余 S 状态格子中随机抽取初始感染者避免覆盖免疫格 candidates np.argwhere(grid S) chosen candidates[rng.choice(len(candidates), sizeINIT_INFECTED, replaceFalse)] grid[chosen[:, 0], chosen[:, 1]] I print(易感人数:, np.sum(grid S)) print(感染人数:, np.sum(grid I)) print(免疫人数:, np.sum(grid R))这段代码有两个细节值得注意。第一状态数组用 np.int8 而不是默认的 int64四万个格子在模拟几百步后会产生大量中间数组int8 能省不少内存读写也更快。第二免疫人群的置入必须先于初始感染者——如果先随机撒感染者再覆盖免疫运气不好时初始感染者会被免疫掩码改成 R你的疫情从一开始就少了一个传染源而且这种 bug 很难肉眼发现。如果你还没装 numpy终端执行 pip install numpy 即可。跑课程设计或论文仿真时固定 SEED 非常重要调试阶段随机种子不固定前后两次运行状态完全不同你就没法判断代码改对了没有。批量做统计分析时再放开种子。3.2 Moore 邻居索引与环形边界向量化统计感染邻居CA 迭代中最频繁的操作是“统计每个格子周围有多少感染者”。最笨的方法是对每个格子写双重 for 循环查邻居200×200 网格一次要查 4 万个格子每一步模拟都要扫一遍跑 300 步就是 1200 万次邻居查询纯 Python 循环会慢到让你怀疑人生。常见的高效做法是利用 numpy 的 roll 函数把整个网格沿行和列各平移一次然后累加比较结果一次性算出所有格子的感染邻居数。这就是向量化的核心思路不逐格处理而是让 numpy 对整块数组做批处理。from numpy import roll def count_infected_neighbors(grid): 统计每个格子周围 8 个邻居中感染者的数量。 count np.zeros_like(grid, dtypenp.int8) for di in (-1, 0, 1): for dj in (-1, 0, 1): if di 0 and dj 0: continue # roll 会把网格整体平移被移出的边界从对面补进来等价于环形边界 count (roll(grid, shiftdi, axis0) I) count (roll(grid, shiftdj, axis1) I) return count要注意循环里每次 roll 只平移单个轴但实际需要的是同时平移两个轴的 8 个偏移。上面的实现先用 roll 沿行平移 di再沿列平移 dj每次调用 roll 的结果和 grid 直接比较统计进入 count。最后 count 数组里每个格子的值就是它周边 8 个邻居中 I 的数量。环形边界是一个常见选择网格最上方的格子它的上方邻居是最下方的格子模拟的是一个“处处等效”的无边界空间。如果你的场景是一栋被封闭的宿舍楼边界外不该有任何邻居这时就得改用 np.pad 给网格包一圈“墙”或者对越界索引做屏蔽。选择哪种边界取决于你想模拟“开放城市”还是“封闭楼栋”没有绝对的对错。3.3 初始感染者放置策略聚集爆发还是分散多点初始感染者怎么放对短时间内的传播形态影响很大。如果 5 个初始感染者集中在网格的一小块区域疫情会在局部先形成密集传播第一波峰值来得早但传播路径相对清晰如果 5 个点随机撒在整张网格各处初期会出现多个独立的小火苗然后同时扩散峰值更宽。最直接的做法是 random.choice 均匀撒点这也是大多数模板代码的默认方式。但如果你想模拟“某栋楼率先爆发”的场景可以在初始化时让感染者落在某个以 (x, y) 为中心的小矩形或圆形区域内。实现并不复杂生成一个初始感染坐标的列表然后循环将对应格子置为 I。我一般建议先跑均匀撒点版本等模型行为正常了再去改聚集分布。初始感染者位置造成的差异本质上属于随机性的一部分后面做批量多次运行时会被平滑掉但如果你只跑单次那这 5 个点落在哪结果可能天差地别这一点在第 5 章避坑部分还会展开。4. 核心迭代逻辑S→E→I→R 状态转换的三条概率规则4.1 同步更新先读旧网格再写新网格进入迭代后最核心的原则是同步更新——从一个旧网格推导出一个新网格绝不能边遍历边改原网格。如果就地更新一个格子在本轮刚变成 I 后立刻被它的邻居看到相当于在同一步里发生了二次传播模拟出的速度会比现实快好几倍这就是很多初学者跑出来“一天之内全图飘红”的元凶。正确的状态推进规则拆成三条S 格子的感染概率取决于周围感染者数量E 格子按概率转 II 格子按概率转 R。注意 R 是终态不再承担任何后续转换。def ca_step(grid, beta, sigma, gamma): 推进一个时间步。入参 grid 为旧状态返回新状态不修改原数组。 next_grid grid.copy() # 逐步写回 next_grid保证同步更新 n_infected count_infected_neighbors(grid) # 基于旧网格统计 # 规则 1易感者被感染。邻居里有 n 个感染者时 # 被感染概率 1 - (1 - beta)^n即至少发生一次有效接触 exposed_mask (grid S) (rng.random(grid.shape) 1 - (1 - beta) ** n_infected) next_grid[exposed_mask] E # 规则 2暴露者发病进入感染态 infectious_mask (grid E) (rng.random(grid.shape) sigma) next_grid[infectious_mask] I # 规则 3感染者康复或隔离进入终态 R recovered_mask (grid I) (rng.random(grid.shape) gamma) next_grid[recovered_mask] R return next_grid这段代码逻辑上最关键的地方是 rng.random(grid.shape) 生成了一张与网格同形状、值在 0 到 1 之间的随机矩阵然后与转移概率比较。这样做比逐格循环快一个数量级。规则 1 里的感染概率公式 1-(1-beta)^n 值得多说一句如果直接写成 beta * n当 n5、beta0.2 时概率会超过 100%不合常理用 1 减去“每次接触都没传染成功”的概率才是多接触的独立事件概率。这也是很多人对照公式后调参踩坑的地方。4.2 时间步与概率换算β、σ、γ 的参数表初始化参数不能随手填。模型里每个概率都必须对应一个明确的时间尺度。这里给一张我在课程设计和论文仿真里常用的参数取值表适用于“一个时间步 一天”的设定参数含义典型取值计算口径betaβ单次有效接触感染概率0.02 ~ 0.05取决于口罩佩戴、社交距离、室内通风sigmaσ潜伏期转阳概率1/7 ≈ 0.143平均潜伏期 7 天则每日转阳概率为 1/7gammaγ感染康复概率1/14 ≈ 0.071平均传染期 14 天则每日康复概率为 1/14细心的读者会发现平均潜伏期 7 天对应的 σ 是 1/7但概率版模型里 E 状态每天以 1/7 概率转 I这意味着实际潜伏期分布是指数型的最短 1 天长则十几天。如果你想严格固定潜伏期正好 7 天就不能用概率σ而要给 E 格子加一个“剩余潜伏期计数器”每步减 1减到 0 再转 I。两种做法各有适用场景简单概率版适合快速跑整体趋势计数器版适合更严谨地复现“隔离 14 天”这类政策效应。这里最大的坑是“一步不等于一天”。如果你希望加快模拟速度让一个时间步代表半天那 σ 就要从 1/7 调整为 1/14γ 调整为 1/28反过来如果你把 300 步理解成“300 天”但其实一步只有半天时间轴就整整少了一半。做可视化时坐标轴上的“天”必须和代码里的参数换算严格对应否则看图的人会得到一个完全错误的疫情周期判断。4.3 加干预口罩、疫苗接种和隔离在 CA 上怎么落地元胞自动机最有价值的用途之一就是模拟干预手段。口罩可以抽象成一个全局乘数把有效感染概率从 beta 调整为 beta * mask_efficiency比如口罩使传染风险下降 50%就把 BETA 变成 0.5 * BETA。这不是精确的流行病学建模但做对比实验时能直观看出“全民戴口罩 vs 不戴”对峰值的影响。隔离措施最粗暴的做法是把部分 I 格子直接置为 R代表“被收治隔离不再传染”但这种处理丢失了时间信息。更常见的处理是给隔离中的格子加标记状态仍然是 I但 set 它的邻居计数不参与传播。简单方式是在 ca_step 里单独生成一个隔离掩码矩阵规则 1 统计邻居时把隔离格子的感染邻居排除掉。疫苗接种则可以完全复用初始化时的免疫格子逻辑。把 INIT_IMMUNE_RATE 从 0.03 提到 0.7你就能看到模拟曲线峰值被显著压低——这正是群体免疫的直观演示。做这类对比时我强烈建议把对照组和实验组放进同一次运行流程里批量跑不要靠肉眼比较两次单次运行的曲线。5. 可视化落地与避坑从“能跑”到“能信”5.1 matplotlib 逐帧渲染imshow 与 FuncAnimation可视化是这个项目标题的交付重点。核心思路分成两路一路是动态的网格空间图用 imshow 把整数状态矩阵渲染成彩色格子另一路是 SEIR 四条曲线的时间序列图。两者结合才能既看到空间爆发过程又看到宏观曲线变化。import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation from matplotlib import colors as mcolors # 四种状态对应的颜色白S 黄E 橙I 蓝R cmap mcolors.ListedColormap([#f5f5f5, #ffd92f, #f28e2b, #4e79a7]) fig, ax plt.subplots(figsize(6, 6)) ax.set_title(SEIR 元胞自动机疫情模拟) im ax.imshow(grid, cmapcmap, vmin0, vmax3) fig.tight_layout() def update(frame): global grid grid ca_step(grid, BETA, SIGMA, GAMMA) # 推进一步 im.set_array(grid) # 更新像素数据 return [im] ani FuncAnimation(fig, update, frames300, interval50, blitTrue) ani.save(seir_ca_simulation.gif, writerpillow)interval50 表示每帧间隔 50 毫秒也就是动画以约 20 帧/秒的速度重绘。要注意 FuncAnimation 的 update 函数里必须用 global grid 引用全局变量否则 grid 每帧都会重新从初始状态开始动画就会变成“永远在原地抖动”。另一个高频问题是网格尺寸过大时渲染卡顿200×200 的 imshow 一般没事如果上到 500×500 就要考虑用 AxesImage.set_array 而不是每帧重建 imshow 对象。曲线绘制相对简单跑完模拟后统计每个时间步的 S/E/I/R 数量再做折线图。这里有一个常见 office 场景的坑就是时间轴刻度太密集导致整条 x 轴糊成一团后面避坑部分专门说。5.2 常见坑过程可复现结果才能被信任下面五条踩坑记录按“现象 → 原因 → 解决”的顺序写基本覆盖了第一次做这类模拟会遇到的典型问题。第一条模拟跑到第 3 步全图一片橙红感染人数瞬间见顶。原因是写迭代时对旧网格就地更新了或者新网格又被同一轮的传染源二次读取。解决方法是严格遵循 4.1 节的写法所有邻居统计发生在旧网格上所有状态写入发生在 next_grid 上最后整体替换。写完后可以打印前两步的感染人数变化做校验如果一步之内翻了几倍基本就是同步更新被破坏了。第二条python 画图横坐标太密集。模拟跑 300 天plt.plot 默认在每个整数点都画一个刻度300 个刻度标签挤在一起完全看不清。解决方法是手动设置刻度稀疏度用 plt.xticks(range(0, T, 30)) 或 plt.MaxNLocator(6) 强制只显示少数几个刻度。这个小事在答辩时特别影响观感一个清晰的横轴比什么美化都重要。第三条I 曲线始终为零疫情根本没爆发。原因大概率是 R0 1也就是 β 太小而 γ 太大传播“烧”不起来。先在纸上算一下Moore 邻域下理论接触数为 8R0 ≈ β×8/γ。比如 β0.005、γ0.07 时R0 0.005×8/0.07 ≈ 0.57传染病自己就消失了。这不是 bug而是没调参。第四条R 数量开局就占了很大比例初始感染者反而没有落在预期位置。原因是初始化顺序错了免疫掩码把感染者覆盖掉了。回到 3.1 节的初始化代码先做免疫置入再选择感染者。第五条单次运行结果和上一次完全不一样怀疑代码有随机 bug。这不是 bugCA 本来就是随机过程单次运行只是其中一个样本。调试时固定 random seed 确保复现做结论时批量跑几十次取中位数和置信区间。不这么做你的模型结论就只是抽了一次签。5.3 调参顺序不要一上来就追求“拟合真实数据”调参的第一个原则是从 R0 反推参数而不是瞎试。先根据你模拟的疾病设定目标 R0比如想模拟新冠早期传播R0 在 2.5 附近然后固定 γ 1/14、邻居数 8反解 ββ R0×γ/8 ≈ 0.044。有了这个基准点再往两边微调曲线形态很快就对上了。第二个原则是先调爆发规模再调时间轴最后调空间分布。爆发规模可以通过峰值感染人数占比观察时间轴对应 E 和 I 的持续步数空间分布则看感染热点是集中还是分散。每一步只动一个参数否则出现差异时你根本说不清是哪一个参数导致的。提示调参时建议把当时的参数组合连同结果图一起保存文件名写成“beta003_sigma014_gamma007_t300.png”这样的格式。参数变量没有记录在图片里的话过两周回来看你根本不知道这张图跑的是什么参数这是做仿真最不值当的返工。6. 让模拟可信的进阶技巧批量运行取中位数并用 R0 验证模型6.1 批量跑 32 次画中位数与四分位带单次 CA 运行是一个随机样本曲线震荡大没法用来下结论。常见做法是固定参数批量跑几十次收集每次的感染人数曲线再按时间步取中位数和四分位区间。中位数比均值更适合这类分布因为疫情大爆发会制造很长的右尾几次极端播报会把均值拉高而中位数能代表“典型场景”。def init_grid(size, n_infected, immune_rate, rng): 整合第 3 节初始化逻辑为可复用的工厂函数。 grid np.full((size, size), S, dtypenp.int8) grid[rng.random((size, size)) immune_rate] R # 先免疫 candidates np.argwhere(grid S) ch candidates[rng.choice(len(candidates), n_infected, replaceFalse)] grid[ch[:, 0], ch[:, 1]] I # 后感染 return grid def simulate(beta, sigma, gamma, T, trials32): 批量运行 32 次返回感染者数量矩阵 shape(trials, T)。 curves [] for trial in range(trials): rng np.random.default_rng(SEED trial) # 每次不同种子 grid init_grid(GRID_SIZE, INIT_INFECTED, INIT_IMMUNE_RATE, rng) i_curve [] for _ in range(T): grid ca_step(grid, beta, sigma, gamma) # 注意 count 内部用到 rng 会冲突 i_curve.append(np.sum(grid I)) curves.append(i_curve) return np.array(curves) arr simulate(BETA, SIGMA, GAMMA, T) median np.median(arr, axis0) q25, q75 np.percentile(arr, [25, 75], axis0) plt.plot(median, color#f28e2b, label感染中位数) plt.fill_between(range(T), q25, q75, alpha0.3, color#f28e2b, label25%-75%区间) plt.xlabel(时间步天); plt.ylabel(感染人数); plt.legend(); plt.show()注意上面代码中 ca_step 内部依赖一个全局 rng 变量批量场景里容易串种子。更干净的做法是让 ca_step 接受 rng 作为参数比如 ca_step(grid, beta, sigma, gamma, rng)每次循环传入当前 rng。我在实际代码里都是直接把 rng 加进函数签名的避免全局状态污染。6.2 用 R0 β × 邻居数 ÷ γ 验证模型是否合理参数设得对不对可以用基本再生数 R0 做一次快速校验。经验公式是 R0 β × 有效接触数 ÷ γ。这里的“有效接触数”不是嘴上说的 8而是邻居真正处于感染状态时的接触总数实际中会比 8 小因为邻居里往往是 S 和 R 居多。一个可操作的做法是在运行时统计“I 格子的 S 邻居接触次数”反推真实有效接触数再和手算的 R0 对一下。不同场景的有效接触数参考Von Neumann 邻域理论最大 4Moore 邻域理论最大 8加隔离措施时再乘上隔离比例。隔离比例越高有效接触数越低R0 相应下降。模型接受度的一条底线是只有当有效传播数经过批量仿真后确实低于 1疫情才应当在模拟中稳步消失。如果你的批量结果显示疫情明明熄灭了R0 却算出来是 2.3那么要么公式里接触数估计不准要么代码里隔离掩码没有真正生效。反过来R0 大于 1 而疫情消失了也要回去查代码。我最早做这个方向的时候把 σ 直接当成“潜伏期天数 7”就填进代码导致 E 状态 7 步之后才被激活曲线形态和真实疫情差了十万八千里。后来才意识到参数应该换算成“每天转移概率 1/天数”。更尴尬的是有一版结果我拿单次运行就下了“封校有效”的结论被导师追问“这个结论在 100 次运行里能复现多少次”我才去补了批量与置信区间。现在我的习惯是固定种子只调 bug批量才下结论。希望帮到你。本文还有配套的精品资源点击获取