2024年5月30日发(作者：翁茂德)
第34卷第6期
2010年12月
物 探 与 化 探
GEOPHYSICAL&GEOCHEMICALEXPLORATION
Vo.l34,No.6
Dec.,2010
汉克尔变换的数值计算与精度的对比
张伟,王绪本,覃庆炎
(成都理工大学地球探测与信息技术教育部重点实验室,四川成都 610059)
摘要:汉克尔变换的数值解法是计算电磁测深正演理论曲线的有效工具。对基于数字滤波法的快速汉克尔变换
算法进行了公式推导,并用Guptasama和Singh给出的线性滤波系数进行了计算,对比分析了该数值算法与理论解
析式的误差分布特点。结果表明,该数值算法的计算结果连续一致逼近其理论值,且不存在振荡现象,计算精度
高,在数值模拟研究中具有较大的实用价值。
关键词:汉克尔变换;快速汉克尔变换;滤波系数
中图分类号:P631 文献标识码:A 文章编号:1000-8918(2010)06-0753-03
在水平层状介质地表的电磁测深数值计算中,
除平面波和偶极场源波区的电磁场可用具有递推性
质的频率特性函数R(X)表示,有限发-收距电磁测
深正演解析式均用含零阶或一阶贝塞尔函数的积分
形式表示出来,这些积分形式是汉克尔变换式,汉克
尔变换的数值解法是电磁测深数值正演计算的基础
问题,因此具有重要的研究意义。
基于数字滤波法的快速汉克尔变换具有编程简
单、运算速度快、执行效率高等优点,该方法广泛地
应用在电磁法(TEM、CASMT、FIP等)正演数值计算
中。笔者采用Guptasama和Singh给出的J
0
和J
1
变
换线性滤波系数进行计算
[2]
式(1),可得
]
G(v)=
-]
F(u)H
Q
M
(v-u)du=F*H
M
。(3)
G是函数F和H
M
的褶积,因此G可看成输入函数F
和线性滤波函数H
M
的系统响应。
在满足抽样定理的情况下,即$<1/|2f
c
|,$为
抽样间隔,f
c
为被抽样函数的最高频率。引入正弦
插值函数
sin[(P/$)(t-n$)]
,并定义函数P(x)=
(P/$)(t-n$)
sin(Px)/Px,则正弦插值函数可表示成P(t/$-
n)。因此,当抽样间隔满足抽样定理时,有
F(u)=
]
n=-]
,将计算结果与理论解
析式进行对比,分析该算法的计算精度与误差分布
特点,并验证该算法的可行性。
E
]
F(n$)P(u/$-n)。
]
(6)
式(6)代入式(3)可得G(v)的近似
G(v)c=
-]
1 快速汉克尔变换
]
[2-3]
H(v-u)
E
F(n$)P(u/$-n)du=
Q
n=-]
n=-]
]
定义v阶汉克尔变换为
g(C)=
0
f(K)KJ
(KC)dK,v>-1,
Q
M
E
]
F(n$)H
M
(M-n$),
P(u/$)H
Q
M
*
(7)
(1)
*
*
其中,J
v
(x)为v阶第一类贝塞尔函数,g(C)为输入
函数f(K)的汉克尔变换,逆汉克尔变换与正汉克尔
变换完全对称。汉克尔变换广泛出现在带源的轴对
称问题中,如可控源大地电磁法(CSAMT)和瞬变电
磁法(TEM)的正演计算。
引入变量替换:K=e
-u
H
M
=
-]
(v-u)du,
其中,H
M
为滤波系数。
Ghos、Koefoed及Anderson等人由已知的含贝
塞尔函数积分恒等式,利用傅里叶变换法得到了汉
克尔变换的滤波系数,滤波抽样间隔为ln10/10,得
到的汉克尔滤波系数H
0
和H
1
分别为60和46个。
后来,一些学者对求取线性滤波系数的计算方法进
行了改进,其中Guptasama和Singh给出了61点和
/r
0
,C=r
0
e,u、MI(-],
x
M
]),r
0
为选定的常数。再定义复合函数F(u)=
f(K)K、G(v)=g(C)C、H
M
(x)=J
M
(x)e
,将其代入
收稿日期:2009-12-25
基金项目:国家863高技术研究发展计划项目(2006AA06Z110)
#754#
物 探 与 化 探34卷
120点的J
0
变换线性滤波系数以及47点和140点
的J
1
变换线性滤波系数。据Guptasarma和Singh测
评,该滤波系数的精度比其他早期的线性滤波器的
计算精度更高,计算公式为
G(v)=
波系数:n=61,a=-5.0825,s=1.E-1。
120点滤波系数选择为n=120,a=-8.3885,s=
9.E-2。图2、表1分别给出了零阶滤
波系数的计算结果和相对误差。
E
F(n$)H
i=1
a+(i-1)s
*
n
*
M
,(8)
n$=(1/r)@10,i=1,2,,,n。
数值计算的精度取决于积分区间的长度n、抽样点
的位置nK和滤波系数H
M
,通常n越大,计算精度
越高。
由于滤波系数在n$y]时呈指数衰减,式(8)
[4]
便可在有限滤波长度下近似求和计算。令v=m$,
m为整数,作变量代换nc=m-n,得到
Gc(m$)=
nc=N
1
N
2
E
F[(m-nc)$]H
M
(nc$)。(9)
*
图2 算例1的61点和120点系数计算结果
表1 零阶滤波系数相对误差一览
lgr
(-2,0.5)
(0.5,0.8)
0.82
0.85
0.88
0.91
0.94
0.97
61点滤波
绝对误差
<4.0E-6
<4.0E-6
1.84846E-6
1.55426E-6
2.53443E-7
8.24982E-7
1.01789E-7
5.11557E-7
相对误差/%
(0,6.3E-3)
(9.3E-3,1.49)
20.288553
85.926025
89.674484
2459.57226
35065.5546
292609.656
120点滤波
绝对误差
<6.0E-7
<3.0E-7
1.10096E-7
2.10929E-8
6.99633E-8
5.61145E-8
1.12878E-8
5.09593E-8
相对误差/%
(0,5.0E-4)
(2.0E-4,0.14)
1.208392
1.166112
24.75473
167.2979
388.8577
29148.62
从上式推导可以看出,由于贝赛尔函数随自变量缓慢
地振荡衰减,不能采用常规的数值积分方法进行计
算
[5]
,而快速汉克尔变换由于在计算中避免了直接数
值积分和Bessel函数的计算,只有简单的乘积和求和
运算,而且滤波系数只须计算一次,故具有很好的计
算精度和运算速度。
2 计算实例与精度验证
利用上述优化快速汉克尔变换理论,编制出一套
相应的计算机程序,分别使用Guptasarma和Singh给
出的0阶、1阶滤波系数进行快速汉克尔数值计算,
将计算结果与理论解析值进行对比。
2.1 零阶滤波系数的计算结果对比
]
2.2 一阶滤波系数的计算结果对比
]
算例2:f(r)=
析解为(r/4)e
-r
2
/4
0
Q
(Ke
2-K
2
)J
1
(rK)dr,其理论解
,如图3所示。
算例1:f(r)=
解为(1/2)e
-r
2
/4
0
Q
(Ke
-K
2
)J
0
(rK)dr,其理论解析
,如图1所示。
图3 算例2的解析式计算结果
Hankel变换数值解:
2-K
2
1
r
E
K(K)C,K(K)
i=1
iii
n
=
图1 算例1的解析式计算结果
Hankel变换数值解:
i
K
i
e,K
i
=(1/r)@10
1
K(K
i
)C
i
,K(K
i
)=
r
E
i=1
,i=1,2,,,n。61点滤
n
i
K
,K,i=1,,,n。47点滤
i
e
i
=(1/r)@10
波系数选取为n=47,a=-3.,s=
a+(i-1)s
1.1E-1。140点滤波系数:n=140,a=
-7.91001919000,s=8.79671439570E-2。表2、
图4给出了一阶滤波系数的计算结果和相对误差。
-K
2
a+(i-1)s
6期张伟等:汉克尔变换的数值计算与精度的对比
表2 一阶滤波系数相对误差一览
lgr
47点滤波
绝对误差
<4.0E-7
(3.7E-9,7.1)E-6
1.26234E-6
2.57783E-6
1.46030E-6
3.77645E-7
1.26715E-6
9.13894E-7
相对误差/%
(0,1.7E-2)
(9.5E-6,1.79)
4.194172
40.261058
136.22622
277.03173
10023.993
112026.34
绝对误差
<3.75E-8
(2.9E-8,6.0E-7)
2.05771E-7
2.00331E-7
3.30527E-8
1.54478E-7
6.55898E-8
6.25419E-8
140点滤波
#755#
相对误差/%
(1.3E-4,1.6E-4)
(1.6E-5,0.1)
0.683681
3.128803
3.083366
113.3219
518.8551
7666.463
(-2,-1)
(-1,0.8)
0.82
0.85
0.88
0.91
0.94
0.97
但在其自变量较小时,尽管其理论解析值较小,该算
法的计算精度也可达到10数量级,因此,该算法
适用性好。例如,该算法在瞬变电磁法的初期到中
晚期的计算具有较高的精度,而在晚期数据的计算
中,由于场值的衰减和仪器数据采集精度的限制,因
而一般舍去不参加计算。此外,试算结果表明,滤波
系数较多,计算精度越高,H
0
的120点较61点、H
1
的140点较47点的计算精度高出一个数量级。
参考文献:
[1] 朴化荣.电测测深法原理[M].北京:地质出版社,1995:112.
[2] itallinearfiltersforHankelJ0andJ1
transforms[J].GeophysicalProspecting,1997,45:745-76.
[3] 程志平.电法勘探教程[M].北京:冶金工业出版社,2007:44.
[4] 阮百尧.均匀水平大地上频率域垂直磁偶源电磁场数值滤波
解法[J].桂林工学院学报,2005,25(1).
[5] 胡俊,聂在平.索末菲尔德积分新方法)快速汉克尔变换[J].
电子学报,1998,26(3).
[6] erprogramnumericalintegrationofrelated
Hankeltransformsoforders0and1byadaptivedigitalfiltering
[J].Geophysic,1979,44(7).
-4
图4 算例2的解析式计算结果
3 结论
通过上述计算结果的对比可以看出,基于数字
滤波的快速汉克尔变换算法连续一致逼近其理论解
析解,不存在振荡现象,且计算精度高,两者的绝对
误差在10数量级以下。在理论解析值较大时,相
-6
对误差在10数量级以下;在理论解析值较小时,
由于其自身解很小,该算法的逼近精度很低,随着自
变量的增长,误差可以迅速地增长到10数量级上;
6
-6
RESEARCHANDAPPLICATIONONNUMERICALINTEGRATION
OFHANKELTRANSFORMSBYDIGITALFILTERING
ZHANGWe,iWANGXu-ben,QINQing-yan
(-Detection&InformationTechniquesofMinistryofEducation,ChengduUniversityofTechnology,Chengdu 610059,China)
Abstract:NumericalintegrationofHankeltransformsiseffectivetoolsforEMSoundings'forwardnumericalsimulation,thispaper
madeoutformuladerivationofnumericalintegrationofHankeltransformsbydigitalfiltering,andusedigitalcoefficientstodonumer-i
calcomputewhichwasputforwardbyGuptasamaandSingh,finallycontrastedtotheoreticalresolveexpressionandanalyzedthisalgo-
rithms'ultsshowthatthecalculationofthisalgorithmcontinuouslyapproximateitstheoreticalsolution,ithas
nooscillation,highprecisionandgreatpracticalvalueinthenumericalsimulationstudy.
Keywords:Hankeltransform,Fasthankeltransform(FHT),Filtercoefficients
作者简介:张伟(1983-),男,成都理工大学信息工程学院博士研究生。
2024年5月30日发(作者：翁茂德)
第34卷第6期
2010年12月
物 探 与 化 探
GEOPHYSICAL&GEOCHEMICALEXPLORATION
Vo.l34,No.6
Dec.,2010
汉克尔变换的数值计算与精度的对比
张伟,王绪本,覃庆炎
(成都理工大学地球探测与信息技术教育部重点实验室,四川成都 610059)
摘要:汉克尔变换的数值解法是计算电磁测深正演理论曲线的有效工具。对基于数字滤波法的快速汉克尔变换
算法进行了公式推导,并用Guptasama和Singh给出的线性滤波系数进行了计算,对比分析了该数值算法与理论解
析式的误差分布特点。结果表明,该数值算法的计算结果连续一致逼近其理论值,且不存在振荡现象,计算精度
高,在数值模拟研究中具有较大的实用价值。
关键词:汉克尔变换;快速汉克尔变换;滤波系数
中图分类号:P631 文献标识码:A 文章编号:1000-8918(2010)06-0753-03
在水平层状介质地表的电磁测深数值计算中,
除平面波和偶极场源波区的电磁场可用具有递推性
质的频率特性函数R(X)表示,有限发-收距电磁测
深正演解析式均用含零阶或一阶贝塞尔函数的积分
形式表示出来,这些积分形式是汉克尔变换式,汉克
尔变换的数值解法是电磁测深数值正演计算的基础
问题,因此具有重要的研究意义。
基于数字滤波法的快速汉克尔变换具有编程简
单、运算速度快、执行效率高等优点,该方法广泛地
应用在电磁法(TEM、CASMT、FIP等)正演数值计算
中。笔者采用Guptasama和Singh给出的J
0
和J
1
变
换线性滤波系数进行计算
[2]
式(1),可得
]
G(v)=
-]
F(u)H
Q
M
(v-u)du=F*H
M
。(3)
G是函数F和H
M
的褶积,因此G可看成输入函数F
和线性滤波函数H
M
的系统响应。
在满足抽样定理的情况下,即$<1/|2f
c
|,$为
抽样间隔,f
c
为被抽样函数的最高频率。引入正弦
插值函数
sin[(P/$)(t-n$)]
,并定义函数P(x)=
(P/$)(t-n$)
sin(Px)/Px,则正弦插值函数可表示成P(t/$-
n)。因此,当抽样间隔满足抽样定理时,有
F(u)=
]
n=-]
,将计算结果与理论解
析式进行对比,分析该算法的计算精度与误差分布
特点,并验证该算法的可行性。
E
]
F(n$)P(u/$-n)。
]
(6)
式(6)代入式(3)可得G(v)的近似
G(v)c=
-]
1 快速汉克尔变换
]
[2-3]
H(v-u)
E
F(n$)P(u/$-n)du=
Q
n=-]
n=-]
]
定义v阶汉克尔变换为
g(C)=
0
f(K)KJ
(KC)dK,v>-1,
Q
M
E
]
F(n$)H
M
(M-n$),
P(u/$)H
Q
M
*
(7)
(1)
*
*
其中,J
v
(x)为v阶第一类贝塞尔函数,g(C)为输入
函数f(K)的汉克尔变换,逆汉克尔变换与正汉克尔
变换完全对称。汉克尔变换广泛出现在带源的轴对
称问题中,如可控源大地电磁法(CSAMT)和瞬变电
磁法(TEM)的正演计算。
引入变量替换:K=e
-u
H
M
=
-]
(v-u)du,
其中,H
M
为滤波系数。
Ghos、Koefoed及Anderson等人由已知的含贝
塞尔函数积分恒等式,利用傅里叶变换法得到了汉
克尔变换的滤波系数,滤波抽样间隔为ln10/10,得
到的汉克尔滤波系数H
0
和H
1
分别为60和46个。
后来,一些学者对求取线性滤波系数的计算方法进
行了改进,其中Guptasama和Singh给出了61点和
/r
0
,C=r
0
e,u、MI(-],
x
M
]),r
0
为选定的常数。再定义复合函数F(u)=
f(K)K、G(v)=g(C)C、H
M
(x)=J
M
(x)e
,将其代入
收稿日期:2009-12-25
基金项目:国家863高技术研究发展计划项目(2006AA06Z110)
#754#
物 探 与 化 探34卷
120点的J
0
变换线性滤波系数以及47点和140点
的J
1
变换线性滤波系数。据Guptasarma和Singh测
评,该滤波系数的精度比其他早期的线性滤波器的
计算精度更高,计算公式为
G(v)=
波系数:n=61,a=-5.0825,s=1.E-1。
120点滤波系数选择为n=120,a=-8.3885,s=
9.E-2。图2、表1分别给出了零阶滤
波系数的计算结果和相对误差。
E
F(n$)H
i=1
a+(i-1)s
*
n
*
M
,(8)
n$=(1/r)@10,i=1,2,,,n。
数值计算的精度取决于积分区间的长度n、抽样点
的位置nK和滤波系数H
M
,通常n越大,计算精度
越高。
由于滤波系数在n$y]时呈指数衰减,式(8)
[4]
便可在有限滤波长度下近似求和计算。令v=m$,
m为整数,作变量代换nc=m-n,得到
Gc(m$)=
nc=N
1
N
2
E
F[(m-nc)$]H
M
(nc$)。(9)
*
图2 算例1的61点和120点系数计算结果
表1 零阶滤波系数相对误差一览
lgr
(-2,0.5)
(0.5,0.8)
0.82
0.85
0.88
0.91
0.94
0.97
61点滤波
绝对误差
<4.0E-6
<4.0E-6
1.84846E-6
1.55426E-6
2.53443E-7
8.24982E-7
1.01789E-7
5.11557E-7
相对误差/%
(0,6.3E-3)
(9.3E-3,1.49)
20.288553
85.926025
89.674484
2459.57226
35065.5546
292609.656
120点滤波
绝对误差
<6.0E-7
<3.0E-7
1.10096E-7
2.10929E-8
6.99633E-8
5.61145E-8
1.12878E-8
5.09593E-8
相对误差/%
(0,5.0E-4)
(2.0E-4,0.14)
1.208392
1.166112
24.75473
167.2979
388.8577
29148.62
从上式推导可以看出,由于贝赛尔函数随自变量缓慢
地振荡衰减,不能采用常规的数值积分方法进行计
算
[5]
,而快速汉克尔变换由于在计算中避免了直接数
值积分和Bessel函数的计算,只有简单的乘积和求和
运算,而且滤波系数只须计算一次,故具有很好的计
算精度和运算速度。
2 计算实例与精度验证
利用上述优化快速汉克尔变换理论,编制出一套
相应的计算机程序,分别使用Guptasarma和Singh给
出的0阶、1阶滤波系数进行快速汉克尔数值计算,
将计算结果与理论解析值进行对比。
2.1 零阶滤波系数的计算结果对比
]
2.2 一阶滤波系数的计算结果对比
]
算例2:f(r)=
析解为(r/4)e
-r
2
/4
0
Q
(Ke
2-K
2
)J
1
(rK)dr,其理论解
,如图3所示。
算例1:f(r)=
解为(1/2)e
-r
2
/4
0
Q
(Ke
-K
2
)J
0
(rK)dr,其理论解析
,如图1所示。
图3 算例2的解析式计算结果
Hankel变换数值解:
2-K
2
1
r
E
K(K)C,K(K)
i=1
iii
n
=
图1 算例1的解析式计算结果
Hankel变换数值解:
i
K
i
e,K
i
=(1/r)@10
1
K(K
i
)C
i
,K(K
i
)=
r
E
i=1
,i=1,2,,,n。61点滤
n
i
K
,K,i=1,,,n。47点滤
i
e
i
=(1/r)@10
波系数选取为n=47,a=-3.,s=
a+(i-1)s
1.1E-1。140点滤波系数:n=140,a=
-7.91001919000,s=8.79671439570E-2。表2、
图4给出了一阶滤波系数的计算结果和相对误差。
-K
2
a+(i-1)s
6期张伟等:汉克尔变换的数值计算与精度的对比
表2 一阶滤波系数相对误差一览
lgr
47点滤波
绝对误差
<4.0E-7
(3.7E-9,7.1)E-6
1.26234E-6
2.57783E-6
1.46030E-6
3.77645E-7
1.26715E-6
9.13894E-7
相对误差/%
(0,1.7E-2)
(9.5E-6,1.79)
4.194172
40.261058
136.22622
277.03173
10023.993
112026.34
绝对误差
<3.75E-8
(2.9E-8,6.0E-7)
2.05771E-7
2.00331E-7
3.30527E-8
1.54478E-7
6.55898E-8
6.25419E-8
140点滤波
#755#
相对误差/%
(1.3E-4,1.6E-4)
(1.6E-5,0.1)
0.683681
3.128803
3.083366
113.3219
518.8551
7666.463
(-2,-1)
(-1,0.8)
0.82
0.85
0.88
0.91
0.94
0.97
但在其自变量较小时,尽管其理论解析值较小,该算
法的计算精度也可达到10数量级,因此,该算法
适用性好。例如,该算法在瞬变电磁法的初期到中
晚期的计算具有较高的精度,而在晚期数据的计算
中,由于场值的衰减和仪器数据采集精度的限制,因
而一般舍去不参加计算。此外,试算结果表明,滤波
系数较多,计算精度越高,H
0
的120点较61点、H
1
的140点较47点的计算精度高出一个数量级。
参考文献:
[1] 朴化荣.电测测深法原理[M].北京:地质出版社,1995:112.
[2] itallinearfiltersforHankelJ0andJ1
transforms[J].GeophysicalProspecting,1997,45:745-76.
[3] 程志平.电法勘探教程[M].北京:冶金工业出版社,2007:44.
[4] 阮百尧.均匀水平大地上频率域垂直磁偶源电磁场数值滤波
解法[J].桂林工学院学报,2005,25(1).
[5] 胡俊,聂在平.索末菲尔德积分新方法)快速汉克尔变换[J].
电子学报,1998,26(3).
[6] erprogramnumericalintegrationofrelated
Hankeltransformsoforders0and1byadaptivedigitalfiltering
[J].Geophysic,1979,44(7).
-4
图4 算例2的解析式计算结果
3 结论
通过上述计算结果的对比可以看出,基于数字
滤波的快速汉克尔变换算法连续一致逼近其理论解
析解,不存在振荡现象,且计算精度高,两者的绝对
误差在10数量级以下。在理论解析值较大时,相
-6
对误差在10数量级以下;在理论解析值较小时,
由于其自身解很小,该算法的逼近精度很低,随着自
变量的增长,误差可以迅速地增长到10数量级上;
6
-6
RESEARCHANDAPPLICATIONONNUMERICALINTEGRATION
OFHANKELTRANSFORMSBYDIGITALFILTERING
ZHANGWe,iWANGXu-ben,QINQing-yan
(-Detection&InformationTechniquesofMinistryofEducation,ChengduUniversityofTechnology,Chengdu 610059,China)
Abstract:NumericalintegrationofHankeltransformsiseffectivetoolsforEMSoundings'forwardnumericalsimulation,thispaper
madeoutformuladerivationofnumericalintegrationofHankeltransformsbydigitalfiltering,andusedigitalcoefficientstodonumer-i
calcomputewhichwasputforwardbyGuptasamaandSingh,finallycontrastedtotheoreticalresolveexpressionandanalyzedthisalgo-
rithms'ultsshowthatthecalculationofthisalgorithmcontinuouslyapproximateitstheoreticalsolution,ithas
nooscillation,highprecisionandgreatpracticalvalueinthenumericalsimulationstudy.
Keywords:Hankeltransform,Fasthankeltransform(FHT),Filtercoefficients
作者简介:张伟(1983-),男,成都理工大学信息工程学院博士研究生。