正演是反演的基础，它的精度直接决定了反演解释的质量.在电(磁)法勘探的正演计算中，经常遇到贝塞尔函数的积分问题，而常规的积分方法计算比较复杂.19世纪70年代，荷兰学者把数字线性滤波方法引入地球物理学领域中(Ghosh，1971)，使得贝塞尔函数积分的计算变得简单且快速.基于数字线性滤波的快速汉克尔变换，在国内外学者的努力下，得到了快速的发展，出现了很多针对电测深的短滤波系数(Koefoed et al., 1972; Das and Ghosh， 1974).随着计算机的发展以及在正演计算中对精度要求的不断提高，采用短滤波系数远远不能满足高精度的计算要求，出现了多种长滤波系数求取方法：Johanson(Johansen and Sorensen， 1979)研究了汉克尔滤波系数的解析计算方法；Anderson(Anderson，1982)研究了用傅里叶变换求取滤波系数的方法，该方法对计算机内存要求较高；Guptasarma(Guptasarma and Singh， 1997)研究了计算滤波系数的维纳-霍普法，该方法的缺点是对优化方法要求比较高；Kong(Kong，2007)提出采用直接法求取的滤波系数，在海底电磁正演计算中取得了较好的效果.长滤波系数的采用，满足了计算精度的要求，但是需要消耗大量的计算时间，这也成为选择系数时的考量.

现有的汉克尔长滤波系数，大部分能满足高精度计算的要求，相对误差可控制在10^－5以上，甚至达到了10^－10以上.在地球物理的实际应用时，需要权衡计算精度和计算耗时(Anderson，1979; Chave，1983; Anderson，1989).现有的汉克尔滤波系数的选择很多，每种选择的计算精度与计算耗时各不相同.近年来，虽然汉克尔变换的理论研究进展缓慢，但是汉克尔线性滤波在实际计算时依然应用广泛(岳安平等，2009；李永兴等，2010；李帝铨等，2011； Key，2012；肖加奇等，2013).Guptasarma等(Guptasarma and Singh， 1997)对几组高精度滤波系数的选择进行了精度对比评价，但目前还没有学者对汉克尔滤波在具体勘探方法的正演计算中进行系统的研究.本文在Guptasarma研究的基础之上，分析了几种高精度长滤波系数的精度，并应用于电阻率测深法、频率域电磁测深法、瞬变电磁测深法中，对比分析了它们的计算精度和运算时间，对汉克尔滤波系数的选择起到了指导作用，对节省计算机资源有着相当的意义.

1 快速汉克尔正演理论

在地球物理场的模拟中，正演公式往往包含零阶或一阶贝塞尔函数的积分，也就是一种汉克尔变换.如果把连 续函数g(r)在(0，∞)上表示为

当然，采样点的位置会因为求取滤波系数的方法、或者说定义方法的不同而不同，但最终快速汉克尔变换，都能表示能这样一种线性形式(朴化荣，1990)，且滤波系数成指数衰减.由于汉克尔滤波算法计算贝塞尔函数积分，较用传统的方法计算速度至少快一个数量级，所以称为快速汉克尔算法.

2 经典滤波系数精度对比

根据近年来各滤波系数被引用的频率，如表 1所示，选择了五组不同长度的高精度长滤波系数，其中第4组和第5组属同一种算法求得的系数.如式(3)~(6)所示是四个有解析解的汉克尔积分式，用来对比不同系数计算时的精度.其中(3)和(4)式为零阶积分式，(5)和(6)式为一阶积分式.图 1、图 2为不同系数计算零阶变换的精度对比，图 3、图 4为计算一阶变换的精度对比.

表 1(Table 1)

表 1 五种高精度长滤波系数 Table 1 The five high-precision long filter coefficients 序号 作者 零阶系数长度 一阶系数长度 1 Anderson 283 283 2 Anderson 801 801 3 Guptasarma and Singh 61 47 4 Guptasarma and Singh 120 140 5 Fong 201 201

表 1 五种高精度长滤波系数 Table 1 The five high-precision long filter coefficients

如图 1可知，对于该汉克尔积分式(3)，在r<10^－2时，801点的滤波系数精度最高，120点的精度次之；在r>10^－2时，201点系数的精度最高，801点系数精度次之.61点滤波系数精度全区间最低.五种滤波系数在r∈[10^－4，10⁰]时，精度最高；而当r>1时，精度都急剧下降.由图 2可知，对于该汉克尔积分式(4)，在r>10^－2时，201点系数精度最高；在r<10^－2时，120点系数精度最高.在r>10^－4时，283点滤波系数精度最低；在r<10^－6时，61点和201点滤波系数精度最低.

图 1 五组系数计算变换式3的精度对比Fig. 1 The relative errors of the five J0 filters for(3)

图 2 五组系数计算变换式4的精度对比Fig. 2 The relative errors of the five J0 filters for(4)

