井間視電阻率的幾何成像方法

沈平，強(qiáng)建科，李永軍，阮百堯

(1. 中南大學(xué) 信息物理工程學(xué)院，湖南 長沙，410083；2. 華北科技學(xué)院 安全工程學(xué)院，北京，065201；3. 桂林理工大學(xué) 資源與環(huán)境系，廣西 桂林，541004)

早在20世紀(jì)80年代末，井間電阻率成像就引起了國外地球物理學(xué)家們的極大興趣。Shima等[1]在1987年提出了“電阻率層析成像”概念，并研究了成像的可能性和反演算法。之后，Daily等[2]用拉普拉斯方程的非線性反演方法進(jìn)行了井間電阻率成像研究。Shima[3-4]提出了用有限元法模擬和α中心法聯(lián)合反演來重建地下電阻率分布圖像，獲得了較好的結(jié)果。國內(nèi)開展這方面的研究工作始于20世紀(jì)90年代中后期。周兵等[5-8]通過求取Jacobi矩陣的格林函數(shù)法和平滑度約束反演方法來實現(xiàn)2.5D井間電阻率成像，得到了較好的井間電阻率成像效果。吳小平等[9-10]將共軛梯度法用于三維電阻率成像，此法在每次反演迭代過程中僅需要1次正演計算，大大提高了計算速度。同時，在反演中將平滑度約束引用到最小二乘準(zhǔn)則中，有效地消除了多余構(gòu)造信息，可得到與實際情況非常接近的、穩(wěn)定可靠的反演結(jié)果。底青云等[11-14]提出了積分三維反演成像法，此法在計算上優(yōu)于有限單元法，與共軛梯度法一樣快捷，可以方便地進(jìn)行大電阻率差異背景下的成像。呂玉增等[15-16]在分析電性不均體對電位影響的基礎(chǔ)上，找到了1個直接成像的參數(shù)，但該參數(shù)的物理概念不是很明確，也沒有按照普遍的物理規(guī)則來運(yùn)算。綜上所述，現(xiàn)有井間成像方法對異常體的縱向位置定位精度較低，需要進(jìn)一步研究。

1 幾何成像方法

本文提出一種幾何成像方法。與傳統(tǒng)的井間電阻率成像方法相比，觀測方式是一致的，它們都是在一口鉆井中布置供電電極 A，另一口鉆井中布置測量電極M，對供電電極依次供電，觀測測量電極的電位，然后，將觀測電位換算成視電阻率。其觀測系統(tǒng)如圖1所示。

圖1 觀測系統(tǒng)示意圖Fig.1 Sketch map of observation system

該方法與傳統(tǒng)的井間成像方法的不同之處在于成像的原理不同。傳統(tǒng)的成像方法借助于雅可比矩陣求出供電點與測量點之間的電位響應(yīng)關(guān)系，從理論上講，這是一種比較精確的方法，但實際反演效果并不理想[17-18]，其原因是反演效果取決于約束條件和迭代次數(shù)。

基于靜電場原理，井間電測方法屬于“透射”工作方式，當(dāng)異常體靠近供電源時，異常體吸引或排斥電流強(qiáng)烈，電場分布畸變嚴(yán)重，測量點處能夠反映這種變化；根據(jù)電場互換原理，當(dāng)異常體靠近接收點時，異常體的電性特征也會嚴(yán)重影響測量點處的電位。為了在空間上充分展現(xiàn)井間電場的這種分布特征，提出井間幾何成像思想：首先，利用鏡像法求出與所有電位對應(yīng)的視電阻率，然后，以觀測井深為縱坐標(biāo)，供電井深為橫坐標(biāo)繪制視電阻率等值線圖。在等值線圖上，那些圈閉就是異常體的位置，這種圖示反映了供電點和測量點之間的映射關(guān)系。

