地球物理学进展  2014, Vol. 29 Issue (5): 2183-2190   PDF    
引用本文
胡瑞华, 林君, 李桐林, 刘长胜, 孙彩堂, 周逢道. 均匀大地交变电偶极子电磁响应直接计算式[J]. 地球物理学进展, 2014, 29(5): 2183-2190, doi: 10.6038/pg20140529  
HU Rui-hua, LIN Jun, LI Tong-lin, LIU Chang-sheng, SUN Cai-tang, ZHOU Feng-dao. Direct calculation of the response of an electric dipole on a homogeneous earth. Progress in Geophysics, 2014, 29(5): 2183-2190, doi: 10.6038/pg20140529   
均匀大地交变电偶极子电磁响应直接计算式
胡瑞华1,2 , 林君1,2, 李桐林3 , 刘长胜1,2, 孙彩堂1,2, 周逢道1,2    
1. 吉林大学仪器科学与电气工程学院, 长春 130061;
2. 吉林大学地球信息探测仪器教育部重点实验室, 长春 130061;
3. 吉林大学地球探测科学与技术学院, 长春 130061
收稿日期: 2013-11-25 修回日期: 2014-1-5 .
基金项目：油页岩勘探开发利用产学研用合作创新研究项目(OSR-02-01)和吉林大学学部青年基金(4305050102RV) 联合资助.
作者简介：胡瑞华,男,1974年生,宁夏中宁人,博士研究生,主要从事电磁法正反演与应用研究.(E-mail: hurh@jlu.edu.cn)
通讯作者：李桐林,男,1962年生,内蒙通辽人,教授,主要从事电磁法理论与应用教学和研究,特别是复电阻率正反演. (E-mail: lilaoshizh@163.com)
摘要：为了能计算出交变电偶极子在均匀大地中的电磁响应,本文在朴化荣编著的《电磁测深法原理》一书中的谢昆诺夫矢量位的表达式和矢量位散度的表达式基础上,用数学手段严格推导出矢量位的偏导数和矢量位散度的偏导数,从而得到不含汉克尔变换的场计算式;由于场计算式中含有贝赛尔函数,为计算方便,选用Matlab编程并调用其提供的贝赛尔函数,然后用场直接计算式求解了在源参数和场点参数已知情况下的8个不同深度处的场分量并以表格的形式给出了电磁分量与总场的计算值.通过总结得出3点结论:(1)本文方法可计算除源点外任一点场,解决了数字滤波法无法计算源点正下方场的问题;(2)本文方法计算场的时间小于1毫秒,计算速度快;(3)本文场计算式有严格的数学推导,其计算结果可以被三维异常体的正反演直接使用.因此,本文的场计算方法具有一定的应用前景.
关键词： 电偶极子     矢量位     电磁场     计算式     Matlab    
Direct calculation of the response of an electric dipole on a homogeneous earth
HU Rui-hua1,2, LIN Jun1,2, LI Tong-lin3 , LIU Chang-sheng1,2, SUN Cai-tang1,2, ZHOU Feng-dao1,2    
1. College of Instrument Science & Electrical Engineering, Jilin University, Changchun 130061, China;
2. Key Lab of Geo-Exploration and Instrumentation, Ministry of Education, Jilin University, Changchun 130061, China;
3. College of Geoexploration Science and Technology, Jilin University, Changchun 130061, China
Abstract: To compute an AC electric dipole's response on a homogeneous earth, the paper uses rigorous mathematical methods deducing the S.A.Schelkunoff vector potential partial derivatives and its divergence's ones based on the vector potential and its divergence's expressions from the book named “The Principle of Electromagnetic Sounding” written by Piao Huarong, through which the fields computation expressions are gained not including Hankel transforms. The expressions, however, still are made up Bessel functions. For convenience of calculation, the Matlab programming platform is chosen to program, in which the Bessel functions supplied by Matlab are called.And then the paper has computed fields components of eight different depth on the conditions known of source parameters and fields positions through the fields direct computation expressions. The computation results are displayed in form of sheets at last. Through summarizing, the paper concludes three points as the following: Firstly, the computation method supplied by the paper can compute any point's fields except the source one, which has solved the problem of the fields' computation from the source right below that the digital filter does not deal with. Secondly, the consume of calculation time is less than 1 millisecond using the paper's method, which declares it has compute rapidly. Lastly, because of strict mathematical deduction, the fields' computation expressions' results can be used directly by the forward and the inversion on 3D anomalous body. As a result, the fields' computational method from the papers has certain application prospects.
Key words: electric dipole     vector potential     electromagnetic     calculation expressions     matlab    
0 引 言

