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. 运行与调参建议
- 层数确定：如未知，可先固定层数（通常 3–5 层），或使用信息准则（AIC/BIC）在循环中比较。
- 参数归一化：若数量级差异大，可将电阻率、厚度做对数变换后在 PSO 中寻优。
- 改进策略：
- 引入 动态惯性权重 或 收缩因子 抑制早收敛；
- 在目标函数中加入 先验约束（如平滑度、已知钻孔信息）；
- 使用 并行计算（Parallel Computing Toolbox）加速适应度评估。
- 扩展到二维：将电阻率剖面网格化，粒子维度为网格单元电阻率，正演采用有限差分/有限元（如使用 MARE2DEM、ModEM 接口）。
6. 结果示例（控制台输出）
真实模型: ρ=[100 20 500], h=[500 1000]
反演模型: ρ=[102.3 19.4 487], h=[521 989]
适应度值: 0.028
视电阻率曲线对比图将直观展示拟合效果。