公司网站维护建设的通知,网站建设 智能建站,站长工具seo综合查询推广,莆田有建设网站的公司码在当今世界#xff0c;传染病的传播模拟已成为公共卫生领域的重要研究工具。传统的传染病模型往往忽略了空间因素#xff0c;将整个群体视为一个均质的整体。然而#xff0c;现实中的疾病传播与地理空间位置密切相关#xff0c;相邻区域之间的传播风险远高于远距离区域。今…在当今世界传染病的传播模拟已成为公共卫生领域的重要研究工具。传统的传染病模型往往忽略了空间因素将整个群体视为一个均质的整体。然而现实中的疾病传播与地理空间位置密切相关相邻区域之间的传播风险远高于远距离区域。今天我们将通过一个结合地理信息系统 (GIS) 和随机模拟算法的案例深入探讨如何在空间维度上模拟传染病的传播过程。
引言当传染病模型遇上地理空间
传染病的传播从来都不是抽象的数字游戏它实实在在地发生在我们生活的城市空间中。想象一下在伦敦这样的大都市中不同区域之间的人口流动、空间邻近性如何影响疾病的传播速度和范围这正是我们今天要探索的主题。
我们将使用经典的 SIR 传染病模型但会为其添加 地理空间 这一重要维度。SIR 模型将人群分为三类易感者 (Susceptible)、感染者 (Infectious) 和康复者 (Recovered)。通过引入空间权重矩阵我们可以量化不同区域之间的邻近关系使模型能够更真实地反映疾病在地理空间上的传播过程。
数据准备加载伦敦的地理空间数据
首先我们需要准备地理空间数据。代码中加载了伦敦市 (City of London) 和卡姆登区 (Camden) 的低级统计区 (LSOA) 边界数据。这些数据是英国国家统计数据的一部分能够帮助我们精确界定不同的地理单元。
# 加载伦敦市和卡姆登区的shapefile
lon gpd.read_file(LSOA_2011_BFE_City_of_London.shp)
cam gpd.read_file(LSOA_2011_BFE_Camden.shp)
# 合并两个GeoDataFrame
df pd.concat([lon, cam], ignore_indexTrue)
地理数据通常有不同的坐标参考系统 (CRS)。为了确保后续计算的准确性我们需要进行坐标转换。我们先保留一个 WGS84 (经纬度) 坐标系的副本用于地图绘制然后将数据转换为等距坐标系 (EPSG:3857) 用于质心计算和距离测量。
# 保留WGS84副本用于绘制经纬度地图
df_ll df.to_crs(epsg4326)# 在进行质心计算前重新投影到等距CRS(EPSG:3857)
df df_ll.to_crs(epsg3857)
坐标系统的转换是地理空间分析中非常关键的一步。想象一下如果我们在经纬度坐标系中直接计算两点之间的距离由于地球曲率的影响结果会有很大误差。而转换到等距坐标系后我们就可以像在平面上一样准确计算距离和面积了 。 接下来我们计算每个区域的质心并获取它们在两种坐标系中的坐标
# 计算两种CRS下的质心
centroids_3857 df.geometry.centroid
coords_3857 np.vstack([centroids_3857.x, centroids_3857.y]).T
centroids_ll df_ll.geometry.centroid
coords_ll np.vstack([centroids_ll.x, centroids_ll.y]).T
空间权重矩阵定义区域间的邻近关系
在空间分析中邻近性 是一个关键概念。如何定义两个区域是否相邻在代码中我们使用了 Queen contiguity 规则这是空间统计中常用的一种邻近性定义。Queen contiguity 认为如果两个区域共享一条边或一个角则它们是相邻的。
# 使用Queen邻接构建空间权重
w libpysal.weights.Queen.from_dataframe(df, use_indexTrue)# 转换为权重矩阵
n len(df)
W np.zeros((n, n), dtypeint)
for i, neighbors in w.neighbors.items():for j in neighbors:W[i, j] 1W[j, i] 1 # 对称矩阵
这种邻接关系定义很有意思它模拟了一种 棋盘式 的邻近性。想象一下国际象棋中的皇后可以沿对角线移动Queen contiguity 就是用这种方式定义邻近性的 chess。与只考虑共享边的 Rook contiguity 相比它考虑了更多的邻近关系。
为了更直观地理解我们构建的权重矩阵我们可以绘制权重分布的直方图
# 正权重的直方图
dist W[W 0]
plt.figure()
plt.hist(dist, bins20)
plt.title(Distribution of Weights 0)
plt.show()
由于我们使用的是二元邻接矩阵 (相邻为 1不相邻为 0)这个直方图应该会呈现出明显的双峰分布其中大部分权重为 0而相邻区域的权重为 1。不过在后续步骤中我们还会引入阈值处理让邻接关系变得更加 模糊。
构建空间网络阈值化邻接关系
有时候我们可能需要考虑更灵活的邻近关系定义。例如不仅考虑直接相邻的区域还可以根据距离或其他因素定义一个 影响范围。在代码中我们通过设置一个阈值将权重矩阵转换为一个二元邻接矩阵
# 阈值邻接
threshold 0.16
adj (W threshold).astype(int)
这里的阈值设置是一个关键参数它决定了哪些区域会被视为 相邻。阈值越高邻接关系越严格网络中的连接越少阈值越低邻接关系越宽松网络中的连接越多。这个参数的设置需要根据具体问题和数据特点来调整 。
我们可以将这些邻接关系可视化看看在地理空间上区域之间是如何连接的
# 绘制高于阈值的边
g_df df_ll.copy()
plt.figure(figsize(6, 6))
g_df.plot(facecolornone, edgecolorgray)
plt.scatter(coords_ll[:, 0], coords_ll[:, 1], s10, colorblack)
for i in range(n):for j in range(i 1, n):if adj[i, j]:plt.plot([coords_ll[i, 0], coords_ll[j, 0]],[coords_ll[i, 1], coords_ll[j, 1]],colorblue, linewidth0.5)
plt.title(Network Edges (thresholded))
plt.show() 看着这些蓝色的连接线在地图上延伸我们仿佛看到了疾病可能传播的 高速公路。这些连接线构成了一个空间网络疾病将在这个网络上进行传播 。
SIR 模型从经典到空间扩展
现在我们来介绍 SIR 模型的基本原理。SIR 模型是一个经典的传染病模型它将人群分为三个状态
S (Susceptible): 易感者尚未感染但可能被感染的人I (Infectious): 感染者具有传染性的人R (Recovered): 康复者已经康复并获得免疫力的人
模型的核心是两个速率参数
β(beta): 感染率表示易感者与感染者接触后被感染的概率γ(gamma): 康复率表示感染者康复的概率
在传统的 SIR 模型中通常假设人群是完全混合的即每个人与其他人接触的概率是相等的。但在我们的空间扩展版本中这个假设被打破了。我们考虑了空间邻近性即一个区域的感染状态会受到其相邻区域感染状态的影响。
在代码中我们设置了感染率和康复率参数
# SIR参数
beta 3.0 # 感染率
gamma 1.0 # 康复率
这些参数的值会显著影响疾病传播的动态。例如较高的感染率会导致疾病更快地传播而较高的康复率会缩短感染期减缓疾病的传播速度 ⚙️。
初始状态设置疫情的起点
每个传染病模拟都需要一个起点。在代码中我们随机选择一个区域作为初始感染区域
# 初始化状态: S, I, R
status np.array([S] * n)
i0 random.randrange(n)
status[i0] I
这个初始感染点的选择虽然是随机的但在现实中疫情的起点可能与人口密度、交通枢纽等因素有关。如果我们有更详细的数据我们可以根据实际情况选择更合理的初始感染点 。
Gillespie 随机模拟算法捕捉传染病的随机性
传染病的传播本质上是一个随机过程因为每个接触事件是否导致感染都是概率性的。为了准确捕捉这种随机性我们使用了 Gillespie 随机模拟算法这是一种在分子尺度上模拟随机过程的强大方法。
Gillespie 算法的核心思想是
计算系统中所有可能事件的发生速率根据这些速率计算每个事件发生的概率随机选择一个事件发生计算到下一个事件发生的时间间隔重复上述过程
在我们的 SIR 模型中有两种可能的事件
感染事件易感者被感染者感染康复事件感染者康复
让我们看看代码中是如何实现这个算法的
# 存储模拟结果
T_max 10.0 # 最大模拟时间
times [0.0] # 时间点列表
states [status.copy()] # 每个时间点的状态列表
adj_list {i: list(np.where(adj[i] 1)[0]) for i in range(n)} # 邻接列表
current_time 0.0 # 当前时间# Gillespie模拟
i 0
while current_time T_max:# 找出所有感染者I_indices np.where(status I)[0]# 计算每个感染者的康复速率rec_rates {idx: gamma for idx in I_indices}# 找出所有易感者S_indices np.where(status S)[0]# 计算每个易感者的感染速率与相邻感染者数量有关inf_rates {idx: beta * sum(status[j] I for j in adj_list[idx]) for idx in S_indices}# 合并所有事件及其速率all_rates {**rec_rates, **inf_rates}total_rate sum(all_rates.values())# 如果没有事件发生结束模拟if total_rate 0:break# 提取事件节点和速率nodes, rates zip(*all_rates.items())# 计算每个事件发生的概率probs np.array(rates) / total_rate# 随机选择一个事件event_node np.random.choice(nodes, pprobs)# 更新状态if status[event_node] S:status[event_node] I # 易感者被感染else:status[event_node] R # 感染者康复# 计算到下一个事件的时间间隔dt np.random.exponential(1.0 / total_rate)current_time dt# 存储当前时间和状态times.append(current_time)states.append(status.copy())i 1
Gillespie 算法的一个重要特点是它精确地模拟了随机过程的时间演化而不是像欧拉法那样使用固定的时间步长。这使得它在处理具有不同时间尺度的事件时特别有效 ⏳。
在我们的实现中每个易感者的感染速率不仅取决于感染率 β还取决于其相邻区域中有多少个感染区域。这正是空间因素被引入模型的关键所在 。
结果分析传染病的时间动态
模拟结束后我们可以分析传染病在时间上的传播动态。首先我们创建一个时间序列数据框记录每个时间点上易感者、感染者和康复者的数量
# 时间序列DataFrame
df_counts [{time: t,S: int((st S).sum()),I: int((st I).sum()),R: int((st R).sum())}for t, st in zip(times, states)]
df_ts pd.DataFrame(df_counts)
然后我们可以绘制这三类人群随时间变化的曲线
plt.figure()
for compartment in [S, I, R]:plt.plot(df_ts[time], df_ts[compartment], labelcompartment)
plt.xlabel(Time)
plt.ylabel(Count)
plt.legend()
plt.title(SIR Dynamics Over Time)
plt.show() 典型的 SIR 模型时间曲线应该呈现出这样的模式
易感者数量 (S) 随时间单调递减感染者数量 (I) 先增加后减少形成一个峰值康复者数量 (R) 随时间单调递增
看着这些曲线的变化我们可以直观地理解传染病的传播过程从少数人感染开始逐渐扩散到高峰然后随着越来越多的人康复感染人数逐渐减少 。
空间传播动画疾病在地理空间上的扩散
最令人兴奋的部分莫过于将疾病的空间传播过程可视化。通过创建动画我们可以直观地看到疾病如何从初始感染点开始沿着空间网络向周围区域扩散。
首先我们准备动画所需的数据
# 准备DataFrame用于动画
records []
for step, st in enumerate(states):for idx, s in enumerate(st):records.append({time_step: step,node: idx,status: s,x: coords_ll[idx, 0],y: coords_ll[idx, 1]})
anim_df pd.DataFrame(records)
然后我们实现动画生成函数
# 使用PillowWriter生成动画
def animate_sir_map(output_filelondon.gif, fps5):fig, ax plt.subplots(figsize(6, 6))def update(frame):ax.clear()df_ll.plot(axax, facecolorlightgray, edgecolorwhite)df_step anim_df[anim_df[time_step] frame]for label, color in zip([S, I, R], [blue, red, green]):sub df_step[df_step[status] label]ax.scatter(sub[x], sub[y], labellabel, s20)ax.set_title(fTime step: {frame})ax.legend(locupper right, fontsizesmall)ax.set_axis_off()writer PillowWriter(fpsfps)anim FuncAnimation(fig, update, frameslen(states), interval200)anim.save(output_file, writerwriter)print(fAnimation saved as {output_file} using PillowWriter)# 生成动画
animate_sir_map() 当我们播放这个动画时可以看到红色的感染区域如何从一个点开始逐渐扩散到周围的蓝色易感区域然后变成绿色的康复区域。这种时空动态直观地展示了空间邻近性如何影响疾病的传播模式 →→。
模型的局限性与改进方向
尽管我们的空间 SIR 模型比传统模型更贴近现实但它仍然有一些局限性 空间权重的简化我们使用了简单的二元邻接矩阵来表示空间关系而现实中的疾病传播可能受到多种因素影响如人口流动、交通网络等。 同质性假设模型假设每个区域的人口特征是相同的而实际上不同区域的人口密度、年龄结构等可能差异很大。 确定性参数我们使用了固定的感染率和康复率而在现实中这些参数可能随时间和空间变化。
针对这些局限性我们可以考虑以下改进方向 更复杂的空间权重可以基于实际交通数据或人口流动数据构建更真实的空间权重矩阵。 异质性建模可以为不同区域设置不同的人口特征和疾病传播参数。 引入行为干预可以模拟社交距离、封锁等干预措施对疾病传播的影响 ️。 考虑人口流动可以引入更复杂的人口流动模型模拟通勤等日常活动对疾病传播的影响。
结语空间视角下的传染病防控
通过这个案例我们看到了地理空间视角在传染病建模中的重要性。空间邻近性、区域连接性等因素显著影响着疾病的传播模式和速度。在公共卫生政策制定中考虑这些空间因素可以帮助我们更精准地分配资源、实施干预措施。
想象一下当我们能够在地图上实时看到疾病的传播路径和高风险区域时公共卫生部门就可以更有针对性地部署医疗资源、实施隔离措施从而最大限度地减少疾病的影响 。
这个模型不仅是一个学术研究工具也为我们理解现实世界中的传染病传播提供了一个窗口。随着地理信息技术和计算能力的不断发展我们有理由相信未来的传染病模型将更加精确、更加贴近现实为全球公共卫生安全提供更强有力的支持 。
如果你对这个模型感兴趣不妨尝试调整参数看看会发生什么变化增加或减少感染率 β看看疾病传播速度如何变化调整空间权重的阈值看看空间网络结构如何影响传播模式。探索不同场景下的疾病传播动态也许你会有新的发现哦
完整代码省流版
import random
import numpy as np
import pandas as pd
import geopandas as gpd
import libpysal
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter# Reproducibility
random.seed(33333334)
np.random.seed(33333334)# Load shapefiles for City of London and Camden
lon gpd.read_file(LSOA_2011_BFE_City_of_London.shp)
cam gpd.read_file(LSOA_2011_BFE_Camden.shp)
# Combine GeoDataFrames
df pd.concat([lon, cam], ignore_indexTrue)# Keep a WGS84 copy for plotting maps in lon/lat
df_ll df.to_crs(epsg4326)# Re-project to a metric CRS before centroid calculations (EPSG:3857)
df df_ll.to_crs(epsg3857)# Compute centroids in both CRS
centroids_3857 df.geometry.centroid
coords_3857 np.vstack([centroids_3857.x, centroids_3857.y]).T
centroids_ll df_ll.geometry.centroid
coords_ll np.vstack([centroids_ll.x, centroids_ll.y]).T# Build spatial weights using Queen contiguity
w libpysal.weights.Queen.from_dataframe(df, use_indexTrue)# Convert to weight matrix
n len(df)
W np.zeros((n, n), dtypeint)
for i, neighbors in w.neighbors.items():for j in neighbors:W[i, j] 1W[j, i] 1 # symmetric# Plot base map and centroids
plt.figure(figsize(6, 6))
df_ll.plot(facecolornone, edgecolorblack)
plt.scatter(coords_ll[:, 0], coords_ll[:, 1], s10)
plt.title(Polygons and Centroids (WGS84))
plt.show()# Histogram of positive weights
dist W[W 0]
plt.figure()
plt.hist(dist, bins20)
plt.title(Distribution of Weights 0)
plt.show()# Threshold adjacency
threshold 0.16
adj (W threshold).astype(int)# Plot edges above threshold
g_df df_ll.copy()
plt.figure(figsize(6, 6))
g_df.plot(facecolornone, edgecolorgray)
plt.scatter(coords_ll[:, 0], coords_ll[:, 1], s10, colorblack)
for i in range(n):for j in range(i 1, n):if adj[i, j]:plt.plot([coords_ll[i, 0], coords_ll[j, 0]],[coords_ll[i, 1], coords_ll[j, 1]],colorblue, linewidth0.5)
plt.title(Network Edges (thresholded))
plt.show()# SIR parameters
beta 3.0
gamma 1.0# Initialize status: S, I, R
status np.array([S] * n)
i0 random.randrange(n)
status[i0] I# Storage
T_max 10.0
times [0.0]
states [status.copy()]
adj_list {i: list(np.where(adj[i] 1)[0]) for i in range(n)}
current_time 0.0# Gillespie simulation
i 0
while current_time T_max:I_indices np.where(status I)[0]rec_rates {idx: gamma for idx in I_indices}S_indices np.where(status S)[0]inf_rates {idx: beta * sum(status[j] I for j in adj_list[idx]) for idx in S_indices}all_rates {**rec_rates, **inf_rates}total_rate sum(all_rates.values())if total_rate 0:breaknodes, rates zip(*all_rates.items())probs np.array(rates) / total_rateevent_node np.random.choice(nodes, pprobs)if status[event_node] S:status[event_node] Ielse:status[event_node] Rdt np.random.exponential(1.0 / total_rate)current_time dttimes.append(current_time)states.append(status.copy())i 1# Time series DataFrame
df_counts [{time: t,S: int((st S).sum()),I: int((st I).sum()),R: int((st R).sum())}for t, st in zip(times, states)]
df_ts pd.DataFrame(df_counts)plt.figure()
for compartment in [S, I, R]:plt.plot(df_ts[time], df_ts[compartment], labelcompartment)
plt.xlabel(Time)
plt.ylabel(Count)
plt.legend()
plt.title(SIR Dynamics Over Time)
plt.show()# Prepare DataFrame for animation
records []
for step, st in enumerate(states):for idx, s in enumerate(st):records.append({time_step: step,node: idx,status: s,x: coords_ll[idx, 0],y: coords_ll[idx, 1]})
anim_df pd.DataFrame(records)# Animate with PillowWriter
def animate_sir_map(output_filelondon.gif, fps5):fig, ax plt.subplots(figsize(6, 6))def update(frame):ax.clear()df_ll.plot(axax, facecolorlightgray, edgecolorwhite)df_step anim_df[anim_df[time_step] frame]for label, color in zip([S, I, R], [blue, red, green]):sub df_step[df_step[status] label]ax.scatter(sub[x], sub[y], labellabel, s20)ax.set_title(fTime step: {frame})ax.legend(locupper right, fontsizesmall)ax.set_axis_off()writer PillowWriter(fpsfps)anim FuncAnimation(fig, update, frameslen(states), interval200)anim.save(output_file, writerwriter)print(fAnimation saved as {output_file} using PillowWriter)animate_sir_map()
