地球物理学报  2020, Vol. 63 Issue (8): 3192-3204   PDF    
基于共形辛Euler算法的非金属地下管道精细化高效探地雷达正演模型
雷建伟1,2,3, 方宏远1,2,3, 李银萍1,2,3, 张娟4, 杨曼1,2,3     
1. 郑州大学水利科学与工程学院, 郑州 450001;
2. 重大基础设施检测修复国家地方联合工程实验室, 郑州 450001;
3. 水利与交通基础设施安全防护河南省协同创新中心, 郑州 450001;
4. 中交第一公路勘察设计研究院有限公司, 西安 710075
摘要:我国城市地下管线种类繁多、规模巨大.近年来,因城市管网家底不清、位置不明导致的管道挖断事故频发.探地雷达(GPR)是目前最为有效的地下设施探测定位设备,通过对探地雷达回波信号进行反演分析,可获取管道埋深、位置等信息.然而,地下管道等地下圆形目标,传统正演模型中阶梯近似方法会产生一定误差.本文结合辛Euler算法与曲面共形技术,建立探地雷达电磁波在含非金属管道地下结构中传播的精细化数值模型,并通过图形处理器(GPU)加速技术进一步提升模型计算效率.数值模拟结果显示,共形网格技术有效降低了由于阶梯近似造成的虚假绕射波,并通过结合GPU并行计算,在不增加网格密度的前提下,实现了对地下管道电磁响应的高效精细化计算.
关键词: 探地雷达      非金属地下管道      辛Euler算法      曲面共形技术      高效数值模型     
GPR forward model of underground non-metallic pipeline based on parallel conformal symplectic Euler algorithm
LEI JianWei1,2,3, FANG HongYuan1,2,3, LI YinPing1,2,3, ZHANG Juan4, YANG Man1,2,3     
1. College of Water Conservancy&Environmental Engineering, Zhengzhou University, Zhengzhou 450001, China;
2. National local joint engineering laboratory of major infrastructure testing and rehabilitation technology, Zhengzhou 450001, China;
3. Collaborative Innovation Center of Water Conservancy and Transportation Infrastructure Safety, Henan Province, Zhengzhou 450001, China;
4. CCCC First Highway Consultants Co., LTD, Xi'an 710075, China
Abstract: There are many types and large scales of urban underground pipelines in China. In recent years, pipeline excavation accidents occur frequently due to unclear urban pipeline network location. Ground Penetrating Radar (GPR) is the most effective tool for detecting and locating underground facilities. The information of pipeline depth and location can be obtained by inverse analysis of GPR echo signal. However, the ladder approximation method in the traditional forward model can produce some errors in the circular target of underground pipelines. Based on the symplectic Euler algorithm and surface conformal technique, a fined numerical model of GPR electromagnetic wave propagation in underground structures containing non-metallic pipeline is established in this paper. The computational efficiency of the model is further improved by Graphics Processing Unit (GPU) acceleration technology. Numerical simulation results show that the surface conformal technique can greatly reduce the pseudo-waves that are generated by the ladder approximation. By combining GPU parallel calculation, the efficient and precise calculation of the electromagnetic response of the underground pipeline is realized without increasing the grid density.
Keywords: Ground penetrating radar    Underground non-metallic pipeline    Symplectic Euler algorithm    Surface conformal technology    Efficient numerical simulation    
0 引言

随着我国城市化进程的加快,城镇地下管网规模不断扩大,供水、排水、燃气、供热、电力、通讯等各类管线密如织网.然而,由于管线权属部门不同,多头管理,施工管理不严等原因,地下管网分布位置信息缺失严重,市政基础设施与建筑基坑施工过程中经常出现挖断既有管网的事故,造成较大经济损失,严重影响居民正常生活.如何快速、准确地确定地下管线的分布位置已经成为城市建设工程中亟待解决的问题.

目前,城市地下管线探测最常用的设备是管线探测仪,但这类设备仅能探测金属管道,不能对非金属管道进行识别.探地雷达是利用电磁波探测地下结构内部物质分布规律的无损检测方法,具有连续、快速、适用范围广等优点,不仅可以应用于地下金属管道探测定位,也可以对地下非金属管道进行检测.通过对探地雷达实测信号进行反演分析,可获取管道位置、埋深等信息,精确高效的地下结构探地雷达电磁波传播数值模型是进行信号反演的前提.

