企业官网建设流程全解析

1. 从扩散时间到扩散参数：一个函数背后的真实物理场景

先把这个需求翻译成人话：你给我一个以秒为单位的扩散时间 t，我返回三个方向的扩散标准差 σx、σy、σz。这通常是大气扩散模型、污染物泄漏模拟、尾气扩散评估或者粒子追踪程序里的一个核心函数。我之前在做工业厂区泄漏风险评估时，就经常要跟这类参数打交道。

为什么要按时间算扩散参数？因为污染物从释放源出来以后，不是老老实实待在一个地方的，它会随着空气湍流在水平方向和垂直方向不断铺开。这个"铺开的程度"在数学上就体现为扩散标准差 σ，单位是米。t 越大，烟团扩散的范围越宽，σ 自然就越大。问题是，这个增长不是简单的线性关系，它跟大气稳定度、地面粗糙度、下风向距离甚至采样时间都有关系，所以才会需要专门的函数来计算。

在实际工程里，这个函数通常服务于两类场景。第一类是高斯烟团模型，比如突然发生的一次性泄漏，污染物像"一口气"那样被释放出来，然后随风扩散，这时候每个时刻的烟团都用一组 σ 来描述它的形态。第二类是拉格朗日粒子扩散模型，大量粒子被释放到流场里，每一步移动都需要根据当地的湍流强度给定一个随机位移，而这个位移幅度就跟 σ 的增长率直接相关。无论哪类场景，σ 的计算精度直接决定了落地浓度预测靠不靠谱。

需要明确的是：σx、σy、σz 分别代表纵向（风向方向）、横向（垂直于风向的水平方向）和垂直方向的扩散尺度。在一些简化模型里 σx 常被忽略或与 σy 相等，但真实的三维场景下三者并不相同，尤其是垂直方向受大气稳定度影响最明显，近地面时还会被地面压缩。所以，别小看这个"返回三个值"的函数，它背后是把大气边界层物理学浓缩成了一行行代码。

注意：这里说的 t 是"扩散时间"，不是"源释放开始之后的墙钟时间"。如果存在风场平流，扩散时间应该是粒子/烟团离开源后经历的时间，二者要区分开。很多新手在这个地方栽过跟头，算出来的 σ 偏大或偏小都找不出原因。

2. 为什么 σ 随时间是"非均匀"增长的

2.1 湍流扩散的物理本质：从分子扩散到涡旋输送

要理解 σ 的时变规律，得先知道扩散的物理机制。分子扩散是布朗运动式的，σ² 与时间成正比，即 σ ∝ √t，这是最简单的费克扩散。但大气边界层里起主导作用的不是分子扩散，而是湍流扩散——由大大小小的涡旋把污染物"搅"开。

如果你站在一个烟囱下风向看烟气，会发现烟羽一会儿向左偏、一会儿向右偏，整体呈扇形展开。这个"展开"的速度取决于边界层里的湍流强度，而湍流强度又和风速、日照、云量、地表状况都有关系。泰勒（Taylor）在1921年提出的单粒子扩散统计理论给出了一个关键结论：在扩散时间很短时，粒子运动是"惯性保持"的，σ 与 t 近似线性增长；在扩散时间足够长后，粒子彻底"忘记"了初始速度，σ 才转为与 √t 成正比。

也就是说，σ 随时间的增长经历了"线性段 → 过渡段 → 平方根段"三个阶段。如果你用一个简单的 σ = a·t 或者 σ = a·√t 去硬套全场，误差会很大。而实际工程里我们要的往往就是覆盖从近源几十米到远源几十公里的结果，所以必须分段处理。

2.2 三个方向为什么不一样

大气湍流在三个方向上是各向异性的。水平方向的湍流主要受大尺度涡旋和地形扰动影响，强度大；垂直方向则受大气层结（稳定度）强烈约束——稳定层结下垂直湍流被压制，不稳定层结下热泡翻腾，垂直扩散显著增强。

具体到参数化方案里：

