常欣莉, 黨峰, 魏寶君, 陳濤, 王成園

(1.中國石油大學(xué)(華東)理學(xué)院, 山東 青島 266580; 2.中國石油集團(tuán)測(cè)井有限公司, 陜西 西安 710077)

0 引 言

井間電磁成像技術(shù)是將發(fā)射器置于發(fā)射井中,在一口或多口鄰近的井中接收電磁信號(hào)并通過對(duì)數(shù)據(jù)進(jìn)行反演成像從而獲得井間油氣藏分布的方法[1-3]。數(shù)值模擬算法研究是研發(fā)井間電磁成像儀器、分析井間電磁探測(cè)特性和對(duì)井間電磁測(cè)量數(shù)據(jù)進(jìn)行成像的前提。三維電磁場(chǎng)數(shù)值模擬算法的研究在20世紀(jì)70年代中期就已開始,積分方程算法[4-6]、有限元素法[7]和有限差分法[8-10]都得到了廣泛應(yīng)用。與需離散求解整個(gè)井間區(qū)域的有限元法和有限差分法相比,體積分方程法只需對(duì)三維異常體區(qū)域進(jìn)行離散,具有快速、方便等特點(diǎn),是計(jì)算三維井間電磁響應(yīng)的一種有效方法。Hohmann和Wannamaker等用體積分方程法計(jì)算了簡(jiǎn)單的三維模型,鮑光淑等[11]采用積分方程法計(jì)算出空間任一點(diǎn)的電場(chǎng)和磁場(chǎng),魏寶君等[12-13]提出了改進(jìn)型局域非線性迭代算法求解積分方程,徐凱軍等[14]利用積分方程法實(shí)現(xiàn)了均勻?qū)щ姲肟臻g三維大地電磁響應(yīng)的數(shù)值模擬。

傳統(tǒng)的關(guān)于體積分方程的應(yīng)用一般只限于各向同性地層[15-17],而陳桂波等[18]則將體積分方程用于各向異性地層,模擬了層狀各向異性地層中三維異常體的散射電磁波,并將層狀各向異性地層作為背景層,但該文獻(xiàn)沒有考慮傾斜井的情況。本文將體積分方程與層狀各向異性介質(zhì)中并矢Green函數(shù)的遞推矩陣方法[19-20]相結(jié)合模擬傾斜井及層狀各向異性地層中的井間電磁的響應(yīng)。體積分方程的求解通過弱化BCGS算法這一完全解法實(shí)現(xiàn),計(jì)算時(shí)將層狀各向異性原狀地層作為背景地層,從而將計(jì)算區(qū)域限制在三維異常體區(qū)域內(nèi),使未知量數(shù)目和計(jì)算量均大大減少,計(jì)算精度和效率同時(shí)得到保證。

1 基本理論

1.1 層狀各向異性地層中的體積分方程

所采用的井間地層模型見圖1。該模型包含了水平層狀各向異性原狀地層和電導(dǎo)率異常區(qū)域。每層背景地層的電導(dǎo)率張量均可表示為σb=diag(σhb,σhb,σvb),其中σhb為背景地層的水平電導(dǎo)率,σvb為背景地層的垂直電導(dǎo)率,不同層的電導(dǎo)率張量不同。將發(fā)射線圈視為磁偶極子并假設(shè)發(fā)射源隨時(shí)間的變化關(guān)系為exp(iωt),其中ω為角頻率,則接收井中接收線圈處的磁場(chǎng)強(qiáng)度矢量可表示為體積分方程形式

Δσ(r′)·E(r′,rT)dr′

(1)

式中,rT、rR分別為發(fā)射線圈和接收線圈的位置坐標(biāo);D為包含電導(dǎo)率異常體的積分區(qū)域;GHJ(rR,r′)為層狀各向異性背景地層中r′處的單位電流元在接收點(diǎn)rR處的磁型并矢Green函數(shù);Δσ(r′)=σ(r′)-σb(r′)=diag(Δσh,Δσh,Δσv)為異常體與背景地層電導(dǎo)率張量之差;E(r′,rT)為發(fā)射線圈在積分區(qū)域內(nèi)產(chǎn)生的總電場(chǎng);Hb(rR,rT)為發(fā)射源在層狀各向異性背景地層中產(chǎn)生的磁場(chǎng)強(qiáng)度。式(1)稱為數(shù)據(jù)方程,若積分區(qū)域內(nèi)的總電場(chǎng)E(r′,rT)已知,由式(1)可獲得接收線圈處的磁場(chǎng)強(qiáng)度。將水平層狀各向異性原狀地層作為背景地層,積分區(qū)域D就是電導(dǎo)率異常區(qū)域,從而使未知量數(shù)目和計(jì)算量均大大減少。

