基于三维数值模拟的激电法差分装置响应特征分析

凌嘉宣，陈志芳，邓维，戴世坤，周印明，陈轻蕊

（1. 桂林航天工业学院计算机科学与工程学院，广西桂林 541004；2. 浙江省钱塘江管理局勘测设计院，浙江杭州 310016；3. 中南大学地球科学与信息物理学院，湖南长沙 410083；4. 湖南科技大学计算机科学与工程学院，湖南湘潭 411201；5. 国防科技大学气象海洋学院，湖南长沙 410073）

0 引言

激发极化法（简称激电法）作为电法勘探重要的分支方法，通常情况下比电阻率法和大地电磁法更适用于目标体与围岩电阻率差异不大的浸染型金属矿床的勘探［1］。根据数据采集方式，激电法测量方式包括地—地、地—井、井—地、井—井［2］，其中后三种方法是利用已有钻孔探测井旁盲矿或矿体空间分布，确定矿体埋深、方位，追踪并圈定矿化带。这些方法勘探成本较高。地—地测量是在地面发射和接收，无需钻孔，成本较低，但是不适用于深部矿产的勘测［3］。不同测量方法各有优缺点，因此开展激发极化法探测应根据矿区（工区）实际情况和地质任务选择合适的测量方法。

正演是反演的基础，为了能更好地进行精细反演，需开展激电法三维数值模拟方法研究。激电法数值模拟与电阻率法数值模拟具有相通性，即激电法正演需要进行两次直流电阻率法的正演计算［4］。目前三维数值模拟方法主要包括积分方程法［5-8］、有限单元法［9-15］及有限差分法［15-18］，其他数值模拟方法均由这三种方法衍生得到。针对不同测量方法的激电三维数值模拟方法主要包括：Hohmann［19］和朴华荣等［20］基于积分方程法研究了地—地测量的三维极化体激发极化效应；张辉等［21］利用体积分方程法对同时存在激电效应和电磁效应的半空间三维体进行研究，并分析了激电效应和电磁效应并存对正演结果的影响；Tsourlos 等［22］和黄俊革［23］利用有限单元法开展了地—地三维激电数值模拟，验证了采用偏导数矩阵与视极化率的关系式计算极化率的可行性；吕玉增［24］和王智［25］利用基于四面体网格交叉剖分的有限单元法开展了地—井、井—地三维激电数值模拟，系统分析了井旁不同形态矿体的激电异常特征；周峰等［26］利用有限差分法开展了井—井激电三维正演，并利用二次异常电位差曲线反演得到矿体的激发极化特征；李长伟［27］利用放射状三棱柱单元结构网格的有限单元法实现了井—井激电法三维数值模拟；高文龙等［28］利用变步长有限差分法实现了井—井激电法三维数值模拟，并分析了侵入带、井径、不同电极系装置及不同井液电阻率对激电效应结果的影响；李静和等［29］为了高效、高精度探测重金属污染的激发极化空间分布，提出了接触式供电—地面观测和接触式供电—直接观测两种观测系统，进行了接触式激发极化法渗漏型目标探测的应用研究，解决了复杂施工环境条件的限制问题，同时提高了观测信号的强度和观测精度。

以上传统的激电法三维数值模拟大多基于空间域，在计算复杂模型时存在计算量大、存储需求高、效率较低等问题。本文基于二维傅里叶变换，将激电法三维数值模拟问题转化为波数域一维数值模拟问题进行求解，可以避免大型稀疏线性方程组的求解，降低计算量和存储需求，明显提高计算效率。 此外，传统激电测深装置主要为二极装置、三极装置、偶极装置等，这些装置对地下介质的分辨率有限。本文通过改善数据采集装置提高对地下电性异常体的分辨率，即引入差分装置［30］和积分装置。差分装置指的是以两组供电电流强度相等、方向相反的偶极发射装置并列或共线作为供电装置，利用偶极电极接收信号；积分装置以两组供电电流大小相等、方向相同的偶极发射装置并列或共线作为供电装置，也是利用偶极电极接收信号。本文基于井—井测量方式，对比传统测深装置、差分装置及积分装置对井旁极化矿体的分辨率。最后，对不同背景介质下二极装置和差分装置得到的视电阻率和视极化率响应特征进行分析，比较各自的优势。

1 理论方法

1.1 边值问题

