作者简介: 李锋平(1989-),男,江西井冈山人,汉族,在读硕士研究生,研究方向为电磁法勘探理论与应用。E-mail: geophysics_Lee@163.com
在瞬变电磁的一维正演响应模拟中,常用的方法是先在频率域中求解,之后将结果转换到时间域。但该方法在晚期的计算精度通常不高,因此,使用5种频-时域转换方法(正、余弦变换的数值滤波算法,G-S逆拉普拉斯算法、正、余弦变换的折线逼近法)进行了计算,与解析解对比,得出余弦变换的数值滤波算法在晚期计算中精度最高的结论,并对这五种转换方法产生误差问题的原因进行了讨论和分析。本研究有利于瞬变电磁一维正演响应的高精度计算,使其在多维计算中得到更好的应用。
In the forward modeling of transient electromagnetic response, the commonly used method is to calculate the response in the frequency domain, and then the results were converted into the time domain. In this paper, five frequency-time conversion methods (Digital filtering algorithms of sine and cosine transform, G-S inverse Laplace algorithm, polygonal approximations method of sine and cosine transform) were used to conduct calculation. it is concluded that the digital filtering method based on cosine transform has the highest computing precision among these five methods, as shown by the comparison with the analytical algorithm, and the reasons of the error caused by these five kinds of conversion methods were discussed and analyzed, which is beneficial to the high precision calculation of one-dimensional forward modeling of transient electromagnetic response, so that these methods can be better applied to the calculation of multi-dimensional frequency-time conversion.
瞬变电磁法正演响应的计算一般都是先求得频率域的解, 然后由频率域转换到时间域得到瞬态解 然后由频率域转换到时间域得到瞬态解[1-5]。目前, 可用于频
Gaver-Stehfest算法的精度与计算机环境字长密切相关, 字长越大精度越高, 因此该算法对计算机硬件和软件要求较高
为此, 笔者针对中心回线装置, 在均匀半空间情况下, 使用多种频 -时域转换方法(正、余弦变换的数值滤波算法, G-S逆拉普拉斯算法, 正、余弦变换的折线逼近法)计算均匀大地表面的感应电压值, 然后将计算结果与均匀半空间的解析解进行了对比分析, 优选算法, 力求对其做出一个全面细致的归纳总结, 从而使其在瞬变电磁多维正演计算中得到更好的应用。
1 几种频-时域转换方法的简介
Gaver-Stehfest逆拉普拉斯算法的主要思想即利用δ 函数的基本积分性质, 让变换函数本身和它跟δ 函数乘积的积分联系在一起, 然后选择恰当的δ 函数使得积分形式转换为和变换函数的拉普拉斯变换相关的形式。这样就可利用熟悉的拉普拉斯变换式求取原函数, 从而达到逆变换的目的[26]。
1982年, Knight及Raiche 首先提出了Gaver-Stehfest 算法。朴化荣教授随后提出了Gaver-Stehfest 逆拉普拉斯算法, 其计算表达式为[27]
在阶跃电流激励的情况下, 中心回线装置感应电动势的计算表达式为[28]
中心回线装置情况下, 二次场磁场分量对时间的导数的可以表示为[29]
因此, 根据法拉第感应定律, 中心回线装置感应电动势的计算表达式还可表述为
对式(2)进行Gaver-Stehfest 逆拉普拉斯变换, 可以得到感应电动势V(t)的表达式
式中:j=1, 2, 3, …; N是抽样点的位置; q是有效接收面积; c=ln2/t, 是采样间隔; W(j, N)是Gaver-Stehfest变换的系数。
折线逼近法的思想是把足够光滑的F(ω )函数曲线分成N段(N的大小取决于函数F(ω )的分布形态和计算精度), 并在每个段内用折线去进行逼近F(ω )。若函数F(ω )的曲线愈是光滑, 分段越是细致, 其逼近误差即可达到足够小[30]。
根据前面所述原理, 其数值计算表达式分别为[31]
其中, a和b为积分的上、下限, ω 为采样角频率, t为采样时间。则中心回线装置的感应电动势的正、余弦折线逼近的计算近似表达式分别为:
正弦变换的折线逼近形式
余弦变换的折线逼近形式
根据美国地球物理学家P.Raab和F.Frisehknecht推出的均匀半空间中心回线装置的感应电压的解析表达式[28]
式中:z=
令发射回线边长为100 m, 接收面积为1 m2, 置于电阻率为100 Ω · m的均匀大地表面上, 供电电流为 1 A, 在t=0时切断, 得到解析解和以上5种频
![]() | 表1 半空间下5种算法计算结果的平均相对误差 % |
从图1可以看出, 除了正弦变换的折线逼近法, 其余四种转换方法的计算结果曲线和解析解的曲线都比较吻合, 同时, G-S逆拉普拉斯算法和正弦变换的折线逼近法的计算结果晚期部分出现了负值现象(如图中所示的负值取绝对值1和负值取绝对值2曲线, 结果为负数, 取绝对值后呈现在图上。), 而其余方法则没有这种情况。从图2和图3以及表1可以看出, 除余弦变换的滤波算法在早期高频段(106~108 Hz)出现了一定相对误差, 其余方法都有很高的转换精度, 但瞬变电磁法实际观测的采样时间范围一般在n× (10-3~102) ms, 所以这个影响可以忽略; 另外, 从以上图表不难看出, 正弦变换的折线逼近法的相对误差是最大的, 特别是在晚期时段, 最高达到了接近380× 100%的相对误差, 而余弦变换的滤波算法的转换精度在10-6~1 s时间段是最好的, 相对误差接近零, 和解析解曲线基本吻合, 在10-8~1 s或是10-6~1 s时间段的平均相对误差也分别只有0.83%和 0.08%, G-S逆拉普拉斯算法、正弦变换的滤波算法、余弦变换的折线逼近法依序次之。
下面对5种准换方法产生误差问题的原因进行讨论。对于G-S逆拉普拉斯算法, 其要求双精度运算, 李建慧等认为G-S算法的精度与滤波系数J的取值、计算机字长有关。计算机字长越大, J的取值相应变大, G-S算法的精度就会越高[34]。通过试算, 发现其计算精度还与大地电阻率值有关, 图4为在滤波系数值J相同的情况下, 大地电阻率分别为10、100、1 000 Ω · m时, G-S逆拉普拉斯算法计算的感应电压相比于解析解的相对误差曲线, 从图中可以看出, 随着大地电阻率值的增大, 其相对误差也相应增大; 另外, 图5为在滤波系数值J不同的情况下, G-S逆拉普拉斯算法计算的感应电压相比于解析解的相对误差曲线, J的取值分别为10、12、14, 当J=12时, 其计算结果的相对误差相对较小。因此, 在实际运用中应根据所用计算机的有效字长选取最适合的滤波系数。
从式(6)和式(7)可以看出, 正、余弦变换的折线逼近法中含有sin(ω t)和cos(ω t)项, 当角频率ω 和采样时间t比较大时, sin(ω t)和cos(ω t)会强烈震荡(如图6所示, 其中ω =105), 这样, 就需要将F(ω )函数曲线分成更多的段来逼近原函数, 分段越多, 其计算精度越高。但是这种分段还取决于F(ω )函数的分布形态, 不是很好确定。在瞬变电磁中, 分段越多, 意味着采样时刻点个数的增多, 其个数可能超过几千甚至几万个。从图7可以看出, 对于正弦变换的折线逼近法, 随着分段数N的增大, 其相对误差并没有太大变化, 尤其是在10-1 s后, 其差别不是很大; 而对于余弦变换的折线逼近法, N从200增大到6 000时, 其整体相对误差的改变还是很可观的, 但当N从6 000增大到 50 000时, 在10-1s后, 其相对误差也没有太大的变化。由于这两种方法需要的采样点数都相对比较大, 在实际运用中不可取。
正、余弦变换的滤波算法都是基于正、余弦变换是汉克尔变换的特例, 运用成熟的汉克尔变换理论推导而得, 从理论上讲, 两者都是等效的, 但两者的转换结果却相差很大, 特别是在低频晚期段, 余弦变换的滤波算法精度很高, 而在高频段出现了一定的相对误差, 正弦变换的滤波算法则反之。原因是在频域到时域的转换时, 正弦变换的滤波算法计算(公式(13))采用的核函数是垂直磁场的虚分量, 即ImHz(ω ), 此核函数在早期高频段(108~106 Hz)函数曲线比较光滑, 斜率变化不大(如图8a所示, 图中虚线为负数取绝对值后的数), 所以, 在早期高频段, 计算精度比较理想; 但随着频率的减小, 在(106~1 Hz)频段, 曲线斜率迅速增大(如图8b所示, 图中虚线为负数取绝对值后的数), 虽然函数曲线依然光滑, 但因是负数, 所以导致计算结果会缓慢大于解析解, 因而在晚期低频段出现了较大的相对误差。反观余弦变换的滤波算法(公式(14)), 采用的核函数是垂直磁场的实分量, 即ReHz(ω ), 此核函数在早期高频段(108~106 Hz)出现了正负交替剧烈震荡的现象(如图8a所示), 随着频率的继续减小, 在(106~1 Hz)频段其曲线变得很光滑(如图8b所示), 104 Hz之后函数曲线的斜率几乎趋于零, 其值不再发生变化。因此, 余弦变换的滤波算法的计算结果在晚期低频段精度很高, 且瞬变电磁法实际观测的采样时间范围一般在n× (10-3~102) ms, 在这个范围内余弦变换的滤波算法的相对误差近乎为零, 所以在高频段出现的相对误差不会影响实际观测。
1) 正弦变换的数值滤波算法由于核函数 Im(Hz(ω ))的性质, 在早期高频段的计算有优势, 精度较高, 但是在晚期低频段计算结果出现了一定的相对误差, 转换精度不是很理想。
2) 余弦变换的数值滤波算法的转换精度在瞬变电磁法实际的采样时间范围内明显高于其它四种方法, 可以满足对转换精度的要求。
3) G-S逆拉普拉斯算法转换精度与滤波系数J的取值、计算机字长、大地电阻率等有关, 但其计算精度总体上来说较好, 仅次于余弦变换的数值滤波算法, 且计算速度非常快, 在实际运用中比较有优势。
4) 正弦变换的折线逼近法在晚期低频段的计算结果相对误差非常大, 且加大采样频点数其转换精度提升得不是很明显, 其计算速度也较慢, 不适合在实际应用中使用。
5) 余弦变换的折线逼近法在晚期低频段的计算结果相对误差明显小于正弦变换的折线逼近法的计算结果, 虽然加大采样频点数其转换精度提升的比较大, 但增加了运算时间, 也不适合在实际运用中使用。
在瞬变电磁一维正演响应计算时, 首先必须确保频率域的计算结果有足够的精度, 然后在此基础上采用更高计算精度, 效率更好的频
The authors have declared that no competing interests exist.
[1] |
|
[2] |
|
[3] |
|
[4] |
|
[5] |
|
[6] |
|
[7] |
|
[8] |
|
[9] |
|
[10] |
|
[11] |
|
[12] |
|
[13] |
|
[14] |
|
[15] |
|
[16] |
|
[17] |
|
[18] |
|
[19] |
|
[20] |
|
[21] |
|
[22] |
|
[23] |
|
[24] |
|
[25] |
|
[26] |
|
[27] |
|
[28] |
|
[29] |
|
[30] |
|
[31] |
|
[32] |
|
[33] |
|
[34] |
|