大地电磁二维正演MATLAB程序
一、MATLAB实现
1.1 主程序
%% 大地电磁二维正演程序
% 功能：实现大地电磁法二维正演模拟
% 方法：有限差分法求解亥姆霍兹方程
clear; clc; close all;
fprintf('=== 大地电磁二维正演程序开始 ===\n');
%% 1. 模型参数设置
fprintf('设置模型参数...\n');
% 网格参数
nx = 100; % x方向网格数
nz = 50; % z方向网格数
dx = 50; % x方向网格间距 (m)
dz = 50; % z方向网格间距 (m)
% 物理常数
mu0 = 4*pi*1e-7; % 真空磁导率 (H/m)
omega = 2*pi*1; % 角频率 (rad/s)，对应1Hz
% 创建地电模型
fprintf('创建地电模型...\n');
rho = create_resistivity_model(nx, nz);
% 可视化模型
figure('Position', [100, 100, 800, 400]);
subplot(1, 2, 1);
imagesc(rho);
colorbar;
xlabel('X网格'); ylabel('Z网格');
title('电阻率模型 (\Omega\cdotm)');
axis xy;
% 计算电导率
sigma = 1 ./ rho;
%% 2. 有限差分法求解
fprintf('开始有限差分求解...\n');
% 构建系数矩阵
fprintf('构建系数矩阵...\n');
[A, b] = build_fd_matrix(nx, nz, dx, dz, sigma, omega, mu0);
% 求解线性方程组
fprintf('求解线性方程组...\n');
tic;
Ey = A \ b;
solve_time = toc;
fprintf('求解完成，耗时 %.2f 秒\n', solve_time);
% 重塑为二维网格
Ey_grid = reshape(Ey, [nz, nx]);
% 可视化电场分布
subplot(1, 2, 2);
imagesc(real(Ey_grid));
colorbar;
xlabel('X网格'); ylabel('Z网格');
title('实部电场分布 (V/m)');
axis xy;
%% 3. 计算视电阻率和相位
fprintf('计算视电阻率和相位...\n');
% 提取地表电场
Ey_surface = Ey_grid(1, :);
% 计算磁场（通过法拉第定律）
Hy_surface = zeros(1, nx);
for i = 2:nx-1
% 中心差分计算dEy/dx
dEy_dx = (Ey_grid(1, i+1) - Ey_grid(1, i-1)) / (2*dx);
Hy_surface(i) = dEy_dx / (1i * omega * mu0);
end
% 计算视电阻率
app_res = 1 ./ (omega * mu0) .* abs(Ey_surface ./ Hy_surface).^2;
% 计算相位
phase_rad = angle(Ey_surface ./ Hy_surface);
phase_deg = rad2deg(phase_rad);
% 可视化结果
figure('Position', [100, 100, 1200, 400]);
subplot(1, 3, 1);
plot(1:nx, real(Ey_surface), 'b-', 'LineWidth', 2);
xlabel('测点编号'); ylabel('电场实部 (V/m)');
title('地表电场实部');
grid on;
subplot(1, 3, 2);
semilogy(1:nx, app_res, 'r-', 'LineWidth', 2);
xlabel('测点编号'); ylabel('视电阻率 (\Omega\cdotm)');
title('视电阻率曲线');
grid on;
subplot(1, 3, 3);
plot(1:nx, phase_deg, 'g-', 'LineWidth', 2);
xlabel('测点编号'); ylabel('相位 (度)');
title('相位曲线');
grid on;
%% 4. 多频点正演
fprintf('进行多频点正演...\n');
% 频率范围
freqs = logspace(-2, 2, 20); % 0.01Hz 到 100Hz
app_res_freq = zeros(length(freqs), nx);
phase_freq = zeros(length(freqs), nx);
for f_idx = 1:length(freqs)
freq = freqs(f_idx);
omega = 2*pi*freq;
fprintf(' 频率 %.2f Hz (%.0f/%.0f)\n', freq, f_idx, length(freqs));
% 重新构建矩阵
[A, b] = build_fd_matrix(nx, nz, dx, dz, sigma, omega, mu0);
% 求解
Ey = A \ b;
Ey_grid = reshape(Ey, [nz, nx]);
% 计算视电阻率和相位
Ey_surface = Ey_grid(1, :);
% 计算磁场
Hy_surface = zeros(1, nx);
for i = 2:nx-1
dEy_dx = (Ey_grid(1, i+1) - Ey_grid(1, i-1)) / (2*dx);
Hy_surface(i) = dEy_dx / (1i * omega * mu0);
end
% 存储结果
app_res_freq(f_idx, :) = 1 ./ (omega * mu0) .* abs(Ey_surface ./ Hy_surface).^2;
phase_freq(f_idx, :) = rad2deg(angle(Ey_surface ./ Hy_surface));
end
% 可视化多频结果
figure('Position', [100, 100, 1000, 400]);
% 视电阻率断面
subplot(1, 2, 1);
imagesc(1:nx, freqs, log10(app_res_freq));
colorbar;
xlabel('测点编号'); ylabel('频率 (Hz)');
title('视电阻率断面 (log_{10}\Omega\cdotm)');
set(gca, 'YScale', 'log');
axis xy;
% 相位断面
subplot(1, 2, 2);
imagesc(1:nx, freqs, phase_freq);
colorbar;
xlabel('测点编号'); ylabel('频率 (Hz)');
title('相位断面 (度)');
set(gca, 'YScale', 'log');
axis xy;
%% 5. 与一维半空间模型对比
fprintf('与一维半空间模型对比...\n');
% 一维半空间解析解
rho_half = 100; % 半空间电阻率
depth = 50; % 深度 (m)
% 计算一维视电阻率
app_res_1d = zeros(length(freqs), 1);
phase_1d = zeros(length(freqs), 1);
for f_idx = 1:length(freqs)
freq = freqs(f_idx);
omega = 2*pi*freq;
% 波数
k = sqrt(1i * omega * mu0 / rho_half);
% 反射系数（自由表面）
R = (k*depth - 1) / (k*depth + 1);
% 视电阻率
app_res_1d(f_idx) = rho_half * abs(1 + R)^2;
% 相位
phase_1d(f_idx) = rad2deg(angle(1 + R));
end
% 对比结果
figure('Position', [100, 100, 800, 400]);
subplot(1, 2, 1);
hold on;
for i = 1:nx
plot(freqs, app_res_freq(:, i), 'b-', 'LineWidth', 0.5, 'Color', [0.7, 0.7, 1]);
end
plot(freqs, app_res_1d, 'r-', 'LineWidth', 3, 'DisplayName', '一维半空间');
xlabel('频率 (Hz)'); ylabel('视电阻率 (\Omega\cdotm)');
title('视电阻率对比');
set(gca, 'XScale', 'log');
legend('Location', 'best');
grid on;
subplot(1, 2, 2);
hold on;
for i = 1:nx
plot(freqs, phase_freq(:, i), 'b-', 'LineWidth', 0.5, 'Color', [0.7, 0.7, 1]);
end
plot(freqs, phase_1d, 'r-', 'LineWidth', 3, 'DisplayName', '一维半空间');
xlabel('频率 (Hz)'); ylabel('相位 (度)');
title('相位对比');
set(gca, 'XScale', 'log');
legend('Location', 'best');
grid on;
%% 6. 保存结果
fprintf('保存结果...\n');
% 保存工作空间变量
save('mt_2d_forward_results.mat', ...
'rho', 'sigma', 'freqs', 'app_res_freq', 'phase_freq', ...
'app_res_1d', 'phase_1d', 'nx', 'nz', 'dx', 'dz');
fprintf('\n=== 正演完成 ===\n');
fprintf('结果已保存到 mt_2d_forward_results.mat\n');
1.2 辅助函数：创建电阻率模型
function rho = create_resistivity_model(nx, nz)
% 创建二维电阻率模型
% 输入：
% nx, nz: x和z方向的网格数
% 输出：
% rho: 电阻率模型 (nz x nx)
% 初始化：背景电阻率 100 Ω·m
rho = ones(nz, nx) * 100;
% 第一层：高阻盖层 (深度0-10层)
rho(1:10, :) = 500;
% 第二层：低阻层 (深度10-30层)
rho(10:30, :) = 50;
% 第三层：高阻基底 (深度30层以下)
rho(30:end, :) = 1000;
% 添加低阻异常体 (深度15-25层，横向30-50网格)
rho(15:25, 30:50) = 10;
% 添加高阻异常体 (深度20-28层，横向70-85网格)
rho(20:28, 70:85) = 5000;
fprintf(' 模型尺寸: %d x %d 网格\n', nx, nz);
fprintf(' 背景电阻率: %.1f Ω·m\n', 100);
fprintf(' 盖层电阻率: %.1f Ω·m\n', 500);
fprintf(' 低阻层电阻率: %.1f Ω·m\n', 50);
fprintf(' 基底电阻率: %.1f Ω·m\n', 1000);
fprintf(' 低阻异常体: %.1f Ω·m\n', 10);
fprintf(' 高阻异常体: %.1f Ω·m\n', 5000);
end
1.3 辅助函数：构建有限差分矩阵
function [A, b] = build_fd_matrix(nx, nz, dx, dz, sigma, omega, mu0)
% 构建有限差分系数矩阵
% 输入：
% nx, nz: 网格数
% dx, dz: 网格间距
% sigma: 电导率 (nz x nx)
% omega: 角频率
% mu0: 真空磁导率
% 输出：
% A: 系数矩阵
% b: 右端项
% 总网格点数
N = nx * nz;
% 初始化稀疏矩阵
A = sparse(N, N);
b = zeros(N, 1);
% 填充矩阵
for iz = 1:nz
for ix = 1:nx
% 当前网格索引
idx = (iz-1)*nx + ix;
% 当前点的电导率
sig = sigma(iz, ix);
% 中心差分系数
coeff_xx = 1 / (dx^2);
coeff_zz = 1 / (dz^2);
coeff_c = 1i * omega * mu0 * sig;
% 主对角线
A(idx, idx) = -2*coeff_xx - 2*coeff_zz + coeff_c;
% x方向相邻点
if ix > 1
idx_left = (iz-1)*nx + (ix-1);
A(idx, idx_left) = coeff_xx;
end
if ix < nx
idx_right = (iz-1)*nx + (ix+1);
A(idx, idx_right) = coeff_xx;
end
% z方向相邻点
if iz > 1
idx_up = (iz-2)*nx + ix;
A(idx, idx_up) = coeff_zz;
end
if iz < nz
idx_down = iz*nx + ix;
A(idx, idx_down) = coeff_zz;
end
% 边界条件：地表电场为1
if iz == 1
b(idx) = 1;
end
end
end
% 应用Dirichlet边界条件
% 左边界
for iz = 1:nz
idx = (iz-1)*nx + 1;
A(idx, :) = 0;
A(idx, idx) = 1;
b(idx) = 1;
end
% 右边界
for iz = 1:nz
idx = iz*nx;
A(idx, :) = 0;
A(idx, idx) = 1;
b(idx) = 1;
end
% 底边界
for ix = 1:nx
idx = (nz-1)*nx + ix;
A(idx, :) = 0;
A(idx, idx) = 1;
b(idx) = 1;
end
end
1.4 扩展：TM模式正演
function [Hx_grid, app_res, phase] = tm_mode_forward(nx, nz, dx, dz, rho, freq)
% TM模式正演（磁场平行于走向）
% 输入：
% nx, nz: 网格数
% dx, dz: 网格间距
% rho: 电阻率模型
% freq: 频率
% 输出：
% Hx_grid: 磁场分布
% app_res: 视电阻率
% phase: 相位
omega = 2*pi*freq;
mu0 = 4*pi*1e-7;
% 计算电导率
sigma = 1 ./ rho;
% 构建系数矩阵（TM模式方程）
A = sparse(nx*nz, nx*nz);
b = zeros(nx*nz, 1);
for iz = 1:nz
for ix = 1:nx
idx = (iz-1)*nx + ix;
sig = sigma(iz, ix);
% TM模式方程：∇×(1/σ ∇×Hx) + iωμ0 Hx = 0
% 简化为中心差分形式
coeff_xx = 1 / (dx^2 * sig);
coeff_zz = 1 / (dz^2 * sig);
coeff_c = 1i * omega * mu0;
A(idx, idx) = -2*coeff_xx - 2*coeff_zz + coeff_c;
% 相邻点
if ix > 1
A(idx, idx-1) = coeff_xx;
end
if ix < nx
A(idx, idx+1) = coeff_xx;
end
if iz > 1
A(idx, idx-nx) = coeff_zz;
end
if iz < nz
A(idx, idx+nx) = coeff_zz;
end
% 边界条件
if iz == 1 || iz == nz || ix == 1 || ix == nx
A(idx, :) = 0;
A(idx, idx) = 1;
b(idx) = 1;
end
end
end
% 求解
Hx = A \ b;
Hx_grid = reshape(Hx, [nz, nx]);
% 计算视电阻率（TM模式）
Hx_surface = Hx_grid(1, :);
% 计算电场（通过安培定律）
Ex_surface = zeros(1, nx);
for i = 2:nx-1
dHx_dz = (Hx_grid(2, i) - Hx_grid(1, i)) / dz;
Ex_surface(i) = dHx_dz / (1i * omega * mu0);
end
% 视电阻率
app_res = 1 ./ (omega * mu0) .* abs(Ex_surface ./ Hx_surface).^2;
% 相位
phase = rad2deg(angle(Ex_surface ./ Hx_surface));
end
二、程序原理详解
2.1 大地电磁法基本原理
大地电磁法（Magnetotelluric, MT）是基于天然交变电磁场探测地球内部电性结构的地球物理方法。
基本方程：
- 麦克斯韦方程组（频域）：
- 亥姆霍兹方程（TE模式）：
2.2 有限差分离散化
将连续偏微分方程离散化为差分方程：
整理得：
2.3 视电阻率计算
TE模式：
其中：
TM模式：
参考代码 大地电磁二维正演程序 www.youwenfan.com/contentcnu/63160.html
三、性能优化与扩展
3.1 使用迭代求解器
%% 使用迭代求解器加速
function Ey = iterative_solver(A, b, tol, max_iter)
% 使用BiCGSTAB迭代求解器
[Ey, flag, relres, iter] = bicgstab(A, b, tol, max_iter);
if flag == 0
fprintf(' 迭代收敛于 %d 次迭代，相对残差 %.2e\n', iter(2), relres);
else
fprintf(' 迭代未收敛，标志: %d\n', flag);
end
end
3.2 自适应网格加密
%% 自适应网格加密
function [fine_grid, fine_rho] = adaptive_refinement(coarse_grid, coarse_rho, refinement_level)
% 在异常体附近加密网格
[nz, nx] = size(coarse_rho);
% 找到需要加密的区域（电阻率变化大的区域）
grad_rho = gradient(coarse_rho);
grad_magnitude = sqrt(grad_rho{1}.^2 + grad_rho{2}.^2);
% 确定加密区域
threshold = 0.5 * max(grad_magnitude(:));
refine_mask = grad_magnitude > threshold;
% 创建精细网格
fine_factor = 2^refinement_level;
fine_nz = nz * fine_factor;
fine_nx = nx * fine_factor;
fine_grid = zeros(fine_nz, fine_nx);
fine_rho = zeros(fine_nz, fine_nx);
% 插值到精细网格
[X_coarse, Y_coarse] = meshgrid(1:nx, 1:nz);
[X_fine, Y_fine] = meshgrid(linspace(1, nx, fine_nx), linspace(1, nz, fine_nz));
fine_rho = interp2(X_coarse, Y_coarse, coarse_rho, X_fine, Y_fine, 'cubic');
end
3.3 并行计算加速
%% 并行计算多频点
function [app_res_freq, phase_freq] = parallel_multi_freq(nx, nz, dx, dz, rho, freqs)
% 使用并行池加速多频点计算
if isempty(gcp('nocreate'))
parpool('local', 4); % 启动4个worker
end
app_res_freq = zeros(length(freqs), nx);
phase_freq = zeros(length(freqs), nx);
parfor f_idx = 1:length(freqs)
freq = freqs(f_idx);
omega = 2*pi*freq;
mu0 = 4*pi*1e-7;
sigma = 1 ./ rho;
% 构建并求解
[A, b] = build_fd_matrix(nx, nz, dx, dz, sigma, omega, mu0);
Ey = A \ b;
Ey_grid = reshape(Ey, [nz, nx]);
% 计算视电阻率和相位
Ey_surface = Ey_grid(1, :);
Hy_surface = zeros(1, nx);
for i = 2:nx-1
dEy_dx = (Ey_grid(1, i+1) - Ey_grid(1, i-1)) / (2*dx);
Hy_surface(i) = dEy_dx / (1i * omega * mu0);
end
app_res_freq(f_idx, :) = 1 ./ (omega * mu0) .* abs(Ey_surface ./ Hy_surface).^2;
phase_freq(f_idx, :) = rad2deg(angle(Ey_surface ./ Hy_surface));
end
end
四、实际应用建议
4.1 模型设计指南
| 地质结构 | 电阻率特征 | 建议网格设置 |
|---|---|---|
| 沉积盆地 | 低阻层（10-100 Ω·m） | 浅部网格加密 |
| 结晶基底 | 高阻（>1000 Ω·m） | 深部网格可放宽 |
| 断裂带 | 低阻异常（<10 Ω·m） | 异常体附近加密 |
| 矿化带 | 高阻异常（>1000 Ω·m） | 异常体附近加密 |
4.2 频率选择策略
- 浅部探测（0-500m）：使用高频（10-1000 Hz）
- 中深部探测（500-5000m）：使用中频（0.1-10 Hz）
- 深部探测（>5000m）：使用低频（0.001-0.1 Hz）
4.3 数值稳定性建议
- 网格比例：保持
，避免各向异性 - Courant条件：
，其中 为趋肤深度 - 边界条件：使用吸收边界或扩展网格减少边界反射
- 预处理：对系数矩阵进行预处理加速收敛
五、结果解释与地质推断
5.1 视电阻率曲线特征
| 曲线形态 | 地质含义 |
|---|---|
| 高阻平台 | 结晶基底或火成岩 |
| 低阻凹陷 | 沉积盆地或含水层 |
| 阶梯状下降 | 层状介质界面 |
| 局部异常 | 矿体或构造破碎带 |
5.2 相位曲线特征
- 相位超前（>45°）：指示低阻层覆盖
- 相位滞后（<45°）：指示高阻层覆盖
- 相位突变：指示电性界面
六、总结
本MATLAB程序实现了大地电磁二维正演的完整流程：
- 模型建立：支持复杂二维地电结构
- 数值求解：有限差分法求解亥姆霍兹方程
- 多频正演：覆盖宽频带探测
- 结果可视化：视电阻率和相位断面图
- 验证对比：与一维解析解对比验证
该程序可直接用于：
- 大地电磁法勘探设计
- 反演算法测试
- 地质结构解释
- 教学演示和科研分析
