企业官网建设流程全解析

做ERT 的朋友应该都有体会：正演好算，反演难调；但如果一开始连灵敏度分布都没吃透，后面不管是设计电极阵列、决定井距，还是给反演加模型权重，都会像闭着眼过河。我写这篇就是想用一个基于 MATLAB 的完整实例，把ERT在 surface（表面）和 XBH（跨井）两类电极配置下的 2D/3D 灵敏度分布计算讲透。代码不依赖任何商业反演软件，核心是解析核函数，改几行坐标就能换阵列、换井距、换测线长度，非常适合用来做观测系统布设前的快速评估。

这篇文章适合谁？做高密度电法、跨孔电阻率 CT 的学生和工程师，想搞懂“为什么跨井中间灵敏度高、两侧是盲区”的人，以及想把灵敏度结果当作先验信息喂给反演算法的MATLAB 玩家。我会尽量把公式、代码、坑都摊开讲，保证你跟着敲一遍就能出图，并能自己改动阵列去对比。

1. 先搞清楚：灵敏度分布到底是什么，为什么值得算

1.1 物理含义：它就是每个网格体素对测量值的“贡献权重”

灵敏度的严格定义是：地下某处电阻率发生微小变化时，某个观测电位差（或视电阻率）随之变化的程度。用数学语言说，它是观测数据对模型参数的导数，也就是反演理论里的 Frechet 导数。如果你没接触过反演，可以打个比方：灵敏度图就像一张“投票权重图”——每次测得的电压数据相当于一次投票，地下每个位置的电阻率变化都按自己的权重参与这次投票。权重大的地方，数据对那里的电阻率变化敏感，反演时那里就好分辨；权重小甚至为零的地方，数据基本“看不到”那里，反演结果就只能靠周围区域插值硬撑。

在均匀半空间背景下，这个权重可以用解析式表达。对一个四电极装置（A、B供电，M、N测电位差），灵敏度核函数本质上等于两组“伪电场”的点积负值：

S(r) = - [ (r - r_A)/|r - r_A|³ - (r - r_B)/|r - r_B|³ ] · [ (r - r_M)/|r - r_M|³ - (r - r_N)/|r - r_N|³ ]

这里第一项是 A、B 两个供电点在介质中产生的电流密度场的空间分布，第二项是假定 M、N 反过来作为供电点时产生的伴随场。两个场在同一点上做点积，就得到了该点的灵敏度。这样的形式非常方便 MATLAB 向量化，完全不需要做数值扰动，也不用建有限元网格。

1.2 为什么“看灵敏度”比“看正演电位”更能指导野外工作

正演电位告诉你某个观测值大概是多少，但没法直接告诉你哪块区域对测量有贡献、哪块区域是“透明的”。灵敏度分布则直接回答三个工程问题：这个配置到底能探多深？能覆盖多宽？井间哪些区域是反演盲区？

举个实际例子：表面测线布设时，很多人都以为“电极距拉大，探测深度就自动变大”。确实，极距增大后灵敏度峰值会向深部移动，但灵敏度的绝对值也在同步衰减。你拉大极距换来的浅层信息会明显丢失，深层虽然能“看到”了，但分辨率下降得很厉害。用灵敏度图一对比，这些权衡一目了然。跨井观测更是如此，A、B 电极在一口井里，M、N 在另一口井里，灵敏度带其实主要沿两根电极之间的连线分布，中间覆盖好，两侧近乎空白。如果目标体恰好偏向井壁附近，可能横跨多条电极深度组合的数据都贡献不大。这些判断，光看测线布设图是看不出来的。

1.3 灵敏度矩阵在反演里的位置

很多人容易忽略一点：灵敏度分布不仅仅是评价观测系统用的，它本身就是反演方程里的核心矩阵。在 Gauss-Newton 类迭代反演中，每次更新模型都要用到雅可比矩阵，而雅可比矩阵的每一行对应的正是某一组电极测量在不同空间位置上的灵敏度。把解析灵敏度分布算出来，实际上就拿到了均匀背景下的雅可比矩阵初值，可以用于：

