航空时间域电磁法回线源有限差分初始场计算

许洋铖 林 君 嵇艳鞠 巩源峰 关珊珊

(吉林大学仪器科学与电气工程学院/地球信息探测仪器教育部重点实验室,吉林 长春130025)

再根据垂直磁场和水平磁场的关系[15],求出磁场水平分量。图1给出了采用回线积分法和纳比积安中心回线解析法计算的圆形线圈均匀半空间模型,回线中心点垂直二次磁场响应曲线。为了对比各种方法,文中各种方法均采用均匀半空间模型,同样的计算参数,如表1所示。

由图1中可见,回线积分法和解析表达式计算结果一致,误差精度能达到0.000 1%,计算精度完全满足有限差分初始场的要求。

1.引 言

直升机时间域电磁系统(airborne time domain electromagnetic)具有效率高、成本低、分辨率高、操控灵活等优点,被广泛应用于矿产资源、地质调查和环境监测等领域。世界上发达国家已经将直升机航空电磁测量系统用于商业化服务[1],我国立项开展吊舱式时间域直升机航空电磁勘查系统的研发。

目前,时间域航空电磁野外数据,多采用层状均匀大地模型进行反演解释,当在电阻率横向变化比较明显的区域,层状均匀大地不能反映真实的地质情况,只有二维、三维异常体模型更符合客观实际地质条件。时间域有限差分方法(a finite-difference approach for time-domain solution,简称 FDTD),比较适合模拟地下复杂形状的三维导电体以及电阻率剧烈变化的模型,具有计算简单、快速等特点。

对于时间域瞬变电磁响应的数值计算,1985年,HOHMANN和ADHIDJAJA[2]在YEE[3]提出交错网格的基础上,实现了二维导电体的显式有限差分瞬变电磁响应正演,1992年,LEPPIN M[4]利用隐格式有限差分,计算了方形线圈三维导电体的瞬变电磁响应,1993年,WANG T和HOHMANN G W[5]采用FDTD方法实现了回线源的三维瞬变电磁正演问题,近年来,魏文博、谭捍东[6]和宋维琪[7]等研究了回线源三维瞬变电磁有限差分正演计算,闫述[8]、李桐林[9]等也研究了二维长导线源有限差分正演,都取得了一定的成果。

采用有限差分法计算时间域电磁响应时,初始场值的计算和边界条件的选择直接决定了计算精度和速度。HOHMAN和ADHIDJAJA等人的计算都是将回线源初始场等效成有限个磁偶极子初始场的和。在吊舱式时间域航空电磁系统,采用圆形回线发射、中心装置,初始场计算含有双重贝塞尔函数,若采用等效成磁偶极子进行近似,会给迭代带来较大误差[10],而且对每个剖分网格都要计算,将导致计算量变得非常庞大,为此,本文研究了圆形回线源的初始场数值计算方法,提出了混合法计算初始电磁场。

2.初始场数值计算

有限差分法中,对于初始时刻的瞬变电磁场的计算,根据电磁场理论[11],在足够短的时间内,假设大地的上部分是均匀半空间,采用均匀半空间电磁响应的积分方程进行计算初始电磁场。选择初始时间[12]为

式中：μ1是大地表层磁导率;σ1是大地表层电导率;Δmin是最小网格步长。在地面及空气中圆回线的二次场电磁响应表达式为

式中 ：核函数 F0=rTEeu0(z-h)λ,u0=(λ2-ω2μ0ε0+iω μ0σ0)1/2,为空气中磁导率,I为发射电流,a 为发射回线半径,h为发射线圈距地面高度,z为接收线圈距地面距离,rTE为地表反射系数,J1为一阶贝塞尔函数,为零阶贝塞尔函数。

上述式(2)～(4)中都含有双重贝塞尔函数,无法用解析式表达出来。由于有限差分方程迭代的计算量大,为此,要求初始场值的计算方法不仅精度高,而且速度快。针对式(2)～(4)采用了回线积分法、宗量近似法、贝塞尔函数展开法计算初始电磁场值,利用汉克尔变换[13]、高斯积分和Guptasarma数字滤波方法[14],实现了时间域电磁响应计算。

2.1 回线积分法

回线积分法是将电偶极子产生的切向电场分量(不含接地项)沿发射回线进行积分,先得到发射回线的感应电压,再根据互易等效原理,计算接收线圈的感应电压,最后根据法拉第电磁感应定律,以及圆形线圈与方形回线面积等效,得到频率域二次磁场表达式,为

