杨 鑫，姜彦南，魏 兵

(1.西安电子科技大学 物理与光电工程学院, 陕西 西安 710071;2.西安电子科技大学 信息感知技术协同创新中心，陕西 西安 710071；3.桂林电子科技大学 信息与通信工程学院, 广西 桂林 541000)



·数理科学·

基于瞬变磁偶极子的正演算法

杨 鑫1,2，姜彦南3，魏 兵1,2

(1.西安电子科技大学 物理与光电工程学院, 陕西 西安 710071;2.西安电子科技大学 信息感知技术协同创新中心，陕西 西安 710071；3.桂林电子科技大学 信息与通信工程学院, 广西 桂林 541000)

以时域有限差分方法为技术手段,研究了瞬态电磁法中磁偶极子瞬态场的初始条件实现方法,提出了一种带边界条件的瞬变磁偶极源正演算法。文中先利用数值积分和Gaver-Stehfest变换等方法,实现了频域瞬态偶极源至时域空间分布的计算，采用高阶差分方法克服二阶差分方程带来的计算不稳定性。然后结合地-空边界和地下边界条件,获得瞬态磁偶极子源的初始场分布,模拟了电波在地下的扩散过程,并给出了地表的电动势响应曲线。算例的对比结果验证了文中算法的有效性。

瞬态场;磁偶极子;时域有限差分方法;正演计算

瞬变电磁法(Transient electromagnetic methods)是以电磁感应原理为基础的时间域人工场源电磁勘探方法。通过对地下响应信号的提取和处理,分析异常体的结构形态、导电性和埋藏深度等特性,达到找矿和解决地质问题等目的。相对于众多物探方法,瞬变电磁法具有探测深度大、分辨率高和适应性强等优点,广泛应用于工程以及环境地质勘探中[1]。

20世纪80年代之前，瞬变电磁法的研究主要限于1维和2维模型。之后，2.5维和3维模型成为研究热点[2-4]。在地电模型的计算中，根据地形或环境的不同常需采用不同形式模拟源，如矩形源、环形源或者阵列源等[5-6]。而针对不同类型源所激发的瞬态场，往往需采用不同的计算方式以获得相应的数学模型。并且，同一物理类型的源因激励方式的不同也需进行不同的数学推导，比如陈曙东等在2012年通过寻找冲激响应函数等方式对自由空间中的回线源在方波和梯形波等不同条件下分别推导出源的瞬变电磁响应[7]；2015年殷长春等人利用积分方程等方法推导出了均匀半空间中环形源在3维空间中的响应[8]；2015年Mark Goldman以及Kinza Haider等人研究了瞬变电磁法在海洋探测中的应用，并指出利用偶极方法在海洋探测中还处于理论探索阶段[9-10]。为了提高探测效率而发展出的混合探测法，比如瞬变电磁法与直流法(Direct current method)相结合[11]等。在这些研究中，瞬变电磁法对于一些简单环境的计算是有效的，但对于更复杂的地电模型如异常非规则体分布和山区地形等，这样的算法就受到了限制[12]。为有效地计算这些复杂的模型，就需要解决诸多问题，比如复杂瞬变源的计算、边界条件的确定等，但国内外目前关于瞬变激励源与边界条件的深入研究还不多。理论上讲，任意复杂源都可视为一定数量偶极子组成，因此，文中对磁偶极子源进行深入研究，提出了一种以频率域磁偶极子作为瞬变激励源并结合吸收边界条件的新方法。从频域偶极源出发，通过一定的数值计算手段实现瞬态场初始条件计算，极大地简化瞬态源计算以及程序实现的复杂性，为复杂源的加入及复杂环境下的瞬态场计算提供了条件。文中采用时域有限差分方法(FDTD)实现地下扩散场的计算。

1 瞬变电磁法基本理论

1.1 瞬态场扩散方程

在用瞬变电磁法处理地球物理问题时,大地可近似为无源、线性和各向同性的媒质。瞬变电磁波在地下扩散过程中,由于位移电流∂D(r,t)/∂t十分微弱,特别是在扩散后期,位移电流与传导电流J(r,t)相比可忽略不计。因此,地下扩散电波方程组可忽略位移电流项，其方程组为

(1)

(2)

(3)

(4)

在用FDTD对电磁场进行数值计算时,一般仅需式(1)和式(2)就可表征电磁波的运动过程。场在扩散晚期表现出似稳性,此时式(1)近似为