• 观测数据加权：灵敏度高的测点给高权重，低灵敏度测点降权，避免反演被冗余数据带偏；
• 模型参数加权：用灵敏度的倒数做模型协方差，让反演在低灵敏度区域更保守，减少虚假异常；
• 正则化参数设计：用灵敏度分布做空间变化的阻尼因子，让深部和盲区不至于被过度拟合。

所以，别看它叫“灵敏度分布”，它贯穿了观测设计和反演两个阶段，值得你花时间把它彻底算明白。

2. 表面和跨井配置的灵敏度形态差异：为什么差的不是一点半点

2.1 表面电极配置：浅部灵敏，深度快速衰减

表面配置，也就是所有电极都放在地表，是高密度电法中应用最广的模式，比如温纳（Wenner）、施伦贝谢（Schlumberger）、偶极-偶极（Dipole-Dipole）等。它们的共同特点是电流从地表注入，在均匀半空间里电流密度向深部迅速扩散衰减，所以灵敏度天然就是浅部高、深部低的形态。

以温纳阵列为例，四个电极在一条直线上等间距排列，A 在最外侧供电，M、N 在中间测量，B 在最外侧回流。这种装置的灵敏度分布有一个非常容易踩坑的特征：浅部灵敏度不是单调的，而是呈正负相间的区域分布，像一只展开的蝴蝶翅膀。负灵敏度并不代表“反着探测”，它只是说明该处电阻率升高反而会引起视电阻率降低，这与测量电位差的定义方式和几何位置有关。这个特征在做反演时尤其重要，如果你忽略符号直接把灵敏度取绝对值去解释，很容易在浅部区域产生错误的覆盖判断。

表面装置的探测深度大致可以用一个经验区间衡量：在极距 a 不太大的情况下，温纳装置灵敏度峰值深度大致在 0.5a~1.0a 左右。极距增大，灵敏度峰值会下移，但峰值本身的幅值会衰减，深层分辨率上限也会被拉低。这也是为什么表面高密度电法越往深部走，解释结论越要谨慎。

2.2 跨井XBH配置：井间覆盖好，但视线狭窄

跨井ERT（Cross-Borehole ERT）把供电电极和测量电极分别放入两个钻孔内，电流路径从一口井穿到另一口井，天然绕开了地表电流密度衰减快的问题，对井间区域的分辨率远远优于表面观测。代价是覆盖范围非常集中于“发射电极到接收电极”之间的窄带区域，而且这个窄带的方向和展布取决于电极深度组合。

常用的跨井布极方式是 AB 在同一口井、MN 在同一口井的另一深度段，或者 AB 在1号井、MN 在2号井。无论哪种，灵敏度形态都有一个共同特点：在连接发射点与接收点的空间带上，灵敏度值最高，在这条带的宽度方向上快速衰减。你可以把它想象成探照灯的光柱——照到的地方亮，没照到的地方完全黑。单组跨井电极组合的探测范围根本无法覆盖整个井间区域，必须靠多组不同深度的发射-接收电极组合把灵敏度“扫”出来。

在实际三维场景中，如果两排井组成一个井阵，跨井灵敏度分布在三维空间里不再只是一条带，而是围绕井间区域的一系列类似“管束”的敏感空间，整体形态更像从一个井眼向另一个井眼辐射的高斯型隧道。这也是跨井ERT做盐碱地监测、微塑料污染扩散追踪、岩溶管道探测时，一定要求三维采测的原因——2D剖面会丢掉大量空间信息。

2.3 两类配置的对比速查表

做完这个对比你就该明白：没有“万能电极配置”，只有“适不适合当前目标”的配置。选配置之前先算一下灵敏度，是最省钱的方案评估手段。

3. MATLAB实现：解析核函数计算灵敏度的核心流程

3.1 为什么选解析法而不是有限元/有限差分

