做气体泄漏扩散模拟和大气环境影响评价的同行应该都绕过不了一个东西高斯扩散模型。前段时间我整理项目代码发现每次都要重复写一段“给扩散时间算三个方向扩散参数”的逻辑干脆收敛成一个函数输入t扩散时间秒返回σx、σy、σz。代码没几行但外人看着简单真正落地的时候坑全藏在参数表、稳定度分级和适用范围里。如果你也在做烟团模型、毒气泄漏后果分析、或者污染物浓度分布计算这篇内容应该能帮你少踩几个雷。我不只讲函数怎么实现还会把背后的模型原理、参数怎么选、边界怎么处理、验证结果多大量级算合理一次性讲清楚。1. 扩散参数到底在算什么高斯模型的三维离散逻辑1.1 σx、σy、σz的物理意义先别急着看代码。你要彻底理解这个函数在算什么才能知道什么时候能用、什么时候要换模型。高斯扩散模型的核心假设是污染物在三维空间中的浓度分布近似服从高斯分布。所谓σx、σy、σz就是三个方向上的标准差单位是米。物理上它们刻画的是“烟团”或者“烟羽”在x方向顺风向、y方向水平侧向、z方向垂直方向上扩散开来的尺度。你可以把它想象成一滴墨水滴进静止的水里刚开始墨迹是一个清晰的小圆点过一会儿变成一团云边缘模糊范围越来越大。这个“范围”没法用一个半径描述因为浓度从中心向外是连续递减的所以就用标准差来衡量。σ越大污染物分布越分散中心浓度越低σ越小污染物集中在轴线上峰值浓度就很高。对瞬时释放更准确的叫法是“烟团”比如储罐瞬间破裂后物料直接抛洒到空气中。这种场景下空间某一点的浓度可以写成C(x, y, z, t) Q / [(2π)^(3/2) · σx · σy · σz] · exp[-((x - xc)²) / (2σx²)] · exp[-(y²) / (2σy²)] · {exp[-((z - H)²) / (2σz²)] exp[-((z H)²) / (2σz²)]}其中Q是瞬时释放量xc是烟团中心顺风向位置H是释放高度。你看浓度直接除以三个σ的乘积。所以σ给错了浓度不是差一点是错一个量级。而对连续泄漏比如管线漏气、烟囱排放用的是烟羽模型。连续公式里只在y、z两个方向有σy、σz顺风向的σx不出现在公式里因为连续源的顺风向扩散相比平流输送可以忽略。这就是为什么很多老代码只算σy和σz而你如果要返回σx、σy、σz三个值基本可以断定做的是瞬时烟团模型。1.2 扩散参数决定了模拟的生死有人可能觉得既然就三个数那随便找个公式套一套不就完了真不是。扩散参数是高斯模型里最敏感、误差最容易放大的输入量。原因很简单浓度和σ的乘积成反比。假如σy和σz每个都大了30%浓度峰值就会变成原来的1/(1.3×1.3)≈0.59直接腰斩。反过来如果σ给小了算出来的浓度可能超出实际几倍应急决策时会把“轻伤事故”误判成“重大伤亡事故”。这个敏感性在实际项目中非常致命。我做化工园区风险评价时同一个泄漏场景用不同稳定度等级算出来下风向1000米处的最大浓度可以差5到10倍。而稳定度等级只是代码里的一个字符串参数容易传错。另外σ还决定了地面反射项的修正位置。烟团碰到地面后浓度分布不再是简单的高斯需要叠加一个镜像源。镜像源的位置和高度也要用到σz尤其是σz较小时地面反射对近地面浓度的影响非常明显不能忽略。1.3 为什么入参是t而不是直接给距离常见的扩散参数经验公式比如Briggs公式都是以“下风向距离x”作为自变量。那我们的函数为什么偏偏要输入扩散时间t因为t是烟团模型里更自然的量。第一瞬时释放的烟团在风场里运动中心位置随时间移动。模拟的时候系统是按时步推进的每个烟团从释放那一刻开始计时它的扩散时间天然就是t而不是x。第二真实风场不会恒定。烟团经过复杂地形风速可能变、风向可能转。如果要用x作参数得算“等效下风向距离”也就是沿轨迹积分风速这个积分在风场不规则时很麻烦。用t就简单t是拉格朗日时间烟团经历了多长时间就是多长的扩散时间。第三多烟团模型里每个烟团的释放时间不同、年龄也不同。这时候对每个烟团单独记录“出生时间”比记录“它走了多远”要可靠得多。远距离扩散时风场、湍流都在变化用时间轴统一管理所有烟团程序也更好写。所以这个函数的设计是合理的外界给的是扩散时间t函数内部再把t换算成下风向距离x然后走经验公式。2. 扩散参数的计算原理与模型选择2.1 从湍流统计理论看扩散本质要理解经验公式为什么长那样得先知道一点湍流扩散理论。单个粒子在湍流里运动它的位移方差随时间怎么增长这个问题的经典回答来自Taylor的统计理论。粒子在某个方向上的速度脉动有一个自相关函数拉格朗日积分时间尺度T_L描述了速度脉动“记住自己方向”的时间。在这个框架下扩散参数随时间有两个极限行为近距离t远小于T_L粒子基本按初始速度直线运动位移标准差近似正比于 tσ ≈ σ_v · t。远距离t远大于T_L粒子运动已经完全随机化和布朗运动一样位移方差近似线性增长σ ≈ sqrt(2 · K · t)。这里K是湍流扩散系数。中间区域比较复杂没有简单的解析解。理解这两个极限很重要。有些项目里直接用σ sqrt(2Kt)针对连续排放、长时间平均是对的但你要是用它去算刚释放几十秒的烟团会把扩散尺度算得明显偏小近源浓度偏高。而σ ∝ t的近场规律在泄漏发生后短时间内的应急评估里非常关键。高斯烟团参数化方案本质上就是用经验公式去近似这条完整的σ-t曲线。不同方案适应不同场景Briggs参数化是工程界用得最多的。2.2 Pasquill-Gifford与Briggs经验公式的坐标系大气扩散里最经典的一族经验方法是Pasquill-GiffordPG曲线。它把大气按稳定度分成A到F六档A是极不稳定、B是不稳定、C是弱不稳定、D是中性、E是较稳定、F是稳定。稳定度不同湍流强度差异很大扩散参数自然天差地别。Briggs在PG曲线基础上整理了便于计算的幂函数形式是目前工程模型使用的主流。最常用的Briggs扩散参数公式以x为下风向距离单位米σy、σz单位米分为乡村和城市两套系数。乡村条件下横向扩散参数σy的形式是统一的σy a · x · (1 b · x)^(-1/2)垂直扩散参数σz则分档A: σz 0.20 · x B: σz 0.12 · x C: σz 0.08 · x · (1 0.0002 · x)^(-1/2) D: σz 0.06 · x · (1 0.0015 · x)^(-1/2) E: σz 0.03 · x · (1 0.0003 · x)^(-1) F: σz 0.016 · x · (1 0.0003 · x)^(-1)城市条件下由于建筑引起的机械湍流更强整体扩散速度快A和B合并成A-BE和F合并成E-F横向公式变成σy a · x · (1 0.0004 · x)^(-1/2)其中A-B取0.32C取0.22D取0.16E-F取0.11。垂直方向也有对应的一套。你要注意这些公式的x是下风向距离最终我们在函数里用x u · t换算得到。下面这张表能帮你直观感受不同稳定度下扩散尺度的差异。取x1000米乡村条件稳定度σy米σz米直观感受A210200烟云又宽又高垂直混合剧烈B152120扩散快云团蓬松C10573中等状态D7638中性垂直扩散明显受限E5723云团扁平F3812烟云窄而矮几乎贴地看着这组数据你就能理解为什么夜间稳定条件下泄漏特别危险——σz只有12米污染物贴着地面不散近地面浓度远高于白天。2.3 t到x的换算风速u取值有讲究函数内部用x u · t把扩散时间换成下风向距离这里u取多少会直接影响计算结果。很多人随手拿一个气象站10米高度风速就用了这不一定对。扩散过程中的输运速度应该取泄漏源释放高度到烟羽中心高度之间的平均风速而不是地面风速。工程上常用幂律风廓线做修正u(z) u_ref · (z / z_ref)^p其中u_ref是参考高度处的风速一般取10米气象站风速。指数p和大气稳定度有关稳定度ABCDEF风廓线指数p0.070.140.210.330.440.54所以你看到这里有个隐含问题稳定度D时p0.33把地面风速换算到100米高处风速会大不少。如果源高是100米你用10米风速直接去算x相当于把扩散时间对应的下风向距离低估了后面所有σ都会被带偏。我一般会给函数再加一个参数传入有效输运风速或者传源高让函数内部做风廓线修正。这样接口更明确不会让大家为了省事直接填地面风速。3. Python实现t入参σx、σy、σz出参3.1 函数签名设计的几个决策写这个函数之前先确定外部传入哪些参数t扩散时间秒必填。stability稳定度等级推荐用字符串“A”“B”“C”“D”“E”“F”城市模型还支持“AB”“EF”。有人喜欢用parse到枚举可以字符串比较直白配合校验就行。terrain下垫面类型只有“rural”和“city”两种。默认用rural保守。wind_speed有效风速米/秒也就是前面说的输运风速不是可选的应该让调用方明确传入。sigma_min扩散参数下限值默认给一个很小的正数就行防止t0时出现除零问题。为什么不直接用API直接接受x因为题目既然定了t是输入而且烟团模型需要按年龄管理扩散尺度所以外部接口保持t。对连续烟羽的用户他可以自己把网格点坐标到源点的顺风向距离除以u变成等效t再调函数一样的。返回值顺序就按题目要求σx、σy、σz。3.2 可直接复用的Python实现下面是我整理的版本没有依赖第三方库纯标准库就能跑import math # 乡村条件 Briggs 系数: (ay, by, py, az, bz, pz) # σy ay * x * (1 by * x) ** py # σz az * x * (1 bz * x) ** pz _BRIGGS_RURAL { A: (0.22, 0.0001, -0.5, 0.20, 0.0, 0.0), B: (0.16, 0.0001, -0.5, 0.12, 0.0, 0.0), C: (0.11, 0.0001, -0.5, 0.08, 0.0002, -0.5), D: (0.08, 0.0001, -0.5, 0.06, 0.0015, -0.5), E: (0.06, 0.0001, -0.5, 0.03, 0.0003, -1.0), F: (0.04, 0.0001, -0.5, 0.016, 0.0003, -1.0), } # 城市条件 Briggs 系数 _BRIGGS_CITY { AB: (0.32, 0.0004, -0.5, 0.24, 0.0010, 0.5), C: (0.22, 0.0004, -0.5, 0.20, 0.0, 0.0), D: (0.16, 0.0004, -0.5, 0.14, 0.0003, -0.5), EF: (0.11, 0.0004, -0.5, 0.08, 0.0015, -0.5), } def diffusion_sigma(t, stabilityD, terrainrural, wind_speed2.0, sigma_min0.1): if t 0: raise ValueError(扩散时间不能为负数) if terrain not in (rural, city): raise ValueError(terrain 仅支持 rural 或 city) if terrain rural: if stability not in _BRIGGS_RURAL: raise ValueError(f乡村条件不支持稳定度: {stability}) ay, by, py, az, bz, pz _BRIGGS_RURAL[stability] else: if stability not in _BRIGGS_CITY: # 城市模型 A、B 合并为 ABE、F 合并为 EF if stability in (A, B): key AB elif stability in (E, F): key EF else: raise ValueError(f城市条件不支持稳定度: {stability}) else: key stability ay, by, py, az, bz, pz _BRIGGS_CITY[key] # 扩散时间转为下风向距离米 x wind_speed * t # 横向扩散工程上取 σx σy sigma_y ay * x * (1 by * x) ** py # 垂直扩散 sigma_z az * x * (1 bz * x) ** pz sigma_x sigma_y return max(sigma_x, sigma_min), \ max(sigma_y, sigma_min), \ max(sigma_z, sigma_min)调用方式很简单sx, sy, sz diffusion_sigma( t600, stabilityD, terrainrural, wind_speed2.0 ) print(sx, sy, sz)稳定度D、风速2米/秒、扩散10分钟后等效x1200米。算出来的σxσy≈90.7米σz≈43.0米。这个量级在工程上是合理的——10分钟烟团已经铺开大约几百米尺度垂直方向受限只有横向的一半。3.3 边界处理与数值稳定性细节代码很短但有几个细节是实际项目里反复遇到问题的地方。第一t0时的处理。如果什么都不做x0σ全部等于0浓度公式里除以σ立即出现inf或nan。现实里泄漏刚刚发生的瞬间烟团不可能是无限小源本身有热浮力、初始动量和羽流抬升初始尺寸不可能为零。所以我在返回值里用sigma_min做了下限保护。工程上更严谨的做法是给源设定一个初始扩散参数σ0比如1到5米然后把公式写成sqrt(σ_model² σ0²)这样既避免除零也符合“源有初始体积”的物理事实。第二适用距离限制。Briggs公式不是万能的很多文献明确它适用于下风向100米到10公里左右。距离太近时风和建筑尾流影响复杂太远时大气边界层对垂直扩散的制约会变得明显单纯的幂函数形式会高估σz。所以函数内部最好加一个warning的阈值if x 10000: print(fWARNING: x{x:.1f}m 超出Briggs公式常用适用范围(10km以内))不要在函数里直接抛异常因为有些模型评估还是会继续跑但要让调用者知道结果可信度在下降。第三稳定度输入校验。字符串作为参数输入最大的风险就是拼写错误。代码里我做了显式校验不支持的组合直接抛ValueError这比错误计算、返回一个明显不合理的结果要负责任得多。4. 验证、测试与工程落地的几个坑4.1 用典型场景验证结果是否合理函数写完了第一件事不是部署而是验证。没有验证条件的项目至少可以做量级验证。我的习惯是选用一组经典场景做基准。最常用的是乡村、D级稳定度、x1000米。按Briggs公式σy 0.08 × 1000 × (1 0.0001 × 1000)^(-0.5) ≈ 76.3米 σz 0.06 × 1000 × (1 0.0015 × 1000)^(-0.5) ≈ 37.9米这个结果的物理含义是中性层结条件下泄漏1公里后烟云的半宽约76米垂直厚度约38米。你去翻大气扩散实验数据类似条件下的横向扩散尺度一般就是几十米到百米的量级说明函数输出没有系统性偏差。再验证一个不稳定条件乡村、A级、x1000米。σy≈210米σz200米。白天强日照、小风天大气湍流旺盛垂直方向上下混合剧烈烟云变得很“胖”这是符合直觉的。还有一个极端稳定条件乡村、F级、x1000米。σy≈38米σz≈12米。夜间晴空、静稳烟云又窄又矮污染物基本趴在地面附近移动近地面浓度可以高到离谱。这也是为什么夜间重气泄漏事故更容易酿成大祸。这几个测试场景一次通过基本可以确认公式实现没写错。如果你算出来D级σz有几百米那肯定公式或者单位有误先回头查x的单位和风速单位。4.2 现实中影响σ的三个隐形因素工程落地时即使公式实现完全正确结果也可能和实测对不上。除了稳定度分级误差还有几个常见的隐形因素。第一个是粗糙度。Briggs只有乡村和城市两档但实际地表千差万别。农田、草原、丘陵、森林、高楼区粗糙度差异很大。城市档的σ比乡村档大不少因为建筑引起的机械湍流强。如果你的源在一个县城工业园区周边建筑不高用城市档可能高估扩散但完全用乡村档又会低估因为工业设备和厂房本身就是扰流器。我的建议是做敏感度分析分别用两套参数跑一遍看结果落在什么区间报告里写清楚模型假设。第二个是热力因素导致的不均匀性。比如泄漏源是高温气体或者比空气重高斯模型本身就不太适用。重气扩散比如液化气大规模泄漏初始阶段会在地面铺开密度分层抑制垂直扩散再用高斯烟团会严重低估近距离地面浓度。这时候要换SLAB、AFTOX或者类似重气模型。第三个是平均时间尺度。高斯模型的扩散参数和气象数据平均时间有关。PG曲线对应的取样时间大约是10分钟Briggs参数也基于类似量级。你要是输入的是1分钟平均风速或者3小时平均风速算出来的σ和实测会有系统性偏差。应急场景里气象数据更新频率尽量控制在10分钟左右不要拿日均风速来跑瞬时释放。4.3 时间参数在完整模拟中的正确理解最后说一个最容易踩的坑就是t在连续烟羽和瞬时烟团里的不同用法。连续烟羽模型里如果你用这个函数去算网格点浓度t不等于模拟总时长而等于“从源到计算点的输运时间”。比如源在原点计算点在下风向500米风速2米/秒那t应该取250秒而不是当前模拟到了第3600秒。很多新手把模拟时间直接当t传入导致远处的σ巨大浓度被稀释到几乎没有这完全是错用。而瞬时烟团模型里每个烟团的t是该烟团自身的年龄也就是从释放那一刻起经过的时间。如果场景是持续泄漏边释放边扩散你会不断产生新的烟团最老的烟团可能有3600秒年龄但新释放的烟团只有10秒年龄。对不同的烟团要分别拿年龄去调函数不能用一个全局t。用函数的正确姿势是先明确场景瞬时释放所有污染物一次性进入环境对所有位置用同一个时间源起点t是释放后经过的时间。连续释放把连续源离散成一系列微小烟团每个烟团有自己的释放时间浓度是所有烟团贡献的叠加。在我自己的项目里我习惯把烟团数据设计成这样class Puff: def __init__(self, x0, y0, z0, release_time, mass): self.x0 x0 self.y0 y0 self.z0 z0 self.release_time release_time self.mass mass def age(self, current_time): return max(current_time - self.release_time, 0.0)计算某个时刻浓度时遍历所有烟团用age调用diffusion_sigma得到每个烟团的σ再按烟团公式叠加。这样做程序结构清晰也方便以后扩展风场、沉降、衰变这些过程。4.4 稳定度等级的选择速查这里把工程里最常用的稳定度判断逻辑整理成一个速查表你可以直接贴在工位上。时间段日照/云量10米风速稳定度等级白天强日照小于2m/sA白天强日照2~3m/sB白天中等日照3~5m/sC白天弱日照任意D夜间云量较多任意D夜间少云2~3m/sE夜间晴空静稳小于2m/sF全天大风大于6m/s任意D注意这张表只是粗略判断真正的环境评价项目需要按标准方法结合太阳高度角、云量、风速选定稳定度。理论上稳定度等级还有更精细的划分超过A-F会引入G但工程应急领域用到A-F就足够了。我实际操作中的体会是别过度信任默认值。函数里把稳定度默认成D级虽然安全但D级只能代表中性状态你处理夜间辐射逆温场景时必须手动改成E或F。如果你懒得判断可以输入一组不稳定、中性、稳定三种稳定度跑出浓度区间给应急决策提供的是范围而不是单点值。收尾一点个人经验这个函数我前后改了很多版最开始的版本只返回σy和σz后来加烟团模型才补上σx。过程中踩过最大的坑是参数表的错误引用——不同文献里Briggs系数的x单位有按米写的也有按公里写的抄错一次计算结果差了好几个量级折腾了整整一天才定位到问题。所以现在我的习惯是任何扩散参数函数写完之后先跑一组手工算过的基准值再进模型。你说函数简单也好说它就二三行也好扩散模型对参数误差的放大效应真的不能小看。接下来你可以在这个基础上继续扩展比如把Moreira、Degrazia这些更适合长时间扩散的参数化方案加进去做成可切换的模式或者接上风场模块让每个烟团的风速随位置变化。扩散模拟的深度没有天花板先把最底层的扩散参数算稳上面的大楼才能盖得住。