这样，式(3)就不能由式(1)导出,就导致独立方程不仅限于式(1)和式(2)。根据霍姆赫兹定理可知,要确定一个矢量,必须同时给定它的散度和旋度。若此时再用式(1)和式(2)对电磁场进行数值计算,就会因为B(r,t)的不确定性而导致计算错误。所以,在扩散场的后期,计算时就必须将式(3)显式地包含在内，才能保证计算的正确性。这样独立方程也就变为式(1)～(3),构成描述地下瞬态电磁波运动的方程组。

FDTD处理电磁计算问题时，其原理是将场方程转化为关于计算时间的显示差分迭代方程。因此，方程中必须存在场量对时间的偏导数。式(1)能满足了这一条件，但式(2)却不然。为解决这一问题，引入了Du Fort-Frankel方法。该方法在式(2)中引入虚拟位移电流[13]Jd(r,t)=γ∂E(r,t)/∂t后，变为

(5)

这样,结合式(1)，式(3)和式(5),就可实现地下电波的模拟计算。

1.2 边界条件

地层一般包含多种媒质(比如矿藏、水和气体等),就必定存在诸多不连续边界。而利用计算机处理的地电模型必定为有限区域,这种有限的区域具有的边界称为截断边界。所以,计算地电模型时将存在两种边界:不同媒质分界面和截断边界。由唯一性定理可知,要确定某区域的电磁场分布,就必须给定边界处的场分布。对于一些简单的地球物理模型,常用的地下边界处理办法是强制加狄利克雷边界条件,比如给定计算域边界处的磁场或者电场。但对于复杂的地下环境,确定地下边界处的场分布就显得比较困难,所以文中采用吸收边界条件方法处理地下截断边界。但对于地空边界条件的设置,由于瞬态源位于地面,若用吸收边界条件对其进行处理,则在地电模型中就相当于将源置位于吸收边界。这种设置方法将导致地表处场的计算错误,所以对于地空边界还需另行设置。在直角坐标系下,式(6)和式(7) 给出了地空边界的一种设置方法。取z轴正向指向自由空间,z=0为地表面,则磁场在地表处沿x和y方向的地空边界条件分别为

(6)

(7)

式(6)和式(7)给出的物理意义为:若知道波数域中地表处磁场垂直分量值,就可外推计算地表上高h处的磁场水平分量,即由内部场确定边界场，从而实现地空边界条件的计算。

2 时域有限差分方法

2.1 二阶差分法

以三维空间麦克斯韦方程组为例给出时域有限差分方法的基本计算思路。首先将电磁波传播的空间域划分为一定数量的方格。通过空间的这种离散划分,将电场和磁场在空间网格上进行交错排列便构成了时域有限差分方法的计算基础。电磁场各个分量的空间离散位置如图1所示。

图1 三维空间离散网格及其场分布Fig.1 3D FDTD grid and fields distribution

再将连续时间流t离散为一定数量的时间间隔(时间步),Δtn=tn+1-tn,n取正整数。tn代表时间的离散点,并约定在该整数离散时刻点计算电场值。用tn+Δtn/2代表每个时间步的中间时刻点,约定用该时刻点计算磁场值。按这种关系对式(1)、式(2)和式(3)进行离散,即可得电磁场的二阶差分迭代关系式。以式(1)的Bx分量为例:

2.2 高阶差分法

在电流源关断的早期,扩散场包含的频率十分丰富,但集中表现为高频段。所以对早期扩散场进行数值计算时,应该视情况而使用合适的高阶差分。文中选用四阶差分,由后面的计算结果可以表明,这种差分可以近似地计算早期电磁场的扩散过程。下面给出四阶差分的实现方法。

以1维可导函数f(x)为例,对f(x)进行关于x网格划分,每个网格编号为Δxi(i=1，2，3，…)。则f(x)在每个网格Δxi中点处Δxi-Δxi/2的一阶导数可由以下关系式求出，

f ′(x)|i+1/2=

a-1f(x)|i-1+a0f(x)|i+

a1f(x)|i+1+a2f(x)|i+2。

其中,参数ai(i=-1，0，1，2)可由如下方程组求出，

其中，yi为

将地下扩散方程按上式离散，可得其四阶差分离散形式。以Bx为例给出离散后的四阶差分方程，

类似地，可得磁场其他分量以及电场分量的四阶差分关系式。得到差分关系后,只要离散参数间满足空间以及时间的计算稳定性要求，就可以实现扩散方程的数值计算。

3 源的验证和分析

3.1 频域激励源

在圆柱坐标系下,位于均匀大地表面的磁偶极子源关断瞬间,电磁场在地下的频域解[14]为

(8)

Er(ω，ρ)=

(9)

(10)

(11)