很多人一看到“计算灵敏度”就想到要建复杂网格、写有限元正演、再求逆，其实完全没必要。题目要算的是均匀半空间背景下的灵敏度核函数，这是有精确解析解的，用插值伪电场点积的方法几十行代码就能搞定。

那什么时候必须上有限元？当地表有起伏、地下有强非均匀体、电极处于钻孔填充材料等复杂介质中时，解析解失效，只能靠数值正演逐点扰动求灵敏度。但做观测系统设计时，我们通常关心的是灵敏度分布的“形态”和“相对量级”，均匀介质的解析核已经足够。等到要精确反演真实数据时，再切换到有限元正演也来得及。先用解析法完成快速扫参，用一个小时筛出三五种候选配置，再挑两三种上有限元细算，效率会高得多。

3.2 从点源电位到灵敏度核的推导主线

在均匀介质中，位于 r_s 的点电流源产生的电位场为 φ(r) ∝ 1/|r - r_s|。对应的电场强度（梯度）为 ∇φ(r) ∝ (r - r_s) / |r - r_s|³，这正是代码里反复出现的核心向量。

对于 A、B 供电，M、N 测量电位差的情况，利用互易原理和 Born 近似，灵敏度核可以写为“AB 伪电场”与“MN 伪电场”的点积。所谓 MN 伪电场，就是把测量电极 M、N 也想象成一对“虚拟供电电极”，让电流从 M 流入、N 流出，再算出该电流场在空间各点的分布。两个伪电场在同一个网格点上做内积，再统一加负号，就得到灵敏度值。

代码实现时不用显式写出背景电阻率系数，因为系数对所有网格点是一个常数，最终在归一化绘图时会被约掉。如果你想得到真实的物理量纲值，再乘回去也不难。重点是把三个坐标系下的向量分量拆对：二维就是 Ex 和 Ez，三维还要加上 Ey。

3.2.1 半空间地表电极的处理

地表电极有“镜像效应”：电流源在地表时，地下部分电位相当于把源加倍。写成格林函数就是 φ = 1/(2πσr) 而不是全空间的 1/(4πσr)。不过在灵敏度核公式里，所有项都含同一个公共系数，归一化后仍然被消去，所以代码里不需要刻意区分“地表源”和“地下源”，至少在做相对分布图时不用。如果你真的需要绝对数值，把表面电极的梯度项再乘一个 2 即可。

3.2.2 关于2D剖面不是严格2D的问题

必须提醒一句：ERT 真实的物理过程是三维的，所谓“2D 剖面灵敏度”其实是一个沿着垂直剖面方向做了某种积分或假设的投影结果。严格的处理叫 2.5D，即地质模型在 y 方向不变，但点源产生的三维场在剖面上采样。很多教材里直接画出“2D 灵敏度剖面”，本质上是把三维核函数在 y=0 处切片显示，或者沿 y 方向做了积分。

我这里两种都可以满足：如果只是想快速看形态，直接在 (x, z) 网格上算切面值就行；如果想更严谨一点，可以把三维灵敏度在 y 方向从 -L 到 L 做数值积分，得到等效的二维灵敏度。代码上差别不大，只是多加一层积分循环或累加。

3.3 网格剖分与精度控制

网格剖分是这类计算最容易被忽略、也最容易出错的地方。原则上，网格范围要覆盖住所有灵敏度非零的区域；步长则决定了分辨率计算的精细程度。

实际操作中我的习惯是：

• 网格横向范围至少向外扩出最大电极距的 2 倍，比如最大极距 3m，x 方向从 -6m 到 6m；
• 网格纵向范围取最大电极距的 1.5~2 倍，太深的地方灵敏度几乎为零，白白增加计算量；
• 步长用测线最小电极距的 1/20~1/10。比如最小电极距 0.5m，步长取 0.025~0.05m 即可。