圖1 地層模型

在積分區(qū)域D內(nèi),總電場(chǎng)強(qiáng)度矢量E(r,rT)滿足

Δσ(r′)·E(r′,rT)dr′,r∈D

(2)

式中,GEJ(r,r′)為層狀各向異性背景地層中r′處的單位電流元在場(chǎng)點(diǎn)r處的電型并矢Green函數(shù);Eb(r,rT)為發(fā)射源在層狀各向異性背景地層中產(chǎn)生的電場(chǎng)強(qiáng)度。式(2)稱為目標(biāo)方程,該式是計(jì)算異常體內(nèi)總電場(chǎng)分布的第二類Fredholm積分方程。

GEJ在r=r′處存在強(qiáng)奇異性,在積分區(qū)域D內(nèi)不能直接求解式(2),而是根據(jù)GEJ與GAJ所滿足的關(guān)系

E(r,rT)=Eb(r,rT)+

Δσ(r′)·E(r′,rT)dr′r∈D

(3)

式中,μb為背景地層的磁導(dǎo)率,由于假設(shè)地層是非磁性的,故取μb為真空中的值μ0;GAJ(r,r′)為層狀各向異性背景地層中r′處的單位電流元在場(chǎng)點(diǎn)r處的磁矢勢(shì)并矢Green函數(shù)。由于GAJ的奇異性大大降低并且對(duì)其含奇異點(diǎn)的積分可以進(jìn)行弱化處理,故求解式(3)可避免奇異值的出現(xiàn)。

將積分區(qū)域D分成N個(gè)立方體形小單元Dj(j=1,2,…,N),其體積為ΔVj。假設(shè)每個(gè)小單元內(nèi)異常體與背景地層的電導(dǎo)率張量之差恒定,設(shè)為Δσj,并假設(shè)在每個(gè)小單元內(nèi)電場(chǎng)強(qiáng)度值恒定,由其中心點(diǎn)的值代替,則式(1)的離散形式為

Δσj·E(rj,rT)ΔVj

(4)

式(3)的離散形式為

E(ri,rT)=Eb(ri,rT)+

Δσj·E(rj,rT)ri∈D,i=1,2,…,N

(5)

式(5)可表示為線性形式

(6)

(7)

式(6)可采用穩(wěn)定型雙共軛梯度(BCGS)方法進(jìn)行迭代求解[21],從而得到總電場(chǎng)強(qiáng)度在異常區(qū)域D內(nèi)的分布。將采用上述算法得到的總電場(chǎng)強(qiáng)度E(ri,rT)代入式(4),即可得到接收線圈處的磁場(chǎng)強(qiáng)度在地層直角坐標(biāo)系中的各個(gè)分量。

根據(jù)發(fā)射線圈和接收線圈方向的不同,磁場(chǎng)強(qiáng)度各分量可表示為矩陣形式

(8)

H∥(rR,rT)=Hxxcosα1cosα2+Hyxcosα1cosβ2+Hzxcosα1cosγ2+Hxycosβ1cosα2+Hyycosβ1cosβ2+Hzycosβ1cosγ2+Hxzcosγ1cosα2+Hyzcosγ1cosβ2+Hzzcosγ1cosγ2

(9)

1.2 層狀各向異性背景地層的并矢Green函數(shù)

設(shè)發(fā)射線圈的磁偶極矩為MT,則式(1)和式(3)中發(fā)射源在背景地層中產(chǎn)生的電場(chǎng)強(qiáng)度和磁場(chǎng)強(qiáng)度可表示為

Eb(r,rT)=GEM(r,rT)·MT

(10)

Hb(rR,rT)=GHM(rR,rT)·MT

(11)

式中,GEM和GHM分別為水平層狀各向異性背景地層中單位磁偶極子源產(chǎn)生的電型和磁型并矢Green函數(shù)。魏寶君等[20]給出了計(jì)算層狀各向異性介質(zhì)中磁偶極子源和電偶極子源并矢Green函數(shù)的遞推矩陣方法,可快速精確地求解式(1)和式(3)所需要的GEM、GHM和GHJ。

(12)

若源點(diǎn)r′和場(chǎng)點(diǎn)ri不在同一地層中,則GAJ只包含散射項(xiàng)SGAJ。若源點(diǎn)r′和場(chǎng)點(diǎn)ri在同一地層中,GAJ除包含散射項(xiàng)外還包含背景項(xiàng)PGAJ。若r′和ri雖在同一地層中但i≠j,即ri不在小單元Dj內(nèi),則對(duì)PGAJ的積分不存在奇異值問題,可近似為

