大地电磁正演计算：从理论到实现

在地球物理勘探中，大地电磁测深（MT）是一种重要的技术手段，而TE（横电）和TM（横磁）模式的正演计算又是其中的核心环节。

大地电磁正演计算：从理论到实现

TE与TM模式的基本原理

大地电磁测深利用天然交变电磁场来探测地下电性结构。由于地下介质通常表现为二维或三维非均匀性，我们需要分别计算TE和TM两种极化模式的电磁场响应，以全面了解地质构造。

极化模式	电场方向	磁场方向	物理特性	分辨率特点
TE模式 (横电模式)	平行于走向	在走向和垂直平面内	电场平行构造走向	纵向分辨率较高
TM模式 (横磁模式)	在走向和垂直平面内	平行于走向	磁场平行构造走向	横向分辨率较高

在实际大地电磁测深工作中，由于地下介质的非均匀性，TE和TM极化模式的视电阻率曲线往往表现出很大差异。理解这两种模式的响应规律，对于正确认识和判断实际地质结构至关重要。

正演计算的数学基础

控制方程

从麦克斯韦方程组出发，可以推导出大地电磁二维正演所满足的偏微分方程边值问题。

对于TE模式，控制方程可以表示为：

\(\nabla \cdot (\mu^{-1} \nabla E_y) + i\omega \sigma E_y = 0\)

对于TM模式，控制方程为：

\(\nabla \cdot (\sigma^{-1} \nabla H_y) + i\omega \mu H_y = 0\)

其中(\(E_y\))和(\(H_y\))分别是垂直于剖面方向的电场和磁场分量，(\(\sigma\))是电导率，(\(\mu\))是磁导率，(\(\omega\))是角频率。

数值离散方法

大地电磁正演主要采用数值模拟方法，常用的有：

• 有限单元法（FEM）：适合复杂地形和复杂地质条件
• 有限差分法（FDTD）：实现相对简单
• 谱元法（SEM）：计算精度高，收敛快

其中，有限单元法应用最为广泛。为提高计算精度，可采用双线性插值或双二次插值技术。

MATLAB实现框架

基于有限单元法的TE和TM模式正演计算MATLAB实现框架：

主程序框架

function [apparent_resistivity_TE, apparent_resistivity_TM] = MT2D_Forward(model, frequencies)
% 二维大地电磁TE和TM模式正演计算
% 输入参数：
%   model: 结构体，包含电导率模型、网格信息等
%   frequencies: 频率向量
% 输出参数：
%   apparent_resistivity_TE: TE模式的视电阻率
%   apparent_resistivity_TM: TM模式的视电阻率

    nFreq = length(frequencies);
    apparent_resistivity_TE = zeros(size(model.x, 1), nFreq);
    apparent_resistivity_TM = zeros(size(model.x, 1), nFreq);
    
    % 并行循环处理不同频率
    parfor iFreq = 1:nFreq
        freq = frequencies(iFreq);
        omega = 2 * pi * freq;
        
        % TE模式正演计算
        E_field_TE = TE_Forward(model, omega);
        
        % TM模式正演计算  
        H_field_TM = TM_Forward(model, omega);
        
        % 计算视电阻率
        apparent_resistivity_TE(:, iFreq) = calculate_ApparentResistivity_TE(E_field_TE, model, omega);
        apparent_resistivity_TM(:, iFreq) = calculate_ApparentResistivity_TM(H_field_TM, model, omega);
    end
end

TE模式正演计算

