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

Hankel变换的数值积分算法及误差分析 张 健 杨紫恒 杨发康 甘肃省地矿局第一地质矿产勘查院 甘肃 天水 741020 摘要: 快速Hankel变换方法是用数值方法计算电磁测深正演理论曲线最有效的工具。本文主要介绍了快速Hankel变换方法在电磁法勘探数值模拟方面的应用，引用Hankel变换的数值积分算法，其基本思路是，采用高斯求积方法对有限区间积分进行计算，进而结合连分式展开方法加快部分积分和的收敛速度；基于Matlab算法平台，进行算例计算分析，验证了Hankel变换数值积分算法的有效性，体现了Hankel变换数值积分算法简单易懂,计算次数少,精度高等特点。 关键词： Hankel变换，电磁法勘探，数值积分算法 0 引言 Hankel变换广泛应用于信号处理各领域，是一个有效的数学工具。在电磁法勘探数值模拟方面，对于有线发--收距电磁测深正演公式均用含零阶或一阶贝塞尔函数的积分形式表示出来。实际上，这些积分形式的正演公式是一种Hankel变换式。 1 Hankel变换应用 从20世纪70年代初开始发展起来的快速Hankel变换方法是用数值方法计算电磁测深正演理论曲线的最有效的工具。电磁测深频率域电磁响应一般可表示为如下含有第1类i 阶Bessel函数的积分形式： 1-1 式中为积分核函数，为与地层有关的层参数，是第1类i阶Bessel函数。对这个无穷积分直接数值求积是困难的，因为积分核中含有振荡的贝塞尔函数，故而多采用数字滤波技术计算该积分。当核函数单调、快速衰减，且发收距合适时，数字滤波确能取得较好的效果。但该方法采用极为有限的滤波系数来进行积分截断，人为的主观臆断性太强，大大损失了其计算精度。 2 Hankel变换数值积分算法 汉克尔变换的数值积分算法，其基本思路是，采用高斯求积方法对有限区间积分进行计算，进而结合连分式展开方法加快部分积分和的收敛速度。一反数字滤波法的常规做法，数值积分算法使用期望误差（绝对误差/相对误差）来控制计算是继续或是终止。为了直接采用数值求积方法计算（1-1）式，可将该积分写成如下部分积分之和 1-2 1-3 式中为第1类的i阶Bessel函数，为的第个零点用距离归一化后的值。为参与求和的部分积分项的个数，取决于计算方法。可采用Gauss求积公式计算 1-4 式（1-4）中是求积节点个数，为求积节点，为求积系数。一般采用点求积公式即可达到计算精度。 采用连分式加速式 1-2 积分项求和的收敛。 1-5 式中可由区间积分项求得。 本文采用数值积分方法，在matlab计算平台中，预先计算了贝赛尔函数的300个零点，最初的几个零点之间的区间积分采用matlab 的自适应Simpsons积分函数 QUAD ，之后的积分区间根据计算精度，采用自适应Lobatto积分函数 QUADL ，积分项的求和采用Euler法加速其收敛，积分的终止条件是最后加入的区间积分项小于设定的容许误差，或者已计算了设定的所有区间数。实际上计算不超过30个过零点区间就可达到所需要的精度。 3 算例分析 下面设计算例验证算法的有效性。 1-6 Matlab计算程序如下 r logspace -4,9,256 ; FUNC @ x exp -x ; for i 1: length r sum i NJCST J1, FUNC, r i ,1e-11, 1e-11, 300 ; end ycx sqrt r.^2+1 -1 ./ r.*sqrt r.^2+1 ; loglog r,sum,r-,r,ycx,b-. ; figure semilogx 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 计算数据对比 序号 ρ 绝对误差 1 6.31E-03 5.83E-13 2 6.98E-03 4.63E-12 3 7.72E-03 3.41E-13 4 8.55E-03 5.53E-14 5 9.46E-03 1.77E-13 6 1.05E-02 5.86E-13 7 1.16E-02 1.69E-13 8 1.28E-02 1.52E-12 9 1.42E-02 1.98E-12 10 1.57E-02 1.46E-14 11 1.73E-02 9.36E-14 12 1.92E-02 2.63E-1