高斯烟羽与烟团模型：从理论假设到GIS空间可视化实战-平芜编程栈

1. 高斯模型：从烟雾到数学的奇妙转化

第一次接触高斯烟羽模型时，我正参与一个化工厂周边空气质量评估项目。站在厂区外看着烟囱冒出的白烟，突然意识到那些看似随意的飘散轨迹，竟然可以用数学公式精确描述。这就像用天气预报来预测烟雾的"舞蹈"，而背后的核心就是高斯分布这个统计学中的经典概念。

高斯模型本质上是用正态分布来描述污染物在空气中的扩散规律。想象把一瓶香水打翻在房间中央，香味会逐渐向四周扩散 - 离得近的地方浓度高，远的地方浓度低，形成一个浓度梯度。高斯模型就是用数学语言来描述这种扩散过程，只不过把房间换成了大气环境，香水换成了工业排放物。

在实际应用中，我们主要使用两种变体：烟羽模型和烟团模型。烟羽模型适合描述持续排放的场景，比如火力发电厂常年运行的烟囱；而烟团模型更适合突发事故，比如化工厂的泄漏事件。这两种模型就像天气预报中的"持续降雨"和"短时雷阵雨"的区别，用的都是降水原理，但具体表现方式不同。

2. 深入理解高斯烟羽模型

2.1 模型的核心假设

让我用一个厨房的例子来解释这些抽象假设：假设你在厨房煮一锅汤，蒸汽持续从锅面升起。如果厨房门窗紧闭（稳定风场），没有抽油烟机干扰（无下垫面影响），蒸汽会均匀地向四周扩散（正态分布）。这就是高斯烟羽模型的理想场景。

具体来说，模型建立在四个关键假设上：

风场稳定不变，就像开着恒定档位的电风扇
污染物在水平和垂直方向的分布都遵循钟形曲线（正态分布）
污染物总量守恒，不会凭空消失或增加
污染源持续稳定排放，不是一阵一阵的

这些假设虽然理想化，但为建模提供了可行的简化方案。就像物理课上假设"无摩擦力"一样，我们先建立理想模型，再逐步加入现实世界的复杂因素。

2.2 模型方程解析

高斯烟羽模型的数学表达式看起来复杂，但拆解后很容易理解：

def gaussian_plume(Q, u, σy, σz, y, z, H): """ Q: 源强 (mg/s) u: 风速 (m/s) σy,σz: 水平和垂直扩散参数 y: 横向距离 (m) z: 高度 (m) H: 有效源高 (m) """ term1 = Q / (2 * np.pi * u * σy * σz) term2 = np.exp(-(y**2)/(2 * σy**2)) term3 = np.exp(-((z-H)**2)/(2 * σz**2)) + np.exp(-((z+H)**2)/(2 * σz**2)) return term1 * term2 * term3

这个公式计算的是下风向某点(x,y,z)处的污染物浓度。其中最后一项的相加考虑了地面反射效应 - 就像声音碰到墙壁会产生回声一样，污染物碰到地面也会"反弹"。

3. 高斯烟团模型：应对突发泄漏

3.1 与烟羽模型的区别

去年参与一次化工厂应急演练时，我深刻体会到烟团模型的价值。当阀门突然爆裂，大量氯气瞬间释放，这时烟羽模型就不适用了 - 因为泄漏不是持续的，而是一次性的"烟团"。

烟团模型把污染物看作一个不断膨胀的气球，计算它在移动过程中浓度的时空变化。想象把一滴墨水瞬间滴入流动的溪水，观察墨水团随水流扩散的过程，这就是烟团模型的直观表现。

3.2 烟团模型方程

烟团模型的数学表达与烟羽模型类似，但增加了时间维度：

def gaussian_puff(Q, t, σx, σy, σz, x, y, z, H): """ t: 时间 (s) σx: 沿风向扩散参数 其他参数同烟羽模型 """ term1 = Q / ((2*np.pi)**1.5 * σx * σy * σz) term2 = np.exp(-0.5 * ((x-u*t)**2/σx**2 + y**2/σy**2)) term3 = np.exp(-0.5 * (z-H)**2/σz**2) + np.exp(-0.5 * (z+H)**2/σz**2) return term1 * term2 * term3

这个公式多了时间t和沿风向扩散参数σx，因为烟团在移动过程中不仅横向和垂直扩散，还会沿风向拉伸变形。

4. 关键参数获取与计算

4.1 扩散参数的确定

扩散参数σy和σz是模型中最棘手的部分，它们取决于大气稳定度和下风向距离。我国《环境影响评价技术导则》推荐使用Pasquill-Gifford曲线来估算这些参数。

实际操作中，我通常使用查表法。比如对于城市地区的B类稳定度，σy和σz可以表示为：

