熊　飛 , 呂玉增, 李　勇

(1.華安奧特(北京)科技股份有限公司，北京　100085;2.桂林理工大學(xué)　地球科學(xué)學(xué)院，桂林　541004；3.廣西礦冶與環(huán)境科學(xué)實(shí)驗(yàn)中心，桂林　541004)

?

井-井電阻率2.5維近似成像研究

熊飛1, 呂玉增2,3, 李勇2

(1.華安奧特(北京)科技股份有限公司，北京100085;2.桂林理工大學(xué)地球科學(xué)學(xué)院，桂林541004；3.廣西礦冶與環(huán)境科學(xué)實(shí)驗(yàn)中心，桂林541004)

摘要：在使用有限單元法解決了井-井2.5維正演問題后，采用矩形剖分方式將井間剖面離散成若干個(gè)點(diǎn)，再將正演所得的視電阻率數(shù)據(jù)采用一種新穎的方式分配重組到井間每個(gè)離散點(diǎn)上，繪制出了井間剖面電阻率近似成像圖。該近似成像方法避開了反演成像中龐大的內(nèi)存需求和計(jì)算量，計(jì)算效率高，并且成像結(jié)果表明該方法具有較好的效果，是一種行之有效的井間視電阻率資料處理手段。

關(guān)鍵詞：井-井電阻率法； 2.5維正演； 有限元； 近似成像

0引言

井-井電阻率法將供電點(diǎn)和測量點(diǎn)都置于鉆孔中，更加接近目標(biāo)體，因而在深部盲礦體的勘探方面相對地面電法更為有效。當(dāng)下現(xiàn)有的淺部礦產(chǎn)資源日益枯竭，急需對深部隱伏礦床進(jìn)行勘探，使得井-井電阻率成像研究具有一定的現(xiàn)實(shí)意義。

在上世紀(jì)90年代，國、內(nèi)外就開始了井間電阻率成像的研究，Daily等[1]首先將有限單元法應(yīng)用于井間二維電阻率成像中；Shima[2]又將有限單元法與α中心法聯(lián)合用于井間電阻率圖像的重建中；董清華[3]系統(tǒng)介紹了井間電阻率層析成像的一些進(jìn)展；王興泰等[4]使用佐迪反演進(jìn)行了電阻率成像研究；底青云等[5]采用積分法實(shí)現(xiàn)了電阻率三維成像；Zhdanov等[6]研究了井間三維成像；Wilkinson等[7]和Chambers等[8]實(shí)現(xiàn)了礦井電阻率成像；A.Arato等[9]實(shí)現(xiàn)了基于交錯(cuò)網(wǎng)格的二維井間電阻率層析成像，提高了井間成像效果。

反演成像過程往往要涉及到大量的計(jì)算，使得成像時(shí)間增加，制約了算法的實(shí)用性，在提高成像速度方面，人們做了大量研究，取得了諸多成果：董清華[10]、劉國強(qiáng)等[11]對井間電阻率反演成像中偏導(dǎo)數(shù)的計(jì)算方法進(jìn)行了對比研究；Knoll[12]、吳小平等[13]采用直接計(jì)算Jacobian矩陣向量積的方法，避開了Jacobian矩陣的計(jì)算和存儲(chǔ)；呂玉增[14]、沈平等[15]直接避開反演方法探討了井間電阻率的近似成像。

這里借鑒前人的思路，從點(diǎn)電流場在地下介質(zhì)的分布情況以及地下異常球體的電位解析式出發(fā)，提出了一種新穎的井-井電阻率近似成像方案。算例表明，該方法只需2 s左右便可完成一次近似成像，并且具有較好的成像效果。

1井-井2.5維有限元數(shù)值模擬

對于復(fù)雜的地下介質(zhì)分布，有限單元法是一種強(qiáng)有力的數(shù)值模擬方法[16]，因而這里的2.5維正演問題采用有限單元法解決。考慮到電位在電源點(diǎn)附近具有奇異性，如果直接采用總電位法求解，包含電源項(xiàng)，勢必會(huì)導(dǎo)致求解的電源附近節(jié)點(diǎn)處的電位值具有較大誤差，為了避免這一問題，作者采用異常電位法求解。