若i=j,則對(duì)PGAJ的積分存在奇異值問題,需進(jìn)行弱化處理。PGAJ各分量的具體表達(dá)式為

(13)

(14)

(15)

(16)

(17)

(18)

對(duì)于式(14)和式(15)的積分,由于被積函數(shù)在小單元Dj內(nèi)是相對(duì)于中心點(diǎn)rj的奇函數(shù),故有

(19)

2 數(shù)值算例及討論

采用體積分方程的弱化BCGS算法對(duì)層狀各向異性背景地層中的多分量井間電磁響應(yīng)進(jìn)行數(shù)值模擬,并分析發(fā)射井和接收井的傾角、背景地層的各向異性等因素對(duì)井間電磁響應(yīng)的影響。所采用的模型為3層各向異性地層模型(見圖1),上、下層電導(dǎo)率分別為σhb1=0.1 S·m-1、σvb1=0.1 S·m-1,σh3=0.1 S·m-1、σvb3=0.1 S·m-1。中間層σhb2=1.0 S·m-1、σvb2大小可變,厚度為10.0 m。上、下層界面垂向坐標(biāo)分別為z1=10 m、z2=20 m。在中間地層中有一個(gè)10 m×10 m×10 m的各向同性電導(dǎo)率異常體,其電導(dǎo)率為σa=0.02 S·m-1。電導(dǎo)率異常體在x、y、z方向劃分的單元數(shù)目均為5,小單元總數(shù)為125,每個(gè)小單元的尺寸為2 m×2 m×2 m。發(fā)射源置于發(fā)射井中垂向坐標(biāo)為15 m處,電導(dǎo)率異常體在井間居中。取發(fā)射源磁偶極矩為MT=1 A·m2,頻率1 kHz。取發(fā)射源與第1個(gè)接收點(diǎn)的井間橫向距離為100 m,若接收井傾斜,井間距則隨接收點(diǎn)垂向位置的增加而增加。

2.1 傾角對(duì)井間電磁響應(yīng)的影響

圖2 接收井傾斜時(shí)散射磁場(chǎng)軸向分量隨接收點(diǎn)垂向位置的變化

圖3 接收井傾斜時(shí)總磁場(chǎng)軸向分量隨接收點(diǎn)垂向位置的變化

假設(shè)背景層各向同性,即σvb2=1.0 S·m-1。取發(fā)射井和接收井在同一個(gè)xoz面上,即β1=β2=90 °。以發(fā)射井垂直、接收井傾斜為例,此時(shí)井間距隨接收點(diǎn)垂向位置的增加而增加。圖2和圖3分別給出了接收井傾斜角度不同時(shí)散射磁場(chǎng)和總磁場(chǎng)軸向分量隨接收點(diǎn)垂向坐標(biāo)的變化?？梢钥闯?傾斜井中磁場(chǎng)的軸向分量與傾斜角度及井間距密切相關(guān),傾斜角度及井間距的變化既改變了磁場(chǎng)峰值的位置也改變了峰值的大小。且傾角越大其虛分量指示地層界面的位置就越準(zhǔn)確。若接收井和發(fā)射井均傾斜,磁場(chǎng)響應(yīng)情況類似,此時(shí)磁場(chǎng)軸向分量變?yōu)?個(gè)分量的投影組合。

圖4 散射磁場(chǎng)xx分量隨接收點(diǎn)垂向位置的變化

圖5 總磁場(chǎng)xx分量隨接收點(diǎn)垂向位置的變化

圖6 散射磁場(chǎng)zz分量隨接收點(diǎn)垂向位置的變化

圖7 總磁場(chǎng)zz分量隨接收點(diǎn)垂向位置的變化

2.2 背景層各向異性對(duì)井間電磁響應(yīng)的影響

圖8 接收井傾斜30 °時(shí)散射磁場(chǎng)軸向分量隨接收點(diǎn)垂向位置的變化

