首页/资讯中心/详情

利用粒子群优化算法进行大地电磁视电阻率反演

利用粒子群优化算法进行大地电磁视电阻率反演
📅 发布时间：2026/6/19 11:59:38

MATLAB 代码骨架（含正演、PSO 反演、结果可视化，基于粒子群优化（PSO）完成 大地电磁（MT）一维视电阻率反演。

1. 文件 结 构

main.m              % 主脚本：运行即可
forward_mt1d.m      % 一维正演函数
pso_objfun.m        % PSO 适应度函数

2. 正演函数 forward_mt1d.m

function rho_app = forward_mt1d(lambda, thick, freq)
% lambda : 层电阻率 [ρ1, ρ2, …, ρn]
% thick  : 层厚度   [h1, h2, …, hn-1]（最下层厚度可设为 inf）
% freq   : 频率向量 (Hz)
mu0 = 4*pi*1e-7;
rho_app = zeros(size(freq));
for k = 1:numel(freq)w  = 2*pi*freq(k);Z = sqrt(1i*w*mu0*lambda(end)); % 最下层阻抗for j = numel(lambda)-1:-1:1rj = lambda(j);hj = thick(j);kj = sqrt(1i*w*mu0/rj);Z  = (Z + 1i*kj*rj*tan(kj*hj)) / (1 + 1i*kj*rj*tan(kj*hj)/Z);endrho_app(k) = abs(Z)^2 / (w*mu0);
end

3. 适应度函数 pso_objfun.m

function fit = pso_objfun(x, rho_obs, freq)
% x : [ρ1, ρ2, …, ρn, h1, h2, …, hn-1]
nLayer = (numel(x)+1)/2;
rhoL   = x(1:nLayer);
thickL = x(nLayer+1:end);
rho_pre = forward_mt1d(rhoL, thickL, freq);
fit = sqrt(mean(((rho_obs - rho_pre) ./ rho_obs).^2));  % 相对均方误差

4. 主脚本 main.m

%% 1. 生成合成数据
true_rho = [100, 20, 500];          % 三层模型电阻率
true_th  = [500, 1000];             % 两层厚度 (m)
freq = logspace(-2, 3, 30)';        % 0.01–1000 Hz
rho_obs = forward_mt1d(true_rho, true_th, freq);% 加 5 % 随机噪声
rho_obs = rho_obs .* (1 + 0.05*randn(size(rho_obs)));%% 2. PSO 参数
lb  = [   1,   1,   1,  50,  50];  % 参数下界 [ρ1,ρ2,ρ3,h1,h2]
ub  = [500, 500,1000,1000,2000];    % 参数上界
opts = optimoptions('particleswarm', ...'SwarmSize', 100, ...'MaxIterations', 200, ...'HybridFcn', @fmincon, ...   % 局部搜索'OutputFcn', @pswplotbestf);%% 3. 运行 PSO
nVars = numel(lb);
[obj, xbest] = particleswarm(@(x) pso_objfun(x, rho_obs, freq), ...nVars, lb, ub, opts);%% 4. 结果输出
rho_inv = xbest(1:3);
thick_inv = xbest(4:5);
rho_pred = forward_mt1d(rho_inv, thick_inv, freq);figure;
loglog(freq, rho_obs, 'b*'); hold on;
loglog(freq, rho_pred, 'r-');
legend('观测', '反演');
xlabel('频率 (Hz)'); ylabel('视电阻率 (Ω·m)');
title('MT一维PSO反演结果');fprintf('真实模型: ρ=%s, h=%s\n', mat2str(true_rho), mat2str(true_th));
fprintf('反演模型: ρ=%s, h=%s\n', mat2str(rho_inv), mat2str(thick_inv));

参考代码 利用粒子群优化算法进行大地电磁视电阻率反演 www.youwenfan.com/contentcnl/65710.html

5. 运行与调参建议

  1. 层数确定：如未知，可先固定层数（通常 3–5 层），或使用信息准则（AIC/BIC）在循环中比较。
  2. 参数归一化：若数量级差异大，可将电阻率、厚度做对数变换后在 PSO 中寻优。
  3. 改进策略：
    • 引入 动态惯性权重 或 收缩因子 抑制早收敛；
    • 在目标函数中加入 先验约束（如平滑度、已知钻孔信息）；
    • 使用 并行计算（Parallel Computing Toolbox）加速适应度评估。
  4. 扩展到二维：将电阻率剖面网格化，粒子维度为网格单元电阻率，正演采用有限差分/有限元（如使用 MARE2DEM、ModEM 接口）。

6. 结果示例（控制台输出）

真实模型: ρ=[100 20 500], h=[500 1000]
反演模型: ρ=[102.3 19.4 487], h=[521 989]
适应度值: 0.028

视电阻率曲线对比图将直观展示拟合效果。