瞬变电磁PSO反演：原理、实现与MATLAB代码

瞬变电磁PSO反演：原理、实现与MATLAB代码

一、瞬变电磁PSO反演基本原理

1.1 瞬变电磁反演问题

瞬变电磁法（TEM）通过观测地下介质对激励电磁场的瞬态响应来推断地下电性结构。反演问题可表述为：

目标函数：

min f(m) = ||d_obs - d_cal(m)||² + λR(m)

其中：

1.2 粒子群优化算法原理

PSO算法模拟鸟群觅食行为，每个粒子代表一个候选解：

v_i^{k+1} = w·v_i^k + c1·r1·(pbest_i - x_i^k) + c2·r2·(gbest - x_i^k)
x_i^{k+1} = x_i^k + v_i^{k+1}

其中：

二、MATLAB实现完整代码

2.1 主程序（main_TEM_PSO.m）

%% 瞬变电磁PSO反演主程序
clc; clear; close all;
warning('off', 'all');

%% 1. 参数设置
fprintf('瞬变电磁PSO反演程序\n');
fprintf('========================================\n');

% 反演参数
params.n_layers = 3;           % 地电模型层数
params.n_particles = 50;       % 粒子数量
params.max_iter = 200;         % 最大迭代次数
params.w = 0.9;                % 惯性权重
params.w_min = 0.4;            % 最小惯性权重
params.w_max = 0.9;            % 最大惯性权重
params.c1 = 2.0;               % 个体学习因子
params.c2 = 2.0;               % 社会学习因子
params.v_max = 0.2;            % 最大速度限制

% 模型参数范围
params.rho_min = 1;            % 最小电阻率 (Ω·m)
params.rho_max = 1000;         % 最大电阻率 (Ω·m)
params.h_min = 10;             % 最小厚度 (m)
params.h_max = 500;            % 最大厚度 (m)

% TEM观测参数
params.times = logspace(-5, -1, 30)';  % 观测时间 (s)
params.loop_size = 100;                % 发射回线边长 (m)
params.current = 10;                   % 发射电流 (A)

%% 2. 生成合成数据（用于测试）
fprintf('生成合成数据...\n');
true_model.rho = [100, 20, 500];      % 真实电阻率 (Ω·m)
true_model.h = [50, 150];             % 真实厚度 (m)

% 正演计算
[d_obs, times] = TEM_forward(true_model, params);
d_obs = d_obs .* (1 + 0.05 * randn(size(d_obs)));  % 添加5%噪声

fprintf('合成数据生成完成\n');
fprintf('  观测时间点数: %d\n', length(times));
fprintf('  时间范围: %.2e ~ %.2e s\n', min(times), max(times));

%% 3. PSO反演
fprintf('\n开始PSO反演...\n');
[best_model, best_fitness, convergence] = TEM_PSO_inversion(d_obs, params);

%% 4. 结果显示
fprintf('\n反演完成！\n');
fprintf('========================================\n');
fprintf('真实模型:\n');
fprintf('  电阻率: %s Ω·m\n', mat2str(true_model.rho));
fprintf('  厚度: %s m\n', mat2str(true_model.h));

fprintf('\n反演模型:\n');
fprintf('  电阻率: %s Ω·m\n', mat2str(best_model.rho));
fprintf('  厚度: %s m\n', mat2str(best_model.h));
fprintf('  最终适应度: %.6f\n', best_fitness);

%% 5. 可视化
plot_TEM_results(true_model, best_model, d_obs, times, convergence, params);

2.2 瞬变电磁正演计算（TEM_forward.m）

function [dBdt, times] = TEM_forward(model, params)
% 瞬变电磁一维正演计算（中心回线装置）
% 输入: model - 地电模型结构体
%       params - 参数结构体
% 输出: dBdt - 感应电动势 (V/m²)
%       times - 观测时间序列

% 提取参数
rho = model.rho;      % 各层电阻率
h = model.h;          % 各层厚度
times = params.times; % 观测时间
a = params.loop_size/2; % 回线半径
I = params.current;   % 发射电流

