直流电阻率法从地表到孔中、低维到高维、粗放到精细快速发展，在有色金属、煤田等矿产资源勘探和地质调查中应用广泛、效果明显^[1-2]。直流电阻率法实测资料的合理解释离不开有效的三维正反演技术，并且直流电阻率法三维正反演一直是地球物理学界的研究热点^[3-8]。

国内外学者对电阻率法正演做了大量研究。20世纪70年代末期，国外学者A. Dey等^[9]将吸收边界条件的思想应用于电阻率法三维数值模拟，有效减少网格数目和边界效应的影响，取得良好的模拟效果；R. C. Fox等^[10]研究了起伏地形对电阻率法和激发极化法观测数据的影响；T. Lowry^[11]、ZHAO Shengkai^[12]等将电位分解成背景值和异常值的叠加，提出了异常电位法，克服场源的奇异性，有效地提高计算精度；C. Rücker等^[13]采用非结构四面体网格有限元法实现了复杂地形下电阻率法三维模拟；强建科等^[14]将三棱柱单元引入起伏地形条件下电阻率法有限元三维正演中，并研究了起伏地形对电阻率法正演响应的影响特征；汤井田等^[15]采用三棱柱网格，推导了无限元三维单元映射函数，取得较高的计算精度。在井–地电阻率法数值模拟方面，王志刚等^[16]基于积分方程法离散泊松方程，实现了井–地电阻率法三维并行正演算法；部分学者^[17-18]基于有限差分法实现了线源井–地电法三维正演，并取得了良好的数值模拟效果；李长伟等^[19]基于三棱柱网格实现了井–地电法有限元三维正演；王智等^[20]对井–地电阻率法三维正演的结果采用归一化总水平导数法，一定程度增强了异常体的响应特征。

上述研究在进行正演模拟时，均未考虑钻孔因素(井液电阻率和钻孔孔径)对视电阻率数据的影响。而考虑钻孔的地电模型，网格剖分势必变得困难。为解决上述问题，采用有限元法，结合混合网格的空间离散模式，解决小孔径条件下“空间离散困难”的问题。对于钻孔所在地层，采用三棱柱网格对钻孔进行快速离散；对于其他地层，采用非结构化四面体网格进行离散，以减少有限元的计算量。针对不同的地电模型，采用合适的网格类型进行离散，并开展高精度、快速的三维正演模拟。

1   混合网格有限元三维正演算法

1.1   有限元边值问题及变分问题

直流电阻率法有限元三维正演中，异常电位$ {u_s} $满足的边界问题^[21]为：

式中：$ \sigma $为地下介质电导率，S/m；$ \sigma ' $为异常电导率，$ \sigma ' = \sigma - {\sigma _0} $，S/m；$ {\sigma _0} $为均匀介质电导率，S/m；$ {u_{\rm{s}}} $为异常电位，V，由地下的电性不均匀体产生；$ {u_0} $为正常电位，V，是点电流源在均匀半空间或全空间产生的电位，可由解析解求得；$ {\boldsymbol{n}} $为模型边界的外法向方向；r和${\boldsymbol{r}}$分别为边界上的点到点电流源的距离和有向距离，m；$ {\varGamma _0} $和$ {\varGamma _\infty } $分别为计算区域$ \varOmega $的地面边界和无穷远边界。

则异常电位满足的变分问题^[21]为：

式中：$ {\text{d}}\varOmega $为体积分；$ {\text{d}}\varGamma $为计算区域的无穷远边界面积分。

1.2   混合网格单元分析

将计算区域$ \varOmega $离散为互不重叠的四面体和三棱柱单元，四面体和三棱柱节点编号规则如图1所示：

图  1  四面体和三棱柱的网格节点编号

Figure  1.  Nodes number of tetrahedral and triangular prism element

下载: 全尺寸图片 幻灯片

则第$e$个单元内任意一点的异常电位为：

