井中磁源瞬變電磁三維時(shí)域有限差分?jǐn)?shù)值模擬

孟慶鑫 ，潘和平

(1.中國地質(zhì)大學(xué)(武漢) 地球物理與空間信息學(xué)院，湖北 武漢，430074；2.石家莊經(jīng)濟(jì)學(xué)院 勘查技術(shù)與工程學(xué)院，河北 石家莊，050031)

井中瞬變電磁法同井組合方式(陣列或偶?偶裝置)以置于井中的脈沖電流作為場源(通常為磁偶源)激發(fā)感應(yīng)電磁場，切斷場源后通過探頭或線圈接收地層介質(zhì)由一次場激發(fā)而產(chǎn)生的瞬態(tài)二次場，通過研究分析二次場特征獲取關(guān)于異常體或地層介質(zhì)的有用信息。由于脈沖信號(hào)一次激發(fā)后在不同時(shí)間點(diǎn)進(jìn)行測量相當(dāng)于在多個(gè)頻率點(diǎn)觀測結(jié)果[1]，故瞬態(tài)場較之單頻域電磁場而言包含了更豐富的信息，在實(shí)際工作中有效率高和場強(qiáng)大等優(yōu)點(diǎn)，得到了普遍的關(guān)注和研究。近年來，國外將瞬變電磁法用于裸眼井測井及地質(zhì)導(dǎo)向鉆井，并應(yīng)用于工礦領(lǐng)域，國內(nèi)也開始了相關(guān)儀器的研制和物理實(shí)驗(yàn)，應(yīng)用前景廣闊。關(guān)于瞬變電磁數(shù)值模擬技術(shù)，前人取得了許多成果[2?5]，國內(nèi)工作者也對(duì)多種方法進(jìn)行應(yīng)用和研討[6?8]，在全空間瞬變場方面，岳建華等[9]實(shí)現(xiàn)了礦井全空間電偶源瞬變場的三維時(shí)域有限差分模擬，楊海燕等[10]探討了礦井 TEM 視電阻率解釋方法，宋汐瑾等[11]對(duì)管套井中瞬變電磁響應(yīng)的理論計(jì)算進(jìn)行了推演，虞兵等[12]采用有限元法分析了過管套井的脈沖瞬態(tài)場特性。時(shí)域有限差分法(FDTD)能夠直觀簡潔地反映瞬態(tài)場傳播及與地質(zhì)體相互作用的物理過程，在相關(guān)問題研究中得到了廣泛使用和認(rèn)可。本文作者基于前人的研究[6?7,9,13?16]，通過應(yīng)用FDTD法對(duì)全空間瞬變電磁場傳播和三維低阻異常體響應(yīng)進(jìn)行模擬并簡析其特征規(guī)律。

1 時(shí)域有限差分法基本原理

1.1 基本關(guān)系式與差分方程

在無源、均勻、非磁性且各向同性的線性介質(zhì)中，忽略位移電流，準(zhǔn)靜態(tài)下近似無源Maxwell微分方程組為：

式中：E為電場強(qiáng)度；H為磁場強(qiáng)度；σ為電導(dǎo)率；μ為磁導(dǎo)率。

在直角坐標(biāo)系下，將求解區(qū)域剖分成均勻立體網(wǎng)格，使空間中連續(xù)電場值離散為各節(jié)點(diǎn)處電場，以任意節(jié)點(diǎn)電場Ei,j,k為中心與周邊6個(gè)節(jié)點(diǎn)所構(gòu)成立方體作為體積元(見圖 1(a))應(yīng)用式(5)進(jìn)行體積分。采用Gauss公式變體積分為面積分，并用中心差商近似替代電場分量對(duì)空間的偏導(dǎo)數(shù)的方法解出式(5)左端的差分格式，利用對(duì)空間電導(dǎo)率的網(wǎng)格剖分與 DuFort Frankel法對(duì)式(5)右端電場分量對(duì)時(shí)間的偏導(dǎo)數(shù)進(jìn)行近似替代的方法解出其差分格式，最終綜合兩端式子可得擴(kuò)散方程的七點(diǎn)差分格式，詳細(xì)推導(dǎo)過程見文獻(xiàn)[9]，這里僅列出最后結(jié)果：；Δxi，Δyj和 Δzk為各軸向網(wǎng)格距離；σi,j,k為節(jié)點(diǎn)電導(dǎo)率。

1.2 場源和初始條件