• σy（横向）主要由水平湍流决定，往往可以近似看成随风速和距离稳定增长，受稳定度影响相对温和；
• σz（垂直）受稳定度影响最大，稳定条件下它可能只有 σy 的几分之一，不稳定条件下两者可能接近甚至 σz 反超；
• σx（纵向）在烟团模型中通常按与 σy 类似的规律处理，但在存在风切变时纵向拉伸效应会很明显。

所以一个合格的"t 转 σ"函数，核心不是套一个万能公式，而是正确处理各向异性和稳定度分段这两件事。

3. 三种主流的工程实现方案对比

要想把"t 秒 → σx, σy, σz"做进代码里，业内其实有好几条路可以走，每条路的精度、适用范围和坑都不一样。我按工程中见到的频率逐一拆开说。

3.1 方案一：Pasquill–Gifford 曲线数字化（最常用，适合中小尺度）

这是帕斯奎尔在1961年基于大量野外扩散实验（著名的草原实验等）总结的经典方法。它按气象条件把大气稳定度分成 A-F 六个等级（A 极不稳定，D 中性，F 极稳定），然后给每个等级配套一组 σy(x) 和 σz(x) 随下风向距离变化的经验曲线。原始数据是曲线图，后来很多人把它拟合成幂函数表达式，变成代码可以直接用的形式。

实现时，你需要先输入风速、日照、云量等气象数据去判定稳定度等级，再把"扩散时间 t"换算成"下风向距离 x"：x = u·t（u 是平均风速）。然后代入对应等级的经验公式去算 σ。这组公式并不复杂，网上流传的 Martin-Turner 拟合、Briggs 拟合都属此类。

优点：简单、标准、被环保部门广泛认可，写报告容易溯源。 缺点：只适用于相对均匀的地形和下垫面，对复杂地形或城市街区无能为力；而且它本身就是从中小尺度实验拟合出来的，超过几十公里后外推可信度明显下降。

3.2 方案二：Briggs 公式（适合开阔乡村与城市，分段明确）

Briggs 在 1973 年前后对 Pasquill 曲线做了改进，给出了一组更便于计算的解析公式，并且分"乡村"和"城市"两套参数。它的自变量同样是下风向距离 x，曲线形态按距离分段（比如 0.1km 以下、0.1-1km、1km 以上），每一段用不同的系数。

我用 Briggs 公式比较多，原因是它在表达式上非常规整，非常适合写进通用函数库里。而且它对 A-F 稳定度全覆盖，参数表查起来方便，特别适合做批量模拟。需要注意：Briggs 公式给的是 σy 和 σz 随 x 的变化，如果你需要在烟团模型里使用，还需要解决"距离 x 怎么对应扩散时间 t"的问题，这个我放在代码实现里细讲。

3.3 方案三：泰勒统计理论的 Lagrangrian 时间尺度法（适合粒子模型）

如果你做的是拉格朗日粒子扩散模型（比如用大涡模拟或诊断风场驱动粒子运动），那么上面的高斯烟羽系经验公式就不够"物理"了。更合理的做法是直接用泰勒扩散统计理论：

σ_i²(t) = 2·σ_vi²·T_Li²·(t/T_Li - 1 + e^(-t/T_Li))

其中 i 代表 x、y、z 三个方向；σ_vi 是第 i 方向脉动风速的标准差；T_Li 是第 i 方向的拉格朗日时间尺度，代表粒子速度"记忆"持续的时间尺度。

这个公式的物理图像非常清楚：t 远小于 T_L 时括号内近似为 t²/(2T_L²)，σ ∝ t，粒子还没忘掉初速度；t 远大于 T_L 时括号里近似为 t/T_L - 1，σ ∝ √t，扩散进入正常扩散区。中间过渡段自动平滑衔接，不需要人为分段，这是它最大的优点。

脉动风速标准差 σ_vi 通常用莫宁-奥布霍夫相似理论（Monin-Obukhov Similarity Theory）来估算，需要输入摩擦速度 u*、奥布霍夫长度 L、边界层高度 zi 等变量；T_Li 也有对应的经验公式。这条路参数多，计算量也大，但物理一致性最好，尤其适合非均匀、非稳态的复杂流场。

