地球物理学进展  2016, Vol. 31 Issue (1): 268-273   PDF    
引用本文
刘晓, 谭捍东. 2016. 瞬变电磁法全空间三维伪谱法模拟[J]. 地球物理学进展, 31(1): 268-273, doi: 10.6038/pg20160130  
LIU Xiao, TAN Han-dong. 2016. 3D pseudo-spectral method for TEM modeling in whole-space. Progress in Geophysics, 31(1): 268-273, doi: 10.6038/pg20160130   
瞬变电磁法全空间三维伪谱法模拟
刘晓, 谭捍东     
中国地质大学(北京)地球物理与信息技术学院, 北京 100083
收稿日期: 2015-03-09 修回日期: 2015-06-16.
基金项目: 国家自然科学基金项目(41374078)和国土资源部地质调查项目(12120113086100,12120113101300)联合资助.
作者简介: 刘晓,男,1987年11月生,河南人,中国地质大学(北京),博士研究生,主要从事瞬变电磁的模拟研究.(E-mail:friendat2010@163.com)
通讯作者: 谭捍东,男,1966年生,中国地质大学(北京),教授,主要从事地球物理的研究与教学工作.(E-mail:thd@cugb.edu.cn)
摘要: 本文根据均匀介质中时间域麦克斯韦方程组,应用伪谱法模拟三维全空间瞬变电磁场的扩散.将电磁场的六分量整理成统一的方程,通过初值条件替代源项将源响应问题转化为初值问题.使用快速傅里叶变换近似空间偏导数,应用切比雪夫多项式近似演化算子来求解方程.通过与全空间均匀介质瞬态线源响应的解析解的比较,验证伪谱法的有效性.计算和分析不同时刻全空间均匀介质中高阻体和低阻体的响应,说明瞬变电磁法对高阻体灵敏度不高,但能较好的分辨低阻体的特点,并反映出时间域伪谱法的优点,即可以连续地计算出电磁场扩散及与低阻体的响应过程.最后讨论伪谱法未来可能应用于矿井瞬变电磁和海底瞬变电磁等领域.
关键词: 瞬变电磁场     时间域电磁法     伪谱法    
3D pseudo-spectral method for TEM modeling in whole-space
LIU Xiao, TAN Han-dong     
School of Geophysics and Information Technology, China University of Geosciences, Beijing 100083, China
Abstract: Based on the 3D Maxwell's equation in time domain, the pseudo-spectral method is used to simulate electromagnetic propagation in whole-space in time domain in this paper. The six components of electromagnetic field are transformed to same form. By means of replacing source with initial condition, the problem of respond with source is transformed to the problem with initial condition. Then the fast Fourier transforming is used to compute spatial derivatives and the Chebyshev expansion to approximate evolution operator. To verify the method, the numerical solutions of response of transient line in whole-space are compared with analytical solutions. The numerical solutions of response of a transient electric dipole in whole-space shows that transient electromagnetic method(TEM) has high resolution to low resistance and low resolution to high resistance. And the advantage of pseudo-spectral method is that it can simulate electromagnetic propagation and the response of the low resistance continuously. Finaly, it discusses the potential application fields of pseudo-spectral method are mine and seafloor TEM.
Key words: transient electromagnetic field     TEM modeling in time domain     pseudo-spectral method    
0 引 言

直接在时间域进行电磁扩散模拟是研究瞬变电磁法的一种策略，相比于先计算频率域的结果再变换到时间域，它避免了变换带来的误差，计算精度高，尤其是晚期.常用的电磁场时域数值模拟方法包括时域有限差分法、时域有限元法、时域积分方程法和伪谱法等.有限差分时间分步迭代法已经成为经典算法(Yee，1966; Wang and Hohmann，1993)，在某一时刻，它先用差分方程离散偏微分方程，计算出当前时刻的场值，然后递推下一时刻的场值(牛之琏，2007; 岳建华等，2007)，它要求苛刻的时间稳定性条件.有限元法将麦克斯韦方程组转化为微分方程，空间上用有限元离散插值，时间偏导数用差分近似来求解方程(Sugen et al.，1993; Stalnaker，2004).有限差分法和有限元法的计算精度不高.时域积分方程法的优点是不需要人为地设置边界条件，缺点是计算复杂，适合于模拟简单的异常体(Bennett，1968;李建慧等，2012).时间域的伪谱法是另一种高效、高精度的数值方法，它利用快速傅里叶变换近似空间导数(Fornberg，1996;Carcione，1999)，和其他方法相比，它的计算精度高、计算速度快，很适合于大规模模型的计算.

