尧图网页设计 
简介一套基于Python的元胞自动机病毒传染SEIR模型可视化实现面向Python学习者、数据科学初学者及对流行病模拟感兴趣的研究者。项目借助Numpy构建二维元胞矩阵将人群划分为易感者、暴露者、感染者、康复者四类状态通过Matplotlib动态展示病毒传播过程可直观理解局部交互与空间因素对疫情发展的影响。资源包共3个文件约8KB包含一个ZIP压缩包、一个核心Python脚本和一份Markdown说明文档。ZIP内应封装完整项目或依赖内容py文件实现SEIR状态转移与元胞自动机更新规则md文档则提供代码结构说明与使用指引便于快速运行和二次开发。目前已有142人学习下载。实际使用中读者可调整感染率、恢复率、潜伏期等参数观察不同干预策略下感染者曲线与空间扩散热力图的变化从而掌握元胞自动机建模、SEIR模型数值模拟以及Python科学计算可视化的完整流程对课程设计、入门科研或公共卫生策略演示均有实用价值。1. 为什么用元胞自动机跑SEIR空间传播是微分方程看不见的先说明一点SEIR模型原本是一组常微分方程S、E、I、R四类人的人数随时间变化画出来就是几条平滑曲线。但曲线不会告诉你“病毒从哪个街区开始传为什么有些小区一直没事有些地方一夜之间全变红”。如果你要的是直观演示、教学课件或者疫情传播空间特征研究元胞自动机是比微分方程更合适的载体。这个标题里的rar工程做的就是把SEIR状态机放进一个个格子组成的世界里用Python模拟每天谁被传染、谁进入潜伏期、谁恢复再把这一过程用动态图展示出来。我在拿到类似项目时第一步不是看可视化而是先把格子世界的更新规则建立在脑子里。因为模拟类项目最容易出的问题不是代码跑不起来而是“跑起来了但结果一眼假”。这篇文章按“状态机→模拟主循环→可视化→干预机制→踩坑→批处理参数扫描”的顺序展开你照着写一份自己的模拟400行以内就能交付。2. SEIR状态机与元胞更新规则先把格子世界搭起来2.1 四个状态与状态转移条件元胞自动机的核心是“每个格子的明天由它自己和邻居的今天决定”。在SEIR模型里每个格子代表一个人状态只有四种S易感、E暴露/潜伏、I感染、R恢复。这个离散化人和空间的做法和真实感染过程非常像——你不会因为“全市平均感染率”被传染而是因为你身边有人感染了。状态转移有三条规则S遇到邻居I有一定概率变成E。这一步是概率事件不是必然事件。E经过潜伏期后变成I。潜伏期是固定天数比如4天到期自动转感染。I经过恢复期后变成R。恢复期结束自动转恢复并假设获得免疫不再感染。关键点是E到I、I到R是确定性转移S到E是随机转移。如果你把S到E也做成确定规则病毒会像潮水一样扫过整个网格真实疫情反倒没这么整齐。随机性来自两个地方一是邻居中有没有I二是每次接触是否真的“中招”。前者由空间结构决定后者由传染概率控制。2.2 最小可运行的时间步函数下面这段代码是最简实现核心是用两个额外数组记录E和I的剩余天数避免把状态和时间耦合在同一个小格子值里。状态数组只存0、1、2、3倒计时数组另存。import numpy as np # 状态编码 S, E, I, R 0, 1, 2, 3 def count_infected_neighbors(grid): 统计每个格子的8邻域中感染者的数量周期性边界由 np.roll 处理 infected (grid I).astype(np.int8) counts np.zeros_like(infected) for di in (-1, 0, 1): for dj in (-1, 0, 1): if di 0 and dj 0: continue counts np.roll(np.roll(infected, di, axis0), dj, axis1) return counts def simulate_step(grid, inc_left, rec_left, prob_infect0.25, latent_days4, recover_days6, rngNone): new_grid grid.copy() new_inc inc_left.copy() new_rec rec_left.copy() # E 潜伏期结束转为 I to_infectious (grid E) (inc_left 1) new_grid[to_infectious] I new_inc[to_infectious] 0 new_rec[to_infectious] recover_days # I 恢复期结束转为 R to_recovered (grid I) (rec_left 1) new_grid[to_recovered] R new_inc[to_recovered] 0 new_rec[to_recovered] 0 # S 被邻居中的 I 感染 infected_counts count_infected_neighbors(grid) susceptible (grid S) # 每个 I 独立接触时合并为一次骰子见 2.3 说明 exposure_prob 1 - (1 - prob_infect) ** infected_counts if rng is None: rng np.random.default_rng() newly_exposed susceptible (infected_counts 0) (rng.random(grid.shape) exposure_prob) new_grid[newly_exposed] E new_inc[newly_exposed] latent_days # 未发生状态转移的 E/I剩余天数减 1 new_inc[(grid E) ~to_infectious] - 1 new_rec[(grid I) ~to_recovered] - 1 # 统计 S E I R 各有多少人 counts np.bincount(new_grid.ravel(), minlength4) return new_grid, new_inc, new_rec, counts代码里最容易被忽略的是同步更新传染概率计算用的grid是旧状态而不是这一步刚刚变成I的人。如果你用更新后的new_grid去统计感染者会造成同一步内“新感染者立刻又传染别人”一个时间步里病毒跳好几代R0 被严重高估。count_infected_neighbors里用np.roll做的周期边界意味着网格最上面一行和最下面一行是邻居。对于传播模拟来说这避免了“病毒撞到边界就停住”的人为假象。如果你希望边界是隔离区比如一座岛的边缘可以把np.roll换成np.pad固定边界后面第4章会细说。2.3 邻居统计与感染概率的两种写法常见做法有两种区别直接影响疫情走势。写法一每个S格子遍历8个邻居凡遇到I就独立掷一次骰子阈值是prob_infect。当I邻居多时相当于一天内多次接触感染概率是1 - (1 - prob_infect)^n也就是代码里的合并写法。但很多初版代码会写成循环里“只要有任何一个I击中就感染”那么同一个S在本轮可能被连续判定多次结果概率被重复计算曲线会偏爆炸。写法二只判断“邻居是否有I”有I就以prob_infect概率感染。这种方式简单但把多个I邻居和一个I邻居同等对待低估了高密度感染区的传播速度。我建议用合并概率1 - (1 - p)^n它介于上面两种极端之间参数含义也更接近流行病学里的“单位时间内每次有效接触的传播概率”。如果你想跑得快可以用卷积替代循环。scipy.ndimage.convolve或np.roll累加都能在百格级别瞬间算完100×100网格用纯Python循环虽然也能跑但每步几十毫秒的时间会拖累后面的动画。3. 把一天的模拟画成动画可视化最小闭环3.1 用imshow同步展示感染地图模拟主循环跑起来后最直观的可视化是把二维数组grid直接丢给matplotlib.imshow。状态编码0~3刚好对应四个颜色但一定记得设置vmin0, vmax3否则matplotlib会按当前帧数据的最大最小值自动缩放导致明明都是E和R颜色却在跳动。动画用FuncAnimation更新函数里只需要调用一次set_data不要重新imshow。重新创建图像对象不仅慢还会让动画越跑越卡。基本框架如下import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation def animate_simulation(grid, inc_left, rec_left, total_steps200, prob_infect0.25, latent_days4, recover_days6, seed42): rng np.random.default_rng(seed) fig, (ax_map, ax_curve) plt.subplots(1, 2, figsize(13, 5)) im ax_map.imshow(grid, cmapviridis, vmin0, vmax3) ax_map.set_title(SEIR spatial spread) ax_a ax_curve.plot([], [], labelS, colorblue)[0] ax_b ax_curve.plot([], [], labelE, colororange)[0] ax_c ax_curve.plot([], [], labelI, colorred)[0] ax_d ax_curve.plot([], [], labelR, colorgreen)[0] ax_curve.legend() ax_curve.set_xlim(0, total_steps) ax_curve.set_ylim(0, grid.size) history np.zeros((total_steps, 4)) def update(step): nonlocal grid, inc_left, rec_left if step 0: grid, inc_left, rec_left, counts simulate_step( grid, inc_left, rec_left, prob_infectprob_infect, latent_dayslatent_days, recover_daysrecover_days, rngrng ) else: counts np.bincount(grid.ravel(), minlength4) history[step] counts im.set_data(grid) xdata np.arange(step 1) ax_a.set_data(xdata, history[:step 1, 0]) ax_b.set_data(xdata, history[:step 1, 1]) ax_c.set_data(xdata, history[:step 1, 2]) ax_d.set_data(xdata, history[:step 1, 3]) ax_curve.set_title(fDay {step} I{int(counts[2])}) return im, ax_a, ax_b, ax_c, ax_d anim FuncAnimation(fig, update, framestotal_steps, interval80, blitTrue) plt.show() return anim这里把rng从外部传入而不是每次在simulate_step里新建是保证可复现的关键。interval80表示每帧间隔80毫秒太快会看不清爆发过程太慢又让人着急。我在参数扫描时一般把interval调低到20先用动画确认逻辑再关掉动画跑批量。3.2 实时曲线S/E/I/R人数变化地图只能看到空间分布判断SEIR模型行为是否正确必须同时看人数曲线。曲线用四个plot对象更新时只改set_data避免每一帧clear()重画。S总数随时间下降E先升后降I峰值延后于ER缓慢上升后趋平。如果你发现E的曲线和I的曲线完全同步、没有错峰说明潜伏期没有真正生效。最常见的实现错误是E格子虽然状态是1但倒计时没有走或者E在更新时被当成I参与传染计算。检查办法很简单在simulate_step里统计E-I的人数应该约等于上一帧E总数除以潜伏期。这条曲线的形状也是调节参数的依据。潜伏期变长I的峰值右移恢复期变长I的平台期变宽传染概率提高峰值高度和总感染人数都会变大。批量扫描这些参数就能量化模型敏感性。3.3 参数配置从哪里来潜伏期和恢复期的现实参考很多自己写的模拟项目把参数拍脑袋定成1、2、3结果曲线一天内走完看不出SEIR的“潜伏期滞后”效果。这里给出我常用的参数映射表参数取值范围模拟含义prob_infect0.05 ~ 0.4每个S与单个I相邻一天内的有效传染概率latent_days2 ~ 7E状态持续天数对应现实潜伏期中位数recover_days4 ~ 10I状态持续天数对应症状期到恢复/隔离初始感染者1 ~ 20打一个或多个随机种子点网格尺寸50 ~ 200每个格子代表一个个体总人口网格数注意距离单位的含义。格子之间是邻居关系一天的接触次数被简化成一次判定。如果你把prob_infect设为0.38个邻居全是I时S一天的感染概率是1 - 0.7^8 ≈ 0.94这意味着一个被I包围的S基本必感染。这是合理的因为在真实场景里8个感染邻居意味着极高密度的病毒环境。参数还有一层换算如果你想让“一天”对应真实世界一天潜伏期就写4或5如果你想让一个时间步代表半天潜伏期得乘2。这个换算直接决定动画的时间轴标签建议在代码开头就写明TIME_UNIT day否则改起来容易忘。4. 边界条件与干预机制让模拟更接近现实4.1 周期性边界与固定边界的选择第2章用np.roll实现的周期边界的优点是“世界是封闭的”没有边界墙病毒不会在边缘聚堆。但周期边界也有副作用左上角的居民和右下角的居民是邻居空间拓扑变成一个环面和你平时看的地图不完全一样。如果做的是城市内局部传播模拟通常希望边界是封闭隔离带也就是网格外面无人居住。做法是给数组四周加一圈固定为S的边框内部才是有效区域或者传播统计时跳过边缘。我一般用np.pad固定边界def pad_grid(grid): return np.pad(grid, 1, modeconstant, constant_valuesS)但要注意加了pad之后邻域统计的索引全部要偏移1。更省事的方法是不改数组而是在count_infected_neighbors里禁用边缘格子的感染传播valid_region np.ones(grid.shape, dtypebool) valid_region[0, :] False valid_region[-1, :] False valid_region[:, 0] False valid_region[:, -1] False newly_exposed valid_region不管选哪种边界必须刻意处理。很多初版模拟翻车就是因为索引[-1]被Python解释成倒数第一个元素本意是右边界越界结果右边界和左边界悄悄连起来了画出来一团乱。我的建议是先明确“这个模型世界里人的活动范围有没有出入口”再决定边界方式。4.2 隔离和封锁怎么建模常规SEIR假设所有人都混在一起元胞自动机的好处是可以建模空间干预。隔离病房是最直接的一种当一个I被识别后把它移入隔离区隔离区里的I不再作为传染源。实现方式是在count_infected_neighbors中额外传入一个isolated_maskdef count_infected_neighbors(grid, isolated_maskNone): infected (grid I).astype(np.int8) if isolated_mask is not None: infected[isolated_mask 1] 0 # ...后续累加逻辑不变这样隔离者虽然状态还是I但对周围S没有传染力。封锁则更简单在全剧时间内把prob_infect临时调低到原来的30%并限制人员流动网格里没有移动所以只体现为接触概率降低。注意封锁不是把prob_infect改成0否则会出现疫情被“瞬间清零”的假象现实里封锁也有漏网之鱼。还有一种常见干预是“易感人群接种”。把部分随机S格子直接置为R观察疫情如何绕过免疫人群。这在模拟中几乎零成本但在展示时很有说服力你会在动画里看到病毒在接种区域周围裂成两条传播链。4.3 用总人数守恒验证模型没有BUG在把所有功能堆上去之前先做一次最小验证初始总人数等于每一步四个状态人数之和并且保持不变。一行断言就够了counts np.bincount(new_grid.ravel(), minlength4) assert counts.sum() grid.size, 人数不守恒状态转移里有人被凭空创建或删除我第一次写这类模拟时就把E转I的格子忘在inc_left里清零导致同一个人在第5天被同时当作E和I处理统计人数超过总人口。这种问题靠看曲线很难发现因为S、E、I、R四条线各自看起来都合理但总和在慢慢膨胀。除了总人数还应该验证单调性S只能减少不能增加R只能增加不能减少I从0开始到峰值后回落。如果R出现下降说明有“恢复者再次感染”的代码路径在没有免疫衰退的模型里这是bug。5. 避坑模拟“瞬间全灭”或“指数爆炸”的5个常见问题5.1 概率叠加方式错误导致R0虚高现象把prob_infect调到0.1结果不到20天全网格变红感染曲线陡得像垂直上升。原因更新时直接用new_grid统计邻居感染数且每个S对每个I邻居分别掷骰子多个独立判定让整体感染概率变成了1-(1-p)^nn8时概率是0.57但问题不在这而在于同一轮里刚变E或刚变I的人又立刻参与传播等于一帧内跑了好几代。解决统计邻居必须基于旧grid并且newly_exposed只从旧S中产生。代码上把count_infected_neighbors(grid)放在任何状态转移之前后续所有判断都用这个旧统计结果。5.2 没有固定随机种子曲线无法复现现象每次运行代码峰值时间、总感染人数都不一样调参时根本分不清是参数影响还是随机波动。原因np.random.random每次使用系统熵没有种子。解决在一开始创建rng np.random.default_rng(42)并一路传给simulate_step。跑多次分析时只对“初始感染者坐标”做随机化其他随机源全部固定。注意不要把np.random.seed和default_rng混用np.random.seed只影响旧接口传参更可控。5.3 帧列表越攒越多保存GIF时内存爆炸现象动画在屏幕上跑没问题但按下保存后内存飙升笔记本风扇狂转几G内存瞬间耗尽。原因用ani.save()保存GIF时matplotlib 会先把每一帧渲染成一幅图像并放到列表里特别是progress_callback的实现200帧以上就非常吃内存。如果是blitTrue的动画保存时还要重新渲染。解决优先保存为 MP4用 FFMpegWriter 流式写入不内存累积如果必须GIF控制total_steps在100以内或者减少interval但用pillow手动按帧写from PIL import Image frames [] # 每个step已经有无标题的fig图保存成RGB数组 rgb np.array(fig.canvas.buffer_rgba())[:, :, :3] frames.append(Image.fromarray(rgb)) frames[0].save(seir.gif, save_allTrue, append_imagesframes[1:], duration80, loop0)手动写法虽然代码多一点但每一帧是即时转换后丢弃不会把所有帧原图留在内存里。5.4 用list逐元素更新100×100网格要跑几十秒现象每个格子用for循环判断邻居logical 逻辑本身没问题但一步耗时几十毫秒200步就是好几秒动画基本卡死。原因Python的嵌套循环在二维网格上是灾难100×100是1万个格子每格访问8个邻居就是8万次操作。加上状态分类和倒计时轻松超过百万次Python级操作。解决能用numpy向量化就用向量化grid E这种布尔数组操作是C级循环再配合np.roll统计邻居性能和for循环是天壤之别。如果网格超过300×300还能用scipy.ndimage.convolve进一步加速from scipy.ndimage import convolve kern np.ones((3, 3), dtypenp.int8) kern[1, 1] 0 infected_counts convolve((grid I).astype(np.int8), kern, modewrap)注意scipy的modewrap就是周期边界。如果不想引入scipy依赖np.roll方案足够用。5.5 潜伏期倒计时和恢复倒计时没同步推进现象E人数曲线在变成I的当天出现“异常的凸起”或者I人数比E人数先达到峰值。原因状态和倒计时不是同一套逻辑。比如先new_grid[to_infectious] I后后面的倒计时递减条件写成new_grid E结果刚变成I的格子不再递减new_inc但它已经用完了潜伏期而新感染E的格子又被重复扣了一天潜伏期比设定少一天。解决倒计时的递减对象必须是旧状态grid不因本轮转移而改变。我在第2章的代码里特意用了(grid E) ~to_infectious而不是(new_grid E)就是为了避免这个坑。建议在simulation_step结束前加打印调试每个类别的到期人数和剩余人数都应该是单调的。6. 批处理与导出不只有动画还能看参数敏感性6.1 把模拟包成一个可复用函数动画适合人眼观察但不适合“量化对比”。我会把模拟过程抽成独立函数只返回每一时刻的S/E/I/R计数完全不绘制图像。这样批量扫描参数时就不用打开窗口一张张看。def run_simulation(grid_size100, initial_infected5, prob_infect0.25, latent_days4, recover_days6, steps200, seed1): rng np.random.default_rng(seed) grid np.full((grid_size, grid_size), S, dtypenp.int8) inc_left np.zeros_like(grid) rec_left np.zeros_like(grid) # 随机洒初始感染者避免集中在中心 idx rng.choice(grid_size * grid_size, sizeinitial_infected, replaceFalse) rows, cols np.unravel_index(idx, grid.shape) grid[rows, cols] I rec_left[rows, cols] recover_days history np.zeros((steps, 4)) for step in range(steps): grid, inc_left, rec_left, counts simulate_step( grid, inc_left, rec_left, prob_infectprob_infect, latent_dayslatent_days, recover_daysrecover_days, rngrng, ) history[step] counts return history这个函数的好处是参数全在入口随机种子可控。需要检查某次运行的具体状态时在循环里加一个if step断点即可不用改动画代码。6.2 批量扫描感染概率画对比曲线用循环调用run_simulation把不同prob_infect的I人数曲线画在同一张图上import matplotlib.pyplot as plt for p in [0.1, 0.2, 0.3, 0.4]: hist run_simulation(prob_infectp, seed42) plt.plot(hist[:, 2], labelfp{p}) plt.xlabel(time step) plt.ylabel(Infectious count) plt.legend() plt.title(Sensitivity of infection probability) plt.show()你很快会发现几个不变规律p越大峰值越高、峰值越早p低于某个阈值时疫情自然消亡。这个“阈值”就是模拟里的流行病学基本再生数R0临界。用扫描结果对比比在动画里猜参数靠谱得多。同样的方法也可以扫latent_days和recover_days。我习惯固定seed一次扫一个维度其他变量都不动这样每条曲线之间的差异完全来自目标参数。如果seed不固定扫描结果会被随机噪声淹没得不出清晰结论。6.3 保存GIF动画时的一个小技巧如果你最终还是要一份GIF放进课件有个细节容易忽略FuncAnimation.save的时长由interval控制和模拟的steps无关。比如你模拟了500步但GIF默认100ms一帧出来就是50秒的慢放。建议保存前检查durationanim.save(seir.gif, writerpillow, fps15)fps和interval的换算关系是fps 1000 / interval。如果你动画是80ms一帧保存时fps12就接近原速。GIF文件大小也和颜色数量有关SEIR只有4个状态可以用optimizeTrue减少体积一般能压到1~2MB。最后说一句我的习惯写这种模拟项目时永远先跑一个不带可视化的批处理脚本只有确认曲线形态正确后才打开动画看空间传播。否则动画里一粒粒红点扩散很酷但你可能花一下午在欣赏一个根本不符合常理的疫情。空间可视化是验证模型的最后一步不是第一优先。希望这套代码和排错思路对你有所帮助也祝你调出一张既好看又可信的传播图。本文还有配套的精品资源点击获取 元胞自动机模拟SEIR疫情：Python实现与可视化
简介一套基于Python的元胞自动机病毒传染SEIR模型可视化实现面向Python学习者、数据科学初学者及对流行病模拟感兴趣的研究者。项目借助Numpy构建二维元胞矩阵将人群划分为易感者、暴露者、感染者、康复者四类状态通过Matplotlib动态展示病毒传播过程可直观理解局部交互与空间因素对疫情发展的影响。资源包共3个文件约8KB包含一个ZIP压缩包、一个核心Python脚本和一份Markdown说明文档。ZIP内应封装完整项目或依赖内容py文件实现SEIR状态转移与元胞自动机更新规则md文档则提供代码结构说明与使用指引便于快速运行和二次开发。目前已有142人学习下载。实际使用中读者可调整感染率、恢复率、潜伏期等参数观察不同干预策略下感染者曲线与空间扩散热力图的变化从而掌握元胞自动机建模、SEIR模型数值模拟以及Python科学计算可视化的完整流程对课程设计、入门科研或公共卫生策略演示均有实用价值。1. 为什么用元胞自动机跑SEIR空间传播是微分方程看不见的先说明一点SEIR模型原本是一组常微分方程S、E、I、R四类人的人数随时间变化画出来就是几条平滑曲线。但曲线不会告诉你“病毒从哪个街区开始传为什么有些小区一直没事有些地方一夜之间全变红”。如果你要的是直观演示、教学课件或者疫情传播空间特征研究元胞自动机是比微分方程更合适的载体。这个标题里的rar工程做的就是把SEIR状态机放进一个个格子组成的世界里用Python模拟每天谁被传染、谁进入潜伏期、谁恢复再把这一过程用动态图展示出来。我在拿到类似项目时第一步不是看可视化而是先把格子世界的更新规则建立在脑子里。因为模拟类项目最容易出的问题不是代码跑不起来而是“跑起来了但结果一眼假”。这篇文章按“状态机→模拟主循环→可视化→干预机制→踩坑→批处理参数扫描”的顺序展开你照着写一份自己的模拟400行以内就能交付。2. SEIR状态机与元胞更新规则先把格子世界搭起来2.1 四个状态与状态转移条件元胞自动机的核心是“每个格子的明天由它自己和邻居的今天决定”。在SEIR模型里每个格子代表一个人状态只有四种S易感、E暴露/潜伏、I感染、R恢复。这个离散化人和空间的做法和真实感染过程非常像——你不会因为“全市平均感染率”被传染而是因为你身边有人感染了。状态转移有三条规则S遇到邻居I有一定概率变成E。这一步是概率事件不是必然事件。E经过潜伏期后变成I。潜伏期是固定天数比如4天到期自动转感染。I经过恢复期后变成R。恢复期结束自动转恢复并假设获得免疫不再感染。关键点是E到I、I到R是确定性转移S到E是随机转移。如果你把S到E也做成确定规则病毒会像潮水一样扫过整个网格真实疫情反倒没这么整齐。随机性来自两个地方一是邻居中有没有I二是每次接触是否真的“中招”。前者由空间结构决定后者由传染概率控制。2.2 最小可运行的时间步函数下面这段代码是最简实现核心是用两个额外数组记录E和I的剩余天数避免把状态和时间耦合在同一个小格子值里。状态数组只存0、1、2、3倒计时数组另存。import numpy as np # 状态编码 S, E, I, R 0, 1, 2, 3 def count_infected_neighbors(grid): 统计每个格子的8邻域中感染者的数量周期性边界由 np.roll 处理 infected (grid I).astype(np.int8) counts np.zeros_like(infected) for di in (-1, 0, 1): for dj in (-1, 0, 1): if di 0 and dj 0: continue counts np.roll(np.roll(infected, di, axis0), dj, axis1) return counts def simulate_step(grid, inc_left, rec_left, prob_infect0.25, latent_days4, recover_days6, rngNone): new_grid grid.copy() new_inc inc_left.copy() new_rec rec_left.copy() # E 潜伏期结束转为 I to_infectious (grid E) (inc_left 1) new_grid[to_infectious] I new_inc[to_infectious] 0 new_rec[to_infectious] recover_days # I 恢复期结束转为 R to_recovered (grid I) (rec_left 1) new_grid[to_recovered] R new_inc[to_recovered] 0 new_rec[to_recovered] 0 # S 被邻居中的 I 感染 infected_counts count_infected_neighbors(grid) susceptible (grid S) # 每个 I 独立接触时合并为一次骰子见 2.3 说明 exposure_prob 1 - (1 - prob_infect) ** infected_counts if rng is None: rng np.random.default_rng() newly_exposed susceptible (infected_counts 0) (rng.random(grid.shape) exposure_prob) new_grid[newly_exposed] E new_inc[newly_exposed] latent_days # 未发生状态转移的 E/I剩余天数减 1 new_inc[(grid E) ~to_infectious] - 1 new_rec[(grid I) ~to_recovered] - 1 # 统计 S E I R 各有多少人 counts np.bincount(new_grid.ravel(), minlength4) return new_grid, new_inc, new_rec, counts代码里最容易被忽略的是同步更新传染概率计算用的grid是旧状态而不是这一步刚刚变成I的人。如果你用更新后的new_grid去统计感染者会造成同一步内“新感染者立刻又传染别人”一个时间步里病毒跳好几代R0 被严重高估。count_infected_neighbors里用np.roll做的周期边界意味着网格最上面一行和最下面一行是邻居。对于传播模拟来说这避免了“病毒撞到边界就停住”的人为假象。如果你希望边界是隔离区比如一座岛的边缘可以把np.roll换成np.pad固定边界后面第4章会细说。2.3 邻居统计与感染概率的两种写法常见做法有两种区别直接影响疫情走势。写法一每个S格子遍历8个邻居凡遇到I就独立掷一次骰子阈值是prob_infect。当I邻居多时相当于一天内多次接触感染概率是1 - (1 - prob_infect)^n也就是代码里的合并写法。但很多初版代码会写成循环里“只要有任何一个I击中就感染”那么同一个S在本轮可能被连续判定多次结果概率被重复计算曲线会偏爆炸。写法二只判断“邻居是否有I”有I就以prob_infect概率感染。这种方式简单但把多个I邻居和一个I邻居同等对待低估了高密度感染区的传播速度。我建议用合并概率1 - (1 - p)^n它介于上面两种极端之间参数含义也更接近流行病学里的“单位时间内每次有效接触的传播概率”。如果你想跑得快可以用卷积替代循环。scipy.ndimage.convolve或np.roll累加都能在百格级别瞬间算完100×100网格用纯Python循环虽然也能跑但每步几十毫秒的时间会拖累后面的动画。3. 把一天的模拟画成动画可视化最小闭环3.1 用imshow同步展示感染地图模拟主循环跑起来后最直观的可视化是把二维数组grid直接丢给matplotlib.imshow。状态编码0~3刚好对应四个颜色但一定记得设置vmin0, vmax3否则matplotlib会按当前帧数据的最大最小值自动缩放导致明明都是E和R颜色却在跳动。动画用FuncAnimation更新函数里只需要调用一次set_data不要重新imshow。重新创建图像对象不仅慢还会让动画越跑越卡。基本框架如下import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation def animate_simulation(grid, inc_left, rec_left, total_steps200, prob_infect0.25, latent_days4, recover_days6, seed42): rng np.random.default_rng(seed) fig, (ax_map, ax_curve) plt.subplots(1, 2, figsize(13, 5)) im ax_map.imshow(grid, cmapviridis, vmin0, vmax3) ax_map.set_title(SEIR spatial spread) ax_a ax_curve.plot([], [], labelS, colorblue)[0] ax_b ax_curve.plot([], [], labelE, colororange)[0] ax_c ax_curve.plot([], [], labelI, colorred)[0] ax_d ax_curve.plot([], [], labelR, colorgreen)[0] ax_curve.legend() ax_curve.set_xlim(0, total_steps) ax_curve.set_ylim(0, grid.size) history np.zeros((total_steps, 4)) def update(step): nonlocal grid, inc_left, rec_left if step 0: grid, inc_left, rec_left, counts simulate_step( grid, inc_left, rec_left, prob_infectprob_infect, latent_dayslatent_days, recover_daysrecover_days, rngrng ) else: counts np.bincount(grid.ravel(), minlength4) history[step] counts im.set_data(grid) xdata np.arange(step 1) ax_a.set_data(xdata, history[:step 1, 0]) ax_b.set_data(xdata, history[:step 1, 1]) ax_c.set_data(xdata, history[:step 1, 2]) ax_d.set_data(xdata, history[:step 1, 3]) ax_curve.set_title(fDay {step} I{int(counts[2])}) return im, ax_a, ax_b, ax_c, ax_d anim FuncAnimation(fig, update, framestotal_steps, interval80, blitTrue) plt.show() return anim这里把rng从外部传入而不是每次在simulate_step里新建是保证可复现的关键。interval80表示每帧间隔80毫秒太快会看不清爆发过程太慢又让人着急。我在参数扫描时一般把interval调低到20先用动画确认逻辑再关掉动画跑批量。3.2 实时曲线S/E/I/R人数变化地图只能看到空间分布判断SEIR模型行为是否正确必须同时看人数曲线。曲线用四个plot对象更新时只改set_data避免每一帧clear()重画。S总数随时间下降E先升后降I峰值延后于ER缓慢上升后趋平。如果你发现E的曲线和I的曲线完全同步、没有错峰说明潜伏期没有真正生效。最常见的实现错误是E格子虽然状态是1但倒计时没有走或者E在更新时被当成I参与传染计算。检查办法很简单在simulate_step里统计E-I的人数应该约等于上一帧E总数除以潜伏期。这条曲线的形状也是调节参数的依据。潜伏期变长I的峰值右移恢复期变长I的平台期变宽传染概率提高峰值高度和总感染人数都会变大。批量扫描这些参数就能量化模型敏感性。3.3 参数配置从哪里来潜伏期和恢复期的现实参考很多自己写的模拟项目把参数拍脑袋定成1、2、3结果曲线一天内走完看不出SEIR的“潜伏期滞后”效果。这里给出我常用的参数映射表参数取值范围模拟含义prob_infect0.05 ~ 0.4每个S与单个I相邻一天内的有效传染概率latent_days2 ~ 7E状态持续天数对应现实潜伏期中位数recover_days4 ~ 10I状态持续天数对应症状期到恢复/隔离初始感染者1 ~ 20打一个或多个随机种子点网格尺寸50 ~ 200每个格子代表一个个体总人口网格数注意距离单位的含义。格子之间是邻居关系一天的接触次数被简化成一次判定。如果你把prob_infect设为0.38个邻居全是I时S一天的感染概率是1 - 0.7^8 ≈ 0.94这意味着一个被I包围的S基本必感染。这是合理的因为在真实场景里8个感染邻居意味着极高密度的病毒环境。参数还有一层换算如果你想让“一天”对应真实世界一天潜伏期就写4或5如果你想让一个时间步代表半天潜伏期得乘2。这个换算直接决定动画的时间轴标签建议在代码开头就写明TIME_UNIT day否则改起来容易忘。4. 边界条件与干预机制让模拟更接近现实4.1 周期性边界与固定边界的选择第2章用np.roll实现的周期边界的优点是“世界是封闭的”没有边界墙病毒不会在边缘聚堆。但周期边界也有副作用左上角的居民和右下角的居民是邻居空间拓扑变成一个环面和你平时看的地图不完全一样。如果做的是城市内局部传播模拟通常希望边界是封闭隔离带也就是网格外面无人居住。做法是给数组四周加一圈固定为S的边框内部才是有效区域或者传播统计时跳过边缘。我一般用np.pad固定边界def pad_grid(grid): return np.pad(grid, 1, modeconstant, constant_valuesS)但要注意加了pad之后邻域统计的索引全部要偏移1。更省事的方法是不改数组而是在count_infected_neighbors里禁用边缘格子的感染传播valid_region np.ones(grid.shape, dtypebool) valid_region[0, :] False valid_region[-1, :] False valid_region[:, 0] False valid_region[:, -1] False newly_exposed valid_region不管选哪种边界必须刻意处理。很多初版模拟翻车就是因为索引[-1]被Python解释成倒数第一个元素本意是右边界越界结果右边界和左边界悄悄连起来了画出来一团乱。我的建议是先明确“这个模型世界里人的活动范围有没有出入口”再决定边界方式。4.2 隔离和封锁怎么建模常规SEIR假设所有人都混在一起元胞自动机的好处是可以建模空间干预。隔离病房是最直接的一种当一个I被识别后把它移入隔离区隔离区里的I不再作为传染源。实现方式是在count_infected_neighbors中额外传入一个isolated_maskdef count_infected_neighbors(grid, isolated_maskNone): infected (grid I).astype(np.int8) if isolated_mask is not None: infected[isolated_mask 1] 0 # ...后续累加逻辑不变这样隔离者虽然状态还是I但对周围S没有传染力。封锁则更简单在全剧时间内把prob_infect临时调低到原来的30%并限制人员流动网格里没有移动所以只体现为接触概率降低。注意封锁不是把prob_infect改成0否则会出现疫情被“瞬间清零”的假象现实里封锁也有漏网之鱼。还有一种常见干预是“易感人群接种”。把部分随机S格子直接置为R观察疫情如何绕过免疫人群。这在模拟中几乎零成本但在展示时很有说服力你会在动画里看到病毒在接种区域周围裂成两条传播链。4.3 用总人数守恒验证模型没有BUG在把所有功能堆上去之前先做一次最小验证初始总人数等于每一步四个状态人数之和并且保持不变。一行断言就够了counts np.bincount(new_grid.ravel(), minlength4) assert counts.sum() grid.size, 人数不守恒状态转移里有人被凭空创建或删除我第一次写这类模拟时就把E转I的格子忘在inc_left里清零导致同一个人在第5天被同时当作E和I处理统计人数超过总人口。这种问题靠看曲线很难发现因为S、E、I、R四条线各自看起来都合理但总和在慢慢膨胀。除了总人数还应该验证单调性S只能减少不能增加R只能增加不能减少I从0开始到峰值后回落。如果R出现下降说明有“恢复者再次感染”的代码路径在没有免疫衰退的模型里这是bug。5. 避坑模拟“瞬间全灭”或“指数爆炸”的5个常见问题5.1 概率叠加方式错误导致R0虚高现象把prob_infect调到0.1结果不到20天全网格变红感染曲线陡得像垂直上升。原因更新时直接用new_grid统计邻居感染数且每个S对每个I邻居分别掷骰子多个独立判定让整体感染概率变成了1-(1-p)^nn8时概率是0.57但问题不在这而在于同一轮里刚变E或刚变I的人又立刻参与传播等于一帧内跑了好几代。解决统计邻居必须基于旧grid并且newly_exposed只从旧S中产生。代码上把count_infected_neighbors(grid)放在任何状态转移之前后续所有判断都用这个旧统计结果。5.2 没有固定随机种子曲线无法复现现象每次运行代码峰值时间、总感染人数都不一样调参时根本分不清是参数影响还是随机波动。原因np.random.random每次使用系统熵没有种子。解决在一开始创建rng np.random.default_rng(42)并一路传给simulate_step。跑多次分析时只对“初始感染者坐标”做随机化其他随机源全部固定。注意不要把np.random.seed和default_rng混用np.random.seed只影响旧接口传参更可控。5.3 帧列表越攒越多保存GIF时内存爆炸现象动画在屏幕上跑没问题但按下保存后内存飙升笔记本风扇狂转几G内存瞬间耗尽。原因用ani.save()保存GIF时matplotlib 会先把每一帧渲染成一幅图像并放到列表里特别是progress_callback的实现200帧以上就非常吃内存。如果是blitTrue的动画保存时还要重新渲染。解决优先保存为 MP4用 FFMpegWriter 流式写入不内存累积如果必须GIF控制total_steps在100以内或者减少interval但用pillow手动按帧写from PIL import Image frames [] # 每个step已经有无标题的fig图保存成RGB数组 rgb np.array(fig.canvas.buffer_rgba())[:, :, :3] frames.append(Image.fromarray(rgb)) frames[0].save(seir.gif, save_allTrue, append_imagesframes[1:], duration80, loop0)手动写法虽然代码多一点但每一帧是即时转换后丢弃不会把所有帧原图留在内存里。5.4 用list逐元素更新100×100网格要跑几十秒现象每个格子用for循环判断邻居logical 逻辑本身没问题但一步耗时几十毫秒200步就是好几秒动画基本卡死。原因Python的嵌套循环在二维网格上是灾难100×100是1万个格子每格访问8个邻居就是8万次操作。加上状态分类和倒计时轻松超过百万次Python级操作。解决能用numpy向量化就用向量化grid E这种布尔数组操作是C级循环再配合np.roll统计邻居性能和for循环是天壤之别。如果网格超过300×300还能用scipy.ndimage.convolve进一步加速from scipy.ndimage import convolve kern np.ones((3, 3), dtypenp.int8) kern[1, 1] 0 infected_counts convolve((grid I).astype(np.int8), kern, modewrap)注意scipy的modewrap就是周期边界。如果不想引入scipy依赖np.roll方案足够用。5.5 潜伏期倒计时和恢复倒计时没同步推进现象E人数曲线在变成I的当天出现“异常的凸起”或者I人数比E人数先达到峰值。原因状态和倒计时不是同一套逻辑。比如先new_grid[to_infectious] I后后面的倒计时递减条件写成new_grid E结果刚变成I的格子不再递减new_inc但它已经用完了潜伏期而新感染E的格子又被重复扣了一天潜伏期比设定少一天。解决倒计时的递减对象必须是旧状态grid不因本轮转移而改变。我在第2章的代码里特意用了(grid E) ~to_infectious而不是(new_grid E)就是为了避免这个坑。建议在simulation_step结束前加打印调试每个类别的到期人数和剩余人数都应该是单调的。6. 批处理与导出不只有动画还能看参数敏感性6.1 把模拟包成一个可复用函数动画适合人眼观察但不适合“量化对比”。我会把模拟过程抽成独立函数只返回每一时刻的S/E/I/R计数完全不绘制图像。这样批量扫描参数时就不用打开窗口一张张看。def run_simulation(grid_size100, initial_infected5, prob_infect0.25, latent_days4, recover_days6, steps200, seed1): rng np.random.default_rng(seed) grid np.full((grid_size, grid_size), S, dtypenp.int8) inc_left np.zeros_like(grid) rec_left np.zeros_like(grid) # 随机洒初始感染者避免集中在中心 idx rng.choice(grid_size * grid_size, sizeinitial_infected, replaceFalse) rows, cols np.unravel_index(idx, grid.shape) grid[rows, cols] I rec_left[rows, cols] recover_days history np.zeros((steps, 4)) for step in range(steps): grid, inc_left, rec_left, counts simulate_step( grid, inc_left, rec_left, prob_infectprob_infect, latent_dayslatent_days, recover_daysrecover_days, rngrng, ) history[step] counts return history这个函数的好处是参数全在入口随机种子可控。需要检查某次运行的具体状态时在循环里加一个if step断点即可不用改动画代码。6.2 批量扫描感染概率画对比曲线用循环调用run_simulation把不同prob_infect的I人数曲线画在同一张图上import matplotlib.pyplot as plt for p in [0.1, 0.2, 0.3, 0.4]: hist run_simulation(prob_infectp, seed42) plt.plot(hist[:, 2], labelfp{p}) plt.xlabel(time step) plt.ylabel(Infectious count) plt.legend() plt.title(Sensitivity of infection probability) plt.show()你很快会发现几个不变规律p越大峰值越高、峰值越早p低于某个阈值时疫情自然消亡。这个“阈值”就是模拟里的流行病学基本再生数R0临界。用扫描结果对比比在动画里猜参数靠谱得多。同样的方法也可以扫latent_days和recover_days。我习惯固定seed一次扫一个维度其他变量都不动这样每条曲线之间的差异完全来自目标参数。如果seed不固定扫描结果会被随机噪声淹没得不出清晰结论。6.3 保存GIF动画时的一个小技巧如果你最终还是要一份GIF放进课件有个细节容易忽略FuncAnimation.save的时长由interval控制和模拟的steps无关。比如你模拟了500步但GIF默认100ms一帧出来就是50秒的慢放。建议保存前检查durationanim.save(seir.gif, writerpillow, fps15)fps和interval的换算关系是fps 1000 / interval。如果你动画是80ms一帧保存时fps12就接近原速。GIF文件大小也和颜色数量有关SEIR只有4个状态可以用optimizeTrue减少体积一般能压到1~2MB。最后说一句我的习惯写这种模拟项目时永远先跑一个不带可视化的批处理脚本只有确认曲线形态正确后才打开动画看空间传播。否则动画里一粒粒红点扩散很酷但你可能花一下午在欣赏一个根本不符合常理的疫情。空间可视化是验证模型的最后一步不是第一优先。希望这套代码和排错思路对你有所帮助也祝你调出一张既好看又可信的传播图。本文还有配套的精品资源点击获取