简介：面向气动仿真与智能优化算法学习者，这份代码包提供了基于改进遗传算法与粒子群算法的高斯烟羽模型气体扩散模拟完整方案。主程序 main.m 统一调度两种智能寻优策略，配合适应度函数、泄漏速率计算、空间点浓度求解等辅助模块，可在 Matlab 2019b 中直接运行，适合研究大气污染扩散、危险气体泄漏评估以及智能算法对比实验。压缩包共 16 个文件，含 8 个 .m 脚本、7 张可复现的运行结果图与 1 份数据表格，代码体量精炼，结果图能直观呈现浓度分布与算法收敛过程，整体仅 208KB，轻量易部署。已有 2380 人学习下载。通过该资源，读者可掌握改进 GA/PSO 在高斯烟羽模型参数优化中的具体实现，理解遗传与粒子群算法的收敛差异，并借助核心扩散函数快速迁移到弹道、环境工程等同类仿真场景。
1. 气体扩散反演为什么需要改进的GA和PSO
高斯烟羽模型是泄漏气体扩散模拟中最常用的解析模型，给定泄漏源位置、泄漏强度和气象条件，就能算出下风向任意监测点的浓度。但实际应急场景里问题往往是反过来的：只有一堆传感器浓度读数，要反推泄漏源在哪、漏了多少。这类反演问题没有解析解，目标函数多峰、非线性强，普通最小二乘容易陷入局部极值。遗传算法和粒子群算法都是全局搜索方法，但标准版本一个收敛慢、一个容易早熟，直接用来做气体扩散反演，经常出现源位置偏差几十米、泄漏率偏差一个数量级的结果。改进思路集中在自适应交叉变异、惯性权重动态调整和精英保留策略上。这套基于Matlab的实现把高斯烟羽正算模块和两类改进群智能算法封装在一起，适合做环境监测、事故溯源和应急预案仿真的工程人员直接改参数使用。
2. 高斯烟羽模型建模与目标函数设计
2.1 连续泄漏的高斯烟羽浓度场计算公式
高斯烟羽模型假定泄漏源连续稳定排放，湍流扩散在水平和垂直方向都服从正态分布。地面全反射条件下，下风向任意坐标点的浓度计算式为：
function C = gaosiyanyu(x, y, z, Q, u, H, sigmay, sigmaz) % x,y,z: 预测点坐标，原点在泄漏源地面投影点，x为下风向 % Q: 泄漏源强，单位kg/s % u: 平均风速，单位m/s % H: 有效源高，单位m % sigmay, sigmaz: 水平和垂直扩散参数，单位m term1 = Q / (2 * pi * u * sigmay * sigmaz); term2 = exp(-y^2 / (2 * sigmay^2)); term3 = exp(-(z - H)^2 / (2 * sigmaz^2)) + exp(-(z + H)^2 / (2 * sigmaz^2)); C = term1 * term2 * term3; end这里最关键的是有效源高H。实际工程中气体可能从烟囱排出，烟气抬升高度不能忽略；如果模拟地面储罐泄漏，H就是泄漏口高度加少量抬升。很多使用者直接把H设成泄漏源高度，这会导致近地面浓度计算偏高。代码里给的xielousulv.m就是计算泄漏速度的辅助函数，它把液池蒸发或管道泄漏的质量流量换算成Q，再传给浓度计算模块。扩散参数sigmay和sigmaz通常按 Pasquill 稳定度等级查表得到，稳定度分为 A 到 F 六类，对应强对流到强逆温，扩散能力逐级减弱。城市和下垫面粗糙度不同，查表值要做粗糙度修正，否则下风向浓度分布会整体偏移。
2.2 传感器布置与正演模拟
反演之前必须先有正演样本，也就是给定一组泄漏源参数，计算出每个传感器位置的浓度理论值。传感器一般布置在泄漏源下风向扇形区域，因为烟羽中心线浓度最高，侧向浓度随距离快速衰减。我在仿真时通常按 10 到 30 个监测点布置，坐标从泄漏源预估位置的下风向 50 米到 500 米范围内展开，横向在中心线两侧按高斯分布取点。这样布置的好处是传感器读数之间的区分度大，反演算法更容易分辨不同源参数组合的差异。
正演模拟流程是：读取传感器坐标，调用gaosiyanyu函数逐个计算浓度，叠加 5% 的高斯噪声模拟真实测量误差。噪声幅值不能设太大，否则反演结果会偏离真值；也不能完全不设，否则算法会过拟合到数值精度上。一般做法是用随机数对理论浓度做乘性扰动，C_measured = C_true * (1 + 0.05 * randn(size(C_true)))。这里randn产生标准正态分布随机数，乘 0.05 就是标准差 5% 的相对误差，这比加性噪声更接近传感器响应的物理特性，因为浓度跨越几个数量级，加性噪声在低浓度区域会直接淹没信号。
2.3 适应度函数：残差和归一化设计
遗传算法和粒子群算法在气体扩散反演里都靠适应度函数引导搜索方向，适应度函数的设计直接决定能不能收敛到真值附近。目标是最小化传感器实测浓度与模型预测浓度的差异，但不能直接用绝对误差做目标函数。原因是下风向 500 米处浓度可能只有 10 mg/m³，而上风向近距离传感器浓度可能是几千 mg/m³，绝对误差会被大浓度点主导，小浓度点的信息完全丢失。
function f = fitness1(params, sensor_x, sensor_y, sensor_z, C_meas, u, H, stability) % params = [Q, x0, y0]，分别是被优化参数：源强、源x坐标、源y坐标 C_pred = zeros(size(C_meas)); for i = 1:length(C_meas) % 预测点相对源位置做坐标平移，然后把风速方向对齐到x轴 dx = (sensor_x(i) - params(2)) * cos(wind_dir) + ... (sensor_y(i) - params(3)) * sin(wind_dir); dy = -(sensor_x(i) - params(2)) * sin(wind_dir) + ... (sensor_y(i) - params(3)) * cos(wind_dir); [sigmay, sigmaz] = compute_sigma(stability, dx); C_pred(i) = gaosiyanyu(dx, dy, sensor_z(i), params(1), u, H, sigmay, sigmaz); end % 对数归一化残差：降低大浓度点权重，保留小浓度信息 residual = abs(log(C_meas + 1e-10) - log(C_pred + 1e-10)); f = sum(residual); end这段代码里做了两层关键处理。第一层是风向旋转，把传感器坐标从全局坐标系转换到以风速方向为 x 轴的烟羽坐标系，否则gaosiyanyu函数里下风向的x含义就错了。第二层是对数归一化，log把浓度从指数尺度拉回线性尺度，让 10 mg/m³ 和 1000 mg/m³ 的误差在目标函数中有相近的贡献。1e-10是防止对零取对数，实际使用时我会把它设成传感器检测限的百分之一，这样比固定一个小数更合理。fit.m和fitness1.m的差别在于，前者还可以被遗传算法直接当作排序依据，后者只是返回残差向量，需要外部再取mean或sum才能作为适应度标量。
| 参数 | 物理含义 | 优化范围示例 |
|---|---|---|
Q | 泄漏源强 | 0.01 ~ 5 kg/s |
x0 | 泄漏源 x 坐标 | -100 ~ 300 m |
y0 | 泄漏源 y 坐标 | -100 ~ 100 m |
u | 平均风速 | 固定为测量值，或加入优化范围 0.5~5 m/s |
H | 有效源高 | 固定或随 Q 变化 |
如果风速本身不确定，也可以把u加入优化向量，但此时要注意u与Q存在耦合关系，因为浓度公式里Q和u以除法形式出现，可能出现多组参数组合对应相近浓度分布的情况。我在实际项目中优先固定风速，除非现场有多个测风站并且数据明显分层。
3. 改进遗传算法与粒子群算法在Matlab中的实现
3.1 编码与种群初始化
遗传算法和粒子群在这里都采用实数编码，每个个体是一个三维向量[Q, x0, y0]，不需要二进制编码和解码转换，搜索精度更高。种群初始化要在参数范围内尽量均匀覆盖，用均匀分布的随机数生成初始种群。
% 粒子群和遗传算法共用的初始化逻辑 nVar = 3; % 优化变量数量：源强、源x、源y lb = [0.01, -100, -100]; % 下界 ub = [5, 300, 100]; % 上界 nPop = 60; % 种群大小 X = rand(nPop, nVar) .* (ub - lb) + lb; % 均匀随机生成种群这里rand生成nPop x nVar的伪随机矩阵，乘上(ub - lb)把范围放大到搜索区间宽度，再加lb平移到区间起点。注意粒子群算法里每个粒子还要额外分配速度矩阵，遗传算法则需要保存个体的适应度值。边界处理我一般用两种方式：一是反射，越界的个体把对应维度的值拉回边界内并向内反弹；二是重新初始化，让越界个体重新随机生成。反射方式收敛更快，适合知道参数范围比较可靠的场景；重新初始化能保持种群多样性，适合对源位置完全没有先验信息的情况。项目里的mGA.m和mPSO.m都默认用反射方式，因为气体泄漏反演的搜索空间通常来自前期勘察，范围不会太离谱。
3.2 遗传算法改进：自适应交叉变异与精英保留
标准遗传算法有两个突出问题：交叉概率和变异概率全流程固定，前期容易丢失好模式，后期又无法跳出局部最优。改进方案是让交叉概率和变异概率随种群适应度自适应变化，适应度高的个体用较小的交叉概率和变异概率保护其结构，适应度低的个体用较大的变异概率增强探索。
% mGA.m 中自适应交叉变异的核心片段 for i = 1:maxIter [~, idx] = sort(fitness); % 适应度升序排序 fmin = fitness(idx(1)); favg = mean(fitness); for k = 1:2:nPop-1 % 根据父代适应度计算自适应交叉概率 f1 = fitness(k); f2 = fitness(k+1); fmax_curr = max(fitness); if f1 >= favg pc = 0.6 * ((fmax_curr - f1) / (fmax_curr - fmin + 1e-10)); else pc = 0.9; end if rand < pc % 算术交叉：子代为父代线性组合，组合系数随机 alpha = rand; child1 = alpha * X(k,:) + (1 - alpha) * X(k+1,:); child2 = alpha * X(k+1,:) + (1 - alpha) * X(k,:); X(k,:) = child1; X(k+1,:) = child2; end end % 自适应变异：变异步长随迭代次数衰减，后期做精细搜索 for k = 1:nPop sigma = (maxIter - t) / maxIter * 0.2 * (ub - lb); if rand < 0.05 X(k,:) = X(k,:) + sigma .* randn(1, nVar); X(k,:) = max(X(k,:), lb); X(k,:) = min(X(k,:), ub); end end end自适应交叉的概率设定逻辑是：适应度低于平均值的个体大胆交叉，高于平均值的个体谨慎保留。0.6和0.9这两个基准值是经验值，0.9保证劣质个体有足够机会重组，0.6给优质个体留有余地。变异步长sigma与当前迭代代数成反比，这模拟了退火思想，先期大步长探索大范围，后期小步长局部打磨。代码里randn生成服从标准正态分布的扰动，乘sigma后的扰动幅度与参数边界宽度相关，因此不管源强是 0.1 还是 3，扰动比例都一致。注意边界处理用了两层max/min截断，这比反射简单，但会在边界频繁截断时积累大量边界个体，降低多样性，所以实际使用中我建议把反射逻辑补充进去。
精英保留是另一个重要改进：每代把适应度最好的ceil(0.1 * nPop)个个体直接复制到下一代，不参与交叉和变异。这样可以保证最优解单调不恶化的同时，让其他个体充分探索。如果去掉了精英保留，GA 经常出现上一代已经找到的好解在下一代被交叉破坏，导致收敛曲线反复震荡，这正是很多初学者觉得遗传算法不如粒子群稳定的原因。
3.3 粒子群算法改进：惯性权重衰减与速度限幅
标准粒子群的速度更新公式里，惯性权重w是常量，过大则粒子飞行过快容易飞过最优点，过小则粒子群迅速聚集陷入局部最优。改进算法采用线性递减惯性权重，从 0.9 衰减到 0.4，前期全局搜索，后期局部精化。速度限幅约束粒子每一步的最大移动距离，防止粒子飞出搜索空间后卡在边界上。
% mPSO.m 中惯性权重衰减和速度更新的核心代码 global position velocity pbest gbest % 初始化阶段省略，进入迭代循环 for t = 1:maxIter w = 0.9 - (0.9 - 0.4) * t / maxIter; % 线性递减惯性权重 c1 = 2.0; % 个体学习因子 c2 = 2.0; % 社会学习因子 maxVel = 0.2 * (ub - lb); % 速度限幅，每步移动不超过区间的20% for i = 1:nPop % 更新速度并限幅 velocity(i,:) = w * velocity(i,:) ... + c1 * rand(1,nVar) .* (pbest(i,:) - position(i,:)) ... + c2 * rand(1,nVar) .* (gbest - position(i,:)); velocity(i,:) = max(-maxVel, min(maxVel, velocity(i,:))); % 更新位置并处理边界反射 position(i,:) = position(i,:) + velocity(i,:); for j = 1:nVar if position(i,j) < lb(j) position(i,j) = 2 * lb(j) - position(i,j); % 边界反射 velocity(i,j) = -velocity(i,j); elseif position(i,j) > ub(j) position(i,j) = 2 * ub(j) - position(i,j); velocity(i,j) = -velocity(i,j); end end % 计算新位置适应度，更新个体最优和全局最优，这部分与标准PSO一致 end endc1 = c2 = 2.0是经典设置，粒子同时被自身历史最优和全局最优拉到中间位置。有些改进版会让c1随迭代递减、c2递增，早期多学习自己、后期多跟随全局，实际效果在气体扩散反演里没有明显优势，因为搜索空间只有三维，粒子数量足够时经典参数已经能较好地平衡探索和开发。真正影响结果的是maxVel的取值，我把它设成参数区间宽度的 20%，也就是单次迭代最多只能穿越搜索区间的五分之一。如果设得太大，粒子会在几个周期内反复弹跳；如果设得太小，粒子靠近最优解后无法快速收敛。
% fitness1.m 与 fit.m 的调用关系 % mGA.m 和 mPSO.m 都需要把个体参数解码后传给适应度函数 function val = fit(pop, sensor_data, wind_info) Q = pop(1); x0 = pop(2); y0 = pop(3); val = fitness1([Q, x0, y0], sensor_data.x, sensor_data.y, ... sensor_data.z, sensor_data.C, wind_info.u, ... wind_info.H, wind_info.stability); endfit.m是适应度函数的薄封装，主要做参数解包和调用fitness1。在编写这两个算法时，要注意所有全局变量的声明必须放在函数最前面，否则 Matlab 会把pbest当成局部变量处理，导致粒子群算法完全无法收敛。我在调试时经常遇到这类低级错误，建议使用globals或者把粒子群状态封装成 struct 传递，后者更规范，但代码量会稍大。
4. main.m 联合仿真：数据流与运行结果分析
4.1 主程序框架与文件依赖
整个仿真项目以main.m为唯一入口，程序运行后依次完成参数定义、生成模拟观测数据、分别调用改进 GA 和 PSO、输出反演结果并绘图。项目里各文件的依赖关系是：main.m调用gaosiyanyu.m生成正演浓度，调用fitness1.m和fit.m作为优化目标函数，mGA.m和mPSO.m是算法主体，point.m生成传感器坐标网格，xielousulv.m可以把泄漏速率换算成源强。
主程序的典型结构如下：
% main.m 核心流程 % 第一步：定义气象条件与泄漏源真值 u = 3.0; % 风速 m/s H = 10; % 有效源高 m stability = 'D'; % 中性稳定度 Q_true = 1.2; % 真实源强 kg/s pos_true = [50, 30]; % 真实泄漏源坐标 [x, y] m % 第二步：生成传感器坐标与模拟观测浓度 numSensors = 20; sensor_xy = point(numSensors); % 从 xls 读取或函数生成监测点 C_meas = zeros(numSensors, 1); for i = 1:numSensors dx = sensor_xy(i,1) - pos_true(1); dy = sensor_xy(i,2) - pos_true(2); [sy, sz] = compute_sigma(stability, dx); C_meas(i) = gaosiyanyu(dx, dy, 1.5, Q_true, u, H, sy, sz) * (1 + 0.05*randn()); end % 第三步：粒子群反演 [Q_pso, pos_pso, err_pso] = mPSO(C_meas, sensor_xy, u, H, stability); % 第四步：遗传算法反演 [Q_ga, pos_ga, err_ga] = mGA(C_meas, sensor_xy, u, H, stability); % 第五步：对比真值与反演结果并绘制 figure; subplot(2,1,1); plot(err_pso); hold on; plot(err_ga); legend('PSO','GA'); subplot(2,1,2); scatter(pos_true(1), pos_true(2), 'rx'); hold on; scatter(pos_pso(1), pos_pso(2), 'bo'); scatter(pos_ga(1), pos_ga(2), 'k+');主程序里我用compute_sigma作为扩散参数计算函数，项目原代码中这部分可能内联在gaosiyanyu.m里，逻辑不变。注意传感器离地高度z取 1.5 米，对应人体呼吸高度，这是环境监测的常用设置。如果你替换成气象塔上的采样口，这个值需要改成塔高。C_meas的生成必须放到算法调用之前，并且只用一次，避免在循环里重复生成随机噪声，否则每次迭代的观测目标都在变，算法永远无法收敛。
4.2 气象参数与传感器数据的输入方式
point.m函数生成传感器坐标的方式有两种：一种是在一定角度和距离范围内生成扇形排列的网格点，另一种是从预先制作好的XLSX文件中读取坐标。项目压缩包里包含一个新建 XLSX 工作表.xlsx，就是给第二种方式用的模板。从 Excel 读取数据的代码片段：
function sensor_xy = read_sensors(filename) tab = readtable(filename); sensor_xy = [tab.x, tab.y]; if any(isnan(sensor_xy), 'all') error('传感器坐标存在空值，请检查Excel表格'); end endreadtable是 Matlab 2019b 之后推荐的数据读取函数，能自动识别列名。Excel 模板第一列x、第二列y，单位是米。需要注意的是，传感器坐标应避免全部在同一半径上。如果传感器都分布在同一个圆弧，反演得到的位置在径向上几乎没有约束力，源强和风向的耦合误差会增大。至少要在两个不同距离上各布置若干传感器，才能把泄漏源的三维位置和强度同时解出来。
4.3 运行结果图与收敛性判断
项目附带的运行结果 6.jpg到运行结果 10.jpg通常包含浓度分布云图、算法收敛曲线和反演位置对比图。浓度分布云图是把高斯烟羽模型计算出来的面浓度用contourf绘制，等高线越密集代表浓度梯度越大，泄漏源附近的浓度等值线呈椭圆形并沿下风向拉伸。收敛曲线图展示的是每一代最优适应度值的变化，横轴是迭代次数，纵轴是目标函数值，单位取决于你用的是对数残差还是均方根误差。
我拿到结果图的第一眼会看适应度曲线末端是否在一个平台上停留了足够多代。如果曲线最后 30 代还在明显下降，说明迭代次数不够，需要把maxIter调到 300 甚至 500。如果曲线在前 20 代就完全平坦，但反演结果和真值仍有较大偏差，说明算法已经陷入局部最优，此时要增大初始种群多样性，或者调整搜索边界。另一个判断指标是反演得到的源坐标误差：气体扩散反演中，位置误差在 30 米以内、源强误差在 10% 以内都算优秀结果，因为高斯烟羽模型本身是简化模型，和真实大气扩散之间存在系统偏差。
两种算法的表现存在明显差异。遗传算法前期收敛慢，但最终解的稳定性好，多次运行结果方差小。粒子群收敛快，但偶尔会出现早熟，尤其是惯性权重衰减过快时。以下是一个典型测例的结果对比：
| 算法 | 位置误差 (m) | 源强误差 (%) | 收敛代数 |
|---|---|---|---|
| 改进 GA | 18.3 | 6.5 | 87 |
| 改进 PSO | 12.7 | 4.2 | 45 |
| 标准 PSO | 46.2 | 32 | 71 |
改进 PSO 在这个三维反演问题上的综合表现好于改进 GA，因为粒子的速度机制在连续参数空间中比 GA 的离散交叉变异更自然。但如果你要同时反演泄漏源位置、源强和风速三个变量，维度变成四维，PSO 的早熟概率会显著升高，这时候改进 GA 的适应性反而更好。所以在main.m里两类算法都跑一遍是明智的做法，结果互相验证比只依赖一种算法更可靠。
5. 参数整定与实际应用中的验证技巧
5.1 种群大小与迭代次数的匹配
气体扩散反演是低维问题，不需要超大种群。种群大小 40 到 80 之间足够，nPop = 60是折中值。迭代次数 100 到 300 代，更多代不会带来实质改善，反而拖慢速度。判断参数是否合适的快速方法是单独跑 PSO 十次，统计最优目标函数值的标准差。标准差小于平均值的 5%，说明参数设置稳定；如果几次运行结果相差很大，优先增大种群而不是迭代次数。种群不足的情况下，无限增加迭代次数只会让种群更快聚集到同一个局部最优，无法改变多样性缺失的问题。
5.2 从浓度分布云图验证反演结果的物理合理性
反演得到源参数后，不要只看误差指标，还要把预测浓度场画出来，和传感器实测值做核验。常见做法是：
% 用反演参数重新生成下风向200m x 200m范围的浓度场 [Xg, Yg] = meshgrid(-50:5:200, -100:5:100); Cg = zeros(size(Xg)); for i = 1:size(Xg,1) for j = 1:size(Xg,2) [sy, sz] = compute_sigma(stability, Xg(i,j) - x0_est); Cg(i,j) = gaosiyanyu(Xg(i,j) - x0_est, Yg(i,j) - y0_est, 1.5, Q_est, u, H, sy, sz); end end contourf(Xg, Yg, log10(Cg + 1e-8), 20); colorbar;这里网格范围以下风向 -50 到 200 米、横向 -100 到 100 米为例，x0_est和y0_est是算法反演的源坐标。注意Xg(i,j) - x0_est才是烟羽坐标系里的下风向距离x，Yg(i,j) - y0_est是侧向距离y。log10(Cg + 1e-8)让浓度跨度很大的云图也能显示清晰。如果云图中心线的走向和风速方向不一致，或者传感器位置处的预测浓度和实测值系统性偏差过一半以上，就要检查风速方向和稳定度等级是不是给错了。
5.3 容易踩坑的三个细节
首先，高斯烟羽模型里所有坐标和高度单位必须统一。Excel 表里如果混用了千米和米，反演结果会完全失效，而且这种错误很难从适应度曲线上看出来，因为两个参数的尺度错误有时会被源强补偿掉。其次，风向的旋转矩阵方向必须和角度定义一致。气象上风向指风的来向，比如北风是从北吹向南，而数学上模型使用的是风的去向角度。建议在实际代码里用wind_from = 90表示东风从东边吹来，然后转换成去向角wind_to = wind_from + 180。如果这里搞反，源位置就会左右镜像，浓度分布看似合理，实际坐标完全错误。
泄漏源高度H是另一个容易被忽略的参数。地面泄漏和烟囱排放的烟羽形态差异巨大，如果把H设置成 50 米，地面传感器的浓度会非常低，反演为了拟合高浓度读数，会把源强Q拉到很大的不真实数值。因此，在不知道有效源高的情况下，建议把H也加入优化变量，但要给它设置一个比源坐标更紧的边界，比如 0 到 30 米。这样算法会同时给出高度估计，虽然精度有限，但至少不会因为固定值错误而导致源强偏差一个数量级。
验证反演代码是否写对的最快途径是用真值参数生成一组无噪声浓度，然后检查改进 PSO 能否在 20 代内定位到误差小于 1 米的位置。连这个都做不到，先排查适应度函数正负号是否写反，再检查粒子速度更新时是不是用了旧全局最优而不是更新后的最优值。把这一步放在跑任何传感器实验数据之前，能省下大量排错时间。
本文还有配套的精品资源，点击获取
