在當(dāng)今世界，傳染病的傳播模擬已成為公共衛(wèi)生領(lǐng)域的重要研究工具。傳統(tǒng)的傳染病模型往往忽略了空間因素，將整個群體視為一個均質(zhì)的整體。然而，現(xiàn)實中的疾病傳播與地理空間位置密切相關(guān)，相鄰區(qū)域之間的傳播風(fēng)險遠高于遠距離區(qū)域。今天，我們將通過一個結(jié)合地理信息系統(tǒng) (GIS) 和隨機模擬算法的案例，深入探討如何在空間維度上模擬傳染病的傳播過程。

引言：當(dāng)傳染病模型遇上地理空間

傳染病的傳播從來都不是抽象的數(shù)字游戲，它實實在在地發(fā)生在我們生活的城市空間中。想象一下，在倫敦這樣的大都市中，不同區(qū)域之間的人口流動、空間鄰近性如何影響疾病的傳播速度和范圍？這正是我們今天要探索的主題。

我們將使用經(jīng)典的 SIR 傳染病模型，但會為其添加 "地理空間" 這一重要維度。SIR 模型將人群分為三類：易感者 (Susceptible)、感染者 (Infectious) 和康復(fù)者 (Recovered)。通過引入空間權(quán)重矩陣，我們可以量化不同區(qū)域之間的鄰近關(guān)系，使模型能夠更真實地反映疾病在地理空間上的傳播過程。

數(shù)據(jù)準(zhǔn)備：加載倫敦的地理空間數(shù)據(jù)

首先，我們需要準(zhǔn)備地理空間數(shù)據(jù)。代碼中加載了倫敦市 (City of London) 和卡姆登區(qū) (Camden) 的低級統(tǒng)計區(qū) (LSOA) 邊界數(shù)據(jù)。這些數(shù)據(jù)是英國國家統(tǒng)計數(shù)據(jù)的一部分，能夠幫助我們精確界定不同的地理單元。

# 加載倫敦市和卡姆登區(qū)的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_index=True)

地理數(shù)據(jù)通常有不同的坐標(biāo)參考系統(tǒng) (CRS)。為了確保后續(xù)計算的準(zhǔn)確性，我們需要進行坐標(biāo)轉(zhuǎn)換。我們先保留一個 WGS84 (經(jīng)緯度) 坐標(biāo)系的副本用于地圖繪制，然后將數(shù)據(jù)轉(zhuǎn)換為等距坐標(biāo)系 (EPSG:3857) 用于質(zhì)心計算和距離測量。?

# 保留WGS84副本用于繪制經(jīng)緯度地圖
df_ll = df.to_crs(epsg=4326)# 在進行質(zhì)心計算前重新投影到等距CRS(EPSG:3857)
df = df_ll.to_crs(epsg=3857)

坐標(biāo)系統(tǒng)的轉(zhuǎn)換是地理空間分析中非常關(guān)鍵的一步。想象一下，如果我們在經(jīng)緯度坐標(biāo)系中直接計算兩點之間的距離，由于地球曲率的影響，結(jié)果會有很大誤差。而轉(zhuǎn)換到等距坐標(biāo)系后，我們就可以像在平面上一樣準(zhǔn)確計算距離和面積了 📐。

?

接下來，我們計算每個區(qū)域的質(zhì)心，并獲取它們在兩種坐標(biāo)系中的坐標(biāo)：

# 計算兩種CRS下的質(zhì)心
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

空間權(quán)重矩陣：定義區(qū)域間的鄰近關(guān)系

在空間分析中，"鄰近性" 是一個關(guān)鍵概念。如何定義兩個區(qū)域是否相鄰？在代碼中，我們使用了 "Queen contiguity" 規(guī)則，這是空間統(tǒng)計中常用的一種鄰近性定義。Queen contiguity 認(rèn)為，如果兩個區(qū)域共享一條邊或一個角，則它們是相鄰的。

# 使用Queen鄰接構(gòu)建空間權(quán)重
w = libpysal.weights.Queen.from_dataframe(df, use_index=True)# 轉(zhuǎn)換為權(quán)重矩陣
n = len(df)
W = np.zeros((n, n), dtype=int)
for i, neighbors in w.neighbors.items():for j in neighbors:W[i, j] = 1W[j, i] = 1  # 對稱矩陣

這種鄰接關(guān)系定義很有意思，它模擬了一種 "棋盤式" 的鄰近性。想象一下國際象棋中的皇后可以沿對角線移動，Queen contiguity 就是用這種方式定義鄰近性的 � chess。與只考慮共享邊的 "Rook contiguity" 相比，它考慮了更多的鄰近關(guān)系。

為了更直觀地理解我們構(gòu)建的權(quán)重矩陣，我們可以繪制權(quán)重分布的直方圖：

