大气扩散模型是一类研究污染物在大气中扩散、运输、转化和沉降等过程的模型, 对解决工业发展带来的污染物排放及危险品泄漏等问题具有重要意义。大气扩散模型通常是独立于 GIS 而在环境科学领域发展起来的, 在环境评价领域的应用较为广泛, 对应用于应急管理中危险气体扩散模拟及辅助决策模型的研究则相对较少, 因此对地形影响、三维扩散及计算速度等应急响应的特殊需求尚缺乏深入的研究。
近年来, 随着我国经济建设的快速发展, 对油气资源的需求越来越大, 伴随着资源开采过程而来的危险事故和环境污染时有发生。我国含硫天然气(体积浓度>1%)的储量占全国天然气储量的 1/4,高含硫气田的开发技术难度大, 风险高, 剧毒重气体含量高, 一旦发生泄露事故, 就有可能带来巨大的人员伤亡和财产损失。因此, 将大气扩散模型应用于毒害性气体泄漏事故的应急管理之中十分重要, 能够为事故的应急预案及灾后处理提供理论依据和决策支持。
目前已有学者对基于大气扩散模型的应急管理系统进行研究和开发, 不同研究中根据具体问题的需求选用了不同的大气扩散模型。按基本原理扩散模型可以分为三大类: 数值模拟模型、高斯模型和拉格朗日模型。
数值模拟模型基于流体力学中纳维尔–斯托克斯方程的大气扩散偏微分方程, 常见的模型包括CMAQ, GEOS-Chem 和 Fluent 等 。 Feng 等 [1]使 用CFD 模型建立地下空间有毒气体泄漏风险评估模型, 并应用于广州某地铁中转站。陈建国等[2]对开县井喷事故进行数值模拟, 成功地模拟了有害气体输运的积聚和导向等现象; 练章华等[3]用 Fluent 软件建立能反映实际硫化氢扩散过程的三维重气流动和扩散的模型, 用 CFD 软件模拟开县事故源周边的污染浓度变化。但是, 数值模拟方法计算精细、耗费大、计算时间较长, 不适用于应急情况下的实时计算。
高斯模型的基本原理是将大气扩散过程遵循的基本偏微分方程——质量守恒方程在合理假设的基础上推导转化成一个简单的数学公式, 称为高斯烟羽分布公式, 常见的模型包括 ADMS[4], AERMOD[5]和 ALOHA 等。郭遵强[6]提出高斯模型与 GIS 组件 MapX 集成的应急系统, 并对突发性大气污染事故进行预测, 输出结果主要用二维等值线和条形统计图的方式表现。
拉格朗日模型的基本原理是追踪并计算污染物粒子在确定性的风场及浮力和不确定性的湍流等共同作用下的轨迹常微分方程, 如 HYSPLIT[7]和 NAME等。学者们还提出一种利用拉格朗日轨迹来改进高斯模型的思路, 称为拉格朗日–高斯烟团模型, 即将污染物按一定体积分割为若干个烟团, 使用拉格朗日方法计算烟团的轨迹, 烟团内部污染物的分布则使用高斯方法计算, 最终各个烟团进行叠加得到总浓度场, 常见的烟团模型包括 CALPUFF[8]和 RIMPUFF[9]等。Sørensen 等[10]提出基于拉格朗日方法的大气应急模型 DERMA, 并与 ARGOS (accident reporting and guidance operational system)核事故决策支持系统进行集成。Sørensen 等[11]基于 RIMPUFF, 提出一种模拟病毒在大气中扩散的应急模型。Li 等[12]开发了针对高含硫气体应急响应的实时 GIS 平台,实现在复杂地形下基于 CALPUFF 模型, 使用实时数据模拟气体扩散的过程。谢颖斯[13]提出大气扩散模型 HYSPLIT 和 CALPUFF 与 GIS 组件 DotSpatial集成的应急系统。本研究选用的 CALPUFF 模型属于拉格朗日–高斯烟团模型, 往往用于较大尺度的实验区域(大于 50 km), 由于该模型综合考虑了复杂地形和气象参数等多个方面, 拥有更高的计算精度[14], 被逐渐应用于应急响应 GIS 系统中。
目前使用大气扩散模型进行案例模拟和应急预案的相关研究中, 对气体扩散的结果输出和可视化大多局限于传统的二维方法, 缺少三维方法的研究。二维表达可用于模拟地表层的污染扩散情形,却无法提供垂直方向的信息, 难以模拟真实大气环境污染的空间层级结构或展现气体在空间中扩散的过程和趋势[15]。虽然也有研究采用将气体浓度计算结果按高度分层叠加显示的方式, 在局部小场景中达到较好的三维可视化效果[16], 但将各层割裂开来, 难以直观地体现气体浓度在整个空间的连续分布情况, 并且对地表高程以及影像的分辨率要求较高, 限制了在大尺度或平原场景中的使用。
针对当前毒害性气体泄漏应急管理中存在的问题和需求, 本研究基于 CALPUFF 模型设计并实现毒害性泄漏气体的三维浓度计算方法。基于面绘制方法实现大气扩散模型的三维动态模拟。“12·23”开县特大井喷事故是我国油气开发历史上伤亡最重的安全事故之一, 造成深远的社会影响。在该案例中, 川东北地区丘陵沟壑众多, 需考虑复杂地形对泄漏气体扩散的影响; 应急救援指挥中需要快速的实时计算, 及时提供参考。本文对开县事故进行模拟和分析, 将计算结果与遥感影像、地形数据等进行集成, 实现扩散的时空过程模拟。
1 研究区概况
罗家寨、滚子坪气田位于四川省达州市宣汉县及重庆市开县境内。该地区高含硫气田的气藏储量丰富, 但地形复杂, 人口稠密, 气藏压力高, 硫化氢含量高, 气藏埋藏深, 勘探开发难度极大, 安全风险和环保风险高。
2003年12月23日, 罗家寨气田开县(现重庆市开州区)境内罗家 16H 井因操作人员违章造成井喷失控事故, 确认 243 人死亡, 4000 多人受伤, 8.3 万多人受灾。罗家 16H 为高含硫(硫化氢)气井, 属于中国石油天然气集团公司西南油气田分公司川东北气矿。该气井为中国石油天然气集团公司作为水平井开采新工艺而设计的钻井, 设计井斜深 4322 m,垂深 3410 m, 水平移位 1586.5 m, 水平段长 700 m。2003年5月23日开钻, 设计日产 100 万 m 3[17]。如图1 所示, 罗家 16H 井所在的高桥镇晓阳村位于重庆市开县西北方向, 距离县城约 80 km, 气井位于黄泥垭口附近被山丘包围着的一个平台上。该平台上还有罗家 14 号、15 号和 16 号井, 矿井周围 30 m范围内有六七户人家, 周围 1000 m 米范围内村舍林立, 约有数百户[18]。
图1 研究区位置
Fig.1 Location of studied area
2 研究方法
2.1 基于 CALPUFF 模型的技术路线
本研究提出的面向应急响应的气体扩散三维浓度场计算和可视化框架分为 3个主要模块: 底层数据库模块、三维计算模块和三维可视化模块, 总体架构如图2 所示。
图2 气体扩散三维计算及可视化流程
Fig.2 Flow chart of three dimensional calculation and visualization of gas diffusion
底层数据库以文件数据库的形式存储输入的地形高程数据(DEM)、土地利用数据(LULC)、污染源参数和气象数据(高空数据 UP.DAT 和地面数据SURF.DAT)。
CALPUFF 模型是由美国西格玛公司研发的新一代非稳态气相和空气质量建模系统, 主要由 4 部分组成: 地理和气象预处理器、气象模块 CALMET、扩散模块 CALPUFF 和后处理器 CALPOST。三维计算模块基于 CALPUFF 模型, 经过数据交换、数据预处理、气象处理、浓度计算和后处理,得到三维气体扩散浓度场。
数据交换接口实现模型与 GIS 之间的数据传输, 进行基本的投影转换、数据检查和单位转换等规范化处理。数据预处理步骤中, 特定格式的地形高程文件、土地利用文件、地表气象数据文件和高
空气象数据文件经过 TERREL 地形高程数据预处理器、CTGPROC 土地利用数据预处理器、MAKEGEO 地理数据最终预处理器、SMERGE 地面站气象数据预处理器和 READ62 探空气象数据预处理器等处理后, 可输出 CALMET 模块所需的格式化文件 GEO.DAT, SURF.DAT 和 UP.DAT 等。气象处理步骤中向 CALMET 子模块输入以上文件及包含定义模型运行信息的控制文件 CALMET.INP, 输出报告文件 CALMET.LST 及包含模型产生的小时网格气象数据的无格式的磁盘文件 CALMET.DAT 或PACOUT.DAT, 建立模拟区域的三维格点风温场。在浓度计算步骤中, 由 CALPUFF 子模块读入 CALMET 生成的气象场及用户输入的污染源参数, 循环修改控制文件 CALPUFF.INP 中的受体(需要计算浓度值的点位)层高度, 计算该层受体点的浓度及干湿沉降通量, 输出 CONC.DAT 文件, 循环结束后得到网格化的三维浓度场。后处理步骤中, 由 CALPOST 子模块对这些输出文件进行处理, 整理为用户指定的格式并进行能见度模拟。本研究将上一步生成的 CONC.DAT 转化为 n个格式化的 grd 文件, n 为系统模拟的扩散时间(分钟级别), 即第 i个 grd 文件描述开始扩散后第 i 分钟的污染气体浓度分布。
三维可视化模块将受体点数据读入, 利用 MarchingCubes 面绘制算法追踪, 生成指定层数和浓度值的等值面时序三维模型。将模型读入三维软件中, 同时叠加高程和地表影像, 实现扩散场景的三维动态可视化。
2.2 改进的气体扩散模型多层计算
CALPUFF 目前只能一次性地输出二维浓度场,因此许多研究只关注近地表的气体浓度分布。实际上, CALPUFF 模型中受体点的高度信息在 CALPUFF 模块控制文件 CALPUFF.INP 中的“非网格(离散)受体信息”组中指定, 通过修改这一参数, 就可以实现对不同高度层的浓度场计算及多层结果的输出。
本研究完善了气体扩散模型的计算流程, 将二维计算变为三维计算。改进后的流程可以根据用户指定的层高及层数, 在一系列空中层设置受体点,并结合气体类型、泄漏速率和泄漏点距地面高度等参数, 在这些层上执行浓度计算, 得到多层浓度值,最终整合生成三维浓度场, 为后续可视化提供基础数据。进行多层计算改进后 CALPUFF 模型的计算流程如图3 所示。
图3 基于CALPUFF的三维浓度场计算模型
Fig.3 Three dimensional concentration field calculation model based on CALPUFF
在 CALPUFF 模块程序中, 通过循环修改 CALPUFF.INP 文件的“非网格(离散)受体信息”组, 指定相应的高度参数, 开展不同层计算, 得到逐层扩散浓度结果, 该结果不断写入网格化的浓度场文件CONC.DAT。在 CALPOST 后处理模块, 将 CONC.DAT 整理成格式化的输出文件, 在输出路径中, 创建文件夹“Layer0”至“Layerx”, 分别保存各层的气体浓度时序结果。每层的文件夹内, 则从气体泄漏开始时刻起, 以 5 分钟为间隔, 按时刻依次创建文本文件, 达到时空数据结构化存储的目的。算法流程如下。
算法1 基于 CALPUFF 的三维浓度场计算算法
输入 网格化逐时风场文件 CALMET.DAT; CALPUFF 控制文件 CALPUFF.INP; CALPOST 控制文件CALPOST.INP; 受体层数N; 受体层高 deltah
输出 三维气体浓度场文件夹Layer0至Layer(N–1)
1.for each n∈[0, N) do
2.计算此受体层高度 h = n * deltah
3.修改 CALPUFF.INP 文件, 更新各受体高度为h
4.向 CALPUFF 模块输入 CALPUFF.INP 和 CALMET.DAT, 执行计算
5.修改 conc.dat 文件
6.end for
7.向 CALPOST 模块输入 CALPOST.INP 和 conc.dat,执行后处理
8.输出格式化后的文件列表Layer0至Layer (N–1)
2.3 气体扩散时序数据三维重建方法
三维可视化的方法分为两大类: 面绘制和体绘制。面绘制是对三维空间目标的表面进行绘制显示, 其数据处理基础是序列二维图像数据, 再结合其空间结构关系, 还原三维空间结构。面绘制主要分为体素级重建方法和切片级重建方法, 其中, 体素级重建方法的本质是在一个三维的数据场中生成等值面的过程[19], 其实现多基于经典的 Marching Cubes (移动立方体)算法。本研究选用 Marching Cubes 算法对气体扩散过程进行三维可视化。
Marching Cubes 算法由 Lorenson [20]等提出, 以体素为单位, 制造出想要重构的部分和复杂背景的分界[21], 然后从体素中抽取三角面片来拟合此边界, 大量的三角面片构成一个等值面。等值面的构造经过两个主要计算过程: 1) 体素中三角面片逼近等值面的计算; 2) 三角面片各顶点法向量的计算[22]。Marching Cubes 算法的基本假设是沿着体素的棱边,数据场呈连续线性变化。构建某个气体浓度值的扩散等值面时, 当搜索到相邻两受体点处的浓度值分别大于和小于等值面值时, 在该条边上必有且仅有一点是这条边与等值面的交点。根据这一原理, 就可以判断所求等值面与哪些体素相交。每个体素包含 8个顶点, 根据顶点处数据值大于或小于等值面的值, 将顶点标记为 1 或 0 两种状态, 故每个体素有 28= 256种组合状态, 根据互补对称性和旋转对称性简化为 15种。得到面片后, Marching Cubes 方法采用中心差分计算出体素各顶点处的梯度, 再次通过线性插值求出三角面片各顶点的梯度, 即各顶点处的法向量, 然后选择适当的光照模型进行计算,生成具有真实感的图形[20]。
但是, Marching Cubes 算法存在连接方式的二义性问题, 相邻体素连接不当时, 可能出现拓扑不一致, 从而使等值面出现孔隙。解决二义性问题的方法主要有二义面平均值判定法、子构型查找表法、梯度一致性准则法和渐近线判别法。本研究使用子构型查找表法, 即对于 15种基本构型中具有二义性的构型, 各自建立一个子查找表, 每个子查找表包含两种三角剖分方式。此外, 还需存储一个表来记录这些三角剖分方式的相容性[23]。
本文模型首先读取 CALPOST 输出的逐时和逐层的气体浓度文本文件, 以点数据结构记录受体网格的三维坐标, 以浮点数组记录对应位置的浓度值,每个受体点即为立方体体素的一个顶点, 进而从三维浓度场追踪提取等值面。算法 2 为使用 Marching Cubes 算法追踪等值面的流程。
算法 2 基于 Marching Cubes 的三维浓度场等值面追踪算法
输入 受体点位置坐标 points; 浓度场数组 conc;受体网格维度 dims[3]; 等值面数 contNum; 等值面值数组 contValue [contNum]; 体素中面片连接方式表triangleCases[256]
输出 等值面集合 contours
1.for each k∈[0, dims[2]) do
2.for each j∈[0, dims[1]) do
3.for each i∈[0, dims[0]) do
4.由(i, j, k), (i+1, j, k), (i+1, j+1, k), (i, j+1,k), (i, j, k+1), (i+1, j, k+1), (i+1, j+1, k+1), (i, j+1, k+1)组成当前受体点体素的 8个顶点, 从 conc 中取出顶点处浓度值, 记入数组s[8]
5.if 该受体点处梯度未计算 then
6.调用中心差分法函数计算 8 顶点梯度, 记入数组gradient[8]
7.for each n∈[0, contNum) do
8.当前追踪的等值面值 value =contValues[n]
9.根据 8个顶点值 s[8]与 value的关系确定当前体素的索引下标值index
10.得到该体素中面片连接方式triangleCases[index]
11.构造该体素中的三角面片 triCase
12.for each edge∈triCase -> edges do
13.由中心差分法计算三角面片各顶点的坐标、浓度值、梯度值
14.将当前三角面片插入等值面集合contours中
15.end for
16.end for
17.end for
18.end for
19.end for
图4(a)为某个时刻原始浓度场的一部分, 为 20行×14 列×11 层的矩阵。在该模拟设置中, 受体网格的水平间隔为 500 m, 垂直间隔为 20 m, 矩阵中心红色点位处表示气体浓度较高, 即污染源位置,蓝色点位处则表示气体浓度较低。图4(b)为由该浓度场重建的等值面模型, 可以看出气体由浓度最高的污染源处向东、向外层扩散, 浓度逐渐下降, 与原始浓度场表现的变化趋势一致。另外, 追踪气体等值面时, 为保证三维表面的封闭性, Marching Cubes 算法会将原始数据范围向外扩展, 最终气体模型顶部达到地表以上约 600 m, 超出受体网格的垂直范围, 因此该扩展部分的数值主要由地表起伏趋势和临近受体点的数值决定。
图4 由原始浓度场追踪得到等值面
Fig.4 Isosurfaces from the original concentration field
为评估等值面追踪的时间效率, 本研究针对同一气体扩散案例, 改变受体网格分辨率, 对同一区域进行扩散过程模拟。选取一个特定浓度值的等值面进行构建, 对相关效率指标进行统计分析。表1选取扩散初期至后期的 9个时刻, 列出当时的非零值受体点数、构成等值面的顶点数、三角面片数以及等值面追踪所需的时长。其中, 网格分辨率为 200, 500 和 1000 m 时的数据体维度分别为 100×100×11、40×40×11 和 20×20×11。程序运行的电脑配置为 Windows10 操作系统, 16 G 内存, CPU采用 Intel(R) Core(TM) i7-1065G7 1.30 GHz, 基于Pycharm Community Edition 2020 平台, 使用 Python语言编写测试程序。
由表1 可知, 对于某一固定分辨率, 每个时刻的等值面构建时长主要受非零值受体点数影响, 随着扩散时间增加, 模拟区域内浓度值超过零值的受体点数增加, 所需的构建时间也增加。当分辨率为200 m 时, 在 20~180 分钟, 非零值受体点数由 110增至3355, 所需时长由 152.36 ms 递增至998.68 ms。三角面片数与等值面具体形态相关, 对构建时长的影响不显著。
表1 不同分辨率下等值面构建效率指标对比
Table 1 Comparison of efficiency indexes of isosurface construction under different resolutions
受体网格效率评估指标模拟时刻非零值受三角面等值面构分辨率/m/min体点数 片数 建时长/ms 20 110 24560 152.36 40 242 18200 177.88 60 484 19432 307.13 80 814 19820 468.80 200100 1309 18788 613.50 120 1837 18596 691.76 140 2453 18424 747.72 160 3256 17424 1019.16 180 3355 17908 998.68 20 33 6752 39.95 40 66 15280 53.18 60 132 16944 55.16 80 176 17176 68.31 500100 220 18792 118.81 120 264 19204 148.48 140 308 19584 184.02 160 418 18868 137.03 180 418 19804 152.49 20 11 25404 22.06 40 11 28652 29.00 60 33 5804 40.43 80 44 10976 26.76 1000100 66 13044 56.31 120 88 14032 80.86 140 99 17056 73.82 160 132 17628 78.02 180 132 17848 94.79
当可用的非零网格点过少时, 得到等值面的形态可能失真, 故越高的网格分辨率越利于三维场景的重建。综合考虑时间效率, 使用 200 m 分辨率时,每个时刻的构建时长基本上控制在 1 秒以内, 可以满足应急场景的需求。故本研究使用 200 m 浓度场网格。
3 气体泄漏事故模拟案例研究
3.1 实验参数设置
事故发生于 2003年12月23日晚 21 时 55 分, 24日 15 时发生井喷的主管道成功封堵, 16 时放喷管线实施点火, 共持续 18 小时。当日该地气温为 4.6~8.0℃, 相对湿度为 94%~99%, 平均风速为 0.13 m/s,最大风速为 0.7 m/s, 风向西北偏西。根据现场调研,罗家 16H 气井地层压力高达 40 MPa 以上, 井口喷出的气体富含 H2S, 预计无阻流量为 400×104~1000×104 m3/d, 天然气中 H2S 含量约为 151 mg/m3[24]。事故造成井场周围居民和井队职工共 243 人死亡, 上万人因硫化氢中毒接受治疗, 九万余人疏散转移, 经济损失上亿元[25]。
本研究关注毒害性气体泄漏后, 对污染源周围人居环境和自然环境的影响。研究区地形起伏多变, 以污染源为中心的 50 km×50 km 范围内海拔为50~1738 m, 平均 996 m, 因而气体扩散情况复杂,垂直方向的变化明显。我们在 z 轴方向设置从地表向上, 以 20 m 为间距的 11 层受体点。另外, 在突发气体泄漏事故的应急场景中, 由于初始喷发速度高、喷发气体量大, 有毒气体浓度在时间和空间上的变化迅速, 常规大气监测中以小时为单位的模拟时间步长无法满足应急场景的需求。一项关于 CALMET 时空分辨率对 CALPUFF 模拟浓度场影响的研究结果[26]显示, 10 min 的时间步长即可达到较好的模拟结果, 能够满足实际应用的需求, 本研究选择更精细的 5 min 时间步长。在 CALPUFF 中输入的污染源、气象条件及模拟设置参数如表2 所示。
根据 H2S 气体的性质, 本研究将追踪浓度分为4 级, 由扩散外层向内依次为 0.01, 6.9, 15 和 100 mg/m3。表3 为气体浓度超标阈限和对人体产生不同程度伤害的阈值。
表3 硫化氢对人体的影响及阈限
Table 3 Effects and threshold of hydrogen sulfide on human body
浓度/(mg·m–3) 说明0.01 嗅觉阈6.9 有明显且浓厚的难闻气味15 有令人讨厌的气味, 眼睛可能受刺激100 出现眼及呼吸道刺激症状
3.2 模拟结果分析
3.2.1 三维时序扩散过程
将表2 中参数输入计算系统, 经过模型计算和三维重建, 得到气体泄漏后快速扩散阶段的动态模拟过程(图5), 分别表示事故发生后 15, 30, 45, 60,75, 90, 120, 150 和 180 分钟的情形。
表2 模拟参数
Table 2 Simulation parameters
参数类型污染源属性 参数名称 取值UTM坐标(km) 237.321, 3474.793烟囱高度(m) 10出口初速度(m/s) 40恒定排放速率(g/s) 6076.39气象参数风速(m/s) 0.3风向(°) 270温度(℃) 7时间步长(min) 5模拟时间设置 总时长(h) 10模拟范围设置X方向长度(km) 20 Y方向长度(km) 20风场网格尺寸(m) 200浓度场网格尺寸(m) 200
在本次事故中, 气井所在的开县高桥镇及相邻的麻柳乡、正坝镇、天和乡 4个乡镇 9.3 万余人受灾。其中, 受灾最重的是高桥镇的晓阳和高旺两个村, 共 2419 人受到影响。90%的受影响群众居住在距离气井 2 km 外, 但由于暴露在高硫化氢环境中,90%的遇难者均位于气井周围 1.5 km 以内。由图5 的模拟结果可以观察到, 由于井喷位置处于低洼地带, 且当天气象条件较稳定, 风速小且空气湿度大, 扩散范围总体上呈现偏心半球形, 有毒气体从污染源水平向外和垂直向上扩散, 气体半球的中心在偏西风作用下位于污染源下风的东南方向。从整体上看, 高浓度硫化氢聚集在污染源附近一定范围内, 这些特点与实际受灾情况基本上相符。
图5 气体扩散过程三维动态可视化
Fig.5 Three dimensional dynamic visualization of gas diffusion process
3.2.2 典型时刻二维和三维对比
将气体等值面三维模型导出, 叠加地形和影像,并与地表二维扩散结果进行逐时对比。
在事故发生 15 分钟时, 污染源东南方向 500 m范围内, 近地表的硫化氢气体浓度已超过 15 mg/m3(图6(a)中橙色区域)。晓阳村位于污染源西南方向约 500 m 处, 此时硫化氢浓度约 3 mg/m3; 高旺村位于其东偏南方向约 800 m 处, 浓度达到约 6 mg/m3。据调查, 离井口最近的晓阳村村民听到井喷刚发生时的轰然巨响后, 不到 10 分钟, 便闻到一股臭鸡蛋味, 即已明显超过嗅觉阈值。在井喷开始后约10 分钟, 晓阳村和高旺村的硫化氢浓度均迅速从 0升至40 mg/m3 [2-3]。考虑到 CALPUFF 模型多用于长时间尺度的环境评价领域, 常规应用中一般采用模式默认的模拟方案和参数设置, 如逐时风场、默认的湍流参数化方案和地表参数等[27], 在应急情景中的模拟结果容易偏小, 但从数值量级的角度看,本方法的短时模拟结果仍然具有一定的参考意义。
事故发生 30 分钟时, 污染源周围 1000 m 以内的大部分区域硫化氢浓度已达到 15 mg/m3 (图6(c)中橙色区域), 该浓度已经能够对人眼产生刺激性,并伴随明显的异味。位于泄漏源东南方向约 500 米的高桥镇已完全处于蓝色区域, 根据基于 Fluent 的数值模拟结果[24], 井喷后半小时的硫化氢浓度已超过人的嗅觉阈值, 与本研究结果一致。
事故 1 小时后, 两座村庄已经完全被浓度约100 mg/m3 的硫化氢气体覆盖。该浓度的硫化氢气体会刺激人的眼睛及呼吸道, 甚至危害生命安全。从二维结果中可以观察到, 近地表气体被划分为以污染源为圆心的两个扇形区域。对应图6(f)的三维地形可知, 在污染源北、西南、东南方向均为狭长山谷, 气体从靠近西侧山丘的山脚处喷发, 同时在偏西风作用下, 沿山谷东侧蔓延。由于周围山体的阻挡, 高浓度硫化氢气体在污染源附近的洼地中不断聚集。另外, 从三维地形观察到, 高旺村位于一座低山的西坡, 相对于污染源的海拔高度更高, 但此处整个山坡都被高浓度气体覆盖, 可见在有风、离污染源较近的情况下, 地势较高处也非安全地带。辛保泉等[28]通过风洞实验也发现, 泄漏源下游为山峰条件时, 气体可能与山体碰撞或在山体背后的山坡上形成高浓度的累积。
图6 事故发生15, 30, 60, 90分钟的气体浓度分布
Fig.6 Gas concentration distribution at 15, 30, 60 and 90 min after accident
事故后 1.5 小时内, 污染源周围 5 km 均已受到浓度超标的硫化氢气体影响(图6(g)中蓝色区域)。该范围也是本次事故中报道受灾最重的区域, 居民均需要被疏散转移。
3.2.3 扩散模型剖面分析
由图6 可以看出, 三维空间模型可以完整、直观地反映事故发展过程和危害影响范围。为了呈现随空间位置变化的更细致的气体扩散特征, 可以配合垂直方向二维剖面图辅助来综合分析。图7 为事故发生 60 分钟时的扩散二维图像, 沿线 AA′, BB′,CC′和 DD′得到下风向 0, 500, 1000 和 1500 m 处横向扩散的气体模型剖面图(图8)。
图7 事故发生60分钟的二维浓度分布图及剖面线位置
Fig.7 Two-dimensional concentration distribution and locations of profile lines at 60 min after accident
可以发现, 在该时刻, 在下风向 1500 m 以内不同距离处, 气体横向扩散的范围大致相同, 但不同浓度等级的圈层范围逐渐变化。0 m 处, 即经过泄漏源点的剖面中, 浓度 100mg/m3 以上的占比很大(图8(a)中红色范围), 表明此处灾害严重, 而沿下风向则逐渐减轻, 橙色范围占比增加。不同位置处剖面的轮廓形态也差异很大, 可以观察到它们受地形起伏影响显著, 高浓度气体易在低洼处聚集。对比图8(a)~(d)的地表剖面线, 每隔 500 m 的地形趋势几乎截然不同, 可见事故发生地地形复杂, 会给灾害预测和应急救援带来很大的不确定性和难度。
图8 事故发生60分钟的下风向不同位置处气体剖面图
Fig.8 Gas profiles at different locations downwind at 60 min after accident
3.2.4 高桥案例总结分析
通过对开县案例扩散过程的模拟分析, 可以得到如下结论。
1) 该案例中, 气体的扩散主要受地形和气象条件影响, 危害区域分布在顺风的山坡前部或山谷洼地, 该区域不适宜居民居住, 一旦发生突发安全事故, 后果不堪设想, 应当在气井选址和民房选址时予以考虑。高排放率的有毒气体在小风速和低洼地形作用下迅速覆盖污染源周边, 并以偏心半球形态不断扩大蔓延, 在高压、高含硫且有风的情况下,污染源附近会受到高浓度毒害气体影响, 向附近高地疏散并非安全的应对措施。
2) 污染源附近是发生事故后气体浓度迅速上升且最大浓度持续出现的区域, 在川东北地区, 常年平均风速较小, 且丘陵地貌不利于气体的水平输送和扩散, 人员撤离半径可以酌情划定为距离污染源 5 km。
4 结论与展望
本研究针对传统的气体扩散计算和可视化问题, 设计并实现基于 CALPUFF 模型的气体浓度计算和三维可视化框架, 将模拟维度从近地表的二维扩展为动态的考虑垂直分布的空间三维加时间维度, 能够更直观、方便地进行应急模拟和事故发生机理的分析。将该方法应用到开县井喷事故案例中, 验证了其可行性和有效性, 并得到如下主要结论。
气体泄漏后在不同高度层的扩散分布不同, 通过三维可视化的方式更易于从多个维度确定气体扩散的范围, 能够对气体的空间分布有连续的和直观的认识; 叠加三维地形场景后, 便于分析各地的受灾程度以及扩散形态的形成原因; 加入时间维度后,从动态扩散过程更能反映灾情的发展变化趋势, 对于应急事件的预演以及宏观把握具有重要的指导意义。
本研究还存在一些不足之处, 未来应进行如下方面的研究。
1) 在模型方面, CALPUFF 模型是为长时间、大范围的大气质量模拟评估所设计, 应用于应急场景时存在一些难以处理的问题: 应急场景具有极强的突发性和危害性, 要求模型必须同时具有高计算速度和高计算精度, 而现有的 CALPUFF 等大气扩散模型空间分辨率为十米级, 气象输入参数为秒级,时空分辨率不能自洽; 进行模拟参数设定时, 气体排放率、温度和风速等输入参数为恒定值, 但在实际场景中, 这些条件往往在不断的变化, 导致模拟结果与现实的差异。
2) 在三维计算方面, 多层受体、高网格分辨率的设定导致计算量增大, 加上 CALPUFF 模型考虑因素的复杂性, 计算时长通常在分钟级。若要提高计算效率, 可以考虑应用 CPU 多核并行计算以及GPU 加速等技术。
3) 在体视化效果方面, Marching Cubes 算法存在不完善之处, 如数据量较小时等值面形态失真、追踪范围不够精确等, 可以通过算法的改进得到修正。根据模拟目的的不同, 如关注事故对近地表居民区的影响, 关注毒害物质对山间动植物的危害、关注气体扩散对高层大气的污染, 等等, 可以选择性地设定受体层的数量和层高, 聚焦于特定的三维空间范围, 以期提高计算模拟效率。
[1]Feng J R, Gai W M, Yan Y B.Emergency evacuation risk assessment and mitigation strategy for a toxic gas leak in an underground space: the case of a subway station in Guangzhou, China.Safety Science, 2021,134: 105039
[2]陈建国, 潘思铭, 刘奕, 等.复杂地形下有害气体泄漏的模拟预测与输运规律.清华大学学报(自然科学版), 2007, 47(3): 381–384
[3]练章华, 周兆明, 王辉, 等.特大井喷 H2S 扩散的数值模拟分析.天然气工业, 2009(11): 112–115
[4]Carruthers D J, Holroy D R J, Hunt J C R, et al.UKADMS: a new approach to modelling dispersion in the earth’s atmospheric boundary Layer.Journal of Wind Engineering and Industrial Aerodynamics, 1994,5(52): 139–153
[5]U.S.Environmental Protection Agency.User’s guide for the AMS/EPA regulatory model (AERMOD)[DB/OL].(2019–08)[2022–01–01].https://www3.epa.gov/ttn/scram/models/aermod/aermod_userguide.pdf
[6]郭遵强.大气污染扩散模型的研究及在环境应急系统中的实现[D].青岛: 中国海洋大学, 2008
[7]Stein A F, Draxler R R, Rolph G D, et al.NOAA’s HYSPLIT atmospheric transport and dispersion modeling system.Bulletin of the American Meteorological Society, 2015, 96(12): 2059–2077
[8]Scire J S, Strimaitis D G, Yamartino R J.A user’s guide for the CALPUFF dispersion model [DB/OL].(2000–01)[2022–01–01].http://www.src.com/calpuff/download/CALPUFF_UsersGuide.pdf
[9]Thykier-Nielsen S, Deme S, Mikkelsen T, et al.Description of the atmospheric dispersion module RIMPUFF [DB/OL][2022–01–01].(1999–04)[2022–01–01].https://resy5.iket.kit.edu/RODOS/Documents/Public/Handbook/Volume3/4_2_6_RIMPUFF.pdf
[10]Sørensen J H, Baklanov A, Hoe S.The danish emergency response model of the atmosphere (DERMA).Journal of Environmental Radioactivity, 2007, 96:122–129
[11]Sørensen J H, Jensen C Ø, Mikkelsen T, et al.Modelling the atmospheric dispersion of foot-and-mouth disease virus for emergency preparedness.Physics and Chemistry of the Earth, Part B: Hydrology,Oceans and Atmosphere, 2001, 26(2): 93–97
[12]Li M, Liu H, Yang C.A real-time gis platform for high sour gas leakage simulation, evaluation and visualization // ISPRS Annals of Photogrammetry,Remote Sensing and Spatial Information Sciences.Fairfax, VA, 2015: 225–231
[13]谢颖斯.基于 GIS 的突发事故核污染物扩散模拟系统研发与应用[D].广州: 华南理工大学, 2014
[14]李梅, 杨冬偶, 何望君.大气扩散模型 AERMOD 与CALPUFF 对比研究及展望.武汉大学学报(信息科学版), 2020, 45(8): 1245–1254
[15]Ridzuan N, Ujang U, Azri S, et al.Visualising urban air quality using AERMOD, CALPUFF and CFD models: a critical review // The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences.Turkey, 2020: 355–363
[16]邬群勇, 黄君毅, 盛玲, 等.CALPUFF 模型模拟结果的时空多维可视化表达.地球信息科学学报,2015, 17(2): 206–214
[17]陈新, 张华东, 陈荣光, 等.开县特大井喷事件及严重后果成因与反思.现代预防医学, 2007(12):2229–2231
[18]阳作明, 白青松, 伍亚光.开县“12.23”天然气井喷硫化氢扩散气象条件分析 // 中国气象学会 2006年年会“大气成分与气候, 环境变化”分会场论文集.成都, 2006: 443–447
[19]Pan Z, Song T, Guo M, et al.Comparison of medical image 3D reconstruction rendering methods for robotassisted surgery // 2nd International Conference on Advanced Robotics and Mechatronics (ICARM).Hefei, 2017: 94–99
[20]Lorensen W E, Cline H E.Marching cubes: a high resolution 3D surface construction algorithm.ACM SIGGRAPH Computer Graphics, 1987, 21(4): 163–169
[21]Farías D I H, Cabrera R G, Fraga T C, et al.Modification of the marching cubes algorithm to obtain a 3D representation of a planar image.Programming and Computer Software, 2021, 47(3): 215–223
[22]Newman T S, Yi H.A survey of the marching cubes algorithm.Computers & Graphics, 2006, 30(5): 854–879
[23]秦绪佳.医学图象三维重建及可视化技术研究[D].大连: 大连理工大学, 2001
[24]李剑峰, 姚晓晖, 刘晓琴.基于 Fluent 的开县井喷事故后果模拟与分析.环境科学研究, 2009, 22(5):559–566
[25]Li J, Zhang B, Yang W, et al.The unfolding of“12.23” Kaixian blowout accident in China.Safety Science, 2009, 47(8): 1107–1117
[26]康凌, 朱好, 黄倩倩, 等.CALMET 时空分辨率对CALPUFF 模拟浓度场的影响.北京大学学报(自然科学版), 2021, 57(6): 1006–1018
[27]朱好, 张宏升, 蔡旭晖, 等.CALPUFF 在复杂地形条件下的近场大气扩散模拟研究.北京大学学报(自然科学版), 2013, 49(3): 452–462
[28]辛保泉, 喻健良, 党文义, 等.复杂地形高含硫天然气风洞扩散实验及安全防护距离.天然气工业,2020, 40(11): 149–158