做一维频域电磁反演时雅可比矩阵——也就是灵敏度矩阵——往往是第一个让新手摔跟头的地方。原因并不复杂正演告诉你给定的层状模型能算出什么样的视电阻率和相位曲线但反演需要回答的是“每一个模型参数往哪个方向动、动多少曲线才会更接近实测点”。去年我在处理一组实际频域测深数据时就踩过这个坑一开始图省事用有限差分去估算灵敏度结果在视电阻率曲线极小值附近怎么迭代都不收敛。后来花了一晚上把解析灵敏度的数学链推清楚同样的数据三到五次迭代就进去了。这篇文章就以大地电磁MT这种最典型的频域电磁方法为例把一维层状模型解析灵敏度矩阵的计算过程完整拆解一遍并给出可以直接拿去跑的Matlab代码。文章适合两类人一类是刚入门电磁反演、想搞清楚灵敏度到底是什么的研究生另一类是手里已有正演程序、想把它扩展成反演程序的工程师。你不需要提前精通反演理论只要熟悉一维正演的递归公式后面每一步都是链式法则的机械应用——但机械也容易漏项所以我会把每个偏导数的来源都标注清楚。1. 解析灵敏度不是炫技是反演迭代的必选项1.1 灵敏度矩阵在频域EM反演里的角色反演问题通常写成最小化目标函数的形式Φ(m) ||d_obs - F(m)||² λ||m - m_prior||²其中 d_obs 是实测数据向量F(m) 是正演算子m 是模型参数向量λ 是正则化权重。在Gauss-Newton这类迭代解法中模型更新方向由一阶导数决定也就是雅可比矩阵 J ∂F/∂m。对于一维频域电磁数据最常见的观测项是若干个频点上的视电阻率 ρa 和阻抗相位 φ模型参数则是每一层的电阻率 ρj 和厚度 hj。J 矩阵中第 i 行第 j 列的元素就是“第 i 条数据对第 j 个模型参数的偏导数”。这个矩阵的质量直接决定反演每一步能不能朝着正确的方向走。很多人觉得灵敏度只是反演的“副产品”其实它本身就是信息量最大的中间产物。通过 J你可以知道某一条视电阻率曲线主要受浅层还是深层电阻率控制、哪些参数在某个频段基本不可见、哪些参数之间存在严重耦合。这些信息在做测深方案设计时同样有用。所以我一直建议即使你不做反演把灵敏度矩阵算出来看看也大有裨益。1.2 数值差分为什么在这种场景下不划算教科书里最常用的灵敏度计算方法是中心差分J_ij ≈ [F_i(m δe_j) - F_i(m - δe_j)] / (2δ)乍看很简单但实际使用时会遇到几个麻烦。首先是步长 δ 的选择。电磁响应对参数并不是线性变化特别是在视电阻率曲线拐点附近δ 取太大会把高阶项误差引入δ 取太小又会被浮点舍入误差主导。我在两层模型里测试过δ 取 10⁻³ 和 10⁻⁶ 得到的灵敏度可以差出几个百分点而没有先验信息时你很难判断哪个值是对的。其次是计算成本。虽然一维正演本身不贵但若干频率乘以若干层数每个参数都要额外调正演两次。如果模型有 20 层、频率有 40 个点一次完整灵敏度计算就需要 40×20×2 1600 次正演。这还没算反演要迭代几十轮。相比之下解析灵敏度是在正演递推过程中“搭车”算出来的几乎不需要额外正演次数。更重要的是精度稳定性问题。相位数据在接近 90° 时变化往往很平缓数值差分得到的微小数值会被多次正演的误差淹没。而解析公式直接给出连续可微的偏导数在同一个频率和模型点上是精确的。所以我的观点非常直接一维场景下解析灵敏度不是“高级选项”而是反演工程中的必选项。1.3 解析公式完成后数值差分依然有大用途我并不是说有了解析公式就要彻底抛弃数值差分。恰恰相反解析推导过程中只要某一个符号传错得到的灵敏度就可能“看起来合理但实际错误”——例如所有元素都反向、公因子差了某个复数因子。所以我现在的固定流程是先把解析代码写出来再用中心差分做一次全参数、全频点对比确认最大相对误差小于 10⁻⁵之后才放心把灵敏度矩阵投入反演循环。这篇文章后面会给出验证脚本的思路你完全可以沿用这个流程。2. 一维MT正演从亥姆霍兹方程到阻抗递推式2.1 模型、符号和单位的第一道坎假设地下模型由 N 层水平均匀介质组成第 j 层电阻率为 ρj单位 Ω·m厚度为 hj单位 m第 N 层为半空间即 hj ∞。地表的天然电磁波可以近似为垂直入射的平面波主要考虑 TE 模式的电场和磁场分量。取直角坐标系x-y 为水平面z 轴向下。频率为 f Hz 的平面波在地下传播定义角频率 ω 2πf磁导率 μ 4π×10⁻⁷ H/m绝大多数陆地岩石的磁导率与真空几乎一致可以直接用真空值电导率 σj 1/ρj。单位混乱是我见过最多的问题。频率到底用周波还是角频率厚度用 km 还是 m电阻率用 Ω·m 还是 S/m 转换我的建议非常简单全部用国际单位制f 用 Hzρ 用 Ω·mh 用 m最后得到的视电阻率天然是 Ω·m。代码里我会把 ω 显式写成 2pifreq避免中间换算出错。2.2 层状介质阻抗递归式的推导脉络在频率域中TE 极化波在均匀层内的电场和磁场满足亥姆霍兹方程解可以表示为向下和向上传播的指数函数叠加也就是双曲函数 tanh 的形式。定义第 j 层内的复波数k_j sqrt(i ω μ σ_j)这里 i 是虚数单位需要取实部为正的平方根根式保证场随深度衰减而不是增长。该层对应的本征阻抗为Z0_j i ω μ / k_j它表示假设该层无限厚时地表阻抗应当是这个值。对于厚度为 hj 的实际层顶面阻抗 Zj 与下一层顶面阻抗 Z_next 的关系为Zj Z0_j · (Z_next Z0_j · tanh(k_j h_j)) / (Z0_j Z_next · tanh(k_j h_j))递推从最底层半空间开始。半空间没有下一层所以 Z_N Z0_N。然后从 j N-1 一直向上计算到 j 1最终得到地表阻抗 Z1。这个公式是大多数一维和二维电磁正演程序的“心脏”。代码实现时需要注意 tanh 的参数是复数Matlab 中直接调用 tanh 即可。物理上当层厚度远大于趋肤深度时tanh(kh) 趋近于 1此时上层对中层与深层之间的耦合几乎没有影响递推公式自动退化到均匀半空间的情况。2.3 视电阻率和阻抗相位的换算关系地表阻抗 Z1 本身不是最终观测数据大地电磁通常记录视电阻率和相位ρa |Z1|² / (ω μ) φ arg(Z1)注意这里的 arg 返回弧度值如果想得到度数乘以 180/π。均匀半空间时ρa ρφ 45°这是一个非常容易记忆的自检标准。相位定义在不同教材里可能有差异常见的是 E/H 比值的辐角。如果代码输出的均匀半空间相位是 -45°说明 k 或 Z0 的符号约定与你所在课题组的标准不同这不是错误只要把符号约定统一即可。但如果你要用现成反演代码必须确认它期望的相位约定是什么。2.4 为什么我说正演是解析灵敏度最好的“载体”从递推式可以看出计算 Z1 的过程中每一层的 k_j、Z0_j、t_j tanh(k_j h_j) 和分母 Den_j 都会被用到。这些中间量在标准正演完成后其实已经全部存在如果你另写一个灵敏度程序又得重新算一遍或者把正演结果存下来。最优雅的做法是让灵敏度计算直接嵌入正演的递推循环里每更新一层就把这一层的导数关系也更新一次。这样不仅省掉重复计算而且天然保证了正演和灵敏度的一致性——同一个函数里算出来的 Z 和 dZ 不会因为中间量定义不同而错位。3. 解析灵敏度推导在递归式上做链式法则3.1 每一步需要维护的量dZ/dm把模型参数按 ρ1, h1, ρ2, h2, ..., ρN-1, hN-1, ρN 的顺序排成一个长向量 m长度 npar 2N-1。递推的过程中我们不仅要知道当前层的顶面阻抗 Zj还要知道它对所有模型参数的偏导行向量 dZj/dm这是一个复数向量。从最底层开始Z_N Z0_N因此 dZ_N/dm 只有一个非零分量对 ρ_N 的导数 dZ_N/dρ_N Z0_N/(2ρ_N)其余全为 0。然后逐层向上更新。关键思路是第 j 层顶面阻抗 Zj 是 Z_next、Z0_j、t_j 这三者的函数而 Z_next 本身又依赖更深层的参数。因此链式法则可以拆成三部分来自 Z0_j 对当前层电阻率的导数、来自 t_j 对当前层电阻率和厚度的导数、以及来自 Z_next 对更深层参数的导数。这一步就是整个灵敏度矩阵计算的核心逻辑。剩下的工作只是把每个偏导数写成显式表达式。3.2 三个核心偏导数 A/B/C 的由来为了表述简洁把阻抗递推式记作Z f(Z0, Z_next, t)其中 t tanh(k h)Den Z0 Z_next·t对任意实数参数 p 求导得到dZ/dp (∂f/∂Z0)(dZ0/dp) (∂f/∂Z_next)(dZ_next/dp) (∂f/∂t)(dt/dp)三个偏导数中∂f/∂Z_next 表示下一层顶面阻抗的变化如何传导到当前层∂f/∂Z0 和 ∂f/∂t 表示本层电性参数变化的影响。经过直接求导可以得到A ∂f/∂Z_next Z0²(1-t²) / Den²B ∂f/∂t Z0(Z0² - Z_next²) / Den²C ∂f/∂Z0 t(Z_next² Z0² 2 Z0 Z_next t) / Den²这三个式子看起来很吓人但它们只是对两个复数的商做商规则偏导展开后每一项都能对上。A 和 C 无量纲B 的量纲与阻抗相同。代码中直接使用这三个表达式即可不需要再化简。如果担心记错可以把它们当作“灵敏度传播因子”来理解A 把深层参数的灵敏度传到上层B 和 C 则把当前层电阻率或厚度的影响注入灵敏度向量。3.3 中间量导数的具体表达式除了 A/B/C还需要几个基础导数。首先是波数对电阻率的导数。由于 k sqrt(iωμ/ρ)所以dk/dρ -k / (2ρ)其次是本征阻抗对电阻率的导数。Z0 iωμ/k代入上式可得非常简洁的形式dZ0/dρ Z0 / (2ρ)然后是双曲函数项 t tanh(kh) 的导数。利用 1-t² sech²(kh)可以写成dt/dh k(1 - t²)dt/dρ h(1 - t²) · (dk/dρ) -k h(1-t²) / (2ρ)最后半空间阻抗 Z_N Z0_N它的灵敏度只与 ρ_N 有关公式同上。这些导数在代码中都是一行就能算完的但符号和系数错一个整个矩阵就会废掉所以建议全部推导一遍再写代码。3.4 从表面阻抗到视电阻率和相位的灵敏度转换表面阻抗 Zsurf 的灵敏度向量 dZdm 是复数行向量。但我们的观测数据是实数的视电阻率和相位因此还需要把复灵敏度转换成实灵敏度。对于任意实数模型参数 p由视电阻率定义可得dρa/dp [Z · conj(dZ/dp) conj(Z) · dZ/dp] / (ωμ) 2 real(conj(Z) · dZ/dp) / (ωμ)这里 real 是取实部。相位 φ angle(Z)其导数为dφ/dp imag(conj(Z) · dZ/dp) / |Z|²如果相位以度输出再乘以 180/π。代码中我会用这种形式实现因为它比拆开实部虚部分别求导更紧凑也不容易漏项。4. Matlab实现与数值验证4.1 代码结构与参数索引设计我写的函数名为mt1d_sensitivity输入频率向量、电阻率向量和厚度向量输出视电阻率、相位以及对应的两个灵敏度矩阵。参数索引规则是ρj 在第 2j-1 个位置hj 在第 2j 个位置这样在递推内层根据 j 就能直接定位 dZ 向量中对应列无需额外映射。函数内部先做正演和灵敏度联合递推最后统一转换为视电阻率和相位灵敏度。代码对每个频率分别循环这是为了把推导体现在清晰的结构中。当模型层数不多且频率点数在几十以内时这个循环开销完全可以接受。若以后需要处理上千频点再考虑把频率循环向量化也不迟。4.2 完整代码下面是可以直接复制运行的Matlab函数。注释里保留了公式出处方便你对照本文内容检查。function [rho_a, phi_deg, drhoa_dm, dphi_dm] mt1d_sensitivity(freq, rho, h) % MT1D_SENSITIVITY 一维大地电磁解析灵敏度矩阵 % % 输入: % freq - 频率向量 (Hz)长度为 Nf可以是行向量或列向量 % rho - 各层电阻率 (Ohm m)长度 Nrho(end) 为半空间电阻率 % h - 各层厚度 (m)长度 N-1最后一层不用给 % % 输出: % rho_a - 视电阻率 (Ohm m)Nf x 1 % phi_deg - 阻抗相位 (度)Nf x 1 % drhoa_dm - 视电阻率对模型参数的偏导Nf x (2N-1) % dphi_dm - 相位对模型参数的偏导Nf x (2N-1) % % 模型参数顺序: % m [rho_1, h_1, rho_2, h_2, ..., rho_{N-1}, h_{N-1}, rho_N] mu0 4 * pi * 1e-7; omega 2 * pi * freq(:); N length(rho); Nf length(freq); npar 2 * N - 1; rho rho(:); h h(:); rho_a zeros(Nf, 1); phi_deg zeros(Nf, 1); drhoa_dm zeros(Nf, npar); dphi_dm zeros(Nf, npar); for i 1:Nf w omega(i); % ---------- 递推起点半空间阻抗及其灵敏度 ---------- sigN 1 / rho(N); kN sqrt(1i * w * mu0 * sigN); Z0N 1i * w * mu0 / kN; Z_next Z0N; dZ_next zeros(1, npar); dZ_next(2 * N - 1) Z0N / (2 * rho(N)); % dZ_N / drho_N % ---------- 从 N-1 层逐层向上递推 ---------- for j N-1 : -1 : 1 rj rho(j); hj h(j); kj sqrt(1i * w * mu0 / rj); Z0j 1i * w * mu0 / kj; t tanh(kj * hj); Den Z0j Z_next * t; Zj Z0j * (Z_next Z0j * t) / Den; sech2 1 - t * t; A Z0j^2 * sech2 / Den^2; B Z0j * (Z0j^2 - Z_next^2) / Den^2; C t * (Z0j^2 Z_next^2 2 * Z0j * Z_next * t) / Den^2; % 深层参数的灵敏度先乘上传播因子 A dZ A * dZ_next; % 当前层电阻率 rho_j 的贡献 dk_drho -kj / (2 * rj); dt_drho hj * sech2 * dk_drho; dZ0_drho Z0j / (2 * rj); idx_r 2 * j - 1; dZ(idx_r) dZ(idx_r) C * dZ0_drho B * dt_drho; % 当前层厚度 h_j 的贡献 dt_dh kj * sech2; idx_h 2 * j; dZ(idx_h) dZ(idx_h) B * dt_dh; Z_next Zj; dZ_next dZ; end Zsurf Z_next; dZdm dZ_next; % ---------- 输出视电阻率和相位 ---------- rho_a(i) abs(Zsurf)^2 / (w * mu0); phi_deg(i) angle(Zsurf) * 180 / pi; % ---------- 灵敏度转换 ---------- factor conj(Zsurf) * dZdm; drhoa_dm(i, :) 2 * real(factor) / (w * mu0); dphi_dm(i, :) imag(factor) / abs(Zsurf)^2 * 180 / pi; end4.3 均匀半空间和两层模型的手工验证拿到代码第一件事不是直接扔进反演而是做两个基本算例。第一个算例是均匀半空间。令 N1rho[100]h[]任意频率下理论上都有 ρa100 Ω·mφ45°。解析灵敏度应当满足 dρa/dρ11dφ/dρ10。运行代码后检查这个结果如果不对说明 Z0 或视电阻率定义有问题而不是灵敏度的问题。第二个算例是两层模型。比如 rho[100, 20]h[500]频率从 0.001 到 10 Hz。用中心差分对比解析灵敏度核心思路是对第 k 个参数做一个小扰动 ±δ分别调用正演得到两条视电阻率曲线差分得到近似灵敏度再与解析灵敏度比较。注意扰动幅度要取得合适我通常取 1e-6 倍的参数值既保证差分精度又不至于被舍入误差干扰。下面这段脚本是数值差分对比的核心片段你可以嵌入自己的测试环境delta rho(1) * 1e-6; [~, ~, J_par_rho1, ~] mt1d_sensitivity(freq, rho delta * [1;0], h); [~, ~, J_minus_rho1, ~] mt1d_sensitivity(freq, rho - delta * [1;0], h); J_num_rho1 (J_par_rho1 - J_minus_rho1) / (2 * delta);对所有参数重复这个过程然后计算max(abs(J_num - J_ana))。如果数值在 1e-6 量级或更小解析代码基本可以信任。注意相位灵敏度如果遇到跨象限反转要先把相位做 unwrap 再比较否则会出现看似很大的“假误差”。4.4 关于代码效率的一点经验上面代码的核心复杂度是频率数 × 层数²因为每一层都要维护长度为 2N-1 的灵敏度向量。对于 N50、频率 100 个点的情况循环次数约 25 万次复数运算Matlab 跑起来也就几秒钟完全够用。如果你追求极致性能可以把 dZ 向量改成稀疏结构因为深层参数在每一层的传播只是整体乘一个 A 因子不需要每次都全列更新。但反演中 J 矩阵最终都是稠密矩阵所以稀疏化意义不大。5. 实际使用中的坑与向其他频域EM方法的延伸5.1 单位、频率范围和厚度离散的坑灵敏度矩阵对单位极其敏感。我在代码里明确使用了 SI 单位但很多既有程序习惯用 km 和周期。如果直接把数据套进这个函数结果会非常离谱。另一个坑是频率范围选取不当导致灵敏度饱和当频率很高时趋肤深度很小最深层半空间的灵敏度会趋近于 0当频率很低时浅层灵敏度趋近线性但可能失去约束。反演前应该先画一张灵敏度随频率和深度变化的伪彩图看看你的观测频段到底覆盖了哪些深度范围。5.2 相位分支与反正切的不连续相位由 angle(Zsurf) 得到取值范围是 (-π, π]。当观测相位跨过 ±180° 时会出现常规相位序列的“跳变”。虽然作为局部导数的 dφ/dp 本身是连续的但如果你用 angle 直接输出而不做处理反演迭代时可能因为相位跳变产生虚假残差。解决办法有两种一是在正演输出时把相位归一到 [0, 90] 区间对 MT 数据通常可行二是在反演中用复阻抗实部和虚部作为数据而不是用视电阻率和相位。后者的灵敏度计算更加平滑但物理上不如视电阻率直观。5.3 深层参数与厚度参数的可分辨性问题解析灵敏度能告诉你什么参数能分辨什么不能。实际经验是高阻薄层经常“隐身”因为低阻层的灵敏度会把信号屏蔽掉厚层的厚度灵敏度往往大于电阻率灵敏度而薄层恰恰相反。如果你在反演中发现某一列灵敏度矩阵数值极小不要强行迭代——那不是数值问题而是数据本身缺乏对该参数的分辨能力。这时候调整正则化权重或者修改层参数化方式比盲目增加迭代次数更有效。5.4 向可控源和频率域电磁测深法扩展标题说的是“频域 EM 数据”MT 只是最典型的一维频域 EM 例子。对于可控源音频大地电磁CSAMT或频率域电磁测深FEM只要正演输出仍是基于层状介质阻抗递推得到的接收点响应灵敏度传播的 A/B/C 结构就完全一样。不同的只是正演表达式里多了源项积分最终的观测数据可能是电场分量或磁场分量。你只需要把改进后的正演结果代入同样的链式法则先在源项部分求导再沿着阻抗递推链传播即可。所以在学透一维 MT 解析灵敏度后再扩展其他频域方法会非常轻松。最后分享一个实用技巧在反演循环里我会把正演和灵敏度合成一个函数避免中间量重复计算。如果模型层数超过 50 层或频率点数超过 200 个建议在循环外面预分配矩阵并把 freq 循环改为向量化表达式。代码中的 A/B/C 任何一个式子和你的推导对不上最稳妥的验证方法就是先做均匀半空间和两层模型的数值差分对比这一步永远不要省。