王亞洲，汪宏年*,莫修文，康莊莊，殷長(zhǎng)春

1 吉林大學(xué)物理學(xué)院計(jì)算方法與軟件國(guó)際中心,長(zhǎng)春 130012 2 吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院,長(zhǎng)春 130026

0 引言

目前,隨鉆電磁波測(cè)井已經(jīng)成為一種重要測(cè)井方法（Bell et al.,2006;Seydoux et al.,2014）,是進(jìn)行地質(zhì)導(dǎo)向、油水層識(shí)別和含油飽和度計(jì)算等工作的重要工具.常規(guī)隨鉆電磁波測(cè)井中采用的線(xiàn)圈系收發(fā)距往往較短（最大收發(fā)距3 m）,且工作頻率較高（100 k、400 k和2 MHz）,因此其探測(cè)范圍往往較小（最大探測(cè)深度5 m左右）（Li et al.,2005;Fang,2011;于蕾等,2021）.隨著大斜度井和水平井在油氣勘探開(kāi)發(fā)中的廣泛應(yīng)用,不僅需要地質(zhì)導(dǎo)向技術(shù)實(shí)現(xiàn)井眼軌跡優(yōu)化,還需要精細(xì)了解井周數(shù)十米范圍內(nèi)的油氣層空間分布情況,以便進(jìn)行深部鉆井的智能“巡航”,在更大范圍內(nèi)優(yōu)化井眼軌跡,提高深層油氣開(kāi)發(fā)效率以降低成本（Kennedy et al.,2009;Omeragic et al.,2005）.因此,近十多年來(lái),在常規(guī)隨鉆電磁波測(cè)井研究應(yīng)用基礎(chǔ)上,新型隨鉆超深前視電阻率測(cè)井技術(shù)也得到了快速發(fā)展（Rabinovich et al.,2011）.該項(xiàng)技術(shù)采用多頻（1、2、4、8 kHz和16 kHz等）、長(zhǎng)收發(fā)距（5～35 m）多分量線(xiàn)圈系組合測(cè)量（Dong et al.,2015;Wu et al.,2018）,并結(jié)合現(xiàn)代三維電磁正反演技術(shù)對(duì)井眼周?chē)?5 m徑向范圍內(nèi)電導(dǎo)率空間分布進(jìn)行精細(xì)成像,實(shí)現(xiàn)鉆井優(yōu)化、地質(zhì)導(dǎo)向以及儲(chǔ)層地質(zhì)構(gòu)造和水文地質(zhì)條件的高分辨率精細(xì)描述.

在隨鉆超深前視電阻率測(cè)井技術(shù)研究、開(kāi)發(fā)和應(yīng)用過(guò)程中,需要進(jìn)行儀器參數(shù)優(yōu)化設(shè)計(jì)、考察復(fù)雜地層條件下探測(cè)特性并建立一套相應(yīng)的資料處理和反演方法,特別需要指出的是,當(dāng)前三維反演效率是滿(mǎn)足不了實(shí)時(shí)處理需要的,所以在隨鉆資料處理中,可以根據(jù)地下電導(dǎo)率空間分布情況,預(yù)先設(shè)計(jì)井眼軌跡,并根據(jù)井眼軌跡和地層電導(dǎo)率分布,事先計(jì)算出一套正演模擬數(shù)據(jù)（理論合成響應(yīng)）,然后將實(shí)測(cè)資料與理論合成響應(yīng)進(jìn)行對(duì)比,達(dá)到地質(zhì)導(dǎo)向的目的,這些工作都需要隨鉆超深前視電阻率測(cè)井三維正演模擬技術(shù).從理論上說(shuō),根據(jù)地層模型以及計(jì)算精度要求可以選擇不同的數(shù)值模擬方法進(jìn)行隨鉆超深前視電阻率測(cè)井響應(yīng)的數(shù)值模擬.例如,對(duì)于分層介質(zhì),可以選用傳輸線(xiàn)法（Michalski and Mosig,1997;Yang et al.,2014）、傳播矩陣法（康莊莊等,2020）、廣義反射系數(shù)方法（Hong et al.,2014,2017）以及數(shù)值模式匹配法（Wang et al.,2008;林藺等,2017;汪宏年等,2021）等算法,而對(duì)于復(fù)雜的三維地層模型,則可以采用三維有限元法（齊彥福等,2020）、三維有限差分法（Wang and Signorelli,2004;楊守文等,2009）、三維有限體積法（Haber et al.,2000;Wang et al.,2020;陳博等,2021）和積分方程法（陳桂波等,2009;湯井田等,2018）等.其中,有限體積法具有較好的守恒特性且易于實(shí)現(xiàn),而引入耦合勢(shì)能夠有效的解決低感應(yīng)數(shù)問(wèn)題（Weiss and Newman,2003）.此外,在對(duì)源進(jìn)行離散時(shí),為了獲得更高的精度,常常將電磁場(chǎng)分解為背景場(chǎng)與二次場(chǎng),從而避免了直接對(duì)源離散（Irons et al.,2012;彭榮華等,2018）,但增加了額外的背景場(chǎng)計(jì)算成本,限制了其在三維數(shù)值方法中的應(yīng)用,目前針對(duì)如何快速計(jì)算背景場(chǎng)的討論較少.

