三维井地电阻率法异常响应特征增强算法模拟研究

王 智,潘和平

(中国地质大学(武汉)地球物理与空间信息学院,湖北武汉430074)

井地电阻率法是在井中供电地面接收电磁场的一类电法勘探方法,最早由前苏联科学家于1958年提出[1],用来圈定煤层或油气藏水平边界。通常情况下，井地电阻率法的结果以视电阻率等值线图或异常电位等值线图表示。国内外对井地电阻率法的正演数值模拟做了很多的研究。Scriba[2]和Spitzer[3]利用边界元实现了三维地电断面的正演问题;Pridmore等[4]用有限元实现了三维正演模拟的计算;徐凯军[5]和王志刚等[6]利用有限差分完成了垂直线源井地电位法的三维正演模拟;柯敢攀等[7]利用有限元实现了垂直线源井地电位法的三维数值模拟;戴前伟等[8]利用有限元法完成了2.5维的井地电位法的正演问题;贾正元等[9]利用有限差分实现了三维井地电阻率法的数值模拟,并对点源、线源、倾斜源三种情况下产生的异常特征进行了详细的总结;刘海飞等[10]采用有限元实现了三维连续电性介质的数值模拟。

在井地电阻率圈定异常体边界的研究方面,张天伦[11-12]用物理水槽实验分别进行在井中无套管和井中有套管情况下确定油气藏边界的研究,证实了利用井地电阻率法确定异常体边界的可行性;并在后来提出了采用非无穷远处三极剖面法确定油气藏边界[13],且应用此方法在新疆等地成功发现了小块油气藏[14-15]。汤井田等[16]利用井地电阻率法歧离率确定高阻油藏边界;戴前伟等[17]利用梯度算子与拉普拉斯算子对异常响应进行增强,取得了一定的效果。戴前伟等[8]利用不同点源深度时的视电阻率异常断面图的中心深度位置来确定异常体的上、下边界;汤井田等[16]利用歧离率确定高阻油气藏的边界取得一定效果,但只限于钻孔刚好打到油气藏上,因此歧离率法发现井旁盲矿的能力很弱。

有时井地电阻率法的异常响应信号非常弱小,特征不明显，特别是当异常体的埋深较深时，这种现象更加突出。因此，在总结前人研究成果的基础上，我们从井地电阻率法满足的边值问题出发，采用点电流源，井中供电、地表接收的观测方式，利用有限元法实现了三维井地电阻率法的数值模拟；计算了多种不同情况下三维异常体模型的异常电位水平导数，分析了异常电位水平导数的极值与异常体边界对应关系的影响因素；引入归一化总水平导数法[18]，以研究提高深部异常响应特征和增强异常响应的边界范围。

1 三维井地电阻率有限元数值模拟算法

1.1 基本原理

在稳恒电流源供电下的研究区域中，各点电位满足如下的边值问题:

(1)

与(1)式对应的变分问题为

(2)

式中:A是与源项有关的坐标向量;σ为介质的电导率;ωA是地面源点对地下区域的立体角,若源点在地下ωA=4π,若源点在地面ωA=2π;I为电流源项;u为总电位;Γs为地面边界;Γ∞为无穷远边界;Ω为Γs与Γ∞组成的封闭区域;r为源点到测点的距离;cos(r,n)是矢径r与外法向方向n的夹角余弦。

1.2 有限单元法

本文利用有限单元法求解上述变分问题,分为4个步骤。

1.2.1 网格剖分

采用四面体剖分,如图1所示。

图1 四面体单元

1.2.2 线性插值

设ui(i=1,2,3,4)是单元中4个节点的电位值,则四面体单元e内任一点p(x,y,z)电位用这4个角点的电位进行线性插值近似得到

(3)

式中：Ni(i=1,2,3,4)是形函数,它是x,y,z的线性函数:

(4)

式中：V是四面体的体积;Vi是插值点p(x,y,z)与四面体其它3个角点所组成的四面体体积;ai,bi,ci,di(i=1,2,3,4)是与四面体单元顶点坐标有关的常数。

1.2.3 单元积分

变分问题(2)式中，第1项对应的单元积分为

