井-井电阻率2.5维近似成像研究

熊　飞 , 吕玉增, 李　勇

(1.华安奥特(北京)科技股份有限公司，北京　100085;2.桂林理工大学　地球科学学院，桂林　541004；3.广西矿冶与环境科学实验中心，桂林　541004)



井-井电阻率2.5维近似成像研究

熊飞1, 吕玉增2,3, 李勇2

(1.华安奥特(北京)科技股份有限公司，北京100085;2.桂林理工大学地球科学学院，桂林541004；3.广西矿冶与环境科学实验中心，桂林541004)

摘要：在使用有限单元法解决了井-井2.5维正演问题后，采用矩形剖分方式将井间剖面离散成若干个点，再将正演所得的视电阻率数据采用一种新颖的方式分配重组到井间每个离散点上，绘制出了井间剖面电阻率近似成像图。该近似成像方法避开了反演成像中庞大的内存需求和计算量，计算效率高，并且成像结果表明该方法具有较好的效果，是一种行之有效的井间视电阻率资料处理手段。

关键词：井-井电阻率法； 2.5维正演； 有限元； 近似成像

0引言

井-井电阻率法将供电点和测量点都置于钻孔中，更加接近目标体，因而在深部盲矿体的勘探方面相对地面电法更为有效。当下现有的浅部矿产资源日益枯竭，急需对深部隐伏矿床进行勘探，使得井-井电阻率成像研究具有一定的现实意义。

在上世纪90年代，国、内外就开始了井间电阻率成像的研究，Daily等[1]首先将有限单元法应用于井间二维电阻率成像中；Shima[2]又将有限单元法与α中心法联合用于井间电阻率图像的重建中；董清华[3]系统介绍了井间电阻率层析成像的一些进展；王兴泰等[4]使用佐迪反演进行了电阻率成像研究；底青云等[5]采用积分法实现了电阻率三维成像；Zhdanov等[6]研究了井间三维成像；Wilkinson等[7]和Chambers等[8]实现了矿井电阻率成像；A.Arato等[9]实现了基于交错网格的二维井间电阻率层析成像，提高了井间成像效果。

反演成像过程往往要涉及到大量的计算，使得成像时间增加，制约了算法的实用性，在提高成像速度方面，人们做了大量研究，取得了诸多成果：董清华[10]、刘国强等[11]对井间电阻率反演成像中偏导数的计算方法进行了对比研究；Knoll[12]、吴小平等[13]采用直接计算Jacobian矩阵向量积的方法，避开了Jacobian矩阵的计算和存储；吕玉增[14]、沈平等[15]直接避开反演方法探讨了井间电阻率的近似成像。

这里借鉴前人的思路，从点电流场在地下介质的分布情况以及地下异常球体的电位解析式出发，提出了一种新颖的井-井电阻率近似成像方案。算例表明，该方法只需2 s左右便可完成一次近似成像，并且具有较好的成像效果。

1井-井2.5维有限元数值模拟

对于复杂的地下介质分布，有限单元法是一种强有力的数值模拟方法[16]，因而这里的2.5维正演问题采用有限单元法解决。考虑到电位在电源点附近具有奇异性，如果直接采用总电位法求解，包含电源项，势必会导致求解的电源附近节点处的电位值具有较大误差，为了避免这一问题，作者采用异常电位法求解。

1.1二极观测

井-井观测一般采用在两口不同的井中布设供电电极和测量电极的工作方式，主要包括二极装置、三极装置和四极装置等。这里近似成像研究基于二极装置。它将供电电极A置于供电井中，供电电极B置于地面无穷远处；将测量电极M置于测量井中，测量电极N置于无穷远处，通过移动供电井中的供电电极A实现全孔供电，并且在每一个供电深度处，都移动测量井中的测量电极M进行全孔测量(图1)。与地面二极观测不同，井-井二极观测中视电阻率计算公式定义如下：

(1)

其中：UM为M处的测量电位；KAM为装置系数。

图1　井-井二极观测示意图Fig．1　Sketch of borehole-borehole pole-pole survey

1.2点源异常电位的变分问题及有限元正演

与地面电法不同，井-井观测供电点置于地下半空间中，因而在正常电位的表达式以及无穷远边界条件的推导中，都应设置关于地表对称的电源镜像点。经过推导，得到井-井2.5维点源异常电位变分问题：

(2)