隨鉆超深前視電阻率測(cè)井儀器的收發(fā)距變化范圍大,本文將基于二次場(chǎng)耦合勢(shì)三維有限體積法研究非均質(zhì)各向異性地層中任意傾角的斜井中隨鉆超深前視電阻率測(cè)井響應(yīng)數(shù)值模擬算法,以便提高正演模擬效率和整個(gè)計(jì)算區(qū)域內(nèi)近場(chǎng)和遠(yuǎn)場(chǎng)的計(jì)算精度.為此,本文將通過(guò)Yee氏非均勻交錯(cuò)網(wǎng)格與三維有限體積法對(duì)二次場(chǎng)耦合勢(shì)Helmholtz方程進(jìn)行離散,并根據(jù)傳輸線(xiàn)算法形成層狀地層中電磁場(chǎng)并矢Green函數(shù)庫(kù),同時(shí)結(jié)合二維Newton插值快速計(jì)算異常體內(nèi)任意位置上的背景電場(chǎng)和散射電流,有效提高體分布散射電流源離散效率.在此基礎(chǔ)上,采用并行直接求解器PARDISO快速穩(wěn)定求解多源正演問(wèn)題.最后,將算法應(yīng)用于隨鉆超深前視電阻率測(cè)井儀器的正演模擬,驗(yàn)證了算法的有效性,并討論了建庫(kù)插值方法的計(jì)算效率.

1 基本原理

本節(jié)首先給出儀器結(jié)構(gòu)和地層模型,并通過(guò)地層模型與電磁場(chǎng)的分解引入二次場(chǎng)與二次場(chǎng)耦合勢(shì)Helmholtz方程,然后重點(diǎn)討論體分布電流源計(jì)算與離散,給出背景場(chǎng)快速計(jì)算方法,以提高二次場(chǎng)方法的計(jì)算效率.

1.1 地層模型與電磁場(chǎng)分解

中國(guó)科學(xué)院地質(zhì)與地球物理研究所目前正在自主研發(fā)隨鉆超深前視電阻率儀器,圖1是儀器結(jié)構(gòu)示意圖,該儀器主要由一個(gè)三分量正交共位發(fā)射線(xiàn)圈Txyz（圖1b）和三個(gè)不同收發(fā)距三分量正交共位接收線(xiàn)圈Rxyz（圖1a）組成,并且發(fā)射器靠近鉆頭.為方便起見(jiàn)本文假定三個(gè)不同收發(fā)距L分別為5 m、15和35 m,而五個(gè)工作頻率則假定為1、2、4、8 kHz和16 kHz （Dong et al.,2015;Wu et al.,2018）.