(5)

其中,

(6)

第2项为电源项,只与电源点的位置有关;第3项为边界单元项,若单元的一个面的1,2,3节点落在边界上,则

(7)

(8)

式中:Δ为单元e在边界上的三角形的面积。

1.2.4 方程组的求解

将单元积分结果扩展成由全体节点组成的矩阵后,进行全部单元相加得到如下的线性方程组:

(9)

式中：K为刚度系数矩阵;U为待求电位矩;P为与电源有关的矩阵。

目前解决这类方程的求解方法有很多,直接法有GS法、奇异值分解、LDLT分解等;迭代法有Newton法、共轭梯度法及其改进方法等。其中吴小平等[19]与柳建新等[20]分别用对称超松弛预条件共轭梯度法(SSOR-PCG)与不完全LU分解预处理的稳定双共轭梯度算法(ILU-BICGSTAB)实现了快速稳定的计算。结合前人的研究,我们采用不完全LU分解预处理的稳定双共轭梯度算法求解线性方程组(9)[20],可得到各个节点的电位。为了消除源点附近的奇异性,我们采取异常电位法计算,对于异常电位的计算过程与总电位计算相似,最后归结为求解大型方程组得到的异常电位,异常电位和正常电位相加得到总电位。

1.3 算法正确性的验证

基于上述理论推导过程,编制了相应的程序。以下用两个模型验证程序的正确性和有效性。

模型一为3层水平层状模型(图2),每层参数

为:第1层电阻率50Ω·m,厚度为10m;第2层电阻率为100Ω·m,厚度为10m;第3层电阻率为10Ω·m。图2坐标轴上A为源点,M为测点。对该模型进行正演计算,并将计算结果与数值滤波法计算的结果进行对比,如图3所示。从图中看出二者吻合得非常好,最大误差为0.03%。

图2 水平层状模型

图3 水平层状模型数值计算与滤波法结果对比

模型二为地下低阻球体模型(图4),源点和测点均在地面,模型参数为:A为源点,坐标为(50m,50m,0)M为测点;围岩电阻率100Ω·m,球体电阻率为10Ω·m,球体半径为5m。球体中心埋深10m,坐标为(50m,50m,10m),观测电极距2m。对该模型进行正演计算,将计算结果与均匀半空间内球体外的电位解析解结果进行对比(图5),从图5中看出二者吻合得非常好,最大误差为2.85%。

图4 低阻球体模型a 三维模型; b 模型XOY平面

图5 低阻球体模型数值计算与解析法结果对比

2 水平导数极值圈定异常体范围的影响因素分析

对于边缘识别方法在重磁位场中的应用研究较多[21],其中数值类边缘识别方法均是利用极大值位置或者零值位置来确定地质体的边缘位置,其理论基础是二度体铅垂台阶模型的重力异常特征,在该模型边缘处重力异常总水平导数和解析信号振幅达到极大值，而垂向导数达到零值,故可利用这些特征位置来确定二度体铅垂台阶的边缘位置[22]。在电法勘探中,通常在异常体的边界范围处,即异常体与围岩横向接触处,地表的电流密度变化同样很大[16]。由此为利用水平导数确定其边界位置提供了可能性。为了解水平导数极值确定异常体边界范围的有效性,下面从源的埋深、异常体的埋深、围岩与异常体的差异、异常体的横向位置对异常电位的水平导数的极值影响进行分析。总水平导数的定义如下[23]:

THD(x,y,0)=

(10)

式中:THD表示原始异常f的总水平导数;f(x,y,0)为某一平面上的数据体,在本文中为地面的异常电位数据。

2.1 总水平导数的求取方式

我们给出两种求取水平导数的方式。

方法一:从推导单元分析的积分方程过程中可得到异常电位法水平导数的解析求法。由式(3)可得:

(11)

其中，

(12)

同理可得

(13)

