1. 从零理解SEIR模型它如何刻画新冠疫情的传播如果你刚开始接触数学建模或者想用Python做一些有意思的仿真新冠疫情SEIR模型绝对是一个经典且富有启发性的起点。它不像一些复杂的模型那样让人望而生畏但又比简单的SIR模型更贴近现实能让你直观地看到隔离、防控措施是如何在数学上影响疫情发展的。简单来说SEIR模型把人群分成了四类易感者S、潜伏者E、感染者I和康复者R。这个“E”的加入就是它比SIR模型更精细的地方因为它考虑到了病毒有潜伏期——一个人被感染后不会立刻具有传染性而是要经过一段时间。想象一下一个封闭的社区。最开始几乎所有人都是易感者S只有极少数的感染者I作为“种子”。当一个易感者接触到一个感染者他有一定的概率被感染但不会马上变成感染者而是先进入潜伏期成为潜伏者E。潜伏者自己还没有症状也不传染别人但身体里的病毒正在悄悄复制。过了一段时间这个时间就是潜伏期的倒数潜伏者才会“转正”成为具有传染性的感染者I。感染者会持续传播病毒直到他康复或死亡进入康复者R阵营并且假设康复后获得了免疫力不会再被感染。整个传播的动态就用一组微分方程来描述这组方程定义了这四类人如何随时间此消彼长。那么这个模型能告诉我们什么最核心的就是预测疫情高峰何时到来、规模有多大以及评估不同防控策略的效果。比如通过调整模型中的“接触率”参数你可以模拟“保持社交距离”或“封控”的效果通过调整“潜伏期转感染率”和“感染者的移出率”包括康复和死亡你可以研究加快检测速度让潜伏者更快被隔离或者提高医疗救治水平的影响。对于Python小白来说实现这个模型并看到曲线随着参数调整而变化是理解微分方程和数值计算威力的绝佳方式。接下来我会手把手带你从原理到代码完整复现这个模型并分析几个关键场景。2. SEIR模型微分方程组的拆解与参数意义要真正玩转SEIR模型不能只停留在“有四类人”的概念上必须深入理解驱动这个系统的微分方程组。每一个方程都代表了一类人群数量的变化率而方程中的参数就是撬动整个疫情发展的“杠杆”。我们假设总人口N是常数不考虑出生和自然死亡即 S E I R N。2.1 核心微分方程组下面是经典的SEIR模型微分方程组我会逐一解释每一项的含义易感者S的变化率 dS/dt -β * S * I / NdS/dt表示易感者数量随时间的变化速度。因为易感者只会减少被感染所以这个值通常是负的。-β * S * I / N这是减少的部分。β是有效接触率表示一个感染者单位时间内能传染的人数。S/N是易感者在总人口中的比例。所以β * I可以理解为单位时间内一个感染者造成的“传染力”这个力作用在易感者比例S/N上再乘以总易感者人数S就得到了单位时间内新感染的人数。负号表示S在减少。潜伏者E的变化率 dE/dt β * S * I / N - σ * Eβ * S * I / N这部分和上面一样是从易感者“流入”潜伏者的新感染人数。- σ * E这是从潜伏者“流出”的人数。σ是潜伏期转感染率。它的倒数1/σ就是平均潜伏期比如如果σ0.2平均潜伏期就是5天。σ * E表示单位时间内有多少潜伏者结束了潜伏期变成了感染者。感染者I的变化率 dI/dt σ * E - γ * Iσ * E从潜伏者“流入”感染者的人数即上一方程流出的部分。- γ * I从感染者“流出”的人数。γ是康复率或移出率。它的倒数1/γ就是平均感染期从具有传染性到康复或死亡的时长。γ * I表示单位时间内有多少感染者康复或死亡了。康复者R的变化率 dR/dt γ * Iγ * I这部分就是从感染者康复过来的人数。康复者数量只增不减。注意这里有一个关键的简化假设即康复者获得永久免疫力不会再变回易感者。对于新冠病毒的某些变体和时间段这可能不完全准确但作为基础模型是合理的。2.2 关键参数的现实意义与取值理解参数是调整模型、拟合现实数据的基础。下面这个表格总结了三个核心参数参数符号参数名称现实意义典型取值参考以新冠原始毒株为例如何影响疫情曲线β有效接触率衡量病毒传播能力综合了人际接触频率和单次接触传播概率。0.5 - 0.8 /天最关键的杠杆。β越大传播越快峰值越高、越早出现。降低β相当于实施社交疏离、戴口罩、减少聚集。σ潜伏期转感染率潜伏期长短的倒数。σ越大潜伏期越短。1/5.2 ≈ 0.192 /天 (潜伏期约5.2天)影响疫情发展的“节奏”。潜伏期越短σ越大感染者更早出现疫情发展更迅速。γ康复率移出率感染期长短的倒数。γ越大感染期越短。1/7 ≈ 0.143 /天 (感染期约7天)影响感染者的存量。感染期越短γ越大感染者被移出的速度越快能有效压低峰值、缩短流行时间。基本再生数R0这是一个衍生出来的、极其重要的流行病学指标。它表示在完全易感人群中一个感染者在其整个传染期内平均能传染多少人。在SEIR模型中R0 β / γ。如果R0 1疫情会扩散R0 1疫情会逐渐消失。防控措施的本质就是通过降低β减少接触或提高γ加快隔离治疗使得有效再生数Rt随时间变化的R0小于1。2.3 为什么是微分方程你可能会问为什么非要用微分方程因为疫情传播是一个连续变化的过程。每一刻都有人从S变成E从E变成I从I变成R。微分方程完美地刻画了这种连续、动态的变化关系。而我们用Python要做的就是通过数值方法比如欧拉法或更常用的龙格-库塔法来“求解”这组方程即计算出从初始时刻开始未来每一天S、E、I、R各有多少人。这个过程叫做“数值仿真”。3. 手把手Python实现从方程到动态曲线理论说得再多不如一行代码。我们使用Python中最常用的科学计算库NumPy和绘图库Matplotlib来实现SEIR模型的求解与可视化。这里我会选择scipy库中的odeint函数来求解微分方程它比手写欧拉法更稳定、更精确。3.1 环境准备与库安装首先确保你的Python环境已经安装了必要的库。打开你的命令行CMD、Terminal或PyCharm的终端输入以下命令进行安装pip install numpy matplotlib scipy如果你使用的是Anaconda这些库通常已经预装了。为了清晰我们将代码写在一个完整的脚本里。3.2 核心代码实现与逐行解读创建一个新的Python文件比如seir_model.py然后输入以下代码。我会在代码块中加入大量注释确保你能看懂每一行的作用。# seir_model.py import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt # 1. 定义SEIR模型的微分方程组 def seir_model(y, t, beta, sigma, gamma, N): 定义SEIR模型的微分方程。 参数: y: 包含[S, E, I, R]当前值的列表 t: 当前时间odeint内部使用我们传入时间序列 beta, sigma, gamma, N: 模型参数 返回: dydt: [dS/dt, dE/dt, dI/dt, dR/dt]的列表 S, E, I, R y # 解包当前状态 # 易感者变化率 dS_dt -beta * S * I / N # 潜伏者变化率 dE_dt beta * S * I / N - sigma * E # 感染者变化率 dI_dt sigma * E - gamma * I # 康复者变化率 dR_dt gamma * I return [dS_dt, dE_dt, dI_dt, dR_dt] # 2. 设置模型参数和初始条件 # 总人口 N 1000000 # 假设一个100万人口的城市 # 初始条件假设最初有10个感染者有100个潜伏者可能来自外部输入其余都是易感者 I0 10 E0 100 R0 0 S0 N - I0 - E0 - R0 # 将初始条件组合成一个列表 initial_conditions [S0, E0, I0, R0] # 流行病学参数参考新冠原始毒株的典型值 beta 0.6 # 有效接触率约合R04.2因为gamma0.143 sigma 1/5.2 # 潜伏期转感染率潜伏期约5.2天 gamma 1/7 # 康复率感染期约7天 # 计算基本再生数R0 R0_basic beta / gamma print(f基本再生数 R0 {R0_basic:.2f}) # 3. 定义时间点模拟200天 days 200 t np.linspace(0, days, days1) # 从第0天到第200天共201个点 # 4. 求解微分方程组 solution odeint(seir_model, initial_conditions, t, args(beta, sigma, gamma, N)) # solution是一个数组每一行对应一个时间点每一列对应[S, E, I, R] S, E, I, R solution.T # 转置后分别赋值给四个变量 # 5. 可视化结果 plt.figure(figsize(12, 8)) plt.plot(t, S, b, label易感者(S), linewidth2) plt.plot(t, E, orange, label潜伏者(E), linewidth2) plt.plot(t, I, r, label感染者(I), linewidth2) plt.plot(t, R, g, label康复者(R), linewidth2) plt.xlabel(时间 (天), fontsize14) plt.ylabel(人数, fontsize14) plt.title(fSEIR模型仿真 (N{N}, β{beta}, σ{sigma:.3f}, γ{gamma:.3f}, R0{R0_basic:.2f}), fontsize16) plt.legend(locbest, fontsize12) plt.grid(True, linestyle--, alpha0.7) # 标记感染者峰值 peak_day t[np.argmax(I)] # argmax找到I最大值的索引 peak_infected np.max(I) plt.annotate(f峰值: {int(peak_infected)}人\n第{int(peak_day)}天, xy(peak_day, peak_infected), xytext(peak_day10, peak_infected*0.8), arrowpropsdict(facecolorblack, shrink0.05, width1.5, headwidth8), fontsize12) plt.tight_layout() plt.show() # 6. 输出一些关键数据 print(f模拟结束第{days}天时) print(f 易感者(S): {int(S[-1])} 人) print(f 潜伏者(E): {int(E[-1])} 人) print(f 感染者(I): {int(I[-1])} 人) print(f 康复者(R): {int(R[-1])} 人) print(f 最终感染比例: {(R[-1]/N*100):.2f}%)运行这段代码你会看到一张绘制了四条曲线的图分别代表了S、E、I、R四类人群在200天内的变化。红色曲线感染者I会先上升后下降形成一个明显的“波峰”这就是疫情高峰。橙色曲线潜伏者E的峰值通常比感染者早几天出现。3.3 代码运行要点与可能遇到的问题关于odeint它是scipy中一个强大的常微分方程求解器。我们只需要定义方程seir_model和初始条件它就会帮我们计算出所有时间点上的解。args参数用于传递方程中除了y和t之外的其他参数beta, sigma等。如果遇到安装错误确保你的pip版本是最新的pip install --upgrade pip或者如果你在使用国内网络可以尝试使用清华镜像源加速安装pip install numpy matplotlib scipy -i https://pypi.tuna.tsinghua.edu.cn/simple。图形不显示如果你在命令行运行且没有弹出图形窗口可能是因为没有图形后端。可以尝试在代码开头加上import matplotlib; matplotlib.use(TkAgg)或者确保你安装了python-tkLinux或使用Jupyter Notebook环境运行。理解输出观察曲线你会发现易感者蓝色从高位逐渐下降最终趋于一个稳定值即始终未被感染的人。康复者绿色单调上升最终趋于稳定。这模拟了一次完整的疫情流行过程。4. 模型应用模拟不同防控策略的效果现在我们已经有了一个可以运行的“疫情沙盘”接下来最有趣的部分就是用它来模拟不同的公共卫生干预措施看看数学上如何预测它们的效果。我们主要通过调整参数β有效接触率来模拟。4.1 场景一无干预基线情况这就是我们上面代码模拟的情况β0.6R0≈4.2。运行后你会发现疫情发展迅猛感染峰值很高最终大部分人可能超过80%都会被感染。这对应了早期没有采取任何措施的情形。4.2 场景二实施中等强度社交疏离假设通过宣传戴口罩、减少聚会将人群的接触率降低了30%。那么新的β 0.6 * (1 - 0.3) 0.42。我们修改参数重新运行求解和绘图部分可以复制一份代码或写一个循环来比较。# 在原有代码基础上增加一个场景对比 beta_intervention 0.6 * 0.7 # 降低30%接触率 solution_int odeint(seir_model, initial_conditions, t, args(beta_intervention, sigma, gamma, N)) S_int, E_int, I_int, R_int solution_int.T # 将两个场景的感染者曲线画在一起对比 plt.figure(figsize(12, 6)) plt.plot(t, I, r-, labelf无干预 (β{beta}, R0{beta/gamma:.2f}), linewidth2) plt.plot(t, I_int, b--, labelf社交疏离 (β{beta_intervention}, R0{beta_intervention/gamma:.2f}), linewidth2) plt.xlabel(时间 (天)) plt.ylabel(感染者人数 (I)) plt.title(不同防控强度下感染者数量对比) plt.legend() plt.grid(True, linestyle--, alpha0.7) plt.tight_layout() plt.show()你会看到蓝色虚线干预后的峰值比红色实线无干预低很多而且高峰到来的时间也推迟了。这就是“压平曲线”的直观数学体现目的是避免医疗资源被瞬间击穿。4.3 场景三早期严格封控与后期常态化管理现实情况更复杂。我们可能先采取一段时间的严格封控大幅降低β然后再放松到常态化管理β小幅回升。这需要让β成为一个随时间变化的函数beta(t)。# 定义随时间变化的beta函数 def beta_function(t): 分段函数模拟防控策略变化 第0-30天严格封控β很低 (0.15) 第31天以后常态化管理β回升但低于基线 (0.4) if t 30: return 0.15 else: return 0.4 # 修改微分方程函数使其能接受随时间变化的beta def seir_model_variable_beta(y, t, sigma, gamma, N): S, E, I, R y beta_t beta_function(t) # 获取当前时间的beta值 dS_dt -beta_t * S * I / N dE_dt beta_t * S * I / N - sigma * E dI_dt sigma * E - gamma * I dR_dt gamma * I return [dS_dt, dE_dt, dI_dt, dR_dt] # 求解 solution_var odeint(seir_model_variable_beta, initial_conditions, t, args(sigma, gamma, N)) S_var, E_var, I_var, R_var solution_var.T # 绘图 plt.figure(figsize(12, 6)) plt.plot(t, I_var, purple, label动态防控策略, linewidth2) plt.axvline(x30, colorgray, linestyle:, linewidth1.5, label封控结束) plt.xlabel(时间 (天)) plt.ylabel(感染者人数 (I)) plt.title(动态防控策略下的疫情发展) plt.legend() plt.grid(True, linestyle--, alpha0.7) plt.tight_layout() plt.show()这个模拟显示早期严格措施迅速压制了疫情曲线几乎没起来但放松后由于仍有大量易感者疫情会出现一个“反弹峰”不过峰值比完全不干预要低得多。这解释了为什么防控需要持续和灵活调整。4.4 场景四提高检测与隔离效率影响σ和γ除了降低接触率β加快病例发现和隔离也能影响模型。这体现在两个方面缩短潜伏期到被隔离的时间这相当于增大了潜伏者转出率σ。原本潜伏者自然转成感染者现在可能提前被检测出并隔离从而不再进入I类具有传染性。在模型中我们可以增加一个“隔离者”类别但更简单的处理是增大σ并假设被“提前转出”的人直接进入R类实际上是被移除了传播链。缩短感染期通过有效治疗和隔离减少感染者传播病毒的时间相当于增大了康复率γ。你可以尝试同时增大σ和γ例如都增加50%观察曲线变化。你会发现疫情峰值更低流行期更短。实操心得在调整参数做模拟时一次最好只改变一个参数并观察曲线的变化这样才能清晰理解每个参数的“杠杆效应”。同时参数的设置要有现实依据例如将感染期缩短到1天γ1是不太现实的。多尝试多对比是学习建模的最好方法。5. 模型局限性与进阶思考SEIR模型是一个强大的教学和初步分析工具但它是对复杂现实的极度简化。认识到它的局限性才能更好地使用它并知道何时需要更复杂的模型。5.1 SEIR模型的主要简化假设均匀混合假设模型假设人群是完全均匀混合的即任何一个易感者接触任何一个感染者的概率相同。这显然忽略了年龄结构、社交网络、地域差异和超级传播事件。常数参数假设β、σ、γ被假设为常数。但实际上它们会随时间变化如季节影响、病毒变异、防控措施强度变化、医疗资源挤兑等。封闭系统假设总人口N不变不考虑出生、死亡非疫情原因、迁入和迁出。对于短期疫情模拟可以接受对于长期或全球模型则需要调整。完全免疫假设康复者获得永久免疫不会再次感染。这与新冠病毒存在再感染和免疫力衰减的情况不符。忽略无症状感染经典的SEIR模型没有区分有症状和无症状感染者而无症状感染者同样具有传染性但其传播能力和被发现隔离的概率与有症状者不同。5.2 如何扩展模型以贴近现实针对以上局限建模者发展出了许多扩展模型SEIRS模型在SEIR基础上考虑康复者免疫力衰减重新变为易感者S。这需要增加一个从R流回S的项。考虑年龄分层的模型将人群按年龄分组如儿童、成人、老人每组设置不同的接触矩阵谁和谁接触和疾病参数如重症率这能更准确地评估疫苗优先接种策略。仓室细分模型将感染者I细分为无症状感染者A、轻症I_m、重症I_s等将康复者细分为康复R和死亡D。这能模拟医疗资源需求。网络模型放弃均匀混合假设用图网络来模拟个体之间的接触关系可以研究超级传播和社区结构的影响。但这需要更复杂的数据和计算。时变参数与随机模型将参数设为随时间变化的函数甚至引入随机噪声以模拟防控措施的动态调整和现实中的不确定性。对于Python进阶学习者尝试实现SEIRS模型是一个很好的挑战。你只需要在原有方程组中为康复者R的变化率增加一项 ξ * R表示免疫力衰减重新变为易感者同时在易感者S的变化率中减去这一项或加上取决于符号约定。然后调整参数ξ免疫力丧失率观察疫情是否会从一次性爆发转变为周期性流行。5.3 从仿真到预测参数估计与数据拟合我们上面的模拟都是“拍脑袋”给定参数。真正的建模挑战在于如何利用真实的疫情数据每日新增确诊、累计康复、死亡等来反推模型中最可能的参数值β, σ, γ这个过程叫做参数估计或模型拟合。常用方法是最小二乘法。思路是调整模型参数让模型输出的曲线比如每日新增感染者σ*E与真实数据曲线之间的差距误差平方和最小。在Python中可以使用scipy.optimize库中的curve_fit或least_squares函数来实现。这是一个更高级的话题需要处理数据噪声、初始条件不确定等问题但它是连接理论模型与现实世界的桥梁。我个人在尝试用早期武汉数据拟合SEIR模型时发现直接拟合往往效果不佳因为真实的防控措施是剧烈变化的导致β不是一个常数。更好的做法是采用分段常数的β或者使用更复杂的时变函数。这提醒我们建模时对现实背景的理解往往比数学技巧本身更重要。6. 给Python小白的项目拓展建议如果你已经成功运行了基础SEIR模型并做了一些参数调整实验那么你已经入门了。为了深化学习我建议你尝试以下几个小项目每一个都能帮你巩固不同的Python和建模技能可视化升级不用matplotlib尝试使用plotly库来制作交互式图表。你可以添加滑块plotly.widgets来动态调整β、σ、γ参数并实时看到曲线变化。这种交互体验能极大地加深你对参数影响的理解。构建SEIRS模型如前所述实现免疫力衰减的SEIRS模型。观察当免疫力持续时间比如1年小于病毒传播周期时疫情如何从爆发转为地方性流行Endemic甚至出现规律的波峰。模拟“疫苗”接种在模型中引入疫苗接种项。假设每天有一定比例的易感者S被接种疫苗并直接转移到康复者R或一个单独的“已接种者V”类别。研究不同的疫苗接种速度对遏制疫情高峰和最终感染规模的影响。简单GUI应用使用tkinter或PyQt如果觉得难可以用gradio快速构建做一个简单的桌面应用。提供几个输入框让用户输入β、σ、γ、初始感染者等参数点击按钮后弹出模拟结果图表。这将让你把模型打包成一个可分享的小工具。尝试其他数值解法不使用odeint自己用欧拉法Euler‘s Method来实现微分方程的求解。欧拉法虽然精度不如龙格-库塔法但原理非常简单S(t1) S(t) dS/dt * dt。自己实现一遍你会对“数值求解”有刻骨铭心的理解。学习数学建模和Python最好的方法就是“做”。从复制代码到修改参数再到增加新功能最后尝试解决自己的问题。SEIR模型就像一把钥匙帮你打开了用计算思维理解复杂世界的大门。当你看到自己写的几行代码能模拟出一场疫情的兴衰并能通过调整参数来评估不同策略时那种成就感就是持续学习的最佳动力。