另外，由于點電源在非半空間比在地表(半空間)時的電流密度要小，因此，井間二極觀測的異常較小。井間電阻率測量與地表測量不同，地表觀測是旁側(cè)影響測量，即在異常體上方測量，遇低阻體時，由于吸引電流而使觀測點電位降低，但井間測量則相反。井間測量是“對測” 測量，或者說是“透射”測量，遇低阻時由于吸引電流使測量點較其他“等源距”位置上的電位增加，由此計算的視電阻率呈高阻反映。這個結(jié)論與模型的實際情況正好相反，數(shù)值模擬結(jié)果也證實了這一點。因此，在新的成像方法中，高視電阻率異常為低阻體所致，低視電阻率異常為高阻體所致，而且由視電阻率異常形態(tài)能夠判斷異常體的空間形態(tài)和精確位置。在測量過程中，當(dāng)供電電極和測量電極位于同一深度時，距離AM最短，電場最強(qiáng)，如果其間正好存在異常體，則對電流場的擾動最強(qiáng)烈，在此處計算的視電阻率變化最大，這就是該方法縱向分辨率高的原因。

下面利用鏡像法求二極裝置視電阻率。當(dāng)在一口井中供電，另一口井中測量，即供電點A和觀測點M分別位于2口井中時，地面邊界會影響點源A的電流分布，從而影響M點的電位。為了消除地面邊界的影響，采用鏡像法計算M點處的視電阻率。

二極裝置計算示意圖如圖2所示。假設(shè)地面上方與A對稱的位置有1個虛電源A′，供電點A與測量點M之間的距離為r，虛電源A′與測量點M之間的距離為r′，x方向與地面平行，z軸垂直于地面向下。由圖2可知：

圖2 二極裝置計算示意圖Fig.2 Calculation sketch map of two-pole array

其中：xM和zM分別為M點在x方向和z方向的坐標(biāo)；xA和zA分別為A點在x方向和z方向的坐標(biāo)。

由直流點電源公式可得觀測點M處的電位為：

其中：ρ為視電阻率；I為電源A的電流。整理式(3)，可得：

2 三維數(shù)值模擬

在三維空間直角坐標(biāo)系中，穩(wěn)定電流場點電源電位滿足基本方程[19-20]：

其中：σ為電導(dǎo)率；u為電位；δ(A)為以A為中心的δ函數(shù)。

三維電場的邊值問題為：

其中：n為電位的外法向；r為電源點至邊界的距離。

等價的變分問題是

其中：∞Γ為無窮遠(yuǎn)邊界；Ω為求解區(qū)域。

有限單元法是求解上述三維問題的有效方法。本文選用六面體為網(wǎng)格剖分單元，電位三線性插值，電導(dǎo)率分塊均勻。有限單元法計算見文獻(xiàn)[19]。

下面計算幾種典型模型來說明這種幾何成像的效果和應(yīng)注意的問題。

模型1：2井間距為100 m，井深1 km，井間的中上部有1個低阻體，其長×寬×高為20 m×40 m×80 m，頂部埋深為300 m，離左右井的距離均為40 m。圍巖電阻率 ρ0=50 Ω·m，低阻體電阻率 ρ1=8 Ω·m，井中電阻率為3 Ω·m，其剖面圖如圖3所示。

圖3 模型1的剖面圖Fig.3 Profile of model one

圖4 所示為低阻長方體的計算結(jié)果，其中：縱坐標(biāo)為觀測井深，橫坐標(biāo)為供電井深(下同)，淺色區(qū)域為標(biāo)準(zhǔn)模型位置。異常形態(tài)向右傾斜 45°，且呈橢圓狀。由于模型位于2井之間的中間部位，當(dāng)在300 m處供電時，接收點在380 m處得到最大異常；在340 m處供電時，接收點在340 m處得到最大異常，即此時A與M受異常影響最大；在380 m處供電時，接收點在300 m處得到最大異常。因此，異常形態(tài)向右傾斜。