由式(12)与式(13)可知,某一点处的异常电位水平导数只与该单元内各个网格节点的坐标和异常电位有关。所以，在单元内电位u的导数是常数。相邻单元的∂u/∂x,∂u/∂y是不连续的。相邻单元的u通过公共边界点上的u连接起来,u在全区域内是连续的,但导数是不连续的。对于相邻单元的公共节点,可用各单元求得的导数的平均值,作为该点的导数值[24]。在正演的过程中需要计算ai,bi,ci,di(i=1,2,3,4)的值,因此在程序中保留原有的值,正演结束后,将地面四面体单元体内的各个节点的异常电位回代到方程(12)与(13)中,即可得到该单元内某一点处的异常电位的水平导数。

方法二:采用中心差分方法代替微积分进行近似计算。详细的计算过程请参见文献[17]。

地下低阻异常体模型的总水平导数试算结果表明，上述两种方法计算的结果是一致的，由此验证了方法一的正确性。

2.2 源的埋深对异常电位X方向水平导数的影响

图6为一个地下低阻异常体的模型,模型的参数如下:A1(100m,100m,30m),A2(100m,100m,50m),A3(100m,100m,100m),A4(100m,100m,130m)为4个不同埋深的源点,M为测点,B为井口位置。井口在模型的中心位置,坐标为(100m,100m,0),设围岩电阻率ρ0=100Ω·m,异常体电阻率ρn=10Ω·m。异常体顶面埋深5m,距井距离为10m,长、宽、高均为20m,正方异常体中心坐标为(80m,80m,15m)。

图6 不同点源深度的低阻正方体模型a 三维模型; b 模型XOY平面

图7是异常电位X方向水平导数随点源埋深的变化规律,从图中可以看到,在异常体的两个边界处水平导数的极值位置基本不随点源的埋深变化,极值之间的宽度保持不变,水平导数的幅值随点源的埋深增加而减小。由此可见，点源的埋深对异常电位的水平导数的极值相对位置几乎没有影响,即点源的埋深对异常体的水平范围的确定没有影响。

图7 在y=80m处不同深度点源计算的X方向水平导数结果对比

2.3 围岩与异常体的电阻率差异对异常电位X方向水平导数的影响

采用图6的低阻异常体模型分别计算不同电阻率下的X方向水平导数,其中源点坐标为(100m,100m,50m)。设围岩的电阻率ρ0=100Ω·m,各异常体的电阻率ρn分别等于5,10,20,50,90,200,500Ω·m。图8为不同电阻率差异下的X方向水平导数结果对比。

从图8中可以看到,在异常体的两个边界处水平导数的极值位置基本不随围岩与异常体电阻率差异不同而变化,极值之间的宽度保持不变,水平导数的幅值随两者的电阻率差异减小而变小。由此可见，围岩与异常体的电阻率差异对异常电位的水平导数的极值相对位置几乎没有影响,即围岩与异常体的电阻率差异对异常体的水平范围的确定没有影响。

图8 在y=80m处不同电阻率差异下X方向水平导数结果对比

2.4 偏离钻井不同位置时对异常电位X方向水平导数的影响

基于图6的低阻异常体模型,分别计算异常体右边界面距离井5,10,15,20,25m时异常电位X方向的水平导数,其中源点坐标为(100m,100m,50m)。设围岩电阻率ρ0=100Ω·m,异常体电阻率ρn=10Ω·m,图9为偏离钻井不同距离时计算X方向水平导数结果对比。

由图9可以看出,在异常体的两个边界处水平导数的极值位置随异常体偏离钻井位置变化而变化,而极值之间的宽度保持不变,水平导数的幅值随着偏离钻井的距离增大而变小。由此可见，异常体偏离钻井的距离对异常电位的水平导数的极值相对位置几乎没有影响，即异常体偏离钻井的距离对异常体的水平范围的确定没有影响。

2.5 不同埋深的异常体对异常电位X方向水平导数的影响

基于图6的低阻异常体模型，分别计算异常体顶面埋深为1,5,15,20,50m处的异常电位X方向的水平导数。设围岩电阻率ρ0=100Ω·m,异常体电阻率ρn=10Ω·m,图10为不同埋深下的异常体的X方向水平导数结果对比,图11为水平导数的极值宽度D随异常体埋深Z的变化规律。

