简介：DESSABNeT 是一个基于 MATLAB 的传染病传播模拟软件包，以基于代理的建模（ABM）方式模拟 COVID-19 等传染病在城市或区域人口中的扩散，适合计算流行病学研究者、公共卫生建模人员及具备一定 MATLAB 基础的仿真爱好者学习与二次开发。资源共 42 个文件，以 33 个 .m 脚本为主，涵盖人口结构生成、家庭/社交/工作网络构建、每周接触日程、社交限制阶段设置、动态传播模块、R0 计算和结果输出等完整流程；另有 7 个 .mat 预置数据文件、1 个悉尼新冠病例 xlsx 数据表和 README 说明文档，压缩包仅 2.7MB，便于快速浏览与调试。目前已有 712 人学习下载。借助这套代码，读者可以系统了解基于代理的传染病模型如何从人口统计特征生成人工社会，并通过模拟每日接触来研究不同社交限制措施对疫情传播的影响，为后续扩展城市级 SEIRD 模型或开展防控策略对比实验提供落地方案。
1. SEIRD模型为什么难调：一份能直接跑的MATLAB仿真框架
如果你拿到过SEIRD模型的MATLAB代码，大概率经历过这种尴尬：换了数据集、调了参数，跑出来的感染曲线要么早早冲顶后“清零”，要么直接爆掉满屏NaN。我最初拆这个项目，也是因为某高校课程设计里需要做疫情传播与消毒干预的模拟，手写了三版SEIR都不收敛，最后才发现问题不出在算法，而出在对模型结构和干预策略的建模方式上。这份名为DESSABNeT的资源，正好补上了这个空档：它把SEIRD的五个状态转移、消毒网络的干预逻辑和仿真输出脚本打包在一起，代码按系统开源的方式组织，每个模块独立成文件，适合正在做传染病建模课程设计、毕业设计，或者需要对消毒策略做快速预演的从业者直接复现。
2. DESSABNeT的模型结构与代码骨架：五个状态如何串联成一条传播链
2.1 SEIRD的状态划分与微分方程组
SEIRD与经典SEIR的区别在于多了一个D（死亡）状态，让模型能够同时刻画康复与死亡两个出口。DESSABNeT的核心是以下六个微分方程：
function dydt = seird_derivs(t, y, beta, sigma, gamma, mu, N) % y = [S; E; I; R; D] S = y(1); E = y(2); I = y(3); R = y(4); D = y(5); dS = -beta * S * I / N; dE = beta * S * I / N - sigma * E; dI = sigma * E - gamma * I - mu * I; dR = gamma * I; dD = mu * I; dydt = [dS; dE; dI; dR; dD]; end逻辑说明：S（易感者）与I（感染者）接触后以betaSI/N的速率转为E（暴露者）；E以sigma速率转为I（发病）；I再分为两个出口，gamma比例康复进入R，mu比例死亡进入D。这里的关键是把I看成“存量”，所有转移项都乘以I，而不是直接给I一个固定的生存时间。
参数说明：beta是有效接触率，直接与消毒干预挂钩，剂量越大beta越低；sigma=1/潜伏期；gamma=1/感染期；mu是死亡速率而非病死率本身，这一点在后文避坑章会展开。N是总人口，实现时建议用常数，而不是让N随时间变化。
2.2 代码文件结构与各模块职责
DESSABNeT的资源内部不是单脚本跑通，而是按典型仿真工程拆成五个文件。拆包后你会看到类似这样的组织方式：
| 文件 | 职责 |
|---|---|
seird_derivs.m | 定义微分方程组，上面的代码就是它 |
run_simulation.m | 主脚本，负责参数设置、求解循环和结果输出 |
intervention.m | 消毒干预策略模块，把干预强度折算成beta的时变系数 |
plot_results.m | 绘图脚本，出S/E/I/R/D五条曲线 |
calibrate_params.m | 参数校准函数，用于拟合实际疫情曲线 |
这种分文件写法的好处是把“模型定义”和“策略定义”分开。比如你要改消毒策略，只动intervention.m，不需要碰微分方程；你要换数据集，也只改run_simulation.m里的参数初始化段。我在课程设计里提交的版本，就是在intervention.m里加了一个周期性消毒强度数组，其余文件完全没动。
2.3 从SEIR到SEIRD的扩展逻辑：死亡率怎么接进去
很多开源的SEIR代码只有四个状态，想加死亡不能简单把mu塞进gamma里。DESSABNeT的处理方式是让I同时流向R和D，核心是gamma与mu的比例关系。如果已知一个感染者的病死率CFR，合理的拆法为：
CFR = 0.02; % 病死率，单位是百分比 recovery_time = 14; % 平均感染期，单位是天 gamma = 1 / recovery_time; mu = gamma * CFR / (1 - CFR);逻辑说明：gamma与mu之和代表感染者的“退出速率”，CFR是退出时死亡的比例，所以mu = (gamma+mu) * CFR。把gamma表达式代进去整理后，就得到上面的mu计算公式。注意这里CFR用的是比例而不是百分数，配置2%时写0.02。
参数说明：如果不做这个拆分，直接把mu写成一个绝对值，那么感染期长度和病死率会被耦合，导致你调高病死率时，平均感染期也会被篡改。这也是DESSABNeT比普通SEIR代码更“可信”的一个细节。
3. 跑通第一个仿真：从初始化到看懂疫情曲线
3.1 最小可运行示例：用RK4手动迭代
run_simulation.m里用的是定步长四阶Runge-Kutta（RK4），好处是步长可控、调试直观，不会像ode45那样因为自适应步长在干预时点产生奇怪跳点。完整的最小主脚本如下：
% run_simulation.m % SEIRD仿真主脚本，使用RK4求解 N = 100000; % 总人口 I0 = 10; % 初始感染者 E0 = 50; % 初始暴露者 S0 = N - E0 - I0; % 初始易感者 y0 = [S0; E0; I0; 0; 0]; % 初始状态向量 beta = 0.35; % 有效接触率 sigma = 1/5.2; % 潜伏期约5.2天 gamma = 1/14; % 感染期约14天 mu = 0.004; % 死亡速率 T = 180; dt = 0.1; % 总时长180天，步长0.1天 t = 0:dt:T; n = length(t); Y = zeros(5, n); Y(:,1) = y0; for k = 1:n-1 yk = Y(:,k); k1 = seird_derivs(t(k), yk, beta, sigma, gamma, mu, N); k2 = seird_derivs(t(k)+dt/2, yk+dt/2*k1, beta, sigma, gamma, mu, N); k3 = seird_derivs(t(k)+dt/2, yk+dt/2*k2, beta, sigma, gamma, mu, N); k4 = seird_derivs(t(k)+dt, yk+dt*k3, beta, sigma, gamma, mu, N); Y(:,k+1) = yk + dt/6*(k1+2*k2+2*k3+k4); end plot_results(t, Y, N);逻辑说明：每一步按RK4的经典格式算出k1到k4，把四个斜率加权平均得到增量。主循环里没有直接调用MATLAB的ode求解器，这样你能清楚地看到每步的状态变化，调试时也容易定位是哪个状态先失真。
参数说明：dt=0.1是我习惯的初始值，如果改成1可能会震荡；T=180适合模拟一轮暴发过程。N保持固定值，五个状态的初始和必须等于N，否则容易触发避坑章里的守恒问题。
3.2 参数设置：beta、sigma、gamma、mu的物理含义与初始范围
第一次跑通不要追求拟合真实数据，用一组“看起来合理”的初值即可。我把常用范围整理如下：
| 参数 | 物理含义 | 常见范围 | 备注 |
|---|---|---|---|
| beta | 单位时间有效接触数 | 0.2-0.6 | 消毒强度越高beta越低 |
| sigma | 潜伏期的倒数 | 1/7到1/3 | 单位：每天 |
| gamma | 康复速率 | 1/21到1/7 | 与感染期互为倒数 |
| mu | 死亡速率 | 0.001-0.01 | 非CFR，按2.3公式转换 |
其中mu最容易配置错误。很多人直接把病死率填进mu，比如写到0.02，这会让感染者退出速率凭空多出一大截，感染期变短，峰值也会提前。我在复现时习惯先把CFR写在注释里，再用公式转成mu，避免回头忘了这个值的口径。
3.3 输出与可视化：从曲线中定位干预窗口
plot_results.m里的绘图代码负责把五个状态画在同一张图上：
function plot_results(t, Y, N) figure('Color','w'); plot(t, Y(1,:), 'b-', 'LineWidth', 1.5); hold on; plot(t, Y(3,:), 'r-', 'LineWidth', 2); plot(t, Y(4,:), 'g-', 'LineWidth', 1.5); plot(t, Y(5,:), 'k--', 'LineWidth', 1.5); legend({'S','I','R','D','Location','NorthEast'}); xlabel('天'); ylabel('人数'); grid on; title('SEIRD疫情曲线'); end逻辑说明：图中最值得关注的是红色I曲线和黑色D曲线。I曲线的峰值位置告诉你暴发高峰在什么时候，D曲线的斜率告诉你医疗资源压力有多陡。消毒干预的成效主要体现在I曲线峰值是否被压低或推迟。
参数说明：如果只想看感染规模，可以只画I和D两条线。我一般会在I峰值附近加一个垂直虚线标注时间，方便后面对比不同消毒策略的位移量。
注意：如果曲线在到达峰值前出现抖动，优先怀疑步长过大或参数mu接近上限，而不是模型本身写错了。
4. 消毒干预模拟的避坑清单：五个最容易翻车的配置细节
4.1 干预强度把beta折成负数
现象：模拟中途S变成负值，或命令行直接出现warning。
原因：intervention.m里把干预强度写成了beta * (1 - intensity)，而intensity在某个时间段超过1。消毒网络覆盖率达到120%的场景在代码里本来不该出现，但参数扫描时有人会把intensity填到1以上做“极限测试”。
解决：在intervention函数里加截断判断，把intensity限制在0到0.95之间。保留上限0.95也是合理的，因为现实中很难做到100%阻断传播，低于0.95的beta折算等价于保留基础传播通道。
4.2 总人口不守恒，五个状态加起来越跑越偏
现象：随着t增长，S+E+I+R+D的结果从N逐渐漂移到N加减几百。
原因：初始状态赋值时写了S0 = N - E0 - I0，这是对的；但又在某处将N设为随时间变化的变量，比如让它等于当前存活人数，导致后面所有比例项都跟着变。
解决：把N固定为常量，放在参数初始化段，后续只用于比例计算，不参与更新。另外可以在主循环最后加一条断言：
assert(abs(sum(Y(:,k+1)) - N) < 1e-6, '总人口不守恒');这条断言能让问题在第一步就跑出来，不需要等到曲线画完再人肉找原因。我在某次模拟中，就是靠它定位到在intervention里误改了N。
4.3 步长太大，RK4解出负的人数
现象：E曲线在某几个时间点出现负值，或者整体呈现锯齿状波动。
原因：dt设得太大。RK4的局部截断误差是O(dt^4)，但稳定性边界不是无限的。当dt=1时，sigma=1/3这类参数组合已经接近震荡边界；如果dt=2，基本必翻车。
解决：把dt降到0.1以下，或者改用MATLAB的ode45让求解器自适应步长。手动RK4的优势是确定性和可调性，代价就是必须无条件接受dt≤0.1的限制。
4.4 潜伏期参数写反，把sigma直接填成天数
现象：I曲线峰值比预期提前了接近一个潜伏期，你调低beta也不管用。
原因：把sigma = 1/5.2写成了sigma = 5.2。5.2是被期望输入的，但单位反了。
解决：在参数表注释里明确写成“sigma = 1/潜伏期天数”，并在代码前面加一个latent_period = 5.2; sigma = 1/latent_period;。这种写法能让后来者一眼看出物理含义，而不是面对一个神奇数字。
4.5 直接套文献参数，不做本地校准
现象：模型曲线与实际疫情数据相差巨大，预测峰值方向都反了。
原因：beta在文献里对应的社会接触模式、口罩佩戴率、消毒执行度与你模拟的场景完全不同。某地区的beta=0.35，换个城市可能只有0.2。
解决：至少要做一次粗略的手动校准。打开calibrate_params.m，把实际累计确诊数或每日新增放进去，先固定sigma和gamma，只调beta和mu，让模型峰值和真实峰值错位不超过一周。这一条是血泪经验，我在某次课程设计里直接照抄某篇论文的参数，结果模拟曲线在20天就封顶，而实际数据60天才到峰，从一开始就不该指望异地参数能直接落地。
5. 从复现到改进：用最速下降法校准模型参数
5.1 目标函数：把“曲线长得像”翻译成数值
手动校准到一定程度后，就该用数值优化来收敛参数。DESSABNeT的calibrate_params.m用的是MATLAB自带fminsearch，目标函数是模拟感染人数与实际感染人数的均方根误差（RMSE）。具体思路是先跑一遍仿真，再计算误差，把误差返回给优化器去调整参数。
function mse = calibration_error(params, real_I, t, N) beta = params(1); mu = params(2); sigma = 1/5.2; gamma = 1/14; dt = 0.1; y0 = [N-60; 50; 10; 0; 0]; Y = run_seird(y0, beta, sigma, gamma, mu, N, t, dt); sim_I = Y(3,:); mse = sqrt(mean((sim_I - real_I).^2)); end逻辑说明：real_I是实际每日在院或每日新增感染的人数序列，sim_I从仿真结果中抽取。RMSE越小代表模拟曲线与真实曲线贴合度越高。固定sigma与gamma是合理的，因为潜伏期和感染期在不同地区相对稳定，真正需要校准的是接触率与死亡速率。
参数说明：fminsearch对初值敏感，beta初值取0.3，mu取0.005，一般能在30次迭代内收敛。如果优化结果出现beta接近0或mu为负，说明真实数据本身波动太大，需要先对原始数据做七天滑动平均。
5.2 校准流程与初始值技巧
校准的完整流程我一般走三步：
第一步，把真实数据做平滑，去掉周末效应和报告延迟造成的锯齿。 第二步，固定sigma和gamma，只优化beta和mu，避免高维参数空间里的局部最优。 第三步，用优化后的参数回代仿真，对比累计死亡和峰值时间两个指标，若累计死亡对不上，再把mu放出来联合优化。
这里有个经验：不要一上来就同时校准四个参数。SEIRD模型里beta、sigma、gamma之间存在明显的相关性，数据量不够时，多参数联合优化会得到一组“数值好看但物理荒谬”的参数，比如潜伏期变成20天。每次只放开一两个参数，是保证结果可信的关键技巧。
5.3 验证模型稳定性：三次运行看方差
校准完毕不等于工作结束。我会把优化得到的参数作为初值，然后对初值做上下10%的扰动，重复跑三次仿真，看峰值位置和累计死亡的变化范围。如果三次结果差异超过15%，说明参数对初值太敏感，可能是数据质量不足，也可能是模型本身对这个场景缺乏辨识度。
从那以后，我每次拿到新的疫情数据做模拟时，都强制走一遍“先固定sigma和gamma，再校准beta，最后联合优化mu”的流程，再也不敢拿论文参数直接开跑。校准出来的结果哪怕和文献不一致，至少曲线趋势是可信的。这份DESSABNeT框架帮我省掉了大量重复造轮子的时间，如果你也在被SEIRD的参数调优折磨，希望帮到你。
本文还有配套的精品资源，点击获取
