Python实现基于SIR模型的瘟疫扩散地图GIF模拟
本文基于SIR模型，利用Python和NumPy对北欧地区的瘟疫扩散进行空间模拟。通过定义易感、感染及移除状态，结合欧拉法求解微分方程，并将人口密度映射为初始条件。最终使用Matplotlib和PIL库生成展示感染过程的GIF动画，直观呈现病毒在网格地图上的传播动态。
受杰森的《Almost Looks Like Work》启发，本文展示一种病毒传播模型的实现方式。需要明确的是，该模型仅用于演示，并不反映现实世界的疫情数据，请勿将其与西非的真实传染病案例混淆。相反，它更适合作为一种虚构的僵尸爆发场景来进行探讨。接下来，我们将进入具体的模型构建环节。
这里采用SIR模型，其中的字母S、I和R分别代表个体在僵尸疫情中可能处于的三种不同状态。
- S代表易感群体（Susceptible），即那些健康且尚未被感染，但有可能转变为感染者的个体数量。
- I代表染病群体（Infected），即当前已转化为僵尸的个体数量。
- R代表移除群体（Removed），包括因死亡而退出模拟的僵尸，或因治愈而恢复为人类的个体。鉴于僵尸题材中通常不存在治愈机制，我们在此不做过度假设（若将此模型应用于流感传播，则需考虑治愈者因素）。
- 接下来介绍参数β（beta）和γ（gamma）的含义：
- β（beta）表示疾病的传染性强度，在此模型中设定为一旦被咬伤即会立即感染。
- γ（gamma）表示僵尸转化为移除状态（即死亡）的速率，该速率取决于僵尸猎人的平均清除效率。需指出这并非完美模型，请保持耐心理解。
- 公式S′=?βIS描述了健康者转化为僵尸的速率，其中S′表示对时间的导数。
- 公式I′=βIS?γI描述了感染者的变化，包括新增感染量以及僵尸进入移除状态的速率（此处为双关语）。
- 公式R′=γI表示移除群体的增加量，该项在前面的等式中为负值，此处转为正数累加。
上述基础模型未考虑S/I/R的空间分布因素，下面我们将对此进行修正。
一种可行的方法是将瑞典及北欧国家划分为网格，每个网格单元可以感染其相邻单元，具体描述如下：
其中对于特定单元，其相邻的四个单元（上下左右）参与计算（无需考虑对角单元，以免增加不必要的计算复杂度）。
首先进行必要的初始化设置。
import numpy as np import math import matplotlib.pyplot as plt %matplotlib inline from matplotlib import rcParams import matplotlib.image as mpimg rcParams['font.family'] = 'serif' rcParams['font.size'] = 16 rcParams['figure.figsize'] = 12, 8 from PIL import Image
设置合适的beta和gamma数值，即可模拟出大规模感染的效果。
beta = 0.010 gamma = 1
回顾导数的定义：当导数已知且假设时间步长Δt非常小时，经过重新整理，该导数可用于近似预测函数在下一个时间点的取值，此前我们已经定义了u′(t)。
再次进行代码初始化。
import numpy as np import math import matplotlib.pyplot as plt %matplotlib inline from matplotlib import rcParams import matplotlib.image as mpimg rcParams['font.family'] = 'serif' rcParams['font.size'] = 16 rcParams['figure.figsize'] = 12, 8 from PIL import Image
配置适当的beta和gamma值，足以模拟出大半区域被感染的场景。
beta = 0.010 gamma = 1
再次强调导数的应用：在已知导数且Δt极小的前提下，通过重新整理公式，可近似预测函数的下一取值，即我们之前声明过的u′(t)。
这种方法被称为欧拉法（Euler's method），相关代码如下：
def euler_step(u, f, dt): return u + dt * f(u)
我们需要定义函数f(u)。NumPy提供了便捷的数组操作功能。我可能会在后续文章中详细回顾其强大之处，因为需要更多解释，但目前的用法已足以实现预期效果：
def f(u):
S = u[0]
I = u[1]
R = u[2]
new = np.array([-beta*(S[1:-1, 1:-1]*I[1:-1, 1:-1] +
S[0:-2, 1:-1]*I[0:-2, 1:-1] +
S[2:, 1:-1]*I[2:, 1:-1] +
S[1:-1, 0:-2]*I[1:-1, 0:-2] +
S[1:-1, 2:]*I[1:-1, 2:]),
beta*(S[1:-1, 1:-1]*I[1:-1, 1:-1] +
S[0:-2, 1:-1]*I[0:-2, 1:-1] +
S[2:, 1:-1]*I[2:, 1:-1] +
S[1:-1, 0:-2]*I[1:-1, 0:-2] +
S[1:-1, 2:]*I[1:-1, 2:]) - gamma*I[1:-1, 1:-1],
gamma*I[1:-1, 1:-1]
])
padding = np.zeros_like(u)
padding[:,1:-1,1:-1] = new
padding[0][padding[0] < 0] = 0
padding[0][padding[0] > 255] = 255
padding[1][padding[1] < 0] = 0
padding[1][padding[1] > 255] = 255
padding[2][padding[2] < 0] = 0
padding[2][padding[2] > 255] = 255
return padding
导入北欧国家的人口密度地图，并进行下采样处理，以加快计算速度。
from PIL import Image
img = Image.open('popdens2.png')
img = img.resize((img.size[0]/2,img.size[1]/2))
img = 255 - np.asarray(img)
imgplot = plt.imshow(img)
imgplot.set_interpolation('nearest')
北欧国家的人口密度地图（注：图中未包含丹麦）。
S矩阵代表易感个体，其初始值应近似于人口密度分布。感染者初始值设为0，并将斯德哥尔摩设定为第一个感染源。
S_0 = img[:,:,1] I_0 = np.zeros_like(S_0) I_0[309,170] = 1 # patient zero
由于初始阶段尚未有人死亡，因此将R矩阵也初始化为0。
R_0 = np.zeros_like(S_0)
接着初始化模拟的总时长及其他相关参数。
T = 900 # final time dt = 1 # time increment N = int(T/dt) + 1 # number of time-steps t = np.linspace(0.0, T, N) # time discretization # initialize the array containing the solution for each time-step u = np.empty((N, 3, S_0.shape[0], S_0.shape[1])) u[0][0] = S_0 u[0][1] = I_0 u[0][2] = R_0
我们需要自定义一个颜色映射表，以便将感染矩阵可视化地显示在地图上。
import matplotlib.cm as cm
theCM = cm.get_cmap("Reds")
theCM._init()
alphas = np.abs(np.linspace(0, 1, theCM.N))
theCM._lut[:-3,-1] = alphas
现在，让我们坐下来欣赏模拟结果……
for n in range(N-1): u[n+1] = euler_step(u[n], f, dt)
最后进行图像渲染处理，将其制作成GIF格式，毕竟每个人都喜欢GIF动画！
from images2gif import writeGif
keyFrames = []
frames = 60.0
for i in range(0, N-1, int(N/frames)):
imgplot = plt.imshow(img, vmin=0, vmax=255)
imgplot.set_interpolation("nearest")
imgplot = plt.imshow(u[i][1], vmin=0, cmap=theCM)
imgplot.set_interpolation("nearest")
filename = "outbreak" + str(i) + ".png"
plt.savefig(filename)
keyFrames.append(filename)
images = [Image.open(fn) for fn in keyFrames]
gifFilename = "outbreak.gif"
writeGif(gifFilename, images, duration=0.3)
plt.clf()
本站声明
本文内容仅供参考，不代表本站观点。转载或引用请注明出处。如发现内容问题，请通过联系我们页面反馈。