n_layers = length(rho);
n_times = length(times);
mu0 = 4 * pi * 1e-7;  % 真空磁导率

% 计算频率域响应
n_freq = 201;
freq = logspace(-3, 6, n_freq);  % 频率范围：1mHz - 1MHz
omega = 2 * pi * freq;

% 计算波数
k = zeros(n_layers, n_freq);
for i = 1:n_layers
    k(i, :) = sqrt(1i * omega * mu0 / rho(i));
end

% 计算表面阻抗（递归算法）
Z = zeros(1, n_freq);
Z(:) = sqrt(1i * omega * mu0 * rho(end));  % 最底层半空间

for i = n_layers-1:-1:1
    ki = k(i, :);
    Zi = sqrt(1i * omega * mu0 * rho(i));
    th = tanh(ki * h(i));
    Z = Zi .* (Z + Zi .* th) ./ (Zi + Z .* th);
end

% 计算频率域磁场
R = sqrt(a^2);  % 观测点距离
k0 = sqrt(1i * omega * mu0 / rho(1));
Hz_freq = (I * a^2) ./ (2 * (R^2 + a^2)^1.5) .* ...
    (1 + k0 .* R) .* exp(-k0 .* R);

% 频率域转时间域（正弦变换）
dBdt = zeros(n_times, 1);
for i = 1:n_times
    t = times(i);
    integrand = imag(Hz_freq .* omega) .* sin(omega * t);
    dBdt(i) = (2/pi) * trapz(omega, integrand);
end

% 转换为感应电动势
dBdt = -mu0 * dBdt;
end

2.3 PSO反演核心算法（TEM_PSO_inversion.m）

function [best_model, best_fitness, convergence] = TEM_PSO_inversion(d_obs, params)
% 瞬变电磁PSO反演
% 输入: d_obs - 观测数据
%       params - 参数结构体
% 输出: best_model - 最优模型
%       best_fitness - 最优适应度
%       convergence - 收敛历史

n_layers = params.n_layers;
n_particles = params.n_particles;
max_iter = params.max_iter;
n_params = 2 * n_layers - 1;  % n个电阻率 + (n-1)个厚度

% 初始化粒子群
particles = initialize_particles(n_particles, n_params, params);
velocities = zeros(n_particles, n_params);
pbest_positions = particles;
pbest_fitness = inf(n_particles, 1);
gbest_position = [];
gbest_fitness = inf;

% 收敛记录
convergence.best_fitness = zeros(max_iter, 1);
convergence.mean_fitness = zeros(max_iter, 1);
convergence.std_fitness = zeros(max_iter, 1);

% 进度条
fprintf('迭代进度: ');
progress_bar = waitbar(0, 'PSO反演进行中...');

% PSO主循环
for iter = 1:max_iter
    % 更新惯性权重（线性递减）
    w = params.w_max - (params.w_max - params.w_min) * iter / max_iter;
    
    % 评估所有粒子
    fitness = zeros(n_particles, 1);
    for i = 1:n_particles
        % 提取模型参数
        model = extract_model(particles(i, :), params);
        
        % 计算适应度
        fitness(i) = calculate_fitness(model, d_obs, params);
        
        % 更新个体最优
        if fitness(i) < pbest_fitness(i)
            pbest_fitness(i) = fitness(i);
            pbest_positions(i, :) = particles(i, :);
        end
    end
    
    % 更新全局最优
    [min_fitness, min_idx] = min(fitness);
    if min_fitness < gbest_fitness
        gbest_fitness = min_fitness;
        gbest_position = particles(min_idx, :);
    end
    
    % 记录收敛信息
    convergence.best_fitness(iter) = gbest_fitness;
    convergence.mean_fitness(iter) = mean(fitness);
    convergence.std_fitness(iter) = std(fitness);
    
    % 更新粒子速度和位置
    for i = 1:n_particles
        r1 = rand(1, n_params);
        r2 = rand(1, n_params);
        
        % 速度更新
        velocities(i, :) = w * velocities(i, :) + ...
            params.c1 * r1 .* (pbest_positions(i, :) - particles(i, :)) + ...
            params.c2 * r2 .* (gbest_position - particles(i, :));
        
        % 速度限制
        velocities(i, :) = min(max(velocities(i, :), -params.v_max), params.v_max);
        
        % 位置更新
        particles(i, :) = particles(i, :) + velocities(i, :);
        
        % 边界处理
        particles(i, :) = enforce_bounds(particles(i, :), params);
    end
    
    % 显示进度
    if mod(iter, 20) == 0
        waitbar(iter/max_iter, progress_bar, ...
            sprintf('迭代 %d/%d, 最优适应度: %.6f', iter, max_iter, gbest_fitness));
        fprintf('#');
    end
