论文-Hankel变换的数值积分算法及误差分析

Hankel变换的数值积分算法及误差分析张健杨紫恒杨发康(甘肃省地矿局第一地质矿产勘查院甘肃天水741020)摘要:快速Hankel变换方法是用数值方法计算电磁测深正演理论曲线最有效的工具。本文主要介绍了快速Hankel变换方法在电磁法勘探数值模拟方面的应用，引用Hankel变换的数值积分算法，其基本思路是，采用高斯求积方法对有限区间积分进行计算，进而结合连分式展开方法加快部分积分和的收敛速度；基于Matlab算法平台，进行算例计算分析，验证了Hankel变换数值积分算法的有效性，体现了Hankel变换数值积分算法简单易懂,计算次数少,精度高等特点。关键词：Hankel变换，电磁法勘探，数值积分算法0引言Hankel变换0)()()(drprJrrfpfnn广泛应用于信号处理各领域，是一个有效的数学工具。在电磁法勘探数值模拟方面，对于有线发--收距电磁测深正演公式均用含零阶或一阶贝塞尔函数的积分形式表示出来。实际上，这些积分形式的正演公式是一种Hankel变换式。1Hankel变换应用从20世纪70年代初开始发展起来的快速Hankel变换方法是用数值方法计算电磁测深正演理论曲线的最有效的工具。电磁测深频率域电磁响应一般可表示为如下含有第1类i阶Bessel函数的积分形式：drJPKPrfi)(),(),(0(1-1)式中),(PK为积分核函数，P为与地层有关的层参数，)(rJi是第1类i阶Bessel函数。对这个无穷积分直接数值求积是困难的，因为积分核中含有振荡的贝塞尔函数，故而多采用数字滤波技术计算该积分。当核函数单调、快速衰减，且发收距合适时，数字滤波确能取得较好的效果。但该方法采用极为有限的滤波系数来进行积分截断，人为的主观臆断性太强，大大损失了其计算精度。2Hankel变换数值积分算法汉克尔变换的数值积分算法，其基本思路是，采用高斯求积方法对有限区间积分进行计算，进而结合连分式展开方法加快部分积分和的收敛速度。一反数字滤波法的常规做法，数值积分算法使用期望误差（绝对误差/相对误差）来控制计算是继续或是终止。为了直接采用数值求积方法计算（1-1）式，可将该积分写成如下部分积分之和NnnPSf0(1-2)drJPKPnnZZin1)(),((1-3)式中iJr为第1类的i阶Bessel函数，nZ为iJr的第n个零点用距离r归一化后的值。N为参与求和的部分积分项的个数，取决于计算方法。nP可采用Gauss求积公式计算1MnjjijjPhKa,PJar(1-4)式（1-4）中M是求积节点个数，ja为求积节点，jh为求积系数。一般采用7M点求积公式即可达到计算精度。采用连分式加速式(1-2)积分项求和的收敛。01231111ndSdddd(1-5)式中012Nd,d,d,,d可由区间积分项nP求得。本文采用数值积分方法，在matlab计算平台中，预先计算了贝赛尔函数01iJr,i,的300个零点，最初的几个零点之间的区间积分采用matlab的自适应Simpsons积分函数(QUAD)，之后的积分区间根据计算精度，采用自适应Lobatto积分函数(QUADL)，积分项的求和采用Euler法加速其收敛，积分的终止条件是最后加入的区间积分项小于设定的容许误差，或者已计算了设定的所有区间数。实际上计算不超过30个过零点区间就可达到所需要的精度。3算例分析下面设计算例验证算法的有效性。2120111keJkdk(1-6)Matlab计算程序如下r=logspace(-4,9,256);FUNC=@(x)exp(-x);fori=1:length(r)sum(i)=NJCST('J1',FUNC,r(i),1e-11,1e-11,300);endycx=(sqrt(r.^2+1)-1)./(r.*sqrt(r.^2+1));loglog(r,sum,'r-',r,ycx,'b-.');figuresemilogx(r,(sum-ycx))NJCST是汉克尔变换函数。sum是式(1-6)左端积分数值计算结果，ycx是式(1-6)右端解析式计算结果。变量r代表，函数随的变化如图1-1(a)所示，图1-1(b)是数值解与解析解的绝对误差随的变化，表1-1给出了25项数据对比，其中绝对误差为解析解与数值解差值的绝对值。图1-1数值计算结果与解析解对比(a)图1-1数值计算结果与解析解绝对误差(b)表1-1计算数据对比序号ρ绝对误差16.31E-035.83E-1326.98E-034.63E-1237.72E-033.41E-1348.55E-035.53E-1459.46E-031.77E-1361.05E-025.86E-1371.16E-021.69E-1381.28E-021.52E-1210-2100102104106108-1-0.500.511.522.53x10-11ρ绝对误差(b)10-410-2100102104106108101010-1010-810-610-410-2100ρf(ρ)数值解解析解91.42E-021.98E-12101.57E-021.46E-14111.73E-029.36E-14121.92E-022.63E-13132.12E-026.43E-13142.35E-028.46E-13152.60E-028.62E-12162.88E-022.41E-13173.18E-024.82E-14183.52E-021.78E-13193.90E-023.70E-13204.31E-021.32E-13214.77E-027.36E-13225.28E-023.74E-13235.84E-027.17E-14246.46E-022.21E-13257.15E-028.39E-13第一幅图中红线是利用Hankel变换数值积分算法计算的结果，蓝线是解析解，第二幅图中是Hankel变换数值积分算法和解析解的误差对比，在logr在10-1到10-2之间最大误差是4*10-11除此之外，在103之后，两者几乎一致，没有差别，很精确。4结论Hankel变换数值积分方法广泛应用于信号处理各领域，是一个有效的数学工具，简单易懂,计算次数少,精度高,具有很好的应用价值。参考文献［1］朴化荣电磁测深法原理地质出版社［2］高俊超李婧Hankel变换的数值积分及应用天津理工大学学报作者简介：张健(1987-)，男，工程师，2010年毕业于长安大学勘查技术与工程专业，从事物探找矿工作通讯地址：甘肃省天水市麦积区马跑泉路54号电话：[phone]