function E_field = TE_Forward(model, omega)
% TE模式正演计算

    % 生成网格（矩形或三角形）
    if strcmp(model.mesh_type, 'rectangular')
        [points, elements] = generate_RectangularMesh(model.x, model.y);
    else
        [points, elements] = generate_TriangularMesh(model.x, model.y);
    end
    
    % 组装总体刚度矩阵和质量矩阵
    nNodes = size(points, 1);
    K = sparse(nNodes, nNodes); % 刚度矩阵
    M = sparse(nNodes, nNodes); % 质量矩阵
    
    for iElem = 1:size(elements, 1)
        % 获取单元节点坐标和电导率
        elemNodes = elements(iElem, :);
        nodeCoords = points(elemNodes, :);
        sigma_elem = get_ElementConductivity(model, elemNodes);
        
        % 计算单元矩阵
        [Ke, Me] = compute_TE_ElementMatrices(nodeCoords, sigma_elem, omega);
        
        % 组装到总体矩阵
        K(elemNodes, elemNodes) = K(elemNodes, elemNodes) + Ke;
        M(elemNodes, elemNodes) = M(elemNodes, elemNodes) + Me;
    end
    
    % 施加边界条件
    [K, M] = apply_TE_BoundaryConditions(K, M, points, model);
    
    % 求解线性方程组: (K + iωμM) * E = source
    source = create_TE_Source(points, model);
    A = K + 1i * omega * mu0 * M;
    
    % 使用迭代法求解大型稀疏线性方程组
    E_field = bicgstab(A, source, 1e-10, 1000);
end

TM模式正演计算

function H_field = TM_Forward(model, omega)
% TM模式正演计算

    % 生成网格
    if strcmp(model.mesh_type, 'rectangular')
        [points, elements] = generate_RectangularMesh(model.x, model.y);
    else
        [points, elements] = generate_TriangularMesh(model.x, model.y);
    end
    
    % 组装总体矩阵（TM模式使用电导率的倒数）
    nNodes = size(points, 1);
    K_TM = sparse(nNodes, nNodes);
    M_TM = sparse(nNodes, nNodes);
    
    for iElem = 1:size(elements, 1)
        elemNodes = elements(iElem, :);
        nodeCoords = points(elemNodes, :);
        sigma_elem = get_ElementConductivity(model, elemNodes);
        rho_elem = 1 ./ sigma_elem; % 电阻率
        
        % 计算TM模式单元矩阵
        [Ke, Me] = compute_TM_ElementMatrices(nodeCoords, rho_elem, omega);
        
        K_TM(elemNodes, elemNodes) = K_TM(elemNodes, elemNodes) + Ke;
        M_TM(elemNodes, elemNodes) = M_TM(elemNodes, elemNodes) + Me;
    end
    
    % 施加TM模式边界条件
    [K_TM, M_TM] = apply_TM_BoundaryConditions(K_TM, M_TM, points, model);
    
    % 求解线性方程组
    source_TM = create_TM_Source(points, model);
    A_TM = K_TM + 1i * omega * mu0 * M_TM;
    
    H_field = bicgstab(A_TM, source_TM, 1e-10, 1000);
end

网格生成与单元矩阵计算

function [points, elements] = generate_RectangularMesh(x, y)
% 生成矩形网格
    [X, Y] = meshgrid(x, y);
    points = [X(:), Y(:)];
    
    nX = length(x); nY = length(y);
    elements = [];
    for i = 1:nY-1
        for j = 1:nX-1
            n1 = (i-1)*nX + j;
            n2 = (i-1)*nX + j+1;
            n3 = i*nX + j+1;
            n4 = i*nX + j;
            elements = [elements; [n1, n2, n3, n4]];
        end
    end
end