通過改變?chǔ)襳b2的大小改變背景層電導(dǎo)率各向異性系數(shù),并模擬背景層電導(dǎo)率各向異性程度不同時(shí)的井間電磁響應(yīng),假設(shè)發(fā)射井和接收井均垂直于地層界面。圖4和圖5分別給出了背景層各向異性系數(shù)不同時(shí)散射磁場(chǎng)和總磁場(chǎng)的xx分量隨接收點(diǎn)垂向坐標(biāo)的變化。由圖4和圖5可以看出,xx分量散射磁場(chǎng)的強(qiáng)度與總磁場(chǎng)的強(qiáng)度相比要小得多,利用該分量不易于探測(cè)出電導(dǎo)率異常體的存在。另外,不同的各向異性系數(shù)對(duì)總磁場(chǎng)Hxx的實(shí)分量影響較大,且在實(shí)分量中出現(xiàn)的“犄角”可以明顯地指示地層界面的位置,而各向異性系數(shù)對(duì)Hxx虛分量的影響則相對(duì)較小。圖6和圖7分別給出了背景層各向異性系數(shù)不同時(shí)散射磁場(chǎng)和總磁場(chǎng)的zz分量隨接收點(diǎn)垂向坐標(biāo)的變化。由圖6和圖7可以看出,zz分量散射磁場(chǎng)的強(qiáng)度數(shù)值較大,利用該分量可探測(cè)電導(dǎo)率異常體的存在,并且該項(xiàng)受背景層各向異性系數(shù)的影響較大。當(dāng)然,在該模型中由于電導(dǎo)率異常體體積有限,磁場(chǎng)散射項(xiàng)強(qiáng)度仍小于背景項(xiàng)強(qiáng)度,所以總磁場(chǎng)受各向異性影響程度相對(duì)較小。

圖9 接收井傾斜30 °時(shí)總磁場(chǎng)軸向分量隨接收點(diǎn)垂向位置的變化

2.3 背景層各向異性和傾角對(duì)井間電磁響應(yīng)的影響

模擬背景層電導(dǎo)率各向異性和井斜同時(shí)存在時(shí)的井間電磁響應(yīng),仍假設(shè)發(fā)射井和接收井同在xoz平面上。首先考慮發(fā)射井垂直、接收井傾斜30 °時(shí)散射磁場(chǎng)和總磁場(chǎng)軸向分量隨接收點(diǎn)垂向坐標(biāo)的變化,計(jì)算結(jié)果見圖8和圖9。由于磁場(chǎng)背景項(xiàng)不受地層各向異性的影響,散射項(xiàng)雖受各項(xiàng)異性的影響但其強(qiáng)度與總磁場(chǎng)相比相對(duì)較弱,所以總磁場(chǎng)軸向分量受電導(dǎo)率各向異性影響仍然較小。然后考慮發(fā)射井和接收井均傾斜30 °的情況,即γ1=γ2=30 °,其總磁場(chǎng)軸向分量隨接收點(diǎn)垂向坐標(biāo)的變化如圖10所示。當(dāng)發(fā)射井和接收井同時(shí)傾斜時(shí),其軸向分量由Hxx、Hzx、Hxz、和Hzz這4個(gè)分量投影而成?？偞艌?chǎng)Hxx的實(shí)分量受地層電導(dǎo)率各向異性影響較大,在此情況下總磁場(chǎng)軸向分量受地層各向異性的影響較為明顯。

圖10 發(fā)射井和接收井均傾斜30 °時(shí)總磁場(chǎng)軸向分量隨接收點(diǎn)垂向位置的變化

3 結(jié) 論

(1) 將層狀各向異性介質(zhì)并矢Green函數(shù)的遞推矩陣方法與體積分方程的弱化BCGS算法有機(jī)結(jié)合,由于將計(jì)算區(qū)域限制在電導(dǎo)率異常體區(qū)域內(nèi),使未知量數(shù)目和計(jì)算量均大大減少,可有效模擬傾斜井及層狀各向異性地層中的井間電磁響應(yīng)。

(2) 利用上述算法的數(shù)值模擬結(jié)果表明,當(dāng)接收井發(fā)生傾斜時(shí),傾角越大磁場(chǎng)軸向分量的虛分量越能準(zhǔn)確指示地層界面的位置;xx分量散射磁場(chǎng)的強(qiáng)度與總磁場(chǎng)的強(qiáng)度相比要小得多,利用該分量不易于探測(cè)出電導(dǎo)率異常體的存在;不同的各向異性系數(shù)對(duì)總磁場(chǎng)Hxx影響較大,且在其實(shí)分量中出現(xiàn)的“犄角”可以明顯指示地層界面的位置;zz分量散射磁場(chǎng)的強(qiáng)度數(shù)值相對(duì)較大,利用該分量可探測(cè)電導(dǎo)率異常體的存在,并且該項(xiàng)受背景層各向異性系數(shù)的影響亦較大。

參考文獻(xiàn):

[1] Wilt M, Morrison H, Becker A, et al. Crosshole Electromagnetic Tomography: A New Technology for Oil Field Characterization [J]. The Leading Edge, 1995, 14(3): 173-177.

[2] Hoversten G, Newman G, Morrison H, et al. Reservoir Characterization Using Crosswell Electromagnetic Inversion: A Feasibility Study for the Snorre Field, North Sea [J]. Geophysics, 2001, 66(4): 1177-1189.