三維時(shí)域有限差分算法中，源可利用均勻全空間的解析式作為初始條件帶入，前提是異常體與源要有一定的距離，以滿足均勻空間的限制，源的形式也有很大自由度，可根據(jù)需要選擇[1]。這里選取磁偶源，其均勻全空間的電場分量和磁場分量解析式及推導(dǎo)見文獻(xiàn)[15]，本文對(duì)瞬態(tài)場特征模擬分析主要基于電場量，僅列出其解析式如下：

1.3 邊界條件和步進(jìn)時(shí)間

選取Mur吸收邊界條件，采用其二階近似式和一階近似式來處理截?cái)噙吔缟厦娴墓?jié)點(diǎn)和棱角節(jié)點(diǎn)(6面12棱8角)，以x=0截面為例，二階近似式和一階近似式分別為[13?14]：

經(jīng)推導(dǎo)可得截面和棱邊節(jié)點(diǎn)的差分格式，以 x=0面和x=y=0棱為例，分別為[13?14]：

2 全空間磁偶源 TEM 場響應(yīng)特征模擬

根據(jù)上文算法，模擬全空間三維磁偶源瞬變場。本文模型(見圖 1(b))參數(shù)的選取盡量保持一致性，三維直角坐標(biāo)系下全空間模型(長×寬×高為 50 m×50 m×50 m)；大地電導(dǎo)率為0.01 S/m，對(duì)整個(gè)介質(zhì)采取均勻網(wǎng)格剖分(100×100×100)，垂直井位于均質(zhì)模型中軸(x=y=0)；發(fā)射線圈Tx(磁偶源)置于井中，線圈半徑0.1 m，匝數(shù)30，電流強(qiáng)度10 A，接收線圈Rx歸一化；為簡便起見Rx與Tx位于同樣位置，該收發(fā)情況下沿井軸運(yùn)動(dòng)接收的響應(yīng)始終同號(hào)。

圖1 三維全空間均質(zhì)模型Fig.1 Full-space 3-D uniform model

圖2 均質(zhì)全空間中瞬變場(電場量E)全剖面等值線圖(0.378 μs)Fig.2 Electric field contour in homogeneous medium (0.378 μs)

2.1 均質(zhì)全空間TEM場的一般響應(yīng)特征與FDTD計(jì)算精度分析

圖2和圖3所示分別為不同時(shí)刻瞬變電磁場的電場量 E等值線圖。Tx位于介質(zhì)中心(x=0，y=0，z=0)點(diǎn)(見圖1(b))，電場量E與其在直角坐標(biāo)系下的2個(gè)電場分量Ey和Ex的關(guān)系見式(7)。由圖2和圖3可知均質(zhì)全空間中磁偶源 TEM 場的擴(kuò)散情況：均勻介質(zhì)中感應(yīng)渦流產(chǎn)生的 TEM 場以源點(diǎn)為中心隨時(shí)間的延遲向外擴(kuò)散傳播，TEM場強(qiáng)度隨時(shí)間的推移而衰減，瞬態(tài)場在三維均質(zhì)空間分布為含有1個(gè)環(huán)帶狀極大值的渦流場，即通常所說呈“煙圈”狀。圖 2(a)和圖 3(a)為y=0切面電場量等值線圖，由式(7)可知：在y=0面上Ex=0，可視電場量為Ey分量(圖中所示為場強(qiáng)，未作變號(hào)標(biāo)注)；圖2(b)和圖3(b)為z=0切面電場強(qiáng)度等值線圖，在此深度段即為該模型參數(shù)下三維空間中不同z切面的電場量極大值(圖中所示為電場強(qiáng)度)。

圖4所示為均質(zhì)全空間下瞬變場電場分量的解析解與時(shí)域差分解的對(duì)比圖。圖4(a)所示為位于點(diǎn)(0,3,0)的 Ex分量和位于點(diǎn)(?5,0,0)的 Ey分量在 0.201～2.509 μs的解析解和FDTD解(解析計(jì)算方法見式(7))；圖4(b)所示為2個(gè)電場分量FDTD解的相對(duì)誤差(以解析解為比對(duì)參考)。由圖4可知：FDTD計(jì)算結(jié)果可以在較早時(shí)段內(nèi)保持模擬計(jì)算結(jié)果的精度，計(jì)算誤差隨時(shí)步迭代次數(shù)的增加(時(shí)間推移)而逐步增大；除差分代替微分所造成的誤差外，吸收邊界條件與差分格式的計(jì)算精度差異也會(huì)造成誤差增大。在實(shí)際模擬過程中，選擇較小的時(shí)間步進(jìn)增量可保證更高的精度，但同時(shí)也增加了計(jì)算量。以試算結(jié)果為參考，本文選取較早時(shí)道和較小時(shí)步增量，以滿足研究工作的需要。

