利用徑向成層介質(zhì)的Green函數(shù)和積分方程模擬含金屬套管井間電磁場的響應(yīng)

魏寶君,陳 濤,侯學(xué)理,馮琳偉,黨 峰,陳章龍

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

井間電磁探測技術(shù)由于具有更大的橫向探測范圍而成為油氣藏研究的重要手段之一，數(shù)值模擬方法和成像方法是該技術(shù)的重要組成部分。傳統(tǒng)的井間電磁數(shù)值模擬方法和成像方法研究均基于無金屬套管的井眼環(huán)境[1-6]，而在很多情況下井間電磁測量可能在發(fā)射井或接收井中有金屬套管的井眼環(huán)境下進(jìn)行[7]，其測量數(shù)據(jù)與無金屬套管環(huán)境下的測量數(shù)據(jù)存在較大差別，因而有必要研究金屬套管對井間電磁場的影響。筆者將求解積分方程的BCGS算法[8-9]與徑向成層介質(zhì)Green函數(shù)的遞推矩陣方法[10]相結(jié)合，將金屬套管作為一層特殊介質(zhì)包含到Green函數(shù)中，計算區(qū)域?yàn)榫g地層，模擬發(fā)射井含金屬套管時軸對稱二維井間電磁場的響應(yīng)，分析金屬套管對井間電磁場的影響規(guī)律。

1 積分方程

所采用的軸對稱二維井間地層模型見圖1，該模型包含了井眼、發(fā)射井金屬套管、井間均質(zhì)背景地層和電導(dǎo)率異常區(qū)域。采用圓柱坐標(biāo)系，忽略位移電流的影響并假設(shè)地層是非磁性的，其磁導(dǎo)率取真空中的數(shù)值μ0，均質(zhì)背景地層的電導(dǎo)率設(shè)為σb。將發(fā)射線圈視為磁偶極子并假設(shè)發(fā)射源隨時間的變化關(guān)系為exp(iωt)，其中ω為角頻率，則接收井中接收線圈處磁感應(yīng)強(qiáng)度的軸向分量可表示為如下積分方程的形式[5]：

式中，rT、rR分別為發(fā)射線圈和接收線圈的位置坐標(biāo)，D為電導(dǎo)率異常區(qū)域；(rR，r′)為徑向成層地層r′處的單位圓電流元在接收點(diǎn)rR處的Green函數(shù)；Δσ(r′)為異常區(qū)域與背景地層電導(dǎo)率之差，Δσ(r′)=σ(r′)-σb(r′)；A(r′，rT)為發(fā)射線圈在電導(dǎo)率異常區(qū)域內(nèi)產(chǎn)生的總磁矢勢(只有φ分量)；Bb(rR，rT)=μ0NTIT(rR，rT)為發(fā)射線圈在徑向成層地層中產(chǎn)生的磁感應(yīng)強(qiáng)度的軸向分量，NT為發(fā)射線圈匝數(shù)，IT為發(fā)射線圈電流強(qiáng)度。式(1)稱為數(shù)據(jù)方程，若電導(dǎo)率異常區(qū)域內(nèi)的總磁矢勢A(r′，rT)已知，由該式可獲得接收線圈處磁感應(yīng)強(qiáng)度的軸向分量。

圖1 含套管軸對稱二維井間地層模型簡圖Fig.1 Schematic configuration of axial symmetric 2-D cross-hole formation model with metal casing

在異常區(qū)域D內(nèi)，總磁矢勢A(r，rT)滿足如下積分方程[5]：

式中，Ab為磁矢勢的背景值；Γ(r，r′)為徑向成層地層磁矢勢Green函數(shù)。

將電導(dǎo)率異常區(qū)域D分成N個小單元Dj，每個小單元內(nèi)地層電導(dǎo)率恒定，設(shè)為σj，并假設(shè)在每個小單元內(nèi)磁矢勢恒定，并由其中心點(diǎn)的值代替，則式(1)可離散為

式(4)可表示為如下線性形式：

其中I為線性算子，定義為

式(6)可采用穩(wěn)定型雙共軛梯度(BCGS)方法進(jìn)行迭代求解[8-9]，從而得到總磁矢勢在區(qū)域D內(nèi)的分布。求解式(6)時先對所有小單元給定任意初始猜測值A(chǔ)0(例如A0=0)，并計算下列殘差：

選擇任一量0滿足(0，r0)≠0，比如0=r0，第k步的迭代順序?yàn)?/p>

當(dāng)小于給定的誤差標(biāo)準(zhǔn)時，迭代終止。

2 徑向成層地層的Green函數(shù)及其對空間坐標(biāo)的積分

2.1 徑向成層地層的Green函數(shù)