从图10与图11 中可以看到,随着异常体的埋深增加,异常体X方向的水平导数的极值位置向外移动且幅值减小,异常电位的水平导数的极值之间的宽度随埋深增加而增大。这遵从电阻率法的基本规律,异常体的埋深越浅,异常体的两个边界与极值的位置对应越精确。因此,异常体的埋深对异常电位的水平导数的影响较大。

通过大量的模型实验研究,得出如下结论:源的埋深、围岩与异常体的差异、异常体的横向位置对异常电位的水平导数的极值位置几乎没有影响,而异常体的埋深对极值的位置,以及极值之间的宽度产生较大的影响,随着异常体埋深的增加,异常体的水平导数的极值位置向外移动且幅值减小。当异常体的埋深较浅时,异常体的水平导数的极值位置与异常体的边界有着很好的对应关系。

图10 在y=80m处不同埋深的异常体计算的X方向水平导数结果对比

图11 水平导数的极值宽度随异常体埋深的变化

3 归一化总水平导数的定义及应用

当存在两个不同埋深的异常体时,往往深部的异常体的异常响应很弱小,特征不明显,为了使深部与浅部的异常体的响应特征及边界范围都能得到有效的显示,引入归一化水平总导数进行异常体的特征识别。

归一化总水平导数法是利用计算点总水平导数与一定窗口范围内总水平导数最大值的比值进行地质体边界的识别,其具体的表达式为[18]

(14)

式中:THD表示原始异常u的总水平导数;NTHD(i,j)表示点(i,j)处归一化总水平导数值;m,n分别表示窗口的尺寸。

为了对归一化总水平导数对异常体的边界识别的有效性进行检验,设计了2个双异常体模型,源点A的坐标为(100m,100m,80m),井口B的坐标为(100m,100m,0),2个异常体的电阻率均为10Ω·m,围岩的电阻率为100Ω·m。异常体的大小均为20m×20m×20m。①第1组模型:第1个异常体和第2个异常体的顶面埋深均为1m;②第2组模型:第1个异常体顶面埋深为1m,第2个异常体的顶面埋深为30m。各项模型参数如图12与图13所示。

图14与图15分别为两个模型的异常电位平面图。当埋深较浅时,异常电位的梯度带都能显示出异常体的范围,但是异常的边界范围较为模糊；当存在深部异常体时,其异常电位的响应很弱小,梯度带变得更加稀疏,边界范围已无法准确地识别出来。图16中,红色虚线框代表真实异常体的边界在地面的投影,红色实线框为本方法识别出的异常体的边界范围。由图16a与图16b可以看出,当异常体的埋深均较浅时,总水平导数法与归一化总水平导数均能精确地反映出异常体的水平边界及范围,异常体的边界得到了显著的增强,这对异常的响应范围的圈定很有利;随着深度的增加,导数的幅值迅速衰减,总水平导数法已无法识别出深部异常体,如图16c所示;图16d为异常电位的归一化总水平导数法识别结果,可以看到,该法能清晰地识别出不同深度的异常体,提高了异常体的识别能力;由于异常体埋深的增加,导致了极值的位置与异常体实际的边界位置有一定差异,这与前面分析得出的结论一致。

图12 第1组双异常体模型示意a XOY平面; b XOZ垂直断面

图13 第2组双异常体模型示意a XOY平面; b XOZ垂直断面

图14 第1组模型的异常电位平面等值线

图15 第2组模型的异常电位平面等值线

图16 双异常体模型总水平导数与归一化总水平导数的识别结果a 第1组模型的总水平导数; b 第1组模型归一化总水平导数; c 第2组模型的总水平导数; d 第2组模型归一化总水平导数

4 结束语

1) 从点源场满足的基本关系式出发,通过有限元实现了三维井地电阻率法正演计算,设计了两个理论模型计算出解析解,证明了程序的正确性；并给出了异常电位水平导数的2种计算方法。

2) 通过大量的模型试算,详细分析了水平导数极值与异常体边界对应关系的影响因素,异常电位的水平导数极值不受源的埋深、异常体偏离钻井的位置、异常体与围岩的电阻率差异的影响;异常电位的水平导数幅值会随着异常体的埋深增加而减小，且极值间的宽度会逐渐变宽。