1986年Tal-Ezer(1986)提出以切比雪夫多项式为基础的伪谱方法来解双曲线方程，随后该方法得到更多的应用，如求线性周期抛物线问题(Tal-Ezer，1989)，解弹性力学方程(Kosloff and Kessler，1990;刘鲁波等，2007;龙桂华等，2009)，求解麦克斯韦方程(Druskin and Knizhnerman，1994;De Raedt et al.，2003; 李展辉等，2009)等.Carcione(2006)把伪谱法应用在二维麦克斯韦方程组的求解，使用显式谱方法来模拟低频信号的在地表的传播，通过交错傅里叶伪谱法来近似空间偏导数，他的计算结果与解析解较好地一致.

本文把伪谱法应用到全空间三维时间域瞬变场的模拟.在研究海洋瞬变电磁和矿井瞬变电磁等方法时(Carcione and Poletto，2003)，分析全空间中电磁源的理论响应具有直接的现实意义，因为这种情况下，大地空气分界面的影响可以忽略(纳比吉安，1992).

1 麦克斯韦方程组

在均匀无源各向同性介质中，麦克斯韦方程组变成如下形式(S and ers and Reed，1986):

其中∂i是时间偏导数，e和h分别表示电场和磁场，μ是磁导率，σ是电导率，▽是哈密顿算子.

方程(1)、(2)展开后变为

其中∂x、∂y和∂z分别为X方向、Y方向和Z方向的空间偏导数，

由磁感应强度的散度为零(▽·b=0)，得到

即∂xhx=-∂yhy，∂yhy=-∂zhz-∂xhx，∂zhz=-∂xhx-∂yhy，

带入方程4推导得到:

对于电场，当电流密度的散度为零时，即▽·J=▽·(σE)=0，做类似的代换，得到:

方程(6)到(11)具有相同的形式:

其中，w将称为场值向量，G=μ-1σ-1(∂x∂x+∂y∂y+∂z∂z)称为空间偏导数算子.

为求解方程(12)，需要考虑初始条件，若给定初始条件w0，则方程(12)的解为

用N=Nx·Ny·Nz个节点均匀离散方程(13)后得到:

其中N为节点总数，Nx、Ny、Nz分别为X、Y、Z方向剖分的节点数.选择合适的多项式近似算子etGN即得到方程(12)的数值解.

2 伪谱法原理2.1 傅里叶变换近似空间偏导数

对均匀离散的函数u(xj)，j=1，2，…，N，使用快速傅立叶变换将其变换到波数域，r=1，2，…，N，乘上波数系数，使用快速傅里叶反变换变换到空间域，即得到u(xj)的空间偏导数∂xu(xj)(Canuto et al.，1987; Fornberg，1996)，步骤如下:

这里FFT和FFT-1分别是快速傅里叶正反变换，，kr= 为空间波数，Δx为函数u(xj)的剖分间隔.

对传播矩阵GN，应用(15)式表示的过程来计算剖分网格间的导数. 如∂z∂zex，先沿Z方向对ey作FFT变换，用ikr乘以变换结果，对相乘结果作反FFT变换，得到∂zey，再对∂zey沿Z方向重复以上的运算，即得到∂z∂zey.

在网格上使用傅里叶伪谱法近似空间偏导数时，理论上具有无限高的精度，使用快速傅立叶变换效率很高(Carcione，2010).

2.2 切比雪夫多项式近似算子etGN