1.1二極觀測

井-井觀測一般采用在兩口不同的井中布設(shè)供電電極和測量電極的工作方式，主要包括二極裝置、三極裝置和四極裝置等。這里近似成像研究基于二極裝置。它將供電電極A置于供電井中，供電電極B置于地面無窮遠(yuǎn)處；將測量電極M置于測量井中，測量電極N置于無窮遠(yuǎn)處，通過移動(dòng)供電井中的供電電極A實(shí)現(xiàn)全孔供電，并且在每一個(gè)供電深度處，都移動(dòng)測量井中的測量電極M進(jìn)行全孔測量(圖1)。與地面二極觀測不同，井-井二極觀測中視電阻率計(jì)算公式定義如下：

(1)

其中：UM為M處的測量電位；KAM為裝置系數(shù)。

圖1　井-井二極觀測示意圖Fig．1　Sketch of borehole-borehole pole-pole survey

1.2點(diǎn)源異常電位的變分問題及有限元正演

與地面電法不同，井-井觀測供電點(diǎn)置于地下半空間中，因而在正常電位的表達(dá)式以及無窮遠(yuǎn)邊界條件的推導(dǎo)中，都應(yīng)設(shè)置關(guān)于地表對稱的電源鏡像點(diǎn)。經(jīng)過推導(dǎo)，得到井-井2.5維點(diǎn)源異常電位變分問題：

(2)

其中：▽為哈密頓算子；σ為地下介質(zhì)電導(dǎo)率；σ′為介質(zhì)異常電導(dǎo)率；k為傅里葉變換中的波數(shù)；U0為正常傅氏電位；U為待求的異常傅氏電位；Ω為計(jì)算區(qū)域；Γ∞為無窮遠(yuǎn)邊界；r1、r2分別為點(diǎn)電源A與其鏡像點(diǎn)A*至無窮遠(yuǎn)邊界?！薜木嚯x；n為?！尥夥ň€方向；K0、K1分別為第二類零階和一階修正貝塞爾函數(shù)。

同樣根據(jù)鏡像法原理，正常傅氏電位U0的表達(dá)式可寫為式(3)。

(3)

其中：I為點(diǎn)源電流強(qiáng)度；σ0為介質(zhì)正常電導(dǎo)率；r為待求點(diǎn)與電源點(diǎn)之間的距離；r′為待求點(diǎn)與電源鏡像點(diǎn)之間的距離。

采用以井間區(qū)域?yàn)橹行闹鸩较蛲庀∈璧钠史址绞剑瑢⒌叵峦ㄟ^兩井的剖面剖分成若干個(gè)矩形單元，在每個(gè)單元內(nèi)進(jìn)行雙線性插值，將式(2)中F(U)分解為各個(gè)矩形單元上的積分，再將各單元積分相加，可得式(4)。

(4)

根據(jù)F(U)取極值的條件有：

KU=-K′U0

(5)

其中：U、U0分別為所有節(jié)點(diǎn)處U、U0組成的列向量；系數(shù)矩陣K、K′為大型稀疏對稱矩陣。

為了節(jié)省內(nèi)存，采用改進(jìn)的行壓縮存儲(chǔ)方法[17](Modified Sparse Row ，MSR)來存儲(chǔ)系數(shù)矩陣?yán)锏姆橇阍兀瑸榱颂岣哂?jì)算精度和效率，采用對稱逐步超松弛預(yù)條件共軛梯度法[18](SSORPCG)法求解線性方程組。最后將求得的異常傅氏電位通過傅里葉反變換(表1為反變換中采用的一組波數(shù)k和系數(shù)g[16])，可得到各節(jié)點(diǎn)處的異常電位，與正常電位相加便得到了各節(jié)點(diǎn)處的總電位，再將電位換算成視電阻率。

表1　傅里葉反變換的k和g

2井-井電阻率近似成像

基于上述正演算法，可以模擬得到井間不同模型二極觀測下的視電阻率數(shù)據(jù)，再以這些數(shù)據(jù)為基礎(chǔ)，依照如下方式，實(shí)現(xiàn)井-井2.5維電阻率近似成像。

