范宜仁 李 虎 胡云云 孫慶濤

（1.中國(guó)石油大學(xué)地球科學(xué)與技術(shù)學(xué)院，山東 青島 266580；2.中國(guó)石油大學(xué)（華東）CNPC測(cè)井重點(diǎn)實(shí)驗(yàn)室，山東 青島 266580）

引 言

傾斜各向異性地層隨鉆電磁波響應(yīng)數(shù)值模擬與特征分析，對(duì)基于隨鉆電測(cè)井資料進(jìn)行實(shí)時(shí)地質(zhì)導(dǎo)向與儲(chǔ)層評(píng)價(jià)具有非常重要的意義［1-2］.目前，常用的隨鉆電磁波數(shù)值模擬方法包括:有限元法（Finite Element Method，F(xiàn)EM）［3-4］、有 限 差 分 法 （Finite Difference Method，F(xiàn)DM）［5-6］、傳輸線法（Transmission Line Matrix Method，TLM）［7］、數(shù)值模式匹配法（Numerical Mode Matching Method，NMM）［8］和快速傅里葉漢克爾變換結(jié)合的共軛梯度迭代法（Conjugate Gradient Method with Fast Fourier and Fast Hankel Transforms，CGFFHT）［9］.其中 FEM與FDM可以精細(xì)模擬隨鉆電磁波在任意傾角各向異性地層響應(yīng)，但頻域有限差分與有限元法都需要求解大型矩陣，計(jì)算過程復(fù)雜耗時(shí)；NMM 與CGFFHT求解相對(duì)較快，但無法實(shí)現(xiàn)對(duì)復(fù)雜儀器結(jié)構(gòu)與地層模型的精細(xì)模擬.時(shí)域有限差分（Finite-Difference Time-Domain，F(xiàn)DTD）方法最早由 K.S.Yee1966年提出，該方法采用時(shí)間迭代的方法避免了大型稀疏矩陣的求解，能更好地適用于三維復(fù)雜介質(zhì)隨鉆電磁波響應(yīng)的模擬.本文采用柱坐標(biāo)系FDTD實(shí)現(xiàn)了傾斜各向異性地層隨鉆電磁波響應(yīng)模擬，并通過對(duì)不同地層條件下響應(yīng)特征分析，理論研究了地層傾角、電導(dǎo)率各向異性、儀器發(fā)射頻率等對(duì)隨鉆電磁波響應(yīng)的影響.

1 方法原理

1.1 圓柱坐標(biāo)系下各向異性電導(dǎo)率獲取

隨鉆地層模型與儀器結(jié)構(gòu)示意圖如圖1，假設(shè)地層層理方向與儀器軸線方向夾角為θ，則可以通過以下變換，將局部直角坐標(biāo)系x′y′z′中電導(dǎo)率張量變換到與儀器軸重合的直角坐標(biāo)系xyz下，有

第二次旋轉(zhuǎn)，將直角坐標(biāo)系中電導(dǎo)率張量變換到圓柱坐標(biāo)系ρφz下，本次變換可以表示為

圖1 隨鉆地層模型與儀器結(jié)構(gòu)示意圖

此時(shí)σ即為圓柱坐標(biāo)系下電導(dǎo)率張量.

1.2 圓柱坐標(biāo)系下各向異性地層FDTD差分格式

電性各向異性介質(zhì)本構(gòu)方程可表示為

介質(zhì)的磁性仍設(shè)為各向同性，即:

此時(shí)，Maxwell方程可表示為

式（7）時(shí)域差分格式與各向同性時(shí)相同，即，

對(duì)式（8）進(jìn)行離散得

整理后得

1.3 介電常數(shù)放大策略

在磁導(dǎo)率為μ，介電常數(shù)為ε，電導(dǎo)率為σ的地層中，F(xiàn)DTD最大時(shí)間步長(zhǎng)Δtc由空間步長(zhǎng)決定，

式中v為電磁波傳播速度，

在隨鉆電磁波測(cè)井?dāng)?shù)值模擬過程中，由于位移電流遠(yuǎn)小于傳導(dǎo)電流，可以通過放大介電常數(shù)增大時(shí)間步長(zhǎng).將式（8）重新整理為

當(dāng)σh（σv）＞＞ωε時(shí)，方程等效為

在模擬100kHz發(fā)射頻率時(shí)，通過放大介電常數(shù)，計(jì)算需要時(shí)間可縮短80%，不同介電常數(shù)放大倍數(shù)下相位差、幅度比計(jì)算結(jié)果見圖2，表1中給出了模型計(jì)算條件與CPU時(shí)間對(duì)比.

表1 不同介電常數(shù)放大倍數(shù)模型計(jì)算時(shí)間表

圖2 介電常數(shù)放大策略效果圖

圖3 模擬結(jié)果與文獻(xiàn)［10］中Hwa計(jì)算結(jié)果對(duì)比

1.4 模擬結(jié)果驗(yàn)證

圖3中給出了不同地層傾角條件下線圈距為［2430］in隨鉆電磁波相位差計(jì)算結(jié)果，并將計(jì)算結(jié)果與文獻(xiàn)［10］中的結(jié)果進(jìn)行對(duì)比，兩者結(jié)果吻合很好，驗(yàn)證了算法的正確性.

2 各向異性地層隨鉆電磁波響應(yīng)特征