3) 当异常体埋深较浅时，总水平导数法与归一化总水平导数法均能较好地识别出异常体的边界；当存在两种不同埋深的异常体时，深部异常体的异常电位信号特征不明显，受浅部异常体的影响，边界范围不能同时显示出来。而应用归一化总水平导数法能够在一定程度上增强深部异常体的响应特征。可见在一定的地质-地球物理条件下，归一化总水平导数法有望为准确圈定异常体的范围提供一种较为有效的手段，以提高井地电阻率法地质解释的正确率和精确度。

4) 归一化总水平导数法尚处于理论研究和数值模拟阶段,还有很多问题有待深入研究,如复杂的不均匀、不规则的异常体对异常电位的水平导数极值的影响;井地电阻率法中的多种装置的研究(三极、四极);异常电位水平导数极值随埋深变化的距离的定量计算等。

参 考 文 献

[1] 黄仲良.石油重磁电法勘探[M].北京:石油大学出版社,1999:1-386

Huang Z L.Grativity magnetism electricity method in petroleum prospecting[M].Beijing:Petroleum University Press,1999:1-386

[2] Scriba H.Computation of the electrical potential in the three-dimension structures[J].Geophysical Prospecting 1981,29(5):790-802

[3] Spitzer K.3-D finite-difference algorithm for DC resistivity modeling using conjugate methods[J].Geophysical Journal International,1995,123(3):903-914

[4] Pridmore D F,Hohmann G W,Ward S H,et al.An investigation of finite-element modeling for electrical and electromagnetic data in three dimensions[J].Geophysics,1981,46(7):1009-1024

[5] 徐凯军,李桐林.垂直有限线源三维地电场有限差分正演研究[J].吉林大学学报:地球科学版,2006,36(1):137-141

Xu K J,Li T L.The forward modeling of three dimensional geoelectric field of vertical finite line source by finite difference method[J].Journal of Jilin University(Earth Science Edition),2006,36(1):137-141

[6] 王志刚,何展翔,魏文博.井地直流电法三维数值模拟中若干问题研究[J].物探化探计算技术,2006,28(4):322-327

Wang Z G,He Z X,Wei W B.Study on some problems upon 3d modeling of dc borehole-ground method[J].Computing Techniques for Geophysical and Geochemical Exploration,2006,28(4):322-327

[7] 柯敢攀,黄清华.井地电法的三维正反演研究[J].北京大学学报:自然科学版,2009,45(2):264-272

Ke G P,Huang Q H.3D Forward and inversion problems of borehole-to-surface electrical method[J].Acta Scientiarum Naturalium Universitatis Pekinensis,2009,45(2):264-272

[8] 戴前伟,侯智超,王洪华.井地电位2.5维有限元数值模拟异常分析[J].物探化探计算技术,2013,35(4):457-463

Dai Q W,Hou Z C,Wang H H.Anomalous analysis of the 2.5D finite element numerical simulation of borehole-fround electrical method[J].Computing Techniques for Geophysical and Geochemical Exploration,2013,35(4):457-463

[9] 贾正元,张刚.三维井地电阻率法有限差分数值模拟研究[J].现代地质,2012,26(6):1199-2005

Jia Z Y,Zhang G.Research on finite-difference numerical simulation for 3D borehole-to-ground resistivity survey[J].Geoscience,2012,26(6):1199-2005

[10] 刘海飞,陈德鹏,戴前伟,等.连续电性介质线源井-地电位三维有限元数值模拟[J].桂林理工大学学报,2011,31(1):28-38

Liu H F,Chen D P,Dai Q W,et al.3D FEM modeling of borehole-surface potential with Line current source in semi-underground space of continuous variation of conductivity[J].Journal of Gu ilin University of Technology,2011,31(1):28-38

[11] 张天伦.用直流电阻率法确定油气藏边界的初步试验[J].石油地球物理勘探,1990,25(5):584-593

Zhang T L.An experiment in locating reservoir boundary using direct-current resistivity method[J].Oil Geophysical Prospecting,1990,25(5):584-593

