# ABM建模实战：用Schelling模型解释城市社区分化（Python版） 如果你曾仔细观察过任何一座大城市的种族或收入分布地图，那些泾渭分明的色块边界可能会让你感到惊讶。芝加哥、纽约、伦敦……这些被视为文化熔炉的国际化都市，其内部社区却常常呈现出高度的同质化聚集。一个直观的解释是：居民存在强烈的排外倾向，只愿与同类为邻。但1971年，经济学家托马斯·谢林（Thomas Schelling）用一个极其简单的思想实验，颠覆了这种认知。他发现，即使个体只表现出轻微的“同类偏好”，在多次互动后，整个系统也会“涌现”出宏观上高度隔离的格局。这个如今被称为“谢林模型”（Schelling Model）的经典研究，不仅为他赢得了诺贝尔经济学奖，更成为了复杂系统科学和基于主体建模（Agent-Based Modeling, ABM）的奠基之作。 今天，我们不再需要用硬币和方格纸来推演这个模型。借助Python强大的科学计算生态，我们可以轻松地将谢林的思想转化为可运行、可调节、可视化的数字实验。本文正是这样一份面向复杂系统研究者和城市规划爱好者的实战指南。我们将从零开始，构建一个完整的谢林模型仿真，并深入探讨其参数空间、涌现现象以及对现实世界的启示。无论你是想理解ABM的建模范式，还是希望获得一个分析社会空间动态的量化工具，这里都有你需要的代码、思路和洞察。 ## 1. 理论基础：从微观偏好到宏观隔离的魔法 在深入代码之前，有必要厘清谢林模型的核心逻辑。它本质上是一个**基于主体的空间互动模型**。我们可以把城市想象成一个巨大的棋盘（比如50x50的网格），每个格子代表一个居住单元。棋盘上有两种类型的居民（例如用红色和蓝色表示），以及一定比例的空置格子。 每个居民（主体）遵循一条极其简单的行为规则： > **满意度检查**：查看自己周围八个相邻格子（摩尔邻域）中，与自己同类的邻居所占的比例。如果这个比例低于某个预设的“满意度阈值”，我就感到“不满意”，并决定搬家。 搬家过程同样简单：随机选择一个当前空置的格子，搬过去。之后，系统更新所有居民的状态，重复上述过程。模型会持续运行，直到没有居民感到不满意，或者达到预设的迭代次数。 谢林的惊人发现在于，**即使满意度阈值设置得相当低（例如，只要求30%的邻居是同类，其余70%可以是异类），经过多轮迁移后，整个棋盘依然会形成高度隔离的“红区”和“蓝区”**。这意味着，宏观上观察到的强烈隔离现象，并不必然源于个体强烈的排外心理。微弱的、甚至看似“包容”的个体偏好，通过空间上的局部互动和迭代，足以催生全局性的分化结构。这种“整体大于部分之和”的特性，正是复杂系统中**涌现**（Emergence）现象的典型体现。 为了后续建模，我们明确模型的几个关键参数： * **网格尺寸**：定义城市的大小。 * **主体类型与比例**：通常是两种类型，比例可以相等或不相等。 * **空置率**：空置格子的比例，为迁移提供空间。 * **满意度阈值**：核心参数，决定主体对多样性的容忍度。 * **邻域定义**：通常使用摩尔邻域（周围8格）或冯·诺依曼邻域（上下左右4格）。 ## 2. 环境搭建与模型初始化 我们选择Python作为实现语言，主要依赖`numpy`进行高效的矩阵运算，用`matplotlib`进行可视化。首先确保你的环境已安装这些库。 ```bash pip install numpy matplotlib ``` 接下来，我们开始构建模型。我们将采用面向对象的方式，创建一个`SchellingModel`类，这样代码更清晰，也便于管理状态和参数。 ```python import numpy as np import matplotlib.pyplot as plt from matplotlib import colors from typing import Tuple, List, Optional import copy class SchellingModel: """ 谢林隔离模型实现类。 """ def __init__(self, width: int = 50, height: int = 50, agent_ratio: Tuple[float, float] = (0.45, 0.45), empty_ratio: float = 0.1, similarity_threshold: float = 0.3, max_iterations: int = 500): """ 初始化模型参数。 参数: width: 网格宽度 height: 网格高度 agent_ratio: 两类主体的比例 (例如 (0.45, 0.45)， 剩余为空格) empty_ratio: 空置格子比例 (如果提供，将覆盖agent_ratio计算出的空置部分) similarity_threshold: 满意度阈值 (0到1之间) max_iterations: 最大迭代次数 """ self.width = width self.height = height self.total_cells = width * height # 计算各类格子数量 self.empty_count = int(self.total_cells * empty_ratio) self.agent1_count = int(self.total_cells * agent_ratio[0]) self.agent2_count = self.total_cells - self.empty_count - self.agent1_count self.similarity_threshold = similarity_threshold self.max_iterations = max_iterations # 状态编码: 0=空置, 1=主体类型A, 2=主体类型B self.grid = None self.iteration = 0 self.history = [] # 用于记录每轮不满意主体比例 self._initialize_grid() def _initialize_grid(self): """随机初始化网格状态。""" # 创建一个包含所有格子索引的列表 all_positions = np.arange(self.total_cells) np.random.shuffle(all_positions) # 分配状态 state_array = np.zeros(self.total_cells, dtype=np.int8) # 前agent1_count个为类型1 state_array[all_positions[:self.agent1_count]] = 1 # 接着agent2_count个为类型2 start = self.agent1_count end = self.agent1_count + self.agent2_count state_array[all_positions[start:end]] = 2 # 剩余部分保持为0 (空置) # 重塑为二维网格 self.grid = state_array.reshape((self.height, self.width)) ``` 这个初始化方法创建了一个指定大小的网格，并按照设定的比例随机放置两类主体和空置格子。我们使用整数`0`、`1`、`2`来分别代表空置、主体类型A和主体类型B，便于后续计算。 ## 3. 核心引擎：满意度计算与迁移逻辑 模型的核心在于计算每个主体的满意度，并让不满意的主体迁移。这里有几个细节需要注意： 1. **边界处理**：对于位于网格边缘的主体，其“邻居”数量会减少。我们采用常见的做法：只计算实际存在的邻居格子。 2. **满意度计算**：对于空置格子，我们不计算其满意度。 3. **迁移策略**：随机选择一个不满意的主体，并将其与一个随机选择的空置格子交换位置。这模拟了“搬入空房”的行为。 4. **同步与异步更新**：在原始谢林模型中，通常采用顺序更新（异步），即逐个检查并迁移主体，一轮结束后再统一更新状态。这更符合现实，但实现稍复杂。另一种是批量同步更新，即先找出所有不满意主体，然后一次性为他们分配新位置（需避免冲突）。我们这里实现更符合原意的顺序更新。 ```python class SchellingModel: # ... 接上文 __init__ 和 _initialize_grid ... def _get_neighbors(self, row: int, col: int) -> np.ndarray: """获取指定单元格周围所有邻居的状态（不包括自身）。""" # 定义8个方向 (摩尔邻域) directions = [(-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)] neighbors = [] for dr, dc in directions: nr, nc = row + dr, col + dc if 0 <= nr < self.height and 0 <= nc < self.width: neighbors.append(self.grid[nr, nc]) return np.array(neighbors) def _calculate_happiness(self, row: int, col: int) -> float: """计算单个主体的满意度分数。""" agent_type = self.grid[row, col] if agent_type == 0: # 空置格子无满意度 return np.nan neighbors = self._get_neighbors(row, col) # 过滤掉空置邻居 non_empty_neighbors = neighbors[neighbors != 0] if len(non_empty_neighbors) == 0: return 1.0 # 如果没有非空邻居，定义为满意（或根据需求调整） same_type_count = np.sum(non_empty_neighbors == agent_type) happiness = same_type_count / len(non_empty_neighbors) return happiness def _is_unhappy(self, row: int, col: int) -> bool: """判断一个主体是否不满意。""" happiness = self._calculate_happiness(row, col) if np.isnan(happiness): # 空置格子 return False return happiness < self.similarity_threshold def _find_unhappy_agents(self) -> List[Tuple[int, int]]: """找出所有不满意主体的坐标列表。""" unhappy = [] for i in range(self.height): for j in range(self.width): if self.grid[i, j] != 0 and self._is_unhappy(i, j): unhappy.append((i, j)) return unhappy def _find_empty_cells(self) -> List[Tuple[int, int]]: """找出所有空置格子的坐标列表。""" empty_cells = np.argwhere(self.grid == 0) return list(map(tuple, empty_cells)) def run_one_iteration(self) -> float: """ 运行一轮迭代。 返回本轮不满意主体的比例。 """ unhappy_agents = self._find_unhappy_agents() empty_cells = self._find_empty_cells() if not unhappy_agents or not empty_cells: unhappy_ratio = 0.0 else: # 随机打乱不满意主体列表，实现随机顺序更新 np.random.shuffle(unhappy_agents) moved_count = 0 for u_row, u_col in unhappy_agents: if not empty_cells: # 如果没有空位了，停止迁移 break # 随机选择一个空置格子 e_idx = np.random.randint(len(empty_cells)) e_row, e_col = empty_cells.pop(e_idx) # 交换位置 self.grid[e_row, e_col], self.grid[u_row, u_col] = \ self.grid[u_row, u_col], self.grid[e_row, e_col] # 由于原空置格子现在被占据，原主体位置变为空置 empty_cells.append((u_row, u_col)) moved_count += 1 unhappy_ratio = len(unhappy_agents) / (self.total_cells - self.empty_count) self.history.append(unhappy_ratio) self.iteration += 1 return unhappy_ratio def run(self, verbose: bool = False) -> None: """运行模型直到稳定或达到最大迭代次数。""" for i in range(self.max_iterations): unhappy_ratio = self.run_one_iteration() if verbose and i % 20 == 0: print(f"Iteration {i}: Unhappy ratio = {unhappy_ratio:.4f}") if unhappy_ratio == 0.0: if verbose: print(f"Model stabilized at iteration {i}.") break if verbose and self.iteration == self.max_iterations: print(f"Reached max iterations {self.max_iterations}.") ``` > **注意**：上述迁移逻辑中，一个主体搬走后，其原位置立即被标记为空置，并加入空置列表。这保证了空置格子池的动态更新，是模型能够收敛的关键。 ## 4. 可视化与结果分析：看见“涌现” 模型跑起来了，但我们更需要“看见”隔离是如何形成的。我们将创建几个可视化函数，并设计实验来探索参数的影响。 首先，是一个绘制当前网格状态的函数： ```python def plot_grid(model: SchellingModel, title: str = "", ax: Optional[plt.Axes] = None): """绘制模型网格状态。""" cmap = colors.ListedColormap(['white', 'red', 'blue']) bounds = [-0.5, 0.5, 1.5, 2.5] norm = colors.BoundaryNorm(bounds, cmap.N) if ax is None: fig, ax = plt.subplots(1, 1, figsize=(8, 8)) im = ax.imshow(model.grid, cmap=cmap, norm=norm, interpolation='nearest') ax.set_xticks([]) ax.set_yticks([]) ax.set_title(title + f" (Iteration {model.iteration})") # 添加图例 from matplotlib.patches import Patch legend_elements = [Patch(facecolor='white', edgecolor='black', label='Empty'), Patch(facecolor='red', label='Type A'), Patch(facecolor='blue', label='Type B')] ax.legend(handles=legend_elements, loc='upper left', bbox_to_anchor=(1.05, 1)) return ax ``` 现在，让我们运行一个基础实验，观察隔离的演化过程。我们将设置一个较低的满意度阈值（0.3），看看会发生什么。 ```python # 实验1：低阈值下的隔离形成 model_low = SchellingModel(width=40, height=40, agent_ratio=(0.4, 0.4), empty_ratio=0.2, similarity_threshold=0.3, max_iterations=200) fig, axes = plt.subplots(2, 3, figsize=(15, 10)) axes = axes.flatten() # 绘制初始状态 plot_grid(model_low, "Initial State", axes[0]) # 运行并记录中间状态 snapshot_iterations = [0, 5, 15, 50, 150] snapshot_idx = 1 for i in range(model_low.max_iterations): unhappy_ratio = model_low.run_one_iteration() if model_low.iteration in snapshot_iterations and model_low.iteration > 0: plot_grid(model_low, f"Iter {model_low.iteration}", axes[snapshot_idx]) snapshot_idx += 1 if unhappy_ratio == 0.0: plot_grid(model_low, f"Final (Iter {model_low.iteration})", axes[-1]) break if model_low.iteration == model_low.max_iterations: plot_grid(model_low, f"Final (Max Iter)", axes[-1]) plt.tight_layout() plt.show() # 绘制不满意比例随时间的变化曲线 plt.figure(figsize=(10, 5)) plt.plot(model_low.history, marker='o', markersize=4, linestyle='-', linewidth=1) plt.xlabel('Iteration') plt.ylabel('Proportion of Unhappy Agents') plt.title('Convergence Dynamics (Threshold=0.3)') plt.grid(True, alpha=0.3) plt.show() ``` 运行这段代码，你将看到一系列图像。初始状态是红蓝随机混合的“马赛克”。仅仅经过几次迭代，一些小的同色簇就开始出现。随着迭代进行，这些簇像滚雪球一样增长、合并，最终形成几个大的、边界清晰的同色区域。而不满意主体的比例曲线则会从初始的较高值迅速下降，最终趋近于零，表明系统达到了稳定状态（所有主体都满意）。 ## 5. 参数空间探索：阈值、空置率与相变 谢林模型的魅力在于其丰富的参数空间。改变关键参数，系统的宏观表现会截然不同。我们设计一个实验来系统性地探索满意度阈值的影响。 ```python def run_experiment_for_thresholds(thresholds: List[float], n_runs: int = 5, grid_size: int = 30) -> dict: """ 对不同阈值进行多次实验，收集统计结果。 返回: dict: 键为阈值，值为字典，包含平均相似度、平均收敛步数等。 """ results = {} for th in thresholds: print(f"Running experiments for threshold = {th:.2f}") all_similarities = [] all_steps = [] for run in range(n_runs): model = SchellingModel(width=grid_size, height=grid_size, agent_ratio=(0.45, 0.45), empty_ratio=0.1, similarity_threshold=th, max_iterations=300) model.run(verbose=False) # 计算最终状态的平均相似度 (隔离程度) similarities = [] for i in range(grid_size): for j in range(grid_size): if model.grid[i, j] != 0: sim = model._calculate_happiness(i, j) if not np.isnan(sim): similarities.append(sim) mean_similarity = np.mean(similarities) if similarities else 0 all_similarities.append(mean_similarity) all_steps.append(model.iteration) results[th] = { 'mean_similarity': np.mean(all_similarities), 'std_similarity': np.std(all_similarities), 'mean_steps': np.mean(all_steps), 'std_steps': np.std(all_steps), } return results # 定义要测试的阈值范围 thresholds_to_test = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8] exp_results = run_experiment_for_thresholds(thresholds_to_test, n_runs=3, grid_size=40) # 将结果整理成表格展示 import pandas as pd df_results = pd.DataFrame.from_dict(exp_results, orient='index') df_results.columns = ['Avg Similarity', 'Std Similarity', 'Avg Steps to Converge', 'Std Steps'] print("\n实验结果汇总：") print(df_results.round(3)) ``` 运行这个实验后，我们可以用图表更直观地展示结果： ```python fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5)) # 图表1：平均相似度 vs 满意度阈值 thresholds = list(exp_results.keys()) mean_sim = [exp_results[t]['mean_similarity'] for t in thresholds] std_sim = [exp_results[t]['std_similarity'] for t in thresholds] ax1.errorbar(thresholds, mean_sim, yerr=std_sim, fmt='-o', capsize=5, linewidth=2) ax1.axhline(y=0.5, color='gray', linestyle='--', alpha=0.5, label='Random Mix (0.5)') ax1.set_xlabel('Satisfaction Threshold') ax1.set_ylabel('Average Similarity (Segregation Level)') ax1.set_title('Threshold vs. Emergent Segregation') ax1.grid(True, alpha=0.3) ax1.legend() # 图表2：收敛速度 vs 满意度阈值 mean_steps = [exp_results[t]['mean_steps'] for t in thresholds] std_steps = [exp_results[t]['std_steps'] for t in thresholds] ax2.errorbar(thresholds, mean_steps, yerr=std_steps, fmt='-s', capsize=5, color='orange', linewidth=2) ax2.set_xlabel('Satisfaction Threshold') ax2.set_ylabel('Average Iterations to Converge') ax2.set_title('Threshold vs. Convergence Speed') ax2.grid(True, alpha=0.3) plt.tight_layout() plt.show() ``` 通过分析这些图表，你会发现一些非线性的、反直觉的现象： | 阈值范围 | 平均相似度 (隔离程度) | 收敛速度 | 宏观表现 | | :--- | :--- | :--- | :--- | | **很低 (如 0.1-0.2)** | 接近随机混合 (~0.5) | 较快 | 融合良好，无明显大集群 | | **中等 (如 0.3-0.5)** | **显著升高 (如 0.7-0.9)** | 中等 | 出现明显的隔离区域 | | **高 (如 0.6-0.7)** | 依然较高，但可能略低于中等阈值 | **变慢，甚至震荡** | 隔离明显，但边界可能更破碎 | | **极高 (如 >=0.8)** | **可能下降** | 极慢，难以收敛 | 主体频繁搬家，无法形成稳定社区 | > **关键洞察**：存在一个**相变**（Phase Transition）区间。当阈值超过某个临界点（例如0.65附近），系统可能无法达到所有人都满意的稳定状态，而是陷入动态平衡——总有一部分人在搬家。这是因为要求过于苛刻，导致“满意”的区位变得极其稀缺。此外，极高的阈值反而可能导致宏观隔离程度下降，因为主体过于挑剔，无法在任何地方安定下来形成大的集群。 空置率是另一个关键参数。更高的空置率就像给城市提供了更多的“缓冲空间”，让不满意的主体有更多选择，通常会**加速收敛**，但有时也可能让隔离模式变得更加**碎片化**，因为主体更容易找到附近的小块满意区域，而不必长途迁徙加入大集群。你可以通过修改上面的实验代码，将`empty_ratio`作为变量进行类似的探索。 ## 6. 超越基础：模型扩展与城市规划启示 基础谢林模型已经揭示了深刻的原理，但我们完全可以对其进行扩展，使其更贴近现实，或探索更复杂的问题。 **扩展一：异质性主体** 现实中的居民并非同质。我们可以让不同类的主体拥有不同的满意度阈值，甚至让同一类中的不同个体也拥有阈值分布。这可以通过在初始化时为每个主体随机分配一个阈值来实现。这种扩展能模拟社会中“宽容者”和“挑剔者”共存的情况。 **扩展二：多维属性与偏好** 主体可能不仅关心邻居的“种族”，还在意其“收入水平”、“教育程度”等。我们可以定义主体的多维标签，并设计更复杂的满意度函数（例如，加权综合评分）。这需要更复杂的数据结构，但Python的面向对象特性使其易于实现。 **扩展三：连续空间与距离衰减** 网格模型是离散的。更现实的模型可以使用连续空间，主体的满意度随邻居距离的增加而衰减。这需要用到更复杂的空间计算库（如`scipy.spatial`），但能研究城市蔓延、郊区化等过程。 **扩展四：引入外部干预（政策模拟）** 这是ABM最具价值的方向之一。我们可以模拟城市规划政策的影响，例如： * **建设混合收入住房**：在网格中特定位置强制放置不同类型的主体。 * **设置社区多样性配额**：修改满意度规则，对低于多样性要求的社区进行“惩罚”或提供“补贴”。 * **改善交通**：改变主体的“可迁移范围”，从全局随机搬迁变为在局部范围内寻找空位。 例如，一个模拟“公平住房政策”的代码片段可能如下： ```python def apply_fair_housing_policy(model, policy_strength=0.1): """ 一个简单的政策干预：随机选择一小部分空置格子， 并将其强制分配给‘少数’群体，以促进融合。 """ empty_cells = model._find_empty_cells() if not empty_cells: return # 确定哪类主体是“少数”群体（数量较少的） count_a = np.sum(model.grid == 1) count_b = np.sum(model.grid == 2) minority_type = 1 if count_a < count_b else 2 # 根据政策强度，决定干预的格子数量 n_to_intervene = max(1, int(len(empty_cells) * policy_strength)) cells_to_fill = np.random.choice(len(empty_cells), size=n_to_intervene, replace=False) for idx in cells_to_fill: r, c = empty_cells[idx] model.grid[r, c] = minority_type print(f"Policy applied: assigned {n_to_intervene} cells to minority group.") ``` 将这些扩展融入模型后，我们就能从一个解释性的工具，转变为一个**政策实验室**，在数字世界中低成本、快速地测试不同干预措施的长期效果和潜在副作用。 谢林模型半个世纪前的洞见，在今天大数据和计算模拟的加持下，依然焕发着强大的生命力。它提醒我们，面对城市社区分化这样的复杂问题，线性思维和简单归因是危险的。宏观格局是无数微观互动在时间中编织出的“涌现”之网。通过构建并探索这样的计算模型，我们并非寻求一个确定的预测，而是获得一种**系统性的思考方式**——一种理解“模式如何形成”以及“干预可能如何奏效”的直觉。这或许，正是ABM建模带给复杂系统研究者最宝贵的礼物。