切比雪夫多项式近似算子etGN的方程为(Abramowitz and Stegun，1972; Tal-Ezer，1989; Carcione，2006)

其中，bk(t)=cke-btIk(bt)，ck是系数，b是矩阵GN的特征值的绝对值，Ik(bt)是k阶贝塞尔函数，FN=(GN+bI)/b，I为单位矩阵，Tk是k阶切比雪夫多项式，Tk(x)=cos(karccosx)，-1≤x≤1，arccosx为反余弦函数.

将方程16带入方程(14)得到

切比雪夫多项式由下列循环来计算，公式为

在波数域，矩阵>GN的特征值为

其中kx，ky，kz为波数.λ是实数，且位于数轴负半轴.对kx，ky，kz而言，其最大值都是奈奎斯特波数，kxmax=π/Δx，kymax=π/Δy，kzmax=π/Δz，这里，Δx，Δy，Δz分别为X，Y，Z方向的剖分间隔.

为保证方程(17)收敛，取b为

实际计算方程(17)，要对无穷累加进行截断，为保证稳定性，多项式的阶M应满足 ，根据Tal-Ezer(1989)的经验，取 .于是方程(17)变为

随M的增大，截断产生的误差以指数关系衰减.伪谱法要求Δt=o(1/N)作为时间稳定性条件.它克服了两个缺点:低精度和苛刻的稳定性条件(Tal-Ezer，1989).

3 源的处理

阶跃源是瞬变电磁法的常用源，它是在导线或线圈中供以稳定电流，在t=0时刻关闭电流，地下会感应出电磁场.在t=t0>0时刻，将初始时刻的电磁场值作为初值条件来替代源，Oristaglio和Hohmann(1984)和Wang和Hohmann(1993)已应用这种思想到瞬变电磁时域三维有限差分法模拟中，形成瞬变电磁时域模拟处理源的基本方式.

在全空间均匀介质中，阶跃源的响应容易计算，用初始时刻的初值条件替代源，将瞬变电磁阶跃源响应问题转化为初值问题，用伪谱法求解.原则上选取的初始时间要适当，既要使得全空间的解有效，又要使得电磁场有足够的采样.常使用如下初始时间:

(Wang and Hohmann，1993)，其中μ1是均匀全空间的磁导率，σ1是均匀全空间的电导率，Δmin是最小的网格间距.此时的烟圈相当于磁偶极子的等效烟圈穿透深度为1.5Δmin时的时间.

4 边界条件

电磁场扩散方程适用于无限介质，在有限的范围内模拟就需要人为地设置边界条件.傅里叶变换具有周期性，计算自身会产生环形伪影现象，于是，在边界处使用吸收衰减的边界条件，来消除这种非人为的干扰(Kosloff and Baysal，1982，Kosloff and Kosloff，1986; Carcione，2006，2007).

选取吸收因子λ，λ>0，在模拟区域的四个边界带上，用GN-λI代替GN，在边界的外侧λ具有最大值，向内侧去，λ逐渐减小，它的作用相当于在矩阵中GN增加了一个负势.使用这种从外到里λ逐渐减小的过程来消除环形伪影现象.所选择的λ既要消除反射影响，又要消除透射影响.

边界吸收衰减的方法虽然消除了环形伪影现象，但在边界上造成了误差.后续的工作是改进边界条件，既要考虑大地空气边界条件的影响，又能提高边界处的计算精度，为更广泛的实际应用提供理论依据.

5 模型算例5.1 全空间均匀介质瞬态线源响应

对于全空间均匀介质，取电阻率为100 Ω·m，磁导率为4π×10-7 H/m，介质中心为坐标原点，垂直方向为Z轴，水平方向为X轴和Y轴，沿Y轴放置无限长导线，供以10 A的电流，稳定后切断电源，介质中产生瞬变电磁场.由于沿Y方向对称，三维问题可以简化为二维问题，在XOZ平面，X和Z方向剖分节点nx、nz都取120，网格间距都取10 m，吸收边界取8个网格.

