传染病建模实战：SEIR仓室模型与参数校准全流程

1. 从数据到决策传染病建模到底在算什么1.1 一个建模项目该从哪里下手传染病建模这件事听起来像是流行病学家的专属领域但实际上它离我们并不远。任何一个做数据分析、做系统仿真、甚至做游戏数值策划的人都能在这套方法论里找到熟悉的味道。我最初接触这类项目的时候第一反应是这不就是一个带状态转移的动力学系统吗把人群分成几个仓室设定转移规则跑时间步看曲线怎么走。后来发现事情比想象中复杂得多但也比想象中有趣得多。这个项目的核心目标很明确用数学模型去回答关于新型冠状肺炎的几个关键问题。比如病毒的基本传播能力有多强如果不管控感染曲线会飙到什么程度如果采取不同强度的干预措施效果差多少医疗资源在什么时间点会被击穿这些问题靠直觉回答不了靠拍脑袋更不行必须用模型去推演。适合谁来参考这篇内容我认为有三类人。第一类是有一定编程基础、想入门传染病建模的数据从业者你不需要有医学背景但需要理解微分方程或者离散仿真的基本逻辑。第二类是公共卫生或政策研究方向的从业者你可能不写代码但需要理解模型输出的含义和局限才能做出合理判断。第三类是对复杂系统仿真感兴趣的工程师传染病模型本质上是一个典型的复杂系统问题涉及非线性动力学、参数敏感性、不确定性量化这些技能在别的领域同样通用。1.2 为什么选仓室模型而不是纯统计拟合建模路线有很多条最主流的两条路是一条是仓室模型比如SIR、SEIR及其变体另一条是纯统计拟合比如用时间序列模型去拟合确诊数据。这两条路各有优劣我最终选了仓室模型为主、统计方法为辅的混合路线原因如下。纯统计拟合的优势是简单直接你有确诊数据直接上ARIMA或者LSTM去拟合预测就行了。但它有个致命问题它只能告诉你“数据长什么样”不能告诉你“为什么长这样”。一旦干预措施发生变化历史数据的模式就被打破了纯统计模型的预测会严重失准。而仓室模型的优势恰恰在于它把传播机制显式地写进了方程里你可以通过调整参数来模拟不同干预场景这是纯统计方法做不到的。当然仓室模型也有它的代价。它对数据的要求更高你需要估计很多参数而这些参数往往没有直接观测值。比如传播率、潜伏期、无症状比例这些量在早期数据里很难准确估计。所以我的做法是用仓室模型搭建仿真框架用统计方法去估计和校准参数两者结合。注意不要一上来就追求最复杂的模型。我见过不少人直接上SEIR年龄分层空间扩散时变参数结果参数多到根本标定不了最后跑出来的结果还不如最简单的SIR靠谱。建模的第一原则是从简单开始逐步增加复杂度每增加一个参数都要问自己——这个参数我能估出来吗1.3 核心思路的拆解整个项目的思路可以拆成四步定义状态变量把人群划分为若干仓室比如易感者、暴露者、感染者、康复者、死亡者。每个仓室的人数随时间变化。建立转移方程定义仓室之间的转移规则和速率。比如易感者以某个速率变成暴露者暴露者经过潜伏期变成感染者感染者以一定概率康复或死亡。参数估计与校准用实际数据去反推模型参数使得模型输出与观测数据尽可能吻合。情景模拟与敏感性分析调整关键参数模拟不同干预策略下的结果分析哪些参数对结果影响最大。这四步看起来简单但每一步都有大量的细节和坑。接下来我会逐一展开把每个环节的关键技术点、实操方法和踩坑经验都讲清楚。2. 模型搭建的核心细节与参数体系2.1 仓室划分SEIR还是更复杂的变体最基础的传染病模型是SIR只有三个仓室易感者、感染者、康复者。这个模型适合那些感染后立刻具有传染性、且康复后获得免疫力的疾病。但新型冠状肺炎有一个显著特点存在潜伏期且潜伏期内可能具有传染性。这就需要用SEIR模型多一个暴露者仓室。SEIR的基本结构是这样的SSusceptible易感人群还没被感染也没有免疫力。EExposed暴露者已经接触了病毒但还没出现症状处于潜伏期。IInfectious感染者有症状且具有传染性。RRecovered康复者获得免疫力或死亡取决于是否单独建模。但实际项目中我很快发现基础SEIR不够用。原因有几个第一无症状感染者也有传染性需要单独建模第二确诊和实际感染之间存在时间差检测能力会影响数据第三死亡和康复需要分开追踪。所以我最终用的是一个扩展版SEIR增加了无症状感染者仓室和死亡仓室。具体来说我的仓室划分是仓室含义关键参数S易感者初始值、总人口E暴露者潜伏期潜伏期时长I_s有症状感染者传染率、症状期时长I_a无症状感染者传染率通常低于有症状者R康复者康复率D死亡者病死率这个划分比基础SEIR复杂但比那些动辄十几个仓室的模型简单得多。我觉得这是一个比较好的平衡点既能捕捉关键特征又不至于参数多到标定不了。实操心得仓室不是越多越好。每增加一个仓室你就多了一组需要估计的参数。如果你的数据只有每日确诊数那你最多能可靠估计两三个参数。仓室多了参数就不可辨识了模型跑出来的结果看着漂亮其实全是过拟合。2.2 转移方程连续还是离散仓室之间的转移可以用连续时间的微分方程来描述也可以用离散时间的差分方程来描述。两者在数学上是等价的当时间步足够小的时候但在实现上有很大区别。连续时间微分方程的形式是这样的dS/dt -beta * S * I / N dE/dt beta * S * I / N - sigma * E dI/dt sigma * E - gamma * I dR/dt gamma * I其中beta是传播率sigma是潜伏期倒数gamma是康复率倒数。离散时间差分方程则是S(t1) S(t) - beta * S(t) * I(t) / N E(t1) E(t) beta * S(t) * I(t) / N - sigma * E(t) I(t1) I(t) sigma * E(t) - gamma * I(t) R(t1) R(t) gamma * I(t)我最终选了离散时间实现原因很实际数据是按天报告的。连续时间模型虽然数学上更优雅但你需要用数值积分去求解而且和按天报告的数据对齐时反而更麻烦。离散模型直接以天为步长和数据结构天然匹配实现起来也简单得多。当然离散模型有个前提时间步要足够小使得每个步长内状态变化不会太大。对于日粒度数据来说这个条件基本满足。但如果你的数据是周粒度或者月粒度那离散模型的误差就会比较大这时候还是用连续模型更稳妥。2.3 关键参数的含义与估算方法模型里最关键的几个参数我逐个说一下含义和估算方法。传播率beta这是最重要的参数表示一个感染者在单位时间内传染给易感者的概率乘以接触人数。beta越大传播越快。估算方法通常有两种一种是从早期指数增长阶段的增长率反推另一种是用优化算法去拟合确诊数据。我用的是后者因为早期数据噪声太大直接反推误差很高。潜伏期倒数sigmasigma等于1除以平均潜伏期。比如平均潜伏期是5天那sigma就是0.2。这个参数通常可以从文献里查到不需要自己估计。但要注意潜伏期的定义在不同研究里可能不一样有的指从暴露到出现症状有的指从暴露到具有传染性用的时候要统一。康复率倒数gammagamma等于1除以平均感染期。这个参数也是从文献里查比较靠谱。但有个坑感染期的定义也有歧义有的指从出现症状到康复有的指从具有传染性到不再具有传染性。我建议在项目文档里明确写清楚每个参数的定义不然后面自己都会搞混。基本再生数R0R0 beta / gamma表示在完全没有免疫力的群体中一个感染者平均能传染多少人。R0是模型输出的一个综合指标不是直接输入的参数。很多人会把R0和beta搞混其实R0是beta和gamma的组合。病死率CFR死亡人数除以确诊人数。这个参数在模型中通常作为一个比例决定感染者中有多少会死亡。但要注意CFR的计算受检测能力影响很大检测不足的时候CFR会被高估。注意参数估计不是一锤子买卖。随着数据积累你需要不断更新参数估计。我建议把参数估计做成一个可重复的流程每次有新数据就重新跑一遍而不是手动改参数。2.4 初始条件被低估的关键因素很多人把注意力全放在参数上却忽略了初始条件的重要性。初始条件包括初始感染人数、初始暴露人数、总人口、初始免疫比例等。这些值如果设得不对模型输出会差很多。初始感染人数尤其关键。在疫情早期实际感染人数往往远高于报告确诊人数因为检测能力有限。如果你直接用报告确诊数作为初始感染人数模型会严重低估传播速度。我的做法是用报告确诊数乘以一个检测覆盖率倒数来估算实际感染人数。比如你估计检测覆盖率是20%那实际感染人数大约是报告数的5倍。总人口这个参数看起来简单但也有讲究。如果你建模的是一个城市那总人口就是城市常住人口。但如果你建模的是一个国家那就要考虑人口流动的影响可能需要引入空间分层。我一开始做的是单城市模型总人口设为常数后来发现人口流动对结果影响很大才加了流入流出项。初始免疫比例在疫情早期通常设为0因为新病毒大家都没有免疫力。但如果你的模型是从疫情中期开始跑的那就要考虑已经感染过的人群初始免疫比例不能忽略。3. 实操过程从数据到仿真到结果解读3.1 数据准备与清洗数据是建模的燃料。没有好的数据再精巧的模型也是空中楼阁。我用的数据主要是每日新增确诊数、每日新增死亡数、每日康复数。这些数据通常来自公开的疫情通报但直接拿来用会有很多问题。第一个问题是报告延迟。确诊数据往往有延迟今天确诊的人可能是几天前就感染的。这个延迟会使得模型拟合出现偏差。我的处理方法是在模型里加一个观测延迟项把模型预测的感染人数经过一个延迟函数映射到报告确诊数再去拟合。第二个问题是周末效应。很多地方的检测和报告在周末会减少导致数据出现周期性波动。如果不处理模型会把这个波动当成真实的传播变化。我的做法是对数据做7天滑动平均平滑掉周末效应。但滑动平均会引入滞后所以只用来做参数估计不用来做实时预测。第三个问题是数据修正。疫情数据经常会被回溯修正今天看到的数字和一周前看到的可能不一样。我建议在项目里保存数据快照每次拉取数据都记录时间戳这样后面回溯的时候知道用的是哪个版本的数据。数据清洗的代码大概长这样import pandas as pd import numpy as np # 读取原始数据 raw pd.read_csv(daily_cases.csv, parse_dates[date]) # 按日期排序 raw raw.sort_values(date).reset_index(dropTrue) # 计算7天滑动平均 raw[new_cases_smooth] raw[new_cases].rolling(window7, centerTrue).mean() # 处理缺失值用前向填充 raw[new_cases_smooth] raw[new_cases_smooth].fillna(methodffill) # 计算累计确诊 raw[cum_cases] raw[new_cases].cumsum() # 剔除异常值单日新增超过均值5倍标准差的点 mean raw[new_cases].mean() std raw[new_cases].std() raw raw[np.abs(raw[new_cases] - mean) 5 * std] print(f清洗后数据量: {len(raw)}) print(raw.head(10))这段代码看起来简单但每一步都有讲究。滑动平均的窗口选7天是因为一周有7天能完整覆盖一个周期。异常值剔除用5倍标准差是因为疫情数据本身波动就大阈值设太小会误删真实信号。3.2 模型实现用Python搭建SEIR仿真模型实现我选的是Python主要因为生态好numpy做数值计算、scipy做优化、matplotlib做可视化一条龙。下面是我用的核心仿真代码。import numpy as np from scipy.integrate import odeint def seir_model(y, t, N, beta, sigma, gamma): 基础SEIR模型的微分方程 y: [S, E, I, R] t: 时间 N: 总人口 beta: 传播率 sigma: 潜伏期倒数 gamma: 康复率倒数 S, E, I, R y dSdt -beta * S * I / N dEdt beta * S * I / N - sigma * E dIdt sigma * E - gamma * I dRdt gamma * I return [dSdt, dEdt, dIdt, dRdt] # 参数设置 N 1000000 # 总人口 beta 0.5 # 传播率 sigma 0.2 # 潜伏期倒数平均潜伏期5天 gamma 0.1 # 康复率倒数平均感染期10天 # 初始条件 I0 10 # 初始感染者 E0 20 # 初始暴露者 R0 0 # 初始康复者 S0 N - I0 - E0 - R0 # 时间范围 t np.linspace(0, 200, 200) # 求解 y0 [S0, E0, I0, R0] solution odeint(seir_model, y0, t, args(N, beta, sigma, gamma)) # 提取结果 S solution[:, 0] E solution[:, 1] I solution[:, 2] R solution[:, 3] # 计算R0 R0 beta / gamma print(f基本再生数 R0 {R0:.2f}) print(f峰值感染人数 {I.max():.0f}) print(f峰值出现时间 第{t[I.argmax()]:.0f}天)这段代码跑出来的结果R0等于5峰值感染人数大约占总人口的20%左右峰值出现在第60天前后。这个结果和很多早期研究的估计是吻合的。但实际项目中我用的是扩展版SEIR增加了无症状感染者和死亡仓室。代码会更长但结构是一样的。关键是要把每个仓室的转移关系理清楚不要搞混。实操心得写仿真代码的时候一定要加单元测试。比如你可以手动算一个简单场景的结果然后和代码输出对比。我一开始没做测试结果有个仓室的转移方向写反了跑出来的曲线看着挺合理实际上完全错了。后来加了测试才发现。3.3 参数校准用优化算法拟合实际数据参数校准是整个项目里最耗时的环节。我的做法是定义一个损失函数衡量模型输出和实际数据的差距然后用优化算法去最小化这个损失函数。损失函数我选的是负对数似然而不是简单的均方误差。原因是疫情数据通常服从泊松分布或负二项分布用似然函数更合理。具体来说如果每日新增确诊数服从泊松分布那对数似然就是log L sum(y_obs * log(y_pred) - y_pred - log(y_obs!))其中y_obs是实际观测值y_pred是模型预测值。优化算法我用的是L-BFGS-B因为它支持边界约束可以限制参数在合理范围内。比如beta必须在0到1之间sigma必须在0到1之间。from scipy.optimize import minimize from scipy.stats import poisson def loss_function(params, data, N): 负对数似然损失函数 params: [beta, sigma, gamma] data: 实际观测的每日新增确诊数 N: 总人口 beta, sigma, gamma params # 参数边界检查 if beta 0 or sigma 0 or gamma 0: return 1e10 # 运行模型 I0 data[0] if data[0] 0 else 1 E0 I0 * 2 S0 N - I0 - E0 y0 [S0, E0, I0, 0] t np.arange(len(data)) try: solution odeint(seir_model, y0, t, args(N, beta, sigma, gamma)) except: return 1e10 I_pred solution[:, 2] # 计算每日新增用感染者的变化率近似 new_infections np.diff(I_pred, prependI_pred[0]) new_infections np.maximum(new_infections, 1e-6) # 避免log(0) # 负对数似然 nll -np.sum(poisson.logpmf(data, new_infections)) return nll # 初始参数猜测 initial_params [0.3, 0.2, 0.1] # 参数边界 bounds [(0.01, 2.0), (0.05, 1.0), (0.01, 1.0)] # 优化 result minimize( loss_function, initial_params, args(observed_data, N), methodL-BFGS-B, boundsbounds ) print(f最优参数: beta{result.x[0]:.4f}, sigma{result.x[1]:.4f}, gamma{result.x[2]:.4f}) print(fR0 {result.x[0]/result.x[2]:.2f})这段代码有几个关键点。第一初始感染人数不能设为0否则模型永远跑不起来我用的是数据第一天的值。第二预测的每日新增用感染者变化率来近似这在离散数据上是合理的。第三参数边界要设得合理太宽会导致优化跑飞太窄会错过最优解。参数校准的结果通常不是唯一的。不同的参数组合可能给出相似的拟合效果这就是所谓的参数不可辨识性。解决方法是固定一些从文献中能查到的参数比如潜伏期只估计最关键的几个参数比如beta和gamma。或者用贝叶斯方法给出参数的后验分布而不是单点估计。3.4 情景模拟不同干预策略的效果对比模型校准好之后就可以做情景模拟了。这是整个项目最有价值的部分因为你可以回答“如果……会怎样”的问题。我设计了四种情景情景干预措施参数调整无干预什么都不做beta保持不变轻度干预减少50%接触beta乘以0.5中度干预减少70%接触beta乘以0.3强力干预减少90%接触beta乘以0.1模拟结果非常直观无干预情况下峰值感染人数占总人口的20%以上医疗系统必然被击穿。轻度干预能把峰值压到10%左右但仍然很危险。中度干预能把峰值压到3%以下基本可控。强力干预能把峰值压到1%以下但需要维持很长时间。这里有个关键概念叫干预时机。同样的干预强度早做和晚做效果差很多。我做了对比模拟发现如果在感染人数还很少的时候就采取强力干预可以把疫情控制在很小的规模。但如果等到感染人数已经很多了再干预效果会大打折扣因为传播已经形成了惯性。注意情景模拟的结果高度依赖模型假设。如果你的模型没有考虑超级传播事件、没有考虑人口异质性、没有考虑行为反馈那模拟结果只能作为参考不能作为决策的唯一依据。我建议在报告结果的时候一定要附上敏感性分析说明哪些假设对结果影响最大。3.5 敏感性分析哪些参数最关键敏感性分析是检验模型稳健性的重要手段。我用的方法是拉丁超立方采样在参数空间里均匀采样然后看每个参数对输出的影响。具体做法是对每个参数在其合理范围内采样100个值其他参数固定在中位数然后跑模型记录峰值感染人数和峰值时间。最后计算每个参数的偏相关系数相关系数越大说明该参数对结果影响越大。我的分析结果是传播率beta对峰值感染人数的影响最大相关系数超过0.9。潜伏期sigma的影响次之康复率gamma的影响相对较小。这个结果符合直觉传播越快峰值越高。但有个反直觉的发现潜伏期越长峰值反而可能越高。原因是潜伏期长意味着感染者被发现和隔离的时间更晚传播窗口更长。这个发现对防控策略有重要启示快速检测和隔离比单纯减少接触更重要。敏感性分析的代码大概长这样from scipy.stats import spearmanr import itertools # 参数范围 beta_range np.linspace(0.1, 1.0, 20) sigma_range np.linspace(0.1, 0.5, 20) gamma_range np.linspace(0.05, 0.3, 20) # 存储结果 results [] for beta, sigma, gamma in itertools.product(beta_range, sigma_range, gamma_range): # 运行模型 y0 [N - 100, 50, 50, 0] t np.linspace(0, 300, 300) solution odeint(seir_model, y0, t, args(N, beta, sigma, gamma)) I solution[:, 2] peak_infection I.max() peak_time t[I.argmax()] results.append({ beta: beta, sigma: sigma, gamma: gamma, peak_infection: peak_infection, peak_time: peak_time }) df pd.DataFrame(results) # 计算Spearman相关系数 for param in [beta, sigma, gamma]: corr, p_value spearmanr(df[param], df[peak_infection]) print(f{param} vs 峰值感染人数: 相关系数{corr:.3f}, p值{p_value:.4f})这段代码跑完你会得到一张相关系数表清楚地看到每个参数的影响力。这个表在写报告的时候非常有用可以让读者快速理解模型的敏感点在哪里。4. 常见问题与排查技巧实录4.1 模型跑出来的曲线不对劲怎么办这是最常见的问题。模型跑出来的曲线要么是直线要么是爆炸式增长要么是负数总之就是不像真实的疫情曲线。我总结了几种典型情况和排查方法。情况一感染人数一直不增长。这通常是因为初始感染人数设得太小或者beta设得太小。排查方法检查初始感染人数是否大于0检查beta是否大于gamma如果beta小于gammaR0小于1疫情会自动消退。情况二感染人数瞬间爆炸一天之内全感染。这通常是因为beta设得太大或者时间步长太大。排查方法减小beta或者减小时间步长。离散模型的时间步长建议不要超过1天。情况三感染人数出现负数。这是数值计算的误差通常发生在感染人数接近0的时候。解决方法在每步更新后把负数截断为0。情况四曲线形状不对没有峰值或者峰值太早太晚。这通常是参数组合的问题。排查方法画相图看S、E、I、R四个变量的关系是否合理。或者手动算几个时间点的值和代码输出对比。实操心得调试模型的时候先用最简单的SIR模型跑通确认基本逻辑没问题再逐步增加复杂度。不要一上来就上最复杂的模型出了问题你都不知道是哪个环节的错。4.2 参数估计不收敛的排查思路参数估计不收敛是另一个高频问题。优化算法跑了几千次迭代损失函数还是很大或者参数跑到了边界上。我遇到过以下几种情况。情况一初始参数猜得太离谱。优化算法对初始值敏感如果初始值离最优解太远可能会陷入局部最优。解决方法用网格搜索先粗筛一遍找到损失函数较小的区域再用优化算法精细搜索。情况二参数不可辨识。不同的参数组合给出相似的拟合效果优化算法在参数空间里来回震荡。解决方法固定一些参数只估计最关键的几个。或者用正则化给参数加先验约束。情况三数据质量太差。如果数据噪声太大或者有系统性偏差优化算法会去拟合噪声而不是真实信号。解决方法对数据做平滑处理或者用稳健的损失函数比如Huber损失。情况四模型结构不对。如果模型本身就不能描述数据的生成过程那再怎么优化参数也没用。解决方法检查模型假设是否合理比如是否应该用SEIR而不是SIR是否需要加无症状仓室。4.3 预测结果和实际不符的原因分析模型预测和实际不符这是常态不是异常。关键是要理解为什么不符而不是简单地否定模型。我总结了几种常见原因。原因一干预措施发生了变化。模型假设beta是常数但实际上干预措施会改变beta。如果实际干预加强了但模型不知道预测就会偏高。解决方法用时变参数模型让beta随时间变化。原因二病毒变异。病毒变异会改变传播率和致病性模型如果用的是旧参数预测就会失准。解决方法定期更新参数或者用多个变异株分别建模。原因三数据报告延迟。实际感染人数总是先于报告确诊人数模型如果拟合的是报告数据预测就会滞后。解决方法在模型里加观测延迟项。原因四人群行为反馈。当疫情严重时人们会自发减少接触这相当于beta自动下降。模型如果没有考虑这个反馈预测就会偏高。解决方法引入行为反馈项让beta随感染人数变化。4.4 常见问题速查表问题现象可能原因排查方法解决方案感染人数不增长初始值太小或beta太小检查初始值和R0增大初始值或beta感染人数爆炸beta太大或步长太大检查beta和步长减小beta或步长出现负数数值误差检查接近0的值截断为0参数不收敛初始值差或不可辨识检查损失函数曲面网格搜索或固定参数预测偏高干预加强或行为反馈对比实际干预措施用时变参数预测偏低干预放松或病毒变异对比病毒特性更新参数拟合很好但预测差过拟合检查参数数量简化模型或加正则化4.5 独家避坑技巧最后分享几个我在项目中踩过的坑以及总结出来的技巧。技巧一保存每次实验的配置和结果。建模是一个反复迭代的过程你会跑很多次实验。如果不记录每次的参数和结果后面根本分不清哪个结果是哪次跑的。我建议用配置文件管理参数每次实验自动保存到一个带时间戳的文件夹里。技巧二可视化要趁早。不要等模型完全跑通了再画图。每改一次代码就画一次图看看曲线形状对不对。可视化是最快的调试手段。技巧三用合成数据验证模型。在真实数据上跑之前先用合成数据验证模型。合成数据的生成过程你完全知道如果模型能从合成数据里恢复出真实参数说明模型实现没问题。技巧四不要迷信单点预测。模型给出的预测是一个点估计但实际的不确定性很大。我建议用集成方法跑多个模型或多个参数组合给出预测区间而不是单一数字。技巧五和领域专家保持沟通。模型是工具不是目的。模型的假设是否合理参数是否可信结果是否有意义这些都需要领域专家来判断。我见过不少建模项目技术做得很漂亮但假设完全脱离实际结果毫无参考价值。这个项目后续还可以这样扩展加入年龄分层因为不同年龄段的接触模式和重症率差异很大加入空间维度因为不同地区的传播条件不同加入疫苗接种模块模拟不同接种策略的效果。每增加一个维度模型的复杂度都会上升但能回答的问题也更精细。关键是找到那个平衡点模型足够复杂以捕捉关键特征又足够简单以可靠标定。