由图 3可知，对于汉克尔积分式(5)，801点的滤波系数精度最高，在r<1时精度平稳；140点的精度次之，而47点的滤波系数精度最低.所有滤波系数，在r∈[10^－4，10⁰]时，精度较高；而当r>1时，精度都急剧下降.由图 4可知，对于汉克尔积分式(6)，在r较小的时候，801点的滤波系数精度最高；在r>10^－2时，201点的精度高于140点和801点系数的精度.47点滤波系数在r<10^－4时精度最低，而283点的滤波系数在r>1时精度最低.

图 3 五组系数计算变换式5的精度对比Fig. 3 The relative errors of the five J1 filters for(5)

图 4 五组系数计算变换式6的精度对比Fig. 4 The relative errors of the five J1 filters for(6)

以响应和视电阻率为研究对象，分析五种滤波系数在三种电(磁)法勘探中的正演计算精度及计算时间.模拟时无法得到解析解，就对比不同系数计算的曲线的形态特征.对于运算时间的对比，都是在个人电脑相同的软硬件条件下进行的.

3.1 在电阻率测深中的应用与对比分析

采用对称四极装置，求取水平层状大地的视电阻率(陈仲候和何昌礼，1980).电阻率转换函数T(λ)的形式是众所周知的，对于任何层状模型，当λ→∞时T(λ)→1.因此，汉克尔积分核函数λT(λ)随着λ的增大而增大.但是283点滤波系数计算时，使用的核函数必须随着参数的增加是连续边界的复函数.其他的滤波系数不需要这个条件，因为设计时的形式不同.所以使用283点系数时，做以下变换.由于下式的关系为

模型1： 二层层状介质D型模型，各层参数为ρ₁=400 Ωm，ρ₂=20 Ωm，h₁=20 m.图 5为模拟得到的视电阻率曲线图.

图 5 二层D型模型电阻率测深视电阻率曲线图Fig. 5 Apparent resistivity sounding curve of two-layer D-type model

模型2： 三层层状介质K型模型，层参数为ρ₁=60 Ωm，ρ₂=400 Ωm，ρ₃=60 Ωm，h₁=20 m，h₂=50 m.图 6为模拟得到的视电阻率曲线图.

图 6 三层K型模型直流电测深视电阻率曲线图Fig. 6 Apparent resistivity sounding curve of three-layer K-type model

模型3： 两层层状介质D型模型，层参数较极端，两个视电阻率差别较大，ρ₁=1 Ωm，ρ₂=0.000025 Ωm，h₁=1 m.图 7为模拟得到的视电阻率曲线图.

图 7 二层D型极端模型直流电测深视电阻率曲线图Fig. 7 Apparent resistivity sounding curve of two-layer D-type model extremely

由电阻率测深的模型1模型2的模拟结果可知，对于一般的模型进行正演时，五种滤波系数基本都能达到精度要求，正演曲线完全重合.由模型3的模拟结果可知，对于比较极端的模型，283点滤波系数在曲线尾部计算不准确，而201点滤波系数计算结果偏小.

3.2 在频率域电磁测深中的应用与对比分析

选择电偶源发射装置，发射电流1 A，收发距为5000 m，电偶极子长为400 m.计算电磁场响应、卡尼亚电阻率、以及用E_x，H_y，Re(E_x/H_z)分量计算的全区视电阻率(汤井田和何继善，1994；刘颖等，2011).

模型4： 三层层状介质H型模型，各层参数为ρ₁=200 Ωm，ρ₂=50 Ωm，ρ₃=200 Ωm，h₁=80 m，h₂=200 m.图 8为模拟得到的电磁场响应的曲线图，图 9为模拟得到视电阻率曲线图.

图 8 三层层状介质H型模型电磁场响应的绝对值(a)Ex曲线；(b)Hy曲线.Fig. 8 The electromagnetic field curve of three-layer H-type model (a)The curve of Ex；(b)The curve of Hy.

由图 8可知，五种滤波系数计算的模型4的电磁响应曲线图完全重合.而由图 9可知，模型4的卡尼亚视电阻率曲线完全重合，而不同分量计算的全区视电阻率有差别，主要表现在：对于E_x计算的全区视电阻率，五种滤波系数计算的曲线基本重合；对于H_y计算的全区视电阻率曲线，120&140、201 点系数的计算结果能够较好的反映出曲线形态，801点系数计算的曲线在尾部不正确；而对于Re(E_x/H_z)计算的全区视电阻率曲线，中间部分差别较大，801点能较好的反映出真实的曲线形态.笔者用五种滤波系数计算了不同收发距、不同发射电流、不同模型的视电阻率曲线，结果基本都与上例模型一致.用不同滤波系数求取的响应与卡尼亚电阻率完全一致，而全区视电阻率有区别，这是因为响应与卡尼亚电阻率是由贝塞尔函数直接积分得到的量，而全区视电阻率是由响应，通过等效法，二分搜索得到，这算法本身就有误差，它可能使滤波系数计算的微小误差积累，也会产生新的误差.所以全区视电阻率的求取不仅与所选的滤波系数有关，还与计算所用响应分量有关.