目前,探地雷达电磁波数值模拟方法主要包括:射线跟踪方法(邓世坤和王惠莲,1993),时域伪谱法(Carcione et al., 1999),时域有限元方法(陈晓光等,1999),快速多极方法(Ginste et al., 2004),积分方程方法(Simsek et al., 2006),时域有限差分方法(FDTD)(Roberts and Daniels, 1997Werner et al., 2013刘四新和曾昭发,2007傅磊等,2014冯德山等,2016; 冯德山和王珣,2017),ADI-FDTD方法(冯德山等, 2010, 2014),辛算法(方宏远和林皋,2013Fang et al., 2019a),无单元Galerkin法(冯德山等,2013)等.虽然这些方法均能较好地模拟探地雷达电磁波在地下结构中的传播,但在效率、精度方面仍存在一些不足.例如,FDTD方法受CFL稳定性条件限制(Taflove and Hagness, 2005),计算效率和精度受到制约;辛算法虽然计算效率优于FDTD方法,但仍属差分算法,无法突破CFL稳定性条件限制,计算效率提升有限;ADI-FDTD方法虽然不再受CFL稳定性条件的约束,但其数值色散性随时间步长增加而增大.探地雷达正演模型的精度和效率除了与算法本身性能有关外,还与建模方法以及计算机硬件加速技术有关.

目前提高模拟地下管道等不规则形状目标电磁响应计算精度的方法主要有三种:一是增加离散网格密度,但计算时间也相应增加,对于复杂地下结构,计算资源难以满足;二是采用亚网格技术,即针对结构不同的部分,采用不同的网格划分标准,即对不规则形状目标区域采用细分网格,对于其他区域采用正常网格尺寸.该技术可以有效地提高计算效率和精度,但在不同大小网格跳变的界面,会产生反射波(Xie et al., 2010);三是采用共形网格技术.这种技术又分为直角和曲面坐标共形网格技术.前者在直角坐标系下处理曲边或曲面物体,适用于处理任意形状的物体,后者是在广义正交坐标系下处理有曲面的物体,适用于处理有规则的曲面物体.共形网格技术目前多用于隐身飞机表面覆盖的薄层吸波材料电磁特性研究(胡晓娟等,2006).

此外,并行计算是提高复杂问题计算效率的有效手段,通用图形处理器(GPGPU)技术是近年来广泛应用的一种并行加速技术,为实现探地雷达数值仿真算法提供了良好的并行计算平台.目前GPU加速方面的研究主要集中于地震勘探(刘国峰等,2009刘红伟等,2010),声学(López et al., 2013),生物医学(Chi et al., 2011),计算电磁学(Dziekonski et al., 2011Frances et al., 2015Cicuttin et al., 2017)等.

本文基于辛Euler算法、曲面共形技术和GPU加速技术,建立非金属地下管道结构探地雷达电磁响应的精细化数值模型.首先通过一个圆形空洞模型对比串行非共形和共形辛Euler算法计算结果与并行共形辛Euler算法计算结果,验证并行共形辛Euler算法的高效性与精确性;然后通过对充水和充空气的管道模型模拟研究,分析不同地下管道模型的探地雷达反射图像特征.

1 基本原理 1.1 Hamilton系统和辛分块龙格库塔方法

