0 引言

五极纵轴电测深法是20世纪70年代提出的一种电探方法。黄兰珍等^[1]率先对该方法进行了系统的理论研究，分析了三个点电流源在地下均匀半空间的电流密度分布，给出了球体模型五极纵轴电测深曲线的解析计算公式，对多种情况下的测深曲线进行了计算分析，并给出了经验解释公式。目前，该方法已在水资源、矿产及工程勘查方面取得了较好的应用效果^[2-5]。然而，由于应用中存在供电排列是否对称、供电负极分流是否均匀、地表电性是否均匀及地形是否平坦等问题，这些观测环境与人为因素直接影响着五极纵轴电测深法的探测效果。因此，研究五极纵轴激电测深三维有限元数值模拟方法有望为后续模拟试验分析与实测数据反演奠定基础。

分析总结三维复杂地电模型的正演结果，对地球物理工作者正确认识野外实测的异常响应至关重要。Dey等^[6]、Scriba^[7]、Spitzer^[8]、吴小平等^[9]采用有限差分法完成了电阻率三维正演模拟；徐世浙等^[10]、Ma Qinzhong^[11]采用边界元法解决了三维地电断面的正演问题；Rucker等^[12]实现了直流电阻率的三维自适应有限元模拟；Zhou B等^[13]基于高斯积分网格的有限元法，实现了2.5维和三维各项异性介质的直流电阻率正演模拟；Plattner等^[14]基于自适应小波基方法实现了三维直流电阻率数值模拟；阮百尧等^[15]、刘海飞等^[16]采用有限元法实现了基于连续电性介质模型的三维激电数据正演模拟计算。 但上述研究工作仅限于单点源或双点源的供电方式，基于五极纵轴供电和接收方式的三维有限元数值模拟鲜有研究。

考虑到实际复杂地质情况以及模拟精度的要求，笔者基于连续电性介质模型及异常电位途径，对五极纵轴激电测深三维有限元数值模拟方法展开了研究。首先，简要地阐述五极纵轴电测深法的基本原理。然后，推导三个点电流源电位场满足的边值问题和变分问题，给出变分方程的有限元求解方法。最后，对几例地电模型进行模拟计算，检验有限元法的模拟精度和计算效率，给出野外供电极距的选择范围，对比不同装置对层状介质的分辨能力，以期为后续复杂模型的模拟试验分析奠定基础。

1 五极纵轴激电测深法工作原理

五极纵轴激电测深法的工作原理如图 1所示^[1]。在地面以测深点A为坐标原点建立直角坐标系，将垂直构造走向方向设为x轴，将平行构造走向方向设为y轴。沿x轴布设供电电极：坐标原点A为供电正极，供电电流强度为I；B₁和B₂为供电负极，供电电流强度均为－I/2。AB₁=AB₂=L，需大于2倍的勘探深度。工作时，供电电极A、B₁和B₂固定不动，测量电极MN=L/30～L/40，M、N沿y轴移动，测量总场电位差ΔV_MN=V_M－V_N(V_M、V_N分别为M、N的电位)。测深点位置记录在A点，视深度H_s=(Y₁+Y₂)/2。根据场的叠加原理，可整理出视电阻率的计算公式为

2 五极纵轴三点源电位场边值问题与变分问题 2.1 边值问题

当位于地表的稳恒点电流源A(+I)、B₁(－I/2)和B₂(－I/2)同时向地下供电时，在混合边界条件下，电位V的边值问题为

其中：∇为哈密顿算子；σ为地下电导率；I为供电电流强度；Ω为研究区域；δ(A)、δ(B₁)和δ(B₂)分别是以A、B₁和B₂为中心的狄拉克函数；Γ_s和Γ_∞分别为地面和无穷远边界；r_A、r_B1和r_B2分别为点源A、B₁和B₂至边界Γ_∞的距离；cos(r_A,n)、cos(r_B1,n)和cos(r_B2,n)分别为矢径r_A、r_B1和r_B2与边界外法线方向n的夹角余弦。