激发极化法中，地下点电流源产生的电位满足微分方程［9］

式中：σ为介质电导率；U为总电位；I为电流强度；δ为狄拉克函数；ω为点电流源A到计算区域的立体角，当点电流源在地面时，ω=2π，当点电流源在地下时，ω=4π。总电位U可分解为背景电位Ub和异常电位Ua，即U=Ub+Ua；总电导率σ可分解为背景电导率σb和异常电导率σa，即σ=σb+σa。因此，式（1）可改写为

其中σb满足

将式（3）代入式（2），可得

因电位U与电场E满足E=-∇U，散射电流密度J=σbE，故式（4）改写为

对于水平层状介质模型，将式（5）沿三个方向展开，可得

式中J=［Jx，Jy，Jz］，Jx，Jy，Jz为空间域三个方向上的散射电流。

对式（6）沿水平方向进行二维傅里叶变换，可得

从式（7）可看出，方程右端Jx、Jy、Jz为未知，该式为的非线性方程，无法直接求出。因此，本文提出以下求解思路：①首先令总电场等于背景场Eb（n表示迭代次数，此时n=1），求出散射电流J；②对散射电流J进行二维正傅里叶变换，得到波数域散射电流的三分量并代入式（7），求得波数域异常电位；③利用波数域异常电位与异常电场的关系式求出波数域异常电场三分量，并对和进行二维反傅里叶变换，得到空间域异常电位Ua和异常电场分量Eax、Eay、Eaz，再与背景场相加可得总电位和总电场；④计算残差，若ε小于期望值，则输出结果，否则，根据迭代算子修改得到新的总电场，并利用总电场计算散射电流，重复步骤②～步骤④直至满足迭代结束条件。

根据前述求解激电法三维数值模拟问题的思路，需给出方程组的边界条件、求解方程组的方法、背景场、波数域电场与电位的关系式、波数选取方法、迭代算子等。

因式（7）沿x和y方向为波数域，由傅里叶变换特性可知，x和y方向的边界条件自动满足，所以只需给出z方向的边界条件。假设z轴向下为正方向，地面为上边界zmin，此处电流密度法向为零，满足

选取异常体下方某位置为下边界zmax，假设供电电极位于该区域之外，总电位U和背景电位Ub均满足拉普拉斯方程

因此

对应的波数域异常电位满足微分方程

式（12）的通解为

式中C和D为常数。因无穷远处（z→+∞）异常电位Ua等于零，有

对式（14）求导，下边界zmax上的异常电位Ua满足

式（7）、式（8）和式（15）即组成三维激电法的边值问题。

1.2 有限单元法

对于三维激电法的边值问题，可采用基于二次插值的一维有限单元法求解。利用伽辽金法［8］得到边值问题等价有限单元方程

式中：M为垂向剖分单元总数；Ni为第e个单元的第i个节点的二次插值形函数，这里i=1，2，3。式中各项的表达式见附录A。式（16）可表示为矩阵形式的线性方程组

从求解三维激电异常电位边值问题的思路可看出，本文利用二维傅里叶变换对复杂数值模拟问题进行降维求解，即将三维数值模拟问题转换成为不同波数域中的一维数值模拟问题进行求解，可大幅度降低计算量及对存储的需求，提高计算效率。

1.3 关键环节分析

1.3.1 背景场求解思路

本文算法基于异常场，因此需先确定背景场。在给出边值问题时，分别假设背景介质为均匀半空间、层状模型、电阻率随深度连续变化等不同情况。电阻率随深度连续变化的特殊情况就是层状介质，即垂向节点电阻率值连续变化。当背景为均匀半空间时，源在任意位置所产生的背景场的求解方法详见文献［1］。当背景为层状介质时，本文推导了点电流源在任意位置产生的背景场，详见附录B。

1.3.2 波数选取方法

本文采用基于标准FFT法的二维傅里叶变换，波数的选取基于采样定理，具体方法详见文献[32-33]。需要说明的是，二维反傅里叶变换可以采用基于标准FFT法或高斯-FFT法［33-34］。

1.3.3 迭代算子和波数域电场的求解

针对求解电位的边值问题，首先根据背景场求出电位和电场的近似解，然后利用算子进行迭代求解，最终获得高精度数值解。本文采用的迭代算子［35-36］为

1.3.4 三维数值模拟算法流程

本次三维数值模拟采用以下流程：

