简介这是一份基于组合社会力模型CSFM与高斯烟团模型的地铁站毒气袭击紧急疏散模拟系统的Python复现与讲解文档面向具有一定编程基础的研究人员、工程师及应急疏散建模学者。资源包内为单个docx文档整包约60KB包含可运行的CSFM类实现、模型初始化、毒气源处理、社会力计算、气体浓度更新与可视化等完整代码并逐段附有解释和论文复现分析。读者可据此模拟毒气源位置与数量、管理人员反应速度、风速等因素对疏散时间及伤亡率的影响评估不同应急响应措施为地铁站设计和应急预案提供量化参考。文档还梳理了模型的数学基础、改进方向、验证方法及多场景应用潜力兼具科研对照与工程参考价值。目前已有80人学习下载适合对行人疏散与气体扩散仿真感兴趣的读者深入研读。1. 为什么突发毒气场景下普通社会力模型不够用某次我跑一个换乘站晚高峰疏散仿真时用经典Helbing社会力模型得出一个漂亮的结论所有人从两端站台向三个出口均匀撤离拥挤点稳定出现在闸机前。这个结果放在普通火灾里还能看但换成毒气泄漏真实监控里出现的画面完全不是这样——有人逆着人群往毒源方向走有人贴着墙根蹲下也有人在出口已经畅通时突然折返去扶同伴。普通社会力模型只描述“我想去出口其他人挡我”它假设每个行人在整个仿真里决策恒定。毒气场景恰恰破坏了这个假设浓度改变人的移动能力剂量主导人的决策。CSFM组合社会力模型把“毒气浓度影响行为权重”做成受力方程的一部分和高斯烟团给的浓度场耦合后才能复现这些反直觉行为。本文讲的这套模拟系统核心就是两个模型怎么各自实现、怎么在一个事件循环里互相喂数据以及参数调到什么量级结果才不骗人。2. 组合社会力模型 CSFM与经典社会力模型的差异落点CSFM不是一个能追溯到某篇固定论文的统一公式工程实现上一般是在Helbing提出的社会力模型骨架上做三件组合把期望速度的方向从“固定朝出口”改成多方向加权把恐慌系数做成随毒剂量变化的动态量把行人状态从单一行走扩成包含中毒倒地、折返救援的状态机。2.1 从Helbing方程出发CSFM在哪些项上动手经典社会力模型按牛顿体系展开单个行人的加速度由三部分叠加m_i * dv_i/dt m_i * (v_i^0 * e_i - v_i) / τ Σ_j f_ij Σ_w f_iw等号右边第一项是自驱动力描述行人希望以期望速度 v_i^0 朝方向 e_i 走但受身体反应时间 τ 限制实际速度 v_i 会滞后。第二项是行人与行人之间的社会排斥力第三项是墙壁等障碍物对行人的排斥力。行人间的物理接触力挤压与摩擦在高密度下再叠加进 f_ij。CSFM在这套体系上做改动时首先要动的就是 e_i。经典模型里 e_i 是一个常量方向即从当前位置指向出口的连线方向CSFM里它被重写为出口方向、避毒梯度方向、从众方向三个分量的归一化加权和。第二个大改动是 τ反应时间在恐慌状态下会骤降行人动作变得“又急又僵硬”计算公式常写成 τ τ0 * exp(-η * D_i)D_i是累计毒剂量。第三个改动是引入状态机毒剂量超过阈值后行人进入中毒减速或倒地状态这是一个离散事件不能靠连续力场平滑表达。2.2 状态机参数和“行为权重”怎么放进受力方程我一般把每个行人设置成四个状态normal、panic、intoxicated、down。normal到panic的切换条件是检测到浓度或他人恐慌比例超阈值panic到intoxicated的切换条件是Haber法则下的累计剂量超阈值。Haber法则指毒性效应近似等于浓度与暴露时间的乘积写成 D_i(t) ∫ C_i(t) dt这里的 C_i(t) 是行人当前位置的毒气浓度。行为权重的优先级也在这个状态机上定义panic状态下避毒梯度权重最高normal状态下出口方向权重最高intoxicated状态所有速度方向的合成会叠加巨大噪声。从众方向不是全局随机而是取以该行人为圆心、半径3米内所有邻居平均速度的方向这个方向存在局部信息传播的语义也是CSFM与仅靠浓度场驱动力模型的最大区别。2.3 最小Pedestrian实现下面是一个可直接用于仿真核心的最小Python实现省略了碰撞邻居搜索的优化但保留了状态切换和期望方向合成的完整逻辑。import numpy as np class Pedestrian: def __init__(self, pos, goal, v01.2, tau00.5, radius0.2, mass80.0): self.pos np.asarray(pos, dtypefloat) self.goal np.asarray(goal, dtypefloat) self.v np.zeros(2) self.v0 v0 self.tau0 tau0 self.radius radius self.mass mass self.rho 1.0 # 恐慌系数初始为1 self.dose 0.0 # 累计毒剂量 self.state normal def set_panicky(self, factor): self.rho max(1.0, factor) def compute_e(self, grad_c, neighbors_v, w_exit1.0, w_avoid0.4, w_follow0.2): e_exit (self.goal - self.pos) e_exit / np.linalg.norm(e_exit) 1e-9 e_avoid np.zeros(2) if np.linalg.norm(grad_c) 1e-6: e_avoid -grad_c / (np.linalg.norm(grad_c) 1e-9) e_follow np.zeros(2) if len(neighbors_v) 0: e_follow np.mean(neighbors_v, axis0) n np.linalg.norm(e_follow) if n 1e-6: e_follow / n e w_exit * e_exit w_avoid * e_avoid w_follow * e_follow n np.linalg.norm(e) return e / (n 1e-9) def step(self, grad_c, neighbors_v, contact_force, dt): if self.state down: self.v[:] 0.0 return tau self.tau0 * np.exp(-0.15 * (self.rho - 1.0)) e self.compute_e(grad_c, neighbors_v) v_desired self.v0 * e if self.state intoxicated: v_desired * 0.4 acc (v_desired - self.v) / tau contact_force / self.mass self.v acc * dt self.pos self.v * dt代码里 compute_e 方法把三个方向加权后归一化然后自驱动力按社会力模型的标准形式合成。注意 panic 状态和 intoxicated 状态分开处理panic 时反应时间 tau 变小行人动作更急intoxicated 时直接把期望速度打折模拟中毒后肌肉乏力。这个顺序不能反。引入接触力时我用的是Helbing原始参数量级的推荐表单位统一为米、千克、秒参数含义典型值校准依据v0期望速度1.2 m/s正常/ 1.8 m/s恐慌站台通道步行实测τ反应时间0.5 s正常/ 0.3 s恐慌Helbing标定区间A社会排斥力强度2000 N原始模型推荐量级B排斥力作用距离0.08 m原始模型推荐量级k身体弹性系数1.2e5 kg/s²接触力线性项κ摩擦系数2.4e5 kg/(m·s)接触力切向项数值更新建议使用Heun方法即先预估再校正显式欧拉在 τ 较小、Δt0.02s 时容易产生速度振荡将时间步控制在 0.01~0.02 秒一秒仿真推进 50~100 步1000人的场景在纯Python下仍能接近实时。3. 高斯烟团模型地铁站毒气浓度场的计算骨架浓度场是这套系统里另一个独立零件。疏散仿真中最常用的是高斯烟羽模型和高斯烟团模型两者都基于湍流扩散的正态分布假设。烟羽模型描述连续泄漏源烟团模型描述瞬时泄漏源。地铁站内毒气袭击更接近突发瞬时释放而且机械通风会让风向在几个通风周期内变化用多个烟团叠加以时间轴推进比拟合一个稳定烟羽更自然。3.1 为什么选烟团而不是烟羽烟羽模型假设源强恒定、气象条件平稳输出一个不随时间变化的浓度场。应急疏散场景里最关心的恰恰是“释放后第30秒站台两端浓度分别多少”这要求浓度场是时间函数。烟团的数学形式自带时间维每个释放脉冲表达为一个以释放点为中心、随风飘移且不断膨胀的三维高斯钟多个脉冲叠加就是非平稳浓度场。另一个原因是地铁站内部结构对风向的扰动太大。烟羽模型要求给出稳定的风速和风向而站厅内气流由通风系统主导风速通常低于0.5 m/s风向不定。烟团模型把空气运动拆成平均风场位移和湍流扩散两个独立部分机械通风导致的空气交换可以等效成一个额外的衰减系数不需要精确刻画具体流场。3.2 浓度公式和地面反射项解释单次瞬时释放的高斯烟团浓度表达式为C(x,y,z,t) Q / [(2π)^(3/2) * σx * σy * σz]exp(-((x-xc)^2)/(2*σx^2))exp(-((y-yc)^2)/(2*σy^2))[exp(-((z-H)^2)/(2σz^2)) exp(-((zH)^2)/(2σz^2))]Q是释放总量单位为毫克(xc,yc)是烟团中心在t时刻的水平位置由释放点位置与平均风场位移相加H是释放高度。最后一项是两个高斯指数的和物理意义是地面作为反射边界直接到达项和地表反射项叠加。计算人员暴露浓度时z通常取呼吸带高度1.2米而不取地面0米因为z0处反射项会让浓度翻倍高估地面附近人员的吸入量。σx、σy、σz是x、y、z三个方向的扩散参数随扩散时间增长。开阔平坦地形下扩散参数按Pasquill-Gifford曲线查表地铁站内部因为存在复杂的梁柱、屏蔽门和通风气流直接套PG系数会偏大工程上常见做法是取中性稳定度D级参数再乘以1.2~1.5的粗糙度修正。3.3 用numpy实现带风向和换气衰减的高斯烟团采样以下实现一个可独立调用的 PuffSource并把采样网格打到站厅平面。实现里额外加了一个空气交换衰减项模拟机械通风从环境中移除污染物的过程换气次数越高衰减越快。import numpy as np class PuffSource: def __init__(self, x0, y0, z0, Q, t0, wind_u0.3, wind_v0.1, a0.08, b0.9, c0.06, d0.9, H1.5): self.x0 x0 self.y0 y0 self.z0 z0 self.Q Q self.t0 t0 self.wind_u wind_u self.wind_v wind_v self.a, self.b, self.c, self.d a, b, c, d self.H H def sigma_x(self, dt): return self.a * dt ** self.b def sigma_y(self, dt): return self.a * dt ** self.b def sigma_z(self, dt): return self.c * dt ** self.d def concentration(self, x, y, z, t, air_exchange6.0): dt max(t - self.t0, 0.1) if dt 0: return 0.0 xc self.x0 self.wind_u * dt yc self.y0 self.wind_v * dt sx self.sigma_x(dt) sy self.sigma_y(dt) sz self.sigma_z(dt) decay np.exp(-air_exchange * dt / 3600.0) base self.Q * decay / ((2 * np.pi) ** 1.5 * sx * sy * sz) term1 np.exp(-((x - xc) ** 2) / (2 * sx * sx)) term2 np.exp(-((y - yc) ** 2) / (2 * sy * sy)) ref np.exp(-((z - self.H) ** 2) / (2 * sz * sz)) ref np.exp(-((z self.H) ** 2) / (2 * sz * sz)) return base * term1 * term2 * ref上面的代码把N个烟团批量计算时只需在外部循环里对这N个源的浓度求和。sigma_x 与 sigma_y 用了同一套系数这是开阔地形下的近似地铁站内若有纵向机械送风应让风速方向上的sigma系数随位移增长更快可以把sigma_x单独替换成随下风距离x变化的拟合式而不是随时间dt变化。扩散参数的取值直接影响浓度量级。以D级稳定度为例工程近似系数取 a0.08、b0.9、c0.06、d0.9在释放30秒、风速0.3 m/s的条件下烟团中心扩散半径约2.1米近源区浓度极高。若换气次数设为6次/小时30秒后衰减因子仍接近0.95说明短时间内机械通风对稀释的帮助很小只能靠时间堆。4. 把浓度场接进行人动力学耦合点、主循环与参数敏感性两个模型各自跑通不是难点真正的难点在耦合。浓度场数据以格点为离散载体行人动力学要求每个行人拿到自己所在位置的浓度值、浓度梯度以及这个浓度对他当前决策的影响。耦合设计的核心是定义一条浓度到行为参数的映射并且保证每个时间步内先更新浓度场再更新行人。4.1 耦合接口浓度如何量化为感知恐慌系数我在模拟系统里用的映射分三段浓度瞬时值 C_i 用来调整行人的避毒方向和速度衰减系数累计剂量 D_i 用来切换状态邻居恐慌比例 p 用来放大从众行为权重。速度衰减写成 v_eff v0 / (1 beta * C_i)beta是浓度敏感系数单位是立方米每毫克。浓度1 mg/m³、beta0.5时速度衰减到原速的67%这比用线性表达式平滑也避免浓度过高时速度变负。避毒梯度的计算用中心差分。二维网格上的C[i][j]对x的偏导近似为 (C[i][j1] - C[i][j-1]) / (2 * dx)对y的偏导同理。行人相对格点的位置用线性插值取浓度和梯度插值开销小且不产生网格尺度的锯齿状受力。上面这些接口要工作邻居恐慌比例还得靠行人之间的通信。每个行人每步把自己的状态广播给所在网格分桶内的其他行人恐慌比例就是桶内处于panic状态的人数除以桶内总人数。这个设计同时服务于两个目的从众方向计算和恐慌传播判定。4.2 逐帧事件循环和空间分桶主循环把两个子系统串在一个离散时间轴上。伪代码顺序如下for step in range(total_steps): t step * dt for puff in puff_list: puff.update(t) partition_pedestrians(grid_buckets, pedestrians) for p in pedestrians: C interpolate_concentration(grid_field, p.pos) grad_C interpolate_gradient(grid_field, p.pos) p.accumulate_dose(C, dt) p.update_state(C, p.dose) neighbors query_bucket(p.pos, radius3.0) contact compute_contact_force(p, neighbors) p.step(grad_C, neighbors_v, contact, dt) record_snapshot(t, pedestrians, grid_field)这段代码的先后顺序有讲究。puff.update 必须先于行人浓度查询执行保证同一时间步内行人读到的浓度场是当前时刻的。网格分桶在行人更新前做一次即可因为行人速度不超过3 m/s、dt不超过0.02s每步位移不超过6厘米不需要每步重建分桶索引可以两三个步长重建一次。接触力计算采用空间分桶剪枝桶尺寸取2倍行人感知半径这样每个行人只需检查相邻9个桶内的邻居能把O(N²)的邻居匹配降到接近O(N)。1000人场景下这一步是从“跑不动”到“能跑完”的关键优化。4.3 敏感参数表和调参顺序耦合模型参数多全部做网格搜索不现实。我实际调参时先固定行人动力学参数只扫浓度相关参数再固定浓度场只扫行为权重。下面的敏感度表记录的是“参数变化50%时最终疏散完成时间的变化幅度”参数含义敏感度调参建议换气次数 air_exchange通风稀释能力高现场实测优先别猜beta浓度对速度的抑制系数高需要毒理学数据支撑出口熟悉度比例知道出口方向的人数占比高用问卷或监控统计w_avoid避毒方向权重中0.3~0.5之间扫描接触力弹性系数k行人挤压形变中对总疏散时间影响小网格分辨率浓度场空间精度中过细只增加计算量注意换气次数设置为0时污染物只靠湍流扩散稀释近源区会持续保持高浓度疏散仿真里会出现近源行人快速倒地但这不代表真实车站情况只代表通风系统失效的极端场景。5. 把代码跑起来运行环境、验证三件事和常见脚本坑先把文本文件变成能运行的脚本。把代码保存为 metro_evac_csfm.py 时必须确认编码是UTF-8无BOMWindows记事本默认的带BOM UTF-8会让Python 3在解析首行时报 SyntaxError。命令行的运行方式很简单cd ./evac_sim python -m venv .venv source .venv/bin/activate # Windows下用 .venv\Scripts\activate pip install numpy scipy matplotlib numba python metro_evac_csfm.py --agents 800 --dt 0.02 --steps 12000 --scenario station_a这段命令把虚拟环境建在当前目录之后所有依赖都装进 .venv 而不会污染系统Python。numba在纯Python版本性能不足时可用来给受力计算加 jit(nopythonTrue)但要注意PuffSource里的字符串属性没法直接编译需要把烟团参数拆成numpy数组再传入。在WSL Ubuntu里写这段代码时终端字体建议换成Fira Code或JetBrains Mono等宽对齐和箭头符号渲染都很正常体验接近macOS下的终端VS Code里装上Pylance后numpy的ndarray类型提示会让acc、grad_c这类变量的维度错误在写代码阶段就暴露代码诊断的红波浪线比运行时追溯堆栈快得多。跑通之后做三件事的验证缺一不可。第一是守恒校验每个时间步检查“站内人数逃出人数中毒倒地人数”是否等于初始总人数误差超过1e-6说明力计算里有人重复更新或漏更新。第二是时间步收敛性把dt从0.02改成0.01重跑一遍疏散完成时间的差异应小于5%差太多说明数值积分在接触力较硬时已经不稳定。第三是剂量无关性把浓度采样网格加密一倍观察倒地人数曲线不应发生跳变。如果疏散曲线的中段斜率变化和换气次数对不上优先检查风场参数和通风衰减系数而不是去调行人腿力——浓度场的时间常数决定高峰拥堵何时消退行人参数只决定拥堵内的挤压细节。本文还有配套的精品资源点击获取
毒气泄漏场景下组合社会力模型与高斯烟团耦合疏散仿真
简介这是一份基于组合社会力模型CSFM与高斯烟团模型的地铁站毒气袭击紧急疏散模拟系统的Python复现与讲解文档面向具有一定编程基础的研究人员、工程师及应急疏散建模学者。资源包内为单个docx文档整包约60KB包含可运行的CSFM类实现、模型初始化、毒气源处理、社会力计算、气体浓度更新与可视化等完整代码并逐段附有解释和论文复现分析。读者可据此模拟毒气源位置与数量、管理人员反应速度、风速等因素对疏散时间及伤亡率的影响评估不同应急响应措施为地铁站设计和应急预案提供量化参考。文档还梳理了模型的数学基础、改进方向、验证方法及多场景应用潜力兼具科研对照与工程参考价值。目前已有80人学习下载适合对行人疏散与气体扩散仿真感兴趣的读者深入研读。1. 为什么突发毒气场景下普通社会力模型不够用某次我跑一个换乘站晚高峰疏散仿真时用经典Helbing社会力模型得出一个漂亮的结论所有人从两端站台向三个出口均匀撤离拥挤点稳定出现在闸机前。这个结果放在普通火灾里还能看但换成毒气泄漏真实监控里出现的画面完全不是这样——有人逆着人群往毒源方向走有人贴着墙根蹲下也有人在出口已经畅通时突然折返去扶同伴。普通社会力模型只描述“我想去出口其他人挡我”它假设每个行人在整个仿真里决策恒定。毒气场景恰恰破坏了这个假设浓度改变人的移动能力剂量主导人的决策。CSFM组合社会力模型把“毒气浓度影响行为权重”做成受力方程的一部分和高斯烟团给的浓度场耦合后才能复现这些反直觉行为。本文讲的这套模拟系统核心就是两个模型怎么各自实现、怎么在一个事件循环里互相喂数据以及参数调到什么量级结果才不骗人。2. 组合社会力模型 CSFM与经典社会力模型的差异落点CSFM不是一个能追溯到某篇固定论文的统一公式工程实现上一般是在Helbing提出的社会力模型骨架上做三件组合把期望速度的方向从“固定朝出口”改成多方向加权把恐慌系数做成随毒剂量变化的动态量把行人状态从单一行走扩成包含中毒倒地、折返救援的状态机。2.1 从Helbing方程出发CSFM在哪些项上动手经典社会力模型按牛顿体系展开单个行人的加速度由三部分叠加m_i * dv_i/dt m_i * (v_i^0 * e_i - v_i) / τ Σ_j f_ij Σ_w f_iw等号右边第一项是自驱动力描述行人希望以期望速度 v_i^0 朝方向 e_i 走但受身体反应时间 τ 限制实际速度 v_i 会滞后。第二项是行人与行人之间的社会排斥力第三项是墙壁等障碍物对行人的排斥力。行人间的物理接触力挤压与摩擦在高密度下再叠加进 f_ij。CSFM在这套体系上做改动时首先要动的就是 e_i。经典模型里 e_i 是一个常量方向即从当前位置指向出口的连线方向CSFM里它被重写为出口方向、避毒梯度方向、从众方向三个分量的归一化加权和。第二个大改动是 τ反应时间在恐慌状态下会骤降行人动作变得“又急又僵硬”计算公式常写成 τ τ0 * exp(-η * D_i)D_i是累计毒剂量。第三个改动是引入状态机毒剂量超过阈值后行人进入中毒减速或倒地状态这是一个离散事件不能靠连续力场平滑表达。2.2 状态机参数和“行为权重”怎么放进受力方程我一般把每个行人设置成四个状态normal、panic、intoxicated、down。normal到panic的切换条件是检测到浓度或他人恐慌比例超阈值panic到intoxicated的切换条件是Haber法则下的累计剂量超阈值。Haber法则指毒性效应近似等于浓度与暴露时间的乘积写成 D_i(t) ∫ C_i(t) dt这里的 C_i(t) 是行人当前位置的毒气浓度。行为权重的优先级也在这个状态机上定义panic状态下避毒梯度权重最高normal状态下出口方向权重最高intoxicated状态所有速度方向的合成会叠加巨大噪声。从众方向不是全局随机而是取以该行人为圆心、半径3米内所有邻居平均速度的方向这个方向存在局部信息传播的语义也是CSFM与仅靠浓度场驱动力模型的最大区别。2.3 最小Pedestrian实现下面是一个可直接用于仿真核心的最小Python实现省略了碰撞邻居搜索的优化但保留了状态切换和期望方向合成的完整逻辑。import numpy as np class Pedestrian: def __init__(self, pos, goal, v01.2, tau00.5, radius0.2, mass80.0): self.pos np.asarray(pos, dtypefloat) self.goal np.asarray(goal, dtypefloat) self.v np.zeros(2) self.v0 v0 self.tau0 tau0 self.radius radius self.mass mass self.rho 1.0 # 恐慌系数初始为1 self.dose 0.0 # 累计毒剂量 self.state normal def set_panicky(self, factor): self.rho max(1.0, factor) def compute_e(self, grad_c, neighbors_v, w_exit1.0, w_avoid0.4, w_follow0.2): e_exit (self.goal - self.pos) e_exit / np.linalg.norm(e_exit) 1e-9 e_avoid np.zeros(2) if np.linalg.norm(grad_c) 1e-6: e_avoid -grad_c / (np.linalg.norm(grad_c) 1e-9) e_follow np.zeros(2) if len(neighbors_v) 0: e_follow np.mean(neighbors_v, axis0) n np.linalg.norm(e_follow) if n 1e-6: e_follow / n e w_exit * e_exit w_avoid * e_avoid w_follow * e_follow n np.linalg.norm(e) return e / (n 1e-9) def step(self, grad_c, neighbors_v, contact_force, dt): if self.state down: self.v[:] 0.0 return tau self.tau0 * np.exp(-0.15 * (self.rho - 1.0)) e self.compute_e(grad_c, neighbors_v) v_desired self.v0 * e if self.state intoxicated: v_desired * 0.4 acc (v_desired - self.v) / tau contact_force / self.mass self.v acc * dt self.pos self.v * dt代码里 compute_e 方法把三个方向加权后归一化然后自驱动力按社会力模型的标准形式合成。注意 panic 状态和 intoxicated 状态分开处理panic 时反应时间 tau 变小行人动作更急intoxicated 时直接把期望速度打折模拟中毒后肌肉乏力。这个顺序不能反。引入接触力时我用的是Helbing原始参数量级的推荐表单位统一为米、千克、秒参数含义典型值校准依据v0期望速度1.2 m/s正常/ 1.8 m/s恐慌站台通道步行实测τ反应时间0.5 s正常/ 0.3 s恐慌Helbing标定区间A社会排斥力强度2000 N原始模型推荐量级B排斥力作用距离0.08 m原始模型推荐量级k身体弹性系数1.2e5 kg/s²接触力线性项κ摩擦系数2.4e5 kg/(m·s)接触力切向项数值更新建议使用Heun方法即先预估再校正显式欧拉在 τ 较小、Δt0.02s 时容易产生速度振荡将时间步控制在 0.01~0.02 秒一秒仿真推进 50~100 步1000人的场景在纯Python下仍能接近实时。3. 高斯烟团模型地铁站毒气浓度场的计算骨架浓度场是这套系统里另一个独立零件。疏散仿真中最常用的是高斯烟羽模型和高斯烟团模型两者都基于湍流扩散的正态分布假设。烟羽模型描述连续泄漏源烟团模型描述瞬时泄漏源。地铁站内毒气袭击更接近突发瞬时释放而且机械通风会让风向在几个通风周期内变化用多个烟团叠加以时间轴推进比拟合一个稳定烟羽更自然。3.1 为什么选烟团而不是烟羽烟羽模型假设源强恒定、气象条件平稳输出一个不随时间变化的浓度场。应急疏散场景里最关心的恰恰是“释放后第30秒站台两端浓度分别多少”这要求浓度场是时间函数。烟团的数学形式自带时间维每个释放脉冲表达为一个以释放点为中心、随风飘移且不断膨胀的三维高斯钟多个脉冲叠加就是非平稳浓度场。另一个原因是地铁站内部结构对风向的扰动太大。烟羽模型要求给出稳定的风速和风向而站厅内气流由通风系统主导风速通常低于0.5 m/s风向不定。烟团模型把空气运动拆成平均风场位移和湍流扩散两个独立部分机械通风导致的空气交换可以等效成一个额外的衰减系数不需要精确刻画具体流场。3.2 浓度公式和地面反射项解释单次瞬时释放的高斯烟团浓度表达式为C(x,y,z,t) Q / [(2π)^(3/2) * σx * σy * σz]exp(-((x-xc)^2)/(2*σx^2))exp(-((y-yc)^2)/(2*σy^2))[exp(-((z-H)^2)/(2σz^2)) exp(-((zH)^2)/(2σz^2))]Q是释放总量单位为毫克(xc,yc)是烟团中心在t时刻的水平位置由释放点位置与平均风场位移相加H是释放高度。最后一项是两个高斯指数的和物理意义是地面作为反射边界直接到达项和地表反射项叠加。计算人员暴露浓度时z通常取呼吸带高度1.2米而不取地面0米因为z0处反射项会让浓度翻倍高估地面附近人员的吸入量。σx、σy、σz是x、y、z三个方向的扩散参数随扩散时间增长。开阔平坦地形下扩散参数按Pasquill-Gifford曲线查表地铁站内部因为存在复杂的梁柱、屏蔽门和通风气流直接套PG系数会偏大工程上常见做法是取中性稳定度D级参数再乘以1.2~1.5的粗糙度修正。3.3 用numpy实现带风向和换气衰减的高斯烟团采样以下实现一个可独立调用的 PuffSource并把采样网格打到站厅平面。实现里额外加了一个空气交换衰减项模拟机械通风从环境中移除污染物的过程换气次数越高衰减越快。import numpy as np class PuffSource: def __init__(self, x0, y0, z0, Q, t0, wind_u0.3, wind_v0.1, a0.08, b0.9, c0.06, d0.9, H1.5): self.x0 x0 self.y0 y0 self.z0 z0 self.Q Q self.t0 t0 self.wind_u wind_u self.wind_v wind_v self.a, self.b, self.c, self.d a, b, c, d self.H H def sigma_x(self, dt): return self.a * dt ** self.b def sigma_y(self, dt): return self.a * dt ** self.b def sigma_z(self, dt): return self.c * dt ** self.d def concentration(self, x, y, z, t, air_exchange6.0): dt max(t - self.t0, 0.1) if dt 0: return 0.0 xc self.x0 self.wind_u * dt yc self.y0 self.wind_v * dt sx self.sigma_x(dt) sy self.sigma_y(dt) sz self.sigma_z(dt) decay np.exp(-air_exchange * dt / 3600.0) base self.Q * decay / ((2 * np.pi) ** 1.5 * sx * sy * sz) term1 np.exp(-((x - xc) ** 2) / (2 * sx * sx)) term2 np.exp(-((y - yc) ** 2) / (2 * sy * sy)) ref np.exp(-((z - self.H) ** 2) / (2 * sz * sz)) ref np.exp(-((z self.H) ** 2) / (2 * sz * sz)) return base * term1 * term2 * ref上面的代码把N个烟团批量计算时只需在外部循环里对这N个源的浓度求和。sigma_x 与 sigma_y 用了同一套系数这是开阔地形下的近似地铁站内若有纵向机械送风应让风速方向上的sigma系数随位移增长更快可以把sigma_x单独替换成随下风距离x变化的拟合式而不是随时间dt变化。扩散参数的取值直接影响浓度量级。以D级稳定度为例工程近似系数取 a0.08、b0.9、c0.06、d0.9在释放30秒、风速0.3 m/s的条件下烟团中心扩散半径约2.1米近源区浓度极高。若换气次数设为6次/小时30秒后衰减因子仍接近0.95说明短时间内机械通风对稀释的帮助很小只能靠时间堆。4. 把浓度场接进行人动力学耦合点、主循环与参数敏感性两个模型各自跑通不是难点真正的难点在耦合。浓度场数据以格点为离散载体行人动力学要求每个行人拿到自己所在位置的浓度值、浓度梯度以及这个浓度对他当前决策的影响。耦合设计的核心是定义一条浓度到行为参数的映射并且保证每个时间步内先更新浓度场再更新行人。4.1 耦合接口浓度如何量化为感知恐慌系数我在模拟系统里用的映射分三段浓度瞬时值 C_i 用来调整行人的避毒方向和速度衰减系数累计剂量 D_i 用来切换状态邻居恐慌比例 p 用来放大从众行为权重。速度衰减写成 v_eff v0 / (1 beta * C_i)beta是浓度敏感系数单位是立方米每毫克。浓度1 mg/m³、beta0.5时速度衰减到原速的67%这比用线性表达式平滑也避免浓度过高时速度变负。避毒梯度的计算用中心差分。二维网格上的C[i][j]对x的偏导近似为 (C[i][j1] - C[i][j-1]) / (2 * dx)对y的偏导同理。行人相对格点的位置用线性插值取浓度和梯度插值开销小且不产生网格尺度的锯齿状受力。上面这些接口要工作邻居恐慌比例还得靠行人之间的通信。每个行人每步把自己的状态广播给所在网格分桶内的其他行人恐慌比例就是桶内处于panic状态的人数除以桶内总人数。这个设计同时服务于两个目的从众方向计算和恐慌传播判定。4.2 逐帧事件循环和空间分桶主循环把两个子系统串在一个离散时间轴上。伪代码顺序如下for step in range(total_steps): t step * dt for puff in puff_list: puff.update(t) partition_pedestrians(grid_buckets, pedestrians) for p in pedestrians: C interpolate_concentration(grid_field, p.pos) grad_C interpolate_gradient(grid_field, p.pos) p.accumulate_dose(C, dt) p.update_state(C, p.dose) neighbors query_bucket(p.pos, radius3.0) contact compute_contact_force(p, neighbors) p.step(grad_C, neighbors_v, contact, dt) record_snapshot(t, pedestrians, grid_field)这段代码的先后顺序有讲究。puff.update 必须先于行人浓度查询执行保证同一时间步内行人读到的浓度场是当前时刻的。网格分桶在行人更新前做一次即可因为行人速度不超过3 m/s、dt不超过0.02s每步位移不超过6厘米不需要每步重建分桶索引可以两三个步长重建一次。接触力计算采用空间分桶剪枝桶尺寸取2倍行人感知半径这样每个行人只需检查相邻9个桶内的邻居能把O(N²)的邻居匹配降到接近O(N)。1000人场景下这一步是从“跑不动”到“能跑完”的关键优化。4.3 敏感参数表和调参顺序耦合模型参数多全部做网格搜索不现实。我实际调参时先固定行人动力学参数只扫浓度相关参数再固定浓度场只扫行为权重。下面的敏感度表记录的是“参数变化50%时最终疏散完成时间的变化幅度”参数含义敏感度调参建议换气次数 air_exchange通风稀释能力高现场实测优先别猜beta浓度对速度的抑制系数高需要毒理学数据支撑出口熟悉度比例知道出口方向的人数占比高用问卷或监控统计w_avoid避毒方向权重中0.3~0.5之间扫描接触力弹性系数k行人挤压形变中对总疏散时间影响小网格分辨率浓度场空间精度中过细只增加计算量注意换气次数设置为0时污染物只靠湍流扩散稀释近源区会持续保持高浓度疏散仿真里会出现近源行人快速倒地但这不代表真实车站情况只代表通风系统失效的极端场景。5. 把代码跑起来运行环境、验证三件事和常见脚本坑先把文本文件变成能运行的脚本。把代码保存为 metro_evac_csfm.py 时必须确认编码是UTF-8无BOMWindows记事本默认的带BOM UTF-8会让Python 3在解析首行时报 SyntaxError。命令行的运行方式很简单cd ./evac_sim python -m venv .venv source .venv/bin/activate # Windows下用 .venv\Scripts\activate pip install numpy scipy matplotlib numba python metro_evac_csfm.py --agents 800 --dt 0.02 --steps 12000 --scenario station_a这段命令把虚拟环境建在当前目录之后所有依赖都装进 .venv 而不会污染系统Python。numba在纯Python版本性能不足时可用来给受力计算加 jit(nopythonTrue)但要注意PuffSource里的字符串属性没法直接编译需要把烟团参数拆成numpy数组再传入。在WSL Ubuntu里写这段代码时终端字体建议换成Fira Code或JetBrains Mono等宽对齐和箭头符号渲染都很正常体验接近macOS下的终端VS Code里装上Pylance后numpy的ndarray类型提示会让acc、grad_c这类变量的维度错误在写代码阶段就暴露代码诊断的红波浪线比运行时追溯堆栈快得多。跑通之后做三件事的验证缺一不可。第一是守恒校验每个时间步检查“站内人数逃出人数中毒倒地人数”是否等于初始总人数误差超过1e-6说明力计算里有人重复更新或漏更新。第二是时间步收敛性把dt从0.02改成0.01重跑一遍疏散完成时间的差异应小于5%差太多说明数值积分在接触力较硬时已经不稳定。第三是剂量无关性把浓度采样网格加密一倍观察倒地人数曲线不应发生跳变。如果疏散曲线的中段斜率变化和换气次数对不上优先检查风场参数和通风衰减系数而不是去调行人腿力——浓度场的时间常数决定高峰拥堵何时消退行人参数只决定拥堵内的挤压细节。本文还有配套的精品资源点击获取