由式(11)理论上能够计算出时空域的场分布，从而求出场的初始条件。但由于存在emz大积分因子,在计算机实际计算时,计算的空间大小一般可达数百或数千米,可见在这种计算尺度下,该积分因子的计算结果将容易溢出。而对一般个人计算机硬件资源而言,计算区域的宽度和深度是十分有限的。因此,以时域解析解为初始条件存在一定的局限性。为克服这一难点,文中采用数值积分法和Gaver-Stehfest方法求解时域初始场分布。

3.2 源的数值计算

源的频域方程式(8)～(10)为典型的Hankel积分变换形式,因此,可将其简写为

(12)

其中，x> 0，Jr为r阶贝塞尔函数。对于这种类型的积分计算，可采用高斯数值积分方法实现。而对于频域源中复杂的拉普拉斯变换，文中选用Gaver-Stehfest[15]变换法进行数值变换。由于Gaver-Stehfest变换方法是纯实数域的运算，且只需要少量的离散拉普拉斯变量值S(一般取12～16个)就可实现反变换的计算，所以该方法广泛应用于求解拉普拉斯逆变换问题。其变换思路大致为：设所求逆变换的时刻为t，先将S域中的拉普拉斯变量离散为Sm= m·ln(2)/t，再把Sm带回拉普拉斯式中，可得离散的拉普拉斯函数F(Sm)。最后，用式(13)表示的Gaver-Stehfest变换式求函数F(s)拉普拉斯逆变换值f(t)，

(13)

其中,Km为Gaver-Stehfest变换系数,n为变换所需S域的离散点数。由式(13)可见,Gaver-Stehfest变换可极大地降低拉普拉斯变换计算的复杂性。这样,通过Hankel积分变换求出频率域场表达式,然后再对频率域解作逆拉普拉斯变换,即可求得时空域的初始场分布,计算结果如图2所示。图中显示了当T=6 ms时电磁场部分分量在地下的纵向分布情况。

由图2可见,在源关断初期,场的能量集中分布在源附近,计算结果与瞬态源关断瞬间的实际场分布相符合。电场随着与源点的距离增加而呈明显的高次衰减分布,且主要分布于地表浅层。此外,电场中所携带的能量明显低于磁场中所携带的能量,并且当远离源时能量主要存在于磁场中,符合电磁能量在导体中的分布特点。图2中还清晰地显示磁场分量分布的“烟圈”效应,这些“烟圈”中心的位置随着剖面与源的距离加大而逐渐加深,体现了瞬变场的扩散特点。

图2 地下场纵向剖面图Fig.2 Longitudinal crossection of underground fields

3.3 激励源数值实现及验证

图3～5给出以频率域磁偶极子源作为初始条件的地下瞬变电磁场的计算结果。图3给出两个不同时刻的地表磁场的垂直分量响应电动势数值解及时域解析解对比,可见数值计算结果与解析解能很好地吻合。由于时域解析源是关于1/ρ(ρ为观察点至源点的距离)的函数[15],可见其在源点处存在奇异性,所以场在源点会出现不稳定性,这种不稳定性随着场的迭代计算而逐渐向周围扩散。而数值解是通过网格离散来实现,就可以避开源点的奇异性,实现稳定的计算。图4给出了在地表x轴上50m和200m处的磁场垂直分量的电动势响应随时间的变化情况。图5给出了不同时刻磁场垂直分量在剖面y=0处的地下分布图,很好地显示了地下场的扩散情况。在扩散初期,磁场垂直分量的能量主要集中在地表源附近。随着时间的增加,场逐步往地下和水平方向扩散,扩散过程显示出了明显的“烟圈”效应。在场计算的后期,磁场垂直分量的“烟圈”效应逐渐消失于整个计算域,瞬态场完成在地下的扩散。

图3 数值解与解析解对比曲线Fig.3 Comparison of numerical solutions with analytical solutions

图4 地表磁场垂直分量响应曲线Fig.4 Ground surface response curve of the vertical magnetic field

图5 不同时刻地下磁场垂直分量分布Fig.5 Vertical component distribution of underground magnetic fields at different time

4 结 语

从磁偶极子的瞬变场频域特性出发,计算电波在地下的时域场分布,提出了一种结合吸收边界条件的瞬变场初始条件计算方法。先分别利用Gaver-Stehfest变换方法和高斯数值积分方法求解频域磁偶极子源在地下的时域场分布,然后由地表场外推自由空间中的场以确定地空边界条件,最后利用FDTD实现了瞬变偶极场的电动势响应，计算并给出了不同时刻地下场的分布图,清晰地显示了地下扩散场的“烟圈”效应。由验证结果可知，文中正确地模拟了瞬变磁偶极子源在地下的扩散过程。该算法可将源的复杂计算过程简化为基本的数值计算,避免了大积分因子的计算,极大地降低了计算的复杂度。文中算法还具有通用性较强的特点,可为复杂瞬变源的模拟提供可靠的参考。