（1）对求解域采用均匀网格剖分，并利用采样定理得到波数域采样点。

（2）计算背景电场，并求出相应散射电流，利用二维傅里叶变换得到波数域散射电流。

（3）令边值问题式（7）等式右端的波数散射电流等于前一个步骤得到的值，求解方程组得到波数域异常电位，进而得到波数域异常电场。

（4）利用二维反傅里叶变换得到空间域异常电位和异常电场，加上背景场得到总场。

（5）将更新后的总电场与先前总电场进行对比，判断是否满足残差要求。若满足，则输出总电位和电场；否则，利用迭代算子修改总电场，并根据修改的总电场计算波数域散射电流，重复步骤（3）～步骤（5），直至满足残差要求。

2 算例分析

2.1 算法正确性验证

为了验证本文算法的正确性，利用均匀半空间背景下的低阻、高极化率球体模型进行测试。模型如图1所示。采用均匀网格剖分计算区域，节点总数为101×101×81，水平方向剖分间隔为1 m，垂直方向剖分间隔为0.5 m。算法利用Gauss-FFT法进行二维傅里叶反变换。采用二极装置进行测量，电流强度为1 A，供电电极置于点A（-10 m，0，0）。测点沿地面x方向布设。

图1 低阻、高极化球体模型示意图

采用解析解［1］作为理论值，两种算法的总电位和相对误差如图2 所示。可以看出，两种算法得到的总电位曲线吻合较好，最大相对误差为0.54%，验证了本文算法的正确性。

图2 图1 模型的电位解析解和数值解（左）及其相对误差（右）

2.2 差分装置激电异常响应特征

为了研究均匀半空间背景下激电差分装置对高极化率异常体的异常响应特征及不同背景条件下差分装置的异常响应特征，分别讨论以下两种情况。

2.2.1 高阻/低阻、高极化矿体模型

对于井—井勘探，为了探测井旁盲矿，为了节约成本，通常在单孔井中进行测量，即供电电极和测量电极均放置在同一口井中［2］。在金属矿产勘探中，黄铁矿、黄铜矿、方铅矿等矿产多表现为低阻、高极化特性，铬铁矿、黑钨矿及侵染状硫化矿等通常表现为高阻、高极化特性［37］。为了研究高阻或低阻的高极化矿体在差分装置下的激发极化响应特征，不失一般性，可利用简化的金属矿模型——两个大小相同的棱柱异常体模型进行测试。

模型如图3 所示，两个异常体中心点在地面上的投影位于直角坐标系原点，钻孔井口G坐标为（10 m，0，0），其他参数详见图中标注。分别利用五种装置（二极装置、三极装置、偶极装置、积分装置、差分装置）沿井孔测量。采用均匀网格剖分计算区域，节点总数为61×51×81，水平方向剖分间隔为2 m，垂直方向剖分间隔为1 m。另外，为了分析电极位于不同位置时盲矿的激电异常特征，分别将供电电极对A-B置于井口（中间点坐标：（10.0 m，0，1.5 m））、井中（中间点坐标（10 m，0，50 m））和井底（中间点坐标（10 m，0，90 m）），使用极距为1 m 的接收电极对M-N沿钻孔自上而下布设。

图3 两个棱柱异常体模型示意图

供电电极位于井口、井中和井底时得到的视电阻率和视极化率曲线如图4 所示。可以看出，供电电极放置于不同位置时，五种装置得到的极化体视电阻率和视极化率曲线形态大体相似，即矿体附近均出现异常。这种现象可依据“电偶极子”叠加原理进行分析，即矿体产生的异常场可近似等效为多个“电偶极子”组合叠加的结果［38］。此外，低阻异常体产生的异常远大于高阻异常体产生的异常，说明低阻异常体对电流更敏感。另外，从图中还看出，对于同一供电位置，差分装置测得的异常幅值远大于其他装置，积分装置和偶极装置产生的异常幅值次之，三极装置再次之，二极装置测得的异常幅值最小，说明通过物理差分可有效提高分辨率。因三极装置观测的是电位差，所以对介质的分辨率高于二极装置；积分装置和偶极装置相当于在三极装置的基础上对供电源进行了一次差分，所以它们的分辨率也高于三极装置。尽管积分装置和偶极装置的分辨率相同，但是积分发射装置是由两组偶极发射装置同向供电组合而成，所以供电电流更大，信噪比更高；差分装置则是在偶极装置的基础上对供电源再次进行差分计算，所以分辨率进一步提高。综上所述，差分装置勘探效果优于其他装置，可为野外勘探测量装置的选择提供参考。2.2.2 不同背景介质下的高极化矿体模型