圖3 均質(zhì)全空間中瞬變場(電場量E)全剖面等值線圖(1.325 μs)Fig.3 Electric field contour in homogeneous medium (1.325 μs)

圖4 均質(zhì)全空間下磁偶源瞬變場電場分量解析解與FDTD 解(0.201～2.509 μs)Fig.4 Electric field contour in homogeneous medium(0.201～2.509 μs)

2.2 三維板狀導(dǎo)體的瞬變場異常響應(yīng)特征

如圖1(b)所示，建立包含三維低阻板狀體的均質(zhì)全空間模型；設(shè) 3種低阻板狀導(dǎo)體參數(shù)情況：(Condition-A) 板狀體電導(dǎo)率為 0.1 S/m，長×寬×高為10 m×10 m×1 m；(Condition-B) 板狀體電導(dǎo)率為0.1 S/m，長×寬×高為 14 m×14 m×1 m；(Condition-C)板狀體電導(dǎo)率為0.5 S/m，長×寬×高為10 m×10 m×1 m。Tx與Rx位于相同位置，沿井軸(z軸)進(jìn)行激發(fā)接收，接收由感應(yīng)渦流場所產(chǎn)生的磁感應(yīng)分量，包括垂直分量Bz和x向水平分量Bx；Tx與三維低阻板狀體的垂直距離為D，單位為m。

圖5(a)和(b)所示分別為均質(zhì)全空間中含1個(gè)水平低阻薄板狀體(圖 5中虛線所示，長×寬×高為 10 m×10 m×1 m，電導(dǎo)率為0.1 S/m；板狀體中心位于(15,0,5)點(diǎn)，Tx位于原點(diǎn)(0,0,0)，模型示意圖見圖1(b)：與正常場相比，圖5所示在低阻導(dǎo)體及附近的電場等值線發(fā)生明顯畸變，等值線密集，表明電磁波擴(kuò)散到低阻導(dǎo)體時(shí)擴(kuò)散速度降低，渦流在低阻導(dǎo)體中密度大，衰變慢，反映出瞬變場對(duì)低阻導(dǎo)體較敏感；圖 6(a)和(b)所示分別為均質(zhì)全空間中含1個(gè)水平薄板狀高阻體(圖6中虛線所示，電導(dǎo)率低于地層介質(zhì)電性參數(shù)，其他參數(shù)同圖5)：在高阻異常體及附近的電場等值線沒有發(fā)生明顯變化，反映出瞬變場對(duì)高阻異常體有較好的穿透性，實(shí)際計(jì)算結(jié)果上與正常場仍有差別。需說明的是：圖5(b)和圖6(b)中所示為電場強(qiáng)度。

圖5 均質(zhì)全空間中包含低阻板狀體的瞬變場(電場量E)全剖面等值線圖(0.497 μs)Fig.5 Electric field contour in homogeneous medium with conductor (0.497 μs)

圖6 均質(zhì)全空間中包含高阻板狀體的瞬變場(電場量E)全剖面等值線圖(0.497 μs)Fig.6 Electric field contour in homogeneous medium with resistor (0.497 μs)

圖7所示為Tx沿井軸移動(dòng)、Rx所測量到的渦流場產(chǎn)生的磁感應(yīng)分量半剖面曲線圖；異常體模型參數(shù)為 Condition-A，中心坐標(biāo)為(15,0,0)，對(duì)應(yīng) 0.497 μs時(shí)刻；曲線按接收響應(yīng)的最大值進(jìn)行歸一化，橫坐標(biāo)D為Tx與異常體的垂直距離。由圖7可知：該地電模型下井中 TEM 收發(fā)裝置接收的響應(yīng)曲線特征：Bz極大值出現(xiàn)在D=0處(激發(fā)源與異常體距離最小處)，響應(yīng)值隨著D(Tx與異常體距離)的增大而減??；Bx響應(yīng)極大值出現(xiàn)在 Tx離異常體一定距離處，出現(xiàn)極大值后響應(yīng)隨距離增加而逐步減?。辉趯?shí)際數(shù)值上垂直分量大于水平分量。

圖7 均質(zhì)全空間下板狀導(dǎo)體的半剖面異常響應(yīng)曲線(0.497 μs)Fig.7 TEM responses for whole-space of homogeneous medium with a conductor (0.497 μs)