其中：▽为哈密顿算子；σ为地下介质电导率；σ′为介质异常电导率；k为傅里叶变换中的波数；U0为正常傅氏电位；U为待求的异常傅氏电位；Ω为计算区域；Γ∞为无穷远边界；r1、r2分别为点电源A与其镜像点A*至无穷远边界Γ∞的距离；n为Γ∞外法线方向；K0、K1分别为第二类零阶和一阶修正贝塞尔函数。

同样根据镜像法原理，正常傅氏电位U0的表达式可写为式(3)。

(3)

其中：I为点源电流强度；σ0为介质正常电导率；r为待求点与电源点之间的距离；r′为待求点与电源镜像点之间的距离。

采用以井间区域为中心逐步向外稀疏的剖分方式，将地下通过两井的剖面剖分成若干个矩形单元，在每个单元内进行双线性插值，将式(2)中F(U)分解为各个矩形单元上的积分，再将各单元积分相加，可得式(4)。

(4)

根据F(U)取极值的条件有：

KU=-K′U0

(5)

其中：U、U0分别为所有节点处U、U0组成的列向量；系数矩阵K、K′为大型稀疏对称矩阵。

为了节省内存，采用改进的行压缩存储方法[17](Modified Sparse Row ，MSR)来存储系数矩阵里的非零元素，为了提高计算精度和效率，采用对称逐步超松弛预条件共轭梯度法[18](SSORPCG)法求解线性方程组。最后将求得的异常傅氏电位通过傅里叶反变换(表1为反变换中采用的一组波数k和系数g[16])，可得到各节点处的异常电位，与正常电位相加便得到了各节点处的总电位，再将电位换算成视电阻率。

表1　傅里叶反变换的k和g

2井-井电阻率近似成像

基于上述正演算法，可以模拟得到井间不同模型二极观测下的视电阻率数据，再以这些数据为基础，依照如下方式，实现井-井2.5维电阻率近似成像。

2.1近似成像原理

根据地下点电源电场的分布规律以及视电阻率的定义，在井-井二极观测中，测量点M处所得到的视电阻率值是供电点A和测量点M之间介质的综合反映，并且主要取决于AM连线附近一小块邻域内介质电阻率。

从异常球体的电位公式[19]可以发现：COS∠AOM直接影响着解析解(A为供电点、O为球心、M为测量点)。由于井间剖面可视为无数个离散的点P，每个点又可以近似当做一个“微小球体”，那么可以用COS∠APM的大小来限定AM连线邻域，并近似作为该A-M观测下视电阻率反映区域。如图2所示，对于任意一条连线AM，存在无数个近似椭圆的区域，记区域边界上任一点(A、M除外)为P*、区域内的任意一点为P，使得∠AP*M为一定值，且COS∠APM

(6)

其中：S=AM/d。

图2　AM邻域限定示意图Fig．2　Sketch of AM neighbourhood limited

由式(6)可以看出，参数S很好地确定了邻域范围；并且当选定一个固定S值后，AM增大时，d也会随之增大，这与电流实际分布情况较为符合。所以作者提出的近似思路在理论上具有一定的合理性。

基于以上思路，给出具体的井间视电阻率近似成像操作方案：

1)采用矩形剖分的方式将井间剖面近似离散成若干个点，横向剖分数为NX，纵向剖分数为NZ(图3)。

图3　近似成像网格化示意图Fig．3　Sketch of meshing for approximate image

2)对于一个格点P，依次试探每一组A-M，求出相应的COS∠APM，并与给定值MAXCOS进行比较，若COS∠APM≤MAXCOS，那么记P点处电阻率为该A-M观测下的视电阻率值，否则记为“0”，将所有的非零值的平均视为P点处电阻率值；如此得到所有点处的电阻率值，便可绘制近似成像图。其中，式(6)为MAXCOS的计算式，COS∠APM的计算式如(7)所示。

(7)

2.2井间模型的近似成像

通过对单一简单模型进行试算表明，当S取9时，会得到相对较好的成像效果，NX、NZ对成像影响较小。采用上述近似成像方法，取S=9、NX=10、NZ=50，对四种有代表性的井间二维模型、模型组合进行近似成像。其中，井深为H=100 m，井间距为D=20 m，背景电阻率为ρ0=100Ω·m。供电极A在供电井中移动范围为0 m～100 m，点距为2 m，共进行51次供电，供电强度为1A；每次供电时，测量极M在测量井中从0 m至100 m每隔2 m依次进行51次测量。

模型一为厚度为4m的倾斜高阻板状体，ρ1=1 000 Ω·m，具体如图4(a)所示。

模型二为厚度为4 m的“L”形高阻板状体，ρ1=1 000 Ω·m，具体如图4(b)所示。