# 正權(quán)重的直方圖
dist = W[W > 0]
plt.figure()
plt.hist(dist, bins=20)
plt.title('Distribution of Weights > 0')
plt.show()

由于我們使用的是二元鄰接矩陣 (相鄰為 1，不相鄰為 0)，這個直方圖應(yīng)該會呈現(xiàn)出明顯的雙峰分布，其中大部分權(quán)重為 0，而相鄰區(qū)域的權(quán)重為 1。不過，在后續(xù)步驟中，我們還會引入閾值處理，讓鄰接關(guān)系變得更加 "模糊"。

構(gòu)建空間網(wǎng)絡(luò)：閾值化鄰接關(guān)系

有時候，我們可能需要考慮更靈活的鄰近關(guān)系定義。例如，不僅考慮直接相鄰的區(qū)域，還可以根據(jù)距離或其他因素定義一個 "影響范圍"。在代碼中，我們通過設(shè)置一個閾值，將權(quán)重矩陣轉(zhuǎn)換為一個二元鄰接矩陣：

# 閾值鄰接
threshold = 0.16
adj = (W > threshold).astype(int)

這里的閾值設(shè)置是一個關(guān)鍵參數(shù)，它決定了哪些區(qū)域會被視為 "相鄰"。閾值越高，鄰接關(guān)系越嚴(yán)格，網(wǎng)絡(luò)中的連接越少；閾值越低，鄰接關(guān)系越寬松，網(wǎng)絡(luò)中的連接越多。這個參數(shù)的設(shè)置需要根據(jù)具體問題和數(shù)據(jù)特點來調(diào)整 🧩。

我們可以將這些鄰接關(guān)系可視化，看看在地理空間上區(qū)域之間是如何連接的：

# 繪制高于閾值的邊
g_df = df_ll.copy()
plt.figure(figsize=(6, 6))
g_df.plot(facecolor='none', edgecolor='gray')
plt.scatter(coords_ll[:, 0], coords_ll[:, 1], s=10, color='black')
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]],color='blue', linewidth=0.5)
plt.title('Network Edges (thresholded)')
plt.show()

?

看著這些藍色的連接線在地圖上延伸，我們仿佛看到了疾病可能傳播的 "高速公路"。這些連接線構(gòu)成了一個空間網(wǎng)絡(luò)，疾病將在這個網(wǎng)絡(luò)上進行傳播 🌐。

SIR 模型：從經(jīng)典到空間擴展

現(xiàn)在，我們來介紹 SIR 模型的基本原理。SIR 模型是一個經(jīng)典的傳染病模型，它將人群分為三個狀態(tài)：

• S (Susceptible): 易感者，尚未感染但可能被感染的人
• I (Infectious): 感染者，具有傳染性的人
• R (Recovered): 康復(fù)者，已經(jīng)康復(fù)并獲得免疫力的人

模型的核心是兩個速率參數(shù)：

• β(beta): 感染率，表示易感者與感染者接觸后被感染的概率
• γ(gamma): 康復(fù)率，表示感染者康復(fù)的概率

在傳統(tǒng)的 SIR 模型中，通常假設(shè)人群是完全混合的，即每個人與其他人接觸的概率是相等的。但在我們的空間擴展版本中，這個假設(shè)被打破了。我們考慮了空間鄰近性，即一個區(qū)域的感染狀態(tài)會受到其相鄰區(qū)域感染狀態(tài)的影響。

在代碼中，我們設(shè)置了感染率和康復(fù)率參數(shù)：

# SIR參數(shù)
beta = 3.0  # 感染率
gamma = 1.0  # 康復(fù)率

這些參數(shù)的值會顯著影響疾病傳播的動態(tài)。例如，較高的感染率會導(dǎo)致疾病更快地傳播，而較高的康復(fù)率會縮短感染期，減緩疾病的傳播速度 ??。

初始狀態(tài)設(shè)置：疫情的起點

每個傳染病模擬都需要一個起點。在代碼中，我們隨機選擇一個區(qū)域作為初始感染區(qū)域：

# 初始化狀態(tài): 'S', 'I', 'R'
status = np.array(['S'] * n)
i0 = random.randrange(n)
status[i0] = 'I'

這個初始感染點的選擇雖然是隨機的，但在現(xiàn)實中，疫情的起點可能與人口密度、交通樞紐等因素有關(guān)。如果我們有更詳細(xì)的數(shù)據(jù)，我們可以根據(jù)實際情況選擇更合理的初始感染點 🦠。

Gillespie 隨機模擬算法：捕捉傳染病的隨機性

傳染病的傳播本質(zhì)上是一個隨機過程，因為每個接觸事件是否導(dǎo)致感染都是概率性的。為了準(zhǔn)確捕捉這種隨機性，我們使用了 Gillespie 隨機模擬算法，這是一種在分子尺度上模擬隨機過程的強大方法。