简单总结三套方案的选型逻辑：做环评等级扩散计算用方案一；写通用计算工具库、需要快且稳的用方案二；做科研级粒子扩散模拟用方案三。如果只是要一个够用且不引入额外气象参数的函数，方案一和方案二都是好选择，下面我给出完整的工程实现。

4. 手把手实现：从理论公式到可直接调用的Python函数

4.1 第一步：确定稳定度等级

稳定度等级是高斯型扩散模型里最重要的输入。工程上最常用的是 Turner 方法，利用地面风速、太阳辐射等级或云量来判定。为了让你快速跑通，我先把稳定度判定的简化版代码写出来。这里用的是"风速+日照强度"的简化规则，适合自动化程序里没有人工观天时使用。

def stability_class(u10, solar_radiation): """ 根据10米风速(m/s)和太阳辐射等级(0=夜间, 1=弱, 2=中, 3=强)返回稳定度等级字符串。 简化版Turner法，仅供自动化程序用，正式环评请用完整观测数据。 """ if solar_radiation == 0: # 夜间 if u10 < 2.0: return 'F' elif u10 < 3.0: return 'E' elif u10 < 5.0: return 'D' else: return 'D' else: # 白天 if solar_radiation == 1: # 弱日照 if u10 < 2.0: return 'A' elif u10 < 3.0: return 'B' elif u10 < 5.0: return 'C' else: return 'D' elif solar_radiation == 2: # 中日照 if u10 < 2.0: return 'A' elif u10 < 3.0: return 'B' elif u10 < 5.0: return 'C' else: return 'D' else: # 强日照 if u10 < 2.0: return 'A' elif u10 < 3.0: return 'B' elif u10 < 5.0: return 'C' else: return 'D'

这个函数虽然简陋，但结构清晰，你完全可以根据项目所在地区的实际气候特征把阈值调得更合理。稳定度判定的精度对最终 σ 影响极大，我试过同一组风速数据，稳定度从 D 变到 F，σz 可能差 3 到 5 倍，落地浓度直接差一个量级，不可大意。

4.2 第二步：在"时间域"直接计算 σ 的 Briggs 连续化形式

接下来就是核心问题：很多模型框架里只有时间 t，没有距离 x。标准的 Briggs 公式是以 x（km）为自变量的经验式，如果你非要把它直接用在时间域，就得在每一小步积分里做"距离-时间同步"处理。我给出一种工程上常用的同步方案，思路是：

1. 以当前时刻 t，用上一次算出的 σ 反推一个"等效下风向距离" x_eff；
2. 用 u·t 得到一个名义距离；
3. 把名义距离代入 Briggs 公式得到这一时刻的 σ。

更简洁的做法是，在模型的时间积分循环里干脆走增量式：

def compute_sigma_t_briggs(t, u, stability): """ Briggs改进式，时间域同步方案。 t: 扩散时间，秒 u: 平均风速，m/s stability: 'A','B','C','D','E','F' 返回 (sigma_x, sigma_y, sigma_z)，单位：米 """ x_km = u * t / 1000.0 # 扩散时间对应的等效距离(km) x_m = u * t # 以米为单位备用 # ---- Briggs 开阔乡村系数表 ---- # 每一项是 (sigma_y系数, sigma_y指数, sigma_z系数1, 指数1, sigma_z系数2, 指数2) # 不同稳定度各不同，此处列出常用A、B、C、D、E、F简化拟合 if stability == 'A': # sigma_y = 0.22*x_km*(1+0.0001*x_km)^-0.5 ; sigma_z = 0.20*x_km sigma_y = 0.22 * x_km * (1.0 + 0.0001 * x_km) ** (-0.5) * 1000.0 sigma_z = 0.20 * x_km * 1000.0 elif stability == 'B': sigma_y = 0.16 * x_km * (1.0 + 0.0001 * x_km) ** (-0.5) * 1000.0 sigma_z = 0.12 * x_km * 1000.0 elif stability == 'C': sigma_y = 0.11 * x_km * (1.0 + 0.0001 * x_km) ** (-0.5) * 1000.0 sigma_z = 0.08 * x_km * (1.0 + 0.0002 * x_km) ** (-0.5) * 1000.0 elif stability == 'D': sigma_y = 0.08 * x_km * (1.0 + 0.0001 * x_km) ** (-0.5) * 1000.0 sigma_z = 0.06 * x_km * (1.0 + 0.0015 * x_km) ** (-0.5) * 1000.0 elif stability == 'E': sigma_y = 0.06 * x_km * (1.0 + 0.0001 * x_km) ** (-0.5) * 1000.0 sigma_z = 0.03 * x_km * (1.0 + 0.0003 * x_km) ** (-1.0) * 1000.0 else: # F sigma_y = 0.04 * x_km * (1.0 + 0.0001 * x_km) ** (-0.5) * 1000.0 sigma_z = 0.016 * x_km * (1.0 + 0.0003 * x_km) ** (-1.0) * 1000.0 # 纵向sigma_x 在本方案中按sigma_y的比例取1.0倍(即水平各向同性) # 如果要考虑风切变的纵向拉伸，可以给一个大于1的修正系数 sigma_x = sigma_y return sigma_x, sigma_y, sigma_z

