当前位置：首页 > news >正文

MATLAB实现直流电法和大地电磁法的一维正演计算

news 2025/10/12 6:55:29

MATLAB实现直流电法和大地电磁法的一维正演计算。这两种方法是地球物理勘探中常用的电法勘探技术。

直流电法一维正演

%% 直流电法一维正演计算
function dc_1d_forward()clear; close all; clc;fprintf('=== 直流电法一维正演计算 ===\n\n');% 模型参数设置model = set_model_parameters();% 计算理论视电阻率[ab2, apparent_resistivity] = calculate_dc_response(model);% 显示结果display_results(model, ab2, apparent_resistivity);% 绘制结果plot_results(model, ab2, apparent_resistivity);
end%% 设置模型参数
function model = set_model_parameters()% 层参数: [厚度(m), 电阻率(Ω·m)]% 最后一层厚度设为Inf（半无限空间）model.layers = [10, 100;    % 第一层: 厚度10m, 电阻率100 Ω·m20, 500;    % 第二层: 厚度20m, 电阻率500 Ω·mInf, 50     % 第三层: 基岩, 电阻率50 Ω·m];% 测量装置参数model.electrode_type = 'Wenner';  % 装置类型: 'Wenner', 'Schlumberger', 'DipoleDipole'model.min_ab = 1;                 % 最小AB/2距离 (m)model.max_ab = 100;               % 最大AB/2距离 (m)model.num_measurements = 20;      % 测点数fprintf('模型参数设置:\n');fprintf('层数: %d\n', size(model.layers, 1));for i = 1:size(model.layers, 1)if isinf(model.layers(i, 1))fprintf('  第%d层: 基岩, 电阻率 = %.1f Ω·m\n', i, model.layers(i, 2));elsefprintf('  第%d层: 厚度 = %.1f m, 电阻率 = %.1f Ω·m\n', i, model.layers(i, 1), model.layers(i, 2));endendfprintf('装置类型: %s\n', model.electrode_type);
end%% 计算直流电法响应
function [ab2, apparent_resistivity] = calculate_dc_response(model)% 生成AB/2距离（对数均匀分布）ab2 = logspace(log10(model.min_ab), log10(model.max_ab), model.num_measurements)';apparent_resistivity = zeros(size(ab2));for i = 1:length(ab2)switch model.electrode_typecase 'Wenner'% 温纳装置apparent_resistivity(i) = wenner_apparent_resistivity(model, ab2(i));case 'Schlumberger'% 施伦贝谢装置apparent_resistivity(i) = schlumberger_apparent_resistivity(model, ab2(i));case 'DipoleDipole'% 偶极-偶极装置apparent_resistivity(i) = dipole_dipole_apparent_resistivity(model, ab2(i));endend
end%% 温纳装置视电阻率计算
function rho_a = wenner_apparent_resistivity(model, AB2)% 温纳装置: AM = MN = NB = AB/3AM = AB2 / 3;rho_a = calculate_potential(model, AM) * 2 * pi * AM;
end%% 施伦贝谢装置视电阻率计算
function rho_a = schlumberger_apparent_resistivity(model, AB2)% 施伦贝谢装置: MN << ABMN = AB2 / 50;  % MN远小于ABAM = (AB2 - MN) / 2;potential_diff = calculate_potential(model, AM) - calculate_potential(model, AM + MN);rho_a = pi * AM * (AM + MN) * potential_diff / MN;
end%% 偶极-偶极装置视电阻率计算
function rho_a = dipole_dipole_apparent_resistivity(model, AB2)% 偶极-偶极装置r = AB2;  % 偶极中心距离rho_a = calculate_potential(model, r) * pi * r^3 / (1 * 1);  % 假设偶极矩为1
end%% 计算电位
function potential = calculate_potential(model, r)% 使用线性滤波法计算层状介质中的电位% 这里使用简化算法，实际应用中可能需要更复杂的数值方法n_layers = size(model.layers, 1);resistivities = model.layers(:, 2);thicknesses = model.layers(1:end-1, 1);% 计算反射系数序列reflection_coeffs = calculate_reflection_coefficients(resistivities);% 计算核函数kernel = calculate_kernel_function(r, thicknesses, reflection_coeffs);potential = kernel / (2 * pi * resistivities(1));
end%% 计算反射系数
function R = calculate_reflection_coefficients(resistivities)n = length(resistivities);R = zeros(n-1, 1);for i = 1:n-1k_i = (resistivities(i+1) - resistivities(i)) / (resistivities(i+1) + resistivities(i));R(i) = k_i;end
end%% 计算核函数（简化版本）
function kernel = calculate_kernel_function(r, thicknesses, reflection_coeffs)% 简化的一维核函数计算% 实际应用中可能需要使用线性滤波法或数值积分kernel = 1/r;  % 均匀半空间项% 添加层状介质修正（简化处理）for i = 1:length(reflection_coeffs)depth = sum(thicknesses(1:i));kernel = kernel + reflection_coeffs(i) / sqrt(r^2 + (2*depth)^2);end
end%% 显示结果
function display_results(model, ab2, apparent_resistivity)fprintf('\n=== 计算结果 ===\n');fprintf('AB/2(m)\t视电阻率(Ω·m)\n');fprintf('-------------------\n');for i = 1:length(ab2)fprintf('%.1f\t%.1f\n', ab2(i), apparent_resistivity(i));end
end%% 绘制结果
function plot_results(model, ab2, apparent_resistivity)figure('Position', [100, 100, 1200, 500]);% 子图1: 理论模型subplot(1, 2, 1);plot_model(model);title('一维地电模型');% 子图2: 视电阻率曲线subplot(1, 2, 2);semilogx(ab2, apparent_resistivity, 'b-o', 'LineWidth', 2, 'MarkerSize', 4);xlabel('AB/2 (m)');ylabel('视电阻率 (\Omega\cdotm)');title([model.electrode_type '装置视电阻率曲线']);grid on;% 添加真实电阻率参考线hold on;colors = ['r', 'g', 'b', 'm', 'c'];for i = 1:size(model.layers, 1)yline(model.layers(i, 2), '--', ...sprintf('第%d层: %.1f Ω·m', i, model.layers(i, 2)), ...'Color', colors(mod(i-1, length(colors))+1));endlegend('视电阻率', 'Location', 'best');
end%% 绘制地电模型
function plot_model(model)depth = 0;max_resistivity = max(model.layers(:, 2));for i = 1:size(model.layers, 1)if i == size(model.layers, 1)% 最后一层rectangle('Position', [0, depth, max_resistivity*1.2, 50], ...'FaceColor', [0.8 0.8 1], 'EdgeColor', 'k');text(max_resistivity*0.6, depth+25, ...sprintf('电阻率: %.1f Ω·m', model.layers(i, 2)), ...'HorizontalAlignment', 'center');else% 中间层rectangle('Position', [0, depth, max_resistivity*1.2, model.layers(i, 1)], ...'FaceColor', [0.8 0.8 1], 'EdgeColor', 'k');text(max_resistivity*0.6, depth+model.layers(i, 1)/2, ...sprintf('厚度: %.1f m\n电阻率: %.1f Ω·m', ...model.layers(i, 1), model.layers(i, 2)), ...'HorizontalAlignment', 'center');depth = depth + model.layers(i, 1);endendset(gca, 'YDir', 'reverse');xlabel('电阻率 (\Omega\cdotm)');ylabel('深度 (m)');xlim([0, max_resistivity*1.2]);ylim([0, depth+50]);
end

