【计算一维频域 EM 数据的解析灵敏度】频域 EM 数据解析灵敏度矩阵的计算附Matlab代码

✅作者简介：热爱科研的Matlab仿真开发者，擅长数据处理、建模仿真、程序设计、完整代码获取、论文复现及科研仿真。

🍎 往期回顾关注个人主页：Matlab科研工作室

🍊个人信条：格物致知,完整Matlab代码及仿真咨询内容私信。

🔥内容介绍

一、核心概念与定义

1.1 频域EM数据与解析灵敏度

频域电磁法（Frequency-Domain EM, FDEM）通过向地下发射不同频率的谐变电磁场，测量地表或空中的电场（E）、磁场（H）响应，反演地下介质的电导率（σ）、磁导率（μ）、介电常数（ε）等物性参数。解析灵敏度描述了**介质物性参数的微小变化对EM响应的影响程度**，是反演过程中权重分配、收敛性判断的核心依据。

对于一维（层状）介质，物性参数仅沿垂直方向（z轴）变化，EM响应可通过解析方法求解（如汉克尔变换、传输线理论），灵敏度矩阵也可通过解析推导获得，避免了数值模拟（如有限元、有限差分）的离散误差。

1.2 灵敏度矩阵的数学表达

设观测数据向量为 \( \mathbf{d} = [d_1, d_2, ..., d_m]^T \)（m为观测点-频率组合数，如不同频率下的阻抗、相位），模型参数向量为 \( \mathbf{m} = [m_1, m_2, ..., m_n]^T \)（n为模型层数对应的物性参数，如各层电导率），则解析灵敏度矩阵 \( \mathbf{J} \in \mathbb{R}^{m \times n} \) 的元素定义为：

\( J_{ij} = \frac{\partial d_i}{\partial m_j} \)

其中，\( J_{ij} \) 表示第j个模型参数的单位微小变化，引起第i个观测数据的变化量。当参数变化量 \( \Delta m_j \ll m_j \) 时，可通过微分近似数据变化：\( \Delta d_i \approx \sum_{j=1}^n J_{ij} \Delta m_j \)，即 \( \Delta \mathbf{d} \approx \mathbf{J} \Delta \mathbf{m} \)。

一维场景中，模型参数通常取各层的电导率（σₖ，k=1,2,...,n）、层厚度（hₖ，k=1,2,...,n-1），磁导率默认取真空磁导率 \( \mu_0 = 4\pi \times 10^{-7} \, \text{H/m} \)（非磁性介质），介电常数在低频段（<10⁶ Hz）影响可忽略，仅高频时需考虑。

二、一维频域EM响应的解析求解基础

计算灵敏度矩阵的前提是获得EM响应的解析表达式，一维层状介质的频域EM响应求解核心是**场分量的边界条件匹配**与**汉克尔变换（对水平波数积分）**。以下以大地电磁法（MT，属于FDEM）为例，简述响应求解逻辑（其他FDEM方法如CSAMT、TEM频域响应可类比推导）。

2.1 基本控制方程

对于谐变电磁场（时间因子 \( e^{i\omega t} \)，ω为角频率），在无源、均匀各向同性介质中，麦克斯韦方程组可简化为亥姆霍兹方程：

\( \nabla^2 \mathbf{E} - i\omega \mu \sigma \mathbf{E} = 0 \)

\( \nabla^2 \mathbf{H} - i\omega \mu \sigma \mathbf{H} = 0 \)

定义复波数 \( k = \sqrt{i\omega \mu \sigma} = (1+i)\sqrt{\pi \omega \mu \sigma} \)，表征电磁场在介质中的衰减与传播特性。

2.2 一维层状介质的场分解与匹配

一维介质中，电磁场可分解为TE极化（电场垂直于测线，E⊥x轴）和TM极化（磁场垂直于测线，H⊥x轴），两类极化的场分量独立求解，边界条件为：

电场切向分量连续：\( E_{x1} = E_{x2} \)（TE极化）、\( E_{z1} = E_{z2} \)（TM极化）；
磁场切向分量连续：\( H_{x1} = H_{x2} \)（TM极化）、\( H_{z1} = H_{z2} \)（TE极化）；
法向分量满足电流连续性：\( \sigma_1 E_{z1} = \sigma_2 E_{z2} \)（TE极化）、\( \mu_1 H_{z1} = \mu_2 H_{z2} \)（TM极化）。

通过汉克尔变换将水平方向（x-y平面）的场分量转化为水平波数域（kₚ域），每层的场分量可表示为上行波与下行波的叠加，结合边界条件递归求解各层的波幅系数，最终通过汉克尔逆变换得到空间域的EM响应（如阻抗 \( Z = E_x / H_y \)、相位 \( \phi = \arg(Z) \)）。

三、解析灵敏度矩阵的计算步骤

基于EM响应的解析表达式，灵敏度矩阵元素 \( J_{ij} = \partial d_i / \partial m_j \) 的计算分为“直接微分法”和“伴随场法”，一维场景中直接微分法更简洁高效，具体步骤如下：

3.1 步骤1：确定观测数据与模型参数