在上述積分方程表達(dá)式中，含有對徑向成層地層 Green 函數(shù)(rR，rT)、Γ(r，rT) 及積分的計算，可采用文獻(xiàn)[10] 的遞推矩陣方法。設(shè)徑向成層地層共有n+1層，由里向外的編號為l=0，1，…，n，各層半徑為r0，r1，…，rn-1，單位圓電流發(fā)射源在第j層，各層電參數(shù)分別為μl、εl。任意第l層內(nèi)的 Green函數(shù)Γl(r，r′) 可表示為如下形式：

式中，ν和ν分別為扣除指數(shù)項(xiàng)后的第一類和第二類ν階復(fù)宗量變型Bessel函數(shù)。將變型Bessel函數(shù)的指數(shù)項(xiàng)單獨(dú)列出計算可避免上溢現(xiàn)象[10]。其中，Al(λ)和Bl(λ)為待定系數(shù)，由圓柱形層界面處電場和磁場的連續(xù)性條件確定，可采用遞推矩陣方法快速求解[5，10]。

由式(7)利用變型Bessel函數(shù)對徑向坐標(biāo)的微分關(guān)系可得到l(r，r′)，若r＞r′，有

利用式(8)可計算井間無電導(dǎo)率異常體時接收線圈處磁感應(yīng)強(qiáng)度的軸向分量Bb(rR，rT)。

2.2 Green函數(shù)對空間坐標(biāo)的積分

設(shè)某一矩形小單元De中心點(diǎn)的坐標(biāo)為(rc，zc)、其徑向和軸向?qū)挾确謩e為Δr、Δz，在該單元內(nèi)對Green 函數(shù)Γl(r，r′) 積分，有

待定系數(shù)由下列方程得到：

其中，，矩陣M的表達(dá)式見文獻(xiàn)[10]。結(jié)合文獻(xiàn)[10]經(jīng)推導(dǎo)可得到的各非零元素為

在該單元內(nèi)對 Green 函數(shù)l(r，r′) 積分，有

求解式(11)亦采用遞推矩陣方法。

3 數(shù)值模擬舉例

3.1 有電導(dǎo)率異常區(qū)塊時的井間電磁場

設(shè)發(fā)射源磁偶極矩MT=1 Am2、井間距為100 m、金屬套管的內(nèi)外半徑分別為8 cm和10 cm、金屬套管的電導(dǎo)率為5×105S/m、發(fā)射井內(nèi)鉆井液電導(dǎo)率為0.05 S/m、發(fā)射頻率為100 Hz。模型1為井間存在一個電導(dǎo)率異常區(qū)塊的情況。固定發(fā)射源垂向坐標(biāo)為15 m，電導(dǎo)率異常區(qū)塊的范圍為徑向30～70 m、垂向10～50 m，小單元數(shù)目為20×20=400，每個小單元的尺寸為2 m×2 m。圖2給出了井間背景地層電導(dǎo)率為0.2 S/m、異常區(qū)塊電導(dǎo)率為2.0 S/m時有金屬套管和無金屬套管情況下的磁感應(yīng)強(qiáng)度隨接收線圈垂向坐標(biāo)zR的變化關(guān)系。圖3給出了井間背景地層電導(dǎo)率為1.0 S/m、異常區(qū)塊電導(dǎo)率為0.02 S/m時有金屬套管和無金屬套管情況下的磁感應(yīng)強(qiáng)度隨接收線圈垂向坐標(biāo)的變化關(guān)系。

圖2 模型1第一種情況下有金屬套管和無金屬套管時總磁感應(yīng)強(qiáng)度隨接收線圈垂向坐標(biāo)的變化關(guān)系Fig.2 Relationship between total magnetic induction intensity with or without metal casing and receiver coil's vertical position for the first case of model 1

圖3 模型1第二種情況下有金屬套管和無金屬套管時總磁感應(yīng)強(qiáng)度隨接收線圈垂向坐標(biāo)的變化關(guān)系Fig.3 Relationship between total magnetic induction intensity with or without metal casing and receiver coil's vertical position for the second case of model 1

由上述算例可以看出，有金屬套管時實(shí)部磁感應(yīng)強(qiáng)度的數(shù)值要小于無金屬套管時的磁感應(yīng)強(qiáng)度值，而有金屬套管時虛部磁感應(yīng)強(qiáng)度的數(shù)值要大于無金屬套管時的磁感應(yīng)強(qiáng)度值。計算發(fā)現(xiàn)，對磁感應(yīng)強(qiáng)度的總幅度值而言，由于金屬套管對電磁信號的吸收，有金屬套管時磁感應(yīng)強(qiáng)度的總幅度值要小于無金屬套管時的磁感應(yīng)強(qiáng)度的總幅度值。另由圖2和圖3發(fā)現(xiàn)，當(dāng)有金屬套管存在時除實(shí)、虛部磁感應(yīng)強(qiáng)度的數(shù)值發(fā)生變化外，其響應(yīng)曲線沿垂向坐標(biāo)出現(xiàn)“平移”，即兩種情況下磁感應(yīng)強(qiáng)度的極值點(diǎn)不再對應(yīng)接收線圈的同一垂向坐標(biāo)，且虛分量的“平移”現(xiàn)象更加明顯。這說明金屬套管的存在還導(dǎo)致了磁感應(yīng)強(qiáng)度相位的改變。