Gillespie 算法的核心思想是：

1. 計算系統(tǒng)中所有可能事件的發(fā)生速率
2. 根據(jù)這些速率計算每個事件發(fā)生的概率
3. 隨機選擇一個事件發(fā)生
4. 計算到下一個事件發(fā)生的時間間隔
5. 重復(fù)上述過程

在我們的 SIR 模型中，有兩種可能的事件：

1. 感染事件：易感者被感染者感染
2. 康復(fù)事件：感染者康復(fù)

讓我們看看代碼中是如何實現(xiàn)這個算法的：

# 存儲模擬結(jié)果
T_max = 10.0  # 最大模擬時間
times = [0.0]  # 時間點列表
states = [status.copy()]  # 每個時間點的狀態(tài)列表
adj_list = {i: list(np.where(adj[i] == 1)[0]) for i in range(n)}  # 鄰接列表
current_time = 0.0  # 當(dāng)前時間# Gillespie模擬
i = 0
while current_time < T_max:# 找出所有感染者I_indices = np.where(status == 'I')[0]# 計算每個感染者的康復(fù)速率rec_rates = {idx: gamma for idx in I_indices}# 找出所有易感者S_indices = np.where(status == 'S')[0]# 計算每個易感者的感染速率，與相鄰感染者數(shù)量有關(guān)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())# 如果沒有事件發(fā)生，結(jié)束模擬if total_rate <= 0:break# 提取事件節(jié)點和速率nodes, rates = zip(*all_rates.items())# 計算每個事件發(fā)生的概率probs = np.array(rates) / total_rate# 隨機選擇一個事件event_node = np.random.choice(nodes, p=probs)# 更新狀態(tài)if status[event_node] == 'S':status[event_node] = 'I'  # 易感者被感染else:status[event_node] = 'R'  # 感染者康復(fù)# 計算到下一個事件的時間間隔dt = np.random.exponential(1.0 / total_rate)current_time += dt# 存儲當(dāng)前時間和狀態(tài)times.append(current_time)states.append(status.copy())i += 1

Gillespie 算法的一個重要特點是它精確地模擬了隨機過程的時間演化，而不是像歐拉法那樣使用固定的時間步長。這使得它在處理具有不同時間尺度的事件時特別有效 ?。

在我們的實現(xiàn)中，每個易感者的感染速率不僅取決于感染率 β，還取決于其相鄰區(qū)域中有多少個感染區(qū)域。這正是空間因素被引入模型的關(guān)鍵所在 🌍。

結(jié)果分析：傳染病的時間動態(tài)

模擬結(jié)束后，我們可以分析傳染病在時間上的傳播動態(tài)。首先，我們創(chuàng)建一個時間序列數(shù)據(jù)框，記錄每個時間點上易感者、感染者和康復(fù)者的數(shù)量：

# 時間序列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], label=compartment)
plt.xlabel('Time')
plt.ylabel('Count')
plt.legend()
plt.title('SIR Dynamics Over Time')
plt.show()

?

典型的 SIR 模型時間曲線應(yīng)該呈現(xiàn)出這樣的模式：

• 易感者數(shù)量 (S) 隨時間單調(diào)遞減
• 感染者數(shù)量 (I) 先增加后減少，形成一個峰值
• 康復(fù)者數(shù)量 (R) 隨時間單調(diào)遞增

看著這些曲線的變化，我們可以直觀地理解傳染病的傳播過程：從少數(shù)人感染開始，逐漸擴散到高峰，然后隨著越來越多的人康復(fù)，感染人數(shù)逐漸減少 📈📉。

空間傳播動畫：疾病在地理空間上的擴散

最令人興奮的部分莫過于將疾病的空間傳播過程可視化。通過創(chuàng)建動畫，我們可以直觀地看到疾病如何從初始感染點開始，沿著空間網(wǎng)絡(luò)向周圍區(qū)域擴散。

首先，我們準(zhǔn)備動畫所需的數(shù)據(jù)：

# 準(zhǔn)備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)

然后，我們實現(xiàn)動畫生成函數(shù)：?

# 使用PillowWriter生成動畫
def animate_sir_map(output_file='london.gif', fps=5):fig, ax = plt.subplots(figsize=(6, 6))def update(frame):ax.clear()df_ll.plot(ax=ax, facecolor='lightgray', edgecolor='white')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'], label=label, s=20)ax.set_title(f"Time step: {frame}")ax.legend(loc='upper right', fontsize='small')ax.set_axis_off()writer = PillowWriter(fps=fps)anim = FuncAnimation(fig, update, frames=len(states), interval=200)anim.save(output_file, writer=writer)print(f"Animation saved as {output_file} using PillowWriter")# 生成動畫
animate_sir_map()

?