模型2：2井間距為100 m，井深1 km，井間的中上部有1個低阻層，其寬×高為30 m×40 m，頂部埋深為220 m。圍巖電阻率ρ0=50 Ω·m，低阻層電阻率ρ1=15 Ω·m，井中電阻率為 3 Ω·m，其剖面圖如圖 5所示。

圖4 三維低阻體的視電阻率等值線圖Fig.4 Apparent resistivity isopleth map of 3D low-resistivity body

圖5 模型2的剖面圖Fig.5 Profile of model two

圖6 所示為水平低阻層的視電阻率等值線圖，其異常形態(tài)近似為“十字星型”。當(dāng)供電點 A與測量點M位于低阻層的上下兩側(cè)時，測量點處的異常比較大。特別當(dāng)供電點A位于220～260 m處，且測量點M也位于220～260 m處時，能夠得到最大異常，按照成圖方式就形成了“十字星型”異常。

模型3：2井間距為100 m，井深1 km，井間的右下部有1個高阻體，其長×寬×高為20 m×40 m×80 m，頂部埋深為700 m，距右井的距離為10 m。圍巖電阻率 ρ0=50 Ω·m，高阻體電阻率 ρ1=800 Ω·m，井中電阻率為3 Ω·m，其剖面圖如圖7所示。

圖8所示為埋深為700 m的高阻模型的計算結(jié)果。如前所述，視電阻率異常中心呈低阻反映，其等值線為水平扁橢圓，原因是高阻體距離測量井太近(10 m)，也就是說大多數(shù)供電點都能在測量井的700～780 m處得到明顯異常。

圖6 水平低阻層的視電阻率等值線圖Fig.6 Apparent resistivity isopleth map of horizontal low-resistivity layer

圖7 模型3的剖面圖Fig.7 Profile of model three

圖8 三維高阻體的視電阻率等值線圖Fig.8 Apparent resistivity isopleth map of 3D high-resistivity body

模型4：2井間距為100 m，井深1 km，井間的中上部和右下部分別有1個低阻層和1個高阻體。低阻層寬×高為30 m×40 m，其頂部埋深為140 m；高阻體長×寬×高為20 m×40 m×80 m，其頂部埋深為700 m，距右井的距離為10 m。圍巖電阻率為50 Ω·m，低阻層電阻率為10 Ω·m，高阻體電阻率為1 kΩ·m，井中電阻率為3 Ω·m。

圖9所示為模型4的模擬結(jié)果。從圖9可以看出：“高阻十字星型”表明是1個水平低阻層，“水平扁橢圓”表明是1個靠近測量井的三維高阻體。組合異常清楚可靠，互不干擾。

圖9 組合模型的視電阻率等值線圖Fig.9 Apparent resistivity isopleth map of combination model

模型5：兩井間距為100 m，井深1 km，井間深度為320～400 m處有2個模型：左邊是1個長×寬×高為20 m×40 m×80 m的高阻體，電阻率為1 kΩ·m，距左井距離為 10 m；右邊是 1個長×寬×高為20 m×40 m×80 m的低阻體，電阻率為10 Ω·m，距右井距離為10 m；兩異常體相距40 m，圍巖電阻率為 50 Ω·m，井中電阻率為 3 Ω·m。

圖10所示為同一深度2個異常體的視電阻率等值線圖。從圖10可以看出：總體異常形態(tài)呈“十字”相交；當(dāng)供電點A位于高阻體附近時，在低阻體附近處上下300 m范圍內(nèi)都有較強(qiáng)異常；同樣，當(dāng)接收點M位于低阻體附近時，在供電點高阻體附近上下300 m范圍內(nèi)也有較強(qiáng)異常。盡管兩異常相交，但仍然能夠清晰區(qū)分。