这段代码最核心的思想是：在时间域模型里，用 u·t 将时间映射到等效距离，再套用距离域的经验公式。它隐含的假设是风速恒定、风向平稳，且烟团质心以风速匀速平流。这个假设在大部分中小尺度模拟里是可接受的，但你要清楚它在强切变或非平稳风场里会失真。

4.3 第三步：基于泰勒统计理论的更物理实现

如果你在写粒子扩散模型，上面这种"只依赖距离"的经验法就不够用了，因为粒子每时每刻受到的湍流作用都和你模拟的风场、温度层结密切相关。这时候更推荐直接用 Taylor 公式，一步到位从 t 算 σ，不需要先转距离。

import math def compute_sigma_taylor(t, sigma_vx, sigma_vy, sigma_vz, T_Lx, T_Ly, T_Lz): """ 基于Taylor单粒子扩散统计理论的sigma(t)计算。 t: 扩散时间, 秒 sigma_vx, sigma_vy, sigma_vz: 三个方向的脉动风速标准差, m/s T_Lx, T_Ly, T_Lz: 三个方向的拉格朗日时间尺度, 秒 返回 (sigma_x, sigma_y, sigma_z), 单位: 米 """ if t <= 0: return 0.0, 0.0, 0.0 def _sigma2(sigma_v, T_L): # sigma^2 = 2 * sigma_v^2 * T_L^2 * (t/T_L - 1 + exp(-t/T_L)) tau = t / T_L if T_L > 0 else 0.0 if tau <= 0: return 0.0 return 2.0 * sigma_v * sigma_v * T_L * T_L * (tau - 1.0 + math.exp(-tau)) sx = math.sqrt(max(_sigma2(sigma_vx, T_Lx), 0.0)) sy = math.sqrt(max(_sigma2(sigma_vy, T_Ly), 0.0)) sz = math.sqrt(max(_sigma2(sigma_vz, T_Lz), 0.0)) return sx, sy, sz

这个函数短小精悍，但它比经验公式"诚实"得多：它把三个方向的湍流强度和时间尺度交给调用方来决定。那 σ_v 和 T_L 怎么来？最常见的做法是在模型初始化时，根据 Monin-Obukhov 相似理论估算边界层湍流参数，然后算出各方向的 σ_v 和 T_L，这个计算量不小，但在粒子模型里属于一次性成本。

这里有一个非常容易踩的坑：σ_v 的单位和量级。σ_v 是脉动风速的标准差，量级通常在 0.2 到 2 m/s 之间，别和平均风速搞混。T_L 的量级在近地面通常是几十秒，在中性层的上部可以达到几百秒。如果你代码里 T_L 填了 1 秒，那 σ 会长期处于线性增长区，扩散尺度会被严重低估；如果填了 10000 秒，那 σ 会非常快地进入 √t 段，结果又会偏大。所以这个方案虽然物理更严谨，但对输入参数敏感性高，一定要做好气象前处理。

4.4 第四步：代入真实场面的效果对比