圖2是傾斜井眼與三維非均質(zhì)各向異性地層模型分解示意圖.其中,井眼傾角為θ,而地層模型由起伏地層界面和幾個(gè)局部異常體組成（圖2a）,各層水平和垂直電導(dǎo)率分別用σH,n和σV,n,n=1,2,3表示.為便于描述傾斜或水平井眼情況下的正演過(guò)程,引入地層坐標(biāo)系x′y′z′和井眼坐標(biāo)系xyz（圖2a）,其中,地層坐標(biāo)系的z′軸垂直于層狀背景地層界面d1和d2,且x′y′平面與背景地層界面平行.井眼坐標(biāo)系的z軸與儀器鉆進(jìn)方向一致,y軸與地層坐標(biāo)系中y′軸平行,井眼坐標(biāo)系可由地層坐標(biāo)系繞y′軸逆時(shí)針旋轉(zhuǎn)θ角得到.

參照?qǐng)D2a中地層電導(dǎo)率空間分布,可以將其分解為層狀背景模型（圖2b）與三維局部異常體的組合（圖2c）,對(duì)電磁場(chǎng)進(jìn)行相應(yīng)的分解,可以得到總電磁場(chǎng)方程（見(jiàn)附錄A中（A1））、背景電磁場(chǎng)方程（A3）以及二次電磁場(chǎng)方程（A4）.其中背景電磁場(chǎng)采用傳輸線(xiàn)方法求解,二次電磁場(chǎng)則需要采用數(shù)值方法求解,并且為了克服低感應(yīng)數(shù)問(wèn)題引入耦合勢(shì)（A5）,整理后得到二次場(chǎng)耦合勢(shì)Helmholtz方程:

圖1 隨鉆超深前視電阻率儀器結(jié)構(gòu)示意圖（a） 三分量陣列接收器;（b） 三分量發(fā)射器.Fig.1 Schematic diagram of the instrument structure of LWD ultra-deep look ahead multicomponent resistivity logging（a） Three-component array receivers;（b） Three-component transmitters.

圖2 復(fù)雜地層模型與分解示意圖（a） 三維各向異性地層模型;（b） 水平層狀背景地層模型;（c） 三維局部異常體模型.Fig.2 Schematic diagram of complex formation model and decomposition（a） 3D anisotropic formation model;（b） Horizontal layered background model;（c） 3D local abnormal bodies model.

=-iωμ0ΔJS,p（r,rs）

（1）

和

（2）

為求解方程（1）和（2）,選擇足夠大計(jì)算區(qū)域Ω并在其外邊界?Ω采用理想電導(dǎo)體邊界條件（PEC）:

（3）

1.2 三維有限體積法

（4）

（5）

借助有限體積法對(duì)方程（5）分別在剖分單元Vi+1/2,j,k、Vi,j+1/2,k和Vi,j,k+1/2上進(jìn)行體積平均,從而得到方程（1）的右端項(xiàng)離散結(jié)果:

（6a）

（6b）

（6c）

（7）

（8）

（9）

將方程（9）的二次場(chǎng)數(shù)值結(jié)果與背景場(chǎng)結(jié)合,可以計(jì)算出各個(gè)發(fā)射源產(chǎn)生的電磁場(chǎng):

（10）

1.3 背景電磁場(chǎng)快速算法

（12）

其中,基函數(shù)為

（13）

將式（12）的計(jì)算結(jié)果代入到附錄（B6）式中就可以得到地層坐標(biāo)系中背景電場(chǎng)并矢Green函數(shù),結(jié)合如下轉(zhuǎn)化公式:

2 數(shù)值結(jié)果

本節(jié)首先利用2.5D數(shù)值模式匹配算法（2.5D NMM）（Wang et al.,2008;林藺等,2017）的數(shù)值結(jié)果對(duì)本文采用的二次場(chǎng)三維有限體積法（3D FV）進(jìn)行對(duì)比驗(yàn)證,然后討論背景電場(chǎng)計(jì)算過(guò)程,并將算法應(yīng)用于隨鉆超深前視電阻率測(cè)井儀器的正演模擬.在數(shù)值結(jié)果中,將以發(fā)射源與接收器中心點(diǎn)作為記錄點(diǎn)的位置,以便于不同收發(fā)距正演結(jié)果的對(duì)比,并將正演模擬結(jié)果轉(zhuǎn)化為視電導(dǎo)率,以降低頻率和收發(fā)距對(duì)正演結(jié)果的影響.視電導(dǎo)率各個(gè)分量為（張燁等,2012）:

（15）