式中：对于四面体单元$t = 4$，对于三棱柱单元$t = 6$；$ u_{\text{s}}^e(x,y,{\textit{z}}) $为第$e$个单元内任意一点的异常电位；$ u_{{{\text{s}}_i}}^e $为第$e$个单元上第$i$个节点的异常电位。

四面体的线性插值基函数$N_i^e(x,y,{\textit{z}})$可表示为：

式中：$i$为四面体节点索引；$ {V^e} $为四面体单元体积；$ {a_i} $、$ {b_i} $、$ {c_i} $、$ {d_i} $为四面体线性插值基函数系数(参照文献[21]计算)。

正三棱柱的线性插值基函数^[14]$ N_i^e(x,y,{\textit{z}}) $可表示为：

其中，${{a_k} = {x_m}{y_n} - {x_n}{y_m}} $，${{a_m} = {x_n}{y_k} - {x_k}{y_n}}，{{a_n} = {x_k}{y_m} }- $$ {x_m}{y_k} $，${{b_k} = {y_m} - {y_n}}，{{b_m} = {y_n} - {y_k}}$，${{b_n} = {y_k} - {y_m}}， {{c_k} = {x_m} }- $$ {x_n}$，${{c_m} = {x_n} - {x_k}}，{{c_n} = {x_k} - {x_m}} $，${{a_k} = {a_{k'}}}，{{a_m} = {a_{m'}}}，{{a_n} }= $$ {a_{n'}} $，$ {{b_k} = {b_{k'}}}，{{b_m} = {b_{m'}}}，{{b_n} = {b_{n'}}} $，$ {{c_k} = {c_{k'}}}，{{c_m} = {c_{m'}}}，{{c_n}} = $$ {c_{n'}} $。

式中：$i$为三棱柱节点索引；$\varDelta $为三棱柱三角形面积；$h$为三棱柱单元高度$(h = 2s)$；${{\textit{z}}_0}$为三棱柱中心${\textit{z}}$坐标；$ {a_i} $、$ {b_i} $、$ {c_i}(i = k,m,n,k',m',n') $为三棱柱线性插值基函数系数。

采用有限元法对式(2)进行离散后，对每个单元进行积分得到单元系数矩阵，最终形成总体系数矩阵：

式中：${N}e$为混合网格单元总数；${{\boldsymbol{K}}_e} = \sigma ({{\boldsymbol{K}}_{1e}} + {{\boldsymbol{K}}_{2e}})$，${{\boldsymbol{K}}_e}$为与$ \sigma $相关的单元系数矩阵；${{\boldsymbol{P}}_e} = \sigma '({{\boldsymbol{K}}_{1e}} + {{\boldsymbol{K}}_{2e}})$，${{\boldsymbol{P}}_e}$为与$ \sigma ' $相关的单元系数矩阵；${{\boldsymbol{K}}_{1e}}$为式(2)中与体积分相关的单元系数矩阵的扩展矩阵；${{\boldsymbol{K}}_{2e}}$为式(2)中与边界积分相关的单元系数矩阵的扩展矩阵；$ {\boldsymbol{K}} $为总电位的系数矩阵；$ {\boldsymbol{P}} $为异常电位的系数矩阵。

对于每个四面体单元，扩展矩阵${{\boldsymbol{K}}_{1e}}$和${{\boldsymbol{K}}_{2e}}$表示为：

式中：$ \varDelta { _1} $为四面体或三棱柱的三角形面积。

对于每个正三棱柱单元，扩展矩阵${{\boldsymbol{K}}_{1e}}$表示为：

当三棱柱侧面为无穷远边界时，扩展矩阵${{\boldsymbol{K}}_{2e}}$为：

式中：$\varDelta{ _2} $为三棱柱外边界四边形的面积。当三棱柱侧面(四边形)位于无穷远边界时，$ {{\boldsymbol{K}}_{2e}} $积分采用式(13)计算；当三棱柱底面和四面体表面(三角形)位于无穷远边界时，$ {{\boldsymbol{K}}_{2e}} $积分采用式(8)计算。至此，四面体与三棱柱单元积分和边界积分完成，由于采用标量有限元，四面体与三棱柱共节点处的积分点重合；如图1所示，当四面体单元2、3、4节点与三棱柱单元$ k' $、$ m' $、$ n' $重合时，积分点重合的单元自动满足标量场的一致性。令式(6)的变分为零，得到$ {\boldsymbol{Ax}} = {\boldsymbol{b}} $形式的线性方程组：

式中：$ {\boldsymbol{A}} = {\boldsymbol{K}} $；$ {\boldsymbol{x}} = {{\boldsymbol{u}}_{\text{s}}} $；$ {\boldsymbol{b}} = - {\boldsymbol{P}}{{\boldsymbol{u}}_{\text{0}}} $。

1.3   线性方程组求解

本文采用改进的对称逐步超松弛预处理共轭梯度迭代算法(SSOR-PCG)^[22]求解式(14)。引入对称逐步超松弛迭代(SSOR预条件)的分裂矩阵作为预处理矩阵：

式中：${\boldsymbol{H}}$为预处理矩阵；${\boldsymbol{D}}$为${\boldsymbol{A}}$的对角阵；${\boldsymbol{L}}$为${\boldsymbol{A}}$的严格下三角矩阵；$ 0 < \omega < 2 $为松弛因子(算法细节参照文献[22])。

1.4   视电阻率计算

在求解方程后得到各节点的异常电位值，从而得到各节点的总电位值：

式中：$ u $为每个节点的总电位。

本文主要研究井–地、地–井工作方式的视电阻率响应特征，图2为井–地、地–井工作方式的示意图。二极装置与三极装置的视电阻率计算方式分别如下：

式中：${\rho _{\rm{s}}}$为视电阻率，$\Omega \cdot {\rm{m}}$；$ I $为电流大小，A；${K_{AM}}$和${K_{AMN}}$分别为二极装置与三极装置的装置系数；$r$和$r'$分别为点电源$A$和其关于地表镜像点到测点$M$或$N$的距离。

图  2  电阻率法井–地、地–井工作方式

Figure  2.  Surface-to-borehole and borehole-to-surface working mode of resistivity method

下载: 全尺寸图片 幻灯片

2   混合网格有限元电阻率法数值模拟分析

2.1   混合网格算法精度验证

采用一维层状地层模型(图3) 验证混合网格算法的正确性。计算区域大小为1 000 m×1 000 m×1 000 m，第一层地层的厚度为6 m，电阻率为$50{\text{ }}\Omega \cdot {\rm{m}}$；第二层地层的厚度为10 m，电阻率为$100{\text{ }}\Omega \cdot {\rm{m}}$；基底电阻率为$50{\text{ }}\Omega \cdot {\rm{m}}$。图4a为混合网格剖分示意图，顶部两层地层均采用三棱柱网格剖分，基底地层采用四面体网格剖分，最终生成170 360个三棱柱，294 820个四面体以及57 948个边界三角形和3 200个边界四边形。图4b为非结构四面体网格剖分示意图，最终生成473 222个四面体，10 348个边界三角形。发射源置于地表(0 m, 0 m, 0 m)处，发射电流为1 A，在x方向1~100 m的范围内布置100个等间隔测点。采用SSOR-PCG算法求解线性方程组，混合网格求解耗时11.02 s，四面体网格求解耗时10.75 s；由解析解和数值解的对比结果(图5)可知：混合网格数值解的最大误差为0.15 %，四面体网格数值解的最大误差为1.56%(误差计算方法如下式)，2种方法均满足计算精度。但混合网格的计算精度要略高于四面体网格，这是由于在较薄的地层离散网格时会产生质量较差的四面体网格，从而影响计算精度，这也证实本文混合网格算法的可行性、有效性。

图  3  一维层状模型

Figure  3.  1D layered model

下载: 全尺寸图片 幻灯片

图  4  层状模型的网格剖分

Figure  4.  Grid discretization of layered model

下载: 全尺寸图片 幻灯片

图  5  不同方法的结果对比

Figure  5.  Comparison of forward modeling results of different methods

下载: 全尺寸图片 幻灯片

式中：${\rho _{\rm{a}}}$为混合网格或四面体网格三维数值解，$\Omega \cdot {\rm{m}}$；${\rho _{\rm{b}}}$为一维解析解，$\Omega \cdot {\rm{m}}$；E为数值解误差，%。

2.2   钻孔模型井–地观测方式数值模拟分析

为研究钻孔因素对井–地电阻率法视电阻率响应的影响特征，设计钻孔模型(图6为钻孔模型的网格剖分示意图)。钻孔孔径为10 cm，钻孔深度为50 m，计算区域大小为1 000 m×1 000 m×1 000 m，钻孔内存在电阻率为$3{\text{ }}\Omega \cdot {\rm{m}}$的井液，围岩电阻率为$100{\text{ }}\Omega \cdot {\rm{m}}$。在钻孔周围采用三棱柱网格进行离散，井底下方采用四面体网格进行离散，混合网格最终生成265 145个节点，497 200个三棱柱，76 587个四面体以及55 496个三角形和6 160个四边形。发射源分别位于井口(0 m, 0 m, 0 m)，井中(0 m, 0 m, 5 m)、(0 m, 0 m, 10 m)、(0 m, 0 m, 15 m)处，发射电流均为1 A，采用井–地观测方式计算地表测点的视电阻率响应。数值模拟结果如图7所示，将视电阻率的色标限制在围岩电阻率的上下5%范围内，超出此范围，即可认为钻孔因素对视电阻率数据的影响不可忽略。由图7a可以看出，当发射源位于地表时，钻孔附近视电阻率受钻孔因素影响最大，但随着测点收发距的增加，视电阻率响应逐渐趋于围岩电阻率。如图7b所示，当发射源位于地下5 m时，地表视电阻率受钻孔因素的影响减小，但仍然不能忽略。当发射源位于地下10 m(图7c)，15 m(图7d)时，地表视电阻率响应和围岩电阻率非常接近，可以忽略钻孔因素对视电阻率数据的影响。因此，在实际井–地电阻率法资料处理解释中，根据发射源位置，合理考虑钻孔因素的影响是必要的。

图  6  钻孔模型的局部网格剖分

Figure  6.  Local grid of borehole models

下载: 全尺寸图片 幻灯片

图  7  发射源位于不同深度的视电阻率响应

Figure  7.  Apparent resistivity response of point current source at different depth

下载: 全尺寸图片 幻灯片

2.3   钻孔模型地–井观测方式数值模拟分析

为研究钻孔因素(井液电阻率和钻孔孔径)对地–井观测方式视电阻率响应曲线的影响，设计钻孔地电模型(图8)。钻孔深度为100 m，围岩采用层状地层模型，该模型的地层厚度及电阻率参数见表1，异常体的大小为30 m×20 m×20 m，电阻率为$10{\text{ }}\Omega \cdot {\rm{m}}$。发射源位于地表井口(0 m, 0 m, 0 m)处，发射电流为1 A，接收电极位于井中2~100 m，电极间距2 m，以三极装置方式进行接收。分别设计如下2个模型：① 钻孔孔径为50 mm时，计算不同井液电阻率(3、10、$50{\text{ }}\Omega \cdot {\text{m}}$)情况下的视电阻率响应曲线；② 井液电阻率为$3{\text{ }}\Omega \cdot {\rm{m}}$时，计算不同钻孔孔径(0、50 、100 、150 mm)情况下的视电阻率响应曲线。

图  8  钻孔模型混合网格离散

Figure  8.  Mixed grid discretization of borehole model

下载: 全尺寸图片 幻灯片

表  1  钻孔模型的地层参数

Table  1.  Formation parameters of borehole model

地层编号	地层厚度l/m	地层电阻率ρ/(Ω·m)
1	20	100
2	10	50
3	40	200
4	10	50
5	20	100
6	∞	150

下载: 导出CSV 

| 显示表格

图9a展示了地–井观测方式在不同井液电阻率情况下的视电阻率响应曲线。由图9a可知，浅部测点的视电阻率数据受井液电阻率影响较大，深部测点的视电阻率数据基本不受井液电阻率影响；而且随着井液电阻率的增大，视电阻率曲线所受影响逐渐减小。实际钻孔中的井液电阻率均为低阻，通过分析井液电阻率为$3{\text{ }}\Omega \cdot {\rm{m}}$时的视电阻率响应可知：当测点收发距为10 m左右时，井液的影响才可以忽略不计。因此，在实际地–井电阻率法勘探中，钻孔因素应加以重视。图9a中的视电阻率曲线很难反映地层的电性特性和界面信息，这是由于地下介质存在明显的电阻率差异，高阻地层对电流具有排斥作用、低阻异常和低阻地层对电流具有吸引作用，地层与异常体之间的电流场相互耦合，使得视电阻率响应变得尤为复杂，很难对地下介质的电阻率信息进行有效判别。因此，在复杂地层环境下，地–井观测方式的视电阻率曲线变化较为复杂，给直流电阻率法资料的处理解释工作带来巨大困难。

图  9  不同钻孔参数下视电阻率响应曲线

Figure  9.  Apparent resistivity response curves of different borehole parameters

下载: 全尺寸图片 幻灯片

现将所有地层的电阻率设为$100{\text{ }}\Omega \cdot {\rm{m}}$，异常体的电阻率为$10{\text{ }}\Omega \cdot {\rm{m}}$，在忽略钻孔的情况下，研究地下低阻异常体对视电阻率曲线的影响特征。图9b展示了地–井观测方式的视电阻率响应曲线，可以看出，由于低阻异常体对电流的吸引作用，低阻异常体与围岩之间的电流场相互耦合，导致视电阻率响应在一定区域内出现大于围岩电阻率的假异常。

图9c展示了地–井观测方式在不同钻孔孔径情况下的视电阻率响应曲线。相比于井液电阻率，视电阻率响应更易受钻孔孔径大小影响。且随着钻孔孔径的增大，视电阻率曲线在浅部畸变更为严重。而图9a和图9c在200 m以下视电阻率曲线基本不受钻孔因素影响，在实际电阻率资料解释中收发距较小测点的视电阻率响应需要加以重视，收发距较大的测点在资料信噪比较高的情况下可以直接进行处理解释。而对于电阻率法测井，发射源与测点均在井中，可以直接借助本文算法确定最佳收发距，以消除钻孔参数的影响，从理论上确保测井数据的有效性。

3   结 论

a. 提出基于四面体和三棱柱混合网格有限元的直流电阻率法三维正演算法。通过一维层状模型的数值验证，混合网格算法可以在保证网格质量的同时，提高计算精度。

b. 为有效处理钻孔网格剖分，采用混合网格算法对井–地、地–井方式下的地电模型进行正演模拟。数值结果表明，井液电阻率或钻孔孔径的变化会对收发距较小测点处的视电阻率响应产生严重的影响，在实际井–地、地–井电阻率法勘探过程中，钻孔参数应加以重视。

c. 对于地–井观测方式，由于测点位于地下，复杂的地层信息和钻孔参数会导致地–井观测方式的视电阻率响应出现假异常。

d. 在实际勘探中，可利用本文算法，结合钻孔信息进行正演模拟分析，选择合适的收发距，以保证实际直流电阻率资料的合理解释。为更好地解决以上因素带来的问题，下一步将展开基于混合网格下的钻孔参数约束的直流电法三维反演工作，为直流电阻率法精细化解释提供理论基础。