简介这份数学建模竞赛资源围绕核电站泄漏后放射性气体浓度分布规律与扩散模型展开适合参加全国大学生数学建模竞赛的本科生、研究生以及从事环境风险评估或应急管理研究的人员。资料以陕西师范大学2011年模拟赛题为载体完整呈现了从问题重述、模型假设、符号说明到模型建立与求解的整套文档并给出了高斯烟羽模型、一维与三维抛物型扩散模型、有限时间泄漏扩散模型及高斯烟团模型的推导过程与MATLAB求解思路。压缩包内共1个doc文件大小约1.18MB内容结构清晰包含公式推导、符号说明和福岛核泄漏案例应用可作为同类赛题建模、论文写作与模型对比分析的重要参考。资源已有204人学习使用适合需要借鉴完整建模流程和扩散模型构建方法的学习者下载。1. 标题指向的到底是什么从竞赛题到可落地的气体扩散建模任务“43367核电站泄漏后放射性气体浓度分布规律和气体扩散模型研究”这类题目是数学建模竞赛里“大气环境应急”方向的代表2025年研究生数学建模竞赛E题也走这条路。它要我们做一件事给定核电站事故释放率、气象观测和监测点数据反推或预测下风向放射性气体浓度分布并给出一个可解释的扩散模型。不是让谁真的去处理核事故而是考察你从物理规律到数值实现再到可信度验证的完整能力。最划算的路线是用高斯烟羽/烟团解析模型而不是大动干戈做CFD。这篇文章就按我自己在类似赛题里跑通的流程把浓度分布规律、参数选型、踩坑点和验证脚本一次说透适合第一次选这类题的参赛队。2. 放射性气体浓度分布规律为什么高斯烟羽模型是竞赛默认起点2.1 浓度分布的核心物理量源项、气象场与地面浓度放射性气体进入大气后浓度分布是平流和湍流扩散共同作用的结果。平流是风把污染物往下风向推扩散是湍流把它往四周稀释。物理量上我们需要三个输入和一个输出源项释放率QBq/s、气象场风速u、风向、大气稳定度和受体坐标下风向距离x、横向距离y、离地高度z输出就是该点的放射性核素活度浓度CBq/m³。竞赛里给的原始数据往往不是这个口径常见的是“事故开始后某小时内释放总量”“10m高度逐时风速风向”“混合层高度”需要先转换成模型参数。浓度分布规律有个直观的几何特征烟羽中心线浓度随距离衰减横风向剖面呈钟形曲线。理论上湍流的随机游走使粒子在横向和垂直方向上的位置分布近似正态所以高斯烟羽模型的解析解长这样C(x,y,z) Q / (2π·u·σy·σz) · exp(-y²/(2σy²)) · [exp(-(z-H)²/(2σz²)) exp(-(zH)²/(2σz²))]式子里x是下风向距离y是横风向距离z是离地高度σy和σz是水平和垂直扩散参数随x增长H是有效源高等于烟囱高度加上烟气抬升高度。式子里的第二个指数项是地面反射项地面不能穿透污染物会在边界累积等效于源下方有个镜像源这在近地面浓度计算里非常重要。为什么竞赛默认用高斯烟羽而不是纳维-斯托克斯方程原因很实际解析解只缺几个参数十分钟能跑完全场评委关注的是参数是否合理、结论是否有依据CFD的网格和湍流模型反而容易变成黑匣子花了三天还调不出一个收敛解。我一般先写高斯烟羽作为基座如果题目气象条件有明显瞬变再升级成烟团模型这是数模竞赛里成本最低的稳妥路线。但高斯模型也有明显假设地形平坦、气象条件定常、释放连续。核电站事故释放往往持续几小时甚至几天风速风向会变释放率也可能随时间变化所以我们要在2.2节里分清楚“什么情况用烟羽什么情况必须换烟团”。2.2 从高斯烟羽到拉格朗日烟团两种模型的适用边界拉格朗日烟团模型把连续释放离散成一系列小烟团每个烟团在释放时刻的位置就是源点之后随风做平移同时扩散变胖。某时刻某点的浓度是所有烟团贡献的总和。这样做的好处是天然支持时间变化风向变了就改变每个烟团的平移向量释放停了就不再发新烟团已发的继续向下游走并稀释。缺点是计算量比烟羽大但在五分钟模拟时长内完全可接受。竞赛里如何选型我常用一张判断表场景特征推荐模型理由连续泄漏、风场稳定、要稳态峰值浓度高斯烟羽公式简单、参数少、容易做敏感性分析短时泄漏30分钟内或释放中途停止高斯烟团能刻画浓度的时间上升/下降段风向在几小时内摆动超过30°高斯烟团烟羽模型会错误地把浓度算在一个最长方向上平坦开阔下风向数百公里高斯烟羽干湿沉积衰减项烟团粒子太多时内存吃紧建筑物尾流、复杂山地峡谷CFD或风洞竞赛不推荐解析模型在这里基本失效这张表是选型的基准。我在处理研究生数学建模竞赛E题这类题目时会先看气象数据如果逐时风向的标准差超过15°就直接放弃稳态烟羽改用时间步进烟团。切换成本不高核心就是3.3节要写的叠加逻辑。有效源高H对浓度分布的影响也值得单独说。H取得偏大地面浓度峰值就离源更远且更小H偏小峰值在近源处提前出现。很多队伍忽略烟气抬升高度直接拿烟囱几何高度计算。常见做法是用霍兰德公式做中性条件抬升估计或者题目如果给了烟囱出口温度和速度用布里格斯公式。但如果题目没给这些信息就说明出题人不想在这里卡你取烟囱高度即可然后在论文里写明这是保守估计。扩散参数σy和σz决定了烟羽的铺展速度。同样在5.1节里给查表法这里说结论σ主要由大气稳定度决定。稳定度A类极不稳定时烟羽快速变宽地面浓度低F类稳定时烟羽窄而高可能造成下风向较远距离的带状高浓度。这个机制直接解释了为什么夜间稳定天气下核电泄漏的地面监测反而可能更危险——污染物不容易垂直混合保持高浓度被风吹到更远。混合层高度也是浓度分布的一个隐形边界。当混合层高度较低时垂直扩散受限烟羽在垂直方向会反复反射形成近似均匀混合。如果H接近混合层高度地面浓度会比开阔地形模型明显增加。处理办法是加多层镜像源但如果混合层高度低于有效源高公式要改成受限扩散形式否则计算出的浓度会违反质量守恒。3. 建扩散模型的完整步骤从读题到可提交的求解代码3.1 数据整理与参数归一化竞赛题给的数据通常很杂逐时气象表、厂区坐标、释放总量、少数监测站浓度。第一步是把它们统一成模型输入。我一般先用pandas读入CSV只保留时间、风速、风向、稳定度字段。然后做三件事一是把风向从角度转成直角坐标中的风矢量二是把风速从10m高度修正到有效源高三是把稳定度映射到扩散参数。风向转换是高频出错点。气象上定义的是“来风方向”比如270°表示西风风从西边来吹向东。模型需要的是烟羽平移方向也就是来风的反方向。如果直角坐标x轴向东、y轴向北来风角度θ从正北顺时针计那么烟羽速度分量是u_x -wind_speed · sin(θ) u_y -wind_speed · cos(θ)注意这里不要用气象学里那种“从x轴正方向逆时针转”的角度两个坐标系混用会让烟羽向反方向跑。我在代码里专门写一个转换函数把角度处理逻辑固定下来。import pandas as pd import numpy as np def wind_angle_to_vector(wind_speed, wind_angle_deg): 把气象风速风向转换为烟羽平移速度向量。 wind_angle_deg: 来风方向角正北为0顺时针增加。 返回 (ux, uy)分别是向东、向北的速度分量(m/s)。 theta np.deg2rad(wind_angle_deg) ux -wind_speed * np.sin(theta) uy -wind_speed * np.cos(theta) return ux, uy # 示例气象报文给出 风速5m/s风向270度(西风) ux, uy wind_angle_to_vector(5.0, 270) print(fux{ux:.2f}, uy{uy:.2f}) # 预期烟羽向东运动ux5, uy≈0逻辑说明把风速和方位角换算成直角分量时θ为0度时烟羽向南uy-v对应北风θ为90度时烟羽向西ux-v对应东风正好符合“来风方向的反方向就是去向”。参数说明wind_speed单位是m/swind_angle_deg是气象角度不是数学极坐标角度如果你看到数据里风向用“东南风”这种文字先查表换成角度再调用。接下来是参数归一化。释放率如果是总量而不是速率要除以释放持续时间距离全部转成米时间步长统一成秒。这样后面网格和扩散参数才不会出现“厘米配千米”的荒唐量级。def normalize_release(total_bq, duration_s, given_rateNone): 把总释放量或给定速率统一成平均释放率 Bq/s。 if given_rate is not None: return given_rate if duration_s 0: raise ValueError(释放时长必须为正) return total_bq / duration_s # 题目给出事故总释放 3.7e16 Bq释放持续2小时 Q normalize_release(3.7e16, 2 * 3600) print(f平均释放率 Q {Q:.3e} Bq/s)这段的逻辑是高斯烟羽模型要求输入源强是单位时间排放活度。竞赛题有时写“排放总量”有时写“释放速率”必须先统一成后者。参数说明duration_s如果题目说的是“开始后第1小时内释放了70%”不要机械除以一小时先按题目给的释放曲线折算。3.2 用Python实现稳态高斯烟羽浓度场核心函数就是2.1节的公式。我用numpy直接对网格计算这样后面画等值线非常快。注意要把稳定度导出的σy、σz作为可调函数传入方便做敏感性分析。def gaussian_plume_grid(X, Y, z, Q, u, H, sigma_y_func, sigma_z_func): 稳态高斯烟羽地面/任意高度浓度场。 X, Y网格坐标单位m风向为x方向。 z受体离地高度地面取1.5m呼吸带高度。 Q源强 Bq/su有效风速 m/sH有效源高 m。 sigma_y_func, sigma_z_func以x为自变量的扩散参数函数。 x np.maximum(X, 1e-6) # 避免源点除零 y Y sy sigma_y_func(x) sz sigma_z_func(x) horizontal np.exp(-(y ** 2) / (2 * sy ** 2)) vertical (np.exp(-(z - H) ** 2 / (2 * sz ** 2)) np.exp(-(z H) ** 2 / (2 * sz ** 2))) C Q / (2 * np.pi * u * sy * sz) * horizontal * vertical return C # 网格下风向0~5km横向±1km分辨率20m x_axis np.linspace(0, 5000, 251) y_axis np.linspace(-1000, 1000, 101) X, Y np.meshgrid(x_axis, y_axis) # 扩散参数先用简单幂律占位5.1节再换成稳定度查表 sigma_y lambda x: 0.22 * x * (1 0.0001 * x) ** (-0.5) # A类 sigma_z lambda x: 0.20 * x C gaussian_plume_grid(X, Y, z1.5, Q1e12, u5.0, H80, sigma_y_funcsigma_y, sigma_z_funcsigma_z)逻辑说明代码严格按照高斯公式实现vertical项中的两个指数分别代表直接烟羽和地面镜像反射。参数说明z1.5是人体呼吸高度不是地面0m因为在辐射剂量计算里关注的是吸入风险H80是烟囱高度加抬升高度如果题目给的是10m风速要先换算到80m高度处的风速否则浓度峰值位置会偏。这里有一个工程细节如果风向不是x方向需要先把整个网格旋转到风向坐标系算完浓度后再旋转回地理坐标。常见做法是构造旋转矩阵def rotate_grid(X, Y, wind_dir_deg): 把地理网格旋转到以平均下风向为x轴的新坐标系。 theta np.deg2rad(wind_dir_deg) Xr X * np.cos(theta) Y * np.sin(theta) Yr -X * np.sin(theta) Y * np.cos(theta) return Xr, Yr调用前先把X、Y都平移成以源为原点再用上面函数旋转。这个步骤漏掉就会发生“风向明明是东北风浓度却出现在源的西南”的翻车现场。3.3 时间维上的烟团叠加与浓度分布动画如果题目要求的是“泄漏后第2小时的浓度分布”或“释放终止后的衰减”稳态烟羽就不够用了。我用拉格朗日烟团叠加把释放时间切成长度相等的Δt每段释放的活度Q·Δt作为一个独立烟团从源出发中心以当时的风速向量平移扩散参数随时间增长。def puff_concentration_2d(times, Q, duration_s, ux, uy, X, Y): 时间步进烟团叠加返回归一化浓度场。 times: 输出时刻数组(s)Q: 释放率Bq/sduration_s: 总释放时长。 ux, uy: 平均风分量(m/s)这里先按常值处理。 dt times[1] - times[0] C_total np.zeros_like(X) sigma0 1.0 # 烟团初始半径 for t_emit in np.arange(0, duration_s, dt): age times - t_emit # 每个输出时刻下该烟团的年龄 age np.maximum(age, 0) # 烟团中心位置 xc ux * age yc uy * age # 扩散参数随时间增长简化幂律 s sigma0 0.3 * age ** 0.8 # 对每个输出时刻累加该烟团贡献 for i, a in enumerate(age): if a 0: continue dis2 (X - xc[i]) ** 2 (Y - yc[i]) ** 2 C_total (Q * dt) / (2 * np.pi * s[i] ** 2) * np.exp(-dis2 / (2 * s[i] ** 2)) return C_total times np.arange(0, 4 * 3600, 300) # 4小时步长5分钟 C_puff puff_concentration_2d(times, Q1e12, duration_s2 * 3600, ux5.0, uy0.2, XX, YY)逻辑说明外层循环按释放时间发射烟团内层循环把每个烟团在所有输出时刻的浓度贡献累加。烟团的扩散参数随年龄增长年龄只取非负值。参数说明dt越小时间分辨率越高但计算量线性增加这里取5分钟对竞赛题足够。初始半径σ0用于避免源点数值奇异一般取1~5m它的影响只在源附近观测点通常都在几百米以外无需过度纠结。这个结果可以配合matplotlib逐帧保存成图片用matplotlib.animation.FuncAnimation做成浓度分布动画。动画在论文答辩里非常加分因为评委能直观看到烟羽随时间的摆头——这是稳态烟羽模型做不到的效果。4. 数学建模竞赛特有的坑让浓度误差高几个量级的4个高频失误4.1 现象稳定度分类用错浓度相差几个量级很多时候队伍用同一个高斯公式算出来的地面浓度和官方参考值差了100倍。原因基本都在稳定度分类上。题目给一句“多云地面弱风”有人就随手选了D类结果σz偏大浓度峰被过度抹平而真实条件可能是F类稳定烟羽窄而集中。解决方法是严格查帕斯奎尔稳定度表先根据太阳高度角、云量、10m风速定级再取对应扩散参数。竞赛里常见的场景是秋季夜间、风速2m/s、云量小于4成那基本就是E或F类。如果实在拿不准就同时算D类和F类两组结果做敏感性区间不要只给一个“拍脑袋”的稳定度。4.2 现象边界反射处理不当导致源附近浓度虚高在公式里只有一次地面镜像反射时如果有效源高70m、混合层高度100m污染物在上下边界之间来回反射。只取一次反射项会漏掉大量的再反射贡献结果源附近浓度比源强每周期的总量还大出现“浓度比源处还高”的荒谬结果。原因是边界条件没有封闭。解决办法是叠加无穷镜像源序列在每个混合层高度上下两侧各设置多个镜像源直到高阶贡献小于主项的1%为止。另一个更工程化的做法是当混合层高度L远大于σz时保留常规反射项当H接近L时改用混合层内完全反射的均匀化公式并检查整个计算域的总活度守恒。4.3 现象风场随高度变化没考虑地面浓度全偏气象站给的是10m高度风速而泄漏源在60m烟囱口。高处的风速通常比地面大很多如果用10m风速代入模型烟羽向下风向推进得慢浓度峰值会出现得更近、更高。原因是忽略风速廓线。常见做法是用幂律修正u(z) u_10 · (z / 10)^pp值随稳定度变化A类约0.07B类0.12C类0.2D类0.28E类0.35F类0.45。我在代码里实现很简单乘以一个系数即可。这个坑最隐蔽因为不报错曲线形状也像模像样只是峰值位置明显偏移评委一眼就能看出来。4.4 现象模型验证只看趋势不看量级有的队伍在论文里画一条浓度曲线和监测点折线叠在一起只写了“趋势一致”就宣称模型有效。这个很危险因为对数坐标下差两个量级也能看起来“趋势一致”。解决方法是定量化把监测点的观测浓度和模拟浓度做散点图横纵轴都取对数计算相对误差的平均值和中位数再算空间相关系数R。R大于0.8且平均相对误差在50%以内才敢写“模型可信”。最后一章我会给现成脚本。5. 参数怎么设才让人信服扩散参数、源强与气象数据的口径5.1 帕斯奎尔稳定度分级与扩散参数的查表法扩散参数σy和σz是浓度分布的核心。竞赛中最推荐的查表法是布里格斯Briggs参数适用于开阔乡村下垫面与核电站厂区地形量级匹配。常见形式是σy和σz取x的幂函数稳定度σy (m)σz (m)A0.22x(10.0001x)^-0.50.20xB0.16x(10.0001x)^-0.50.12xC0.11x(10.0001x)^-0.50.08x(10.0002x)^-0.5D0.08x(10.0001x)^-0.50.06x(10.0015x)^-0.5E0.06x(10.0001x)^-0.50.03x(10.0003x)^-1F0.04x(10.0001x)^-0.50.016x(10.0003x)^-1注意表格里x的单位是米。用这个表写出的σ函数可以直接替换第3章的lambda函数。如果题目明确是城市地形或山区参数要换成城市型Briggs城市系列或按《环境影响评价技术导则 大气环境》附录附表取值。竞赛中不指定下垫面时默认乡村开阔地形是合理的但论文里要写清楚这个假设。查表法的关键是先确定帕斯奎尔稳定度级别。我按这个顺序白天根据太阳辐射等级强、中、弱和风速查表夜晚根据云量多云/少云和风速查表。如果题目没有给云量只给了“夜间”和“阴天”可以直接按D类中性作为基准因为阴天夜间辐射状况接近中性。稳定度级别写进敏感性分析比费半天劲猜一个“真实值”更诚实。5.2 源项估算与释放时长对浓度分布的影响源项是最大的不确定性来源。题目经常不给具体释放率只给“堆芯损坏程度”“安全壳完好性”。此时常见做法是参考事故分级估算放射性惰性气体释放份额或者更直接地把源项当作待反演参数先用多个假设源强计算浓度再和监测站数据做拟合反推出“等效源强”。这个反演思路比试图读一份现实中根本拿不到的机组源项清单更靠谱。释放时长同样影响分布量级。连续释放时长τ如果远大于迁移时间x/u可以按稳态烟羽处理峰值浓度与τ无关如果τ很短烟羽前段和后段的浓度分布明显不对称。判断方法很简单计算迁移时间t_travel x/u和释放时长τ比较。τ 10·t_travel时稳态误差小于10%放心用烟羽否则用3.3节的烟团叠加。我在竞赛里会把τ和t_travel的比值写进模型选择一节评委看到这种逻辑就很踏实。源高以上说的都是几何高度修正。如果题目给了烟气温度80°C、烟囱出口速度12m/s可以用布里格斯抬升公式算ΔH然后把有效源高H几何高度ΔH代入模型。没给就写“不考虑浮力抬升保守估计”可以有效避免为了不用公式而显得知识盲区。5.3 敏感性与不确定性分析的简洁做法竞赛时间紧张不需要蒙特卡洛。我用局部扰动法以基准参数算一组浓度场然后把每个关键参数上下扰动合理范围看目标点浓度的变化倍数。参数范围取工程惯例源强±30%风速±20%稳定度上下各移一类有效源高±20%。参数基准值扰动范围峰值浓度变化释放率Q1e12 Bq/s±30%±30%线性风速u5 m/s4~6约-17%~25%稳定度D类C~E可能变化2~4倍有效源高H80m60~100峰值距离偏移约20%这张表的价值在于告诉评委你的结论对哪个参数敏感。通常稳定度最敏感这说明你尽到责任了。敏感性分析不只用于写论文也用于反演源强时的迭代步长选择。6. 一个能抄的验证技巧用浓度等值线和误差指标反推模型参数6.1 把计算结果画成浓度等值线并加观测点对照模型算完第一件事就是画图。等值线比彩色云图更适合竞赛论文因为能直接读浓度数值。我会用matplotlib的contourf画填充图再用contour叠加等值线标出源位置和监测站位。方向一定要和厂区地图一致否则监测点全错位。import matplotlib.pyplot as plt def plot_concentration_map(X, Y, C, monitor_points, save_path): 绘制浓度等值线并叠加监测点。 monitor_points: [(x_km, y_km, label), ...] fig, ax plt.subplots(figsize(8, 6)) # 浓度取对数分级因为放射性浓度跨越多个量级 levels np.logspace(np.log10(C[C 0].min()) if C.max() 0 else 0, np.log10(max(C.max(), 1e-6)), 10) cf ax.contourf(X / 1000, Y / 1000, C, levelslevels, cmapReds, extendboth, normplt.LogNorm()) cs ax.contour(X / 1000, Y / 1000, C, levelslevels, colorsblack, linewidths0.5) ax.clabel(cs, fmt%.1e) ax.plot(0, 0, k^, markersize10, label源点) for x_km, y_km, label in monitor_points: ax.plot(x_km, y_km, bo, markersize6) ax.annotate(label, (x_km, y_km), textcoordsoffset points, xytext(5, 5), fontsize9) ax.set_xlabel(下风向距离 (km)) ax.set_ylabel(横向距离 (km)) ax.set_title(放射性气体浓度分布 (Bq/m³)) plt.colorbar(cf, labelBq/m³) plt.savefig(save_path, dpi300, bbox_inchestight)逻辑说明浓度跨量级分布等值线按对数间隔画图上每个点都能读出数值。监测点用蓝色圆点叠上去可以直观看到“模拟是否落在观测值附近”。参数说明LogNorm需要C所有值非负如果有零值先用np.maximum(C, 1e-6)压低levels的个数根据网格范围调整10层对5km尺度通常够用。6.2 用相对误差与空间相关系数判断模型可信度画完图还要定量验证。我一般把监测站处的实测浓度和模拟浓度取出来做三个指标中位相对误差、空间相关系数、以及观测/模拟散点图和1:1线。这三个指标足够撑起论文里的“模型验证”小节。def evaluate_model(C_sim, points_xy, C_obs): 在监测点处评估模拟浓度与观测浓度的误差。 points_xy: [(x, y), ...]单位m与C_sim同坐标系。 C_obs: 对应监测点的实测浓度数组Bq/m³。 C_sim_pts [] for (x, y) in points_xy: ix int((x - x_axis[0]) / (x_axis[1] - x_axis[0])) iy int((y - y_axis[0]) / (y_axis[1] - y_axis[0])) C_sim_pts.append(C_sim[iy, ix]) C_sim_pts np.array(C_sim_pts) C_obs np.array(C_obs) # 相对误差非负便于取中位数 rel_err np.abs(C_sim_pts - C_obs) / np.maximum(C_obs, 1e-12) med_rel_err np.median(rel_err) # 空间相关系数对浓度取对数后计算避免个别高值主导 corr np.corrcoef(np.log10(C_sim_pts 1e-12), np.log10(C_obs 1e-12))[0, 1] return med_rel_err, corr med_err, corr evaluate_model(C, [(1000, 0), (2000, 200), (3000, -150)], [1e5, 2e4, 8e3]) print(f中位相对误差{med_err:.2%}, 对数相关系数R{corr:.3f})逻辑说明如果模拟浓度和观测浓度分布规律一致散点会聚在1:1线附近相关系数接近1。中位相对误差比均值更抗异常点一个远离烟羽的监测站数值很小会把均值拉爆中位数更稳定。参数说明取对数后再算相关是因为放射性浓度从源到几公里外跨越5~6个量级直接算线性相关会被最大值主导对数让每个量级都有同等的发言权。这是我踩过坑之后改成的标准做法。我在最后提醒一个反向使用这类验证的套路如果监测点足够多可以把“假设源强”当作可调参数反复跑2~3组源强观察哪个假设让相关系数最高、相对误差最小这组源强就是反演源强。这在竞赛题里往往比题目给的名义释放量更接近真实值也会成为论文里“基于监测数据的源项修正”亮点的来源。模型验证不只是给评委看的它也是帮你找到参数错误的第一线索。上次比赛我就是因为相关系数只有0.3回去查了代码发现风向旋转少乘了个负号修好之后相关系数跳到0.87最后拿奖的就是那个版本。把这个验证流程写进论文才是真正让模型可信的办法希望帮到你。本文还有配套的精品资源点击获取