地-井与井-地及井-井时域激电异常响应特征分析

陈汉波， 熊　彬

(桂林理工大学　地球科学学院，桂林　541004)



地-井与井-地及井-井时域激电异常响应特征分析

陈汉波， 熊彬*

(桂林理工大学地球科学学院，桂林541004)

摘要：井中激电是寻找深部金属矿产的有效方法之一，这里基于有限差分和异常电位算法实现了时间域井中激发极化法的三维数值模拟。首先直接给出点源三维异常电位场满足的基本微分方程；然后采用一维非零元素行压缩存储模式对所形成的大型稀疏矩阵进行存储，节约了计算所需内存空间；同时引入不完全Cholesky分解稳定化双共轭梯度(ICBG)法求解有限差分线性方程组，提高了求解的效率。构建典型的地电模型，采用地-井、井-地、井-井的观测模式分别进行正演计算，结果表明，时间域井中激电视极化率异常响应较视电阻率异常响应灵敏，综合考虑这两种视参量的空间分布特征，在一定程度上可以提高对相关异常体的分辨率，为今后野外的实际数据资料的解释提供依据。

关键词：有限差分法； 激发极化； IBCG法； 一维非零元素压缩

0引言

激发极化法是根据岩、矿石的激发极化效应解决地质问题的一类地球物理勘探方法。而井中激发极化法(Induced polarization well logging)则是在钻孔中运用激发极化法，并与其他井中物探方法结合运用，可以有效地解决在勘探金属矿产中所遇到的一系列问题。

从七十年代开始，我国就比较系统地对井中激电进行理论研究和实验应用，经过半个世纪的发展，取得了丰硕的成果。吴小平等[1-4]采用有限差分法对直流点电源作用下的三维地电体进行数值模拟；潘和平等[5]结合解析法和数值法实现了地-井时域激电拟牛顿反演研究，提高了反演的精度；吴冰[6]进行了地-井激电法探测盲矿体的数字模拟研究；周峰等[7-8]以球体为研究对象，采用有限差分法实现了井中激电地-井方式井旁球体正反演研究；吕玉增[9-10]实现了地-井、井-地IP三维快速正反演研究，对异常特征进行了分析，同时编制了相应的解释软件；李长伟[11]进行井中激发极化法正反演，同时探讨了相应的快速迭代求解技术；张艳国[12]等进行了地-井断面IP正演模拟研究，增大了地-井IP反演的数据量，提高了反演的精度；郭刚等[13]对井中激电地-井方式下的极化率异常特征进行了计算研究，同时结合实际应用加以说明。以上研究成果对于井中激电的应用具有很大的指导意义。然而，这些研究大多仅考虑井-地或者地-井情形，没有考虑井-井观测模式，也没能同时处理地-井、井-地、井-井情形下的井中激电正演问题，同时计算的效率也有待进一步的提高。作者在前人的基础上，在井-地、地-井、井-井的观测模式下对时间域井中激发极化法三维有限差分进行正演研究；探讨了时间域井中激发极化法正演求解的原理和方法；给出了相应的迭代求解算法过程，对井-地、地-井、井-井IP异常空间分布特征进行分析，以期对井中激电的异常规律作进一步揭示，进而指导野外的实际资料的地质解释。

1边值问题

稳定电流场所满足的微分方程：

·[σ(x,y,z)u(x,y,z)]=

-Iδ(x-x0)δ(y-y0)δ(z-z0)∈Ω

(1)

相应的边界条件

(2)

(3)

式中：σ为电导率；I为电流强度；u为电位；(z,y,z)、(x0,y0,z0) 分别为观测点和点电流源的位置；δ为狄拉克函数；Γs为地面边界；Γ∞为地下无穷边界；n为边界的外法向矢量；r为源点到边界上节点的矢径；θ为外法向矢量n和矢径r之间的夹角。为了消除源点的奇异性，假设电流源附近电导率为δm的均匀半空间或者均匀全空间电位为u0，由无源区介质电阻率分布引起的电位为us，于是

u=us+u0

(4)

显然有源区电位满足式(5)。

σm2u0=-Iδ(x-x0)δ(y-y0)δ(z-z0)(5)

结合式(1)和式(5)可得

·[σus)=-[(σ-σm)u0]

(6)

其相应的边界条件

(7)

(8)

对整个计算区域单元进行差分离散，具体推导过程见文献[14]。最后可得线性方程组，见式(9)。

Kus=b

(9)

式中:K为总体系数矩阵；b为右端项。求解此方程组，便可得us，利用式(4)即可求得总电位值u,根据文献[15]即可计算出极化率。

2线性方程组的求解