再根据垂直磁场和水平磁场的关系[15],求出磁场水平分量。图1给出了采用回线积分法和纳比积安中心回线解析法计算的圆形线圈均匀半空间模型,回线中心点垂直二次磁场响应曲线。为了对比各种方法,文中各种方法均采用均匀半空间模型,同样的计算参数,如表1所示。

由图1中可见,回线积分法和解析表达式计算结果一致,误差精度能达到0.000 1%,计算精度完全满足有限差分初始场的要求。

图2为采用回线积分法计算地面垂直二次磁场分布。计算时间为71 898.547 32 s。虽然沿回线积分的精度很高,可是它的计算速度太慢,不能满足FDTD对速度的要求。

图1 回线积分法中心点垂直磁场

表1 几种方法的计算参数

图2 回线积分法地面垂直磁场

由于发射回线在空间各处电磁场响应没有解析式(地面中心点除外),而回线积分法在中心点的计算值,同与其他方法相比,更接近理论值,所以下面在地面其他点将采用回线积分法作为比较标准。

2.2 宗量近似法

宗量近似法是将贝塞尔函数根据求解区域,进行大、小宗量近似。在回线外,当a≪r时,采用小宗量近似,有以垂直磁场为例,磁场响应表达式(2)写为

电磁响应表达式简化为含一重贝塞尔函数积分,计算变得简单。图3为采用小宗量近似法计算的地面各点垂直磁场响应。

图3 小宗量近似计算的地面垂直磁场

计算所用时间为0.323 043 s,从图3中可看出,小宗量近似只有在a≪r和a≫r时才与回线积分法计算的垂直磁场有很好近似,在a和r相差不大的情况下有较大误差,可达5%,此时不适合有限差分迭代计算。

采用大宗量近似时[16],利用Hankel函数的渐进特性,将磁场积分式分为两部分,即

由Hankel函数的大宗量近似性质和折线逼近算法有

利用贝塞尔函数的导数关系以及差商代替导数,有

式中：

此时将包含双重贝塞尔函数的无限积分转换为只含有一个贝塞尔函数积分的有限个有限区间的和,在每个有限区间上采用7点高斯积分。图4为采用大宗量近似计算的地面电磁响应分布。

图4 大宗量近似法地面垂直磁场

计算时间为223.977549 s。由图4中可以看出,大宗量近似与回线积分的误差精度能达到0.01%,但是由于有有限个区间的积分,虽然计算速度比回线积分法快,但是还是不适合计算大量的初始场,所以不适合FDTD。

2.3 贝塞尔函数展开法

贝塞尔函数展开法是将含有发射半径的贝塞尔函数J0(λ a)项写入核函数中,核函数写为：F1=rTE eu0(z-h)λ J1(λ a)。将 J1(λ a)进行级数展开 ,以垂直磁场为例：

将贝塞尔函数展式(17)代入表达式(16)中,直接计算即可得到电磁场响应。图5为贝塞尔函数展开法计算的电磁场响应分布。

图5 贝塞尔函数展开法地面垂直磁场

计算时间为1.941 0311 s,从图5中可看出,用贝塞尔函数展开直接计算的垂直磁场在地面各点的响应曲线,在r≤30 a时同回线积分法计算误差精度能达到0.01%,但是在r＞30 a时开始出现震荡,这是由于核函数中引入了贝塞尔函数的缘故,所以只能计算r≤30 a的初始电磁场。

2.4 几种算法对比分析

将回线积分法、宗量近似法、贝塞尔函数展开法进行对比、分析,各种方法的计算精度、速度、适用区间,如表2所示。

表2 几种方法比较

回线积分法精度较高,可以计算全空间的电磁场,但是计算速度慢,不适合大规模计算电磁场。

小宗量近似法计算速度较快,但是只能是在极端的情况下,精度才能满足要求;在收发距与圆线圈半径相差不大的情况下,误差比较大,不能计算全部空间的电磁场;大宗量近似虽然精度高,可以计算全部空间的电磁场,但是也受到计算速度的限制。

贝塞尔函数展开法计算较为简单,计算时间快,但是只能在r≤30 a时精度满足要求,在r＞30 a时出现震荡。

3.混合法计算初始电磁场

采用有限差分计算时,需要计算大地表层的初始电磁场,网格划分越细,计算精度越高,计算量也越大,为此,提出了采用混合法进行计算初始电磁场,在r≤30 a采用贝塞尔函数展开法、在r＞30 a时采用小宗量近似法计算FDTD的初始场。

图6为采用混合法计算均匀大地表面的时电场和时磁场分布。由图中可见,初始场基本没有截断误差。