end

close(progress_bar);
fprintf('\n');

% 提取最优模型
best_model = extract_model(gbest_position, params);
best_fitness = gbest_fitness;
end

function particles = initialize_particles(n_particles, n_params, params)
% 初始化粒子群
% 使用拉丁超立方采样确保空间均匀性

particles = zeros(n_particles, n_params);
n_layers = params.n_layers;

for i = 1:n_particles
    % 电阻率参数（对数均匀分布）
    rho_params = 10.^(log10(params.rho_min) + ...
        rand(1, n_layers) * log10(params.rho_max/params.rho_min));
    
    % 厚度参数（对数均匀分布）
    h_params = 10.^(log10(params.h_min) + ...
        rand(1, n_layers-1) * log10(params.h_max/params.h_min));
    
    % 组合参数
    particles(i, :) = [rho_params, h_params];
end

% 添加一些特殊初始化
particles(1, 1:n_layers) = params.rho_min;  % 全最小电阻率
particles(2, 1:n_layers) = params.rho_max;  % 全最大电阻率
particles(3, 1:n_layers) = sqrt(params.rho_min * params.rho_max);  % 几何平均
end

function fitness = calculate_fitness(model, d_obs, params)
% 计算适应度函数
% 使用加权最小二乘 + 模型平滑约束

% 正演计算
d_cal = TEM_forward(model, params);

% 数据拟合项
data_misfit = sum(((d_obs - d_cal) ./ (abs(d_obs) + 1e-10)).^2);

% 模型平滑约束（一阶差分）
rho = model.rho;
h = model.h;

% 电阻率平滑
if length(rho) > 1
    rho_smooth = sum(diff(log10(rho)).^2);
else
    rho_smooth = 0;
end

% 厚度平滑
if length(h) > 1
    h_smooth = sum(diff(log10(h)).^2);
else
    h_smooth = 0;
end

% 总适应度
alpha = 0.1;  % 正则化参数
beta = 0.05;  % 厚度约束权重
fitness = data_misfit + alpha * rho_smooth + beta * h_smooth;
end

function model = extract_model(particle, params)
% 从粒子位置提取模型参数

n_layers = params.n_layers;
rho = particle(1:n_layers);
h = particle(n_layers+1:end);

% 确保物理合理性
rho = max(min(rho, params.rho_max), params.rho_min);
h = max(min(h, params.h_max), params.h_min);

model.rho = rho;
model.h = h;
end

function particle = enforce_bounds(particle, params)
% 强制边界约束

n_layers = params.n_layers;

% 电阻率边界
particle(1:n_layers) = max(min(particle(1:n_layers), params.rho_max), params.rho_min);

% 厚度边界
particle(n_layers+1:end) = max(min(particle(n_layers+1:end), params.h_max), params.h_min);
end

2.4 改进的PSO算法（改进策略）

function [best_model, best_fitness] = improved_PSO_TEM(d_obs, params)
% 改进的PSO算法（包含多种改进策略）
% 1. 自适应惯性权重
% 2. 动态学习因子
% 3. 混沌初始化
% 4. 变异操作

n_particles = params.n_particles;
max_iter = params.max_iter;
n_params = 2 * params.n_layers - 1;