这里尤其要注意：表面配置的灵敏度浅部变化剧烈，如果步长取太大，会把正负相间的区域混在一起，看起来像一团噪声；适当加密一下，图形会干净得多。而三维计算受限于内存，步长要适当放宽，一般取最小电极距的 1/10 左右，先用粗网格看形态，再在关键区域局部细化。

4. 实操一：表面2D灵敏度分布的计算与解读

4.1 电极坐标与网格设计

这一节我以温纳阵列为例，四个电极从 0m 开始排，间距 a=0.7m，所以 A=0m，M=0.7m，N=1.4m，B=2.1m。网格x方向从 -1.4m 到 3.5m，z方向从 0m 到 2.8m，步长 dx=0.05m。这里 z 以地表为 0、向下为正，画图时会比较自然。

网格确定后，依次计算每个网格点到 A、B、M、N 的距离。这段代码的每一步都对应 3.2 节的解析式，没有神秘操作。

4.2 核心代码：温纳阵列2D灵敏度

clear; clc; a = 0.7; % 温纳极距 xA = 0; zA = 0; % 供电电极 A xB = 3*a; zB = 0; % 供电电极 B xM = a; zM = 0; % 测量电极 M xN = 2*a; zN = 0; % 测量电极 N % 计算网格（地表为z=0，向下为正） dx = 0.05; xg = -2*a:dx:5*a; zg = 0:dx:4*a; [XX, ZZ] = meshgrid(xg, zg); % 网格点到四个电极的距离 rA = sqrt((XX - xA).^2 + (ZZ - zA).^2); rB = sqrt((XX - xB).^2 + (ZZ - zB).^2); rM = sqrt((XX - xM).^2 + (ZZ - zM).^2); rN = sqrt((XX - xN).^2 + (ZZ - zN).^2); % 伪电场分量（A、B和M、N产生的空间场） ExA = (XX - xA)./rA.^3; EzA = (ZZ - zA)./rA.^3; ExB = (XX - xB)./rB.^3; EzB = (ZZ - zB)./rB.^3; ExM = (XX - xM)./rM.^3; EzM = (ZZ - zM)./rM.^3; ExN = (XX - xN)./rN.^3; EzN = (ZZ - zN)./rN.^3; % AB伪电场与MN伪电场的点积 ExAB = ExA - ExB; EzAB = EzA - EzB; ExMN = ExM - ExN; EzMN = EzM - EzN; Sen = -(ExAB.*ExMN + EzAB.*EzMN); % 绘图：用上下限截断颜色条，防止正负极值把图形压平 figure('Color','w'); imagesc(xg, zg, Sen); axis xy; axis equal tight; colormap(jet); colorbar; clim([-max(abs(Sen(:)))/4, max(abs(Sen(:)))/4]); xlabel('x / m'); ylabel('z / m'); title('Surface Wenner 2D Sensitivity');

这段代码跑完，你会看到一张正负相间的四瓣型图案，中心区域靠近地表，灵敏度绝对值较大，深部颜色逐渐变浅。核心规律是：灵敏度正负区域的分界大致在电极排列方向的垂直面上，负区往往出现在测线两端外侧的浅部区域。如果图上颜色范围不对、全是刺眼的红蓝尖刺，多半是 colorbar 范围设置问题——我上面用了 max 的四分之一作为截断值，可以根据波形调整。

4.3 结果解读和参数扫描

拿到这图后有件事一定值得做：改变极距 a，看灵敏度峰值深度怎么移动。把 a 从 0.3m 改到 1.5m，你会直观地看到峰值深度和展宽的变化。这个试验比只看探测深度的经验公式有用得多，因为它连分辨率的横向展宽都一起展示了。

我做参数扫描时习惯把多个 a 的灵敏度图保存成 subplot 拼在一张图里，看哪个极距能把目标深度区域包进“高灵敏度带”内，就选哪个极距作为野外采集配置。整个过程只需要改一行 a，重跑一次，几十秒出结果，非常顺手。

5. 实操二：跨井XBH的2D与3D灵敏度分布

