高阶时域有限差分法的多偏移距电磁波数值成像
汪科1,2
1.广西科技大学 电气与信息工程学院,广西 柳州 545006
2.湖南铁道职业技术学院,湖南 株洲 412001

作者简介: 汪科(1978-),男,讲师,主要研究方向为电磁场及信号分析处理。

摘要

数值模拟是数据处理和反演解释的重要环节,为提高高频电磁波层析成像技术在工程勘察中的解译准确性,采用高阶时间域有限差分法(FDTD(2,4))模拟电磁波在溶洞、断层和孤石等3种地电模型下的成像过程,并采用多偏移距模式获取电磁波道集。通过分析电磁波的波场快照结果和成像结果,总结了电磁波在不同地质体中的传播规律,为实际工程勘察工作提供可靠依据。试验结果表明,FDTD(2,4)能高精度模拟电磁波在复杂介质中的数值成像,同时也验证了多偏移距采集模式具有高效性和灵活性。

关键词: 高阶时域有限差分; 多偏移距; 电磁波层析; 数值成像; 工程勘察
中图分类号:P631 文献标志码:A 文章编号:1000-8918(2017)03-0489-07
High-order FDTD method for multi-offset electromagnetic numerical imaging
WANG Ke1,2
1. College of Electrical and Information Engineering, Guangxi University of Science and Technology, Liuzhou 545006, China
2. Hunan Railway Professional Technology College, Zhuzhou 412001, China
Abstract

Numerical simulation is an important part of data processing and inversion interpretation. In order to improve the efficiency and the accuracy of electromagnetic wave imaging technology in engineering survey, it is important to understand the propagation characteristics of electromagnetic waves by forward modeling. In this study, the finite-difference time-domain FDTD (2, 4) was used to simulate the numerical imaging of electromagnetic waves in three kinds of geoelectric models, i.e., karst caves, faults and boulders. The multi-offset mode was used to acquire the electromagnetic wave gathers. Based on the analysis of the electromagnetic wave wavefield and the waveform profile, the authors summarized the propagation law of electromagnetic waves in different geological bodies and provided reliable basis for practical engineering investigation. The experimental results show that FDTD (2, 4) can simulate the numerical imaging of electromagnetic waves in complex media with high accuracy, and also prove that the multi-offset acquisition mode has high efficiency and flexibility.

Keyword: high-order FDTD; multiple offset; electromagnetic wave tomography; numerical imaging; engineering investigation
0 引言

电磁波层析成像技术(电磁波CT)是工程勘察的重要方法之一, 它具有无损、高精度的特点[1, 2, 3]。高频电磁波CT通常采用分离式天线, 即发射天线和接收天线分离布置。电磁波CT的数据处理方法包括走时层析成像、振幅层析成像和相位层析成像, 其中走时层析成像应用最广[4, 5]。直达波初至时的拾取直接影响走时层析成像的反演结果, 而实际生产中初至时的拾取往往具有一定的主观性, 并且地球物理反演中普遍存在多解性, 因此解译者必须掌握电磁波的传播规律才能得到正确的勘探结果。笔者采用数值成像的方法模拟了高频电磁波CT在溶洞、断层、孤石等3种典型地电模型中的应用, 利用多偏移距数据采集模式, 提高了探测效率和目标针对性。在此基础上, 分析、总结了电磁波的特征响应和传播规律, 为实际勘探工作提供可靠依据。

本文数值模拟方法采用高阶时域有限差分法(FDTD)。FDTD是1966年由K.S.Yee[6, 7]提出的一种求解电磁场问题的数值计算方法, 它已广泛应用于各类电磁法的数值模拟中。FDTD的基本原理是在Yee元胞下将Maxwell旋度方程进行中心差分近似, 进而得到时间和空间上相互交替的电场分量E和磁场分量H。传统FDTD在时间和空间上都采用一阶偏导数2阶精度的中心差分近似(FDTD(2, 2)), 近似处理会导致数值色散, 并且数值色散程度与时间和空间的离散步长有关。FDTD(2, 2)通常采取减小空间和时间步长的方法减小数值色散, 但是这样会大大增加网格数目, 从而降低计算效率。为提高数值计算的精度, 同时保证计算效率, Fang等[8]提出了高阶算法, 高阶FDTD的基本思路是时间上采用与传统FDTD相同的2阶精度近似, 空间上采用高阶(2M)精度近似, 并且大量试验已证明高阶FDTD具有良好的数值色散特性[9, 10, 11, 12]。文中推导了FDTD(2, 4)算法的差分公式, 并完成了多偏移距电磁波CT的数值成像。