2.1近似成像原理

根據(jù)地下點(diǎn)電源電場的分布規(guī)律以及視電阻率的定義，在井-井二極觀測中，測量點(diǎn)M處所得到的視電阻率值是供電點(diǎn)A和測量點(diǎn)M之間介質(zhì)的綜合反映，并且主要取決于AM連線附近一小塊鄰域內(nèi)介質(zhì)電阻率。

從異常球體的電位公式[19]可以發(fā)現(xiàn)：COS∠AOM直接影響著解析解(A為供電點(diǎn)、O為球心、M為測量點(diǎn))。由于井間剖面可視為無數(shù)個(gè)離散的點(diǎn)P，每個(gè)點(diǎn)又可以近似當(dāng)做一個(gè)“微小球體”，那么可以用COS∠APM的大小來限定AM連線鄰域，并近似作為該A-M觀測下視電阻率反映區(qū)域。如圖2所示，對于任意一條連線AM，存在無數(shù)個(gè)近似橢圓的區(qū)域，記區(qū)域邊界上任一點(diǎn)(A、M除外)為P*、區(qū)域內(nèi)的任意一點(diǎn)為P，使得∠AP*M為一定值，且COS∠APM

(6)

其中：S=AM/d。

圖2　AM鄰域限定示意圖Fig．2　Sketch of AM neighbourhood limited

由式(6)可以看出，參數(shù)S很好地確定了鄰域范圍；并且當(dāng)選定一個(gè)固定S值后，AM增大時(shí)，d也會(huì)隨之增大，這與電流實(shí)際分布情況較為符合。所以作者提出的近似思路在理論上具有一定的合理性。

基于以上思路，給出具體的井間視電阻率近似成像操作方案：

1)采用矩形剖分的方式將井間剖面近似離散成若干個(gè)點(diǎn)，橫向剖分?jǐn)?shù)為NX，縱向剖分?jǐn)?shù)為NZ(圖3)。

圖3　近似成像網(wǎng)格化示意圖Fig．3　Sketch of meshing for approximate image

2)對于一個(gè)格點(diǎn)P，依次試探每一組A-M，求出相應(yīng)的COS∠APM，并與給定值MAXCOS進(jìn)行比較，若COS∠APM≤MAXCOS，那么記P點(diǎn)處電阻率為該A-M觀測下的視電阻率值，否則記為“0”，將所有的非零值的平均視為P點(diǎn)處電阻率值；如此得到所有點(diǎn)處的電阻率值，便可繪制近似成像圖。其中，式(6)為MAXCOS的計(jì)算式，COS∠APM的計(jì)算式如(7)所示。

(7)

2.2井間模型的近似成像

通過對單一簡單模型進(jìn)行試算表明，當(dāng)S取9時(shí)，會(huì)得到相對較好的成像效果，NX、NZ對成像影響較小。采用上述近似成像方法，取S=9、NX=10、NZ=50，對四種有代表性的井間二維模型、模型組合進(jìn)行近似成像。其中，井深為H=100 m，井間距為D=20 m，背景電阻率為ρ0=100Ω·m。供電極A在供電井中移動(dòng)范圍為0 m～100 m，點(diǎn)距為2 m，共進(jìn)行51次供電，供電強(qiáng)度為1A；每次供電時(shí)，測量極M在測量井中從0 m至100 m每隔2 m依次進(jìn)行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)所示。

模型四為在同一深度賦存兩個(gè)尺寸都為4 m×8 m的異常體，左邊的為低阻異常：ρ1=1 0 Ω·m右邊為高阻異常，ρ2=1 000 Ω·m，具體如圖4(d)所示。

模型一至模型四的近似成像圖如圖5所示，圖5中正方形、長方形等對應(yīng)異常體截面實(shí)際形態(tài)和位置。