function [Ke, Me] = compute_TE_ElementMatrices(nodeCoords, sigma, omega)
% 计算TE模式单元刚度矩阵和质量矩阵
    nNodes = size(nodeCoords, 1);
    Ke = zeros(nNodes, nNodes);
    Me = zeros(nNodes, nNodes);
    
    % 获取形函数及其导数
    [N, dN_dxi, dN_deta, weights, gp] = get_ShapeFunctions(nodeCoords);
    
    % 数值积分（高斯积分）
    for iGauss = 1:length(weights)
        % 计算雅可比矩阵
        J = [dN_dxi(iGauss, :)' * nodeCoords(:, 1), ...
             dN_dxi(iGauss, :)' * nodeCoords(:, 2);
             dN_deta(iGauss, :)' * nodeCoords(:, 1), ...
             dN_deta(iGauss, :)' * nodeCoords(:, 2)];
         
        detJ = det(J);
        invJ = inv(J);
        
        % 形函数全局导数
        dN_dx = invJ(1,1) * dN_dxi(iGauss, :) + invJ(1,2) * dN_deta(iGauss, :);
        dN_dy = invJ(2,1) * dN_dxi(iGauss, :) + invJ(2,2) * dN_deta(iGauss, :);
        
        % 计算单元矩阵
        for i = 1:nNodes
            for j = 1:nNodes
                Ke(i, j) = Ke(i, j) + (1/mu0) * (dN_dx(i)*dN_dx(j) + dN_dy(i)*dN_dy(j)) * weights(iGauss) * detJ;
                Me(i, j) = Me(i, j) + (1i * omega * sigma) * N(iGauss, i) * N(iGauss, j) * weights(iGauss) * detJ;
            end
        end
    end
end

视电阻率计算

function apparent_resistivity = calculate_ApparentResistivity_TE(E_field, model, omega)
% 计算TE模式视电阻率
    mu0 = 4 * pi * 1e-7;
    
    % 通过数值微分计算磁场分量
    [dE_dx, dE_dy] = gradient_Efield(E_field, model.x, model.y);
    
    % 由法拉第定律：H_x = (1/(iωμ)) * ∂E_y/∂z ≈ 0 (二维假设)
    % H_z = (-1/(iωμ)) * ∂E_y/∂x
    H_z = (-1/(1i * omega * mu0)) .* dE_dx;
    
    % 波阻抗 Z_xy = E_y / H_z
    Z_xy = E_field ./ H_z;
    
    % 视电阻率 ρ_xy = |Z_xy|^2 / (ωμ)
    apparent_resistivity = abs(Z_xy).^2 / (omega * mu0);
end

function apparent_resistivity = calculate_ApparentResistivity_TM(H_field, model, omega)
% 计算TM模式视电阻率
    mu0 = 4 * pi * 1e-7;
    
    % 通过数值微分计算电场分量
    [dH_dx, dH_dy] = gradient_Hfield(H_field, model.x, model.y);
    
    % 由安培定律：E_x = (1/σ) * ∂H_y/∂z ≈ 0 (二维假设)  
    % E_z = (1/σ) * ∂H_y/∂x
    sigma = model.sigma;
    E_z = (1./sigma) .* dH_dx;
    
    % 波阻抗 Z_yx = E_z / H_y
    Z_yx = E_z ./ H_field;
    
    % 视电阻率 ρ_yx = |Z_yx|^2 / (ωμ)
    apparent_resistivity = abs(Z_yx).^2 / (omega * mu0);
end

参考代码 实现大地电磁TE和TM模式的正演计算 www.youwenfan.com/contentcnm/81343.html

模型验证与实例分析

模型验证

为了验证正演程序的正确性，可以设计经典地电模型进行测试：

function test_MT2D_Forward()
% 测试二维大地电磁正演程序

    % 创建均匀半空间模型
    model_resistivity = 100; % 欧姆·米
    frequencies = logspace(-3, 3, 20); % 10^-3 到 10^3 Hz
    
    % 创建模型结构
    model.name = 'UniformHalfSpace';
    model.x = linspace(-5000, 5000, 50); % 水平范围
    model.y = linspace(0, 5000, 30); % 垂直范围（深度）
    model.sigma = ones(length(model.y), length(model.x)) / model_resistivity;
    model.mesh_type = 'rectangular';
    
    % 运行正演计算
    [rho_xy, rho_yx] = MT2D_Forward(model, frequencies);
    
    % 与解析解比较
    theoretical_rho = model_resistivity * ones(size(frequencies));
    
    % 绘制结果
    figure;
    loglog(frequencies, rho_xy(25, :), 'ro', 'DisplayName', 'TE模式数值解');
    hold on;
    loglog(frequencies, rho_yx(25, :), 'bs', 'DisplayName', 'TM模式数值解'); 
    loglog(frequencies, theoretical_rho, 'k-', 'DisplayName', '理论解');
    xlabel('频率 (Hz)');
    ylabel('视电阻率 (\Omega\cdot m)');
    legend('show');
    title('均匀半空间模型验证');
end

典型地质构造模拟

你可以使用上述框架模拟典型地质构造，如地垒、地堑、断层等。这些模型的正演结果能够揭示TE和TM模式在不同地质条件下的响应规律，为实际资料解释提供依据。