此外,在如下二次場(chǎng)耦合勢(shì)三維有限體積法數(shù)值模擬中,網(wǎng)格節(jié)點(diǎn)數(shù)為（Nx,Ny,Nz）=（52,52,122）,并且x和y方向中間取30個(gè)步長(zhǎng)為1 m的等間距網(wǎng)格,在中間網(wǎng)格外面再取10個(gè)步長(zhǎng)為2 m等間距網(wǎng)格,最外圍是12個(gè)漸變網(wǎng)格;而z方向的中間是40個(gè)步長(zhǎng)為1 m等間距網(wǎng)格,在其外部則是40個(gè)步長(zhǎng)為2 m等間距網(wǎng)格,然后是30個(gè)步長(zhǎng)為3 m的等間距網(wǎng)格,最外圍則是12個(gè)漸變網(wǎng)格,漸變網(wǎng)格根據(jù)Lebedev網(wǎng)格（Davydycheva et al.,2003）進(jìn)行選取.文中所有數(shù)值結(jié)果都是在CPU型號(hào)為Intel（R） Xeon（R） Platinum 8269CY、主頻2.50 GHz的工作站上計(jì)算得到的.

2.1 算法檢驗(yàn)

為便于用2.5D NMM對(duì)二次場(chǎng)耦合勢(shì)有限體積法加以檢驗(yàn),假定在水平和垂直電導(dǎo)率分別為0.2和0.1 S·m-1的均勻各向異性地層中包含著一個(gè)垂直柱狀異常體（見(jiàn)圖3）,該柱狀異常體的水平和垂直電導(dǎo)率分別為0.125和0.05 S·m-1、其半徑和高度分別為10和20 m,而井眼軌跡正好穿過(guò)圓柱的軸線(xiàn).此外,還給出總場(chǎng)耦合勢(shì)有限體積法的計(jì)算結(jié)果,以對(duì)比二次場(chǎng)方法與總場(chǎng)法的精度.

圖4是儀器工作頻率為1 kHz,收發(fā)距分別是5、15和35 m時(shí)三種不同正演算法計(jì)算得到的三個(gè)主分量的結(jié)果對(duì)比,其中,圖4（a1）、（a2）和（a3）,圖4（b1）、（b2）和（b3）以及圖4（c1）、（c2）和（c3）分別對(duì)應(yīng)收發(fā)距為5、15和35 m時(shí)的主分量視電導(dǎo)率σa,xx、σa,yy和σa,zz,黑色實(shí)線(xiàn)是NMM得到的結(jié)果,綠色星點(diǎn)是二次場(chǎng)耦合勢(shì)有限體積法得到的結(jié)果,藍(lán)色虛線(xiàn)是總場(chǎng)耦合勢(shì)有限體積法得到的結(jié)果.二次場(chǎng)耦合勢(shì)有限體積法與NMM得到的結(jié)果吻合很好,其相對(duì)誤差（絕對(duì)值）均小于3%（見(jiàn)圖4（a1′）、（a2′）和（a3′）,圖4（b1′）、（b2′）和（b3′）以及圖4（c1′）、（c2′）和（c3′）中綠色星點(diǎn)）,而總場(chǎng)耦合勢(shì)有限體積法的誤差明顯更大（見(jiàn)圖4（a1′）、（a2′）和（a3′）,圖4（b1′）、（b2′）和（b3′）以及圖4（c1′）、（c2′）和（c3′）中藍(lán)色虛線(xiàn)）.此外,對(duì)比二次場(chǎng)法和總場(chǎng)法的相對(duì)誤差曲線(xiàn)可以發(fā)現(xiàn),二次場(chǎng)法的相對(duì)誤差雖然都較小,但在異常體所在位置相對(duì)誤差值會(huì)增大,產(chǎn)生這一現(xiàn)象的主要原因是這些深度點(diǎn)對(duì)應(yīng)的發(fā)射源均位于柱狀異常體內(nèi)部,導(dǎo)致發(fā)射源所在單元中背景電場(chǎng)插值誤差增大,從而影響了計(jì)算精度.而總場(chǎng)耦合勢(shì)有限體積法相對(duì)誤差卻與二次場(chǎng)法相對(duì)誤差變化特征完全不同,在異常體所在位置,總場(chǎng)法的相對(duì)誤差明顯降低,而在異常體外部的相對(duì)誤差卻明顯增大,產(chǎn)生這一現(xiàn)象的主要原因是總場(chǎng)法中磁偶極子源的離散精度受網(wǎng)格大小影響非常明顯,在異常體所在區(qū)域,因?yàn)榫W(wǎng)格尺寸較?。? m）導(dǎo)致其離散誤差下降,而在異常體外部采用了較大的網(wǎng)格尺寸（2 m）,引起離散誤差明顯增大.因此可以得出結(jié)論,在相同網(wǎng)格尺寸下,總場(chǎng)法的計(jì)算精度明顯低于二次場(chǎng)方法,隨鉆超深前視電阻率測(cè)井的數(shù)值模擬需要考慮較大的計(jì)算區(qū)域,采用二次場(chǎng)方法有利于放寬對(duì)網(wǎng)格大小的要求.