從圖5中可以看出，①當(dāng)異常體為傾斜板狀時(shí)，異常區(qū)域顯示為一傾斜橢圓，橢圓位置、范圍以及傾斜方向和角度都與異常體實(shí)際對應(yīng)較好；②當(dāng)異常體為“L”形時(shí)，異常區(qū)域近似為“芒果“形態(tài)，大致上反映了異常體的位置以及在橫向和縱向上的展布；③當(dāng)兩個(gè)異常體一上一下，形態(tài)不同、電阻率也不同時(shí)，近似成像圖無論在位置、形態(tài)還是電阻率上都能加以反映；④當(dāng)兩個(gè)形態(tài)相同、電阻率不同的異常體處于同一深度，并且一左一右相互干擾時(shí)，近似成像圖仍能較明顯地辨別。

在AMD Phenom(tm) II N660 Dual-Core ProcessorCPU，3 GHz，2 GB內(nèi)存，Windows 7平臺(tái)下，每次成像時(shí)間只需2 s左右。除此之外，還應(yīng)說明以下兩點(diǎn)：

1)由于近似成像圖是根據(jù)視電阻率數(shù)據(jù)重組而來，因而繼承著視電阻率數(shù)據(jù)本身的特點(diǎn)——高阻區(qū)域?qū)?yīng)低阻異常、低阻區(qū)域?qū)?yīng)高阻異常，這點(diǎn)與反演成像有所區(qū)別。

2)近似成像圖上下邊界處存在明顯的形如“山峰”的假異常，“峰頂”指向異常體。這時(shí)由于成像過程中采用的近似方式引起的，在圖像解釋中應(yīng)注意區(qū)別。

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

3結(jié)語

作者采用有限單元法實(shí)現(xiàn)了井-井電阻率2.5維正演模擬，其中采用MSR法存儲(chǔ)超大型稀疏系數(shù)矩陣，運(yùn)用SSORPCG法求解線性方程組，保證了正演計(jì)算的效率和精度。根據(jù)地下點(diǎn)源電場的分布規(guī)律，提出了一種新穎的井-井電阻率2.5維近似成像方法，據(jù)此進(jìn)行了四種典型模型的近似成像。結(jié)果顯示，這里的井間電阻率近似成像方法能定性判斷異常體電阻率、較準(zhǔn)確反映異常體的位置、形態(tài)，并且具有極快的成像速度，因此可用于井-井觀測數(shù)據(jù)的快速預(yù)處理，得到井間剖面較為直觀的構(gòu)造。

參考文獻(xiàn)：

[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ìn)展[J]. 地球物理學(xué)進(jìn)展, 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]王興泰, 李曉芹. 電阻率圖像重建的佐迪反演及其應(yīng)用效果[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]. 地球物理學(xué)報(bào), 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導(dǎo)數(shù)的算法比較[J]. 物探化探計(jì)算技術(shù), 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]劉國強(qiáng), 陶果, 柯式鎮(zhèn), 等. 攝動(dòng)求解Frchet導(dǎo)數(shù)的電阻率層析成像方法[J]. 地質(zhì)地震, 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]吳小平, 徐果明. 電阻率三維反演中偏導(dǎo)數(shù)矩陣的求取和分析[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]. 物探化探計(jì)算技術(shù), 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]沈平, 強(qiáng)建科, 李永軍,等. 井間視電阻率的幾何成像方法[J]. 中南大學(xué)學(xué)報(bào)：自然科學(xué)版, 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]. 北京: 科學(xué)出版社, 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]. 地球物理學(xué)報(bào), 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]. 北京: 地質(zhì)出版社, 2005.

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

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

基金項(xiàng)目:國家“十二五”國家重大科學(xué)儀器設(shè)備開發(fā)專項(xiàng)(2012YQ030126)；廣西自然科學(xué)基金(2013GXNSFBA019212)；廣西高?？茖W(xué)技術(shù)研究重點(diǎn)項(xiàng)目(2013ZD029)

作者簡介：熊飛(1990-)男，碩士，主要從事電磁法數(shù)值模擬研究，E-mail：xiongfei357@163.com。

文章編號(hào)：1001-1749(2016)03-0308-06

中圖分類號(hào)：P 631.3

文獻(xiàn)標(biāo)志碼：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