一维大地电磁反演MATLAB程序

一维大地电磁反演MATLAB程序。大地电磁反演是通过地表观测的电磁场数据来推断地下电性结构的重要方法。

1. 一维大地电磁正演模拟

function [app_res, phase] = MT1D_forward(rho, thickness, periods)
    % 一维大地电磁正演模拟
    % 输入：
    %   rho - 各层电阻率 (Ω·m) [n_layers x 1]
    %   thickness - 各层厚度 (m) [(n_layers-1) x 1]，最后一层为半无限空间
    %   periods - 周期数组 (s) [n_periods x 1]
    % 输出：
    %   app_res - 视电阻率 (Ω·m)
    %   phase - 相位 (度)
    
    mu0 = 4 * pi * 1e-7; % 真空磁导率
    n_layers = length(rho);
    n_periods = length(periods);
    
    app_res = zeros(n_periods, 1);
    phase = zeros(n_periods, 1);
    
    for i_period = 1:n_periods
        T = periods(i_period);
        omega = 2 * pi / T;
        
        % 计算波数
        k = zeros(n_layers, 1);
        for i_layer = 1:n_layers
            k(i_layer) = sqrt(1i * omega * mu0 / rho(i_layer));
        end
        
        % 从底层开始向上递推计算阻抗
        Z = zeros(n_layers, 1);
        Z(n_layers) = 1i * omega * mu0 / k(n_layers); % 底层阻抗
        
        for i_layer = n_layers-1:-1:1
            % 层间递推关系
            kh = k(i_layer) * thickness(i_layer);
            Z_down = Z(i_layer + 1);
            Z_up = 1i * omega * mu0 / k(i_layer);
            
            % 阻抗递推公式
            numerator = Z_up .* (Z_down + Z_up .* tanh(kh));
            denominator = Z_up + Z_down .* tanh(kh);
            Z(i_layer) = numerator ./ denominator;
        end
        
        % 地表阻抗和视电阻率
        Z_surface = Z(1);
        app_res(i_period) = (1 / (omega * mu0)) * abs(Z_surface)^2;
        phase(i_period) = atan2(imag(Z_surface), real(Z_surface)) * 180 / pi;
    end
end

2. 一维OCCAM反演算法

function [rho_inv, misfit_history, model_history] = MT1D_OCCAM_inversion(...
    periods, app_res_obs, phase_obs, app_res_err, phase_err, ...
    rho_initial, thickness, lambda_initial, target_misfit)
    % 一维大地电磁OCCAM反演
    % 输入：
    %   periods - 周期数组
    %   app_res_obs - 观测视电阻率
    %   phase_obs - 观测相位
    %   app_res_err - 视电阻率误差
    %   phase_err - 相位误差
    %   rho_initial - 初始模型电阻率
    %   thickness - 层厚度
    %   lambda_initial - 初始正则化参数
    %   target_misfit - 目标拟合差
    
    max_iterations = 50;
    n_layers = length(rho_initial);
    n_periods = length(periods);
    
    % 数据向量
    data_obs = [log10(app_res_obs); phase_obs];
    data_err = [app_res_err ./ app_res_obs / log(10); ...
                phase_err * pi / 180]; % 转换为弧度
    
    % 数据协方差矩阵
    Wd = diag(1 ./ data_err);
    
    % 模型协方差矩阵（平滑约束）
    L = compute_smoothness_matrix(n_layers);
    Wm = L' * L;
    
    rho_current = rho_initial;
    lambda = lambda_initial;
    
    misfit_history = zeros(max_iterations, 1);
    model_history = zeros(max_iterations, n_layers);
    
    for iter = 1:max_iterations
        % 正演计算当前模型响应
        [app_res_calc, phase_calc] = MT1D_forward(rho_current, thickness, periods);
        data_calc = [log10(app_res_calc); phase_calc];
        
        % 计算残差和拟合差
        residual = data_obs - data_calc;
        misfit = sqrt(mean((residual ./ data_err).^2));
        misfit_history(iter) = misfit;
        model_history(iter, :) = rho_current;
        
        fprintf('迭代 %d: 拟合差 = %.4f, lambda = %.2e\n', iter, misfit, lambda);
        
        % 检查收敛
        if misfit <= target_misfit
            fprintf('达到目标拟合差，反演收敛！\n');
            break;
        end
        
        if iter > 1 && abs(misfit_history(iter) - misfit_history(iter-1)) < 1e-6
            fprintf('拟合差变化很小，反演收敛！\n');
            break;
        end
        
        % 计算雅可比矩阵
        J = compute_jacobian(rho_current, thickness, periods, data_calc, data_err);
        
        % OCCAM反演核心方程
        A = J' * (Wd' * Wd) * J + lambda * Wm;
        b = J' * (Wd' * Wd) * residual;
        
        % 模型更新
        delta_rho = A \ b;
        rho_new = rho_current .* exp(delta_rho); % 在对数空间更新
        
        % 检查新模型
        rho_new = max(rho_new, 0.1);   % 电阻率下限
        rho_new = min(rho_new, 10000); % 电阻率上限
        
        % 计算新模型的拟合差
        [app_res_new, phase_new] = MT1D_forward(rho_new, thickness, periods);
        data_new = [log10(app_res_new); phase_new];
        residual_new = data_obs - data_new;
        misfit_new = sqrt(mean((residual_new ./ data_err).^2));
        
        % 更新lambda和模型
        if misfit_new < misfit
            rho_current = rho_new;
            lambda = lambda * 0.8; % 减小正则化参数
        else
            lambda = lambda * 2.0; % 增大正则化参数
        end
        
        % 防止lambda过小或过大
        lambda = max(lambda, 1e-10);
        lambda = min(lambda, 1e10);
    end
    
    rho_inv = rho_current;
    misfit_history = misfit_history(1:iter);
    model_history = model_history(1:iter, :);