阶跃线源产生的瞬态电场Y分量为

其中θ2=μσ/4t，ρ2=x2+z2 .

根据方程(22)计算的t0为1.42 μs，为保证瞬变场场有足够的采样，将20 μs时刻的线源瞬变电磁场值作为初始条件，应用伪谱法计算得到不同时刻的电磁场值.图 1为40 μs时刻的XOZ断面的电场Y分量等值线图，电场等值线环绕线源对称分布.图 2为40 μs、100 μs和300 μs时刻，沿OY方向上的数值解与解析解的对比图，两者一致性非常好，均方误差分别为1.959×10-9、7.9694×10-9和7.9694×10-9，说明伪谱方法正确有效.

图 140 μs时刻的电场Y分量等值线图
Fig. 1 Y component of Electrical field contour at 40 μs

图 2 X轴方向上伪谱法计算结果与解析解的对比图
 Fig. 2 Comparison of PSM and the analytical solution of electrical field along X axis

在电场扩散的早期，从40 μs到100 μs，电场随时间衰减较快. 100 μs后，电场逐渐扩散到边界带上，受到吸收边界影响，到300 μs时，边界带上(-600，-520)和(520，600)的数值解小于解析解.

5.2 全空间均匀介质高阻体和低阻体的特征响应

全空间均匀介质电阻率为10 Ω·m，磁导率为4π×10-7 H/m，中心点取为坐标原点，磁偶极源位于YOZ平面上且放在原点，磁距为400 A·m2，供以稳定电流，在t=0时刻切断电流.模型如图 3，原点正上方100 m处对称放置正方体，边长为80 m，电阻率为1 Ω·m，下方对称位置放置形状相同电阻率为100 Ω·m的高阻体.以80×80×80的网格剖分研究区域，网格间距为10 m，吸收边界取8个网格.

图 3 全空间模型图
 Fig. 3 3Dmodel in whole-space

根据方程(22)计算的初始时间t0为14.2 μs，将100 μs时刻的电磁场值作为初值条件，计算不同时刻的电磁场值.图 4为100 μs时刻YOZ平面的瞬变电场Y分量等值线图，感应出的电场集中在源的附近，呈对称分布，没有扩散进入异常体.

图 4 初始时刻t0=100 μs电场Y分量等值线图
 Fig. 4 Electrical field contour at initial instant 100

图 5为不同时刻电场Y分量的等值线图，图a到图f分别为0.2 ms、0.3 ms、0.6 ms、0.9 ms、1.7 ms、3.3 ms、4.9 ms、8.1 ms时刻的等值线图.在扩散的早期(图a、b、c)，低阻体阻碍场的扩散，使电场发生畸变，电场扩散速度变低，衰减变慢.图a、b中，电场不能穿透低阻体，电场等值线绕开低阻体，高阻体对电磁场的扩散影响不大.图c中，电场极值等值线已经完全穿过高阻体而正在穿过低阻体.这验证了瞬变场穿透高阻能力强，对高阻体灵敏度不高，但能较好的分辨低阻体的特点.图d、e、f、g、h中，电磁场扩散逐渐过渡到晚期，低阻体与电磁场发生相互作用，感应出二次电场.图d中，低阻体底端先有明显的响应产生，也就是产生了涡旋电流. 图e中，响应产生的电场向低阻体中心移动. 图f、g、h中，响应产生的电场逐渐主导整个电场的分布，这部分电场主要是低阻体感应出的二次场.随时间的推移，图f中感应出的电场完全主导了整个电场的分布.图g中，高阻体中的扩散电场已经衰减掉. 图h中，低阻体响应产生的二次电场开始向高阻体扩散.与初始场相比，低阻体响应产生的二次场虽然幅值相对小，但持续的时间较长，瞬变电磁法正是研究它随时间的变化规律.

图 5 全空间介质中不同时刻电场Y分量的等值线图
 Fig. 5 Y component of Electrical field contour at different instants in whole-space medium