大地电磁法一维正演

%% 大地电磁法一维正演计算
function mt_1d_forward()clear; close all; clc;fprintf('=== 大地电磁法一维正演计算 ===\n\n');% 模型参数设置model = set_mt_model_parameters();% 计算MT响应[periods, apparent_resistivity, phase] = calculate_mt_response(model);% 显示结果display_mt_results(model, periods, apparent_resistivity, phase);% 绘制结果plot_mt_results(model, periods, apparent_resistivity, phase);
end%% 设置MT模型参数
function model = set_mt_model_parameters()% 层参数: [厚度(m), 电阻率(Ω·m)]model.layers = [1000, 100;   % 第一层: 沉积层5000, 10;    % 第二层: 低阻层Inf, 1000    % 第三层: 高阻基底];% 频率范围 (Hz)model.min_freq = 1e-4;   % 0.0001 Hzmodel.max_freq = 1e3;    % 1000 Hzmodel.num_periods = 30;  % 周期点数fprintf('MT模型参数设置:\n');fprintf('层数: %d\n', size(model.layers, 1));for i = 1:size(model.layers, 1)if isinf(model.layers(i, 1))fprintf('  第%d层: 基底, 电阻率 = %.1f Ω·m\n', i, model.layers(i, 2));elsefprintf('  第%d层: 厚度 = %.1f m, 电阻率 = %.1f Ω·m\n', i, model.layers(i, 1), model.layers(i, 2));endend
end%% 计算MT响应
function [periods, apparent_resistivity, phase] = calculate_mt_response(model)% 生成周期序列（对数均匀分布）frequencies = logspace(log10(model.max_freq), log10(model.min_freq), model.num_periods);periods = 1 ./ frequencies;apparent_resistivity = zeros(size(periods));phase = zeros(size(periods));for i = 1:length(periods)[apparent_resistivity(i), phase(i)] = calculate_impedance(model, periods(i));end
end%% 计算阻抗
function [apparent_resistivity, phase] = calculate_impedance(model, T)% 使用递推算法计算层状介质中的MT阻抗omega = 2 * pi / T;  % 角频率mu0 = 4 * pi * 1e-7; % 真空磁导率n_layers = size(model.layers, 1);resistivities = model.layers(:, 2);thicknesses = model.layers(1:end-1, 1);% 从最底层开始递推Z = sqrt(1i * omega * mu0 * resistivities(end));for k = n_layers-1:-1:1rho_k = resistivities(k);h_k = thicknesses(k);% 波数k_k = sqrt(1i * omega * mu0 / rho_k);% 层阻抗Z_k = sqrt(1i * omega * mu0 * rho_k);% 递推公式Z = Z_k * (Z * cosh(k_k * h_k) + Z_k * sinh(k_k * h_k)) / ...(Z_k * cosh(k_k * h_k) + Z * sinh(k_k * h_k));end% 计算视电阻率和相位apparent_resistivity = (abs(Z)^2) / (omega * mu0);phase = atan2(imag(Z), real(Z)) * 180 / pi;
end%% 显示MT结果
function display_mt_results(model, periods, apparent_resistivity, phase)fprintf('\n=== MT计算结果 ===\n');fprintf('周期(s)\t视电阻率(Ω·m)\t相位(°)\n');fprintf('--------------------------------\n');for i = 1:length(periods)fprintf('%.2e\t%.1f\t\t%.1f\n', periods(i), apparent_resistivity(i), phase(i));end
end%% 绘制MT结果
function plot_mt_results(model, periods, apparent_resistivity, phase)figure('Position', [100, 100, 1200, 800]);% 子图1: 理论模型subplot(2, 2, 1);plot_mt_model(model);title('一维MT地电模型');% 子图2: 视电阻率曲线subplot(2, 2, 2);loglog(periods, apparent_resistivity, 'b-o', 'LineWidth', 2, 'MarkerSize', 4);xlabel('周期 (s)');ylabel('视电阻率 (\Omega\cdotm)');title('MT视电阻率曲线');grid on;% 子图3: 相位曲线subplot(2, 2, 3);semilogx(periods, phase, 'r-s', 'LineWidth', 2, 'MarkerSize', 4);xlabel('周期 (s)');ylabel('相位 (°)');title('MT相位曲线');grid on;% 子图4: 视电阻率和相位的联合图subplot(2, 2, 4);yyaxis left;loglog(periods, apparent_resistivity, 'b-o', 'LineWidth', 2);ylabel('视电阻率 (\Omega\cdotm)');yyaxis right;semilogx(periods, phase, 'r-s', 'LineWidth', 2);ylabel('相位 (°)');xlabel('周期 (s)');title('MT响应曲线');grid on;legend('视电阻率', '相位', 'Location', 'best');
end%% 绘制MT地电模型
function plot_mt_model(model)depth = 0;max_resistivity = max(model.layers(:, 2));for i = 1:size(model.layers, 1)if i == size(model.layers, 1)% 最后一层rectangle('Position', [0, depth, max_resistivity*1.2, max(1000, depth*0.2)], ...'FaceColor', [1 0.8 0.8], 'EdgeColor', 'k');text(max_resistivity*0.6, depth+max(1000, depth*0.2)/2, ...sprintf('电阻率: %.1f Ω·m', model.layers(i, 2)), ...'HorizontalAlignment', 'center');else% 中间层rectangle('Position', [0, depth, max_resistivity*1.2, model.layers(i, 1)], ...'FaceColor', [1 0.8 0.8], 'EdgeColor', 'k');text(max_resistivity*0.6, depth+model.layers(i, 1)/2, ...sprintf('厚度: %.1f m\n电阻率: %.1f Ω·m', ...model.layers(i, 1), model.layers(i, 2)), ...'HorizontalAlignment', 'center');depth = depth + model.layers(i, 1);endendset(gca, 'YDir', 'reverse');xlabel('电阻率 (\Omega\cdotm)');ylabel('深度 (m)');xlim([0, max_resistivity*1.2]);
end%% 主程序 - 运行两个正演计算
function main()fprintf('地球物理电法一维正演计算程序\n\n');% 运行直流电法正演dc_1d_forward();fprintf('\n\n');% 运行大地电磁法正演mt_1d_forward();
end% 运行主程序
main();

特点：

直流电法部分：

支持多种装置：温纳、施伦贝谢、偶极-偶极
层状模型：可定义多层地电结构
可视化：显示地电模型和视电阻率曲线
理论计算：基于反射系数法的简化正演

大地电磁法部分：

频域响应：计算不同周期的MT响应
阻抗计算：使用递推算法计算层状介质阻抗
完整参数：输出视电阻率和相位
多图显示：分别显示模型、视电阻率、相位曲线

参考代码

使用：

修改模型参数：在set_model_parameters函数中调整层参数
选择装置类型：直流电法支持三种常用装置
调整测量参数：设置合适的测量范围和点数
运行分析：程序会自动计算并显示结果

扩展：

添加数据反演功能
支持更复杂的二维、三维模型
添加噪声模拟和数据处理
集成实际野外数据导入

参考

直流电测深正演报告. 百度文库, 2023. wenku.baidu.com/view/f9f53fc8adf8941ea76e58fafab069dc51224775.html
代码直流电法和大地电磁一维正演计算 www.youwenfan.com/contentcsi/65321.html
大地电磁一维正演理论. 中国知网, 2025.

查看全文

http://www.dtcms.com/a/469536.html