end

function L = compute_smoothness_matrix(n_layers)
    % 计算平滑度矩阵（一阶差分）
    L = zeros(n_layers-1, n_layers);
    for i = 1:n_layers-1
        L(i, i) = -1;
        L(i, i+1) = 1;
    end
end

function J = compute_jacobian(rho, thickness, periods, data_calc, data_err)
    % 计算雅可比矩阵（数值差分法）
    n_layers = length(rho);
    n_periods = length(periods);
    n_data = 2 * n_periods;
    
    J = zeros(n_data, n_layers);
    delta = 1e-6; % 差分步长
    
    for i_layer = 1:n_layers
        rho_perturbed = rho;
        rho_perturbed(i_layer) = rho(i_layer) * (1 + delta);
        
        [app_res_pert, phase_pert] = MT1D_forward(rho_perturbed, thickness, periods);
        data_pert = [log10(app_res_pert); phase_pert];
        
        J(:, i_layer) = (data_pert - data_calc) / (rho(i_layer) * delta);
    end
    
    % 考虑数据误差
    for i = 1:n_data
        J(i, :) = J(i, :) / data_err(i);
    end
end

3. 完整的反演演示程序

function MT1D_inversion_demo()
    % 一维大地电磁反演演示程序
    
    close all; clear; clc;
    
    fprintf('=== 一维大地电磁反演演示 ===\n');
    
    % 1. 生成理论模型
    fprintf('生成理论模型...\n');
    [periods, app_res_true, phase_true, rho_true, thickness] = create_synthetic_model();
    
    % 2. 添加噪声生成观测数据
    fprintf('生成观测数据（添加噪声）...\n');
    [app_res_obs, phase_obs, app_res_err, phase_err] = ...
        add_noise_to_data(app_res_true, phase_true);
    
    % 3. 设置初始模型
    fprintf('设置初始模型...\n');
    rho_initial = create_initial_model(length(rho_true));
    
    % 4. 执行反演
    fprintf('开始反演...\n');
    lambda_initial = 10;
    target_misfit = 1.0;
    
    [rho_inv, misfit_history, model_history] = MT1D_OCCAM_inversion(...
        periods, app_res_obs, phase_obs, app_res_err, phase_err, ...
        rho_initial, thickness, lambda_initial, target_misfist);
    
    % 5. 显示结果
    plot_inversion_results(periods, app_res_true, phase_true, ...
        app_res_obs, phase_obs, rho_true, rho_inv, thickness, ...
        misfit_history, model_history);
    
    % 6. 计算反演统计量
    calculate_statistics(rho_true, rho_inv, app_res_true, phase_true, ...
        app_res_obs, phase_obs);
end

function [periods, app_res_true, phase_true, rho_true, thickness] = create_synthetic_model()
    % 创建理论模型（三层模型）
    
    % 模型参数：低阻层夹在两个高阻层之间
    rho_true = [100, 10, 1000]; % Ω·m
    thickness = [500, 1000];    % m
    
    % 周期范围
    periods = logspace(-3, 3, 30)'; % 0.001s 到 1000s
    
    % 正演计算理论响应
    [app_res_true, phase_true] = MT1D_forward(rho_true, thickness, periods);