[3] 曾文沖, 趙文杰, 臧德福. 井間電磁成像系統(tǒng)應(yīng)用研究 [J]. 地球物理學(xué)報(bào), 2001, 44(3): 411-420.

[4] Hohmann G W. Three-dimensional Induced Polarization and Electromagnetic Modeling [J]. Geophysics, 1975, 40(2): 309-324.

[5] Wannamaker P E, Hohmann G W, SanFilipo W. Electromagnetic Modeling of Three-dimensional Bodies in Layered Earths Using Integral Equations [J]. Geophysics, 1984, 49(1): 60-74.

[6] Newman G A, Hohmann G W, Anderson W L. Transient Electromagnetic Response of a Three-dimensional Body in a Layered Earth [J]. Geophysics, 1986, 51(8): 1608-1627.

[7] Pridmore D F, Hohmann G W, Ward S H, et al. An Investigation of Finite-element Modeling for Electrical and Electromagnetic Data in Three Dimensions [J]. Geophysics, 1981, 46(7): 1009-1024.

[8] Alumbaugh D L, Newman G A, Prevost L, et al. Three-dimensional Wideband Electromagnetic Modeling on Massively Parallel Computers [J]. Radio Science, 1996, 31(1): 1-23.

[9] 沈金松. 用交錯(cuò)網(wǎng)格有限差分法計(jì)算三維頻率域電磁響應(yīng) [J]. 地球物理學(xué)報(bào), 2003, 46(2): 281-288.

[10] 沈金松. 用有限差分法計(jì)算各向異性介質(zhì)中多分量感應(yīng)測(cè)井的響應(yīng) [J]. 地球物理學(xué)進(jìn)展, 2004, 19(1): 101-107.

[11] 鮑光淑, 張碧星, 敬榮中. 三維電磁響應(yīng)積分方程法數(shù)值模擬 [J]. 中南工業(yè)大學(xué)學(xué)報(bào), 1999, 30(5): 472-474.

[12] 魏寶君, 張庚驥, 胡波. 用改進(jìn)的局域非線性迭代方法計(jì)算三維井間電磁場(chǎng) [J]. 石油大學(xué)學(xué)報(bào), 2002, 26(3): 124-127.

[13] 魏寶君, 張庚驥. 三維井間電磁場(chǎng)的正反演計(jì)算 [J]. 地球物理學(xué)報(bào), 2002, 45(5): 735-743.

[14] 徐凱軍, 李桐林, 張輝, 等. 利用積分方程法的大地電磁三維正演 [J]. 西北地震學(xué)報(bào), 2006, 28(2): 104-107.

[15] Wei B J, Simsek E, Liu Q H. Improved Diagonal Tensor Approximation (DTA) and Hybrid DTA/BCGS-FFT Method for Accurate Simulation of 3D Inhomogeneous Objects in Layered Media [J]. Waves in Random and Complex Media, 2007, 17(1): 55-66.

[16] Wei B J, Simsek E, Yu C, et al. Three-dimensional Electromagnetic Nonlinear Inversion in Layered Media by a Hybrid Diagonal Tensor Approximation: Stabilized Biconjugate Gradient Fast Fourier Transform Method [J]. Waves in Random and Complex Media, 2007, 17(2): 129-147.

[17] 魏寶君, LIU Q H. 水平層狀介質(zhì)中基于DTA的三維電磁波逆散射快速模擬算法 [J]. 地球物理學(xué)報(bào), 2007, 50(5): 1595-1605.

[18] 陳桂波, 汪宏年, 姚敬金, 等. 用積分方程法模擬各向異性地層中三維電性異常體的電磁響應(yīng) [J]. 地球物理學(xué)報(bào), 2009, 52(8): 2174-2181.

[19] Wei B J, Zhang G J, Liu Q H. Recursive Algorithm and Accurate Computation of Dyadic Green's Functions for Stratified Uniaxial Anisotropic Media [J]. Science in China: Series F, 2008, 51(1): 63-80.

[20] 魏寶君, 王甜甜, 王穎. 用磁流源并矢Green函數(shù)的遞推矩陣方法計(jì)算層狀各向異性地層中多分量感應(yīng)測(cè)井響應(yīng) [J]. 地球物理學(xué)報(bào), 2009, 52(11): 2920-2928.

[21] 魏寶君, LIU Q H. 層狀介質(zhì)中計(jì)算體積分方程的弱化BCGS-FFT算法 [J]. 中國石油大學(xué)學(xué)報(bào): 自然科學(xué)版, 2007, 31(1): 49-55.