井-井電阻率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