前阵子整理一套野外电阻率数据甲方想在两孔之间找低阻含水破碎带。地面测线扫了整整一天反演出来的孔间深层区域几乎是一条平滑过渡带异常体轮廓完全糊掉。后来把电极放进两个钻孔做跨井采集同样的目标体一下子就立体了。这个差别本质上就是灵敏度分布sensitivity distribution在起作用。这篇文章把基于Matlab计算电阻率层析成像ERT表面和跨井XBH电极配置的2D和3D灵敏度分布的思路完整整理一遍内容包括灵敏度公式的来源、伴随法为什么高效、网格和边界怎么处理、结果图怎么读以及我实际跑代码时踩过的坑。标题虽然带着“电磁”标签但方法本身属于直流电阻率法ERT“电磁”在地电勘探板块通常是个大类泛称。适合正在做ERT反演、设计电极排列、或者刚接触地球物理数值模拟的同学参考。1. 灵敏度分布到底刻画了什么反演和测线设计的地基1.1 从反演不适定性说起为什么不能只看数据ERT反演的任务是把地表或井中测到的电位差/视电阻率转换成地下电阻率的分布。一个测区往往有几千个测量数据而地下网格单元动辄上万甚至几万个方程数量远少于未知数。这种情况下反演结果天然是非唯一的必须靠正则化、平滑约束和先验信息来“压住”不合理的波动。那问题来了哪些位置的电阻率变化对哪些测量值影响大哪些区域实际上处于“数据盲区”这就是灵敏度J_ij ∂V_i / ∂ρ_j要回答的问题。一句话理解灵敏度第j个网格单元电阻率变化一个单位第i个测量值会跟着变化多少。灵敏度高的区域测量数据对这个位置的电阻率“看得清”灵敏度接近零的区域反演基本靠相邻单元插值或者先验模型撑着数据本身提供不了约束。把每个单元的这个数值画出来得到的图就是测量系统的“视野地图”。这张图比你想象的重要。设计测线时扫一眼灵敏度剖面能提前判断电极间距、钻孔布置是否覆盖目标体反演时给灵敏度高的区域和低区域设置不同的阻尼或加权也能明显改善收敛。换句话说灵敏度计算虽不直接产出电阻率图却是反演和测线设计共同的地基。1.2 灵敏度图是“视野地图”2D与3D看的不是同一张图同一个电极配置在2D和3D模型下算出来的灵敏度形态并不一样。2D模型默认地下结构在垂直剖面外的方向y方向无限延伸电流沿y方向也当作无限长线源来处理因此灵敏度是沿y方向积分后的等效平面分布。这种模型计算快、普及率高传统ERT反演软件里的灵敏度分析多数基于2D。3D模型则把点源电流的真实扩散效应完整算出来能看到横向边界、聚焦效应、边缘衰减灵敏度图在空间上的细节更接近物理实际。代价是计算量和内存明显上升网格也从几千个单元涨到几万甚至几十万。标题里强调“2D和3D”都要算正是因为这个原因——野外先跑二维快速判断再对关键剖面或三维目标体做三维灵敏度验证是性价比很高的流程。2. 表面测线与跨井XBH测线灵敏度形态为什么会完全不同2.1 表面测量的固有短板灵敏度随深度快速衰减表面布置的电极都在z0的地表电流从电极注入后向地下半空间扩散电流密度随深度的增加迅速减小。电场强度在深部变弱测量结果对深部电阻率的敏感程度自然也就降低。典型的Wenner配置高灵敏度区域主要集中在测量排列中心下方、深度与电极间距a同量级的范围内超过大约0.5a到1a的深度后灵敏度往往只剩浅层峰值的几分之一甚至更低。这意味着表面ERT对浅层结构成像很有效对深部目标体尤其是埋深大于电极排列长度一半的目标体数据和模型之间的联系非常弱。反演时如果想靠表面数据约束150米以下的目标基本是自欺欺人——要么加大电极间距让电流穿透更深要么就得换观测方式。加大电极间距确实能增加探测深度但横向分辨率同步下降浅层灵敏度也被架空这个矛盾是表面测量的固有短板。2.2 XBH跨井测量把“视线”搬到井间跨井XBHCross-Borehole测量把电极放到两个钻孔中电流在两个井之间建立起横向通路。此时电流不再只能从地表垂直向下传播而是可以直接从一井的电极横向流到另一井井间区域的电场强度显著提升。从灵敏度角度看两孔之间的区域会形成一条相对均匀的“高灵敏度通道”水平方向上的分辨能力远强于表面测量。我自己的经验是同样的低阻破碎带埋深在井深中部位置时XBH配置的灵敏度通常比表面Wenner配置高一个数量级以上。代价是跨井配置的盲区也明显钻孔外侧、钻孔底部以下、地表浅薄层往往不在灵敏度高值区内。尤其是钻孔上部和下部超出电极覆盖范围的深度基本没有数据约束。因此跨井布孔时电极的深度范围要尽可能覆盖目标体而不是只看孔间水平距离。2.3 两种配置的灵敏度特征对比对比项表面配置Surface跨井配置XBH高灵敏度区位置地表浅层、测线中心下方两钻孔之间的井间区域主要探测深度与电极间距a直接相关由钻孔深度与电极覆盖范围决定横向分辨率中浅层较好深层快速退化井间整体较好尤其对垂直走向目标典型盲区深层、测线两侧边缘钻孔外侧、钻孔顶底之外计算规模网格纵向延伸即可规模小需要同时覆盖双孔及之间区域规模中等适用场景浅层调查、区域普查、二维剖面孔间精细探查、含水层追踪、灌浆检测这张表是我做项目时经常翻的对照表。测线设计前先拿表里的规律大致判断一下再跑灵敏度图确认一般不会出现“测完才发现覆盖不到目标体”的尴尬。3. 伴随法公式 正演求解Matlab里灵敏度计算的核心引擎3.1 先有正演后有灵敏度稳定电流场的数值求解灵敏度不是凭空画出来的它依赖地下背景电阻率电导率σ分布下的电场。ERT正演对应的是稳定电流场方程 [ abla\cdot(\sigma abla V) -I\delta(\mathbf{r}-\mathbf{r}_s) ] 其中V是电位I是点电源强度r_s是电流电极位置。给定电导率分布和边界条件求出整个空间的电位V再对电极位置采样得到测量电位这就是一次正演。边界条件处理很关键。地表当作绝缘边界∂V/∂n 0侧面和底面不能简单设V0因为截断边界会让电流在边界附近异常聚集产生边界伪影。标准做法是Robin混合边界模拟无穷远的衰减条件让边界处的电位梯度与电位本身满足一个近似比例关系这样有限网格就能近似无限半空间。3.2 伴随法的直觉一次正演换所有视角有了正演结果灵敏度怎么算最直接的想法是扰动法把第j个单元的电阻率调高1%重新正演一次看测量值变化多少这个数值就是灵敏度的一列。但网格单元动辄上万个每个单元都扰动一次等于跑一万次正演在Matlab里根本跑不动。伴随法adjoint method能把这个成本压到几乎最低。对某一测量配置A、B电流电极M、N电位电极灵敏度核函数满足一个很漂亮的对称关系用A、B注入电流时在M、N测到的电位对地下单元j的导数等于“A、B注入电流时的全场电位”和“M、N注入电流时的全场电位”这两个电场的梯度点积在单元j上的积分前面再加一个负号 [ \frac{\partial V_{MN}^{(AB)}}{\partial \sigma_j} -\int_{\Omega_j} abla U_{AB}\cdot abla U_{MN},dV ] 这里U_AB是AB供电时的合成电位场U_MN是MN供电时的合成电位场。MN本来是测量电极但在数学上把它当作虚拟电流源再解一次正演就能得到所有单元对这一次测量的灵敏度而不需要逐个单元扰动。实际代码里还有一个更省的做法每个电极位置只正演一次记录该电极单独供电时的电位场之后任意AB-MN组合的灵敏度都从预先算好的电极电位场里梯度组合出来。一个测区如果有30个电极位置只需要解30次正演就能合成所有上千组四电极测量的灵敏度。互易定理在这里相当于免费送了一半计算量。3.3 扰动法为什么在Matlab里跑不动扰动法的代码逻辑很简单for j 1:N_cells sigma_j sigma; sigma_j(j) sigma_j(j) * 1.01; V_new fwd_solve(mesh, sigma_j, config); J(:,j) (V_new - V_old) / (0.01 * sigma(j)); end问题是N_cells通常有一万个循环里面每次都要重新组装矩阵、重新求解稀疏线性方程组。单次正演可能只要几十毫秒乘以一万次就是十几分钟甚至几个小时而且矩阵反复分解内存分配也乱。伴随法每套测量配置只需处理两个电位场的组合计算量主要花在正演求解上而正演次数又可以通过电极复用大幅缩减。两相对比伴随法几乎成了ERT灵敏度计算的标准方案。4. 2D与3D网格、电极布置和边界条件代码实现的关键细节4.1 网格和边界求解结果的“底座”决定灵敏度精度网格剖分决定了正演电场的分辨率。如果网格太粗电位梯度被平均掉灵敏度图的细节就丢了网格在电极附近太粗还会让点源附近的强电场产生数值畸变。我习惯的做法是非均匀网格电极附近、目标体区域用最小间距向外按1.15到1.3的系数渐次放大远处网格尺寸可以是目标区的好几倍。这样既能保住目标区精度又不会让网格总数膨胀到难以求解。边界条件要特别检查。地表设成绝缘边界是天然正确的但侧面和底面如果用Dirichlet条件强制电位为0电流在边界处会被“憋”出强烈的伪反射灵敏度图在模型边缘会出现很奇怪的团块。换成Robin混合边界后边界伪影基本消失。我判断边界合适不合适的方法很简单算一个均匀半空间模型的灵敏度看边缘区域的值是否平滑过渡如果边缘出现条带或螺旋状高值十有八九是边界条件的问题。4.2 电极位置与测量配置的建模方式把电极位置写进代码时要先统一坐标约定。表面电极一般放在z0处x等间距排列跨井电极则分布在两个钻孔的垂线上每个钻孔有自己固定的x坐标z从孔口向下递增。% 表面测线24个电极间距5m x_surf 0:5:115; z_surf zeros(size(x_surf)); % 跨井XBH左孔x10右孔x30各15个电极间距2m x_bh1 10 * ones(15,1); z_bh1 (0:-2:-28); x_bh2 30 * ones(15,1); z_bh2 (0:-2:-28);测量配置用一个矩阵记录每一列对应一组四电极[A,B,M,N]的电极序号。跨井测量中既包括“左孔供电、右孔测量”的井间组合也包括同孔内的供电-测量组合。这两种组合的本质区别在灵敏度形态上非常明显建模时只需按电极序号组合即可代码层面不需要区分。电极落在网格节点之外时要做点源分配。最省事的做法是让网格生成时把电极坐标强制对齐到最近的节点更精细的做法是把电流按距离权重分配到相邻几个节点。实际误差通常在可接受范围但要提醒一句电极位置离网格节点太远会引入位置误差做3D精细模拟时尽量用第二种分配方式。4.3 2D和3D源码模块的代码差异2D和3D的正演模块差异是源码里最大的一块。2D网格是Nx×Nz的矩形网格未知量数量级在几千到几万Matlab里直接用稀疏矩阵加反斜杠求解很快。3D网格是Nx×Ny×Nz的三维体网格未知量数量级通常在几万到几十万内存和求解时间明显上升。一个关键优化是矩阵分解复用。无论2D还是3D正演时网格结构不变背景电导率不变变化的只是右侧的电流源项。因此可以把系数矩阵A做一次分解之后所有电极的供电问题都复用这个分解结果A assemble_stiffness_matrix(mesh, sigma); L decomposition(A); % 或 [L,U,p]lu(A) U L \ rhs; % rhs是多列矩阵每列对应一个电极这样30个电极的右侧向量可以一次回代解出比循环30次重新求解快一个数量级。3D网格的稀疏矩阵规模虽大但对称正定的刚度矩阵配合Cholesky分解或分解对象在R2020b之后的Matlab版本里表现很稳定。2D和3D还有一个本质区别2D严格来说要做y方向的积分或2.5D波数域变换才能对应真实点源如果源码直接把点源当成线源做二维正演灵敏度剖面反映的是“无限长异常体”的等效响应这个近似在解读结果时心里要有数。5. 运行源码的完整流程与三张结果图怎么读5.1 从参数表到主脚本三分钟跑通第一个算例拿到源码包后先别急着改代码。我建议按这个顺序过一遍第一步打开主脚本找到参数设置区第二步把背景电阻率、网格范围、电极坐标填好第三步调用正演和灵敏度两个子函数第四步跑绘图部分。典型的源码包通常会包含这几个模块网格生成脚本根据测区范围、电极位置自动生成非均匀网格正演求解脚本组装刚度矩阵、加载荷、解方程返回电极位置电位灵敏度计算脚本调用电位场结果对每个配置组合计算灵敏度核绘图脚本画2D剖面、3D切片、累计灵敏度图。我常用的一组测试参数是背景电阻率100欧姆米表面Wenner测线24个电极、间距5米两个钻孔相距20米、每孔15个电极、间距2米。参数填好后先跑一个表面2D灵敏度验证流程通顺再依次跑跨井XBH和3D。第一次跑通建议关闭所有并行池避免并行参数干扰问题定位。5.2 单条灵敏度剖面、累计覆盖图和3D切片怎么读跑完之后会得到三类图读法各有侧重。第一张是单条测量配置的2D灵敏度剖面。单个配置的灵敏度往往有正有负高值区集中在电流路径和测量路径的公共“通视区域”。这张图的意义是理解“某一个测量对到底在看哪里”对研究数据权重和剔除坏数据很有帮助。颜色条建议用对称色标保留正负信息。第二张是累计灵敏度图。把所有配置的灵敏度绝对值叠加起来常用Σ|J_i|或RMS合成得到整个测量系统对地下空间的“照明强度”。表面测线累计图一般是浅层亮、深层暗呈一个透镜状的高值区跨井XBH累计图会在两个钻孔之间的区域形成亮带这个亮带的上下边界基本就是电极覆盖范围的投影。这张图是测线设计评审时最常用的依据。第三张是3D灵敏度切片。通常用slice命令取xy、xz、yz三个正交切面或者用isosurface取某个阈值对应的三维包络。3D切片能看出电流在三维空间里的聚焦效应尤其是孔间区域的三维轮廓比2D图直观很多。由于灵敏度动态范围大画3D前建议先取对数log10(|S_cum|)或log10(1|S_cum|)否则颜色条会被浅层的大值拉爆。5.3 调试技巧和解析解对一遍再画图运行结果异常时最快定位问题的方法是利用均匀半空间模型下的解析灵敏度做基准。比如表面四极法在均匀半空间中有解析的灵敏度核函数形式先用少量网格算数值灵敏度再和解析值对比看误差是否随网格加密而收敛。如果趋势不一致问题多半出在边界条件或梯度计算上如果趋势一致但数值差常数多半是系数约定或符号约定不一致需要回头检查灵敏度公式里的负号和电导率/电阻率参数化方式。这个“先对解析解再上复杂模型”的习惯帮我省了大量排查时间。不要相信任何一张没经过验证的灵敏度图——它看起来漂亮可能只是边界条件假出来的。6. 灵敏度计算的实战经验与易被忽视的陷阱6.1 负灵敏度带来的两个坑抵消和伪振荡灵敏度的正负号在单条配置图里很常见尤其是偶极-偶极、跨井交叉组合。负灵敏度意味着该处电阻率升高时测量电位反而下降。直接把这些值累加会互相抵消累计灵敏度图出现一片灰暗误导“这里没有数据约束”的判断。解决方法是累计时取绝对值或者RMS并且单独保留正负输出备查。另外反演时如果初始模型离真实模型太远高负灵敏区的线性响应会导致迭代方向出现伪振荡需要在阻尼因子里给予考虑。6.2 电极邻域效应与绘图颜色条调节点源电极附近的电流密度理论上趋于无穷大灵敏度数值也会疯涨几个量级。如果直接画图颜色条会被电极附近几个网格完全占据其余区域的分布细节全被压平。我的处理办法是绘图时把电极周围半径等于2~3倍最小网格步长的区域屏蔽掉或者对灵敏度做阈值截断。这在代码里就一行S_plot S_cum; S_plot(dist_to_electrode 1.5 * h_min) NaN;屏蔽半径不能太大否则会把电极间原有高灵敏区也抹掉一般取最小网格尺寸的2倍左右比较平衡。6.3 一个绕不开的模型依赖背景电阻率会影响灵敏度灵敏度不是只由电极几何决定的“通用分布”它依赖背景电导率场。均匀半空间下算出来的灵敏度图和存在高阻屏蔽层或低阻通道时的灵敏度图会有可分辨的差异。低阻体像一个“电流汇”会把电流吸引过去导致其后方区域灵敏度下降高阻体则像屏蔽墙会让电流绕行绕行路径上的灵敏度反而增强。所以灵敏度分析要用接近真实背景的模型来做尤其在做跨井探测时两孔之间的地层差异会直接影响XBH的成像能力。这一点常常被忽略以为灵敏度只跟电极位置有关。6.4 用数值扰动验证灵敏度代码的正确性最后给一个通用的验证方法。在模型里挑三五个有代表性的单元浅层、深层、孔间、边界附近把该单元电阻率改变一个小量δ重新正演得到测量电位差变化ΔV然后和代码输出的灵敏度数值比较 [ J_{num} \frac{V(\rho \delta\rho_j) - V(\rho - \delta\rho_j)}{2\delta\rho_j} ] 数值差分解与伴随法结果一致代码就可以放心用。这个验证脚本不要删掉后续改网格、改边界条件、扩展3D时都跑一遍几秒钟的代价能省下几天排查错误的时间。我在实际项目里体会最深的一件事是真正影响成果质量的因素往往不是反演算法原本的迭代公式而是测线设计阶段有没有把灵敏度分布想清楚。表面测线觉得深部“反演不出来”先别急着加正则化或者调权重先看一眼灵敏度图数据是不是真的穿到了那个深度跨井项目里两孔之间成像糊排查顺序也应该从灵敏度覆盖开始。这套Matlab源码把表面和XBH配置的2D、3D灵敏度计算流程打通之后我几乎每个项目开工前都要跑一轮灵敏度图相当于在野外数据采集之前先给测量系统做一次“视力表检查”。如果你要在这个基础上继续扩展可以考虑把计算出的灵敏度作为反演Jacobi矩阵的直接输入或者把它做成测量配置优化的目标函数——这一步的收益通常比调平滑系数大得多。