Python复现CSFM与高斯烟团：地铁毒气疏散模拟实战

简介一份基于组合社会力模型CSFM与高斯烟团模型的地铁站毒气袭击紧急疏散模拟系统文档面向具备一定编程基础的研究人员、工程师及对应急疏散模拟感兴趣的学者。文档完整复现了论文核心模型从模型初始化、毒气源处理、社会力计算到动态更新与可视化等环节均给出可运行的Python代码及逐步解释并深入探讨毒气源位置与数量、管理响应速度、风速等因素对乘客疏散和伤亡情况的影响。压缩包内含1个docx文档包体约60KB内容紧凑便于查阅。目前已有81人学习浏览。除代码实现外文档还系统整理了模型的数学基础、改进方向、验证方法及多场景应用潜力直击毒气扩散与人员疏散耦合建模的核心难点适合复现实验与二次开发可帮助读者全面掌握CSFM工作原理并据此评估不同应急响应措施效果为地铁站设计和应急预案制定提供科学依据。1. 毒气袭击下的地铁疏散为什么值得把 CSFM 复现一遍论文里有个反直觉结论地铁站毒气袭击时风速越大受伤乘客反而越少。这不是拍脑袋推出来的而是组合社会力模型CSFM加高斯烟团模型跑出来的结果——风把毒气团吹散局部浓度下降乘客暴露量跟着降。这份资源把整套逻辑用 Python 完整复现成了紧急疏散模拟系统毒气源扩散、个体竞争、群体跟随、伤亡估算全在一个可视化的二维地铁站里跑起来。适合做应急疏散研究、地铁站预案评估的工程师和研究生。它不解决所有问题但能让你在平面图纸之外多一个量化推演工具。2. 先看懂四种力再动参数CSFM 的力学拆解与高斯烟团耦合2.1 为什么是社会力模型从元胞自动机到组合模型的选型逻辑先摆出论文里的模型对比结论元胞自动机计算效率高规则也简单但个体行为过度简化Agent-based 模型个体异质性表现好参数设置却复杂到让人头大社会力模型物理意义明确行为真实性最好被论文认定为模拟行人行为的最佳模型。CSFM 的做法是在经典社会力框架上把毒气影响当成一种额外力叠加进去所以它既能保留社会力模型那种“人能感受到推搡、避让、恐慌”的真实感又能回答“毒气环境下多少人会受伤”这种传统行人流模型答不了的问题。从我的复现经验看选社会力模型还有一层现实原因它的力是连续向量符合直觉调试时能直接看某个人的受力分量不像元胞自动机那样只能看到 0/1 状态迁移。对要做应急决策支撑的人来说模型里每个力都有物理含义汇报时也说得清。论文给的对比表信息量不小我按复现时的关注点整理成下面这张表模型类型优点局限性适用场景元胞自动机计算效率高规则简单个体行为过于简化大规模人群简单行为模拟社会力模型物理意义明确行为真实计算复杂度较高中规模人群复杂行为模拟Agent-based个体异质性表现好参数设置复杂需要个体差异化行为的场景组合模型综合多种模型优势模型整合难度大复杂多因素应急场景这张表决定了后文所有参数调试的方向既然选了社会力模型你的工作量就集中在力的系数上而不是在格子和规则上。2.2 四种力的数学表达与代码映射CombinedSocialForceModel 里的 social_force 方法把每个人的受力拆成四部分注释里写得清楚但代码细节值得再抠一遍。首先是目标驱动力nearest_exit min(self.exits, keylambda x: np.linalg.norm(pos_i - x)) e_i (nearest_exit - pos_i) / (np.linalg.norm(nearest_exit - pos_i) 1e-5) driving_force (self.desired_speed * e_i - vel_i) / self.tau这段的逻辑是拿“期望速度方向 e_i”和“当前速度 vel_i”的差除以反应时间 tau。1e-5 是防除零的当人站在出口上时 norm 为 0直接除会崩。desired_speed 默认 1.34 m/s这个值来自行人流文献里自由行走速度的常见值你要是换场景改成 1.0 会更保守。然后是人与人之间的排斥力r_ij pos_i - pos_j d_ij np.linalg.norm(r_ij) if d_ij 0: n_ij r_ij / d_ij people_force self.A * np.exp(-d_ij / self.B) * n_ijA2000 表示接触时的排斥强度B0.08 是力的作用范围。这个形式是经典社会力模型里的指数衰减距离超过 B 的 3 倍左右约 0.24 米力就接近零了所以它本质上是近距避让不是远程规划。障碍物力在代码里被简化成边界力离墙 1 米以内用 k 除以距离近似推离。k1.2e5 默认值对标 Helbing 当年论文里的量级实际跑的时候你会发现边界力只在小范围内起作用所以行人不会在场地中间感受到墙。最后一股力是气体影响它直接取浓度网格的梯度然后乘 10 的系数往外推gas_force -np.array([gas_grad_x, gas_grad_y]) * 10这里系数 10 是手动调的论文没给标准值。我一般先让 gas_force 的量级和 driving_force 可比否则要么气体影响不痛不痒要么行人直接原地打转。这个系数在后面的避坑章里会专门展开。2.3 高斯烟团模型怎么叠进社会力框架毒气浓度用的是高斯烟团模型update_gas_concentration 方法遍历所有毒气源对每个源生成二维高斯分布再叠加风向偏移。关键代码是dx - self.wind_direction[0] * self.wind_speed dy - self.wind_direction[1] * self.wind_speed distance_sq dx**2 dy**2 concentration source[strength] * np.exp(-distance_sq / (2 * self.gas_diffusion_rate**2))wind_speed 和 wind_direction 相当于把高斯中心沿风向平移gas_diffusion_rate 控制烟团半径。这里有个容易误读的点代码里 wind_speed 是直接加到 dx、dy 上的单位其实是“每帧偏移量”不是 m/s。所以论文里“风速越大受伤越少”的结论在复现时是通过这个偏移量体现的你调 wind_speed 时要结合 dt 看实际移动距离。浓度最后被归一化到 [0,1]范围是整个 50×50 网格后面伤亡判定用的 toxicity_threshold0.7 就是相对于归一化浓度说的。这一章把力的来源讲清楚了下一章直接落地跑代码。3. 把论文跑成动画CombinedSocialForceModel 的初始化、更新与可视化3.1 环境与最小运行骨架我用的环境是 Python 3.9依赖只有 numpy、matplotlib 和 scipy.stats其中 scipy 其实只用来导入了 norm核心计算没用到它。跑最小示例只要几行model CombinedSocialForceModel(num_agents50, width40, height40, num_exits2) model.add_gas_source(position[10, 20], strength1.0) model.wind_speed 0.5 model.wind_direction np.array([1, 0.5]) model.visualize(steps200)这段对应论文里“两个毒气源、东北风”的基准实验。num_agents 我建议从 50 开始别一上来就 1000因为 social_force 里有个双重循环人数从 50 涨到 200单步耗时涨的是平方级我刚开始用 500 人跑 200 帧等了十分钟没出动画差点以为是死循环。提示第一次跑先把 num_agents 降到 30确认动画能出再往上加不然双重循环会卡到你怀疑人生。visualize 里用的是 FuncAnimation每帧调用 self.update() 再刷新散点图和热图。blitTrue 是性能关键只有变化的绘图对象会被重绘否则每帧重画整幅图动画会卡成 PPT。3.2 初始化参数出口分布、毒气源、行人位置构造函数里有个细节值得单独说出口位置是均匀分布在左右边界上的不是集中在一侧。代码逻辑是偶数索引出口放在左边界、奇数放在右边界所以 num_exits2 时一个在左边中间一个在右边中间。这个布局模拟的是地铁站两侧都有通道口的常见情况。你要是想改出口位置直接在 exits 数组上赋值就行。行人初始位置是 np.random.rand(num_agents, 2) 乘站厅尺寸完全随机。这会让动画开头出现人挤在毒气源旁边的情况不影响结果但会拉长疏散时间。常见做法是在远离毒气源的区域采样或者像论文那样按入口排队我一般会改成positions np.random.rand(num_agents, 2) * np.array([width, height]) positions[:, 0] np.clip(positions[:, 0], 15, width - 5)把 x 方向限制在场地中后部让前排行人更有“从两端逃向出口”的观感也更接近真实候车人群分布。3.3 主循环与可视化从 update 到 FuncAnimationupdate 方法每帧做三件事更新气体浓度、逐个计算社会力并更新速度位置、做边界裁剪。速度上限写在循环里speed np.linalg.norm(self.velocities[i]) if speed self.desired_speed * 1.5: self.velocities[i] self.velocities[i] / speed * self.desired_speed * 1.5这个 1.5 倍上限是为了防止行人被恐慌加速推到不合理的速度。我在调参时发现如果不加这个限制gas_force 调大后行人会一帧飞出半个站厅动画看起来像瞬移。desired_speed * 1.5 大约是 2 m/s对应跑步速度符合常识。可视化的核心是 update 回调函数def update(frame): self.update() scat.set_offsets(self.positions) if len(self.gas_sources) 0: self.update_gas_concentration() gas_img.set_array(self.gas_concentration.T) return scat, gas_img if len(self.gas_sources) 0 else scat注意 set_array 之前要重新调用 update_gas_concentration因为浓度场每帧都在被风吹着移。gas_img 用 imshow 叠加在散点图下层alpha0.3 保证红热图不遮住行人点。3.4 关键参数对照表参数代码位置默认值作用调参建议A / Bsocial_force2000 / 0.08人际排斥强度与范围A 过大容易人挤人振荡k / kappasocial_force1.2e5 / 2.4e5边界排斥与摩擦kappa 在本实现中未实际调用tauinit0.5反应时间越小行人响应越快也越容易振荡desired_speedinit1.34期望速度恐慌场景可提到 1.8gas_diffusion_rateinit0.1烟团扩散半径决定毒气影响范围wind_speedinit0.0风平移量论文结论关键变量注意单位是帧偏移toxicity_thresholdEnhancedCSFM0.7伤亡判定阈值需配合浓度归一化理解这张表是我按复现顺序整理的kappa 在基础版里定义了却没在 update 里用增强版也没用到属于论文代码里的“预留参数”。遇到这种参数心里有数就行别浪费时间调它。4. 增强版 CSFM伤亡判定、恐慌传染与统计输出怎么落地4.1 健康状态与毒气暴露伤害基础版只算出每个人的位置和气体浓度没回答“多少人受伤”。EnhancedCSFM 用 health_status 数组1 健康、0 受伤/死亡和 exposure_time 累计暴露帧数把毒气伤害拆成三项加权damage ( self.gas_impact_params[concentration_weight] * concentration self.gas_impact_params[exposure_weight] * (self.exposure_time[i] / 100) self.gas_impact_params[fitness_weight] * np.random.normal(0.5, 0.2) ) individual_resistance np.random.normal(1.0, 0.2) damage / individual_resistance if damage self.toxicity_threshold: self.health_status[i] 0 self.velocities[i] * 0.3三项权重分别是浓度 0.8、暴露时长折算 0.5、个体体质随机项 0.3最后除以一个随机抵抗力。这套设计比“浓度过线就受伤”要细短时间高浓度和长时间低浓度都能造成伤害只是路径不同。individual_resistance 用正态分布模拟个体差异np.random.normal(1.0, 0.2) 表示大多数人抵抗力在 0.8 到 1.2 之间少数人特别敏感。受伤后速度直接乘以 0.3对应论文里“受伤者行动受限”的假设。有个实现细节要注意浓度检查只在 concentration 0.1 时才累计暴露时间这是为了防止归一化后接近 0 的数值噪声污染统计。毒性阈值 0.7 是相对值你如果把气体源 strength 调大了一倍归一化后阈值对应的实际物理浓度会变所以每次改 strength 都要重新标定阈值这是我踩过的坑后面避坑章会展开。4.2 恐慌水平的社会传染与速度耦合恐慌在增强版里是一个独立的状态量叠加在基础社会力之上。update_panic_level 做了三件事浓度驱动的基础恐慌上涨、邻域恐慌传染、恐慌影响速度。核心代码for j in range(self.num_agents): if i j: continue if (np.linalg.norm(self.positions[i] - self.positions[j]) 5 and self.panic_level[j] self.panic_level[i]): self.panic_level[i] self.panic_spread_rate * (self.panic_level[j] - self.panic_level[i])距离 5 米以内的恐慌者会把恐慌水平传染给邻居传染系数 panic_spread_rate0.05。这里 5 米是硬编码的感知半径对应人在站厅里能注意到周围五米内有人奔跑的自然反应。恐慌对速度的影响是概率性的70% 概率加速到 1panic*0.5 倍30% 概率因为混乱而减速。这个设计比“恐慌加速”更真实因为高恐慌下人群容易发生堵塞和推搡局部有人反而跑不起来。更新顺序也值得注意EnhancedCSFM.update 先更新浓度再更新健康再更新恐慌最后 super().update(dt) 走基础社会力。如果你把顺序反过来伤亡判断用的浓度就是上一帧的旧值第一帧会漏判。panic_level 大于 0.7 时还会以 10% 概率给速度加一个随机旋转模拟恐慌者的无序转向——这个随机项让动画看起来更“乱”但注意它可能让行人绕远路调参时别把概率调太高。4.3 统计输出与实验对比get_stats 在原文里只给出了统计项定义输出格式是我补的def get_stats(self): healthy np.sum(self.health_status 0) injured np.sum(self.health_status 0) avg_panic np.mean(self.panic_level[self.health_status 0]) return { healthy: healthy, injured: injured, avg_panic: avg_panic, total_agents: self.num_agents }用这个接口做对照实验很方便跑三次模拟分别统计两个毒气源、一个毒气源、无风速三种配置下的受伤人数就能复现论文的结论。我一般会把每次模拟的随机种子固定下来np.random.seed(42) 放在创建模型前否则两次实验的初始位置不同统计出来的差异会混入随机噪声。受伤害度里又用了一次 np.random.normal所以光固定种子还不够建议把 individual_resistance 的生成也固定或者把多次运行的均值当结果。统计输出可以配合可视化里的热图一起看热图变红的位置对应高浓度区散点消失的位置对应伤亡者。用这两者对照能直观确认伤亡是不是集中在毒气扩散路径上。5. 参数调试避坑阈值、速度上限与边界力最容易翻车5.1 四个高频踩坑记录先说结论这套代码从头跑到尾不难难的是改完参数后结果还能自圆其说。以下四个坑我都实际翻过车按频率排序。坑一毒性阈值跟着浓度归一化一起漂移现象改了毒气源 strength 从 1.0 改成 3.0其他参数没动结果受伤人数反而从 40 人降到 10 人。 原因update_gas_concentration 最后做了 max 归一化strength 变大后整个浓度场被等比放大了但 toxicity_threshold 还是 0.7归一化后的相对阈值没变实际伤的人反而少了。 解决把 toxicity_threshold 改成相对值理解或者在 update_gas_concentration 里去掉归一化直接用绝对浓度判伤。我一般保留归一化但每次改 strength 后就跑一次单帧浓度统计把最大浓度和平均浓度打出来再按 0.7 的百分比重新设定阈值。坑二gas_force 系数拍脑袋设 10行人原地鬼打墙现象把气体影响力系数从 10 调到 50行人开始在高浓度区附近绕圈不往出口走了。 原因gas_force 的方向是浓度梯度反方向系数调大后这股的力超过了目标驱动力行人被推离毒气源后又被期望速度拉回去形成振荡。 解决别单独调 gas_force 系数。我一般先跑一帧把 driving_force 和 people_force 的量级打印出来再让 gas_force 的最大值控制在 driving_force 的 30%~50% 范围。这套代码里系数 10 配合 tau0.5 是个稳定组合改一个就得配套改另一个。坑三风速直接平移高斯中心风太大时浓度场“跳帧”现象wind_speed 从 0.5 调到 5.0热图开始一帧一帧地跳不再平滑。 原因浓度更新是每帧直接把 dx 减 wind_speed风速 5 相当于每帧搬了 5 个网格高斯烟团在离散网格上显示为跳跃移动。 解决把风速按帧折算wind_velocity_per_frame wind_speed * dt然后让偏移量不超过两个网格。论文里的风速是本构场风速不是代码里的帧偏移量所以做敏感性分析时建议把单位换算关系写进注释免得后面对不上。坑四边界力只在 1 米内生效行人会贴在墙上滑动现象行人走到墙边后不反弹反而沿着墙一路滑到出口。 原因social_force 里墙力只在 pos_i[0] 1 或 width-1 时触发1 米以外没有任何导向障碍物的力。 解决如果你要模拟立柱、闸机等内部障碍物光有边界力不够。常见做法是把障碍物也建模成 repulsion 源给每个障碍物一个位置和排斥半径在 social_force 里叠加一层。论文里没提障碍物复现时可以不动。5.2 调试时我一般盯哪几个指标除了直接看动画我每轮调参都会输出三组数字存活人数曲线健康人数随帧数下降的斜率、平均恐慌水平、每帧平均速度。存活曲线陡降说明毒性阈值太低平均速度持续大于 2 m/s 说明速度上限失效平均恐慌涨到 0.9 说明恐慌扩散太快得调小 panic_spread_rate。这三个指标互相印证比单看受伤总数更能定位问题。还有一个容易被忽略的细节simulate 和 visualize 里的 steps 不是同一个概念。visualize(steps200) 内部每帧都会调用 update而 update 里 dt0.05 是物理步长所以 200 帧对应的模拟时长是 10 秒。你如果要模拟 5 分钟的疏散过程steps 得设 6000这时双重循环的 O(n²) 会成为瓶颈建议先把 num_agents 降到 100 以内。6. 给模型做体检疏散率曲线与参数扫描的三种做法验证模型不只是看动画顺不顺眼得拿量化曲线说话。我常用的三个做法按成本从低到高排。第一个是固定种子三连跑。np.random.seed(42) 放在模型创建前同一配置跑三次取均值受伤人数、平均恐慌都按均值报。这是最低成本的消噪手段。第二个是画疏散率曲线。把健康人数比例按帧记录下来画成横轴为模拟时间、纵轴为剩余健康人数的折线。这个曲线一眼能看出毒气源放在哪里伤亡来得最快。代码很直接health_ratio [] for _ in range(steps): model.update() stats model.get_stats() health_ratio.append(stats[healthy] / stats[total_agents])第三个是参数扫描。我一般扫风速和毒气源数量两个变量风速从 0 到 2 步长 0.5每个风速下记录最终受伤人数。扫完画一张柱状图论文里“风速越大受伤越少”的结论就直接可视化出来了。扫描时有个细节风速每次改变风位移量要按 dt 折算成帧偏移同时把 gas_diffusion_rate 保持 0.1 不变这样对照才是单一变量。参数扫描看起来笨但它是把黑匣子模型变成可解释结果的必经之路。从那以后我每次复现论文模型都会强制走一遍“固定种子→跑曲线→扫单参数”这个流程不扫参数就不敢拿结论去汇报。希望帮到你。本文还有配套的精品资源点击获取