5.1 2D跨井配置的坐标与代码改动

跨井配置的坐标设定就四个字：把电极挪到井里去。假设 1 号井在 x=-1.5m，2 号井在 x=1.5m，两口井深度深入到 3m。供电电极 A、B 都放在 1 号井内，分别位于深度 0.6m 和 1.4m；测量电极 M、N 都放在 2 号井内，也分别位于 0.6m 和 1.4m。

代码主体和上一节几乎完全相同，只需要替换电极坐标。注意 z 方向依然是向下为正，所以电极坐标是正深度值。

% 跨井XBH 2D灵敏度 xA = -1.5; zA = 0.6; xB = -1.5; zB = 1.4; xM = 1.5; zM = 0.6; xN = 1.5; zN = 1.4; xg = -3:0.05:3; zg = 0:0.05:3; [XX, ZZ] = meshgrid(xg, zg); rA = sqrt((XX - xA).^2 + (ZZ - zA).^2); rB = sqrt((XX - xB).^2 + (ZZ - zB).^2); rM = sqrt((XX - xM).^2 + (ZZ - zM).^2); rN = sqrt((XX - xN).^2 + (ZZ - zN).^2); ExA = (XX - xA)./rA.^3; EzA = (ZZ - zA)./rA.^3; ExB = (XX - xB)./rB.^3; EzB = (ZZ - zB)./rB.^3; ExM = (XX - xM)./rM.^3; EzM = (ZZ - zM)./rM.^3; ExN = (XX - xN)./rN.^3; EzN = (ZZ - zN)./rN.^3; ExAB = ExA - ExB; EzAB = EzA - EzB; ExMN = ExM - ExN; EzMN = EzM - EzN; Sen = -(ExAB.*ExMN + EzAB.*EzMN); figure('Color','w'); imagesc(xg, zg, Sen); axis xy; axis equal tight; colormap(jet); colorbar; clim([-max(abs(Sen(:)))/5, max(abs(Sen(:)))/5]); xlabel('x / m'); ylabel('z / m'); title('Cross-Borehole 2D Sensitivity');

跑完你会看到一条从 (x=-1.5, z=0.6) 指向 (x=1.5, z=0.6) 的高灵敏度条带，在 z=1.4 附近还有一条大体平行的条带，两条带之间有时会出现较弱的正负过渡区。这个形态说明跨井测量对“发射-接收对”连线的覆盖非常好，但两条连线之间的区域覆盖率明显不足。所以实际采集中要加密不同深度组合的电极对，才能保证井间区域没有覆盖空洞。

这只是一种简化展示。更严格的 2D 剖面灵敏度可以通过 3D 核函数沿 y 方向积分得到。做法是把上面的 Sen 计算放到两层循环里，对不同的 y 值各算一次三维切面，再乘以 dy 累加即可。在 MATLAB 里这种累加几十个 y 层也就几秒钟的事，并不会很慢。

5.2 3D跨井灵敏度：核心代码与可视化

三维的计算和二维的差别只在坐标从 (x,z) 变成了 (x,y,z)，电场向量增加了一个 Ey 分量。这时必须用 ndgrid 生成三个方向的三维网格。注意不要用 meshgrid 处理三维，meshgrid 在三维下的维度排列规则非常容易搞错，ndgrid 更直观。

