基于混合网格有限元的直流电阻率法三维正演研究

2022-06-07王新宇毛玉蓉严良俊高文龙

煤田地质与勘探 2022年5期

王新宇，王程，毛玉蓉，严良俊，周磊，高文龙

(1.长江大学油气资源与勘探技术教育部重点实验室，湖北武汉 430100；2.非常规油气省部共建协同创新中心，湖北武汉 430100；3.中煤科工集团西安研究院有限公司，陕西西安 710077)

直流电阻率法从地表到孔中、低维到高维、粗放到精细快速发展，在有色金属、煤田等矿产资源勘探和地质调查中应用广泛、效果明显[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 有限元边值问题及变分问题

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

式中：σ为地下介质电导率，S/m；σ′为异常电导率，σ′=σ−σ0，S/m；σ0为均匀介质电导率，S/m；us为异常电位，V，由地下的电性不均匀体产生；u0为正常电位，V，是点电流源在均匀半空间或全空间产生的电位，可由解析解求得；n为模型边界的外法向方向；r和r分别为边界上的点到点电流源的距离和有向距离，m；Γ0和Γ∞分别为计算区域 Ω的地面边界和无穷远边界。

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

式中：dΩ为体积分；dΓ为计算区域的无穷远边界面积分。

1.2 混合网格单元分析

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

图1 四面体和三棱柱的网格节点编号Fig.1 Nodes number of tetrahedral and triangular prism element

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

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

四面体的线性插值基函数(x,y,z)可表示为：

式中：i为四面体节点索引；Ve为四面体单元体积；ai、bi、ci、di为四面体线性插值基函数系数(参照文献[21]计算)。

正三棱柱的线性插值基函数[14]可表示为：

式中：i为三棱柱节点索引；Δ为三棱柱三角形面积；h为三棱柱单元高度(h=2s)；z0为三棱柱中心z坐标；ai、bi、ci(i=k,m,n,k′,m′,n′)为三棱柱线性插值基函数系数。

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

式中：Ne为混合网格单元总数；Ke=σ(K1e+K2e)，Ke为与 σ相关的单元系数矩阵；Pe=σ′(K1e+K2e)，Pe为与σ′相关的单元系数矩阵；K1e为式(2)中与体积分相关的单元系数矩阵的扩展矩阵；K2e为式(2)中与边界积分相关的单元系数矩阵的扩展矩阵；K为总电位的系数矩阵；P为异常电位的系数矩阵。

对于每个四面体单元，扩展矩阵K1e和K2e表示为：

式中：Δ1为四面体或三棱柱的三角形面积。

对于每个正三棱柱单元，扩展矩阵K1e表示为：

当三棱柱侧面为无穷远边界时，扩展矩阵K2e为：

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

1.3 线性方程组求解

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

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

1.4 视电阻率计算

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

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

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

图2 电阻率法井-地、地-井工作方式Fig.2 Surface-to-borehole and borehole-to-surface working mode of resistivity method

式中：ρs为视电阻率，Ω·m；I为电流大小，A；KAM和KAMN分别为二极装置与三极装置的装置系数；r和r′分别为点电源A和其关于地表镜像点到测点M或N的距离。

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

2.1 混合网格算法精度验证

采用一维层状地层模型(图3) 验证混合网格算法的正确性。计算区域大小为1 000 m×1 000 m×1 000 m，第一层地层的厚度为6 m，电阻率为50 Ω·m；第二层地层的厚度为10 m，电阻率为100Ω·m；基底电阻率为50 Ω·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 个等间隔测点。采用SSORPCG 算法求解线性方程组，混合网格求解耗时11.02 s，四面体网格求解耗时10.75 s；由解析解和数值解的对比结果(图5)可知：混合网格数值解的最大误差为0.15 %，四面体网格数值解的最大误差为1.56%(误差计算方法如下式)，2 种方法均满足计算精度。但混合网格的计算精度要略高于四面体网格，这是由于在较薄的地层离散网格时会产生质量较差的四面体网格，从而影响计算精度，这也证实本文混合网格算法的可行性、有效性。

图3 一维层状模型Fig.3 1D layered model

图4 层状模型的网格剖分Fig.4 Grid discretization of layered model

图5 不同方法的结果对比Fig.5 Comparison of forward modeling results of different methods

式中：ρa为混合网格或四面体网格三维数值解，Ω·m；ρb为一维解析解，Ω·m；E为数值解误差，%。

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

为研究钻孔因素对井-地电阻率法视电阻率响应的影响特征，设计钻孔模型(图6 为钻孔模型的网格剖分示意图)。钻孔孔径为10 cm，钻孔深度为50 m，计算区域大小为1 000 m×1 000 m×1 000 m，钻孔内存在电阻率为3 Ω·m的井液，围岩电阻率为100 Ω·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 钻孔模型的局部网格剖分Fig.6 Local grid of borehole models

图7 发射源位于不同深度的视电阻率响应Fig.7 Apparent resistivity response of point current source at different depth

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

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

图8 钻孔模型混合网格离散Fig.8 Mixed grid discretization of borehole model

表1 钻孔模型的地层参数Table 1 Formation parameters of borehole model

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

图9 不同钻孔参数下视电阻率响应曲线Fig.9 Apparent resistivity response curves of different borehole parameters

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

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