1 FDTD(2, 4)的基本原理
1.1 差分近似公式

FDTD基于Maxwell旋度方程进行差分离散, 然后沿时间轴逐步推进求解电磁场, Maxwell 旋度方程为[13, 14]:

其中:E为电场强度, H为磁场强度, ε 为介电常数, μ 为磁导率, σ 为电导率, σ m为等效磁阻率。

本文电磁波数值模拟基于TE波, 也就是只有ExEyHz3个分量, 则直角坐标系中Maxwell的分量:

Hzy=εExt+σEx, (3)-Hzx=εEyt+σEy, (4)Eyx-Exy=-μHzt-σmHx(5)

在上式离散中, 电场强度B和磁场强度H遵循Yee元胞排布原则。令f(x, y, z, t)代表EH在直角坐标系的某一个分量, 那么它在时间域和空间域的离散表示为

f(x, y, z, t)=f(iΔx, jΔy, kΔz, nΔt)=fn(i, j, k), (6)

f(x)进行泰勒级数展开:

fx±Δx2=f(x)±f'(x)Δx2+12!f″(x)Δx22±13!f‴(x)Δx23+14!f4(x)Δx24+O[Δx)5], (7)fx±3Δx2=f(x)±f'(x)3Δx2+12!f″(x)3Δx22±13!f‴(x)3Δx23+14!f4(x)3Δx24+O(Δx)5), (8)

式中, O[(Δ x)n] 代表与括号中的量具有相同数量级。由式(7)、(8)可推导出一阶偏导的2阶和4阶精度近似式:

f(x)x=fx+12Δx-fx-12ΔxΔx+O(Δx)2), (9)f(x)x=98fx+12Δx-fx-12ΔxΔx- 124fx+32Δx-fx-32ΔxΔx+O(Δx)4)(10)

FDTD(2, 4)在时间上采用一阶偏导数 2 阶近似精度, 在空间上采用一阶偏导数 4 阶精度, 由式(9)、(10)得:

f(x, y, z, t)t=fn+12(i, j, k)-fn-12(i, j, k)Δt, (11)f(x, y, z, t)x=98fni+12, j, k-fi-12, j, kΔx-124fni+32, j, k-fi-32, j, kΔx, (12)f(x, y, z, t)y=98fni, j+12, k-fi, j-12, kΔx-124fni, j+32, k-fi, j-32, kΔx(13)

设观测点Ex(i+1/2, j)、Ey(i, j+1/2)、Hz(i+1/2, j+1/2)、t=(n+1/2)Δ t, 由上式得TE模式下FDTD(2, 4)公式:

Exn+1i+12, j=CA(m)·Exni+12, j+CB(m)·98[Hzn+12i+12, j+12Δy-Hzn+12i+12, j-12Δy]-CB(m)·124[Hzn+12i+12, j+32Δy-Hzn+12i+12, j-32Δy], (14)Eyn+1i, j+12=CA(m)·Eyni, j+12-CB(m)·98[Hzn+12i+12, j+12Δx-Hzn+12i-12, j+12Δx]+CB(m)·124[Hzn+12i+32, j+12Δx-Hzn+12i-32, j+12Δx], (15)Hzn+12i+12, j+12=CP(m)·Hzn-12i+12, j+12-CQ(m)·[98Eyni+1, j+12-Eyni, j+12Δx-124Eyni+2, j+12-Eyni+1, j+12Δx-98Exni+12, j+1-Ezni+12, jΔy+124Exni+12, j+2-Ezni+12, j+1Δy], (16)

观测点的坐标表示为m=(x, y), 其中

CA(m)=1-σ(m)Δt2ε(m)/1+σ(m)Δt2ε(m), CB(m)=Δtε(m)/1+σ(m)Δt2ε(m), CP(m)=1-σm(m)Δt2μ(m)/1+σm(m)Δt2μ(m), Q(m)=Δtμ(m)/1+σm(m)Δt2μ(m)

1.2 数值稳定条件和边界条件