% 混沌初始化（Logistic映射）
particles = chaotic_initialization(n_particles, n_params, params);

% 自适应参数
w_max = 0.9; w_min = 0.4;
c1_max = 2.5; c1_min = 1.5;
c2_max = 2.5; c2_min = 1.5;

% 改进PSO主循环
for iter = 1:max_iter
    % 动态调整参数
    w = w_max - (w_max - w_min) * iter / max_iter;
    c1 = c1_max - (c1_max - c1_min) * iter / max_iter;
    c2 = c2_min + (c2_max - c2_min) * iter / max_iter;
    
    % 引入变异操作（防止早熟）
    if rand() < 0.1  % 10%概率变异
        mutation_idx = randi(n_particles);
        particles(mutation_idx, :) = particles(mutation_idx, :) + ...
            0.1 * randn(1, n_params) .* (params.rho_max - params.rho_min);
        particles(mutation_idx, :) = enforce_bounds(particles(mutation_idx, :), params);
    end
    
    % 其他PSO步骤...
end
end

function particles = chaotic_initialization(n_particles, n_params, params)
% 混沌初始化（提高种群多样性）

particles = zeros(n_particles, n_params);
mu = 4;  % Logistic参数

for i = 1:n_particles
    % 生成混沌序列
    chaotic_seq = zeros(1, n_params);
    chaotic_seq(1) = rand();
    
    for j = 2:n_params
        chaotic_seq(j) = mu * chaotic_seq(j-1) * (1 - chaotic_seq(j-1));
    end
    
    % 映射到参数空间
    rho_params = params.rho_min + chaotic_seq(1:params.n_layers) * ...
        (params.rho_max - params.rho_min);
    
    h_params = params.h_min + chaotic_seq(params.n_layers+1:end) * ...
        (params.h_max - params.h_min);
    
    particles(i, :) = [rho_params, h_params];
end
end

2.5 可视化函数（plot_TEM_results.m）

function plot_TEM_results(true_model, inv_model, d_obs, times, convergence, params)
% 绘制反演结果

figure('Position', [100, 100, 1400, 800]);

% 子图1：数据拟合
subplot(2,3,1);
loglog(times, abs(d_obs), 'bo', 'MarkerSize', 8, 'LineWidth', 1.5);
hold on;
d_cal = TEM_forward(inv_model, params);
loglog(times, abs(d_cal), 'r-', 'LineWidth', 2);
grid on;
xlabel('时间 (s)', 'FontSize', 12);
ylabel('dB/dt (V/m²)', 'FontSize', 12);
title('数据拟合曲线', 'FontSize', 14);
legend('观测数据', '反演数据', 'Location', 'best');
set(gca, 'FontSize', 11);

% 子图2：收敛曲线
subplot(2,3,2);
semilogy(convergence.best_fitness, 'b-', 'LineWidth', 2);
hold on;
semilogy(convergence.mean_fitness, 'r--', 'LineWidth', 1.5);
fill([1:length(convergence.best_fitness), ...
    length(convergence.best_fitness):-1:1], ...
    [convergence.mean_fitness' + convergence.std_fitness', ...
    fliplr(convergence.mean_fitness' - convergence.std_fitness')], ...
    'r', 'FaceAlpha', 0.2, 'EdgeColor', 'none');
grid on;
xlabel('迭代次数', 'FontSize', 12);
ylabel('适应度值', 'FontSize', 12);
title('PSO收敛曲线', 'FontSize', 14);
legend('最优适应度', '平均适应度', '标准差范围', 'Location', 'best');
set(gca, 'FontSize', 11);

% 子图3：电阻率-深度模型
subplot(2,3,3);
plot_model_comparison(true_model, inv_model);
title('电阻率-深度模型对比', 'FontSize', 14);
xlabel('电阻率 (Ω·m)', 'FontSize', 12);
ylabel('深度 (m)', 'FontSize', 12);
grid on;
set(gca, 'FontSize', 11);