圖8(a)所示為D(激發(fā)源Tx與導(dǎo)電板體的垂向距離)不同時(shí)，Rx接收的Bz響應(yīng)分量隨時(shí)間延遲而衰減的情況(為了對(duì)比，圖8中曲線是按各最后時(shí)刻的響應(yīng)值進(jìn)行了歸一化處理，采樣時(shí)間段為0.260～2.509 μs)：隨著激發(fā)源與異常體距離的增加，Rx接收的響應(yīng)衰減速度加快，說明磁偶源產(chǎn)生的瞬態(tài)渦流場的衰減速度與位置有關(guān)，距離場源點(diǎn)越遠(yuǎn)的場強(qiáng)衰減越快；同時(shí)反映出異常體距場源越遠(yuǎn)耦合情況越差，故其所激發(fā)的二次場較小，衰減較快。需說明的是：Bx分量的曲線形態(tài)較復(fù)雜，但大致規(guī)律相近。

圖8 不同情況激發(fā)的瞬變場響應(yīng)的時(shí)間衰減曲線Fig.8 TEM responses for homogeneous medium with conductor at different conditions

圖 8(b)所示為異常體的電性和尺寸參數(shù)不同時(shí)(中心坐標(biāo)為(15,0,0)，Tx位于(0,0,4))，Rx所接收的響應(yīng)隨時(shí)間延遲而衰減的情況(為了對(duì)比，圖8中曲線是按各最后時(shí)刻的響應(yīng)值進(jìn)行了歸一化處理，采樣時(shí)間段為0.260～2.509 μs)：均質(zhì)全空間條件下Rx所接收的響應(yīng)衰減最快，含低阻異常體時(shí)響應(yīng)的衰減速度變慢，這反映出相同激發(fā)條件下，導(dǎo)電板狀體的異常場與地層圍巖產(chǎn)生的場相互作用影響到總場的衰減速度；對(duì)比3種物性參數(shù)板狀體的異常響應(yīng)衰減情況可知，異常體的電導(dǎo)率越大，響應(yīng)的衰減越慢；異常體的尺寸越大，響應(yīng)的衰減越慢；說明異常體本身受激發(fā)產(chǎn)生的二次場越大(響應(yīng)與電導(dǎo)率及尺寸成正比)，其與地層圍巖所形成的場相互作用的過程就越長，故衰減越慢。需要注意的是：針對(duì)不同異常體，Bx分量的特征規(guī)律較復(fù)雜，受很多因素影響。

3 結(jié)論

(1) FDTD法是模擬分析瞬變電磁場特征規(guī)律的有效方法，在一定條件下(選取較小時(shí)步增量和較早時(shí)道)保證模擬的精度和效果。

(2) 模擬并簡述了均質(zhì)全空間磁偶源TEM場響應(yīng)情況；計(jì)算了該地電模型中不同參數(shù)導(dǎo)電板狀體的響應(yīng)情況，分析其特征規(guī)律。

(3) 均質(zhì)全空間下不同位置的瞬態(tài)場衰減情況不同，距場源越遠(yuǎn)，其二次場衰減速度越快；對(duì)于不同參數(shù)的異常體，異常體受激發(fā)產(chǎn)生二次場越大，其與地層圍巖所產(chǎn)生場的相互作用過程越長，整個(gè)響應(yīng)衰減速度越慢。

(4) 關(guān)于全空間磁偶源瞬變場問題的研究，以“渦流場”(如擴(kuò)散“煙圈”)為出發(fā)點(diǎn)，從電場分量的角度進(jìn)行解釋分析更為直觀簡便。

(5) 全空間磁偶源TEM響應(yīng)特征規(guī)律復(fù)雜，針對(duì)實(shí)際應(yīng)用存在很多問題(如不同地電條件下水平磁感分量特征規(guī)律、管套影響下對(duì)復(fù)雜地層介質(zhì)的勘查等)，采用的不同階Mur邊界條件存在計(jì)算精度差異，經(jīng)多次步進(jìn)迭代后模擬效果較差，需改進(jìn)算法；此外，在軸對(duì)稱介質(zhì)的井中 TEM 研究中，選用柱坐標(biāo)進(jìn)行模擬分析更為適用。

[1]牛之璉.時(shí)間域電磁法原理[M].長沙: 中南大學(xué)出版社,2007: 57?64.NIU Zhilian.The theory of time domain electromagnetic methods[M].Changsha: Central South University Press, 2007:57?64.