圖3 圓柱異常體模型示意圖Fig.3 Schematic diagram of cylinder model

圖4 圓柱模型中3D FV與NMM計(jì)算結(jié)果對(duì)比與誤差分析Fig.4 The comparison of results obtained by 3D FV and NMM and error analysis in a cylinder model

2.2 長(zhǎng)方體異常體模型

圖5 長(zhǎng)方體異常體模型示意圖Fig.5 Schematic diagram of cuboid model

圖5是雙層各向異性層狀背景地層中含有單個(gè)長(zhǎng)方體異常體模型示意圖.雙層背景模型中,頂層水平和垂直電導(dǎo)率分別為0.5和0.2 S·m-1,底層水平和垂直電導(dǎo)率分別為0.1和0.05 S·m-1,異常體位于上部地層且水平和垂直電導(dǎo)率分別為0.02和0.005 S·m-1.在地層坐標(biāo)系o′x′y′z′中,層界面位于z′=0 m處,異常體中心位置坐標(biāo)為（0,0,6）m,在x′、y′和z′三個(gè)方向上的長(zhǎng)度分別為30、30和6 m,且上下邊界與背景地層界面平行.此外,假定傾斜井眼位于x′z′平面內(nèi)且傾角為60°,井眼軌跡經(jīng)過(guò)地層坐標(biāo)系中的點(diǎn)（0,0,-20）m,并選擇該點(diǎn)作為井眼坐標(biāo)系oxyz的原點(diǎn).整個(gè)井眼中共包含81個(gè)等間距測(cè)點(diǎn),每一個(gè)測(cè)點(diǎn)都需要對(duì)三個(gè)不同位置發(fā)射源產(chǎn)生的背景電磁場(chǎng)進(jìn)行計(jì)算,但不同測(cè)點(diǎn)對(duì)應(yīng)的發(fā)射源位置可能會(huì)重復(fù),所以在正演過(guò)程中實(shí)際只需要計(jì)算167個(gè)位于不同位置的三分量發(fā)射源產(chǎn)生的背景電磁場(chǎng).

表1 計(jì)算背景電磁場(chǎng)CPU耗時(shí)對(duì)比Table 1 Comparison of CPU time for background electric field by different algorithms

為了對(duì)插值法計(jì)算精度進(jìn)行考察,分別用插值法與傳統(tǒng)TLM算法計(jì)算了傾斜井眼軌跡上背景視電導(dǎo)率（圖6）,選定的收發(fā)距和工作頻率分別為35 m和1 kHz.圖6（a）、（b）、（c）、（d）和（e）分別對(duì)應(yīng)五個(gè)背景視電導(dǎo)率分量σa,xx、σa,yy、σa,zz、σa,zx和σa,xz,其中黑色實(shí)線(xiàn)是TLM算法得到的結(jié)果,而綠色圓點(diǎn)是插值法得到的結(jié)果,兩種方法得到的曲線(xiàn)吻合的非常好;圖6f給出了對(duì)應(yīng)的相對(duì)誤差,各個(gè)分量的相對(duì)誤差（絕對(duì)值）均不超過(guò)1%.