% 子图4：相对误差
subplot(2,3,4);
relative_error = abs(d_cal - d_obs) ./ abs(d_obs) * 100;
semilogx(times, relative_error, 'k-', 'LineWidth', 2);
hold on;
plot([min(times), max(times)], [5, 5], 'r--', 'LineWidth', 1.5);
grid on;
xlabel('时间 (s)', 'FontSize', 12);
ylabel('相对误差 (%)', 'FontSize', 12);
title('数据拟合相对误差', 'FontSize', 14);
legend('误差', '5%阈值', 'Location', 'best');
ylim([0, max(20, max(relative_error)*1.2)]);
set(gca, 'FontSize', 11);

% 子图5：参数搜索空间
subplot(2,3,5);
if length(inv_model.rho) >= 2
    scatter(inv_model.rho(1), inv_model.rho(2), 100, 'r', 'filled');
    hold on;
    scatter(true_model.rho(1), true_model.rho(2), 100, 'b', 'filled');
    xlabel('第一层电阻率 (Ω·m)', 'FontSize', 12);
    ylabel('第二层电阻率 (Ω·m)', 'FontSize', 12);
    title('参数搜索空间', 'FontSize', 14);
    legend('反演解', '真实解', 'Location', 'best');
    grid on;
end
set(gca, 'FontSize', 11);

% 子图6：残差分布
subplot(2,3,6);
residuals = d_cal - d_obs;
histogram(residuals, 20, 'FaceColor', [0.2, 0.6, 0.8]);
hold on;
xline(0, 'r-', 'LineWidth', 2);
xline(mean(residuals), 'g--', 'LineWidth', 2);
xlabel('残差', 'FontSize', 12);
ylabel('频数', 'FontSize', 12);
title('残差分布直方图', 'FontSize', 14);
legend('残差分布', '零线', sprintf('均值: %.2e', mean(residuals)), 'Location', 'best');
grid on;
set(gca, 'FontSize', 11);

% 保存图片
saveas(gcf, 'TEM_PSO_inversion_results.png');
end

function plot_model_comparison(true_model, inv_model)
% 绘制模型对比图

% 计算深度坐标
depth_true = [0, cumsum(true_model.h)];
depth_inv = [0, cumsum(inv_model.h)];

% 创建阶梯图
stairs([true_model.rho, true_model.rho(end)], [depth_true, depth_true(end)], ...
    'b-', 'LineWidth', 2);
hold on;
stairs([inv_model.rho, inv_model.rho(end)], [depth_inv, depth_inv(end)], ...
    'r--', 'LineWidth', 2);

set(gca, 'YDir', 'reverse');  % 深度向下增加
xlim([min([true_model.rho, inv_model.rho])*0.5, ...
    max([true_model.rho, inv_model.rho])*1.5]);
ylim([0, max([depth_true(end), depth_inv(end)])*1.1]);

legend('真实模型', '反演模型', 'Location', 'best');
grid on;
end

参考代码 瞬变电磁PSO反演 www.youwenfan.com/contentcsu/54815.html

三、PSO改进策略与最新研究

3.1 经典PSO存在的问题及改进

问题 改进策略 效果
易陷入局部最优 Lévy flight搜索策略 提高全局搜索能力
收敛速度慢 自适应惯性权重 平衡探索与开发
种群多样性不足 混沌初始化 提高搜索空间覆盖率
早熟收敛 变异算子引入 跳出局部最优

3.2 混合算法策略

  1. PSO-DLS组合算法

    • PSO全局搜索 + 阻尼最小二乘局部优化
    • 无需人工给定初始模型

  2. L-PSO算法

    • 引入Lévy flight长步长搜索
    • 提高跳出局部最优概率

  3. 重心反向学习PSO

    • 动态生成反向解并择优选择
    • 提高全局搜索能力

  4. PSO-FADBN网络

    • PSO优化深度信念网络
    • 因子分析降维 + DBN特征提取