end

function [app_res_obs, phase_obs, app_res_err, phase_err] = ...
    add_noise_to_data(app_res_true, phase_true)
    % 添加噪声模拟观测数据
    
    noise_level = 0.05; % 5%噪声
    
    % 视电阻率噪声（对数正态分布）
    app_res_noise = lognrnd(0, noise_level, size(app_res_true));
    app_res_obs = app_res_true .* app_res_noise;
    app_res_err = app_res_true * noise_level;
    
    % 相位噪声（高斯分布）
    phase_noise = normrnd(0, noise_level * 5, size(phase_true)); % 相位噪声稍大
    phase_obs = phase_true + phase_noise;
    phase_err = ones(size(phase_true)) * noise_level * 5;
end

function rho_initial = create_initial_model(n_layers)
    % 创建初始模型（均匀半空间）
    rho_initial = ones(n_layers, 1) * 100; % 100 Ω·m 均匀模型
end

function plot_inversion_results(periods, app_res_true, phase_true, ...
    app_res_obs, phase_obs, rho_true, rho_inv, thickness, ...
    misfit_history, model_history)
    % 绘制反演结果
    
    figure('Position', [100, 100, 1200, 900]);
    
    % 1. 数据拟合图
    subplot(2, 3, 1);
    loglog(periods, app_res_true, 'k-', 'LineWidth', 2, 'DisplayName', '理论值');
    hold on;
    errorbar(periods, app_res_obs, app_res_obs * 0.05, 'ro', ...
        'MarkerSize', 4, 'DisplayName', '观测值');
    [app_res_inv, phase_inv] = MT1D_forward(rho_inv, thickness, periods);
    loglog(periods, app_res_inv, 'b--', 'LineWidth', 2, 'DisplayName', '反演结果');
    xlabel('周期 (s)');
    ylabel('视电阻率 (\Omega\cdot m)');
    title('视电阻率拟合');
    legend('show');
    grid on;
    
    subplot(2, 3, 2);
    semilogx(periods, phase_true, 'k-', 'LineWidth', 2, 'DisplayName', '理论值');
    hold on;
    errorbar(periods, phase_obs, ones(size(phase_obs)) * 5, 'ro', ...
        'MarkerSize', 4, 'DisplayName', '观测值');
    semilogx(periods, phase_inv, 'b--', 'LineWidth', 2, 'DisplayName', '反演结果');
    xlabel('周期 (s)');
    ylabel('相位 (度)');
    title('相位拟合');
    legend('show');
    grid on;
    
    % 2. 模型对比图
    subplot(2, 3, 3);
    depth = calculate_depth_profile(thickness);
    plot_model_comparison(rho_true, rho_inv, depth);
    
    % 3. 拟合差收敛曲线
    subplot(2, 3, 4);
    plot(1:length(misfit_history), misfit_history, 'b-o', 'LineWidth', 2);
    xlabel('迭代次数');
    ylabel('拟合差');
    title('拟合差收敛曲线');
    grid on;
    
    % 4. 模型演化
    subplot(2, 3, 5);
    plot_model_evolution(model_history, depth);
    
    % 5. 参数灵敏度分析
    subplot(2, 3, 6);
    plot_sensitivity_analysis(rho_inv, thickness, periods);
end

function depth = calculate_depth_profile(thickness)
    % 计算深度剖面
    n_layers = length(thickness) + 1;
    depth = zeros(n_layers, 1);
    for i = 2:n_layers
        depth(i) = depth(i-1) + thickness(i-1);
    end
end

function plot_model_comparison(rho_true, rho_inv, depth)
    % 绘制模型对比图
    semilogx(rho_true, depth, 'k-', 'LineWidth', 3, 'DisplayName', '真实模型');
    hold on;
    semilogx(rho_inv, depth, 'r--', 'LineWidth', 2, 'DisplayName', '反演模型');
    set(gca, 'YDir', 'reverse');
    xlabel('电阻率 (\Omega\cdot m)');
    ylabel('深度 (m)');
    title('电性结构对比');
    legend('show');
    grid on;
end

function plot_model_evolution(model_history, depth)
    % 绘制模型演化过程
    n_iter = size(model_history, 1);
    colors = parula(n_iter);
    
    for i = 1:n_iter
        if i == 1
            semilogx(model_history(i, :), depth, '--', 'Color', colors(i, :), ...
                'LineWidth', 1, 'DisplayName', '初始模型');
        elseif i == n_iter
            semilogx(model_history(i, :), depth, '-', 'Color', colors(i, :), ...
                'LineWidth', 2, 'DisplayName', '最终模型');
        else
            semilogx(model_history(i, :), depth, ':', 'Color', colors(i, :), ...
                'LineWidth', 0.5);
        end
        hold on;
    end
    
    set(gca, 'YDir', 'reverse');
    xlabel('电阻率 (\Omega\cdot m)');
    ylabel('深度 (m)');
    title('模型演化过程');
    grid on;