[12] 张天伦.井中有套管情况下用直流电电阻率法确定油气藏边界的试验与研究[J].石油地球物理勘探,1993,28(3):314-324

Zhang T L.Experimental research on reservoir boundary determination using DC resistivity method at a cased borehole[J].Oil Geophysical Prospecting,1993,28(3):314-324

[13] 张天伦,张白林,聂荔.采用钻井套管作电极的非无穷远三极剖面法寻找剩余油[J].西南石油学院院报,1999,21(1):29-34.

Zhang T L,Zhang B L,Nie L.Looking for residual oils with nonboundless trielectrode section method using the drilling casing as electrodes[J].Journal of Southwest Petroleum Institute,1999,21(1):29-34

[14] 张天伦,张白林,聂荔.用地-井工作方式的三极梯度法寻找小块油气藏[J].石油地球物理勘探,1997,32(1):520-531

Zhang T L,Zhang B L,Nie L.Small reservoir discovery using ground borehole trielectrode gradient method[J].Oil Geophysical Prospecting,1997,32(1):520-531

[15] 张天伦,张白林,聂荔,等.非无穷远三极剖面法在新疆F油田的试验效果[J].西南石油学院学报,2001,23(1):5-10

Zhang T L,Zhang B L,Nie L,et al.The test result of nonboundless trielectrode profiling method in F field[J].Journal of Southwest Petroleum Institute,2001,23(1):5-10

[16] 汤井田,张继峰,冯兵,等.井地电阻率法歧离率确定高阻油气藏边界[J].地球物理学报,2007,50(3):926-931

Tang J T,Zhang J F,Feng B,et al .Determination of borders for resistive oil and gas reservoirs by deviation rate using the hole-to-surface resistivity method[J].Chinese Journal of Geophysics,2007,50(3):926-931

[17] 戴前伟,陈德鹏,陈勇雄.电法勘探中异常响应特征的增强算法及其实现[J].煤田地质与勘探,2013,41(3):75-78

Dai Q W,Chen D P,Chen Y X.The enhanced algorithms and its implementation for the abnormal response characteristics in electrical exploration[J].Coal Geology & Exploration,2013,41(3):75-78

[18] Ma G,Li L.Edge detection in potential fields with the normalized total horizontal derivative[J].Computers&Geosciences,2012,41(1):83-97

[19] 吴小平,汪彤彤.利用共轭梯度算法的电阻率三维有限元正演[J].地球物理学报,2003,46(3):428-432

Wu X P,Wang T T.A 3-D finite-element resistivity forward modeling using conjugate gradient algorithm[J].Chinese Journal of Geophysics,2003,46(3):428-432

[20] 柳建新,蒋鹏飞,童晓忠,等.不完全LU分解预处理的BICGSTAB算法在大地电磁二维正演模拟中的应用[J].中南大学学报:自然科学版,2009,40(2):484-491

Liu J X,Jiang P F,Tong X Z.Application of BICGSTAB algorithm with incomplete LU decomposition preconditioning to two dimensional magnetotelluric forward modeling[J].Journal of Central South University(Science and Technology),2009,40(2):484-491

[21] 王万银,邱之云,杨永,等.位场边缘识别方法研究进展[J].地球物理学进展,2010,25(1):196-210

Wang W Y,Qiu Z Y,Yang Y,et al.Some advance in the edge recognition of the potential field[J].Progress in Geophysics,2010,25(1):196-210

[22] 王万银.位场总水平导数极值位置位置空间变化规律研究[J].地球物理学报,2010,53(9):2258-2270

Wang W Y.Spatial variation law of the extreme value positions of total horizontal derivative for potential field data[J].Chinese Journal of Geophysics,2010,53(9):2258-2270

[23] Cordell L.Gravimetric expression of graben faulting in Santa Fe Country and the Espanola Basin[C]∥New Mexico Geological Society.30thFiled Conference.New Mexico:[s.n.],1979:59-64

[24] 徐世浙.地球物理中的有限单元法[M].北京:科学出版社,1994:120-122

Xu S Z.The finite element method in geophysics[M].Beijing:Science Press,1994:120-122