观测数据 \( d_i \)：选取可解析求解的EM响应量，如大地电磁阻抗幅值 \( |Z| \)、相位 \( \phi \)，可控源音频大地电磁法（CSAMT）的电场幅值 \( |E| \)、磁场幅值 \( |H| \) 等；
模型参数 \( m_j \)：核心参数为各层电导率 \( \sigma_k \)（j对应k时，\( m_j = \sigma_k \)）、层厚度 \( h_k \)（j对应n+k-1时，\( m_j = h_k \)），忽略磁导率和介电常数时，n为层数（厚度参数为n-1个）。

3.2 步骤2：推导EM响应对模型参数的微分表达式

以MT阻抗 \( Z \) 对第k层电导率 \( \sigma_k \) 的灵敏度为例，推导核心逻辑如下：

由EM响应解析式，将 \( Z \) 表示为各层波幅系数 \( A_k、B_k \)（依赖 \( \sigma_k、h_k、k \)）的函数：\( Z = f(A_1,B_1,...,A_n,B_n) \)；
对 \( \sigma_k \) 求偏导，应用链式法则：\( \frac{\partial Z}{\partial \sigma_k} = \sum_{p=1}^n \left( \frac{\partial Z}{\partial A_p} \cdot \frac{\partial A_p}{\partial \sigma_k} + \frac{\partial Z}{\partial B_p} \cdot \frac{\partial B_p}{\partial \sigma_k} \right) \)；
波幅系数 \( A_p、B_p \) 通过边界条件递归获得，其对 \( \sigma_k \) 的微分需结合每层的复波数 \( k = \sqrt{i\omega \mu \sigma} \) 的微分（\( \partial k / \partial \sigma = k/(2\sigma) \)），以及边界条件方程组的隐函数微分求解。

对于层厚度 \( h_k \) 的灵敏度，核心差异在于波幅系数对 \( h_k \) 的微分（\( \partial A_p / \partial h_k \)、\( \partial B_p / \partial h_k \)）与层厚度引起的相位延迟相关，推导时需保留波数与厚度的乘积项（\( k \cdot h \)）的微分。

3.3 步骤3：汉克尔变换域与空间域的转换

EM响应的解析表达式通常在汉克尔变换域（kₚ域）推导，灵敏度的微分也需在变换域完成后，通过汉克尔逆变换转换至空间域，得到观测点处的灵敏度值。

汉克尔逆变换公式（以零阶为例）：\( f(r) = \int_0^\infty k_p J_0(k_p r) F(k_p) dk_p \)，其中 \( J_0 \) 为零阶贝塞尔函数，\( F(k_p) \) 为变换域函数。对灵敏度而言，需满足：

\( \frac{\partial f(r)}{\partial m_j} = \int_0^\infty k_p J_0(k_p r) \frac{\partial F(k_p)}{\partial m_j} dk_p \)

实操中，汉克尔积分通过数值积分（如自适应高斯积分、滤波法）求解，需注意积分区间的截断误差（通常取 \( k_p \in [10^{-4}, 10^4] \, \text{rad/m} \) 覆盖主要波数成分）。

3.4 步骤4：灵敏度矩阵的归一化处理

由于模型参数（如σ：S/m，h：m）与观测数据（如Z：Ω，φ：rad）的量纲差异，直接计算的灵敏度元素量级差异极大，需进行归一化处理，常用方法包括：

模型归一化：\( J'_{ij} = J_{ij} \cdot \frac{m_j}{d_i} \)，表征参数相对变化引起的数据相对变化，适合对数反演；
数据归一化：\( J'_{ij} = J_{ij} / d_i \)，消除数据量级影响；
标准化：\( J'_{ij} = \frac{J_{ij} - \bar{J}_j}{\sigma_{J_j}} \)，按列（同一参数对不同数据的灵敏度）标准化，使各列灵敏度量级一致。

四、关键难点与应对策略

4.1 高维积分的数值稳定性

汉克尔积分的收敛性直接影响灵敏度精度，尤其是低频段（小kₚ）和高频段（大kₚ）的积分贡献。应对策略：采用**滤波法汉克尔变换**（如数字滤波系数逼近贝塞尔函数积分），或自适应积分算法，结合EM响应的波数域衰减特性（大kₚ时场衰减快，可缩短积分上限）。

4.2 边界条件微分的复杂性

多层介质中，某一层参数变化会通过边界条件连锁影响所有上层的波幅系数，导致微分表达式冗长。应对策略：采用**递归法推导**，从底层（半无限介质，仅上行波）向上逐层递推波幅系数的微分，避免直接求解高维方程组。

4.3 高频段介电常数的影响

当频率高于10⁶ Hz时，介电常数ε的贡献不可忽略，复波数需修正为 \( k = \sqrt{i\omega \mu (\sigma + i\omega \varepsilon)} \)，灵敏度需额外考虑对ε的微分。应对策略：将ε纳入模型参数，按电导率灵敏度的推导逻辑扩展，或通过复介电常数 \( \varepsilon^* = \varepsilon + i\sigma/(\omega) \) 统一表征。