Hamilton系统对偶方程(Huang and Wu, 2006D'Ambrosio et al., 2011)可以表示为

(1)

式中,H代表Hanmilton函数.可分Hamilton系统中,Hamilton函数表示为如下形式:

(2)

则对偶方程(1)可简化为

(3)

式(3)中的两个偏微分方程分别采用不同的龙格库塔方法进行计算,两种不同的龙格库塔方法具有不同的Butcher表达式:

(4)

将式(4)代入式(3)可以得到如下关系式:

(5)

式中,dt为时间增量,pnqn表示pqn时间步的离散值,QiPi为中间值.

如果(4)式的系数满足如下关系:

(6)

则该分块龙格库塔方法是辛的(Sanz-Serna,1988).

在所有辛算法中,辛Euler算法计算量最小,该算法的Butcher表式如下:

(7)

1.2 控制方程

各向同性有耗介质中,Maxwell方程组表示为

(8)

其中,HE分别代表磁场和电场向量,σεμ分别表示电导率,介电常数和磁导率.

引入矢量磁位H=▽×A,并令E=-U,则有耗体系中广义Hamilton函数可以写为

(9)

利用式(1),Maxwell方程组可以表示为如下正则方程形式:

(10)

对于二维TM波情形,则式(10)可化简为

(11)

式中,UzAz为场组分UA沿z方向的分量,▽2为拉普拉斯算子,文中采用二阶中心差分对拉普拉斯算符进行离散.

TM波磁场和电场分量HxHyEz可按下式得到

(12)

FDTD方法采用Yee网格离散Maxwell旋度方程组中的偏微分算子,其中电场和磁场在时间上和空间上彼此间隔半个时间和空间步长,即电场位于元胞中央,磁场环绕在电场周围,如图 1a所示.辛Euler算法则是电场和磁场离散定义在相同的空间网格节点和相同的时间步长上,如图 1b所示,这种网格离散方式更适合处理复杂边界问题(Fang and Lin, 2012; Fang et al., 2013, 2019b).

图 1 二维TM波场组分空间分布 (a) FDTD方法场组分空间分布; (b)辛分块Runge-Kutta方法场组分空间分布. Fig. 1 Sketch map of spatial distribution of field components in 2D TM case (a) Sketch map of field components in FDTD method; (b) Sketch map of field components in SPRK method.

将(7)式应用到(11)式,可得如下辛Euler方法迭代式:

(13)

化简后可得一阶辛Euler方法的迭代格式:

(14)

式中,Ai, jnUi, jn分别表示场组分AzUzndt时刻,空间网格节点(i, j)上的离散值.

1.3 边界条件

探地雷达电磁波在地下结构中的传播是开域问题,因此在计算区域截断边界处需设置合理的吸收边界条件.本文采用吸收效果较好的坐标伸缩完全匹配层(CPML)边界(Berenger,1998Roden and Gedney, 2000).在有损耗介质FDTD算法中,二维TM波的CPML频域方程可以写为(张鲁新等,2010余翔等,2017)

(15)

将公式(12)代入到公式(15)得到辛Euler算法二维TM波的CPML频域方程如下:

(16)

式中si是坐标伸缩因子,定义如下:

(17)

式中κi为有效延伸因子,αi为场值在i方向的坐标延伸因子的自由度,σi为CPML边界区域中的电导率.

将公式(16)第一式从频域转换到时域,转换到时域后公式右边会存在卷积形式:

(18)

式中的sxsy分别为1/sx和1/sy的逆拉普拉斯变换,公式右侧的卷积使用递归卷积进行运算.

根据拉普拉斯变换理论,si的脉冲响应公式如下:

(19)

式中u(t)和δ(t)分别为单位阶跃函数和单位冲击函数.

将公式(19)代入到公式(18)得到:

(20)

为了提高卷积计算效率,将ζi(t)离散冲击响应定义为

(21)

式中ai定义如下:

(22)

本文采用一阶向前差商的向后差商来近似二阶微商,即用

(23)

来近似同理.

根据公式(21)和(22),按照辛分块Runge-Kutta方法将式(20)进行空间和时间离散得

(24)

式中的离散卷积计算过程较为复杂,由于Z0i(m)是简单的指数形式,从而它们的和可以用递归卷积递归得到(方思南等,2016).引入一组新的辅助表达式ψi,式(24)改写为

(25)

式中ψzxn+1(i, j)和ψzyn+1(i, j)可由如下公式计算

(26)

(27)

式中bi定义如下:

(28)

公式(25)经整理后得到

(29)

式中

其余场量可类似推出.CPML的PML层内的各参数κασ单调变化,以x方向为例,通常取

(30)

式中d为PML层厚度,x0为PML与模拟区域分界处的网格编号.本文模拟中取d=10,m=4,κmax=5,αmax=0.008,,其中δ模拟网格大小.

2 共形辛Euler算法及其GPU并行计算 2.1 曲面共形技术

本文采用一种基于有效介质参数的共形网格方法来模拟地下管道圆形边界.图 2a显示了圆管的实际细分网格,白色网格是普通网格点,绿色网格是圆管的实际细分网格.图 2b显示了一个圆管实际模拟计算的常规阶梯近似的细分网格,橙色网格是圆管的常规阶梯近似网格.图 2c显示了圆管的共形网格方法的细分网格,紫色网格是共形网格,其他网格为非共形网格.

图 2 圆管的细分网格 (a)圆管的实际细分网格; (b)圆管的常规阶梯近似的细分网格; (c)圆管的共形网格细分. Fig. 2 Schematic view for the subdivisions of a circle pipeline (a) An actual model for a circle pipeline; (b) The subdivision diagram using the conventional ladder approximation; (c) The conformal grid for a circle pipeline.

下面从图 2c中提取共形网格单元,说明二维TM波中共形网格点的等效介质参数关系.在辛Euler算法中,场分量AU是在相同的时间步长和相同的空间网格节点上定义的.在直角坐标下介质共形网格点的计算如图 3所示,其中F是场分量AU的取样点,Δx和Δy是网格的宽度和高度,Sxy1Sxy2分别是介质1和2区域的面积.假设介质1和2的电磁参数分别为ε1σ1μ1ε2σ2μ2.场分量UA位于网格单元的中心.介电常数、导电率和磁导率系数的有效值可以通过对不同介质所占网格区域面积的加权平均得到.

图 3 共形网格点介质参数等效示意 (a)实际的细分网格; (b)共形网格. Fig. 3 Schematic view for the equivalence of the medium parameters for a conformal grid point (a) An actual subdivision grid; (b) A conformal grid.

图 3b所示,场分量采样点F处的等效介电常数、导电率和磁导率系数如下:

(31)

式中eff表示参数的等效值.

将方程(31)代入方程(14),得到共形辛Euler算法的迭代微分方程:

(32)

2.2 共形辛Euler算法的GPU并行计算

CUDA(Compute Unified Device Architecture)能够利用英伟达GPU并行计算引擎高效解决许多复杂计算任务(NVIDIA Corporation,2017),非常适合解决并行计算问题(Luebke and Humphreys, 2007Lindholm et al., 2008李博等,2009Demir and Elsherbeni, 2010陈召曦等,2012石颖等,2013).在GPU上实现的并行共形辛算法与传统在CPU上实现的串行辛算法方法的区别在于并行共形辛Euler算法计算部分在不同的设备上执行,并具有相似的执行过程.二维辛算法问题最初被划分为粗糙的子问题,如系统初始化和建模、磁场和电场的更新和数据输出.在系统初始化和建模中,字段值设置为零.通过读取输入文件来分配系统参数,例如正演数值模型和激励源.数据输出子问题是程序的串行部分,在CPU上执行.在GPU上,场分量AU的更新是相互独立的,U的更新在A更新之后.这两个计算密集型的子问题(内核函数)是在GPU上执行的.GPU中实现二维并行共形辛Euler算法的流程如图 4所示.

图 4 并行共形辛算法程序流程 Fig. 4 The flowchart of parallel conformal symplectic algorithm

并行辛算法中场分量AU的计算是在二维xy方域内进行的, 线程的数量由计算域决定,也就是辛分块单元的数量.在处理的二维共形辛算法问题时,需两个CUDA内核函数.第一个内核函数计算A,第二个更新U.如图 5所示,采用的策略是每个线程(thread)计算一个辛分块单元格,每个线程块(block)计算一组连续的辛分块单元格,按照以上策略对整个计算域进行并行计算.

图 5 二维并行共形辛算法的线程安排方式 Fig. 5 Thread arrangement of two-dimensional parallel conformal symplectic algorithm

采用二维线源传播模型验证并行算法的准确性以及CPML边界条件的吸收效果.开发工具是Visual Studio 2010和CUDA toolkit 7.5,中央处理器为Intel Core i7-6700K搭载NVIDIA GeForce GTX 1070.模拟区域为200 cm×200 cm的正方形,内部介质为空气,在模拟区域的中心加入一个中心频率为1 GHz单位振幅Ricker子波(图 6).空间步长和时间步长分别取0.5 cm和0.01 ns.图 7为不同时刻波场的快照图,由图可知CPML吸收边界条件吸收效果良好.

图 6 中心频率为1 GHz的单位波幅Ricker子波波形 Fig. 6 A Ricker waveform of unity amplitude and central frequency of 1 GHz
图 7 不同时刻波场快照 (a) T=200 ns; (b) T=300 ns; (c) T=400 ns; (d) T=500 ns. Fig. 7 Wave field snapshots at different times
3 数值算例 3.1 含有圆形空洞的两层介质模型

首先,通过含有圆形空洞的两层介质模型验证算法的精度与效率.如图 8所示,该模型的第一层为空气层,相对介电常数εr=1,电导率σ=0;第二层介质代表沙层,相对介电常数εr=12,电导率σ=2 mS·m-1;第二层中存在一个半径为5cm的圆形空洞,其内部介质相对介电常数εr=30,电导率σ=0,假定所有材料的相对磁导率μ都为1.激励源选取中心频率为1 GHz的单位振幅Ricker子波,发射器Tx和接收器Rx之间的间隔为10 cm.

图 8 含有圆形空洞的两层介质模型 Fig. 8 Sketch map of two-layer dielectric model with circular void

为了验证并行共形辛Euler算法的精确性和高效性,分别采用并行共形、并行非共形以及串行非共形辛Euler算法模拟该模型,在并行非共形辛Euler算法模拟中,用不同的网格剖分介质并进行模拟.模拟得到单道雷达数据对比图如图 9所示,并行非共形辛Euler算法模拟中不同网格模拟时间如表 1所示.

图 9 不同算法模拟得到的单道雷达数据对比 Fig. 9 Comparison diagram of single radar data simulated by different algorithms
表 1 不同网格的模拟时间 Table 1 Simulation time of different grid

图 9表 1可知,并行非共形和并行共形辛Euler算法在模拟网格都为800×800时,并行非共行辛Euler算法所需时间相对较少,但其计算精度较低.通过细化网格可提升其模拟精度,当网格为4000×4000时,其精度接近并行共形辛Euler算法的精度,但其所需时间为202.258 s,耗时是并行共形辛Euler算法的17.12倍.

图 10-12分别为串行非共形辛Euler算法、并行非共形辛Euler算法以及并行共形辛Euler算法获得的GPR剖面图.图 10中的串行非共形辛Euler算法需要5674.89 s,图 11中的并行非共形辛Euler算法需要257.635 s,图 12中并行共形辛Euler算法需要316.797 s.与串行非共形辛Euler算法相比,基于GPU实现的并行共形辛Euler算法节约94.42%的计算时间.由图 9-12以及表 1可知,在圆形空洞模拟中,并行共形辛Euler算法在提高计算效率的同时能够显著地减少由常规阶梯近似引起的虚假反射波,提高模拟的精确性.

图 10 串行非共形辛Euler算法得到的GPR剖面 Fig. 10 GPR profile obtained by non-conformal symplectic Euler algorithm
图 11 并行非共形辛Euler算法得到的GPR剖面 Fig. 11 GPR profile obtained by parallel non-conformal symplectic Euler algorithm
图 12 并行共形辛Euler算法得到的GPR剖面 Fig. 12 GPR profile obtained by parallel conformal symplectic Euler algorithm
3.2 地下管道模型

本模型为1:1的实际地下管道模型,如图 13所示.该模型的第一层为空气层,相对介电常数εr=1,电导率σ=0;第二层介质代表土层,相对介电常数εr=12,电导率σ=2 mS·m-1;第二层中存在一个直径为100 cm的水泥混凝土管,管壁厚度为10 cm,其材料的相对介电常数εr=6,电导率σ=0.所有材料都假定为非磁性的,即相对磁导率都为1.空间步长和时间步长分别取0.5 cm和0.01 ns,模拟时间步为5000步.本例中,分别模拟管内充水和管内为空气的情况,其中管内为空气时,相对介电常数εr=1,电导率σ=0;其中管内充水时,其对介电常数εr=81,电导率σ=0.激励源选取中心频率为1 GHz的单位波幅Ricker子波,从x=10 cm起,每隔10 cm发射器Tx和接收器Rx发射接收一次.

图 13 地下管道模型 Fig. 13 Schematic map of underground pipe model

图 14图 15分别为并行共形与非共形辛Euler算法获得的管内充满空气和管内充水的地下管道模型的单道雷达数据对比图.图 1617为串行非共形与共形以及并行共形辛Euler算法获得管内充满空气和管内充水的地下管道模型GPR剖面图.图 16a中串行非共形辛Euler算法需要110047.52 s,图 16b中串行共形辛Euler算法需要117860.57 s,图 16c中并行共形辛Euler算法获得管内充满空气的地下管道模型GPR剖面图需要8163.76 s,与串行非共形辛Euler算法相比,基于GPU实现的并行共形辛Euler算法节约大约92.58%的计算时间.图 17a中串行非共形辛Euler算法需要115152.72 s,图 17b中串行共形辛Euler算法需要119274.63 s,图 17c中并行共形辛Euler算法获得管内充水的地下管道模型GPR剖面图需要8212.88 s,与串行非共形辛Euler算法相比,基于GPU实现的并行共形辛Euler算法节约大约92.87%的计算时间.

图 14 并行共形辛Euler算法获得的地下管道模型单道雷达数据对比 Fig. 14 Comparison diagram of single channel GPR data of underground pipeline model obtained by parallel conformal symplectic Euler algorithm
图 15 并行非共形辛Euler算法获得的地下管道模型单道雷达数据对比 Fig. 15 Comparison diagram of single channel GPR data of underground pipeline model obtained by parallel non-conformal symplectic Euler algorithm
图 16 管内充满空气地下管道模型GPR剖面 (a)串行非共形辛Euler算法得到的GPR剖面;(b)串行共形辛Euler算法得到的GPR剖面;(c)并行共形辛Euler算法得到的GPR剖面. Fig. 16 GPR profile of underground pipeline model filled with air in pipe (a) GPR profile obtained by serial non-conformal symplectic Euler algorithm; (b) GPR profile obtained by serial conformal symplectic Euler algorithm; (c) GPR profile obtained by parallel conformal symplectic Euler algorithm.
图 17 管内充水的地下管道模型GPR剖面 (a)串行非共形辛Euler算法得到的GPR剖面;(b)串行共形辛Euler算法得到的GPR剖面;(c)并行共形辛Euler算法得到的GPR剖面. Fig. 17 GPR profile of underground pipeline model filled with water in pipe (a) GPR profile obtained by serial non-conformal symplectic Euler algorithm; (b) GPR profile obtained by serial conformal symplectic Euler algorithm; (c) GPR profile obtained by parallel conformal symplectic Euler algorithm.

图 14-17可得,不同状况即管内为空或是管内充水的地下管道模型的GPR反射图像特征差别很大.管内为空气的地下管道模型的GPR反射图像中有特别明显的绕射双曲线特征.通过分析探地雷达数据中的绕射双曲线形状特征以及时间可以精确地定位管道位置以及管内状况.

4 结论

本文结合辛Euler算法,曲面共形技术及GPU加速方法,建立了探地雷达电磁波在含非金属管道地下结构中传播的精细化数值模型,实现了对非金属地下管道电磁响应的精确高效计算.模拟结果表明,结合曲面共形网格技术的辛Euler算法,能大大减少正演数值计算中由常规阶梯近似圆形边界引起的误差,有效降低了虚假绕射波的产生;此外,通过GPU加速技术大幅提高了模型计算效率.与串行辛Euler算法相比,基于GPU的并行辛Euler算法可节约计算时间92%以上.该模型可为下一步反演计算提供精确高效的正演数值模型.

References
Berenger J P. 1998. An effective PML for the absorption of evanescent waves in waveguides. IEEE Microwave and Guided Wave Letters, 8(5): 188-190. DOI:10.1109/75.668706
Carcione J M, Lenzi G, Valle S. 1999. GPR modelling by the Fourier method:improvement of the algorithm. Geophysical Prospecting, 47(6): 1015-1029. DOI:10.1046/j.1365-2478.1999.00151.x
Chen X G, Jin Y Q, Nie Z P. 1999. A two-dimensional finite element time-domain method for electromagnetic responses in inhomogeneous media. Acta Geophysica Sinica (in Chinese), 42(2): 268-276.
Chen Z X, Meng X H, Guo L H, et al. 2012. Three-dimensional fast forward modeling and the inversion strategy for large scale gravity and gravimetry data based on GPU. Chinese Journal of Geophysics (in Chinese), 55(12): 4069-4077. DOI:10.6038/j.issn.0001-5733.2012.12.019
Chi J R, Liu F, Weber E, et al. 2011. GPU-accelerated FDTD modeling of radio-frequency field-tissue interactions in high-field MRI. IEEE Transactions on Biomedical Engineering, 58(6): 1789-1796. DOI:10.1109/TBME.2011.2116020
Cicuttin M, Codecasa L, Kapidani B, et al. 2017. GPU accelerated time-domain discrete geometric approach method for Maxwell's equations on tetrahedral grids. IEEE Transactions on Magnetics, 54(3): 7203004.
D'Ambrosio R, Ixaru L G, Paternoster B. 2011. Construction of the ef-based Runge-Kutta methods revisited. Computer Physics Communications, 182(2): 322-329. DOI:10.1016/j.cpc.2010.10.009
Demir V, Elsherbeni A Z. 2010. Compute unified device architecture (CUDA) based finite-difference time-domain (FDTD) implementation. Applied Computational Electromagnetics Society Journal, 25(4): 303-314.
Deng S Q, Wang H L. 1993. The forward synthesizing and migration processing for GPR image. Acta Geophysica Sinica (in Chinese), 36(4): 528-536.
Dziekonski A, Lamecki A, Mrozowski M. 2011. GPU acceleration of multilevel solvers for analysis of microwave components with finite element method. IEEE Microwave and Wireless Components Letters, 21(1): 1-3. DOI:10.1109/LMWC.2010.2089974
Fang H Y, Lin G. 2012. Symplectic partitioned Runge-Kutta methods for two-dimensional numerical model of ground penetrating radar. Computers & Geosciences, 49: 323-329.
Fang H Y, Lin G, Zhang R L. 2013. The first-order symplectic euler method for simulation of GPR wave propagation in pavement structure. IEEE Transactions on Geoscience and Remote Sensing, 51(1): 93-98. DOI:10.1109/TGRS.2012.2202121
Fang H Y, Lin G. 2013. Simulation of GPR wave propagation in complicated geoelectric model using symplectic method. Chinese Journal of Geophysics (in Chinese), 56(2): 653-659. DOI:10.6038/cjg20130229
Fang H Y, Lei J W, Yang M, Li Z W. 2019a. Analysis of GPR Wave Propagation Using CUDA-Implemented Conformal Symplectic Partitioned Runge-Kutta Method. Complexity.
Fang H Y, Lei J W, Zhang J, et al. 2019b. Modelling ground-penetrating radar wave propagation using graphics processor unit parallel implementation of the symplectic Euler method. Near Surface Geophysics, 17(4): 417-425.
Fang S N, Pan H P, Du T, et al. 2016. Three-dimensional cross-well electromagnetic modeling considering numerical dispersion in convolutional perfectly matched layers. Chinese Journal of Geophysics (in Chinese), 59(5): 1888-1897. DOI:10.6038/cjg20160531
Feng D S, Chen C S, Dai Q W. 2010. GPR numerical simulation of full wave field based on UPML boundary condition of ADI-FDTD. Chinese Journal of Geophysics (in Chinese), 53(10): 2484-2496. DOI:10.3969/j.issn.0001-5733.2010.10.022
Feng D S, Wang H H, Dai Q W. 2013. Forward simulation of Ground Penetrating Radar based on the element-free Galerkin method. Chinese Journal of Geophysics (in Chinese), 56(1): 298-308. DOI:10.6038/cjg20130131
Feng D S, Chen J W, Wu Q. 2014. A hybrid ADI-FDTD subgridding scheme for efficient GPR simulation of dispersion media. Chinese Journal of Geophysics (in Chinese), 57(4): 1322-1334. DOI:10.6038/cjg20140429
Feng D S, Yang L Y, Wang X. 2016. The unsplit convolutional perfectly matched layer absorption performance analysis of evanescent wave in GPR FDTD forward modeling. Chinese Journal of Geophysics (in Chinese), 59(12): 4733-4746. DOI:10.6038/cjg20161232
Feng D S, Wang X. 2017. Convolution perfectly matched layer for the finite-element time-domain method modeling of ground penetrating radar. Chinese Journal of Geophysics (in Chinese), 60(1): 413-423. DOI:10.6038/cjg20170134
Frances J, Otero B, Bleda S, et al. 2015. Multi-GPU and multi-CPU accelerated FDTD scheme for vibroacoustic applications. Computer Physics Communications, 191: 43-51. DOI:10.1016/j.cpc.2015.01.017
Fu L, Liu S X, Liu L B, et al. 2014. Airborne ground penetrating radar numerical simulation and reverse time migration. Chinese Journal of Geophysics (in Chinese), 57(5): 1636-1646. DOI:10.6038/cjg20140526
Ginste D V, Rogier H, Olyslager F, et al. 2004. A fast multipole method for layered media based on the application of perfectly matched layers-The 2-D case. IEEE Transactions on Antennas and Propagation, 52(10): 2631-2640. DOI:10.1109/TAP.2004.834396
Hu X J, Ge D B, Wei B. 2006. Study on MCFDTD for 3-D coated targets by using effective parameters. Systems Engineering and Electronics (in Chinese), 28(11): 1652-1667.
Huang Z X, Wu X L. 2006. Symplectic partitioned Runge-Kutta scheme for Maxwell's equations. International Journal of Quantum Chemistry, 106(4): 839-842. DOI:10.1002/qua.20852
López J J, Carnicero D, Ferrando N, et al. 2013. Parallelization of the finite-difference time-domain method for room acoustics modelling based on CUDA. Mathematical and Computer Modelling, 57(7-8): 1822-1831. DOI:10.1016/j.mcm.2011.11.075
Li B, Liu G F, Liu H. 2009. A method of using GPU to accelerate seismic pre-stack time migration. Chinese Journal of Geophysics (in Chinese), 52(1): 245-252.
Lindholm E, Nickolls J, Oberman S., et al. 2008. NVIDIA Tesla:A unified graphics and computing architecture. IEEE Micro, 28(2): 39-55.
Liu G F, Liu H, Li B, et al. 2009. Method of prestack time migration of seismic data of mountainous regions and its GPU implementation. Chinese Journal of Geophysics (in Chinese), 52(12): 3101-3108. DOI:10.3969/j.issn.0001-5733.2009.12.019
Liu H W, Li B, Liu H, et al. 2010. The algorithm of high order finite difference pre-stack reverse time migration and GPU implementation. Chinese Journal of Geophysics (in Chinese), 53(7): 1725-1733. DOI:10.3969/j.issn.0001-5733.2010.07.024
Liu S X, Zeng Z F. 2007. Numerical simulation for Ground Penetrating Radar wave propagation in the dispersive medium. Chinese Journal of Geophysics (in Chinese), 50(1): 320-326.
Luebke D, Humphreys G. 2007. How GPUs work. Computer, 40(2): 96-100.
NVIDIA Corporation. 2017. NVIDIA CUDA C Programming Guide, version 9.1 edition.https://docs.nvidia.com/cuda/archive/8.0/pdf/CUDA_C_Programming_Guide.pdf
Roberts R L, Daniels J J. 1997. Modeling near-field GPR in three dimensions using the FDTD method. Geophysics, 62(4): 1114-1126.
Roden J A, Gedney S D. 2000. Convolution PML (CPML):an efficient FDTD implementation of the CFS-PML for arbitrary media. Microwave and Optical Technology Letters, 27(5): 334-339. DOI:10.1002/1098-2760(20001205)27:5<334::AID-MOP14>3.0.CO;2-A
Sanz-Serna J M. 1988. Runge-Kutta schemes for Hamiltonian systems. BIT Numerical Mathematics, 28(4): 877-883.
Shi Y, Wang W H, Li Y, et al. 2013. 3D surface-related multiple prediction approach investigation based on wave equation. Chinese Journal of Geophysics (in Chinese), 56(6): 2023-2032. DOI:10.6038/cjg20130623
Simsek E, Liu J G, Liu Q H. 2006. A spectral integral method and hybrid SIM/FEM for layered media. IEEE Transactions on Microwave Theory and Techniques, 54(11): 3878-3884. DOI:10.1109/TMTT.2006.883647
Taflove A, Hagness S C. 2005. Computational Electrodynamics:The Finite-Difference Time-Domain Method. Boston: Artech House.
Werner G R, Bauer C A, Cary J R. 2013. A more accurate, stable, FDTD algorithm for electromagnetics in anisotropic dielectrics. Journal of Computational Physics, 255: 436-455. DOI:10.1016/j.jcp.2013.08.009
Xie J, Yang Y, Liu S B, et al. 2010. Huygens subgrid technique for FDTD analysis of plasma problems.//2010 International Conference on Microwave and Millimeter Wave Technology(ICMMT). Chengdu, China: IEEE, 929-932.
Yu X, Wang X B, Li X J, et al. 2017. Three-dimensional finite difference forward modeling of the transient electromagnetic method in the time domain. Chinese Journal of Geophysics (in Chinese), 60(2): 810-819. DOI:10.6038/cjg20170231
Zhang L X, Fu L Y, Pei Z L. 2010. Finite difference modeling of Biot's poroelastic equations with unsplit convolutional PML and rotated staggered grid. Chinese Journal of Geophysics (in Chinese), 53(10): 2470-2483. DOI:10.3969/j.issn.0001-5733.2010.10.021
陈晓光, 金亚秋, 聂在平. 1999. 轴对称有耗介质电磁问题的时域有限元方法. 地球物理学报, 42(2): 268-276. DOI:10.3321/j.issn:0001-5733.1999.02.015
陈召曦, 孟小红, 郭良辉, 等. 2012. 基于GPU并行的重力、重力梯度三维正演快速计算及反演策略. 地球物理学报, 55(12): 4069-4077. DOI:10.6038/j.issn.0001-5733.2012.12.019
邓世坤, 王惠濂. 1993. 探地雷达图像的正演合成与偏移处理. 地球物理学报, 36(4): 528-536. DOI:10.3321/j.issn:0001-5733.1993.04.012
方宏远, 林皋. 2013. 基于辛算法模拟探地雷达在复杂地电模型中的传播. 地球物理学报, 56(2): 653-659. DOI:10.6038/cjg20130229
方思南, 潘和平, 杜婷, 等. 2016. 考虑卷积完全匹配层数值色散的井间电磁三维正演. 地球物理学报, 59(5): 1888-1897. DOI:10.6038/cjg20160531
冯德山, 陈承申, 戴前伟. 2010. 基于UPML边界条件的交替方向隐式有限差分法GPR全波场数值模拟. 地球物理学报, 53(10): 2484-2496. DOI:10.3969/j.issn.0001-5733.2010.10.022
冯德山, 王洪华, 戴前伟. 2013. 基于无单元Galerkin法探地雷达正演模拟. 地球物理学报, 56(1): 298-308. DOI:10.6038/cjg20130131
冯德山, 陈佳维, 吴奇. 2014. 混合ADI-FDTD亚网格技术在探地雷达频散媒质中的高效正演. 地球物理学报, 57(4): 1322-1334. DOI:10.6038/cjg20140429
冯德山, 杨良勇, 王珣. 2016. 探地雷达FDTD数值模拟中不分裂卷积完全匹配层对倏逝波的吸收效果研究. 地球物理学报, 59(12): 4733-4746. DOI:10.6038/cjg20161232
冯德山, 王珣. 2017. 基于卷积完全匹配层的非规则网格时域有限元探地雷达数值模拟. 地球物理学报, 60(1): 413-423. DOI:10.6038/cjg20170134
傅磊, 刘四新, 刘澜波, 等. 2014. 机载探地雷达数值模拟及逆时偏移成像. 地球物理学报, 57(5): 1636-1646. DOI:10.6038/cjg20140526
胡晓娟, 葛德彪, 魏兵. 2006. 基于等效参数的三维涂层目标MCFDTD分析. 系统工程与电子技术, 28(11): 1652-1667. DOI:10.3321/j.issn:1001-506X.2006.11.011
李博, 刘国峰, 刘洪. 2009. 地震叠前时间偏移的一种图形处理器提速实现方法. 地球物理学报, 52(1): 245-252.
刘国峰, 刘洪, 李博, 等. 2009. 山地地震资料叠前时间偏移方法及其GPU实现. 地球物理学报, 52(12): 3101-3108. DOI:10.3969/j.issn.0001-5733.2009.12.019
刘红伟, 李博, 刘洪, 等. 2010. 地震叠前逆时偏移高阶有限差分算法及GPU实现. 地球物理学报, 53(7): 1725-1733. DOI:10.3969/j.issn.0001-5733.2010.07.024
刘四新, 曾昭发. 2007. 频散介质中地质雷达波传播的数值模拟. 地球物理学报, 50(1): 320-326. DOI:10.3321/j.issn:0001-5733.2007.01.040
石颖, 王维红, 李莹, 等. 2013. 基于波动方程三维表面多次波预测方法研究. 地球物理学报, 56(6): 2023-2032. DOI:10.6038/cjg20130623
余翔, 王绪本, 李新均, 等. 2017. 时域瞬变电磁法三维有限差分正演技术研究. 地球物理学报, 60(2): 810-819. DOI:10.6038/cjg20170231
张鲁新, 符力耘, 裴正林. 2010. 不分裂卷积完全匹配层与旋转交错网格有限差分在孔隙弹性介质模拟中的应用. 地球物理学报, 53(10): 2470-2483. DOI:10.3969/j.issn.0001-5733.2010.10.021