理想情况下, 离散的空间步长(Δ x、Δ y)和时间步长(Δ t)越大, 计算速度就越快。但是过大的数值空间步长会引起数值色散, 过大的时间步长会导致计算不稳定。因此为减小数值色散, 空间稳定性条件满足

Δxλmin10, (17)

其中λ min为介质中的最小波长。另外, 为保证计算稳定性, 时间步长上满足Courant稳定性条件, 其中c是真空中的电磁波速度:

Δt=1c1Δx)2+Δy)2(18)

边界条件方面, FDTD(2, 4)采用完全匹配层(PML)。PML具有两个较为突出的优点:①PML吸收边界条件只需要小部分计算网格就可以实现边缘部分反射波的吸收; ②PML的坐标只需要在吸收边界区域内进行拉伸就可以, 避免在整个计算区域内进行复拉伸。

文中使用的数值稳定条件和PML边界条件已广泛应用于地质雷达的正演模拟中, 故不作详细介绍, 其推理过程可参考文献[15]。

2 数值成像与分析

基于matlab编程软件, 利用FDTD(2, 4)对溶洞、断层、孤石等3种典型不良地质体进行数值模拟。通过分析不同时刻的波场快照图和多偏移距波形剖面图, 总结了电磁波在3种典型地电模型中的响应特征和传播规律。需要强调的是, 本节所有数值模拟基于TE模式下y方向上的电场强度Ey

2.1 溶洞

溶洞模型大小为6 m× 11 m, 溶洞为直径1 m的圆形区域。模型背景的相对介电常数为15, 无水溶洞的相对介电常数设为5, 有水溶洞的相对介电常数为55。发射天线水平距离x=0.5 m, 天线中心频率为100 MHz。

图1为发射天线置于5.5 m, 接收天线水平距离X=5.5 m且每次垂直向上移动0.25 m的波形剖面。横坐标Rx为接收天线的偏移距(下同)。偏移距为3.5~7.5 m段道集的初至波能量较其他道集有明显衰弱, 这是因为电磁波发生了反射。发射天线垂直深度为5.5 m时, 其波场快照显示:由于无水溶洞中介电常数相对背景值小, 电磁波传播快, 因此在t=40 ns和t=60 ns时, 波前呈凸状; 相反, 有水溶洞中介电常数相对背景值大, 电磁波传播慢, 因此在t=40 ns和t=60 ns时, 波前呈凹状。由波场快照图还能清晰地观测到电磁波对溶洞左右边界的反射响应。可见, FDTD(2, 4)的模拟精度高。

图1 溶洞模型及其波场快照

2.2 断层

图2模型为6 m× 11 m的倾斜状断层, 模型背景的相对介电常数为15, 断层的相对介电常数设为25。发射天线水平距离x=0.5 m, 天线中心频率为100 MHz。

图2 溶洞模型波形剖面

发射天线深度为1.5 m时, 波场快照如图3所示:t=55 ns时, 电磁波遇到断层, 产生绕射波并向周围传播; t=75 ns时, 波前面呈凹状, 表明该处传播速度慢; t=115 ns时, 绕射波、反射波向四周传播, 波形较为复杂, 且直达波能量减弱。

图3 断层模型波场快照

图4a为发射天线置于1.5 m, 接收天线x=5.5 m且每次垂直向上移动0.5 m的波形剖面, 后半部波形道集有明显的脱节, 此为倾斜断层导致。在实际勘察工作中, 可根据原始采集波形图推测可疑区, 进而针对可疑区加密射线, 这样可以提高电磁波CT的准确性。

图4 断层模型波形剖面

图4b为针对断层适当加密的覆盖射线, 接收天线每0.25 m移动一次, 可提高层析成像的精度。

因为大偏移距处的电磁波能量衰减大, 并且电磁波在断层中传播速度慢, 所以电磁波波峰出现延迟。图5是接收天线置于9、9.25、9.5 m时的单道电磁波, 峰值出现较迟。在拾取直达波时, 应以第一个起跳点作为初至时(箭头所示)。

图5 接收天线不同位置时的单道波形(箭头为初至)

2.3 孤石

如图6所示, 模型大小为6 m× 11 m, 两个孤石皆为1 m× 3 m, 模型背景的相对介电常数为15, 孤石的相对介电常数设为6。发射天线水平距离x=0.5 m, 天线中心频率为100 MHz。

图6 孤石模型的波长快照