图4　井间模型示意图Fig．4　Sketch of models between boreholes

模型三为井间上方赋存一4 m×4 m的低阻异常体，ρ1=1 0 Ω·m，下方赋存12 m×4 m的高阻异常体，ρ2=1 000 Ω·m。具体如图4(c)所示。

模型四为在同一深度赋存两个尺寸都为4 m×8 m的异常体，左边的为低阻异常：ρ1=1 0 Ω·m右边为高阻异常，ρ2=1 000 Ω·m，具体如图4(d)所示。

模型一至模型四的近似成像图如图5所示，图5中正方形、长方形等对应异常体截面实际形态和位置。

从图5中可以看出，①当异常体为倾斜板状时，异常区域显示为一倾斜椭圆，椭圆位置、范围以及倾斜方向和角度都与异常体实际对应较好；②当异常体为“L”形时，异常区域近似为“芒果“形态，大致上反映了异常体的位置以及在横向和纵向上的展布；③当两个异常体一上一下，形态不同、电阻率也不同时，近似成像图无论在位置、形态还是电阻率上都能加以反映；④当两个形态相同、电阻率不同的异常体处于同一深度，并且一左一右相互干扰时，近似成像图仍能较明显地辨别。

在AMD Phenom(tm) II N660 Dual-Core ProcessorCPU，3 GHz，2 GB内存，Windows 7平台下，每次成像时间只需2 s左右。除此之外，还应说明以下两点：

1)由于近似成像图是根据视电阻率数据重组而来，因而继承着视电阻率数据本身的特点——高阻区域对应低阻异常、低阻区域对应高阻异常，这点与反演成像有所区别。

2)近似成像图上下边界处存在明显的形如“山峰”的假异常，“峰顶”指向异常体。这时由于成像过程中采用的近似方式引起的，在图像解释中应注意区别。

图5　井间模型近似成像图Fig．5　The approximate image of models between borehole

3结语

作者采用有限单元法实现了井-井电阻率2.5维正演模拟，其中采用MSR法存储超大型稀疏系数矩阵，运用SSORPCG法求解线性方程组，保证了正演计算的效率和精度。根据地下点源电场的分布规律，提出了一种新颖的井-井电阻率2.5维近似成像方法，据此进行了四种典型模型的近似成像。结果显示，这里的井间电阻率近似成像方法能定性判断异常体电阻率、较准确反映异常体的位置、形态，并且具有极快的成像速度，因此可用于井-井观测数据的快速预处理，得到井间剖面较为直观的构造。

参考文献：

[1]DAILY W, OWEN E. Cross-borehole resistivity tomography[J]. Geophysics, 1991, 56( 8): 1228-1235.

[2]HIROMASA SHIMA. 2D and 3D resistivity image reconstruction using crosshole data[J]. Geophysics, 1992, 57( 10): 1270-1281.

[3]董清华. 井间电阻率层析成像的某些进展[J]. 地球物理学进展, 1997, 12(3): 77-89.

DONG Q H. The latest developments of crosshole resistivity tomography[J]. Progress in Geophysics, 1997, 12(3): 77-89. (In Chinese)

[4]王兴泰, 李晓芹. 电阻率图像重建的佐迪反演及其应用效果[J]. 物探与化探, 1996, 20(3): 228-233.

WANG X T, LI X Q. The image reconstruction of electric resistivity by Zohdy inversion and it’s application effects[J]. Geophysical&Geochemical Exploration, 1996, 20(3): 228-233. (In Chinese)

[5]底青云, 王妙月. 积分法三维电阻率成像[J]. 地球物理学报, 2001, 44(6): 844-851.

DI Q Y, WANG M Y. 3-D resistivity tomography by integral method[J]. Chinese J. Geophys，2001, 44(6): 844-851. (In Chinese)

[6] ZHDANOW M S, YOSHIOKA K. Cross-well electromagnetic imaging in three dimesions[J]. Exploration Geophysics, 2003,34:34-40.

[7]WILKINSON P B, CHANBERS J E,LELLIOTT M，et al. Extreme sensitivity of crosshole electrical resistivity tomography measurements to geometric errers[J]. Geophys, 2008,173:49-62.

[8]CHAMBERS J E, WILKINSON P B,WELLER AL，et al. Mineshaft imaging using surface and crosshole 3D electrical resistivity tomography :A case history from the East Pennine Coalfield[J]. UK. Journal of Applied-Geophysics, 2007,62(4):324-337.