從對上述模型1～5的計算結(jié)果可以看出：本文所采用的井間幾何成像效果較好。由于是非半空間，電流密度較小，視電阻率異常差異不是太大。在視電阻率等值線圖中，地質(zhì)模型位置清楚、準(zhǔn)確，水平低阻層的異常形態(tài)近似為“十字星型”，直立三維高阻體的異常形態(tài)近似為“水平扁橢圓”。但是，在所成圖像中，高視電阻率反映低阻體，低視電阻率反映高阻體，與地面電法結(jié)論相反，必須引起注意。

圖10 同一深度2個異常體的視電阻率等值線圖Fig.10 Apparent resistivity isopleth map of two abnormal body at the same depth

3 結(jié)論

(1) 提出一種井間幾何成像方法，即采用鏡像法原理計算非全空間視電阻率，以測量井為縱坐標(biāo)，供電井為橫坐標(biāo)繪制視電阻率等值線圖，成像效果較好，消除了地面影響。

(2) 幾何成像方法縱向分辨率較高，橫向分辨率較低，能準(zhǔn)確確定異常體中心埋深。

(3) 存在的問題是電阻率為視電阻率，目前不能求出真電阻率，這有待進(jìn)一步研究。

[1] Shima H, Sakayama T. Resistivity tomography: An approach to 2-D resistivity inverse problems[C]//57th Society of Exploration Geophysicists, Expanded Abstracts. Beijing: Petroleum Industry Press, 1989: 204-207.

[2] Daily W, Owen E. Gross-borehole resistivity tomography[J].Geophysics, 1991, 56(8): 1228-1235.

[3] Shima H. Two-dimensional automatic resistivity inversion technique using alpha centers[J]. Geophysics, 1990, 55(6):682-694.

[4] Shima H. 2D and 3D resistivity image reconstruction using cross-hole data[J]. Geophysics, 1992, 57(10): 1270-1281.

[5] 周兵, 曹俊興. 井間電阻率成像數(shù)值模擬[J]. 物探化探計算技術(shù), 1995, 17(4): 9-17.ZHOU Bing, CAO Jun-xing. Numerical modeling for crosshole resisitivity image[J]. Computing Techniques for Geophysical and Geochemical Exploration, 1995, 17(4): 9-17.

[6] 董清華, 田憲謨, 嚴(yán)忠瓊. 井間電阻率層析成象的平滑度約束反演[J]. 物探化探計算技術(shù), 1997, 19(2): 122-127.DONG Qing-hua, TIAN Xian-mo, YAN Zhong-qiong. A smothness-constrained inversion for crosshole resistivity tomography[J]. Computing Techniques for Geophysical and Geochemical Exploration, 1997, 19(2): 122-127.

[7] 董清華, 曹俊興. 井間電阻率層析成像的幾個問題研究[J].地球物理學(xué)進(jìn)展, 1998, 13(4): 84-89.DONG Qing-hua, CAO Jun-xing. Study on some problem for crosshole resistivity tomography[J]. Progress in Geophysics,1998, 13(4): 84-89.

[8] 董清華. 井間電阻率層析成像的某些進(jìn)展[J] . 地球物理學(xué)進(jìn)展, 1997, 12(3): 77-89.DONG Qing-hua. The latest development of cosshole resisitivity tomography[J]. Progress in Geophysics, 1997, 12(3): 77-89.

[9] 吳小平, 徐果明. 利用共軛梯度法的電阻率三維反演研究[J].地球物理學(xué)報, 2000, 43(3): 420-427.WU Xiao-ping, XU Guo-ming. Study on 3-D resistivity inversion using conjugate gradient method[J]. Chinese Journal of Geophysics, 2000, 43(3): 420-427.

[10] 黃俊革, 阮百堯, 鮑光淑. 基于有限單元法的三維地電斷面電阻率反演[J]. 中南大學(xué)學(xué)報: 自然科學(xué)版, 2004, 35(2):295-299.HUANG Jun-ge, RUAN Bai-yao, BAO Guang-shu. Resistivity inversion on 3-D section based on FEM[J]. Journal of Central South University: Science and Technology, 2004, 35(2):295-299.