圖7是模型5中三個(gè)不同收發(fā)距（L=5、15 m和35 m）和五個(gè)不同工作頻率（1、2、4、8和16 kHz）情況下總視電導(dǎo)率和背景視電導(dǎo)率正演模擬結(jié)果,其中,圖7（a1—a5）、圖7（b1—b5）和圖7（c1—c5）分別是收發(fā)距為5 m、15 m和35 m時(shí)五個(gè)不同頻率下五個(gè)視電導(dǎo)率分量σa,xx、σa,yy、σa,zz、σa,zx和σa,xz的曲線(xiàn),圖中的實(shí)線(xiàn)表示總視電導(dǎo)率,虛線(xiàn)為背景視電導(dǎo)率.從圖可以明顯看出,收發(fā)距越長(zhǎng)異常體對(duì)視電導(dǎo)率曲線(xiàn)的影響范圍越大,說(shuō)明收發(fā)距增加能夠明顯提高儀器的探測(cè)范圍;此外,也可以看到工作頻率變化對(duì)視電導(dǎo)率的影響,在源距不變的情況下，頻率越低異常體對(duì)視電導(dǎo)率曲線(xiàn)的影響也越大,說(shuō)明頻率改變也能夠有效調(diào)節(jié)儀器的探測(cè)范圍.因此,通過(guò)多頻率多源距的組合測(cè)量方式,能夠獲得儀器周?chē)煌綔y(cè)距離內(nèi)的地層電導(dǎo)率分布信息,有利于進(jìn)行井眼周?chē)S電導(dǎo)率反演成像.

圖6 插值法與傳輸線(xiàn)法得到的背景磁場(chǎng)響應(yīng)對(duì)比Fig.6 The comparison of results obtained by interpolation method and conventional TLM

圖7 長(zhǎng)方體異常體模型中不同收發(fā)距與頻率下儀器響應(yīng)Fig.7 The instrument responses at different source distances and frequencies in cuboid model

2.3 二維起伏界面模型

圖8 二維起伏界面模型示意圖Fig.8 Schematic diagram of 2D undulating interface model

下面將考察二維起伏界面模型（圖8）中的測(cè)井響應(yīng),模型中背景介質(zhì)為雙層模型,頂層的水平和垂直電導(dǎo)率分別為0.1和0.05 S·m-1,底層的水平和垂直電導(dǎo)率分別為0.01和0.005 S·m-1.在地層坐標(biāo)系x′y′z′中,層界面位于z′=100 m處,起伏界面中心位置位于x′=0 m處,寬D取40 m,高H分別取4、8和12 m,以便考察起伏界面高度對(duì)儀器響應(yīng)的影響,起伏界面在x′z′垂直平面上的截線(xiàn)方程為z′=H（x′+D/2）（x′-D/2）/（D/2）2+100.儀器位于底層地層中,井眼傾角為89°,并且井眼坐標(biāo)系以地層坐標(biāo)系中的點(diǎn)（0,0,105）m為原點(diǎn).

首先給出起伏界面周?chē)艌?chǎng)的空間分布,為簡(jiǎn)潔起見(jiàn),這里僅給出H為8 m時(shí)的結(jié)果.圖9是單位發(fā)射源Mz在位于井眼坐標(biāo)系的坐標(biāo)原點(diǎn)時(shí)產(chǎn)生的磁場(chǎng)（虛部）在y′=0垂直截面上的分布情況.其中,圖9（a1—a5）、圖9（b1—b5）和圖9（c1—c5）分別為背景磁場(chǎng)、二次磁場(chǎng)和總磁場(chǎng)在五個(gè)不同頻率下的計(jì)算結(jié)果,結(jié)果顯示,隨著頻率增加,發(fā)射源附近的背景磁場(chǎng)、二次磁場(chǎng)和總磁場(chǎng)的場(chǎng)值會(huì)逐漸變大,不同頻率對(duì)應(yīng)的磁場(chǎng)空間分布形態(tài)相似.在起伏界面周?chē)?背景磁場(chǎng)方向沿逆時(shí)針?lè)较颍▓D9（a1—a5））,而二次磁場(chǎng)方向沿順時(shí)針?lè)较颍▓D9（b1—b5））,因此起伏界面的存在使得總磁場(chǎng)的場(chǎng)值變小,從總磁場(chǎng)中可以明顯看出起伏界面對(duì)磁場(chǎng)分布形態(tài)的影響（圖9（c1—c5））.