[9]A. ARATO, A. GODIO,L.SAMBUELLI. Staggered grid inversion of cross hole 2-D resistivity tomography[J]. Journal of Applied-Geophysics, 2014,107:60-70.

[10]董清华, 田宪谟. 井间电阻率成像中Frchet导数的算法比较[J]. 物探化探计算技术, 1999, 19(1): 41-45.

DONG Q H, TIAN X M. A comparison of methods for calculating Frchet derivatives in crosshole resistivity image[J]. Computing Techniques for Geophysical and Geochemical Exploration, 1999, 19(1): 41-45. (In Chinese)

[11]刘国强, 陶果, 柯式镇, 等. 摄动求解Frchet导数的电阻率层析成像方法[J]. 地质地震, 2001, 23(2): 314-320.

LIU G Q, TAO G, KE S Z, et al. Improved reconstruction algorithms for resistivity tomography[J]. Seismology and Geology, 2001, 23(2): 314-320. (In Chinese)

[12]KNOLL D A, KEYES D E. Jacobian-free Newton-Krylov methods: a survey of approaches and applications[J]. Journal of Computational Physics, 2004, 193(2): 357-397.

[13]吴小平, 徐果明. 电阻率三维反演中偏导数矩阵的求取和分析[J]. 石油地球物理勘探, 1999, 34(4): 363-372.

WU X P, XU G M. Derivation and analysis of partial derivative matrix in resistivity 3-D inversion[J]. Oil Geophysical Prospecting, 1999, 34(4): 363-372. (In Chinese)

[14]吕玉增, 阮百尧, 黄俊革. 直流电井间三维直接成像[J]. 物探化探计算技术, 2003, 25(1): 60-64.

LV Y Z, RUAN B Y, HUANG J G. The 3-D immediate crosshole tomography with direct current[J]. Computing Techniques for Geophysical and Geochemical Exploration, 2003, 25(1): 60-64. (In Chinese)

[15]沈平, 强建科, 李永军,等. 井间视电阻率的几何成像方法[J]. 中南大学学报：自然科学版, 2010, 41(3): 1079-1084.

SHEN P. QIANG J K, LI Y J, et al. Geometry image method of crosshole apparent resistivity[J]. Journal of Central South University (Science and Technology), 2010, 41(3): 1079-1084. (In Chinese)

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

XU S Z. The finite element method in Geophysics[M]. Beijing: Science Press, 1994.(In Chinese)

[17]SAAD Y. Iterative methods for sparse linear systems[M]. Siam, 2003.

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

WU X P, WANG T T. A 3-D finite-element resistivity forward modeling using conjugate gradient algorithm[J]. Chinese J. Geophys, 2004, 46(3): 428-432. (In Chinese)

[19]李金铭. 地电场与电法勘探[M]. 北京: 地质出版社, 2005.

LI J M. Geoelectric field and electrical exploration[M]. Beijing: Geological Publishing House, 2005. (In Chinese)

收稿日期：2015-04-10改回日期：2015-05-21

基金项目:国家“十二五”国家重大科学仪器设备开发专项(2012YQ030126)；广西自然科学基金(2013GXNSFBA019212)；广西高校科学技术研究重点项目(2013ZD029)

作者简介：熊飞(1990-)男，硕士，主要从事电磁法数值模拟研究，E-mail：xiongfei357@163.com。

文章编号：1001-1749(2016)03-0308-06

中图分类号：P 631.3

文献标志码：A

DOI：10.3969/j.issn.1001-1749.2016.03.03

A study on approximate imaging of borehole-borehole 2.5D resistivity method

XIONG Fei1, LV Yu-zeng2，3, LI Yong2

(1.HUAAN AOTE(Beijing) Technology Co.,Ltd,Beijing100085，China；2.Guilin University of Technology,College of Earth Science,Guilin541004,China;3.Guangxi Scientific Experiment Center of Mining,Metallurgy and Environment,Guilin541004,China)

Abstract:After solving the borehole-borehole 2.5 dimensional forward problem by using the finite element method, we can use the rectangular partition discretized the section between boreholes into a number of points. According to the distribution of the electric field of a point source in the underground space, the apparent resistivity data can be re-allocated on each discrete points between boreholes through a novel way. Then we can draw the approximate resistivity image across hole profile. The approximate imaging method has high computational efficiency as avoiding the large memory requirements and calculations. And the imaging results show that the method is a kind of effective processing method of apparent resistivity data between boreholes because of it’s good effect.

Key words:borehole-borehole resistivity method; 2.5D forward; finite element method; approximate imaging