简介DESSABNeT 是一个基于 MATLAB 编写的开源传染病传播模拟软件包采用基于代理的模型ABM构建人工社会以真实人口统计数据为输入可模拟任何城市或区域中传染病的传播动态最初专门用于澳大利亚城市的 COVID-19 疫情建模。面向流行病学研究者、公共卫生决策者及计算社科领域的 MATLAB 用户可用于评估社交限制措施、模拟不同疾病爆发场景。资源包共 42 个文件以 33 个 .m 格式 MATLAB 脚本为核心覆盖人口设置、社交网络构建、每周时间表配置、疾病传播动态模块及社交限制阶段模拟另有 7 个 .mat 预生成数据文件、1 个 Excel 数据表和 1 个 README 说明文档压缩包整体 2.7MB按模块分目录组织。已有 712 人学习下载。读者可获得完整的 ABM 建模框架代码包括悉尼人口结构设置、家庭/工作/学校网络构建、感染病例注入与传播计算、R0 矩阵计算、隔离与社交限制模拟等并提供可直接运行的示例数据便于移植到其他传染病或城市场景开展仿真实验。1. 这份MATLAB模型为什么值得你先跑一遍再说做传染病仿真的人多数有一个共识SEIR已经不够用了。真实世界每一轮疫情都会出现死亡病例而标准SEIR只把人群分成易感、暴露、感染、恢复四类死亡和康复只占一个R槽模拟出来反而比实际结果乐观不少。这套名为DESSABNeT的SEIRD模型MATLAB代码正好补上了这段缺口——它把死亡D单独拆出来同时把消毒网Disinfection Network的影响以参数形式接入传播链路跑出来的曲线更贴合实际决策场景。适合谁用一类是做疾控策略模拟的研究生和工程师另一类是想在MATLAB里快速做多情形对比的仿真开发者。不需要额外装任何工具包打开主脚本就能看到感染峰值从哪一天开始放缓、消毒强度变化会带来多大差别。2. 从SEIR到SEIRD模型演进与技术栈选型2.1 标准SEIR丢掉的死亡D在这里单独成仓先厘清一个概念SEIRD模型里五个字母分别代表Susceptible易感、Exposed暴露、Infectious感染、Recovered恢复、Deceased死亡。传统的SEIR把R视为“移除”康复和死亡混在一起这会带来两个问题第一死亡病例无法参与后续决策策略评估尤其是在计算床位占用率、病亡率时完全失真第二当恢复率高而死亡少时还好一旦疫情严重R的累加值会掩盖真实致死趋势。DESSABNeT在这套实现里把D配成独立仓室并用“消毒网”参数把外部干预动作嵌进微分方程组。简单说消毒不再只是模型外部的一个注释而是直接作用在传染系数β上的削减因子。细看代码能发现传染率的大小除了受基础传播系数影响还乘了一个随时段改变的消毒强度序列这样每一次消毒动作都会被模型记录成传染曲线上的一个拐点。2.2 技术栈与仓库结构两个文件夹分工十分明确DESSABNeT的MATLAB工程结构做得很轻没有复杂的包管理器依赖核心就两个目录sim/和plot/。前者放模型主体函数后者放可视化脚本职责清楚不会出现“改了绘图却要连带翻模型代码”的情况。文件夹大致是这样DESSABNeT/ ├── sim/ │ ├── main_seird.m # 主运行脚本定义参数 │ ├── seird_ode.m # SEIRD微分方程核心函数 │ ├── disinfection.m # 消毒策略生成函数 │ └── contact_matrix.m # 人群接触矩阵构建可选 ├── plot/ │ ├── plot_curves.m # 绘制五仓室曲线 │ ├── plot_peak.m # 感染峰值对比绘图 │ └── plot_heatmap.m # 消毒强度-峰值热力图 └── data/ └── baseline_params.mat # 默认参数文件main_seird.m是入口所有参数集中在文件头部适合做批量仿真seird_ode.m是模型本体没有任何界面代码方便二次开发时抽出来单测。这套结构对于研究者来说很友好因为它遵循的是一条黄金法则参数、方程、可视化三者分离。如果某天要把模型迁移到Python或R只需要重写seird_ode.m其他部分都可以原样保留。3. 核心模块拆解SEIRD微分主循环与DESSABNeT消毒策略接入3.1 主模型函数五仓室ODE的底层逻辑打开seird_ode.m核心代码并不长但每一行的物理含义都需要说清楚。这个文件是整套模型的运算心脏MATLAB的ode45会反复调用它来计算每一步的仓室变化量。function dydt seird_ode(t, y, beta0, sigma, gamma, delta, dis_amp) % y 列向量: [S; E; I; R; D] S y(1); E y(2); I y(3); R y(4); D y(5); N S E I R D; % 总人口保持恒定 % 消毒网作用传染率随时间被削减 beta beta0 * (1 - dis_amp * exp(-0.1 * t)); dS -beta * S * I / N; dE beta * S * I / N - sigma * E; dI sigma * E - (gamma delta) * I; dR gamma * I; dD delta * I; dydt [dS; dE; dI; dR; dD]; end这段代码最值得细看的是beta那一行。beta0是不做任何干预时的基础传染率dis_amp是消毒强度系数exp(-0.1*t)表示消毒效果随时间的衰减。这种设计隐含了一个公共卫生常识消毒的边际效应是逐日递减的不可能连续三十天保持同样的消杀效率。如果直接把消毒设为常量模拟后期会出现曲线异常偏低、与真实疫情趋势脱节的问题。参数sigma代表潜伏期结束进入感染期的速率gamma是恢复速率delta是因病死亡速率。需要特别留意的是gamma delta共同决定感染者仓室的迁出速度两者此消彼长会直接影响感染峰值的形状。如果调高delta而不动gammaI仓室会更快清空但D仓室会累积出一个令人警醒的死亡曲线。3.2 消毒策略接入disinfection函数怎么生成干预计划disinfection.m的作用是生成一段可读的消毒强度时间表再把这张表喂给主模型。单独拆这个函数出来是有好处的因为实际场景中消毒策略不是连续函数而是分段操作——今天全城消杀明天可能只有重点区域消毒。function dis_table disinfection(base_strength, schedule, tspan) % schedule: [start_day, end_day, strength_scale] % 返回与 tspan 等长的消毒强度向量 dis_table zeros(length(tspan), 1); for i 1:size(schedule, 1) start_day schedule(i, 1); end_day schedule(i, 2); scale schedule(i, 3); idx tspan start_day tspan end_day; dis_table(idx) base_strength * scale; end end这个函数的输入schedule是一个矩阵每一行代表一次消毒行动。比如[10, 20, 1.5]表示第10天到第20天之间执行1.5倍强度的消杀作业。输出向量与仿真时间轴等长方便直接点乘到传染率上。我一般会把schedule从Excel里读进来这样非编程背景的同事也能直接修改干预计划不用动代码。二次开发时也可以把一个真实的消毒日志转成这个格式——只要日期和强度对齐就行。contact_matrix.m是另一个可选模块。当总人口N被分成多个子群时比如学生、在职人员、老年人每个子群的接触频率不同用单一β值就不够了。这时传入一个接触矩阵C传染率会变成beta0 * C * S * I / N更适合做异质性人群的传播模拟。默认工程只跑单一人群但结构上已经预留了这个扩展点。3.3 主入口main_seird串联全流程的关键脚本main_seird.m是直接运行的文件它的任务是把参数定义、ODE求解、结果返回整合成一条流水线。很多人在这个文件上翻车因为改完参数立刻运行结果还在看上一次的原因是MATLAB的工作区缓存没有清理。clear; clc; % 基础参数 N 100000; % 总人口 I0 10; % 初始感染者 E0 0; % 初始暴露者 sigma 1/5.2; % 潜伏期约5.2天 gamma 1/10; % 恢复期约10天 delta 0.015; % 死亡率 beta0 0.32; % 基础传染率 % 消毒策略 schedule [5, 15, 1.2; 20, 30, 0.8]; dis_amp 0.25; % 时间跨度 tspan 0:0.5:100; y0 [N - I0; E0; I0; 0; 0]; % 求解 [t, y] ode45((t, y) seird_ode(t, y, beta0, sigma, gamma, delta, dis_amp), ... tspan, y0); % 留存结果 result.t t; result.y y; save(data/run_result.mat, result);注意这里tspan用步长0.5而不是默认的整数步长原因在于ode45是自适应步长算法输出节点不一定均匀如果后面要做傅里叶分析或频谱分析不均匀时间轴会带来很多额外麻烦。统一成0.5天一步既保持曲线光滑又方便和真实日报告数据做对齐。参数I0设置为10意味着第一天就有10个感染者进场。实际使用中应该把I0设成1或者更小因为社区传播的起点往往很低设置过大会让峰值前置、显得疫情爆发太快。4. 从克隆到跑通运行配置与三大参数组的调优方法4.1 环境准备与两种运行方式MATLAB环境基本没有额外要求只要是R2016b之后的版本都能直接跑。仓库不依赖Simulink也不依赖任何第三方工具箱普通的MATLAB基础授权就够。运行方式有两种。第一种是图形界面操作打开main_seird.m在编辑器中直接按F5运行结果保存在data/run_result.mat然后运行plot/plot_curves.m绘图。第二种是命令行批处理适合批量跑参数扫描matlab -batch run(main_seird.m); plot_curves用-batch模式可以避开MATLAB启动GUI的额外开销连续跑几十组参数时体感快很多。如果机器上装了并行计算工具箱还能用parfor把多组参数同时送进求解器。需要注意的一点是-batch模式下工作区不保留变量必须在脚本内部显式用save把结果写到磁盘否则跑完什么都没留下。4.2 三大参数组传播强度、病程时长、消毒强度这组参数可以划分为三组理解每组对应模型行为的一个方面表格核心参数调节方向参数组参数名含义调高后的影响传播强度beta0基础传染率峰值更高、峰值时间提前病程时长sigma,gamma潜伏期/恢复期时长病程越长感染平台期越久消毒干预dis_amp,schedule消毒强度与时间表峰值下降曲线出现明显转折beta0的取值可以参考基本再生数R0来推算近似公式是R0 beta0 / (gamma delta)。比如希望模拟R03.5的传染病假设gamma0.1,delta0.015那么beta0应该设在3.5 * 0.115 ≈ 0.4附近。这里特别提醒算出来的值只是一个参考起点实际拟合时还需要根据疫情数据反向校准。sigma写的是1/5.2意思是潜伏期平均5.2天这对应指数分布的速率参数。如果改成1/3就表示潜伏期缩短到3天感染峰值会明显左移。调整时注意不要把恢复期和潜伏期设成同一数量级否则E和I仓室会同时膨胀曲线之间互相遮掩后续可视化很难解读。4.3 验证模型行为曲线不收敛时怎么看模型跑通之后先别急着做策略对比第一步是检查五条曲线是否符合常识。正常情况下S曲线缓慢下降E和I先后出现一个单峰R曲线单调上升D曲线缓慢累积。如果看到I曲线剧烈震荡或者出现负值说明数值求解出了问题。绝大多数震荡来自ode45对刚性方程的不适应。SEIRD模型在参数极端时比如beta0很大、delta接近0会呈现刚性特征此时换成ode15s求解器往往能立刻解决[t, y] ode15s((t, y) seird_ode(t, y, beta0, sigma, gamma, delta, dis_amp), ... tspan, y0);另外一个常见验证手段是计算总人口守恒。每行SEIRD理论上恒等于初始总人口N如果数值上出现漂移大概率是因为某个仓室的更新率在积分中丢失。可以在求解结束后跑一行sum(y, 2)查看随时间的变化偏差超过总人口的1%就得回头查微分方程里的正负号。plot_curves.m会自动把五条曲线画在同一张图上但默认的纵轴范围可能让D和R曲线因为数值太小而压扁。我会习惯性把D曲线单独画一个对数坐标的子图否则D的累积趋势很难看清。5. 避坑指南MATLAB模型调试中的5个实际问题5.1 仿真结果全是NaN现象跑完main_seird.m后y矩阵里出现大量NaN。原因初始感染者I0设为0导致感染项beta * S * I / N在第一步就是0所有仓室保持初始值不变但数值求解器的自适应步长尝试把时间步调大后某个仓室计算出现除以零的情况触发NaN连锁扩散。解决I0至少设为1。还要检查beta0是否过大比如超过1.5在ode45默认容差下容易数值爆炸。把beta0回退到0.3附近NaN现象基本就会消失。5.2 改了参数但输出结果完全不变现象在编辑器中修改beta0后运行绘图结果和上一次一模一样。原因MATLAB的脚本运行默认会沿用上一次工作区中残留的变量。如果某次运行后main_seird.m里的赋值语句被跳过比如代码报错后继续执行旧参数就留在工作区里被重复使用。解决在每个主脚本开头都写clear; clc;并且使用run(main_seird.m)强制刷新。注意clear会清空所有变量如果依赖外部导入的Excel数据需要放在clear之后重新读取。5.3 消毒强度调得很大曲线却几乎没变化现象把dis_amp从0.2调到0.9感染峰值只下降了不到5%。原因消毒效果被指数衰减项exp(-0.1*t)主导当t超过30天后衰减系数已经降到0.05以内消毒强度再大也推不动曲线。解决缩短消毒作用窗口把衰减系数从0.1改到0.05或者直接把指数衰减换成矩形窗函数这样消毒强度在作用窗口内保持恒定效果才会真实反映在曲线上。还有一种做法是把dis_amp从标量改成向量配合schedule使用分段常数控制。5.4 绘图时D曲线被压成一条直线现象plot_curves.m生成的图中D和R曲线贴在横轴上毫无辨识度。原因D曲线数值比I曲线小两个数量级线性坐标下几乎不可见。这是SEIRD可视化的通病不是模型出错。解决用subplot把D单独画到第二个子图并给D曲线一个独立的纵轴。也可以改用yyaxis left和yyaxis right把I和D分别放在左右两侧。示例如下figure; yyaxis left; plot(t, y(:,3), LineWidth, 2); ylabel(Infectious); yyaxis right; plot(t, y(:,5), LineWidth, 2); ylabel(Deceased);5.5 多次运行每次结果都有微小差异现象同样参数运行两次峰值时间有时差一天曲线尾部有细微抖动。原因如果代码里用了rand或randn生成初始扰动而随机种子没有固定结果自然无法复现。这个问题在接触矩阵版中尤其常见。解决在main_seird.m开头固定随机种子rng(2025);这一行加上之后任何随机因素都会被锁死多轮实验才能公平对比。6. 怎么改出第二个模型二次开发的三个入手方向第一个入手方向是把单一人群改成多群体模型。现在contact_matrix.m只是预留了接口没有真正启用。要改成多群体需要给每个子群单独分配一组beta和gamma。以学生和老年人两个子群为例学生群体的接触频率高但病程短老年人群体反之。将接触矩阵写成2x2的数值表每个元素表示两个群体的接触次数然后主模型中的传染项从beta * S * I / N改成beta * C(i,j) * S_i * I_j / N需要同时更新所有仓室的下标引用。完成后可以观察到一个典型现象感染的峰值不再是单峰而是出现拖尾的双峰形态。第二个方向是引入时变参数。真实疫情中防控措施的强度是逐周调整的而不是一沉不变。可以把beta0从一个标量改成向量向量的每个元素对应一天的传染率。比如beta_daily linspace(0.32, 0.2, 101);用interp1在求解过程中动态插值。这个改法特别适合复盘某次特定疫情把每日真实传染率曲线直接喂进模型再对比模型输出的感染人数和真实报告数据使用最小二乘法拟合出最优参数组合。第三个方向是把输出从曲线改成指标表。现在的输出只是mat文件真正写报告时还需要峰值日期、累计病例、总死亡数、第90天残余感染数五个指标。我习惯在主脚本末尾加一段汇总代码peak_I max(y(:,3)); peak_day t(y(:,3) peak_I); total_R y(end, 4); total_D y(end, 5); fprintf(峰值感染数: %.0f 第%.1f天\n, peak_I, peak_day); fprintf(累计恢复: %.0f 累计死亡: %.0f\n, total_R, total_D);这段代码让每次仿真的结果一眼就能读取不用再打开绘图窗口肉眼找峰值。跑完一组参数扫描后把这些指标收集成一个表格消毒策略的优化问题就变成了一张可以直观对比的Excel表格。最后说一个我反复踩过的坑。以前我总喜欢在seird_ode.m里直接改参数做试验结果每次临时改完忘记还原后续仿真全都被带偏。后来我强制规定方程文件只保留模型结构所有参数必须由main_seird.m传入任何实验性修改一律写在主脚本里跑完立刻save。从那以后每次做参数扫描我都在运行前先clear清场、固定随机种子、校验N的守恒确认全绿才去看曲线。这套流程虽然多花十秒但能省掉后面至少两小时的排查时间。希望这篇解析能帮你顺利跑通这套SEIRD模型并把消毒网策略真正嵌进你的仿真流程里。本文还有配套的精品资源点击获取