戴世坤 曾 铃* 周印明③ 李 昆 陈轻蕊 凌嘉宣

(①中南大学地球科学与信息物理学院,湖南长沙 410083; ②中南大学有色金属成矿预测与地质环境监测教育部重点实验室，湖南长沙 410083； ③东方地球物理公司综合物化探处，河北涿州 072751)

0 引言

在电磁场地球物理方法中，均匀层状介质电磁场解析表达式含有汉克尔积分[1-2]，其核函数为0阶和1阶Bessel函数。随着其宗量的增大，Bessel函数呈现快速振荡和慢衰减的特征，造成汉克尔积分难以实现高效、高精度的计算，特别是高频和大收发距情况下难度更大。为了解决这个问题，学者们提出了许多方法，主要分为两大类： 第一类为数字滤波方法； 第二类为直接积分法(Direct integration method，DIM)。Ghosh[3]首先将数字滤波法引入地球物理数值计算； 针对电磁法，Johansen等[4]提出了汉克尔数字滤波系数的解析计算方法。之后，很多研究者给出了不同优化滤波系数[5-6]。蔡盛[7]对比了5组高精度快速汉克尔长滤波系数的计算精度，研究结果表明在一定的收发距和频率范围内，计算精度较高，但耗时太大，随着滤波系数的增大而增长。总体来说，数字滤波法适用于中、低频和适中的收发距情况，在高频和大收发距的情况下计算精度很低。考虑到数字滤波法的局限性，众多学者利用多项式精确拟合核函数的思想，达到求解解析解的目的，提出了直接积分类方法。代表性方法包括高斯积分法[8]、基于快速正弦余弦算法的汉克尔变换算法[9]、直接法[10]、离散复镜像法[11-12]、谱域法[13]、基于高阶窗函数的结合非线性变换的算法[14]、外推积分方法(QWE)[15]、三次样条插值法[16]等。其中最后一种方法采用三次样条函数对核函数进行插值，得到一系列核函数为多项式的简单Sommerfeld积分，再利用Bessel函数的递推公式和Lommel公式的渐进展开求解。这些方法虽能较好地应用于低频和适中收发距的均匀层状介质电磁场的计算，但随着频率和收发距的增大，这些方法的计算精度逐渐降低，甚至失效。针对高频电磁场的计算，郑圣谈等[17]提出了高密度采样数字滤波法，但当介质导电率变小、介电常数变大、收发距变大、频率增加时，该方法计算精度会大幅度降低，不具普适性。

针对电磁场汉克尔积分核函数复杂的特点，本文提出一种高效、高精度直接积分方法，用于频率和收发距广普适用的电磁场高效、高精度的数值计算。其基本思想是Bessel函数可以在两个区间分别用不同的多项式展开，每一区间的汉克尔积分可离散为多个单元积分之和，每个单元被积函数采用三次样条插值函数表示，由此可求得积分的解析解，并通过求和获得汉克尔积分的数值解。在此基础上，利用均匀全空间电偶极子电磁场的解析解，正确选取积分范围，并合理剖分积分单元。数值解与解析解对比表明，本文算法正确、可靠。与数字滤波类算法的比较表明，本文算法广泛适用于不同频率和不同收发距电磁场的计算，具有较强的普适性。

1 方法原理

均匀层状介质电磁场的解析解可用汉克尔积分表示

(1)

式中：g(m)为大地响应函数，m表示积分参数；r为发射源到接收点的水平距离，mr为宗量；Jn(·)为n阶Bessel函数。

式(1)中，Bessel函数为0阶和1阶，根据特殊函数手册[18]，可以分别在两个区间用不同的多项式展开。第一区间为0≤mr≤4，Bessel函数可展开为简单的代数多项式，每个积分单元内代数多项式与响应函数的乘积可用三次样条插值函数表示。第二区间为4

在第一个区间(0≤mr≤4)有

J0(mr)=-0.0005014415t7+0.0076771853t6-

0.0709253492t5+0.4443584263t4-

1.7777560599t3+3.9999973021t2-

3.9999998721t+1.0+e0(mr)

(2)

式中截断误差e0≤1×10-8。