模型2為井間存在兩個電導(dǎo)率異常區(qū)塊的情況，取井間背景地層電導(dǎo)率值為1.0 S/m，固定發(fā)射源垂向坐標(biāo)為40 m。第一個電導(dǎo)率異常區(qū)塊的范圍為徑向30～70 m、垂向10～30 m，電導(dǎo)率為5.0 S/m。第二個電導(dǎo)率異常區(qū)塊的范圍為徑向30～70 m、垂向50～70 m，電導(dǎo)率為0.02 S/m。每個小單元的尺寸仍為2 m×2 m。圖4給出了該模型有金屬套管和無金屬套管時的磁感應(yīng)強(qiáng)度隨接收線圈垂向坐標(biāo)的變化關(guān)系。由該算例可得到類似結(jié)論。

圖4 模型2有金屬套管和無金屬套管時總磁感應(yīng)強(qiáng)度隨接收線圈垂向坐標(biāo)的變化關(guān)系Fig.4 Relationship between total magnetic induction intensity with or without metal casing and receiver coil's vertical position for model 2

3.2 井間地層電導(dǎo)率改變時金屬套管的影響

假設(shè)井間無電導(dǎo)率異常塊，并假設(shè)發(fā)射線圈和接收線圈具有相同的垂向坐標(biāo)，其余參數(shù)同上。圖5和圖6給出了有金屬套管和無金屬套管兩種情況下頻率分別為50 Hz和500 Hz時接收線圈處的磁感應(yīng)強(qiáng)度隨井間均質(zhì)地層電導(dǎo)率的變化關(guān)系。由圖5和圖6可以看出，在頻率為50 Hz情況下，金屬套管的影響較小，當(dāng)頻率為500 Hz時，金屬套管的影響明顯增大。就磁感應(yīng)強(qiáng)度的幅度而言，由于金屬套管的存在導(dǎo)致穿過套管的電磁信號變?nèi)?，故有金屬套管情況下磁感應(yīng)強(qiáng)度的幅度變小。但金屬套管對磁感應(yīng)強(qiáng)度幅度值的影響是近似恒定的，即對于不同的地層電導(dǎo)率，有金屬套管和無金屬套管情況下磁感應(yīng)強(qiáng)度幅度的比值近似恒定。這個比值取決于金屬套管的參數(shù)和頻率，與地層電導(dǎo)率無關(guān)。就磁感應(yīng)強(qiáng)度的實(shí)分量和虛分量而言，變化關(guān)系要復(fù)雜，主要是存在曲線的“平移”且虛分量更加明顯，這同樣說明金屬套管亦引起磁感應(yīng)強(qiáng)度相位的變化。

圖5 頻率為50 Hz時磁感應(yīng)強(qiáng)度隨井間均質(zhì)地層電導(dǎo)率的變化關(guān)系Fig.5 Relationship between magnetic induction intensity and homogeneous cross-hole formation's conductivity for the frequency of 50 Hz

圖6 頻率為500 Hz時磁感應(yīng)強(qiáng)度隨井間均質(zhì)地層電導(dǎo)率的變化關(guān)系Fig.6 Relationship between magnetic induction intensity and homogeneous cross-hole formation's conductivity for the frequency of 500 Hz

3.3 金屬套管電導(dǎo)率對井間電磁信號的影響

假設(shè)井間均質(zhì)地層電導(dǎo)率為0.1 S/m，其余參數(shù)不變，圖7給出了頻率分別為50 Hz和500 Hz時接收線圈處的磁感應(yīng)強(qiáng)度隨金屬套管電導(dǎo)率σc的變化關(guān)系。由圖7可以看出，在套管電導(dǎo)率很小時，其影響可以忽略不計。當(dāng)套管電導(dǎo)率增大到某一臨界值時導(dǎo)致磁感應(yīng)強(qiáng)度快速變化。當(dāng)套管電導(dǎo)率增加到該臨界值時，磁感應(yīng)強(qiáng)度的實(shí)分量值快速減小，此后逐漸穩(wěn)定到一個較低數(shù)值；而磁感應(yīng)強(qiáng)度的虛分量先增加再快速減小，此后逐漸穩(wěn)定到一個數(shù)值。套管電導(dǎo)率的臨界值與頻率有關(guān)，頻率越高臨界值越小。

圖7 磁感應(yīng)強(qiáng)度隨套管電導(dǎo)率的變化關(guān)系Fig.7 Relationship between magnetic induction intensity and casing's conductivity