[11] 底青云, 王妙月. 積分法三維電阻率成像[J] . 地球物理學(xué)報,2001, 44(6): 844-853.DI Qing-yun, WANG Miao-yue. 3-D resistivity tomography by integral method[J]. Chinese Journal of Geophysics, 2001, 44(6):844-853.

[12] 底青云, 王妙月. 電流線追蹤電位電阻率層析成像方法初探[J]. 地球物理學(xué)進(jìn)展, 1997, 12(4): 27-34.DI Qing-yun, WANG Miao-yue. The resistivity tomography by tracing the electric ray[J]. Progress in Geophysics, 1997, 12(4):27-34.

[13] 白登海, 于晟. 電阻率層析成像理論和方法[J]. 地球物理學(xué)進(jìn)展, 1995, 10(1): 56-75.BAI Deng-hai, YU Sheng. Theory and methods of resistivity tomography[J]. Progress in Geophysics, 1995, 10(1): 56-75.

[14] 楊華, 李金銘, 桂峰. 電阻率層析成像方法技術(shù)近年發(fā)展概況[J]. 地球物理學(xué)進(jìn)展, 1998, 13(4): 90-96.YANG Hua, LI Jin-ming, GUI Feng. The central situation of development for methods and technologies on resistivity tomography in recent years[J]. Progress in Geophysics, 1998,13(4): 90-96.

[15] 呂玉增, 阮百堯, 黃俊革. 直流電井間三維直接成像[J]. 物探化探計算技術(shù), 2003, 25(1): 60-64.Lü Yu-zeng, RUAN Bai-yao, HUANG Jun-ge. The 3-D immediate crosshole tomography with direct current[J].Computing Techniques for Geophysical and Geochemical Exploration, 2003, 25(1): 60-64.

[16] 曾文沖, 趙文杰, 臧德福. 井間電磁成像系統(tǒng)研究[J]. 地球物理學(xué)報, 2001, 44(3): 411-420.ZENG Wen-chong, ZHAO Wen-jie, ZANG De-fu. Application research of crosshole electromagnetic tomography[J]. Chinese Journal of Geophysics, 2001, 44(3): 411-420.

[17] 董清華, 田憲謨. 井間電阻率成像中Frèchet導(dǎo)數(shù)的算法比較[J]. 物探化探計算技術(shù), 1997, 19(1): 41-45.DONG Qing-hua, TIAN Xian-mo. A comparision of methods for calculating Frèchet derivatives in crosshole resistivity image[J].Computing Techniques for Geophysical and Geochemical Exploration, 1999, 19(1): 41-45.

[18] 董清華, 嚴(yán)忠瓊. 井間電阻率成像在工程勘察中的應(yīng)用[J].工程勘察, 1998(1): 73-75.DONG Qing-hua, YAN Zhong-qiong. The application of crosshole resistivity tomography in engineering prospecting[J].Engineering Prospecting, 1998(1): 73-75.

[19] 強(qiáng)建科, 羅延鐘. 三維地形直流電阻率有限元法模擬[J]. 地球物理學(xué)報, 2007, 50(5): 1606-1613.QIANG Jian-ke, LUO Yan-zhong. The resistivity FEM numerical modeling on 3-D undulating topography[J]. Chinese Journal of Geophysics, 2007, 50(5): 1606-1613.

[20] 強(qiáng)建科, 何繼善. 橢球體上雙頻激電法的正演與反演算法[J].中南大學(xué)學(xué)報: 自然科學(xué)版, 2007, 38(6): 1199-1205.QIANG Jian-ke, HE Ji-shan. Algorithm of forward and inversion of dual-frequency induced polarization method on elliptical sphere[J]. Journal of Central South University: Science and Technology, 2007, 38(6): 1199-1205.