def get_sigma_yz(stability_class, x): # x为下风向距离(km) if stability_class == 'B': a_y, b_y = 0.16, 0.82 a_z, b_z = 0.11, 0.86 σy = a_y * x ** b_y # 水平扩散参数(m) σz = a_z * x ** b_z # 垂直扩散参数(m) return σy, σz

4.2 烟气抬升高度的计算

烟气从烟囱排出后，由于热力和动力作用会继续上升，这个抬升高度△H对模型结果影响很大。我常用Briggs公式计算：

def briggs_plume_rise(Qh, u, Fb): """ Qh: 热排放率(kJ/s) u: 风速(m/s) Fb: 浮力通量参数 """ if Fb < 55: Δh = 21.425 * Qh**0.75 / u else: Δh = 38.71 * Qh**0.6 / u return Δh

记得有次评估电厂项目时，忽略了这个抬升高度，结果预测的地面浓度比实测高出一个数量级，闹了个大笑话。

5. GIS空间可视化实战

5.1 数据准备与坐标转换

将高斯模型与GIS结合时，首先要解决坐标统一问题。我常用PyProj库进行坐标转换：

from pyproj import Transformer def lonlat_to_utm(lon, lat): transformer = Transformer.from_crs("EPSG:4326", "EPSG:32650") # WGS84转UTM50N x, y = transformer.transform(lat, lon) return x, y

转换后，我们可以创建规则网格来计算每个格点的浓度值。比如生成1km×1km的网格：

import numpy as np def create_grid(x0, y0, size=1000, resolution=100): x = np.arange(x0, x0+size, resolution) y = np.arange(y0, y0+size, resolution) return np.meshgrid(x, y)

5.2 浓度场计算与可视化

有了网格后，就可以计算每个点的浓度值。以烟羽模型为例：

def calculate_concentration_grid(Q, u, H, stability, grid_x, grid_y): concentrations = np.zeros_like(grid_x) for i in range(grid_x.shape[0]): for j in range(grid_x.shape[1]): x = grid_x[i,j] - x0 y = grid_y[i,j] - y0 σy, σz = get_sigma_yz(stability, x/1000) # km转换 concentrations[i,j] = gaussian_plume(Q, u, σy, σz, y, 0, H) # 地面浓度z=0 return concentrations

使用Matplotlib或GIS软件可以将结果可视化。我更喜欢用QGIS的栅格渲染功能，它能创建更专业的专题图：

import matplotlib.pyplot as plt def plot_concentration(concentration, grid_x, grid_y): plt.figure(figsize=(10,8)) plt.contourf(grid_x, grid_y, concentration, levels=20, cmap='jet') plt.colorbar(label='浓度 (mg/m³)') plt.xlabel('东向坐标 (m)') plt.ylabel('北向坐标 (m)') plt.title('污染物浓度分布') plt.show()

6. 实际应用中的注意事项

6.1 模型局限性

高斯模型虽然简单实用，但也有明显局限。记得有次在山谷地区做项目，模型预测结果与实测相差甚远 - 因为地形对气流的影响太大。这种情况下需要更复杂的CFD模型。

主要局限包括：

不适合复杂地形（山地、城市峡谷）
假设条件在长时间（>1小时）模拟中可能不成立
对化学反应和沉降过程考虑不足

6.2 参数敏感性分析

模型结果对某些参数非常敏感。建议进行敏感性分析，比如风速变化±20%对结果的影响：

def sensitivity_analysis(Q, H, stability, grid_x, grid_y): base_u = 3.0 # 基准风速3m/s concentrations = {} for u in [base_u*0.8, base_u, base_u*1.2]: concentrations[u] = calculate_concentration_grid(Q, u, H, stability, grid_x, grid_y) return concentrations

通过这种分析，可以识别关键参数，在数据收集时给予更多关注。

7. 完整案例：化工厂泄漏模拟

去年参与的某化工厂氯气泄漏应急演练，完整展示了从模型计算到GIS可视化的全流程：

获取泄漏参数：
- 源强Q：通过泄漏孔径和压力计算，约5000mg/s
- 泄漏高度H：15m（设备离地高度）+3m（抬升高度）=18m
- 气象条件：风速2.5m/s，稳定度D类
建立模拟区域：
- 以泄漏点为中心，2km×2km范围
- 网格分辨率50m

计算浓度场并可视化：

Q = 5000 # mg/s u = 2.5 # m/s H = 18 # m stability = 'D' grid_x, grid_y = create_grid(x0, y0, size=2000, resolution=50) conc = calculate_concentration_grid(Q, u, H, stability, grid_x, grid_y) plot_concentration(conc, grid_x, grid_y)