[2]Goldman M M, Stoyert C H.Finite-difference calculations of the transient field of an axially symmetric earth for vertical magnetic dipole excitation[J].Geophysics, 1983, 48(7): 953?963.

[3]Oristaglio M L, Hohmann G W.Diffusion of electromagnetic fields into a two dimensional earth: A finite-difference approach[J].Geophysics, 1984, 49(7): 870?894.

[4]Leppin M.Electromagnetic modeling of 3-D sources over 2-D inhomo-geneties in the time domain[J].Geophysics, 1992, 57(8):994?1003.

[5]Wang T, Hohmann G W.A finite difference time domain solution for three dimensional electromagnetic modeling[J].Geophysics,1993, 58(6): 797?817.

[6]宋維琪, 仝兆歧.3D瞬變電磁場有限差分正演計(jì)算[J].石油地球物理勘探, 2000, 35(6): 751?756.SONG Weiqi, TONG Zhaoqi.Forward finite differential calculation for 3D transient electromagnetic field[J].Oil Geophysical Prospecting, 2000, 35(6): 751?756.

[7]閆述, 陳明生, 傅君眉.瞬變電磁場的直接時(shí)域數(shù)值分析[J].地球物理學(xué)報(bào), 2002, 45(2): 275?284.YAN Shu, CHEN Mingsheng, FU Junmei.Direct time-domain numerical analysis of transient electromagnetic fields[J].Chinese J Geophys, 2002, 45(2): 275?284.

[8]熊彬, 羅延鐘.電導(dǎo)率分塊均勻的瞬變電磁2.5維有限元數(shù)值模擬[J].地球物理學(xué)報(bào), 2006, 49(2): 590?597.XIONG Bin, LUO Yanzhong.Finite element modeling of 2.5-D TEM with block homogeneous conductivity[J].Chinese J Geophys, 2006, 49(2): 590?597.

[9]岳建華, 楊海燕.礦井瞬變電磁法三維時(shí)域有限差分?jǐn)?shù)值模擬[J].地球物理學(xué)進(jìn)展, 2007, 22(6): 1904?1909.YUE Jianhua, YANG Haiyan.3D finite difference time domain numerical simulation for in-mine TEM[J].Progress in Geophysics, 2007, 22(6): 1904?1909.

[10]楊海燕, 鄧居智, 張華, 等.礦井瞬變電磁法全空間視電阻率解釋方法研究[J].地球物理學(xué)報(bào), 2010, 53(3): 651?656.YANG Haiyan, DENG Juzhi, ZHANG Hua, et a1.Research on full-space apparent resistivity interpretation technique in mine transient electromagnetic method[J].Chinese J Geophys, 2010,53(3): 651?656.

[11]宋汐瑾, 黨瑞榮, 郭寶龍, 等.井中磁源瞬變電磁響應(yīng)特征研究[J].地球物理學(xué)報(bào), 2011, 54(4): 1122?1129.SONG Xijin, DANG Ruirong, GUO Baolong, et al.Research on transient electromagnetic response of magnetic source in borehole[J].Chinese J Geophys, 2011, 54(4): 1122?1129.

[12]虞兵, 劉迪仁, 于新娟, 等.管套井脈沖源磁場特性的有限元分析[J].物探與化探, 2011, 35(5): 689?691.YU Bing, LIU Diren, YU Xinjuan, et al.The finite element analysis of magnetic field characteristics of the cased hole with pulse source[J].Geophysical & Geochemical Exploration, 2011,35(5): 689?691.

[13]葛德彪, 閆玉波.電磁波時(shí)域有限差分方法[M].西安: 西安電子科技大學(xué)出版社, 2002: 8?128.GE Debiao, YAN Yubo.Finite difference time domain electromagnetic method[M].Xi’an: Xidian University Press,2002: 8?128.

[14]王長清, 祝西里.電磁場計(jì)算中的時(shí)域有限差分法[M].北京:北京大學(xué)出版社, 1994: 47?95.WANG Changqing, ZHU Xili.Calculation of electromagnetic fields finite difference time domain method[M].Beijing: Beijing University Press, 1994: 47?95.

[15]納比吉安.勘查地球物理電磁法: 第1卷[M].趙經(jīng)祥, 譯.北京: 地質(zhì)出版社, 1992: 155?351.Nabighian M N.Electromagnetic methods in applied geophysics:Volume 1[M].ZHAO Jinxiang, trans.Beijing: Geological Publishing House, 1992: 155?351.

[16]Mur G.Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic field equations[J].IEEE Trans Electromagn Compat, 1981, 23:377?382.