图4 供电电极AB 位于不同位置时的电阻率（上）和极化率（下）曲线

为了进一步研究不同背景下矿体的激发极化特征，建立图5所示低阻、高极化异常体模型。异常体中心位置在地面的投影位于坐标原点。为了便于对比，对此模型分别利用二极装置和差分装置进行测量。点源（两组偶极源中心位置）位于（-20 m，0，0），接收电极M、N极距为1 m，沿钻孔向下布设。假设背景介质的极化率为0，电阻率分为三种情况：（a）均匀半空间，电阻率为1000 Ω ⋅m；（b）三层水平介质，自上而下各层介质电阻率为100、200、100 Ω·m，第一、第二层的厚度分别为10、20 m；（c）连续介质，电阻率随深度变化满足关系式计算区域大小为120 m×100 m×50 m，剖分网格总数均为61×51×51，水平方向剖分间隔均为2 m，垂直方向剖分间隔为1 m。三种不同背景下得到的电阻率和极化率见图6。

图5 低阻、高极化异常体模型示意图

图6 不同背景介质下的电阻率（上）和极化率（下）曲线

对比不同背景介质下两种装置的视电阻率曲线，可见差分装置能更好地反映背景介质和异常体的垂向电性特征和激发极化特征。

需要说明的是，不同背景介质下同一种装置得到的视极化率曲线形态相似，说明无激电效应的背景介质对极化率影响较小。在实际生产中，对复杂背景介质下极化矿体的勘查可综合考虑电阻率和极化率两个参数，以便更好地分析矿体及围岩的空间分布。

3 结论

（1）将空间域与波数域相结合能有效求解激电法三维数值模拟问题，为探索激电法数据的高效、高精度的数值模拟提供一条新思路。

（2）与其他装置相比，差分装置对异常体的视电阻率/视极化率异常响应更明显，分辨率更高。但因该装置是两组发射电流大小相同、方向相反的偶极发射装置组合而成，相对于其他装置，向地下供电能量更低，需增大供电电流，所以在实际应用中，需要综合考虑地质任务及经济成本，选择合适的测量装置及参数。

（3）视电阻率能有效反映目标体及围岩的电性特征，而视极化率能反映目标体的激电效应特征，综合分析这两个参数可更精确地解释异常体的空间分布及地电异常特征。

附录A 有限元单元分析

z轴方向可利用节点二次插值函数表示单元e的电位、电流密度及电导率，其形函数表达式为

式中j、p、m表示单元内节点编号。因此，式(16)的积分项为

式中

其中l为单元长度。

式(16)中的最后一项表示z=zmax处节点的值，此时

附录B 任意位置点源层状介质电位表达式推导

令笛卡尔坐标系原点位于地面，点电流源置于层状介质中任意位置，则源所在地层内测点和非源所在地层内测点的电位和电场分别为

式（B-1）和式（B-2）中：j表示源所在地层序号；i表示非源层序号；系数A和B的表达式根据点电流源在不同位置分为五种情况。假设地下包含n层介质，则点电流源位于第1 层、第2 层、第3～第n-2 层、第n-1层、第n层五种情况下A和B的取值分别见式（B-3）～式（B-7）。

式中

式中

式中

式中

式中

上面公式适用于地下介质至层数n>5 的情况。当n<5时，可将某层分解为电导率相等的多层，即可采用以上公式。若需计算多源下任意位置的电位和电场，可分别计算每个点电流源产生的场值，进行矢量叠加即可。求解点源在层状介质任意位置产生的电位和电场可采用滤波法。传统的滤波法需要逐个测点递推计算层系数，对于层状和节点较多的情况，计算速度慢。因此，本文采用优化的滤波法求解，即先计算电位和电场表达式中与z轴深度方向的相关系数，然后计算所有波数下与z轴深度方向无关的层系数和某一平面节点的电场／电位系数并进行存储，后续计算任意深度上测点的场值时，只需调用存储数据，将其与z轴相关系数相乘并进行累加即可，如此可快速计算背景场值。模型测试结果表明，采用此方法比优化前计算速度提高约20倍。