为消除场源奇异性的影响，将地下介质的电导率σ分解为场源处电导率σ₀与异常电导率σ′之和，即σ=σ₀+σ′，将电位V分解为正常电位U₀与异常电位U之和，即V=U₀+U，再将分解后的σ和V代入式(2)、(3)及(4)，便得到五极纵轴三点源异常电位的边值问题：

其中：

式中：r_AP、r_B1P和r_B2P分别为点源A、B₁和B₂到地下空间任一点P的距离。

2.2 变分问题

本文采用变分法推导五极纵轴三点源异常电位的变分问题。根据能量最小原理，对偏微分方程(5)构造泛函^[17]：

将式(8)两端对U求变分，得

根据场论中∇算子的运算规则

∇·(Aφ)=∇·Aφ+A·∇φ

及奥高公式

有

成立，其中：A和φ分别为任意矢量和标量；Γ=Γ_s+Γ_∞，表示研究区域Ω所包围的边界。则方程(9)变为

将方程(5)、(6)和(7)代入式(10)，得

即有

成立。这样便得到五极纵轴测深三点源异常电位边值问题对应的变分问题：

3 有限单元法

采用有限单元法求解五极纵轴测深异常电位的变分问题(12)，具体求解过程如下。

3.1 单元剖分

用六面体单元对区域Ω进行剖分(图 2a)，将式(12)中对区域Ω的积分分解为对各单元e的积分之和：

3.2 线性插值

在单元e内，电位和电导率均采用三线性插值^[15]：

式中：U_0i和U_i(i=1,2,…,8)分别为单元各节点的正常电位和异常电位；σ_i和σ′_i分别为单元各节点的电导率和异常电导率；n_i为形函数，且有

其中：ξ_i,η_i,ζ_i为点i的坐标；ξ,η,ζ与x,y,z的关系为

式中：x₀、y₀、z₀为子单元的中点；a、b、c为子单元的边长，如图 2b所示。

3.3 单元积分

首先，对式(13)右端的体单元进行积分得：

其中：k_1e=[k_1ij]；k′_1e=[k′_1ij]；U_e=[U_i]^T；U_0e=[U_0i]^T；i,j=1,2,…,8。

将式(14)、(15)代入式(16)，可得到体单元的刚度矩阵元素。将图 2c中8个网格节点的坐标代入式(19)中，便可以得到相应的系数。

接着，对式(13)右端的边界单元进行积分。假定单元e的一个面1234位于无穷远边界上，点源到边界面上任意一点的距离可视为相等，则可将式(13)中的D看作常数，提到积分号外。那么，边界单元Γ_e的积分可整理为

其中：k_2e=[k_2ij]；k′_2e=[k′_2ij]；U_e=[U_i]^T；U_0e=[U_0i]^T。

将式(14)、(15)代入式(20)，可得到边界单元的刚度矩阵元素。将图 2c中面的4个节点坐标代入式(23)便可得到相应的系数。

3.4 总体合成

在单元内，将式(16)和式(20)的积分结果相加，再扩展成由全体节点组成的矩阵和列阵：

其中：u与u₀分别是全体节点的u与u₀组成的列向量；k_e和k′_e分别是k_1e+k_2e和k′_1e+k′_2e的扩展矩阵。由全部单元的F_e(u)相加，得

令式(24)的变分为0，得到线性方程组

考虑到系数矩阵K具有大型、稀疏及对称等特点，故采用超松弛预条件共轭梯度法(SSORCG)求解线性方程(25)^[18-19]。将解方程得到的异常电位U与正常电位U₀相加即可得到总电位V=U+U₀，再结合式(1)便可计算出五极纵轴装置的视电阻率ρ_s。

当地下介质存在体极化特性时，其等效电导率为σ^*=σ(1－η)，η为介质的极化率。重复前述视电阻率的计算过程模拟出等效视电阻率_sρ^*，进而可计算出视极化率：

4 算例分析 4.1 精度与效率分析

为测试分析有限元方法的模拟精度，将文献^[20]中水平两层介质模型的单点源电位解析计算公式转化为五极纵轴激电测深的电位计算公式：