我拿一个典型的中性层结（D 类）场景做个计算对比：风速 u=5 m/s，扩散时间 t=600 秒（10分钟），Briggs 方案：等效距离 x=3km，算得 σy≈80m，σz≈40m 上下（取决于具体系数表细节）。Taylor 方案如果取 σ_v=0.5 m/s、T_L=100s，则 σy≈0.5×sqrt(2×0.5²×100²×(...))≈上百米，量级是基本吻合的。这说明在参数合理的情况下，两套方案结果能互相印证。

但一旦进入强不稳定（A 类）或强稳定（F 类）条件，两套方案的差异会被拉大。尤其 F 类稳定条件下，近地面垂直扩散被强烈抑制，σz 可能只有几米，这时如果你的模型里 σz 算出了几十米，那基本可以断定气象输入或者稳定度判定出了问题。

5. 实操中的三个经典坑与排查技巧

这个函数看起来只有几行，实际用起来问题不少。我把这些年遇到的高频问题整理成速查表，可以帮你快速定位。

一个我特别想强调的坑：在时间步进模型中，σ 不应该被反复重新初始化。很多人把计算 σ 的代码放在每个时间步内，每次都用当前累计 t 去算，这没问题；但如果你用了"增量式 σ += Δσ"的写法，就要注意 Δσ 必须根据当前 σ 所在的增长阶段来给，不能用线性增长近似。我见过一个案例，程序把 σ 当成了随时间均匀增长的量，每步加一个固定值，短时间模拟还行，时间一长 σ 直接爆表，落地浓度负值都出来了。

另一个经验：务必强制加一个最小扩散尺度下限。数值模型中，σ 在 t=0 附近趋向于 0，这会导致扩散系数趋近于 0，模型数值上容易产生振荡甚至发散。常规做法是给 σ 设一个下限，例如 σ_min = 0.1 米，或者把初始 t 设为一个大一点的等效时间（如 1 秒），保证起步阶段数值稳定。

6. 再进一步：从均一湍流到非均匀边界层的扩展思路

前面所有公式都隐含假设整个扩散路径上的湍流是均匀的、定常的。这在地形平坦、气象平稳的假设下还能成立，但真实环境里，从烟囱排出的烟羽要穿过不同高度的风切变、温度层结变化，甚至在混合层顶被"盖帽"反射。这时 σ 的时变函数就不能再靠一个全局统一的公式了。

工程上有几个渐进式改进方案：

第一，把边界层按高度分层，每一层用不同的 σ_v、T_L 或 Briggs 系数，粒子穿过层界面时重新计算当地参数。这个方法实现简单，适合中小尺度精细模拟。

第二，用随机游走模型中的"反射边界"条件来处理混合层顶和地面的限制。地面和混合层顶的反射会让 σz 的实际增长比自由扩散慢，这种效应在近源（t 小）时不明显，在远源（t 大）时能显著影响地面浓度。

第三，从诊断风场数据里实时提取脉动速度标准差，替代固定的相似理论估算。这个方法对数据要求高，但结果最接近真实演进，适合做事故应急模拟。

我在实际项目中通常这样组织代码：核心扩散模块提供一个 compute_sigma(t, u, stability, terrain) 接口，内部根据场景选择 Briggs、Pasquill 或 Taylor 算法；对于三维粒子模式，则封装一个 TurbulenceField 类，按空间位置插值输出当地的 σ_v 和 T_L，再调用 Taylor 函数。这样接口统一，底层算法可替换，业务代码不用大改。

回到最初的需求："t 是扩散时间（秒），返回 σx, σy, σz"。这个函数固然短小，但真正让它在工程里立得住脚的，是它背后那套稳定度判定、参数选型、范围约束和数值稳定处理。如果你只是要一个能出数的函数，上面代码直接抄回去就能跑；如果你想在大项目里把它用得稳、用得准，建议把稳定度判定、参数表校验和边界条件处理这三件事一并做扎实。我自己的体会是，这种"小函数"反而是整个扩散模拟链条里最值得花时间打磨的环节——源头参数错了，后面再精细的流场和化学反应模块都是白搭。