實(shí)際測(cè)井過程中，通過對(duì)兩接收線圈的相位差與幅度比進(jìn)行刻度轉(zhuǎn)化為相位／幅度電阻率，相位差（DP）與幅度比（RA）定義為

式中:V為接收線圈電位；arg為取相位角.對(duì)各向異性介質(zhì)，定義各向異性系數(shù)

2.1 各向異性系數(shù)的影響

為了研究各向異性系數(shù)對(duì)相位（幅度）電阻率的影響，模擬了不同各向異性系數(shù)條件下傾斜地層隨鉆電磁波測(cè)井響應(yīng)，模型中水平電阻率分別選取為1Ω·m、10Ω·m，模擬結(jié)果見圖4.

地層傾角較小時(shí)（小于30°），不同水平電阻率條件下，各向異性系數(shù)對(duì)視電阻率影響較小，隨地層傾角增大，各向異性系數(shù)對(duì)視電阻率影響增大，且地層水平電阻率越低，響應(yīng)受地層各向異性影響越大，當(dāng)?shù)貙觾A角較大時(shí)，隨各向異性系數(shù)增大，視電阻率甚至?xí)^垂直電阻率.

圖4 不同各向異性系數(shù)條件下隨鉆電磁波測(cè)井響應(yīng)

2.2 發(fā)射頻率的影響

傳統(tǒng)隨鉆電磁波儀器發(fā)射頻率一般為500kHz與2MHz，新型的隨鉆方位電磁波增加了100kHz頻率，為了研究不同發(fā)射頻率對(duì)各向異性系數(shù)的敏感性，模擬了地層各向異性系數(shù)為，水平電阻率為0.5Ω·m時(shí)不同地層傾角條件下隨鉆電磁波響應(yīng)（如圖5），隨發(fā)射頻率增大，視電阻率受各向異性影響增強(qiáng)，當(dāng)?shù)貙觾A角較大時(shí)，2MHz頻率下視電阻率甚至?xí)h(yuǎn)遠(yuǎn)超過垂直電阻率.

圖5 不同發(fā)射頻率條件下隨鉆電磁波測(cè)井響應(yīng)

3 結(jié) 論

基于柱坐標(biāo)系時(shí)域有限差分實(shí)現(xiàn)了傾斜各向異性地層隨鉆電磁波測(cè)井響應(yīng)數(shù)值模擬，通過與已有結(jié)果對(duì)比，驗(yàn)證了算法的正確性.通過放大介電常數(shù)，在模擬100kHz發(fā)射頻率時(shí)，計(jì)算時(shí)間可縮短80%.

地層傾角小于30°時(shí)，隨鉆電磁波視電阻率主要反映地層水平電阻率，地層傾角越大、發(fā)射頻率越高，地層電阻率各向異性對(duì)視電阻率影響越大，且相位電阻率比幅度電阻率更加敏感，當(dāng)發(fā)射頻率較高、地層傾角較大時(shí)，隨鉆電磁波視電阻率甚至?xí)h(yuǎn)遠(yuǎn)超過地層垂直電阻率.

［1］ADERSON B，BONNET S，ROSTHAI R，et al.Response of 2MHz LWD resistivity and wireline induction tools in dipping beds and laminated［J］.J Log Analyst，1992，33（5）:461-475.

［2］BITTAR M S，HU G Y.The effects of rock anisotropy on LWD toroidal resistivity sensors［C／OL］／／SPWLA 46th Annual Logging Symposium，Noordwijk，The Netherlands，June 6-9，2004［2012-10-09］.http:／／www.spwla.org／publications／view／item／2785.

［3］CHANG S K，ANDERSON B.Simulation of induction logging by the finite-element method［J］.Geophysics，1984，49（11）:1943-1958.

［4］孫向陽(yáng)，聶在平，李愛勇，等.用高階疊層矢量有限元法計(jì)算隨鉆測(cè)井的三維電磁響應(yīng)［J］.電波科學(xué)學(xué)報(bào)，2009，24（2）:273-279.SUN Xiangyang，NIE Zaiping，LI Aiyong，et al.The modeling of logging-while-drilling tool's threedimensional electromagnetic response using the high order hierarchical vector finite element method［J］.Chinese Journal of Radio Science，2009，24（2）:273-279.（in Chinese）

［5］HOLLAND R.Finite-difference time-domain（FDTD）analysis of magnetic diffusion［J］.IEEE Trans on Electromagnetic Compatibility，1994，36（1）:32-39.

［6］YEE K S.Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media［J］.IEEE Trans on Antennas and Propagation，1966，14（3）:302-307.

［7］LI J，CE L.A three-dimensional transmission line matrix（TLM）method for simulations of logging tools［J］.IEEE Transaction on Geoscience and Remote Sensing，2000，38（4）:1522-1529.

［8］LIU Q H，CHEW W C，TAHERIAN M R，et al.A modeling study of electromagnetic propagation tool in complicated borehole environments［J］.Log Analyst，1989，30（6）:424-436.

［9］ZHANG Z Q，LIU Q H.Simulation of induction logging response using conjugate gradient method with non-uniform fast Fourier and fast Hankel transforms［J］.Radio Science，2001，36（4）:59-608.

［10］LEE H O，TEIXEIRA F L.Cylindrical FDTD analysis of LWD tools through anisotropic dipping layered earth media［J］.IEEE Transaction Geoscience and Remote Sensing，2007，45（2）:383-388.