在第二个区间(4

(3)

式中

(4)

其中截断误差e1、e2≤1×10-8。

对第一区间0≤mr≤4，将积分范围均匀剖分为N个单元。由式(2)可知，Bessel函数的展开式为简单的代数多项式，因此设f(m)=g(m)J0(mr)，这个区间的积分F1可表示为

(5)

式中mi(i=1,2,…,N,N+1)为离散采样点。每个单元内被积函数fi用三次样条插值函数[18-19]表示为

fi=di(m-mi)3+ci(m-mi)2+

bi(m-mi)+ai

(6)

式中ai、bi、ci、di为单元i的三次样条插值函数的系数。设Li=mi+1-mi为单元i的长度，单元i的积分解析表达式为

(7)

因此，第一区间各单元积分累加的结果为

(8)

对于第二区间4

(9)

设

(10)

(11)

(12)

式中

(13)

因此，第二区间各单元积分累加求和结果为

(14)

最终汉克尔积分为

F=F1+F2

(15)

2 积分参数的选取

为了提高计算精度和效率，需合理选取积分范围，恰当地进行单元剖分，而这些参数的选取依赖于响应函数的特征和变化规律[19-20]。本文以均匀全空间电偶极子源电磁场为例，通过研究不同响应函数的变化规律，确定积分范围，并选定单元剖分方案。

在均匀全空间中，将x方向上电偶极子源产生的电、磁场用汉克尔积分表示，其解析解的汉克尔积分表达式分别为

(16a)

(17)

(18)

取电导率σ=0.0001S/m，相对介电常数ε=1，电偶极子源中心坐标为(0,0,0)，偶极矩为1C·m，频率v=1×108Hz， 积分变量m在[1×10-12,1×108]区间按对数等间距均匀采样，每一个对数单位内采样点为5000个。计算上述三种响应函数的值，结果如图1所示。

图1 三种响应函数随m的变化曲线

由图1可以看出，m从1×10-6增大到1×10-1时，三种响应函数的对数值呈线性增大； 当m为1×10-1～1×101时，三种响应函数值随m变化较快，尤其在1×100附近出现尖脉冲；m>102时三种响应函数随着m的增大按照对数以近似线性规律迅速衰减； 在m接近1×103时，三种函数的值均衰减至低于10-20。根据上述分析结果，m的积分范围可选取为[0,1×104]，并整体上按照对数等间距规律进行单元剖分。特别地，在m为1×10-1～1×102时，因响应函数值变化较快，单元剖分可适当加密，尤其在m接近1×100时，响应函数变化剧烈，单元剖分需格外加密。

根据以上原则，对均匀全空间电偶极子源产生电磁场用汉克尔积分表达式进行数值计算，并与解析解进行比较，以验证算法的正确性。下文分别计算均匀全空间中1×104、1×108Hz的电场响应。

对v=1×104Hz，计算zr=40m所在平面的电磁场，计算范围：x方向为-500～500m，y方向为-500～500m，m的采样区间为0～1×104。在区间1×10-3～1×101加密积分单元，总采样点数为373。图2为电场分量的数值解、解析解及相对误差平面图。由图可见，数值解与解析解曲线形态基本一致，相对误差都很小，最大为8×10-6(Ex分量)，即便在源附近，其误差也不大于1%，表明本文算法的正确性和可靠性，且精度较高。

图2 v=1×104Hz时电场分量Ex(a)、Ey(b)、Ez(c)的解析解(左)、数值解(中)及其相对误差(右)

对v=1×108Hz，计算zr=10m所在平面的电场。计算的平面范围x方向为-15～15m，y方向为-15～15m，m的采样区间为0～1×105，在区间1×10-2～1×100加密采样点，总采样点数为811。图3为该模型的电场分量数值解、解析解及相对误差平面图。从图3可看出，数值解与解析解曲线形态基本一致，相对误差很小，其中Ey分量的误差最大，最大为4×10-6，即便如此，在源附近的误差也都不大于1%，这表明算法在高频时计算的场值是正确、可靠的，且精度较高。

图3 v=1×108Hz时电场分量Ex(a)、Ey(b)、Ez(c)的解析解(左)、数值解(中)及其相对误差(右)

上述两个频率的计算结果证明了本文算法正确、可靠，且能适应高频电磁场的计算。

3 模型算例

设计均匀全空间模型和均匀层状介质模型，分别使用本文DIM算法和常用的数字滤波类算法模拟电磁场，以验证本文算法对不同频率、不同收发距电磁场模拟计算的广普适用性。

3.1 均匀全空间模型

首先验证本文算法对于频率的广普适用性。

在均匀全空间中，设发射源到接收点的水平距离R=100m，计算频率范围为1×10-2～1×1010Hz，电导率σ=0.0001S/m，相对介电常数ε=1，电偶极子源中心坐标为(0,0,0)，偶极矩为1C·m。分别使用本文算法与四种经典的数字滤波算法计算该模型的数值解与解析解的相对误差，其中四种数字滤波算法的滤波系数选取见表1。数字滤波方法①和②参见文献[21]，方法③和方法④分别参见文献[22]和文献[23]。这四种方法及本文提出的直接积分法计算结果如图4所示，其中黑色直线表示相对误差为1%。从图4可看出，在频率v<1×106Hz的情况下，这四种数字滤波类算法的计算误差均小于1%；随着频率的升高，四种滤波系数类算法的误差越来越大，而本文直接积分法的计算误差小于1%。表明本文提出的直接积分法计算精度高，在广泛的频率范围内适用。

表1 四种常用的数字滤波算法的滤波系数长度取值

为验证该方法对收发距的广普适用性，在均匀全空间中，取频率v=1×107Hz， 发射源到接收点的水平距离r范围为1×10-3～1×103m，其他参数不变。使用本文算法与四种经典的数字滤波算法分别计算不同收发距时的数值解与解析解的相对误差，其中四种数字滤波算法的滤波系数选取同上(表1)。计算结果见图5。可以看出，在收发距r<2m时，直接积分法和四种数字滤波类算法的计算误差均小于1%； 但随着收发距的增大，四种滤波系数类算法的误差越来越大，而直接积分法在任意r时，计算误差均不大于1%。表明本文提出的直接积分法计算精度高，在广泛的收发距范围内适用。

图4 均匀全空间模型不同频率下直接积分法与四种

图5 均匀全空间模型不同收发距下直接积分法与四种

3.2 层状介质模型电磁场不同算法对比

设计一个三层层状介质模型，模型参数见图6。电偶极子源中心坐标为(0,0,0)，偶极矩为1C·m。

图6 三层层状模型断面示意图

计算频率v=1×108Hz、z=10m所在平面的电场分量，计算平面范围x方向和y方向均为-15～15m。

图7和图8分别为表1中第④种数字滤波法与本文直接积分法计算的Ex、Ey、Ez实部与虚部平面等值线图。由图可见，数字滤波法计算的Ex、Ey、Ez等值线圆滑度较低，表明计算结果有一定误差，而直接积分算法计算的电场等值线圆滑度较高，表明计算结果更准确、精度更高。

图7 三层层状介质模型直接积分法(上)与数字滤波法(下)计算的电场分量Ex(a)、Ey(b)和Ez(c)实部平面等值线图

图8 三层层状介质模型直接积分法(上)与数字滤波法(下)计算的电场分量Ex(a)、Ey(b)和Ez(c)虚部平面等值线图

4 结论

本文提出一种高效、高精度的直接积分方法，适用于不同频率和不同收发距的电磁场模拟计算，尤其适合高频电磁场的计算，具体取得以下结论。

(1)利用Bessel函数可以在两个积分区间分别用不同的多项式展开。将每一区间的汉克尔积分离散成多个单元积分之和，每个单元被积函数采用三次样条插值函数表示，由此可求得积分的解析解。通过叠加，求得汉克尔积分的数值解。

(2)响应函数的特征和变化规律是积分参数选取的关键因素。为此，以均匀全空间电偶极子源电磁场为例，研究了三种典型的响应函数随积分变量的变化规律。计算结果表明，三种响应函数随积分变量呈对数规律变化，据此按对数规律确定积分范围和单元剖分：在响应函数值变化快的地方，单元剖分可适当加密，尤其在尖脉冲附近，单元剖分需格外加密；在响应函数值变化慢的地方，单元剖分可适当稀疏。这样合理剖分积分单元可提高直接积分法的计算精度和效率。均匀全空间电偶极子源电磁场数值解与解析解对比表明本文算法更正确、可靠。

(3)分别设计均匀全空间模型和层状介质模型，对本文提出的直接积分法与数字滤波方法电磁场数值计算结果进行对比，结果表明本文积分法对频率和收发距具有广普适用性，尤其对高频电磁场，其优势更加突出，为高频电磁场的计算提供了重要的方法基础。

附录A 1阶Bessel函数汉克尔积分的解析表达式推导过程

根据特殊函数手册[18]，对J1(mr)函数分两个区间进行展开。

在第一个区间(0≤mr≤4)

0.0236616773t5+0.1777582922t4-

0.8888839648t3+2.6666660544t2-

3.9999999710t+1.9999999998)+e0(mr)

(A-1)

其中截断误差e0≤1×10-8。

在第二个区间(4

(A-2)

式中

(A-3)

0.000999941t-3+0.000266891t-2-

0.001601836t-1+0.093749994)+e2(mr)

(A-4)

其中截断误差e1、e2≤1×10-8。

对于第一区间(0≤mr≤4)，将其剖分为N个单元，由式(A-2)可知，Bessel函数的展开式为简单的代数多项式，因此设f(m)=g(m)J1(mr)，这部分的积分为

(A-5)

设mi(i=1,2,…,N,N+1)为离散采样点，每个单元内被积函数fi可采用三次样条插值函数表示为

fi=di(m-mi)3+ci(m-mi)2+

bi(m-mi)+ai

(A-6)

(A-7)

第一区间各单元积分累加表达式为

(A-8)

对第二区间(4

(A-9)

设

(A-10)

在每个单元内，Fp、Fq变化规律简单，可采用三次样条插值函数近似表示为

(A-11)

(A-12)

式中

(A-13)

可求得第二区间各单元的积分为

(A-14)

据此，汉克尔积分为

F=F1+F2

(A-15)