4.结 论

在计算有限差分法初始电磁场时,需要计算大量的电磁场,对计算精度和速度都有要求。通过对回线积分法、宗量近似法、贝塞尔函数展开法三种方法计算初始场进行对比,研究表明,为了计算全空间的电磁场初值,采用任何单一的方法,都存在缺陷,只有采用混合法,将贝塞尔函数展开与小宗量近似法结合,不仅计算速度快,而且精度高,满足FDTD对计算时间和计算精度的要求,适合FDTD的初始场计算。

[1]雷 栋,胡祥云,张素芳.航空电磁法的发展现状[J].地质找矿论丛,2006,21(1)：40-53.LEI Dong,HU Xiangyun,ZHANG Sufang.Development status quo of airborne electromagnetic[J].Contributions to Geology and Mineral Resources Research,2006,21(1)：40-53.(in Chinese)

[2]ADHIDJIAJA J I,HOHM ANN G W,ORISTAGLIO M L.Two-dimensional transient electrom-agnetic responses[J].Geophsics,1985,50(12)：2849-2861.

[3]YEE K S.Numerical solution of initial boundary value problems involvingMaxwell equations in isotropic media[J].IEEE Trans.Antennas propa-gat,May 1966,AP-14(3)：302-307.

[4]LEPPIN M.Electromagnetic modeling of 3-D Sources over 2-D inhomogeneities in the time do-Main[J].Geophsics,1992,57(8)：994-1003.

[5]WANG T and HOHMANN G W.A finite-difference,time-domain solution for three-dimensiona lelectromagnetic modeling[J].Geophsics,1993,58(6)：797-809.

[6]谭捍东,余钦范,John Booker,等.大地电磁法三维交错采样有限差分数值模拟[J].地球物理学报,2003,46(5)：705-711.TAN Handong,YU Qinfan,John Booker,et al.Magnetotelluric three-dimensional modeling using the staggered-grid finite difference method[J].Chinese Journal of Geophysics,2003,46(5)：705-711.(in Chinese)

[7]宋维琪,仝兆歧.3D瞬变电磁场的有限差分正演计算[J].石油地球物理勘探,2000,35(6)：751-757.SONG Weiqi,TONG Zhaoqi.Forward finite differential calculation for 3-D transient electromagnetic field[J].Oil Geophysical Prospecting,2000,35(6)：751-757.(in Chinese)

[8]闫 述,陈明生,傅君眉.瞬变电磁场的直接时域数值分析[J].地球物理学报,2002,45(2)：275-283.YAN Shu,CHEN Mingsheng,FU Junmei.Direct time-domain numerical analysis of transient electromagnetic fields[J].Chinese Journal of Geophysics,2002,45(2)：275-283.(in Chinese)

[9]徐凯军,李桐林.时域瞬变场电磁场有限差分法[J].世界地质,2004,23(3)：301-305.XU KAijun,LI Tonglin.A finite-difference time-domain solution for transient electromagnetic fields[J].World Geology,2004,23(3)：301-305.(in Chinese)

[10]朱凯光,林 君,程德福.基于MAT LAB的高频谐变电磁场的数值计算[J].电波科学学报,2002,17(3)：224-228.ZHU Kaiguang,LIN Jun,Cheng Defu.Numerical computation of high frequency electromagnetic fields under MA TLAB[J].Chinese Journal of Radio Science,2002,17(3)：224-228.(in Chinese)

[11]NABIGHIAN M N.勘查地球物理-电磁法[M].北京：地质出版社,1992：204-207.

[12]ORISTASGLIO M L and HOHMANN W.Diffusion of electromagnetic field into a two-dimentional earth：a finite-difference approach[J].Geophsics,1984,49(7)：870-894.

[13]GUPT ASARMA D,SINGH B.New digital linear filters for Hankel J0 and J1 transforms[J].Geophysical Prospecting,1997,45(5)：745-762.

[14]GUPTASA RM A D.Computation of the time-domain response of a polarizable ground[J].Geophysics,1982,47(11)：953-963.

[15]NABIGHIAN M N.Toward a three-dimensional automatic interpreta-tion of potential field data generalized Hilbert transforms：fundamental relations[J].Geophsics,1984,49(6)：780-786.

[16]华 军,蒋延生,汪文秉.双重贝塞尔函数积分的数值计算[J].煤田地质与勘探,2001,29(3)：58-62.HUA Jun,JIANG Yansheng,WANG Wenbing.The numerical integration of dual hankel transformation[J].Coal Geology&Exploration,2001,29(3)：58-62.(in Chinese)