综上，伪谱法可以连续计算出电磁场扩散及低阻体的响应全过程，给出非常丰富的时间域的信息.

6 结 论

6.1     本文推导三维的麦克斯韦方程组，应用伪谱法模拟全空间介质中瞬变电磁场的三维响应.瞬变电磁法属于时间域电磁法，由伪谱法计算结果分析可知，时间域电磁法可以给出丰富的时域信息，连续观测电磁场扩散及低阻体的响应过程.

6.2     伪谱法由傅里叶变换计算空间偏导数，优点是计算效率高、精度高和低的稳定性条件，其缺点是难以加入大地空间边界条件，因为傅里叶变换具有全局性.它未来的潜在应用领域是矿井瞬变电磁和海底瞬变电磁.如果解决大地空气边界条件的问题，就可以应用到计算地表瞬变电磁场.

致 谢   感谢Tel-Aviv大学的Carcione J. M教授先前的研究，感谢各位审稿专家和各位编辑对本文的宝贵意见和建议.

参考文献
[1] Abramowitz M, Stegun I A. 1972. Handbook of mathematical functions[M]. New York:Dover Publications, 71-72.
[2] Bennett C L Jr. 1968. A technique for computing approximate electromagnetic impulse response of conducting bodies[Ph. D. thesis]. West Lafayette:Purdue University.
[3] Canuto C, Hussaini M Y, Quarteroni A, Zang T A Jr. 1987. Spectral methods in fluid dynamics[M]. New York:Springer-Verlag, Inc.
[4] Carcione J M. 1999. Staggered mesh for the anisotropic and viscoelastic wave equation[J]. Geophysics, 64(6):1863-1866.
[5] Carcione J M. 2001. Wave fields in real media:Wave propagation in anisotropic, anelastic and porous media[A].//Helbig K, Treitel S, eds. Handbook of Geophysical Exploration:Seismic Exploration, vol. 31[M]. Pergamon:Pergamon Press Inc.
[6] Carcione J M, Poletto F. 2003. Electric drill-string telemetry[J]. Journal of Computational Physics, 186(2):596-609.
[7] Carcione J M. 2006. Geophysical software and algorithms a spectral numerical method for electromagnetic diffusion[J]. Geophysical, 71(1):I1-I9.
[8] Carcione J M. 2007. Wave fields in real media:Wave propagation in anisotropic, anelastic, porous and electromagnetic media[A].//Handbook of Geophysical Exploration:Seismic Exploration[M]. Amsterdam:Elsevier Ltd, 38:380-426.
[9] Carcione J M. 2010. Simulation of electromagnetic diffusion in anisotropic media[J]. Progress in Electromagnetics Research B(PIER B), 26:425-450.
[10] De Raedt H, Michielsen K, Kole J S, et al. 2003. Solving the Maxwell equations by the Chebyshev method:A one-step finite-difference time-domain algorithm[J]. IEEE Transactions on Antennas and Propagation, 51(11):3155-3160.
[11] Druskin V, Knizhnerman L. 1994. Spectral approach to solving three-dimensional Maxwell's diffusion equations in the time and frequency domains[J]. Radio Science, 29(4):937-954.
[12] Fornberg B. 1996. A practical guide to pseudospectral methods[M]. Cambridge:Cambridge University Press.
[13] Kosloff D D, Baysal E. 1982. Forward modeling by a Fourier method[J]. Geophysics, 47(10):1402-1412.
[14] Kosloff R, Kosloff D. 1986. Absorbing boundaries for wave propagation problems[J]. Journal of Computational Physics, 63(2):363-376.
[15] Kosloff D, Kessler D, Filho A Q, et al. 1990. Solution of the equations of dynamic elasticity by a Chebychev spectral method[J]. Geophysics, 55(6):734-748.
[16] Li J H, Zhu Z Q, Zeng S H, et al. 2012. Progress of forward computation in transient electromagnetic method[J]. Progress in Geophysics(in Chinese), 27(4):1393-1400, doi:10.6038/j.issn.1004-2903.2012.04.013.
[17] LI Z H, Huang Q H, Wang Y B. 2009. A 3-D staggered grid PSTD method for borehole radar simulations in dispersive media[J]. Chinese J. Geophys.(in Chinese), 52(7):1915-1922, doi:10.3969/j.issn.0001-5733.2009.07.027.
[18] Liu L B, Chen X F, Wang Y B. 2007. Chebyshev pseudospectral method for seismic wavefield modeling[J]. Northwestern Seismological Journal(in Chinese), 29(1):18-25.
[19] Long G H, Li X F, Zhang M G. 2009. The application of staggered-grid Fourier pseudo-spectral differentiation operator in wavefield modeling[J]. Chinese J. Geophys.(in Chinese), 2009, 52(1):193-199.
[20] Nabighian M N. 1992. Electromagnetic methods in applied geophysics:Voume 1.(in Chinese)[M]. Zhao J X, Wang Y J Trans. Beijing:Geological Publishing House,(1):302-348.
[21] Oristaglio M L, Hohmann G W. 1984. Diffusion of electromagnetic fields into a two-dimensional earth:A finite-difference approach[J]. Geophysics, 49(7):870-894.
[22] Sanders K F, Reed G A L. 1986. Transmission and propagation of electromagnetic waves[M]. 2nd ed. Cambridge:Cambridge University Press.
[23] Stalnaker J L. 2004. A finite element approach to the 3D CSEM modeling problem and applications to the study of the effect of target interaction and topography[Master Degree Thesis]. Texas, USA:Texas A & M University.
[24] Sugeng F, Raiche A, Rijo L. 1993. Comparing the time-domain EM response of 2-D and elongated 3-D conductors excited by a rectangular loop source[J]. Journal of Geomagnetism and Geoelectricity, 45(9):873-885.
[25] Tal-Ezer H. 1986. Spectral methods in time for hyperbolic equations[J]. SIAM Journal on Numerical Analysis, 23(1):11-26.
[26] Tal-Ezer H. 1989. Spectral methods in time for parabolic problems[J]. SIAM Journal on Numerical Analysis, 26(1):1-11.
[27] Wang T, Hohmann G W. 1993. A finite-difference time-domain solution for three-dimensional electromagnetic modeling[J]. Geophysics, 58(6):797-809.
[28] Yue J H, Yang H Y, Hu B. 2007. 3D finite difference time domain numerical simulation for TEM in-mine[J]. Progress in Geophysics(in Chinese), 22(6):1904-1909, doi:10.3969/j.issn.1004-2903.2007.06.036.
[29] 李建慧, 朱自强, 曾思红,等. 2012. 瞬变电磁法正演计算进展. 地球物理学进展, 27(4):1393-1400, doi:10.6038/j.issn.1004-2903.2012.04.013.
[30] 李展辉, 黄清华, 王彦宾. 2009. 三维错格时域伪谱法在频散介质井中雷达模拟中的应用[J]. 地球物理学报, 52(7):1915-1922, doi:10.3969/j.issn.0001-5733.2009.07.027.
[31] 刘鲁波, 陈晓非, 王彦宾. 2007. 切比雪夫伪谱法模拟地震波场[J]. 西北地震学报, 29(1):18-25.
[32] 龙桂华, 李小凡, 张美根. 2009. 错格傅里叶伪谱微分算子在波场模拟中的应用[J]. 地球物理学报, 52(1):193-199.
[33] 纳比吉安 米萨克 N. 1992. 勘查地球物理电磁法(第一卷)[M]. 赵经祥, 王艳君译. 北京:地质出版社, 302-348.
[34] 牛之琏. 2007. 时间域电磁法原理[M]. 长沙:中南大学出版社, 1-200.
[35] 岳建华, 杨海燕, 胡搏. 2007. 矿井瞬变电磁法三维时域有限差分数值模拟[J]. 地球物理学进展, 22(6):1904-1909, doi:10.3969/j.issn.1004-2903.2007.06.036.