clear; clc; % 井位置和电极深度 xA = -1.5; yA = 0; zA = 0.6; xB = -1.5; yB = 0; zB = 1.4; xM = 1.5; yM = 0; zM = 0.6; xN = 1.5; yN = 0; zN = 1.4; % 三维网格：注意用ndgrid xg = -3:0.15:3; yg = -3:0.15:3; zg = 0:0.15:3; [XX, YY, ZZ] = ndgrid(xg, yg, zg); % 距离计算 rA = sqrt((XX-xA).^2 + (YY-yA).^2 + (ZZ-zA).^2); rB = sqrt((XX-xB).^2 + (YY-yB).^2 + (ZZ-zB).^2); rM = sqrt((XX-xM).^2 + (YY-yM).^2 + (ZZ-zM).^2); rN = sqrt((XX-xN).^2 + (YY-yN).^2 + (ZZ-zN).^2); % 伪电场三分量 ExA = (XX-xA)./rA.^3; EyA = (YY-yA)./rA.^3; EzA = (ZZ-zA)./rA.^3; ExB = (XX-xB)./rB.^3; EyB = (YY-yB)./rB.^3; EzB = (ZZ-zB)./rB.^3; ExM = (XX-xM)./rM.^3; EyM = (YY-yM)./rM.^3; EzM = (ZZ-zM)./rM.^3; ExN = (XX-xN)./rN.^3; EyN = (YY-yN)./rN.^3; EzN = (ZZ-zN)./rN.^3; ExAB = ExA - ExB; EyAB = EyA - EyB; EzAB = EzA - EzB; ExMN = ExM - ExN; EyMN = EyM - EyN; EzMN = EzM - EzN; Sen = -(ExAB.*ExMN + EyAB.*EyMN + EzAB.*EzMN); % 可视化：取一个阈值画等值面 thresh = 0.3 * max(Sen(:)); figure('Color','w'); p = patch(isosurface(xg, yg, zg, Sen, thresh)); isonormals(xg, yg, zg, Sen, p); p.FaceColor = 'interp'; p.EdgeColor = 'none'; view(3); camlight; lighting gouraud; xlabel('x / m'); ylabel('y / m'); zlabel('z / m'); title('Cross-Borehole 3D Sensitivity');

三维结果中，你会看到灵敏度等值面像一个从 1 号井延伸到 2 号井的喇叭形管束，在井周围较粗，在两井中间部位收敛。原因是每个点源的场强在源附近急剧增大，所以离电极越近灵敏度越高。这个形状对三维反演解释很重要：如果井间距很大，井间远端区域的灵敏度等值面会非常稀疏，显示反演结果时要格外小心。

5.2.1 三维计算的内存控制

三维计算最大的敌人是网格太密导致内存暴涨。假设网格是 41×41×21，大约是 3.5 万个节点，每个节点需要至少 7~8 个 double 数组，每个数组又是 3.5 万个 8 字节数据，总共几 MB，很轻松。但如果步长加密到 0.05m，网格变成 121×121×61，接近 90 万个节点，每个变量数组就是 7MB，整段代码的临时变量加起来轻轻松松超过 200MB，低配电脑就开始卡了。

我的经验是：先用 0.15m 或 0.2m 的粗步长跑通全流程，确认形态没有问题后，再用 0.08m 左右的步长做局部细化。如果非要全区域细网格，可以分段计算，比如沿 y 方向切几片分别算完再拼起来。同时养成用 clear 清掉临时变量的习惯，特别是 rA、rB 这类用完就再也不用的距离数组，越早释放越好。

5.3 灵敏度归一化与绘图技巧

灵敏度值的动态范围非常大，源附近可以高出远处几十倍。直接绘图通常会看到一个刺眼的光点，而其他区域全是深色，掩盖掉真实形态。我的处理方法是先看量级再决定显示方式：

• 如果只是看形态，直接对灵敏度做 clip，把高于 max 的 20%~30% 压平，或者把低于 min 同样压平，这样正负区域的边界能看得更清楚；
• 如果需要对比不同配置的相对强度，不要用 log，因为灵敏度存在正负号，log 会把符号吃掉。可以改用 sign(Sen).*log10(1+abs(Sen)) 这种变换，保留正负的同时压缩动态范围；
• 三维可视化时，阈值不要取得太低，否则整个空间都会飘满等值面碎片；实际调试时先从最大值的 0.3 倍开始往下调，看哪个阈值能把目标区域的形态勾勒出来。

6. 常见问题与排查技巧实录

6.1 灵敏度出现负值，是不是代码写错了

