井中磁源瞬變電磁三維時域有限差分數值模擬

孟慶鑫 ，潘和平

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

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

1 時域有限差分法基本原理

1.1 基本關系式與差分方程

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

式中：E為電場強度；H為磁場強度；σ為電導率；μ為磁導率。

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

1.2 場源和初始條件

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

1.3 邊界條件和步進時間

選取Mur吸收邊界條件，采用其二階近似式和一階近似式來處理截斷邊界上面的節點和棱角節點(6面12棱8角)，以x=0截面為例，二階近似式和一階近似式分別為[13?14]：

經推導可得截面和棱邊節點的差分格式，以 x=0面和x=y=0棱為例，分別為[13?14]：

2 全空間磁偶源 TEM 場響應特征模擬

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

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

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

2.1 均質全空間TEM場的一般響應特征與FDTD計算精度分析

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

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

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

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

2.2 三維板狀導體的瞬變場異常響應特征

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

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

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

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

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

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

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

圖8 不同情況激發的瞬變場響應的時間衰減曲線Fig.8 TEM responses for homogeneous medium with conductor at different conditions

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

3 結論

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

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

(3) 均質全空間下不同位置的瞬態場衰減情況不同，距場源越遠，其二次場衰減速度越快；對于不同參數的異常體，異常體受激發產生二次場越大，其與地層圍巖所產生場的相互作用過程越長，整個響應衰減速度越慢。

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

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

[1]牛之璉.時間域電磁法原理[M].長沙: 中南大學出版社,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].石油地球物理勘探, 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]閆述, 陳明生, 傅君眉.瞬變電磁場的直接時域數值分析[J].地球物理學報, 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]熊彬, 羅延鐘.電導率分塊均勻的瞬變電磁2.5維有限元數值模擬[J].地球物理學報, 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]岳建華, 楊海燕.礦井瞬變電磁法三維時域有限差分數值模擬[J].地球物理學進展, 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].地球物理學報, 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]宋汐瑾, 黨瑞榮, 郭寶龍, 等.井中磁源瞬變電磁響應特征研究[J].地球物理學報, 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]葛德彪, 閆玉波.電磁波時域有限差分方法[M].西安: 西安電子科技大學出版社, 2002: 8?128.GE Debiao, YAN Yubo.Finite difference time domain electromagnetic method[M].Xi’an: Xidian University Press,2002: 8?128.

[14]王長清, 祝西里.電磁場計算中的時域有限差分法[M].北京:北京大學出版社, 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].趙經祥, 譯.北京: 地質出版社, 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.