3.3 最新研究进展

  1. 模型-数据混合驱动

    • 结合模型驱动与数据驱动优势
    • 计算效率提升75%，边界刻画更清晰

  2. 自适应差分优化

    • 差分进化 + 高斯牛顿法
    • 电阻率-极化率多参数同时提取

  3. 混沌PSO算法

    • Bernoulli映射混沌粒子群
    • 适用于2.5维水下模型反演

四、应用实例与代码扩展

4.1 矿井瞬变电磁反演

% 矿井全空间TEM反演
function mine_TEM_inversion()
    % 矿井特殊参数
    params.mine_depth = 500;      % 矿井深度
    params.tunnel_size = 5;       % 巷道尺寸
    params.water_content = 0.1;   % 含水率
    
    % 全空间正演修正
    function response = mine_TEM_forward(model, params)
        % 考虑全空间效应
        response = TEM_forward(model, params);
        % 全空间修正因子
        correction = calculate_fullspace_correction(params);
        response = response .* correction;
    end
    
    % 反演执行
    [best_model, fitness] = improved_PSO_TEM(d_obs, params);
end

4.2 多参数反演扩展

% 电阻率-极化率多参数PSO反演
function multi_parameter_PSO()
    % 扩展模型参数
    n_layers = 4;
    n_particles = 100;
    
    % 参数：电阻率 + 极化率 + 厚度
    n_params = 3 * n_layers - 1;
    
    % 多目标适应度函数
    function fitness = multi_objective_fitness(particle, d_obs, params)
        % 提取参数
        rho = particle(1:n_layers);
        m = particle(n_layers+1:2*n_layers);  % 极化率
        h = particle(2*n_layers+1:end);
        
        % 计算复电阻率响应
        response = complex_resistivity_forward(rho, m, h, params);
        
        % 多目标加权
        w1 = 0.7;  % 电阻率拟合权重
        w2 = 0.3;  % 极化率拟合权重
        
        fitness = w1 * norm(d_obs.rho - response.rho) + ...
                  w2 * norm(d_obs.m - response.m);
    end
end

五、实践建议与注意事项

5.1 参数设置建议

参数 推荐值 说明
粒子数量 30-100 问题复杂度越高，粒子数越多
最大迭代次数 100-500 根据收敛情况调整
惯性权重 0.4-0.9 线性递减策略
学习因子 c1=c2=2.0 平衡个体与社会学习
速度限制 0.1-0.3 防止粒子飞出搜索空间

5.2 收敛判断准则

  1. 相对变化率：|f_{k} - f_{k-1}|/f_{k-1} < ε
  2. 种群多样性：粒子位置标准差小于阈值
  3. 最大迭代次数：达到预设上限
  4. 计算时间限制：实际应用中的时间约束

5.3 常见问题与解决方案

问题 现象 解决方案
早熟收敛 适应度过早稳定 增加变异操作、混沌初始化
收敛缓慢 迭代数百次未收敛 调整惯性权重、学习因子
过拟合 训练好但验证差 增加正则化项、交叉验证
参数敏感 结果波动大 多次运行取统计结果

六、总结

瞬变电磁PSO反演作为一种全局优化方法，在解决非线性、多极值反演问题中表现出色。通过合理的算法改进和参数设置，可以显著提高反演精度和效率。

关键优势：

  1. 无需初始模型：相比线性化方法，PSO不依赖初始猜测
  2. 全局搜索能力：能够找到全局最优或近似最优解
  3. 并行计算友好：粒子评估可并行进行
  4. 灵活易扩展：易于与其他算法结合形成混合策略

应用领域：

参考文献：

  1. 徐正玉等. 瞬变电磁法非线性优化反演算法对比[J]. 吉林大学学报, 2022
  2. 矿井瞬变电磁PSO-DLS组合算法反演研究, 2019
  3. 瞬变电磁数据L-PSO反演方法, 2022
  4. 张继锋. 基于改进粒子群算法的SOTEM电场分量Ex反演, 2025
  5. 基于改进PSO-FADBN网络的瞬变电磁反演方法, 2025

专注于matlab/simulink,电子电路,编程

关于我们湘ICP备19023172号-4