交变电偶极子源是电法勘探中一类十分重要的激发源.野外探测时，通常将长导线与两个接地电极(常用A，B字母表示两个电极的名称)相连接，在发射系统的控制下经长导线向地下供一定强度的谐变电流，从而实现人工偶极源频域的激励模式.当接地电极AB的长度小于其中心点到观测点(也称场点)之间距离的5~10倍时，观测点处的场可近似认为是偶极子场.根据场的叠加原理，接地长导线场源可以看成偶极场源的组合，因此长导线场源的电磁场求解问题实质上以偶极场源电磁场的积分方式来解决(Nabighian，1992；朴化荣，1990).实际计算时，常把长导线源切分为多个电偶极子源，然后把每个电偶极子场叠加从而实现接地长导线源的电磁场计算.在地球电磁探测应用中，交变电偶极子场源在全空间、半空间以及分层介质模型中激发的正常电磁场(也被称为一次场)的研究和计算方法成为电法勘探者们必须面对和首要解决的问题.

众所周知，任何电磁现象的规律都遵从麦克斯韦微分方程组.从数学的角度看，计算介质电磁响应的本质就是在分界面条件和物理条件(如离源无限远处场衰减为零)的约束下求解微分方程组.由于电磁场本身是三维场，通过直接求解麦克斯韦微分方程组获取场解存在诸多的数学难题.为此，谢昆诺夫定义和引入了一组矢量势(位)和标量势(位)函数，为在由多个均匀段组成的空间中求解波动方程提供了便利，从而间接有效地解决了电磁场难于求解的数学难题(Nabighian，1992).地球电法勘探领域里的学者们应用谢昆诺夫势函数及其理论正逐步解决着复杂的地学问题(Hohmann，1975；Wannamaker et al., 1984a；Zhang et al., 2006；张辉等，2006；Rob，2007； Weidelt，2007; 贾定宇等，2013).根据频域地球探测中激励信号的频率特点，谢昆诺夫位函数满足齐次或非齐次的亥姆霍兹方程(Nabighian，1992).实际证实，由于矢量位的表达式中含有无穷大积分上限的汉克尔变换，且被积函数中还含有震荡的贝赛尔函数，因此很难直接计算矢量位.为了得到矢量位的数值解，一些学者采用常规的计算方法——数字滤波方法实现汉克尔变换的快速计算(Johansen et al., 1979；Guptasarma，1982；Chave，1983；Anderson，1989；Christensen，1990； Anderson，1991；Guptasarma et al., 1997； Kong，2007).由于汉克尔变换中核函数的差异，业界并没有统一的满足领域内所有应用的滤波系数，这导致采用数字滤波方法计算的场结果并不完全一致，从而使解释结果难以完全被人信服.此外，特别是在基于一次场的三维电磁探测的正反演中，三维异常体的响应必须通过已知的一次场来求解(Ting et al., 1981；Wannamaker et al., 1984b；Wannamaker，1991；Xiong，1992；Pankratov et al., 1995；Zhdanov and Fang， 1996；Pankratov et al., 1997；鲍光淑等，1999；陈桂波等，2010；范翠松，2012)，因此一次场计算的精确程度就显得格外重要.为了能绕开用数字滤波计算场的方法，作者寻找到了并推导出均匀大地交变电偶极子矢量位及相关的直接计算表达式，为一次场的计算提供了新算法.

