企业官网建设流程全解析

用Matlab把电阻率层析成像（ERT）的灵敏度分布算清楚，这件事看着偏理论，却是决定反演结果可信度的关键一步。最近我把表面电极和跨井电极（cross-borehole，XBH）配置下的2D/3D灵敏度分布完整跑了一遍，从电极坐标、网格剖分到正演求解，再到灵敏度矩阵组装和可视化，算是把整条链路打通了。这篇内容就是一次实操复盘，重点讲清楚灵敏度分布到底在算什么、Matlab里怎么搭这套计算流程、以及表面ERT和XBH在灵敏度覆盖上的本质差异，刚接触ERT数据模拟或者做反演前设计的朋友可以直接参考。

1. 认识ERT、跨井测量与灵敏度分布

1.1 电阻率层析成像到底在测什么

电阻率层析成像（Electrical Resistivity Tomography，ERT）是一种以地下岩土体电阻率差异为成像目标的地球物理方法。原理简单说就是：在地表或者钻孔中布置金属电极，通过两个供电电极向地下注入电流，然后在另外两个测量电极之间记录电位差，再根据一系列不同电极排列下的测量值反演得到地下电阻率的二维或三维分布。

实际工程里，ERT被用在很多地方：滑坡体滑动面探测、堤坝渗漏通道识别、地下水污染羽追踪、冻土监测，甚至考古遗址边界圈定。它的核心优势是对低阻体特别敏感，比如含水破碎带、黏土夹层、盐水污染区，这些目标往往比其他物性差异更容易被发现。

但注意，我们直接测到的只是有限数量的视电阻率或电位差值，并非地下电阻率的直接照片。要从这些离散测量中还原出电阻率分布，必须解决反演问题。而反演问题能不能稳定解出来，很大程度上取决于一个前置指标——灵敏度分布。它就像是地下介质对每一组测量的“响应权重图”，权重高的区域反演可信度高，权重低的区域即使反演出结果也往往不可靠。

1.2 表面电极配置与跨井电极配置的差异

表面ERT大家比较熟悉，就是把电极沿地表布成一条测线，或者布成二维网状。测量深度受电极排列长度控制，通常有效探测深度约为排列长度的三分之一到四分之一。表面ERT的优点是施工简单，适合浅层目标；缺点是对深部目标的灵敏度衰减非常快，一旦目标深度超过排列长度的一半，基本就很难有稳定的分辨能力。

跨井ERT（cross-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)); % 跨井电极：两井位于x=0和x=20，深度从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。自编有限元的好处是能精确控制关于灵敏度的积分，我通常采用这种路线，流程如下：

1. 根据网格节点坐标生成单元局部刚度矩阵。
2. 组装全局稀疏刚度矩阵K，同时处理边界条件。
3. 构建源项向量q，位置对应供电电极节点。
4. 求解线性方程组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，高值区从两个井轴向外发散，中间形成相对低值带。可以生成如下表所示的对比：

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的人走一遍。