當(dāng)我們播放這個動畫時，可以看到紅色的感染區(qū)域如何從一個點開始，逐漸擴散到周圍的藍色易感區(qū)域，然后變成綠色的康復(fù)區(qū)域。這種時空動態(tài)直觀地展示了空間鄰近性如何影響疾病的傳播模式 🌍→🦠→💊。

模型的局限性與改進方向

盡管我們的空間 SIR 模型比傳統(tǒng)模型更貼近現(xiàn)實，但它仍然有一些局限性：

1. 空間權(quán)重的簡化：我們使用了簡單的二元鄰接矩陣來表示空間關(guān)系，而現(xiàn)實中的疾病傳播可能受到多種因素影響，如人口流動、交通網(wǎng)絡(luò)等。

2. 同質(zhì)性假設(shè)：模型假設(shè)每個區(qū)域的人口特征是相同的，而實際上不同區(qū)域的人口密度、年齡結(jié)構(gòu)等可能差異很大。

3. 確定性參數(shù)：我們使用了固定的感染率和康復(fù)率，而在現(xiàn)實中這些參數(shù)可能隨時間和空間變化。

針對這些局限性，我們可以考慮以下改進方向：

1. 更復(fù)雜的空間權(quán)重：可以基于實際交通數(shù)據(jù)或人口流動數(shù)據(jù)構(gòu)建更真實的空間權(quán)重矩陣。

2. 異質(zhì)性建模：可以為不同區(qū)域設(shè)置不同的人口特征和疾病傳播參數(shù)。

3. 引入行為干預(yù)：可以模擬社交距離、封鎖等干預(yù)措施對疾病傳播的影響 🛡?。

4. 考慮人口流動：可以引入更復(fù)雜的人口流動模型，模擬通勤等日?；顒訉膊鞑サ挠绊憽?/p>

結(jié)語：空間視角下的傳染病防控

通過這個案例，我們看到了地理空間視角在傳染病建模中的重要性?？臻g鄰近性、區(qū)域連接性等因素顯著影響著疾病的傳播模式和速度。在公共衛(wèi)生政策制定中，考慮這些空間因素可以幫助我們更精準(zhǔn)地分配資源、實施干預(yù)措施。

想象一下，當(dāng)我們能夠在地圖上實時看到疾病的傳播路徑和高風(fēng)險區(qū)域時，公共衛(wèi)生部門就可以更有針對性地部署醫(yī)療資源、實施隔離措施，從而最大限度地減少疾病的影響 🏥。

這個模型不僅是一個學(xué)術(shù)研究工具，也為我們理解現(xiàn)實世界中的傳染病傳播提供了一個窗口。隨著地理信息技術(shù)和計算能力的不斷發(fā)展，我們有理由相信，未來的傳染病模型將更加精確、更加貼近現(xiàn)實，為全球公共衛(wèi)生安全提供更強有力的支持 🌟。

如果你對這個模型感興趣，不妨嘗試調(diào)整參數(shù)看看會發(fā)生什么變化：增加或減少感染率 β，看看疾病傳播速度如何變化；調(diào)整空間權(quán)重的閾值，看看空間網(wǎng)絡(luò)結(jié)構(gòu)如何影響傳播模式。探索不同場景下的疾病傳播動態(tài)，也許你會有新的發(fā)現(xiàn)哦！ 🧪

完整代碼（省流版）

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_index=True)# Keep a WGS84 copy for plotting maps in lon/lat
df_ll = df.to_crs(epsg=4326)# Re-project to a metric CRS before centroid calculations (EPSG:3857)
df = df_ll.to_crs(epsg=3857)# 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_index=True)# Convert to weight matrix
n = len(df)
W = np.zeros((n, n), dtype=int)
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(facecolor='none', edgecolor='black')
plt.scatter(coords_ll[:, 0], coords_ll[:, 1], s=10)
plt.title('Polygons and Centroids (WGS84)')
plt.show()# Histogram of positive weights
dist = W[W > 0]
plt.figure()
plt.hist(dist, bins=20)
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(facecolor='none', edgecolor='gray')
plt.scatter(coords_ll[:, 0], coords_ll[:, 1], s=10, color='black')
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]],color='blue', linewidth=0.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, p=probs)if status[event_node] == 'S':status[event_node] = 'I'else:status[event_node] = 'R'dt = 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], label=compartment)
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_file='london.gif', fps=5):fig, ax = plt.subplots(figsize=(6, 6))def update(frame):ax.clear()df_ll.plot(ax=ax, facecolor='lightgray', edgecolor='white')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'], label=label, s=20)ax.set_title(f"Time step: {frame}")ax.legend(loc='upper right', fontsize='small')ax.set_axis_off()writer = PillowWriter(fps=fps)anim = FuncAnimation(fig, update, frames=len(states), interval=200)anim.save(output_file, writer=writer)print(f"Animation saved as {output_file} using PillowWriter")animate_sir_map()