图 9 三层H型模型频率域电磁测深视电阻率曲线图(a)卡尼亚视电阻率曲线图；(b)Ex分量计算的全区视电阻率； (c)Hy分量计算的全区视电阻率； (d)Re(Ex/Hz)分量计算的全区视电阻率. Fig. 9 frequency electromagnetic sounding apparent resistivity curve of three-layer H-type model(a)Cagniard apparent resistivity；(b)resistivity calculated by Ex；(c)Resistivity calculated by Hy； (d)Resistivity calculated by Re(Ex/Hz).

选择电偶源发射装置，下阶跃激发，发射电流为1 A，收发距为5000 m，电偶极子长为400 m.瞬变响应采用G-S变换法，用垂直磁场的导数计算早晚期视电阻率，用垂直磁场来计算全区视电阻率(殷长春和李吉松，1994；陈明生和田小波，1999).

模型5： 两层水平层状介质D型模型，层参数为 ρ₁=200 Ωm，ρ₂=20 Ωm，h₁=200 m.图 10为模拟得到的响应和视电阻率曲线图.

图 10 二层D型模型瞬变电磁测深的响应和视电阻率曲线图(a)垂直磁场导数曲线图；(b)早晚期定义的视电阻率曲线图；(c)全区定义的视电阻率曲线图.Fig. 10 transient electromagnetic sounding curve of two-layer D-type model(a)The curve of dhz/dt;(b)The curve of early and late; (c)The curve of whole area.

模型6： 三层水平层状介质H模型，其中ρ₁=200 Ωm，ρ₂=20 Ωm，ρ₃=200 Ωm，h₁=200 m，h₂=200 m.图 11为模拟得到的响应和视电阻率曲线图.

图 11 三层H型模型瞬变电磁测深的响应和视电阻率曲线图(a)垂直磁场导数曲线；(b)早晚期定义的视电阻率曲线图； (c)全区定义的视电阻率曲线图.Fig. 11 transient electromagnetic sounding curve of three-layer H-type model(a)The curve of dhz/dt;(b)The curve of early and late; (c)The curve of whole area.

由模型5模型6可知，对于垂直磁场导数和早晚期视电阻率而言，在关断时间t<10^-6较小时，283点滤波系数计算结果有振荡，47点计算结果偏大，其它时间区间五条曲线完全重合.对于全区视电阻率而言，在关断时间t<10^-5时，除140点系数计算的曲线外，其它曲线都偏大或者偏小，其它时间区间曲线完全重合.不过，实际生产时很难用到过小的关断时间，如果取t>10^-4，各曲线都能反映出正确的地电形态. 3.4 各模型计算耗时对比

各滤波系数计算所耗时间与积分的核函数有关，在地球物理中，不同模型的积分核函数不同.六个模型的五种滤波系数计算所耗时间对比如图 11所示.由图可知，耗时基本随着滤波系数个数的增加而增加，47&61点系数耗时最小，801点系数耗时最大，801点是47&61点系数的好几十倍.图 11a与c为同为电阻率测深的二层模型，只是层参数不同，但所耗时间基本一致，说明改变模型层参数对计算耗时的影响不大.

图 12 各模型计算所耗时间对比(a)模型1所耗时间；(b)模型2所耗时间；(c)模型3所耗时间；(d)模型4所耗时间；(e)模型5所耗时间；(f)模型6所耗时间.Fig. 12 The time consuming of the five filter coefficients for different models(a)The time consuming of model 1；(b)The time consuming of model 2；(c)The time consuming of model 3；(d)The time consuming of model 4；(e)The time consuming of model 5；(f)The time consuming of model 6.

4.1     每组滤波系数对于不同的积分核函数，有着不同的计算精度；对于同一积分内核，计算区间不同，计算精度不一样；不同滤波系数的精度并不是系数个数越多就越高，它与权值求取方法、采样间隔、系数长度等都有关系，当然，相同算法求得的系数，系数越长精度越高.

4.2    五种滤波系数的在电法勘探中的不同表现，主要体现在：在电阻率测深中，对于极端模型，283点滤波系数在尾部计算结果不准确，201点系数计算结果偏小，一般情况下五种系数都能满足计算要求；在频率域电磁测深中，各响应计算的曲线完全重合，卡尼亚电阻率曲线也完全重合，但各分量计算的全区视电阻率有区别；在瞬变电磁测深中，对于垂直磁场的导数和早晚期视电阻率，283点系数稍有振荡，47点系数结果在关断时间很小时偏大，而对于全区视电阻率除140点系数外，各系数计算的曲线在关断时间很小时都偏大或者偏小.所有模型计算的耗时，都是与系数个数成正比，801点系数计算耗时最大.

4.3    综合来看，283点系数，必须在核函数有连续边界时才能使用，增加了计算的复杂性；201点系数在某些计算区间计算精度较其它系数要高；140&120点系数一般情况下计算精度都较高，但其受求值区间影响较大；801点系数计算精度较高，在不同的求值区间精度也较稳定，但其耗时太大.在具体计算时，需权衡计算精度与计算效率，本文的研究对滤波系数的选取有一定的指导作用.

致 谢    本文在修改过程中，得到中南大学地球科学与信息物理学院郭荣文老师和刘海飞老师的大力帮助，在此表示感谢.