井中磁源瞬变电磁三维时域有限差分数值模拟

孟庆鑫 ，潘和平

(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.