圖9 Mz發(fā)射源激發(fā)的磁場(chǎng)在y′=0垂直平面的空間分布圖Fig.9 Spatial distribution of magnetic field due to Mz source in y′=0 plane

圖10 89° 斜井中雙層背景介質(zhì)的儀器響應(yīng)Fig.10 The instrument responses for a double-layer background medium in 89° deviated well

圖11 89° 斜井中不同二維起伏界面高度下的儀器響應(yīng)Fig.11 The instrument responses for different heights of 2D undulating interface in 89° deviated well

此外,在考察二維起伏界面對(duì)儀器響應(yīng)的影響之前,先討論圖8模型中雙層背景介質(zhì)的測(cè)井響應(yīng),以便于了解在89°的大斜度井中不同源距的響應(yīng)特征以及探邊能力.為此,儀器工作頻率選定為1 kHz,圖10（a）、（b）、（c）、（d）和（e）給出了三個(gè)不同源距五個(gè)背景視電導(dǎo)率分量σa,xx、σa,yy、σa,zz、σa,zx和σa,xz的正演結(jié)果,整個(gè)背景響應(yīng)沿水平方向x′的測(cè)量范圍為-200～1500 m,即記錄點(diǎn)從距離層邊界1.5 m逐步增加到31.2 m.先觀(guān)察并比較圖10（a）、（b）和（c）中三個(gè)不同源距的主分量背景測(cè)井曲線(xiàn),可以看出,短源距L=5 m三個(gè)主分量測(cè)井曲線(xiàn)在整個(gè)的計(jì)算范圍內(nèi)均是隨變量x′單調(diào)變化的（也是層邊界距離的單調(diào)函數(shù)）,這將十分有利于確定層邊界位置,但當(dāng)變量x′增加到大約600 m以后視電導(dǎo)率的變化變緩,說(shuō)明其探邊能力較弱;源距L=15 m兩個(gè)主分量σa,xx和σa,yy與L=35 m主分量σa,yy在整個(gè)計(jì)算范圍內(nèi)也是單調(diào)變化的,并且視電導(dǎo)率在更遠(yuǎn)處才趨于平緩,說(shuō)明其探邊能力明顯增強(qiáng).此外,不難看出,L=15 m的主分量σa,zz和L=35 m的主分量σa,xx和σa,zz曲線(xiàn)不是單調(diào)變化,即在不同位置可能有相同的測(cè)量結(jié)果,顯然不利于層界面反演.不過(guò),與主分量不同的是,三個(gè)不同源距的交叉分量σa,zx和σa,xz均是隨變量x′單調(diào)變化的,但其衰減較快,探測(cè)深度比主分量小.不同源距不同分量測(cè)井曲線(xiàn)的不同特征對(duì)于層界面反演是十分有利的,可以用不同曲線(xiàn)分別進(jìn)行層界面反演,根據(jù)反演結(jié)果是否一致能夠判斷出層界是否起伏或?qū)咏缑娓浇欠裼挟惓ｓw,因此采用多源距多分量組合測(cè)量是十分必要的,能夠起到相互補(bǔ)充信息的目的.

