用Matlab把电阻率层析成像ERT的灵敏度分布算清楚这件事看着偏理论却是决定反演结果可信度的关键一步。最近我把表面电极和跨井电极cross-boreholeXBH配置下的2D/3D灵敏度分布完整跑了一遍从电极坐标、网格剖分到正演求解再到灵敏度矩阵组装和可视化算是把整条链路打通了。这篇内容就是一次实操复盘重点讲清楚灵敏度分布到底在算什么、Matlab里怎么搭这套计算流程、以及表面ERT和XBH在灵敏度覆盖上的本质差异刚接触ERT数据模拟或者做反演前设计的朋友可以直接参考。1. 认识ERT、跨井测量与灵敏度分布1.1 电阻率层析成像到底在测什么电阻率层析成像Electrical Resistivity TomographyERT是一种以地下岩土体电阻率差异为成像目标的地球物理方法。原理简单说就是在地表或者钻孔中布置金属电极通过两个供电电极向地下注入电流然后在另外两个测量电极之间记录电位差再根据一系列不同电极排列下的测量值反演得到地下电阻率的二维或三维分布。实际工程里ERT被用在很多地方滑坡体滑动面探测、堤坝渗漏通道识别、地下水污染羽追踪、冻土监测甚至考古遗址边界圈定。它的核心优势是对低阻体特别敏感比如含水破碎带、黏土夹层、盐水污染区这些目标往往比其他物性差异更容易被发现。但注意我们直接测到的只是有限数量的视电阻率或电位差值并非地下电阻率的直接照片。要从这些离散测量中还原出电阻率分布必须解决反演问题。而反演问题能不能稳定解出来很大程度上取决于一个前置指标——灵敏度分布。它就像是地下介质对每一组测量的“响应权重图”权重高的区域反演可信度高权重低的区域即使反演出结果也往往不可靠。1.2 表面电极配置与跨井电极配置的差异表面ERT大家比较熟悉就是把电极沿地表布成一条测线或者布成二维网状。测量深度受电极排列长度控制通常有效探测深度约为排列长度的三分之一到四分之一。表面ERT的优点是施工简单适合浅层目标缺点是对深部目标的灵敏度衰减非常快一旦目标深度超过排列长度的一半基本就很难有稳定的分辨能力。跨井ERTcross-borehole简称XBH则是把电极分别放在两个或更多钻孔里在其中一个钻孔的电极供电在另一个钻孔的电极测量。这样电流必须穿过两孔之间的地下空间测量结果天然对井间区域敏感。XBH特别适合需要查明两孔之间深部目标的情况比如堤基渗漏通道在两孔之间的走向、地热储层裂隙连通性、污染羽在井间深部的扩散范围。两种配置在几何上互补表面ERT看浅层横向变化XBH看井间深部纵向连通。但问题在于XBH的灵敏度分布并不是均匀覆盖整个井间区域的不同的供电-测量电极组合会形成不同形状的灵敏度“扇形区”有些区域可能重叠严重有些区域则灵敏度接近零。这一步如果不提前算清楚后面反演时很容易出现假异常。1.3 灵敏度分布为什么值得单独算灵敏度分布不是反演的最终产物却决定了反演的“视力范围”。它的定义非常直接某个地下单元电阻率发生微小变化时某一条测量道记录到的电位差会变化多少。这个变化量越大说明该测量道对该单元越“敏感”。把整个灵敏度分布算出来至少有三个实际用途第一电极排列设计。在做外业之前先在计算机里模拟不同电极间距、不同电极数量、不同测线布置下的灵敏度分布就能评估现有电极配置对目标区域的覆盖程度避免花大成本施工后却发现目标区域落在盲区里。第二反演约束和加权。反演的目标函数里通常会加入数据加权矩阵这个矩阵本质上就来自灵敏度信息。灵敏度低的区域天然欠约束如果不做处理反演迭代很容易在某些灵敏度极低的单元里产生无意义的电阻率跳动。第三解释结果的可靠性评价。反演完成后可以对照灵敏度分布图判断哪些异常体是数据真正约束住的哪些只是反演算法为了拟合数据而硬凑出来的边界。所以花点时间把2D和3D灵敏度分布算清楚远比直接套用一个反演软件更值得。下面我从物理和数学层面拆开看这个计算到底在算什么。2. 灵敏度分布背后的物理与数学2.1 什么是灵敏度电位变化如何映射电阻率变化设想地下被剖分成很多小单元每个单元有自己的电阻率。测量得到的是某一电极对之间的电位差记为V。如果把第j个单元的电阻率ρj稍微扰动一个微小量δρj那么电位差V就会跟着变化δV。两者之比∂V/∂ρj就是这个单元对这条测量的灵敏度也就是灵敏度矩阵中的元素J[i,j]。这个定义看起来简单但要直接计算却很费劲。一种最朴素的办法是扰动法把每一个单元的电阻率依次加上一个小扰动重新求一遍正演然后看测量值的变化。假设网格有十万个单元就要做十万次正演计算这在3D问题里完全不现实。所以实际项目中几乎不用扰动法而是用伴随场法或者解析灵敏度公式一次正演就能得到某条测量对所有单元的灵敏度。这里可以打一个比方。灵敏度分布就像一张“地下广播信号覆盖图”电流源是广播发射塔测量电极是收音机地下每个位置的电阻率变化相当于某个地方出现了信号干扰源。有的地方干扰源动一下收音机声音立刻变化说明这个地方灵敏度高有的地方无论怎么改收音机都没反应说明这个地方对这条测量来说是“盲区”。2.2 从正问题到灵敏度矩阵ERT正问题的控制方程是泊松方程∇·(σ∇φ) -Iδ(r-r) Iδ(r-r-)其中σ是电导率φ是电位I是供电电流r和r-分别是正、负供电电极位置。σ与电阻率ρ互为倒数。求解这个方程得到每个节点上的电位值φ再从节点电位差值中得到测量道的响应。若想得到灵敏度比较经典的路径是利用互易性和伴随场原理。对于一条“电流电极A、B测量电极M、N”的测量灵敏度可以写成某个积分形式核心是对“A-B供电情形下的场分布”和“M-N等效供电情形下的场分布”做体积分S[i,j] -∫Ω (∇φ_AB · ∇φ_MN) dΩ / I^2这里的φ_AB是真实供电电流产生的电位场φ_MN是人为“把测量电极当作供电电极”、注入同样大小电流得到的伴随电位场。正因为有互易定理这个伴随场可以通过额外一次正演得到而不是对所有单元逐一扰动。在Matlab实现时通常的做法是先组装好有限元刚度矩阵K然后对每一对供电位置求解Kn×1方程得到电位分布u_A再对每一对测量电极位置也求解一次得到伴随场u_M。这样遍历所有测量道后灵敏度矩阵的规模为“测量道总数 × 网格单元总数”正好是后续反演Jacobian矩阵的形状。2.3 2D和3D在方程与离散上的差别2D和3D的差别不只是“多一个维度”这么简单。2D计算中通常假设测线方向x和深度方向z构成计算平面垂直测线的y方向电阻率不变。此时电流源在数学上被看作一条无限长的“线源”控制方程退化为二维泊松方程求解速度快网格剖分规模通常只有几万个节点。2D灵敏度分布的结果是一个x-z剖面适合分析测线下方某个断面上的探测盲区和分辨率差异。3D计算则必须处理真正的点电流源控制方程保持三维形式。网格规模直接从二维的几万膨胀到几十万甚至上百万。3D灵敏度分布的优点是能把横向y方向的变化也包含进来尤其适合地表网状测线、井间三维阵列这些真正的体积测量。缺点是内存和计算时间呈量级上升需要稀疏矩阵求解器、并行计算等策略配合。在灵敏度计算公式的形式上2D和3D也有差异3D点源情况下场量随距离按1/r衰减积分要乘上单元体积2D线源情况下场量随距离的对数衰减计算公式中包含与测线垂直方向上的等效长度项。很多人直接拿3D公式改个网格就去算2D结果灵敏度量纲和数值都会出错。明白这些区别后Matlab实现的思路就很清晰了先确定测线维度再决定控制方程和积分形式最后再去拼灵敏度矩阵。3. Matlab实现的核心思路3.1 网格剖分与电极位置处理不管是2D还是3D第一步都是剖分网格。Matlab本身没有特别通用的人机交互剖分工具但可以通过PDE Toolbox生成基本网格也可以用distmesh这类开源工具箱做带边界的三角形或四面体网格。实际项目中我的建议是自己写结构化网格生成控制起来更灵活。比如2D剖面可以这样生成% 2D 表面ERT网格 x 0:0.5:40; % 测线方向单位m z [0:0.25:2, 2.5:0.5:10, ... % 近地表加密 12:2:40]; % 深部放宽 [X, Z] meshgrid(x, z); % X为横向坐标Z为纵向坐标向下为正 rho0 100 * ones(size(X)); % 背景电阻率 100 Ω·m网格剖分的核心原则是电极附近加密远离电极的区域逐步放宽。因为在点电流源附近电位梯度变化非常剧烈如果网格太粗灵敏度峰值区域会被严重平滑掉算出的灵敏度分布会丢失局部细节。电极位置建议单独用一个矩阵保存不要直接写在网格函数里% 表面电极沿测线每2m一个电极 electrode_x 0:2:40; electrode_z zeros(size(electrode_x)); % 跨井电极两井位于x0和x20深度从0到30m间距2m borehole1_x zeros(1,16); borehole2_x 20 * ones(1,16); borehole_z 0:2:30;电极一般放在节点上如果电极落在单元内部需要把电流源分配到周围节点上。这是容易出问题的地方后面避坑环节我会专门展开。3.2 正问题求解有限元或解析近似ERT正问题求解有两种常见路线有限元法和解析/半解析近似法。解析近似只适用于均匀半空间或者层状介质快速验证灵敏度公式还可以但实际项目里地形起伏、复杂地下结构都需要有限元。Matlab里可以用自编有限元也可以用PDE Toolbox。自编有限元的好处是能精确控制关于灵敏度的积分我通常采用这种路线流程如下根据网格节点坐标生成单元局部刚度矩阵。组装全局稀疏刚度矩阵K同时处理边界条件。构建源项向量q位置对应供电电极节点。求解线性方程组K*u q得到电位场。关键点在于刚度矩阵K与电阻率有关。默认情况下每个单元有自己的电阻率值单元内电导率取常数积分时用单元中心值。解大型稀疏方程组时Matlab自带的直接求解器比如K\q对于二维问题足够三维问题节点数多了以后建议改用共轭梯度法等迭代求解器配合不完全Cholesky预处理。边界条件处理尤其要注意顶部地表一般设为绝缘边界即电流不能穿过地表侧面和底部是截断边界理论上应该延伸到无穷远实际计算中通常把网格范围取得足够大或者使用Robin型混合边界来吸收外行能量。最简单稳妥的做法是把网格向外扩展使得目标测区距离边界至少5倍以上的电极排布尺寸。3.3 灵敏度矩阵的组装与可视化在有限元正演结果基础上灵敏度矩阵可以逐条测量道计算。核心是拿到两个电位场梯度一个来自供电电极另一个来自测量电极的“伴随供电”。下面给一段Matlab伪代码展示2D情况下灵敏度计算的骨架nMeas size(measurement_pairs, 1); % 测量道总数 nCell size(elements, 1); % 有限元单元总数 S zeros(nMeas, nCell); % 灵敏度矩阵 I 1; % 供电电流归一化为1A for k 1:nMeas % 从测量道信息中提取供电电极和测量电极编号 A measurement_pairs(k, 1); B measurement_pairs(k, 2); M measurement_pairs(k, 3); N measurement_pairs(k, 4); % 求解供电电位场 qA sparse(nodeIndex(A), 1, I, nNode, 1); qB sparse(nodeIndex(B), 1, -I, nNode, 1); phi_AB K \ (qA qB); % 求解测量电极对应的伴随电位场 qM sparse(nodeIndex(M), 1, I, nNode, 1); qN sparse(nodeIndex(N), 1, -I, nNode, 1); phi_MN K \ (qM qN); % 计算每个单元两个场的梯度点积 [gradX_AB, gradZ_AB] gradientField(phi_AB, nodes, elements); [gradX_MN, gradZ_MN] gradientField(phi_MN, nodes, elements); S(k, :) -(gradX_AB .* gradX_MN gradZ_AB .* gradZ_MN) ... * eleArea / (I^2); end这段代码看起来简单实际运行起来有两个地方要特别注意。第一梯度计算不能用Matlab内置的gradient函数直接对节点电位做差分那样精度不够应该按有限元形函数在单元内做积分先算单元内积分点上的梯度再聚合到单元中心。第二灵敏度量纲不同方案差异很大有人用电阻率灵敏度有人用电导率灵敏度还有人用视电阻率灵敏度。计算前先明确自己需要的是哪一种。大多数ERT反演软件最终使用的是对数电阻率灵敏度即∂logρa/∂logρ这样动态范围更均匀。可视化部分二维用pcolor或者contourf画x-z剖面三维用slice函数切不同深度的水平切片也可以用来回旋转的volumetric图展示“花瓣形”灵敏度覆盖。使用imagesc时要注意坐标方向通常把z轴翻转让深度向下显示否则看图习惯会完全反掉。3.4 计算参数与性能取舍Matlab里计算灵敏度最耗时的部分不是灵敏度积分本身而是大量正演求解。比如一条测线上有24个电极四极法测量组合可能有几百上千条每个测量道要解两次正演方程乘以求解一次线性系统的时间总耗时很可观。所以第一个参数建议是不要对所有可能的四极组合都无脑算。先做一次道筛选结合电极排列方式和目标深度挑出彼此重叠度低、且对目标区域灵敏度较高的测量组合。这一步能减少一半以上的正演次数。第二个参数是背景电阻率值。灵敏度积分结果与背景电阻率有关在均匀半空间模型中灵敏度值与背景电导率成反比。实际计算时要用你预计的地下平均电阻率作为初始背景不要用随意给定的数值否则算出来的灵敏度绝对量级会对后续加权矩阵产生误导。第三个参数是并行。Matlab的parfor循环可以很好地套在测量道遍历上因为每个测量道的两次正演相对独立。在四核笔记本上用parfor能把总耗时压到原来的三分之一以下。如果每一条测量道都调用同一个刚度矩阵K记得在循环外把K的稀疏结构固定好避免重复分解内存。4. 表面与跨井配置的灵敏度分布对比4.1 表面ERT的灵敏度特征表面ERT当电极沿测线布设时灵敏度分布有一个非常典型的特点浅层区域的灵敏度值高且集中随深度增加快速衰减。如果测线长度为L电极间距为a那么大约在深度L/3到L/4以下灵敏度通常已经衰减到浅层峰值的百分之几。不同测量排列也会改变灵敏度形态。以温纳Wenner排列为例它的灵敏度在剖面中心位置形成一个较宽的“穹顶状”高值区中心点深度约为电极间距a的水平长度对应中等深度浅部和深部灵敏度相对均匀。偶极-偶极排列的灵敏度则分裂成几个瓣状高值区横向分辨率高但由于供电偶极和测量偶极之间的距离增大浅层噪声容易影响深部测量。斯伦贝谢排列介于两者之间适合兼顾横向和纵向分辨率的场景。这给表面ERT带来一个实际约束如果你只用一个电极间距做测量灵敏度覆盖就像一把“近地表望远镜”——浅层很亮深层全黑。解决办法是使用多个电极间距离散布置将不同间距的灵敏度分布叠合起来才能把深部目标区域照亮。4.2 跨井ERT的灵敏度特征跨井ERT由于电极分布在两个钻孔中灵敏度分布形态和表面ERT完全不同。以一个A井供电、B井测量的单极-单极pole-pole或单极-偶极pole-dipole排列为例灵敏度以供电电极和测量电极之间的连线为轴形成一个沿井轴方向延伸的“香蕉形”或“扇形”高值区。多条不同深度电极对的灵敏度相互叠加后井间深部的覆盖明显改善。但跨井XBH也有一个容易忽视的盲区——两井连线的中间区域。尤其当供电电极和测量电极在各自井内对应位置相近时高灵敏度区域偏向两井附近井间中点处灵敏度相对偏低。如果两井之间的距离过大比如超过目标深度井间中部的灵敏度会趋于微弱反演结果在那一带极不稳定。实际项目里我常用的跨井电极模式有三种单极-单极供电电极在一个井中测量电极在另一个井中。覆盖范围广深度方向分辨率偏低。偶极-偶极两个井中分别取相邻电极对供电和测量。横向灵敏度分辨率高但数据量爆炸需要筛选。跨井表面联合在一个井内供电在另一个井内和地表同时测量。用于补足井间浅部与深部之间的过渡区域。把表面ERT和XBH的灵敏度叠加到同一张图上会发现它们的高灵敏度区往往是斜交的表面ERT高值区靠近地表中线XBH高值区沿井轴呈扇形发散。联合反演能显著缩小盲区这也是近年很多项目采用“井地联合ERT”的原因。4.3 2D与3D可视化解读2D灵敏度分布可视化时常以x为横轴、z为纵轴将每条测量道的灵敏度叠加后取绝对值或归一化值。叠加后的图像会明显呈现出“倒三角形”或者“扇形”的高值区。对于表面ERT高值区顶点朝下中心位置灵敏度最高对于XBH高值区从两个井轴向外发散中间形成相对低值带。可以生成如下表所示的对比配置类型高灵敏度区形态主要盲区位置适合目标表面ERTWenner浅层穹顶状深部、测线两端下方浅层分层、空洞表面ERT偶极-偶极多瓣状横向分辨好瓣间低值区、深部横向不连续体跨井ERTpole-pole井间扇形、沿井轴延伸两井连线中部、深部外侧井间目标深部追踪跨井ERT偶极-偶极沿井轴强聚焦、瓣状交叉井间中心敏感度偏低、远井端高分辨率井间成像3D可视化则更复杂。对井间XBH三维布置可以用slice函数同时切出xz和yz剖面观察灵敏度在三维空间的分布。结论上3D灵敏度分布通常比2D更“瘦”因为点源不像线源那样在横向无限延伸实际影响范围集中在一个三维“苹果形”体内。这个差异直接决定了3D反演对横向约束能力比2D强但也意味着反演问题更病态需要更多正则化约束。从Matlab输出的角度建议3D灵敏度分布离散保存为一个ncells×1的向量每个单元索引对应中心坐标。可视化时先创建三维体数据vol再把灵敏度值赋到相应体素用slice或者volshow展示。注意灵敏度值正负差异很大取绝对值前先观察原始符号分布往往会发现负值区域同样包含强烈信息比如某些测量排列在供电和测量电极之间存在负灵敏度带直接取pcolor绝对值会把这种结构抹掉。5. 实操中的常见问题与避坑5.1 网格密度与电极大小这个问题我踩过几次坑。一是网格剖分过于均匀电极附近网格不够细。当电流源落在粗网格节点上时计算得到的电位场在电极周边会产生较大误差进而导致灵敏度在近电极区域出现奇异的正负交替图案。解决办法是电极周围局部加密可以用高斯分布形式的节点密度离电极越近节点间距越小。例如表面ERT沿测线方向电极附近节点间距取0.1~0.2m远离电极处放宽到1~2m然后过渡到深部指数放大。二是电极单元尺寸不能为零。在有限元里点电流源被分配到节点上如果网格单元面积太小但求解精度不够电流密度会异常高。通常把电极所在的单元尺寸控制在电极间距的十分之一到二十分之一既能保证精度又不至于让单元数量爆炸。5.2 边界条件与无穷远边界边界条件对灵敏度分布的影响非常大。如果网格截断边界离测量区域太近边界上电位回弹会造成灵敏度分布出现“伪高值”边界效应。一个快速检查方法把灵敏度分布图画出来观察边缘是否有一条沿模型边界的明显亮线如果有几乎可以确定是边界反射。我常用的处理顺序是先扩大网格范围让模型四周比测区外扩3~5倍再对侧面和底面施加混合边界Robin边界以模拟无穷远半空间。Matlab中PDE Toolbox支持这类边界条件自编有限元时也可以通过边界积分项添加。如果只想快速验证最简单的做法是把模型底部和侧面设为“零电位边界”但要确保网格边界足够远远到边界电位已经衰减到可忽略。5.3 归一化与矩阵病态灵敏度矩阵数值跨度过大是另一个常见坑。同一个测量道对不同单元的灵敏度可能相差几个数量级如果不做归一化反演时高灵敏度单元会主导目标函数低灵敏度单元几乎没有更新机会。处理方法是采用对数灵敏度。具体来说如果用电阻率ρ作为参数灵敏度矩阵可以改写成∂V/∂logρ ρ·∂V/∂ρ这样灵敏度数值在空间上的分布更稳定且天然与电阻率尺度无关。视觉呈现上也可以对灵敏度图做log10取对数后再画等值线否则最大峰值附近的色标会把整个图渲染成一片亮白。另外注意灵敏度矩阵的行和列方向。在反演程序中雅可比矩阵的行对应测量道列对应单元。如果你把自己算出的灵敏度矩阵转置一下再塞进别人写的反演代码会得到完全错误的更新方向。我在联调时用了一整天才发现是这个低级问题。5.4 内存和速度优化3D灵敏度的存储量很容易失控。假设网格单元有50万个测量道有2000条那么灵敏度矩阵大小为2000×500000如果用double型存储需要8GB内存这还不包括正演求解过程中的临时变量。三个优化建议用稀疏矩阵存储灵敏度矩阵中的非零元通常不是全部但比例也不低。可以先判断哪些单元距离所有电极对过远灵敏度衰减到可以忽略把这些单元剔除或标记为“低灵敏度冻结区”只保留主要影响区域的灵敏度。分块计算与磁盘映射把测量道分成若干块每块算完灵敏度子矩阵后先存到.mat文件最后再合并。这样运行内存可以被压到很小。用parfor并行测量道之间独立直接在循环外层用parfor。共享的刚度矩阵K可以在循环前用decomposition()函数做矩阵分解循环内重复调用分解后的对象求解速度提升明显。我测试过一个典型的3D XBH算例30m井距、两井各16个电极、80万网格单元、1500条测量道。单核循环跑了约3小时加parfor后并行到6个worker缩短到40分钟内存峰值控制在6GB以内。如果进一步用低灵敏度冻结区剔除还能再压缩20%。最后的实操体会把灵敏度分布算完以后你会发现它比反演结果本身更能说明问题。我个人的习惯是在正式反演前先做一张灵敏度覆盖图叠加测区内的工程地质剖面凡是目标深度落在低灵敏度区的要么重新设计电极排列要么直接在报告中标注为“不可靠区域”。这一步做完反演迭代里许多无意义的电阻率抖动其实都能提前规避。另外还有一个非常实用的小技巧把灵敏度矩阵单独抽出一列来画图这一列代表“地下某个固定单元电阻率变化对所有测量道的响应”。它直接揭示了数据里到底有多少条测量道与该单元相关适合用来排查测线端部效应和孤立电极的异常影响。把每一列都扫一遍你就能找到哪些电极对是“冗余”的去掉它们几乎不影响反演质量却能省下不少计算时间。这套Matlab流程跑通之后后续还可以扩展把2D灵敏度结果沿y方向复制成伪3D用来评估测线网络的横向覆盖也可以把灵敏度矩阵输入带正则化的反演框架实现从“灵敏度分析”到“反演成像”的平滑过渡。先算清灵敏度再谈反演这条路线值得每个做ERT的人走一遍。