4 結(jié) 論

(1)金屬套管的存在導(dǎo)致井間電磁信號的強(qiáng)度變?nèi)?，且金屬套管對井間電磁信號強(qiáng)度的衰減作用隨頻率的升高而增大。

(2)金屬套管對磁感應(yīng)強(qiáng)度幅度的影響近似恒定，即有金屬套管和無金屬套管情況下磁感應(yīng)強(qiáng)度幅度的比值近似恒定。這個比值取決于金屬套管的參數(shù)和頻率，與地層電導(dǎo)率無關(guān)。

(3)當(dāng)有金屬套管存在時，井間磁感應(yīng)強(qiáng)度的實(shí)分量和虛分量除數(shù)值發(fā)生變化外，還存在“平移”現(xiàn)象且虛分量更加明顯，說明金屬套管亦引起井間磁感應(yīng)強(qiáng)度相位的變化。

(4)套管電導(dǎo)率很小時，其對井間電磁信號的影響可以忽略不計。當(dāng)套管電導(dǎo)率增大到某一臨界值時導(dǎo)致井間磁感應(yīng)強(qiáng)度快速變化。在該臨界值處磁感應(yīng)強(qiáng)度的實(shí)分量值快速減小，此后逐漸穩(wěn)定到一個較低數(shù)值，而虛分量則先增加后減小，此后逐漸穩(wěn)定到一個數(shù)值。該臨界值與頻率有關(guān)，頻率越高臨界值越小。

[1] ALUMBAUGH D L，MORRISON H F.Theoretical and practical considerations for cross-well electromagnetic tomography assuming a cylindrical geometry[J].Geophysics，1995，60(3)：846-870.

[2] LIU C，SHEN L C，ZHANG Y C.Monitoring saltwater injection using conductivity images obtained by electromagnetic cross-hole measurements[J].Radio Science，1995，30(5)：1405-1415.

[3] ZHANG Y C，LIU C，SHEN L C.An iterative algorithm for conductivity image reconstruction from crosswell EM measurements[J].Geophysics，1996，61(4)：987-997.

[4] 魏寶君.三維井間電磁場的聯(lián)合反演方法[J].石油大學(xué)學(xué)報：自然科學(xué)版，2005，29(3)：20-26.WEI Bao-jun.Combined inversion algorithm of three-dimensional cross-hole electromagnetic fields[J].Journal of the University of Petroleum，China(Edition of Natural Science)，2005，29(3)：20-26.

[5] 魏寶君.井間電磁場的一維、二維聯(lián)合反演方法[J].地球物理學(xué)報，2006，49(1)：264-274.WEI Bao-jun.A combined 1-D/2-D inversion algorithm of cross-hole electromagnetic fields[J].Chinese Journal of Geophysics，2006，49(1)：264-274.

[6] 沈金松，孫文博，趙文杰，等.井間電磁測量數(shù)據(jù)的模擬與處理解釋[J].測井技術(shù)，2007，31(5)：417-424.SHEN Jin-song，SUN Wen-bo，ZHAO Wen-jie，et al.Modeling and processing of the cross-hole electromagnetic and measurements and their interpretation[J].Well Logging Technology，2007，31(5)：417-424.

[7] BEER R，ZHANG P，ALUMBAUGH D，et al.Anisotropy modeling and inversions of DeepLook-EM data from Brazil[C/CD].2010 SEG Annual Meeting，October 17-22，2010，Denver，Colorado，USA.

[8] 魏寶君，LIU Q H.水平層狀介質(zhì)中基于DTA的三維電磁波逆散射快速模擬算法[J].地球物理學(xué)報，2007，50(5)：1595-1605.WEI Bao-jun，LIU Q H.Fast algorithm for simulating 3-D electromagnetic inverse scattering in horizontally stratified medium via DTA[J].Chinese Journal of Geophysics，2007，50(5)：1595-1605.

[9] 魏寶君，LIU Q H.層狀介質(zhì)中計算體積分方程的弱化BCGS-FFT算法[J].中國石油大學(xué)學(xué)報：自然科學(xué)版，2007，31(1)：49-55.WEI Bao-jun，LIU Q H.Weak-form BCGS-FFT algorithm for volume integral equations in stratified medium[J].Journal of China University of Petroleum(Edition of Natural Science)，2007，31(1)：49-55.

[10] 魏寶君，王穎，王甜甜.電磁波電阻率儀器的基本理論及其在隨鉆測量中的應(yīng)用[J].地球物理學(xué)進(jìn)展，2009，24(2)：774-781.WEI Bao-jun，WANG Ying，WANG Tian-tian.The basic theory of electromagnetic wave resistivity instrument and its application to measurement of while-drilling[J].Progress in Geophysics，2009，24(2)：774-781.