其中，V_M、V_N 、V_M^*和V_N^*分别为测量电极M和N处的电位及等效电位。结合视电阻率和视极化率的计算公式(1)和(27)—(30)，可计算出水平两层模型五极纵轴激电测深的理论曲线，以对比分析三维有限元数值解的模拟精度。

根据五极纵轴激电测深的三维有限元理论推导公式，我们基于C++开发了相应的计算程序。对如图 3a所示的两层水平层状模型模拟五极纵轴激电测深曲线，供电极距L为40 m，测量极距MN为2 m，点距为2 m，O是MN的中点，测点从2 m到76 m，应用不规则长方体单元将模型剖分成264 708(81×86×38)个未知节点。图 3b为有限元和解析方法模拟的视电阻率和视极化率曲线，可以看出数值解与解析解的曲线形态完全一致，并且两者拟合得很好；说明算法是正确的。而且从图 3c数值解与解析解的相对误差曲线可看出，视电阻率和视极化率的最大相对误差均控制在0.25%以内，满足模拟精度要求。模拟五极纵轴激电测深(包括视电阻率和视极化率)曲线需要两次正演计算，而且正演中合成刚度矩阵及右端项所需的时间较少(本例仅需1.3 s)，时间主要耗费在解大型稀疏线性方程组上，本文采用对称超松弛预条件共轭梯度法进行求解，终止误差为1×10^-10，两次正演在Intel双核CPU 2.67 GHz的微机上总耗费时间小于18 s，计算效率较高，如图 3d所示。

4.2 供电极距选择

假定地下均匀半空间中存在一大小为5 m×5 m×5 m的低阻异常体。围岩电阻率为100 Ω·m，极化率为1%；异常体电阻率为5 Ω·m，极化率为5%。异常体位于A极正下方11 m处(中心埋藏深度H=13.5 m)，供电极距L=10~100 m，测量极距MN为2 m，点距为2 m，测点从2 m到76 m。地电模型如图 4a所示。图 4b、c为模拟的视电阻率和视极化率曲线，可以看出：当L ＞2H时，视电阻率和视极化率测深曲线出现极值点，并且随着L的增大，测深曲线的极值特征越明显；当L＞5H时，曲线形态几乎趋于极限状态，视电阻率极值点位置S几乎无明显变化。出于经济和勘探效果的考虑，供电极距L选择(3~4)H较佳。

4.3 探测效果对比

假定地下为水平三层介质模型，从地表到深度8 m为第一层，电阻率为100 Ω·m ，极化率为2%；从深度10~14 m为第二层，电阻率为5 Ω·m ，极化率为5%；从深度16 m向下为第三层，电阻率为500 Ω·m，极化率为1%；相邻层之间电性为线性连续变化，具体如图 5中的折线所示。针对五极纵轴装置(FLS)、固定点源装置(FPS)及对称四极装置(VES)，供电极距(分别为L、AO和AB/2)和测量极距(MN或MN/2)如表 1所示。利用三维有限元数值模拟方法，模拟了三种装置的激电测深曲线，如图 5所示。从图 5中可看出，三种装置的曲线形态类似，特别是首支曲线基本完全重合，并且几乎在同一视深度位置出现低阻高极化异常；说明五极纵轴装置与另外两种装置具有同样的纵向分辨率，并都有视深度偏大的特点。对比三种观测装置的布极特点可知，五极纵轴装置的野外工作效率较高。

5 结论

1) 基于连续电性介质模型，推导了五极纵轴激电测深三个点电流源异常电位的边值问题、变分问题及有限元求解方法，编制了三维有限元模拟程序。通过对比水平两层介质模型五极纵轴激电测深的数值解和解析解，发现数值解的最大相对误差控制在0.25%以内，在CPU 2.67 GHz的PC机上模拟一条激电测深曲线的耗费时间为18 s。因此，从模拟精度和计算效率的角度来看，本文正演模拟算法是正确的，编制的模拟程序能够用于理论模型分析。

2) 从经济和勘探效果的角度考虑，五极纵轴激电测深方法在野外施工时，供电极距选择为3~4倍中心埋藏深度较佳。对于水平层状介质模型，五极纵轴装置与固定点源装置和对称四极装置的纵向分辨率基本相当，但野外布极方式相对经济。