在給出了層狀背景響應(yīng)的基礎(chǔ)上,將進(jìn)一步考察層界面起伏情況下的儀器響應(yīng),由于起伏地形僅有40 m的寬度,所以?xún)H給出x′變化范圍為-80～80 m（記錄點(diǎn)與水平層界面距離為3.6～6.4 m）內(nèi)的正演結(jié)果,在該范圍外測(cè)井響應(yīng)受起伏地形的影響已非常小,并且為簡(jiǎn)潔起見(jiàn),這里僅給出工作頻率為1 kHz的結(jié)果.圖11（a1—a5）、圖11 （b1—b5）和圖10 （c1—c5）分別是源距為L(zhǎng)=5、15和35 m時(shí)五個(gè)視電導(dǎo)率分量σa,xx、σa,yy、σa,zz、σa,zx和σa,xz的曲線(xiàn),其中黑、紅和藍(lán)三條實(shí)線(xiàn)分別對(duì)應(yīng)4、8和12 m三個(gè)不同高度的起伏界面,綠色實(shí)線(xiàn)則是層狀背景介質(zhì)中的背景視電導(dǎo)率曲線(xiàn).起伏界面的存在直接增加了儀器與層界面的距離,起伏界面高度越大,各個(gè)分量的視電導(dǎo)率受底部地層的影響也越大,因此從圖11中可以看出,三個(gè)不同源距測(cè)井響應(yīng)與背景響應(yīng)間的偏離程度隨著起伏界面高度的增加而變大,然而,源距L=5 m的探測(cè)深度相對(duì)較小,起伏界面高度8 m和12 m間的測(cè)井響應(yīng)已經(jīng)相差較小（圖11（a1—a5））,并且其各分量視電導(dǎo)率曲線(xiàn)形態(tài)與起伏界面形態(tài)勻是單峰曲線(xiàn),從而有助判斷界面形態(tài).此外,由于源距L=15和35 m的探測(cè)范圍較大,三個(gè)不同界面起伏高度H對(duì)應(yīng)的五個(gè)視電導(dǎo)率分量間的差異十分明顯（圖11（b1—b5）和圖10（c1—c5））,除了圖11（b4）和（b5）對(duì)應(yīng)的源距L=15 m交叉分量σa,zx和σa,xz以外,其他所有視電導(dǎo)率曲線(xiàn)均出現(xiàn)多峰現(xiàn)象,非常不利于層界面深度的反演.

3 結(jié)論

本文基于二次場(chǎng)耦合勢(shì)三維有限體積法建立了一套隨鉆超深前視電阻率儀器數(shù)值模擬算法,有效保證了在大計(jì)算區(qū)域內(nèi)電磁場(chǎng)的計(jì)算精度,與NMM的對(duì)比顯示其相對(duì)誤差小于3%,而在相同的剖分網(wǎng)格上,總場(chǎng)耦合勢(shì)三維有限體積法的計(jì)算精度明顯較差.此外,采用建庫(kù)插值方法計(jì)算背景電場(chǎng),在給出的算例中使背景電場(chǎng)的計(jì)算效率提高25倍,并能保證較好的計(jì)算精度,大大降低了在多發(fā)射源、大異常體情況下二次場(chǎng)方法的計(jì)算成本,提高了二次場(chǎng)方法在三維數(shù)值算法中的應(yīng)用價(jià)值.

數(shù)值結(jié)果進(jìn)一步證實(shí)了源距變大、工作頻率變小能增加儀器的探測(cè)范圍,探測(cè)范圍變大使得在層狀背景介質(zhì)或含有異常體、起伏界面情況時(shí),相應(yīng)的測(cè)井響應(yīng)變化特征會(huì)更為復(fù)雜,而短源距測(cè)量結(jié)果由于其探測(cè)范圍較小,在儀器離界面距離較近時(shí),其響應(yīng)往往是層邊界距離的單調(diào)變化函數(shù),有利于用于層邊界距離的反演;而不同源距交叉分量響應(yīng)均具有隨層邊界距離的單調(diào)變化特征,也有利于層邊界反演.

最后需要說(shuō)明的是,在三維地層條件下,整個(gè)儀器的響應(yīng)十分復(fù)雜,要想準(zhǔn)確確定遠(yuǎn)離井眼的地層中電導(dǎo)率的空間分布,須研究開(kāi)發(fā)出相應(yīng)的三維反演算法.

附錄A 電磁場(chǎng)分解

在井眼坐標(biāo)系下,三維各向異性地層中多分量隨鉆超深前視電阻率測(cè)井響應(yīng)實(shí)質(zhì)上可以表示為求解如下Maxwell方程的問(wèn)題（時(shí)諧因子取e-iω t）:

（A1）

其中,

（A2）

（A3）

p=x,y,z

（A4）

（A5）

附錄B 水平層狀TI介質(zhì)中并矢Green函數(shù)計(jì)算

（B1）

（B2）

（B3）

（B4）

為對(duì)方程（B3）和（B4）進(jìn)行二維Fourier逆變換確定頻率-空間域中的并矢Green函數(shù),定義如下形式Sommerfeld型積分:

（B5）

（B6）