end

function plot_sensitivity_analysis(rho, thickness, periods)
    % 参数灵敏度分析
    n_layers = length(rho);
    n_periods = length(periods);
    
    % 计算雅可比矩阵
    [app_res, phase] = MT1D_forward(rho, thickness, periods);
    data_calc = [log10(app_res); phase];
    data_err = [ones(n_periods, 1) * 0.05; ones(n_periods, 1) * 5 * pi / 180];
    
    J = compute_jacobian(rho, thickness, periods, data_calc, data_err);
    
    % 计算灵敏度
    sensitivity = sqrt(sum(J.^2, 1))';
    
    depth = calculate_depth_profile(thickness);
    plot(sensitivity, depth, 'g-', 'LineWidth', 2);
    set(gca, 'YDir', 'reverse');
    xlabel('灵敏度');
    ylabel('深度 (m)');
    title('参数灵敏度分析');
    grid on;
end

function calculate_statistics(rho_true, rho_inv, app_res_true, phase_true, ...
    app_res_obs, phase_obs)
    % 计算反演统计量
    
    [app_res_inv, phase_inv] = MT1D_forward(rho_inv, thickness, periods);
    
    % 模型误差
    model_rmse = sqrt(mean((log10(rho_inv) - log10(rho_true)).^2));
    
    % 数据拟合误差
    data_rmse_appres = sqrt(mean((app_res_inv - app_res_obs).^2)) / mean(app_res_obs);
    data_rmse_phase = sqrt(mean((phase_inv - phase_obs).^2)) / mean(phase_obs);
    
    fprintf('\n=== 反演统计结果 ===\n');
    fprintf('模型相对误差: %.4f\n', model_rmse);
    fprintf('视电阻率拟合相对误差: %.4f\n', data_rmse_appres);
    fprintf('相位拟合相对误差: %.4f\n', data_rmse_phase);
    fprintf('最终模型电阻率: %s Ω·m\n', mat2str(rho_inv', 3));
end

% 运行演示程序
if __name__ == "__main__"
    MT1D_inversion_demo();
end

4. 实际数据应用示例

function process_real_data()
    % 处理实际大地电磁数据示例
    
    % 读取实际数据（假设为TXT格式）
    % 格式：周期(s) 视电阻率(Ω·m) 视电阻率误差 相位(度) 相位误差
    data = load('mt_data.txt');
    periods = data(:, 1);
    app_res_obs = data(:, 2);
    app_res_err = data(:, 3);
    phase_obs = data(:, 4);
    phase_err = data(:, 5);
    
    % 设置反演参数
    n_layers = 20; % 反演层数
    thickness = ones(n_layers-1, 1) * 100; % 等厚度层
    rho_initial = ones(n_layers, 1) * 100; % 初始模型
    
    lambda_initial = 1;
    target_misfit = 1.0;
    
    % 执行反演
    [rho_inv, misfit_history, model_history] = MT1D_OCCAM_inversion(...
        periods, app_res_obs, phase_obs, app_res_err, phase_err, ...
        rho_initial, thickness, lambda_initial, target_misfit);
    
    % 显示结果
    depth = calculate_depth_profile(thickness);
    figure;
    semilogx(rho_inv, depth, 'b-', 'LineWidth', 2);
    set(gca, 'YDir', 'reverse');
    xlabel('电阻率 (\Omega\cdot m)');
    ylabel('深度 (m)');
    title('实际数据反演结果');
    grid on;
end

参考代码 一维大地电磁反演程序 www.3dddown.com/cna/65755.html

特点

1. 完整的正演模拟：基于传输矩阵法的精确正演
2. 稳健的反演算法：OCCAM反演方法，带正则化平滑约束
3. 全面的可视化：数据拟合、模型对比、收敛曲线等
4. 误差分析：灵敏度分析和统计评估
5. 实际应用：支持实际数据处理

建议

1. 初始模型选择：从均匀半空间开始反演
2. 正则化参数：根据数据质量调整lambda值
3. 层状模型：可根据探测深度需求调整层数和厚度
4. 数据预处理：确保观测数据的质量和误差估计准确