[1] 丁艳飞, 白登海, 许诚. 均匀半空间表面大定源瞬变电磁响应的快速算法[J].地球物理学报,2012,55(6):2087-2096.

[2] SANFILIPO W A, EATON P A, HOHMANN G W. The transient EM response of a prism in a conductive half-space[J].Geophysics, 1985, 50(2): 272.

[3] ENDO M, NOGUCHI K.Three-dimensional modeling considering the topography for the case of the time-domain electromagnetic method[C].2nd International Symposium on 3-Dimensional Electromagnetics, 1999.

[4] MAAO F A. Fast finite-difference time-domain modeling for marine-subsurface electromagnetic problems[J].Geophysics, 2007, 72: 19-23.

[5] UME S, HARRIS J M, ALUMBAUGH D L. A finite element algorithm for 3-D transient electromagnetic modeling[C]∥2009 SEG Annual Meeting. Society of Exploration Geophysicists, 2009.

[6] 孙怀凤, 李貅, 李术才，等. 考虑关断时间的回线源激发TEM三维时域有限差分正演[J].地球物理学报, 2013, 56(3): 1049-1064.

[7] 陈曙东, 林君, 张爽. 发射电流波形对瞬变电磁响应的影响[J]. 地球物理学报, 2012, 55(2): 709-716.

[8] 殷长春, 任秀艳, 刘云鹤. 航空瞬变电磁法对地下典型目标体的探测能力研究[J]. 地球物理学报, 2015, 58(9): 3370-3379.

[9] GOLDMAN M, MOGILATOV V, HAROON A,et al. Signal detectability of marine electromagnetic methods in the exploration of resistive targets[J]. Geophysical Prospecting, 2015, 63(1):192-210.

[10] HAIDER K, ENGESGAARD P, SONNENBORG T O,et al. Numerical modeling of salinity distribution and submarine groundwater discharge to a coastal lagoon in Denmark based on airborne electromagnetic data[J]. Hydrogeology Journal, 2015, 23(2): 217-233.

[11] CHENG Jiulong, LI Fei, PENG Suping, et al. Joint inversion of TEM and DC in roadway advanced detection based on particle swarm optimization[J]. Journal of Applied Geophysics, 2015, 123: 30-35.

[12] 魏明君, 丁云河, 李冰,等. 线源瞬变电磁法在熊耳山区多金属找矿中的应用研究[J].地球物理学进展, 2015, 30(1): 146-152.

[13] WANG T, HOHMANN G W. A finite-difference time-domain solution for three-dimensional electromagnetic modeling[J]. Geophysics, 1993, 58(6): 797-809.

[14] 米萨克·纳比吉安. 勘查地球物理电磁法第一卷理论[M]. 北京: 地质出版社, 1992.

[15] 宋汐瑾, 党瑞荣, 郭宝龙,等. 井中磁源瞬变电磁响应特征研究[J]. 地球物理学报, 2011, 54(4):1122-1129.

(编 辑 李 静)

The transient properties of magnetic dipole source

YANG Xin1,2, JIANG Yannan3, WEI Bing1,2

(1.School of Physics and Optoelectronic Engineering,Xidian University, Xi′an 710071, China;2.Collaborative Innovation Center of Information Sensing and Understanding, Xidian University, Xi′an 710071， China；3.School of Information and Communication Engineering, Guilin University of Electronic Technology, Guilin 541000, China)

The calculation method of initial condition for transient magnetic dipole is studied by means of FDTD, and a forward algorithm for its transient fields is proposed with a boundary condition. The numerical integration and Gaver-Stehfest transform are applied to achieve the magnetic dipole fields from frequency domain to time domain. And the high-order difference method is used to overcome the calculation instability of second order differential equations. Combinating with the ground-air boundary and subsurface boundary conditions, the initial distributed field of transient magnetic dipole is obtained. The diffusion of the EM wave of the underground is simulated and the EMF response curve at the earth′s surface is presented as well. Finally, the numerical results show the validation of the algorithm.

transient field; magnetic dipole; FDTD method; forward calculation

2015-09-30

国家自然科学基金资助项目(61231003，61401344，61571348)；广西自然科学基金资助项目(2014GXNSFAA118283)

杨鑫，男，重庆人，博士生，从事计算电磁学研究。

P631

A

10.16152/j.cnki.xdxbzr.2016-06-005