本文首先直接给出基于谢昆诺夫势的均匀大地交变电偶极子的电磁场计算式，然后在已知矢量位计算式的基础上，推导出场表达式所需的矢量位的偏导数和矢量位散度的偏导数，接着用Matlab编程计算指定参数的场，最后对全文进行总结.

1 基于谢昆诺夫势的均匀大地交变电偶极子的电磁场计算式

将交变电偶极子布置在均匀大地表面，界面以上为空气，空气以下为大地，见图 1.图中采用右手笛卡尔坐标系，电偶源AB与x轴一致，坐标原点0与AB中点重合.AB的长度为dl；均匀大地的电导率为σ1；负谐变激励电流I=I0e－i ω t；k1为电磁波在下半空间的传播常数；M(x，y，z)为任一场点.

图 1 电偶极子在半空间分界面上 Fig. 1 Electric dipole on the half-space surface

根据文献(朴化荣，1990)中谢昆诺夫势理论，对电偶源而言，当使用 A 1表达下半空间矢量位时，则可用

计算下半空间的电磁场，公式中的 E是电场矢量，H 是磁场矢量；k21=-iωμ0σ1，i为虚数单位，ω为激励信号角频率，μ0是真空磁导率； ∇ 、 ∇ ·、 ∇ ×分别是梯度、散度和旋度算子.由于水平电偶源不存在与偶极子方向(这里指x方向)垂直的矢量位水平分量，故可以用 A1= A x1，0，A z1 分量形式表达 A 1，且

式中，J0和J1为第一类实变量Bessel函数.这样，由式(1)分解出电磁场分量形式

公式中 E x、 E y、 E z和H x、 H y、 H z分别是电场 E和磁场H 的三个分量；为偏导数算子.

2 矢量位的偏导数和矢量位散度的偏导数计算

由(3)式可知，若能找到比(2)式更易计算的 A x1，A z1，∇ · A1的表达式以及A x1对y，z的偏导数，A z1对x，y的偏导数，∇ · A 1对x，y，z的偏导数，便能计算出地下任一点的电磁场.文献(朴化荣，1990)(86页，87页)借助索末菲和福斯特积分已经推导出了与(2)式等价的 A x1，A z1及 ∇ · A 1的表达式，但由于其只关注了地表场的计算，没有继续推导地下场的计算公式，也没有给出相应的偏导数表达式.为使文献更加完善，本文的重点是推导出相应的偏导数.

2.1 矢量位偏导数和矢量位散度偏导数的分解

通过观察 A x1，A z1，∇ · A 1的表达式，将它们按Bessel函数的交叉项分解为若干个子项，便于书写及供同行参考和验证.将 A x1分解为5个子项的和，记为

其中，

将 A z1分解为5个子项的和，记为

其中，

将 ∇ · A 1分解为4个子项的和，记为

其中，

以上式中的x，y，z为场点的3个坐标分量，R为收发距(即图 1中OM的长度)且，I0和I1为第一类虚宗量(宗量取为)贝赛尔函数，K0和K1为第二类虚宗量(宗量取为 )贝赛尔函数，偶极距 P E=Idl.进一步将Ax1对y，z的偏导数，记为

A z1对x，y的偏导数，记为

∇ · A 1对x，y，z的偏导数，记为

2.2 矢量位的偏导数和矢量位散度的偏导数的推导

(1)辅助公式

为了推导 A x1，A z1，∇ · A 1对坐标的偏导数，需要以下的辅助公式(何淑芷和陈启流， 1997，叶其孝等，2011)，其中：

①收发距R对x，y，z的偏导数

②虚宗量 k1 2(R±z)对x，y，z的偏导数

结合虚宗量Bessel函数的微分性质(朴化荣，1990；Nabighian，1992)

③I0对x，y，z的偏导数

④I1对x，y，z的偏导数

⑤K0对x，y，z的偏导数