不是。四电极装置（比如温纳）的灵敏度本来就存在负值区，这是电位差观测方式带来的必然结果。负灵敏度区域通常出现在测线端部外侧和较深部位，代表该处电阻率升高反而会导致观测视电阻率下降。代码排查时唯一要注意的是确认符号约定：如果 A、B 供电、M、N 取电位差，且公式采用“AB 伪电场点积 MN 伪电场再加负号”，那么中心浅部的正区和边缘负区的布局是稳定的。如果你调换 M、N 的顺序，图中正负区域会镜像反转，这不算错，但解释时要保持一致。

6.2 为什么跨井灵敏度在井壁附近特别大，中间反而变弱

这是点源场的几何衰减决定的。电流源附近的场强按 1/r² 增长，所以电线越贴近电极的区域灵敏度天然高。注意这不是“井壁处的探测效果好”，它只代表测量数据对井壁附近的电阻率变化极其敏感，而真正需要分辨的井间目标区，灵敏度是相对较低的。实际工程中，井壁附近常常还有泥浆侵入带、井孔套管等干扰体，这部分高灵敏度反而容易让反演把异常“吸”到井壁附近，形成伪分层。所以跨井反演通常要对井壁周围做额外的阻尼约束，或者直接挖掉井孔位置的模型单元。

6.3 坐标Z方向怎么统一：地表0还是向下正值

灵敏度点积公式在两种坐标定义下数值不变，因为点积不随坐标反射改变符号。但麻烦出在绘图和电极坐标录入时。我的建议是所有电极深度一律写成正数，z 网格从 0 开始向下递增，绘图用 axis xy 加 YDir reverse 把深度显示在下方向。如果混用“向上为正”的坐标，再加镜像法时特别容易把电极符号搞反，导致灵敏度分布整体翻转，看起来像目标从左边跑到右边了。

6.4 网格加密后灵敏度图反而出现放射状条纹

这通常不是公式错，而是由于源点处奇异性在细网格下被放大。理论上灵敏度核在电极附近趋于无穷，细网格会让这些极值点变成尖锐的刺。遇到这种情况，绘图前做一次中值滤波或者对灵敏度做小尺度高斯平滑即可。计算时完全不必要追求源点附近的精确值，因为这些区域的真实响应还被电极几何尺寸、接触电阻等因素主导，解析核只是一个近似。

6.5 三维计算太慢，有没有加速方案

有几个立竿见影的办法：一是去掉数组中间的重复计算，比如先算 rA.^3，再在后续出现时直接引用计算结果，避免反复开方；二是用单精度 single 存储临时距离数组，三维计算时内存占用直接减半，速度也会快不少；三是如果 MATLAB 版本支持，可以在三重网格上开启 gpuArray，把距离和点积计算挪到 GPU 上。对于 41×41×21 这种小网格，GPU 加速提升不明显，但对百万级网格就很可观了。

7. 我自己的使用套路与后续扩展

写到最后，分享一个我日常用得最多的组合：先在 2D 切面模式下跑完所有候选阵列，把不同极距、不同井距的灵敏度图截成对比图，这一步基本不花时间。选定一两个候选配置后，再上 3D 模式看井间覆盖的完整形态，重点检查是否有明显的覆盖空洞。等野外数据采完要反演时，我会把均匀介质下算出的灵敏度核改造成雅可比矩阵初值，再交给反演引擎去迭代。这样一整条流程下来，从观测设计到反演初始模型，全部在同一个 MATLAB 框架里完成，不用来回导格式。

顺带提一个扩展方向：这套代码稍微改一下网格和电极输入，就能批量生成整条测线每个测点的灵敏度叠加图，用来评估“哪些区域被多个测点覆盖、哪些区域只被极少数测点覆盖”。这种覆盖度热图比单组灵敏度图更有决策价值，适合野外施工前定测点间距。如果你想再进一步，可以把灵敏度核的输出直接做成矩阵存储，作为反演系统的模型分辨率矩阵，那又完全是另一个层面的用法了。