1. 项目概述从现实问题到数学模型空气质量问题离我们并不遥远无论是城市里偶尔出现的雾霾天还是工厂周边居民对排放的担忧本质上都是空气污染物在特定气象和地理条件下的扩散与累积。作为一名长期与数据和模型打交道的从业者我处理过不少环境相关的项目。我发现很多人对“空气质量模拟”抱有敬畏之心觉得它高深莫测涉及复杂的流体力学和化学方程。但实际上它的核心思想非常直观搞清楚污染源排放了多少东西源强这些东西在风中怎么跑扩散以及跑的过程中会发生什么变化转化与沉降。数学建模就是为这个直观的过程建立一套可计算、可预测的数学语言。“空气质量模拟的数学建模与实战案例”这个标题精准地概括了从理论到实践的全链路。它适合环境科学、大气科学专业的学生深化理解适合从事环境评估、城市规划的工程师解决实际问题也适合对数学建模感兴趣、想找一个有现实意义的课题来练手的朋友。无论你是想弄明白一篇学术论文里的模型还是需要为自己的项目评估环境影响这篇文章都将带你走完从基本原理到代码实操的完整过程。我们将以最经典的高斯扩散模型作为切入点因为它结构清晰、物理意义明确是绝大多数复杂模型的基石。随后我们会用一个完整的实战案例手把手教你如何用Matlab将模型“跑”起来并分析结果。2. 核心模型解析高斯扩散模型及其数学内涵空气污染物在大气中的扩散可以类比为一滴墨水在静水中的晕染。但大气不是静水它有风有湍流情况复杂得多。高斯模型做了一个非常聪明且有效的简化它假设污染物在水平和垂直方向的浓度分布都服从我们熟悉的正态分布高斯分布。这个假设基于统计学中的“中心极限定理”即大量随机、微小的湍流运动共同作用的结果会使污染物粒子的位置分布趋近于正态分布。这是整个模型的灵魂所在。2.1 模型公式与参数物理意义对于一个连续排放的点源比如一根烟囱在下风向任意一点(x, y, z)处的污染物浓度C可以用以下公式表示C(x, y, z) (Q / (2π * u * σy * σz)) * exp(-y²/(2σy²)) * [exp(-(z-H)²/(2σz²)) exp(-(zH)²/(2σz²))]初看有点吓人但我们拆开看每一个部分都有明确的物理意义Q(源强)污染源单位时间内排放的污染物的质量单位通常是克/秒或千克/小时。这是整个模型的“输入”是所有计算的起点。如果源强估错了后续所有模拟都是空中楼阁。u(平均风速)排放口高度处的平均风速。风是污染物输送的“传送带”风速越大稀释作用越强下风向浓度越低。σy,σz(水平和垂直扩散参数)这是模型的核心参数表征了湍流扩散的强弱。它们不是常数而是随着下风向距离x的增加而增大的函数。σy和σz越大表示污染物在横向和垂直方向散得越开中心浓度就越低。它们的取值依赖于大气的稳定度是平静的、中性的还是湍流强烈的和地面粗糙度。通常使用 Briggs 或 Pasquill-Gifford 经验公式来查表或计算。H(有效排放高度)这不是烟囱的物理高度而是烟囱高度加上烟气因初始动量和热浮力产生的抬升高度。抬升高度计算本身就是一个子模型对于热烟气来说非常重要因为它直接影响了污染物在垂直方向的初始分布中心。公式后半部分的两个指数项第一个exp(-(z-H)²/(2σz²))代表从有效高度H处扩散的烟羽。第二个exp(-(zH)²/(2σz²))代表烟羽接触到地面后的反射作用假设地面完全反射污染物。这个镜像项保证了地面处的质量守恒是模型能准确模拟近地面浓度的关键。注意这个公式是连续点源稳态模型。它假设风场稳定、源强恒定、地形平坦并且污染物在输送过程中没有发生化学反应或沉降惰性气体。这些假设是它的局限性也是我们选择它作为入门的原因——先理解理想情况再逐步增加复杂度。2.2 大气稳定度判定模型参数的钥匙前面提到σy和σz依赖于大气稳定度。如何判定稳定度最常用的方法是Pasquill-Gifford (P-G) 分类法。它根据地面风速、日间太阳辐射强度或夜间云量将大气稳定度分为 A-F 六个等级A (极不稳定)晴朗夏日午后微风。B (不稳定)晴朗日中。C (略不稳定)多云白天或晴朗夜晚有微风。D (中性)阴天全天或夜间风速较大时。这是最常见也最“标准”的状态。E (略稳定)晴朗夜晚微风。F (稳定)晴朗夜晚静风。有了稳定度等级和距离x我们就可以通过查 P-G 曲线表或拟合的经验公式得到对应的σy(x)和σz(x)。在编程实现时我们通常会将查表数据拟合成幂函数形式如σy a * x^b方便计算。3. 实战案例模拟工业园区烟囱对下风向的影响理论说得再多不如亲手算一遍。我们设计一个典型的应用场景某工业园区内有一座燃煤锅炉的烟囱我们需要评估其排放的二氧化硫SO₂对下风向敏感点如居民区的浓度贡献。3.1 案例场景与参数设定假设我们获得以下基础数据源参数烟囱几何高度Hs 80 m烟囱出口内径d 2.5 m烟气排放速度Vs 15 m/s烟气温度Ts 150 °C环境气温Ta 20 °CSO₂ 排放速率源强Q 200 g/s气象参数烟囱出口高度处风速u 4.5 m/s环境条件晴朗秋日下午风速适中。根据 P-G 分类判定为B类不稳定大气。受体点我们关心从烟囱下风向 100m 到 5000m地面z0以及中心线y0上的浓度分布。同时重点关注下风向 2000m 处的一个具体受体点。3.2 关键计算步骤分解整个模拟流程可以分解为以下几个关键步骤我们会在Matlab中逐一实现。步骤一计算烟气抬升高度ΔH与有效源高H对于有热浮力和动量的烟羽常用Briggs 公式计算抬升。在不稳定B类大气中浮力抬升占主导。 首先计算烟气热释放率Qh粗略估算Qh ≈ 0.35 * (烟气与环境密度差相关的项) * Vs * d² * (Ts - Ta)/Ts更实用的工程简化是使用国标或行业导则中的公式。这里我们采用一个常见的 Briggs 浮力抬升公式ΔH 1.6 * Fb^(1/3) * x^(2/3) / u其中Fb是浮力通量Fb g * Vs * (d/2)² * (Ts - Ta)/Tsg为重力加速度。 但x是达到最终抬升的距离我们需要迭代或使用最终抬升公式。一个更直接的最终抬升公式为ΔH 1.6 * Fb^(1/3) * (3.5 * x*)^(2/3) / u其中x*是达到最终抬升的下风向距离对于不稳定大气x*取49 * Fb^(5/8)和实际关心距离的较小值。 计算过程略繁但Matlab编程可以轻松处理。假设我们计算出ΔH ≈ 45 m。 则有效源高H Hs ΔH 80 45 125 m。步骤二确定扩散参数 σy 和 σz对于 P-G B类稳定度我们可以使用 Briggs 给出的城市扩散参数公式幂函数形式σy 0.32 * x * (1 0.0004 * x)^(-0.5)单位mσz 0.24 * x * (1 0.001 * x)^(0.5)单位m 其中x是下风向距离米。这个公式在数公里范围内有较好的精度。我们将它直接写入Matlab函数。步骤三实现高斯模型浓度计算函数这是核心代码块。我们将编写一个Matlab函数gaussian_plume输入参数(x, y, z, Q, u, H, stability_class)返回浓度C。function C gaussian_plume(x, y, z, Q, u, H, stability) % 根据稳定度类别和距离x计算扩散参数sigma_y和sigma_z [sigma_y, sigma_z] calc_sigma(x, stability); % 高斯模型公式实现 term1 Q / (2 * pi * u * sigma_y * sigma_z); term2 exp(-0.5 * (y / sigma_y).^2); % 考虑地面反射的两个指数项 term3 exp(-0.5 * ((z - H) ./ sigma_z).^2) exp(-0.5 * ((z H) ./ sigma_z).^2); C term1 .* term2 .* term3; end function [sigma_y, sigma_z] calc_sigma(x, stability) % 以B类稳定度为例使用Briggs城市参数化方案 switch stability case B % 不稳定 sigma_y 0.32 * x .* (1 0.0004 * x).^(-0.5); sigma_z 0.24 * x .* (1 0.001 * x).^(0.5); case D % 中性 sigma_y 0.22 * x .* (1 0.0004 * x).^(-0.5); sigma_z 0.16 * x .* (1 0.001 * x).^(0.5); % 可以添加其他稳定度类别的参数化公式 otherwise error(Stability class not supported.); end end步骤四空间网格化计算与可视化为了看到整个污染烟羽的形态我们需要在x-y平面上创建一个网格计算每个网格点上的浓度并绘制等值线图或三维曲面图。% 定义计算范围 x_vec linspace(100, 5000, 100); % 下风向距离100个点 y_vec linspace(-500, 500, 80); % 横向距离80个点 z_receptor 1.5; % 受体高度通常取人的呼吸带高度1.5-2米 % 创建网格 [X, Y] meshgrid(x_vec, y_vec); C_grid zeros(size(X)); % 参数设定使用之前案例的数据 Q 200; % g/s u 4.5; % m/s H 125; % m stability B; % 遍历网格点计算浓度向量化操作效率更高 for i 1:numel(X) C_grid(i) gaussian_plume(X(i), Y(i), z_receptor, Q, u, H, stability); end % 绘制地面浓度等值线图 figure(Position, [100, 100, 800, 600]); contourf(X, Y, C_grid, 30, LineStyle, none); % 30条填充等值线无线条 colorbar; colormap(jet); % 使用jet色图直观显示浓度高低 xlabel(下风向距离 (m)); ylabel(横向距离 (m)); title(SO2地面浓度分布 (\mug/m^3) - B类稳定度); hold on; plot([0, max(x_vec)], [0, 0], k--, LineWidth, 1.5); % 画出中心线步骤五特定受体点浓度计算与评估现在我们来计算下风向2000米、中心线y0处地面z1.5m的浓度。x_target 2000; y_target 0; z_target 1.5; C_target gaussian_plume(x_target, y_target, z_target, Q, u, H, stability); fprintf(在下风向%d米中心线地面处的SO2浓度为%.2f μg/m³\n, x_target, C_target*1e6); % 注意我们模型中的Q单位是g/s计算出的C单位是g/m³乘以1e6得到μg/m³。假设计算结果是C_target 45.67 μg/m³。我们需要将这个结果与环境空气质量标准进行对比。例如中国《环境空气质量标准》(GB 3095-2012)中SO2的1小时平均浓度一级标准为150 μg/m³二级标准为500 μg/m³。45.67 μg/m³远低于标准限值表明在该特定气象条件下该单一源对2000米处的影响在可接受范围内。实操心得模型计算出的浓度是一次浓度即污染物直接扩散后的浓度。在实际环境评估中还需要考虑背景浓度、其他污染源的叠加以及化学转化如SO2转化为硫酸盐颗粒物。高斯模型是评估单个源贡献的利器但做总体环境评估时需要更复杂的模型或进行多源叠加。4. 模型进阶从理想走向现实基础的高斯模型解决了“有没有”的问题但要回答“准不准”我们必须考虑更多现实因素。这部分是建模工作的深化也是体现专业性的地方。4.1 复杂地形与建筑物的处理平坦地形假设在山区或城市中会失效。常用的处理方法是地形修正使用有效源高H减去受体点地形高度h_t(x,y)即H_eff H - h_t。如果烟羽低于地面H_eff 0则认为该点浓度为0烟羽被山体阻挡。这被称为“烟羽路径”模型虽然粗糙但实用。建筑物下洗当烟囱高度低于附近建筑物高度的2.5倍时烟气可能被卷吸到建筑物背风面的涡流区导致地面浓度急剧升高。处理这种情况需要更专业的计算流体力学CFD模型如AERMOD、ADMS等法规模型内置了相关算法。在简单评估中一个保守的做法是直接将烟囱高度视为0地面源来估算最大可能影响。4.2 非稳态与化学转化风速风向变化高斯稳态模型假设风向风速恒定。现实中风向会摆动这会导致横向扩散增强。一个经验方法是适当增大σy的值或者采用分段稳态模拟将长时间序列的风场数据输入进行逐时计算后再取平均或统计分布。干湿沉降颗粒物或可溶性气体会因重力或雨水冲刷从大气中移除。可以在高斯模型公式中增加一个衰减项exp(-λ * x/u)其中λ是沉降系数与污染物性质和气象条件有关。化学转化像SO2转化为硫酸盐是一个相对缓慢的过程时间尺度数小时到数天。对于近场模拟几公里内通常可以忽略。对于区域尺度模拟则需要耦合箱式化学模型或使用分段线性衰减来近似。4.3 面源与线源的处理实际污染源除了点源烟囱还有面源整个厂区无组织排放和线源公路机动车排放。面源可以将面源划分为多个小点源的集合或者使用虚拟点源法。虚拟点源法的思路是假设污染物从面源中心的上风向某个“虚拟点”释放经过一段初始距离x0的扩散后其扩散参数刚好等于面源初始的尺度。这样就把面源问题转化成了一个位于上风向的等效点源问题。线源对于公路线源通常将其离散化为一系列紧密排列的点源然后对每个点源在下风向受体点的贡献进行积分或求和。有专门的高斯线源模型公式但原理相通。5. 在Matlab中构建完整的模拟与分析工作流一个完整的空气质量模拟项目不仅仅是调用一个函数。它应该是一个可重复、可调整、可分析的工作流。下面我们构建一个更健壮的Matlab脚本框架。5.1 数据准备与参数管理将所有输入参数源、气象、受体组织在结构体或表格中便于管理和修改。% 定义源参数结构体 source.Hs 80; % 烟囱高度 (m) source.d 2.5; % 出口直径 (m) source.Vs 15; % 出口流速 (m/s) source.Ts 150; % 烟气温度 (°C) source.Q 200e-3; % 源强转换为 kg/s (200 g/s - 0.2 kg/s) % 定义气象参数结构体 meteo.u 4.5; % 风速 (m/s) meteo.Ta 20; % 环境温度 (°C) meteo.stability B; % 稳定度等级 % 可以扩展如加入风向、湿度等 % 定义受体网格 receptor.x_range [100, 5000]; receptor.y_range [-500, 500]; receptor.z 1.5; % 受体高度 receptor.nx 100; receptor.ny 80; % 计算有效源高调用一个独立的函数 source.H_eff calculate_effective_height(source, meteo);5.2 批处理与情景分析我们常常需要分析不同气象条件或排放情景下的结果。这可以通过循环来实现。% 定义不同的风速情景 wind_speeds [2.0, 4.5, 7.0]; % m/s % 定义不同的稳定度情景 stability_classes {A, B, D, F}; max_concentrations zeros(length(wind_speeds), length(stability_classes)); for i 1:length(wind_speeds) for j 1:length(stability_classes) meteo_current meteo; meteo_current.u wind_speeds(i); meteo_current.stability stability_classes{j}; % 重新计算有效源高风速影响抬升 H_eff_current calculate_effective_height(source, meteo_current); % 计算整个网格的浓度 C_grid calculate_concentration_grid(receptor, source, meteo_current, H_eff_current); % 记录最大地面浓度及其位置 max_concentrations(i, j) max(C_grid(:)); [idx] find(C_grid max_concentrations(i, j)); % ... 可以存储位置信息 end end % 将结果可视化例如绘制热图 figure; imagesc(max_concentrations); colorbar; set(gca, XTick, 1:length(stability_classes), XTickLabel, stability_classes); set(gca, YTick, 1:length(wind_speeds), YTickLabel, wind_speeds); xlabel(大气稳定度); ylabel(风速 (m/s)); title(不同情景下最大地面浓度 (kg/m^3));5.3 结果可视化与专业出图除了基本的等值线图还可以创建更丰富的可视化来展示结果。浓度剖面图绘制沿中心线y0浓度随下风向距离变化的曲线直观显示浓度衰减过程。三维表面图展示整个浓度场的三维形态。动画如果模拟了不同时间步如逐时变化可以制作浓度场演变动画非常直观。% 绘制中心线浓度剖面 x_line linspace(receptor.x_range(1), receptor.x_range(2), 200); y_line zeros(size(x_line)); C_line zeros(size(x_line)); for k 1:length(x_line) C_line(k) gaussian_plume(x_line(k), y_line(k), receptor.z, ... source.Q, meteo.u, source.H_eff, meteo.stability); end figure; plot(x_line, C_line*1e6, b-, LineWidth, 2); % 浓度转换为μg/m³ grid on; xlabel(下风向距离 (m)); ylabel(SO_2 浓度 (μg/m^3)); title(烟羽中心线地面浓度分布); % 标记最大浓度点 [maxC, idx] max(C_line); hold on; plot(x_line(idx), maxC*1e6, ro, MarkerSize, 10, MarkerFaceColor, r); text(x_line(idx), maxC*1e6*1.05, sprintf(Max: %.1f μg/m³ %.0f m, maxC*1e6, x_line(idx)), ... HorizontalAlignment, center);6. 常见问题、模型校验与避坑指南在实际操作中你会遇到各种预料之外的问题。下面是我从多次项目中总结的一些经验。6.1 模型校验你的模拟结果可信吗一个模型如果无法验证其价值就大打折扣。校验通常有以下几种方法与解析解或经典案例对比对于标准高斯模型可以手动计算几个特定点的浓度与程序输出对比。或者找一篇使用了相同模型和参数的学术论文对比其结果。与监测数据对比这是最理想但也最困难的方法。需要获取模拟时段内、下风向监测点的实际浓度数据。由于背景浓度和其他源的存在直接对比往往差异很大。通常的做法是在背景浓度较低的时段如夜间模拟单个主导源的影响并与监测数据的增量进行对比。相关系数和标准化平均偏差是常用的评价指标。敏感性分析检查模型输出对关键输入参数如风速u、有效源高H、扩散参数方案变化的敏感程度。如果浓度结果对某个参数极其敏感而该参数又非常不确定那么模型预测的不确定性就很大。6.2 典型问题排查表问题现象可能原因排查思路与解决方法计算浓度全部为0或NaN1. 数组运算维度不匹配。2. 除法中分母如u,σy,σz为0。3. 指数运算溢出如距离为负。1. 使用size()函数检查所有参与运算的变量维度确保能进行逐元素运算使用.运算符。2. 在计算前加入判断确保u0x0。对于σy和σz检查其计算公式确保在x很小时不会为0可设置一个最小值如1e-5。3. 检查输入的距离、坐标是否为负。浓度值异常高如远超源强1. 单位不一致。这是最常见错误2. 风速u输入值过小如误用 km/h 代替 m/s。3. 有效源高H计算错误可能为负值或极小值。1.统一使用国际单位制SI长度-m时间-s质量-kg。将Q从 g/s 转为 kg/s浓度结果单位是 kg/m³再根据需要转为 mg/m³ 或 μg/m³。在代码开头用注释明确所有变量的单位。2. 检查风速单位换算1 m/s 3.6 km/h。3. 打印中间变量H的值检查抬升高度计算函数逻辑。浓度分布图不对称或形状怪异1. 受体网格定义有误x和y向量顺序不对。2. 在计算σy和σz时对网格X直接使用了矩阵而公式期望输入是标量或向量。3. 地面反射项计算错误。1. 使用meshgrid或ndgrid时明确X和Y的维度。绘制网格点scatter(X(:), Y(:))检查网格是否正确。2. 确保calc_sigma函数能处理向量输入使用逐元素运算符.*和.^。3. 检查地面反射项exp(-(zH)²/(2σz²))是否已正确加入。最大浓度出现位置与理论不符理论最大地面浓度通常出现在x ≈ H/√2附近对于地面源。对于高架源位置更远。1. 检查扩散参数公式是否适用于当前的稳定度类别和距离范围。不同公式在不同距离上差异很大。2. 检查有效源高H是否计算准确它直接决定了最大浓度点的距离。3. 受体网格分辨率可能不够无法捕捉到准确的峰值点。尝试加密网格。6.3 实操中的关键技巧与心得从简单开始逐步复杂化不要一开始就试图模拟最复杂的场景。先用一组标准参数在平坦地形、稳态条件下把基础模型调通画出合理的浓度分布图。然后再依次加入地形修正、风速变化、多源叠加等模块。每加一个功能都要验证其结果是否物理合理。做好量纲单位管理这是建模中最容易出错的地方。我的习惯是在脚本的最开头将所有输入参数从原始单位如 g/s, km/h统一转换为 SI 单位kg, m, s。在计算结果的最后再根据输出需求转换如 kg/m³ 转 μg/m³。在变量名或注释中注明单位。可视化是调试的最佳工具当结果不对劲时别光看数字。把中间变量都画出来看看σy和σz随距离变化的曲线对吗有效源高H随风速变化的趋势合理吗一张图往往比一堆数字更能揭示问题。理解模型的局限性高斯模型适用于尺度在几十公里以内、地形相对平坦、污染物化学惰性的情况。对于城市街谷、复杂山区、光化学烟雾等问题它的误差会很大。知道模型在哪里会失效和知道它在哪里适用同样重要。代码模块化将计算有效源高、计算扩散参数、计算浓度的函数分别独立编写。这样不仅代码清晰易于调试也方便你未来替换不同的子模型比如换一种抬升高度公式或换一套扩散参数方案。最后我想分享的一点体会是数学建模的魅力在于它用简洁的数学语言刻画了复杂的现实世界。空气质量模拟的模型从高斯点源开始可以扩展到面源、线源可以耦合化学、考虑地形甚至可以与气象模型在线耦合。这个从简到繁的过程正是我们认识问题、解决问题的典型路径。当你用自己写的代码模拟出污染物扩散的“烟羽”图并与理论规律吻合时那种成就感是无可替代的。希望这个从理论到实战的完整拆解能为你打开这扇门。