⑥K1对x，y，z的偏导数

这样就可以进一步推导出 A x1，A z1，∇ · A 1各子项对x，y，z的偏导数(共计32个)，通过适当的化简合并，得到最终的表达式如下.

(2)A x1的5个子项分别对y，z的偏导数

① A x11对y，z的偏导数

② A x12对y，z的偏导数

③ A x13对y，z的偏导数

④ A x14第四项对y，z的偏导数

⑤ A x15对y，z的偏导数

(3)A z1的5个子项分别对x，y的偏导数

① A z11对x，y的偏导数

② A z12对x，y的偏导数

③ A z13对x，y的偏导数

④ A z14对x，y的偏导数

⑤ A z15对x，y的偏导数

(4)∇ × A 1的4个子项分别对x，y，z的偏导数

①div A 11对x，y，z的偏导数

②div A 12对x，y，z的偏导数

③div A 13对x，y，z的偏导数

④div A 14对x，y，z的偏导数

这样，首先将 A x1的各子项分别对y和z的偏导数(式11a~15b)代入(3)式中，便可得到 A x1分别对y和z完整的偏导数，同理得到 A z1分别对x和y完整的偏导数，∇ × A 1对x，y，z完整的偏导数，然后将它们代入(3)式中，便得到计算场的表达式.显然，这样的表达式已不含有如(2)式中的无穷大积分限的积分符号.另外，本文的场计算公式仅无法计算源点处(R=0)的场，而(2)式除了无法计算源点处的场外还无法计算在地表投影为源点的地下场点(r=0)，因此本文的场计算公式扩大了场点的计算范围.

3 场计算实例

虽然在第2小节中已经推导出场的计算式，但由于计算式中含有I0，I1，K0，K1贝赛尔函数，因此能否直接有效地求解这些贝赛尔函数就成为直接计算场的关键，文献(朴化荣，1990)中指出对不同范围的虚宗量可利用贝赛尔函数的解析式或渐进展开式来计算.考虑到编程的便利，选取Matlab2012平台编制程序.Matlab软件是全球公认的优秀的数值计算平台，而且提供了直接可调用的贝赛尔函数(王正林，2009)，与本文中的贝赛尔函数对应关系见表 1，因此(2)式场的表达式是完全可计算的.

根据电偶极子的条件并参考图 1中的观测装置，选取以下参数：σ1=1/100(Ωm)-1，dl=1200 m，I=20 A，ω=2π*10 rad/s，x=12000 m，y=12000 m，磁导率μ0=4π*10－7H/m，对第2节中的场计算式用Matlab2012编程，并分别计算z取0 m，200 m，500 m，800 m，1000 m，1200m，1500m，2000m共计8个不同深度的场分量，计算结果见表 2，各深度处场的计算时间分别为0.533 ms，0.330 ms，0.304 ms，0.295 ms，0.286 ms，0.292 ms，0.291 ms，0.321 ms.表 2中的数据和计算时间是在一台PC机(3.29 GHz的CPU，4G内存，WindowsXP操作系统)上运行的结果.

表 1 文中的贝赛尔函数与Matlab中的贝赛尔函数对应关系 Table 1 Correspondence relationship between the paper’s and Matlab’s

表 2 不同z处的场分量 Table 2 Field components for different z

由表 2中的场分量，进一步计算出总场，计算公式为

“| |”为取模运算符号，见表 3.

表 3 由表 2计算的总场 Table 3 Total fields according to Table 2

从表 2中各点场分量的计算值可以看出：除地表处的 E y分量为实数外，其余分量都为复数；对比电场分量发现，E z远小于 E x和E y；从计算表 2各点场的时间来看，花费时间大约在0.3~0.6毫秒之间；表 3反映了总场随着深度的增加而递减，符合场的传播规律.值得注意的是表 2和表 3的数据是与本文给定的参数相对应的计算结果，若参数发生变化则需重新计算.

4 结 论