经过有限差分法离散形成的系数矩阵是对称正定的，共轭梯度(CG)法是求解此类问题优先选择的方法之一。共轭梯度法在求解条件数很小的线性方程组时，收敛速率快，若线性方程组的条件数很大，则不易收敛。这里所形成的线性方程组条件数很大，为了提高其收敛速度，要对其进行预处理，改善其条件数。作者采用不完全分解求解预条件矩阵，之后结合稳定化双共轭梯度(BICGSTAB)法求解差分方程组。根据Meijerink 等[16]提出的不完全Cholesky分解公式

A=MMT

(10)

其中

M=UD(-1/2)

(11)

及

(12)

式中: U为下三角矩阵；Uii=Dii。将方程(9)改写成方程(13)。

[M-1A(MT)-1](MTu)=M-1b

(13)

将M、MT分别看做稳定化双共轭梯度算法中左右预条件矩阵，则可形成不完全Cholesky分解稳定化双共轭梯度算法，其具体迭代过程[17]如式(14)所示。

P(i)=r{i-1)+β(i-1)(P(i-1)-ωi-1ν(i-1))

s=r(i-1)-αiν(i)

(14)

ωi=(tTs)/(tTt)

r(i)=s-ωit

βi=(ρi/ρi-1/(ωi/αi)

enddo

3数值算例

3.1算例1

为了验证算法的正确性，设计三层水平大地模型。第一层电阻率为50 Ω·m，厚度为5 m，第二层电阻率为100 Ω·m，厚度为 10 m，第三层电阻率为20 Ω·m。供电电极在地面，采用对称四极进行测量，解析解采用数字滤波算法进行计算，数值解与解析解对比如图1所示。由图1可知,数值解与解析解基本一致，这里采取的算法解点源三维地电场是可靠有效的。针对模型1，采用不同的算法，不同的网格大小在CPU Intel(R) Core(TM) i3 M380内存2G的计算机上进行了试算。表1为ICCG和ICBG这两种计算方法在相同的网格大小下的计算时间和迭代次数的对比表。由表1可知，无论在计算时间上，还是在迭代次数方面，不完全稳定化双共轭梯度算法(ICBG)，优于不完全Cholesky共轭梯度法(ICCG)。

3.2算例2

如图2所示，井深500 m,井口坐标(0,0,0)，井旁有一板状异常体，长、宽、高分别为20 m、100 m、50 m，其中心位于(110，0,225)处。采用地-井测量模式进行观测，将供电点分别置于A1=(0,0,0)、A2=(-20,0,0)、A3=(20，0，0)处供电，采用二极装置沿井进行测量。计算结果分别如图3和图4所示。

图1　数值解与解析解Fig．1　Comparison between analytical and 　　numerical solutions

图2　模型过XOZ断面图Fig．2　The XOZ section figure of model

方法50×50×50100×100×100时间/s迭代次数时间/s迭代次数不完全Cholesky共轭梯度法(ICCG)711260183不完全Cholesky稳定化共轭梯度法(ICBG)5624095

图3　不同供电点视电阻率曲线图Fig．3　The apparent resistivity curves

图4　不同供电点视极化率曲线图Fig．4　The apparent chargeability curves

从图3、图4可以看出，视电阻率与视极化率均可以反映出低阻体的存在，视电阻率在异常体上边界达到极小值，视极化率在异常体上边界达到极大值。

3.3算例3

如图5所示，有一竖井，井口位于处。井旁有两个相同大小的异常体：①低阻高极化体，长、宽、高分别为30 m、30 m、20 m，异常体中心为(-25,-25,60)；②高阻高极化体，异常体中心为(25,25,40)。采用二极装置进行井-地测量，供电电极分别定在地下(0,0,20)、(0,0,0)、(0，0，70)处供电，在地表上沿着测线逐点测量。

地表视电阻率和视极化率平面等值线如图6所示。

图5　模型平面图Fig．5　The section figure of model

图6　不同供电点的地表视电阻率/极化率等值线图Fig．6　Ground apparent resistivity and chargeability anomalycontour section of different point sources model(a)供电点(0,0,20)地表视电阻率等值线图;(b)供电点(0,0,20)地表视极化率等值线图;(c)供电点(0,0,25)地表视电阻率等值线图;(d)供电点(0,0,25)地表视极化率等值线图;(f)供电点(0,0,70)地表视电阻率等值线图;(g)供电点(0,0,70)地表视极化率等值线图

从图6(a)、图6(c)可知，地表视电阻率等值线，表现为一个高阻异常和一个低阻异常形态，与异常体在地面的投影位置基本吻合。由于高阻异常体比低阻异常体更靠近地表和供电点，故高阻异常体所形成的视电阻率异常较高于低阻异常体所引起的异常响应。随着供电点的深度加大，高阻体的地表视电阻率响应依然明显，然而，低阻体异常响应却变弱。从图6(f)可以看出，随着供电点位于异常体下面，低阻体在地表投影位置视电阻率呈高阻响应，而高阻体在地表投影位置对应的视电阻率呈低阻响应，这是低阻体吸引电流，高阻排斥电流所致。图6(b)为地表视极化率等值线图，表现为两个高极化异常形态，显示出两个高极化体的存在。高阻体的异常响应要比低阻体更加明显，低阻体的异常响应中心往右上角方向偏移。由图6(d)可以看出，随着供电点的深度加大，异常中心由两个变为一个，对高阻体在地表投影位置比较吻合，主要反映的是高阻体。由图6(g)可看到，随着供电点处于两个异常体下面，地表的视极化率等值线在两个异常体在地表投影位置之间的中心呈高极化响应，此时的高极化响应是两个异常体的综合反映。

3.4算例4

如图7所示，有两竖井，井b1井口位于(0,0,0)处，作为测量井。井b2井口位于(-100,0,0)处，作为供电井。井旁有一异常体，长、宽、高均为50 m。其中心位于(0,0,250)处，井b1垂直穿过异常体中心。

采用井-井测量方式，设供电点分别位于A1=(-100,0,250)、A2=(-100,0,250)、A3=(-100,0,275)处供电，采用二极装置沿井进行测量，计算结果如图8和图9所示。由图8、图9可以看出，视电阻率曲线和视极化率曲线同样能够反映出低阻体的存在，在异常体的上下边界视电阻率和视极化率均达到极值，两者均很好地反映出异常体的大致位置。

图7　模型过XOZ断面图Fig．7　The XOZ section figure of model

图8　不同供电点的视电阻率率曲线Fig．8　The apparent resistivity curves

图9　不同供电点的视极化率曲线Fig．9　The apparent chargeability curves

4结 论

作者研究了时域井中激电的三维正演问题，运用异常电位法消除了源点奇异性，提高源附近解的精度。引入不完全分解稳定化双共轭梯度法迭代法求解线性方程组，快速实现了线性方程组的求解。值得注意的是，该方法同样可以求解复数域的大型线性方程组。

通过分析不同测量模式下的异常空间分布特征可知，采用井-地面积测量方式，可以大致了解异常体在地面的投影位置，但这只能定性地分析。若要精确定位，还需进行反演研究。地-井测量方式可知异常体的大小范围；井-井测量模式有助探测井旁盲矿体，了解异常体与已知矿体之间的连接关系。总体而言，时域井中激电视极化率异常响应较视电阻率响应灵敏，综合此两种视参量，有助于提高对异常体的分辨率。算例结果验证了方法的可靠性和高效性，能够很好反映出异常场的特征，从而为进一步进行井中激电反演和解释研究奠定基础。

参考文献：

[1]吴小平,王彤彤. 利用共轭梯度方法的电阻率三维反演若干问题研究[J].地震地质，2001，23(2)：321-327.

WU X P, WANG T T. Stoudy on some problems for 3D resistivity inversion using conjugate gradient[J]. Seimology and geology,2001,23(2):321-327. (In Chinese)

[2]邓正栋,关洪军,万乐,等. 稳定点电流源场三维有限分正演模拟[J].解放军理工大学学报,2000，1(3)：46-50.

DENG Z D ,GUANG H J,WAN L,et al.Three -Dimensional finite differential front demonstrates simulation of stable point electrical source field[J]. Journal of PLA Universi ty of Science and Technolog y,2000,1(3):46-50.(In Chinese)

[3]宋维琪,刘仕友. 三维地下埋藏点电流源有限差分正演计算[J].石油物探，2006，45(3)：324-328.

SONG W Q,LIU S Y.A new limited difference forward method on 3D geo-electrical medium with embedding point current source[J].Geophysical Prospectinc for Petroleum.2006,45(3):324-328. (In Chinese)

[4]胡朝俊. 三维直流电法共轭梯度反演算法研究[D].北京：中国地质大学，2006.

HU C J. Study of three-dimensional conjugate gradients inversion of DC resistivity data.[D].Master dissertation[D].Beijing：China University of Geosciences，2006. (In Chinese)

[5]潘和平,孟庆鑫. 地-井时域激电拟牛顿反演研究[J].电报科学学报，2013，28(5):987-993.

PAN H P,MENG Q X. The research on quasi-newton methods for inverse problem of surface-borehole induced polarization measurement data[J]. Chinese journal of radio science，2013,28(5):987-993.(In Chinese)

[6]吴 冰. 地-井激电法探测盲矿体的数字模拟研究[D].西安：长安大学，2008.

WU B . Study on digital simulation of surface-borehole IP method detection the blind ore body[D].Xi an: Chang’an University,2008.(In Chinese)

[7]周峰,潘和平,吴国军,等. 井中激电地-井方式井旁球体正反演研究[J].物探与化探,2008,32(3):321-325.

ZHOU F，PAN H P，WU G J,et al. Forward and inversion of sphere besidewell in the method of ground-well induced polarization[J]. Geophysical& geochemical exploration,2008,32(3):321-325.(In Chinese)

[8]周峰,潘和平,吴国军,等. 基于有限差分法的井中激发极化法正演[J].电报科学学报，2010,25(4):785-791.

ZHOU F，PAN H P，WU G J, et al. Three-dimensional forward modeling of induced of polarizationborehole using finite difference method[J]. Chinese journal of radio science,2010，25(4)：785-791. (In Chinese)

[9]吕玉增.地-井，井-地三维快速正反演研究[D].长沙:中南大学，2008.

LV Y Z. The research on 3-D fast forward and inversion of surface—borehole and borehole-surface IP methods[D].Changsha: Central South University,2008. (In Chinese)

[10]吕玉增,阮百尧, 彭苏萍.地-井方位激电观测异常特征研究[J].地球物理学进展，2012，27(1),202-216.

LV Y Z,RUAN B Y,PENG S P.A study on anomaly of surface-borehole direction induced polarization survey[J].Progress in Geophys,2011,26(6):201-216. (In Chinese)

[11]李长伟. 井中激发极化法正反演及快速迭代求解技术研究[D].长沙：中南大学，2011.

LI C W. Study on forward and inversion of induced-polarization well logging and fast iterative solvers[D].Changsha: Central South University,2011. (In Chinese)

[12]张艳国，吕玉增，刘 正，等. 地-井IP断面正演模拟研究[J].工程地球物理学报，2013，10(1)：46-50.

ZHANG Y G,LV Y Z,LIU Z, et al. Surface-borehole cross-section IP Method numersical modeling and data explanation[J]. Chinese journal of engineering geophysics,2013,10(1)：46-50. (In Chinese)

[13]郭 刚，鄂博军. 井中激电地-井方式极化率特征异常正演计算及应用[J].矿产勘查，2014，5(2)：307-311.

GUO G,E B J.Forward calculation and application of the polarization anomalies of ground- well induced polarization method[J]. Mineral exploration,2014,5(2):307-311. (In Chinese)

[14]周熙囊．电法勘探数值模拟技术[M]．四川：四川科学技术出版社，1986.(In Chinese)

ZHOU X Y. Electrical prospecting numericalsimulation technology[M].Sichuang: Sichuan Science and Technology Press ,1986. (In Chinese)

[15]SEIGEL H O. Mathematical formulation and type curves for induced polarization[J].Geophysics，1959, 24(3)：547-565.

[16]MEIJER INK J A, VANDERVORST H A. An iterative solution method for linear system of which the coefficient matrix is a symmetric M - Matrix [J]. Math comp，1977，31：148.

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

WANG Z G,HE Z X,WEI W B. Study on some problems upon 3D modeling of dc borehole -grund methods[J]. Computing technqiques geophysical and geochemical exploration,2006,28(4):323-327. (In Chinese)

收稿日期：2015-03-12改回日期：2015-04-02

基金项目:国家自然科学基金项目(201300541164004)；广西自然科学基金项目(2011GXNSFA018003，2013GXNSFAA019277)；桂林市“漓江学者”专项(2013005)

作者简介：陈汉波(1990-)，男，硕士，主要从事电磁场数值模拟与反演研究，E-mail：563179536@qq.com 。 *通信作者：熊彬(1974-)，男，博导，主要从事电磁法数据模拟与反演成像研究，E-mail：xiongbin@glute.cdu.cn。

文章编号：1001-1749(2016)03-0314-07

中图分类号：P 631.2

文献标志码：A

DOI：10.3969/j.issn.1001-1749.2016.03.04

The anomaly spatial distribution of surface-borehole,borehole-surface and borehole-borehole IP method

CHEN Han-bo， XIONG Bin

(College of Earth Sciences, Guilin University of Technology, Guilin541004, China)

Abstract:Induced polarization well logging is a preferred method for surveying deep metal mineral resources. The main task in this paper is to solve the forward problem of induced polarization well logging by using the finite difference method and the abnormal field Method. First, the 3D anomalous potential basic differential equation of point source is derivated in detailed. The technique of the one-dimensional non-zero element row compression is used to store the sparse coefficient matrix. The result shows that it can save much memory space and reduce compute capacity. It is high efficiency to solve the system of linear equations with the method. After the forward of the surface-borehole, borehole-surface and borehole-borehole IP method. The results display that the apparent chargeability anomaly is more sensitive than the apparent resistivity. Considering the distribution space characteristics of the apparent parametric, it can, to a certain extent, improve the resolution of abnormal body and provide a basis to analyze actual field data in future.

Key words:finite difference method; induced polarization; ICBG method; one-dimensional non-zero element row compression