发射天线置于深度5.5 m处时, 波场快照图如图7所示:t=20 ns时, 电磁波到达孤石1的左侧, 并产生了向后传播的反射波; t=60 ns时, 电磁波经过孤石2; t=100 ns时, 电磁波遇到孤石1和孤石2发生的反射波和绕射波向周围传播, 波形复杂。

图7 断层模型波形剖面

图7为发射天线置于深度5.5 m, 接收天线x=5.5 m且每次移动0.25 m的波形剖面。由图7可知, 电磁波在孤石1和孤石2之间发生多次反射。在实际探测工作中, 与孤石地电模型相似的地质情况很常见, 在走时层析成像过程中两个异常体速度差异效应被反映到同一个初至时矩阵, 因此它的电磁波成像结果通常为单个异常区。考虑到电磁波CT横向分辨率不高、难以分辨水平多个异常体的情况, 分析波形剖面图在数据解译中显得尤其重要。

3 结论

FDTD(2, 4)在空间上采用更高阶精度的中心差分近似, 其具有更好的数值色散特性。本文对溶洞、断层和孤石等地电模型进行了数值模拟成像, 由不同时刻的波场快照图可清晰观测到直达波、反射波和绕射波的传播特性以及电磁波能量的衰减情况, 表明了FDTD(2, 4)具有高模拟精度。

结合波场快照图和波形剖面图, 分析、总结了电磁波在不同介质中的响应特征和传播规律, 为多偏移距电磁波CT的初至时拾取以及成像结果的解译提供有利依据。此外, 文中还展现了多偏移距采集模式的高效性和灵活性。

The authors have declared that no competing interests exist.

参考文献
[1] 赵连锋, 朱介寿, 严忠琼, . 岩溶勘察层析成像试验: 成像效果分析[J]. 物探与化探计算技术, 2002, 24(4): 332-336. [本文引用:1]
[2] 陈昌彦, 沈小克, 苏兆锋, . 电磁波层析成像技术在复杂地质边坡工程勘察中的应用研究[J]. 地球物理学进展, 2012, 27(2): 0796-0803. [本文引用:1]
[3] 冯锐. 电磁波层析成象[J]. 地震学报, 1997, 19(5): 524-534. [本文引用:1]
[4] 张曼舳, 师学明. 电磁波层析成像技术进展[J]. 工程地球物理学报, 2009, 6(4): 418-425. [本文引用:1]
[5] 刘昌铨, 孙振国, 李长法. 层析成像技术在核反应堆地基勘察中的应用[J]. 物探与化探, 2000, 24(1): 69-75. [本文引用:1]
[6] YEE K S. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media[J]. IEEE Transactionon Antennasand Propagation, 1966, 14(4): 302. [本文引用:1]
[7] 林树海, 赵立英. 井问电磁波层析成像中的高精度时间域正演计算[J]. 中国地质大学学报: 地球科学, 2007, 32(4): 469-473. [本文引用:1]
[8] Fang J. Time domain infinite difference computation for Maxwell’s equations [D]. Berkeley: Univ California, 1989. [本文引用:1]
[9] Lan K, Liu Y W, Lin W G. A higher order (2, 4) scheme for reducing dispersion in FDTD algorithm[J]. IEEE Trans Electromagnetic Compatibility, 1999, 41(2): 160-165. [本文引用:1]
[10] 吴丰收, 曾昭发, 黄玲, . 探地雷达信号的高阶时间域有限差分模拟[J]. 物探与化探计算技术, 2009, 32(4): 308-313. [本文引用:1]
[11] 孔繁敏, 李康, 郭毅峰. 高阶FDTD法分析电-大尺寸光波导器件[J]. 光电子·激光, 2004, 15(6): 672-6789. [本文引用:1]
[12] 王健. 高阶FDTD方法及其在散射问题的应用[D]. 安徽: 安徽大学, 2007. [本文引用:1]
[13] 曾昭发, 刘四新, 王者江, . 探地雷达方法原理及应用[M]. 北京: 科学出版社, 2006. [本文引用:1]
[14] 冯德山, 戴前伟, 何继善, . 探地雷达GPR正演模拟的时域有限差分实现(英文)[J]. 地球物理学进展, 2006, 02: 630-636. [本文引用:1]
[15] 葛德彪, 阎玉波. 电磁场时域有限差分法[M]. 西安: 西安电子科技大学出版社, 2002. [本文引用:1]