本文在已有的均匀半空间交变电偶极源矢量位及矢量位散度表达式的基础上，通过数学手段推导出矢量位的偏导数和矢量位散度的偏导数公式，从而解决了均匀半空间交变电偶极源场的直接计算问题.为了方便场计算式中的贝赛尔函数的计算，采用先进的数值计算平台Matlab编程，并计算了指定参数对应的场分量.通过总结得出以下3点结论：

(1)采用直接计算式计算均匀半空间交变电偶极源的激发场的是对用数字滤波法计算场方法的补充和完善，因为本文的场计算式可以计算地下场点在地表投影为源点的区域；

(2)采用直接计算式计算均匀半空间交变电偶源的某点激发场，计算时间小于1毫秒，表明计算速度快；

(3)由于均匀半空间交变电偶极源的激发场的直接计算式有严格的数据推导，因此计算出的场值令人信服，可以被三维异常体的正反演直接使用.

因此，本文的场计算方法具有一定的应用前景.

致 谢 感谢贵刊编辑和专家对本文的评审及提出的修改意见; 感谢吉林大学公教中心数学教研室李立明副教授对本文公式的校对.

参考文献
[1] Anderson W L. 1989. A hybrid fast Hankel transform algorithm for electromagnetic modelling[J]. Geophysics, 54(2): 263-266.
[2] Anderson W L. 1991. Comment on 'Optimized fast Hankel transform filter' by Niels Boie Christensen[J]. Geophysical Prospecting, 39(3): 445-447.
[3] Bao G S, Zhang B X, Jing R Z, et al. 1999. Numerical modeling for 3D EM responses using integral equation solution [J]. Cent. South Univ. Technol (in Chinese), 30(5): 472-474.
[4] Chave A D. 1983. Numerical integration of related Hankel transforms by quadrature and continued fraction expansion[J]. Geophysics, 48(12): 1671-1686.
[5] Cheng G B, Wang H N, Yao J J, et al. 2010. Three-dimension modeling of frequency sounding in layered anisotropic earth using integral equation method[J]. Chinese Journal of Computational Physics. (in Chinese), 27(2): 274-280.
[6] Christensen N B. 1990. Optimized fast Hankel transform filters[J]. Geophysical Prospecting, 38(5): 545-558.
[7] Fan C S, Li T L, Yan J Y. 2012. Research and application experiment on 2. 5D SIP inversion[J].J. Geophys. (in Chinese), 55(12): 4044-4050, doi:10.6038/j.issn.0001-5733.2012.12.016.
[8] Guptasarma D. 1982. Optimisation of shorter digital filters for increasing accuracy[J]. Geophysical Prospecting, 30(4): 501-514.
[9] Guptasarma D, Singh B. 1997. New digital linear filters for Hankel J0 and J1 transforms[J]. Geophysical Prospecting, 45(5): 745-762.
[10] He S Z, Chen Q L. 1997. Math -Physical Method (in Chinese)[M]. Nanjing: Huanan Science and Technology University Press.
[11] Hohmann G W. 1975. Three dimensional induced polarization and electromagnetic modeling[J]. Geophysics, 40(2): 309-324.
[12] Jia D Y, Weng A H, Liu Y H, et al. 2013. Propagation of electromagnetic fields from a horizontal electrical dipole buried in ocean[J]. Progress in Geophys. (in Chinese), 28(1): 0507-0514, doi: 10. 6038/pg20130158.
[13] Johansen H K, Sorensen K I. 1979. Fast Hankel transforms[J]. Geophysical Prospecting, 27(4): 876-901.
[14] Kong F N. 2007. Hankel transform filters for dipole antenna radiation in a conductive medium[J]. Geophysical Prospecting, 55(1): 83-89.
[15] Nabighian M N. 1992. Electromagnetic Method in Applied Geohpysics[M]. Volume 1, Theory. Society of Exploration Geophysicists. Beijing: Geological Publishing House.
[16] Pankratov O V, Avdeev D B, Kuvshinov A V. 1995. Electromagnetic field scattering in a heterogeneous earth: a solution to the forward problem[J]. Physics Solid Earth, 31: 201-209.
[17] Pankratov O V, Kuvshinov A V, Avdeev D B. 1997. High performance three-dimensional electromagnetic modeling using a modeified neumann series: anisotropic case[J]. J. Geomag. Geoelectr, 49: 1541-1547.
[18] Piao H R. 1990. The Principle of Electromagnetic Sounding (in Chinese) [M]. Geological Publishing House,
[19] Rob L E. 2007. Using CSEM techniques to map the shallow section of seafloor: From the coastline to the edges of the continental slope[J]. Geophysics, 72(2): 105-116.
[20] Ting S C, Hohmann G W. 1981. Integral equation modeling of three dimensional magnetotelluric response[J]. Geophysics, 46(2): 182-197.
[21] Wang Z L, Gong C, He Q. 2009. Master MATLAB Science Computing (in Chinese)[M]. Beijing: Publishing House of Electronics Industry.
[22] Wannamaker P E, Hohmann G W, SanFilipo W A. 1984a. Electromagnetic modeling of three dimensional bodies in layered earth s using integral equations[J]. Geophysics, 49(1): 60-74.
[23] Wannamaker P E, Hohmann G W, Ward S H. 1984b. Magnetotelluric responses of three-dimensional bodies in layered earth using integral equations[J]. Geophysics, 49(9): 1517-1533.
[24] Wannamaker P E. 1991. Advances in three-dimensional magnetotelluric modeling using integral equations[J]. Geophysics, 56(11): 1716-1728.
[25] Weidelt P. 2007. Guided waves in marine CSEM[J]. Geophys, 171(1): 153-176.
[26] Xiong Z. 1992. electromagnetic modeling three-dimensional structures by the method of system iterations using integral equations[J]. Geophysics, 57(12): 1556-1561.
[27] Ye Q X, Shen Y H. 2011. Practical Mathematical Manual (in Chinese)[M]. Beijing: Science Press.
[28] Zhang H, Li T L, Dong R X. 2006. Modeling 3-D electromagnetic responses of the electric dipole using volume integral equation method[J].a href="http://www.progeophys.cn/CN/abstract/abstract2109.shtml" target=_blank> Progress in Geophysics (in Chinese), 21(2): 386-390.
[29] Zhdanov M S, Fang S. 1996. Quiasi-linear approximation in 3-D electromagnetic modeling[J]. Geophysics, 61(3): 646-665.
[30] 鲍光淑, 张碧星, 敬荣中,等. 1999. 三维电磁响应积分方程法数值模拟[J]. 中南工业大学学报, 30(5): 472-474.
[31] 陈桂波, 汪宏年, 姚敬金,等. 2010. 利用积分方程法的各向异性地层频率测深三维模拟[J]. 计算物理, 27(2): 274-280.
[32] 范翠松, 李桐林, 严加永. 2012. 2.5维复电阻率反演及其应用试验[J]. 地球物理学报, 55(12): 4044-4050, doi:10.6038/j.issn.0001-5733.2012.12.016.
[33] 何淑芷, 陈启流. 1997. 数学物理方法[M]. 南京: 华南理工大学出版社, 1997.
[34] 贾定宇, 翁爱华, 刘云鹤,等. 2013. 海洋环境中水平电偶极子电磁场特征分析[J]. 地球物理学进展, 28(1): 0507-0514, doi:10.6038/pg20130158.
[35] 朴化荣. 1990. 电磁测深法原理[M]. 北京: 地质出版.
[36] 王正林, 龚纯, 何倩. 2009. 精通MATLAB科学计算[M]. 北京: 电子工业出版社.
[37] 叶其孝, 沈永欢. 2011. 实用数学手册[M]. 北京: 科学出版社.
[38] 张辉, 李桐林, 董瑞霞. 2006. 体积分方程法模拟电偶源三维电磁响应[J]. 地球物理学进展, 21(2): 386-390.