# Analysis of GPR Wave Propagation Using CUDA-Implemented Conformal Symplectic Partitioned Runge-Kutta Method

Authors: Hongyuan Fang; Jianwei Lei; Man Yang; Ziwei Li
Journal: Complexity
 (2019)

---

## Abstract

Accurate forward modeling is of great significance for improving the accuracy and speed of inversion. For forward modeling of large sizes and fine structures, numerical accuracy and computational efficiency are not high, due to the stability conditions and the dense grid number. In this paper, the symplectic partitioned Runge-Kutta (SPRK) method, surface conformal technique, and graphics processor unit (GPU) acceleration technique are combined to establish a precise and efficient numerical model of electromagnetic wave propagation in complex geoelectric structures, with the goal of realizing a refined and efficient calculation of the electromagnetic response of an arbitrarily shaped underground target. The results show that the accuracy and efficiency of ground-penetrating radar (GPR) forward modeling are greatly improved when using our algorithm. This provides a theoretical basis for accurately interpreting GPR detection data and accurate and efficient forward modeling for the next step of inversion imaging.


---

## Body

## 1. Introduction

Ground-penetrating radar (GPR) is a nondestructive testing method that uses the penetrating ability of electromagnetic waves to determine the distribution of substances in an underground structure [1–5]. Forward modeling of GPR is used to simulate and analyze the electromagnetic response of an underground target by a numerical method, in order to obtain the reflected wave-shape of an underground target on the surface and to understand the propagation laws of the electromagnetic wave in an underground structure [6–8]. Moreover, it can deepen the understanding of a GPR-measured profile, and thus improve the processing and interpretation accuracy of the GPR detection data [9–13]. In addition, forward modeling is still the foundation for determining the inversion of underground target structural parameters based on GPR-measured signals. Accurate and highly effective forward modeling is of great significance for improving the inversion accuracy and speed [14, 15].At present, the forward modeling methods of GPR mainly include the finite element method (FEM) [16, 17], ray tracing method [18, 19], the pseudo-spectral time domain method (PSTD) [20], the fast multipole method [21], the finite difference time domain method (FDTD) [22–24], the alternating direction-implicit FDTD (ADI-FDTD) [25], and the symplectic Euler algorithm [26, 27]. Although these methods can accurately simulate the propagation of a GPR electromagnetic wave in underground structures, these algorithms have some shortcomings with respect to efficiency and precision. For example, the FEM is suitable for solving problems with complex boundary conditions or fixed solution problems with complex media, but its computational program is complex, and the solution time is long, which may result in the phenomenon of a pseudo-solution. The PSTD method adopts the fast Fourier transform (FFT) and the inverse fast Fourier transform (IFFT) to calculate the spatial difference quotient, which theoretically has infinite-order accuracy. However, the calculation accuracy is low when simulating metal object targets and cases with discontinuous medium distribution problems because of the inherent periodicity of FFT. The FDTD method is the most widely used time domain electromagnetic field numerical simulation technology, due to its simple implementation, wide application, and suitability for parallel operation. However, the computational efficiency of this method is limited by the Courant–Friedrichs–Lewy (CFL) stability condition [28], and it is not efficient for large-sized or fine-structure targets. For two-dimensional problems, the symplectic Euler algorithm can describe the electromagnetic field by using only two equations. The computational efficiency and computer occupancy rate of the symplectic algorithm are better than those of the FDTD method; however it is still a difference algorithm and thus is limited by the CFL stability condition. Although the ADI-FDTD method is free from the constraints of CFL stability, the numerical dispersion of the algorithm is increased with a large time step.The accuracy and efficiency of a GPR forward model are not only related to the performance of the algorithm itself, but also related to the modeling method and computer hardware acceleration technology. There are three main methods for improving the accuracy when determining the electromagnetic response of a simulated irregular shape target: one is to increase the density of the discrete grids; however the computational memory and time also increase correspondingly. For complex problems, the computational resources are difficult to meet. Another approach is to adopt subgrid technology [29], such that the irregular shape target area uses a subdivision grid, and a normal grid is adopted for other areas. The subgrid technology can greatly improve the calculation accuracy and efficiency, but it is also easy to generate reflection waves at the interface of the coarse and fine grids, and the calculation is also unstable. A last approach is conformal grid technology [30]. This technology is further divided into surface coordinate and Cartesian coordinate conformal grid technologies. The former is suitable for dealing with objects with regular surfaces, such as balls and columns. The latter can deal with objects with arbitrary surfaces, and the formula deductions and calculations are simple. Conformal grid technology is currently used to study the electromagnetic characteristics of the thin-layer absorbing materials for stealth aircraft surface coverage. The application of conformal grid technology has not been developed in accurate simulations of complex underground structural boundaries.In a two-dimensional simulation of complex irregular shape targets, the grid that is used is very dense, and the time and memory needed for the simulation calculation are accordingly increased. It is difficult to meet these computing requirements when using traditional central processing unit (CPU) operations. Parallel computing is an effective method to improve computational efficiency. General purpose graphics processing unit (GPGPU) technology has higher computing power and lower equipment cost [31–33], which provides a good parallel computing platform for the numerical simulation of GPR with complex underground structures. At present, research on GPU acceleration is mainly focused on acoustics [34], biomedical applications [35], computational electromagnetics [36], and other applications. A parallel symplectic algorithm based on GPU has not been developed.Using GPU acceleration techniques, a two-dimensional numerical model of GPR electromagnetic wave propagation, based on a parallel conformal SPRK algorithm, is established in this paper. A refined numerical simulation of the electromagnetic response for any complex underground structure is realized by the means of curved conformal grid technology. This provides a theoretical basis for the accurate interpretation of GPR detection data and provides an efficient and accurate forward model for the next step of GPR two-dimensional inversion imaging.In this paper, the simulation model of GPR wave propagation in two-dimensional subsurface structure is developed by using conformal and nonconformal SPRK methods extended into a programmable GPU-based parallel-computing platform. A line source is considered as the incident wave and the transmitting boundary condition is adopted to avoid boundary reflections. We verified the accuracy and efficiency of the proposed algorithm by three numerical simulations. One simulated the plane wave reflection from an embankment model, and the other simulated the GPR trace profile of a structural damage pavement model, the last one simulated GPR wave propagation in the underground pipes model. We compared the reflection waveforms simulated by CUDA-implemented conformal and nonconformal SPRK methods and the time simulated by CUDA-implemented and traditional conformal symplectic partitioned Runge-Kutta methods. The results indicate that the proposed CUDA-implemented method requires approximately one-tenth the time of the traditional method and the CUDA-implemented conformal SPRK algorithm can greatly reduce the pseudo-waves that are generated by the ladder approximation, in comparison to the nonconformal SPRK algorithm.This paper is organized as follows: Section2 describes the Hamiltonian Systems and the SPRK methods, the Curved Conformal Grid Technology, and the GPU implementation. Section 3 is devoted to presenting the numerical simulations results analysis. Conclusions and future research are provided in Section 4.


## 2. Methodology

### 2.1. Hamiltonian Systems and the Symplectic Partitioned Runge-Kutta Methods [37–39]

A special Hamilton system can be expressed as(1)dpdt=-∂H∂f=pf,gdqdt=∂H∂g=qf,g,where H(p,q) represents the Hamiltonian function, f and g are independent variables, and p and q are dependent variables, such as E and H. The above equations are often solved by the s-stage-partitioned Runge-Kutta methods, which can be specified by the following Butecher table [40]:(2)Applying (2) in (1) yields the Hamiltonian system relations, which are formulated as in the following equations:(3)Pi=pfn+τ∑j=1saijpQjQi=qgn+τ∑j=1sAijqPjfn+1=fn+τ∑i=1sbifQign+1=gn+τ∑i=1sBigPi,i=1,2,…,s,where Qi and Pi are the internal stages corresponding to the variables g and f, and τ denotes the increment of time. If the coefficients of (2) satisfy(4)biAij+Biaji-biBj=0bi=Bii,j=1,…,Lthen the s-stage partitioned Runge-Kutta method is considered symplectic [41, 42]. For example, the Butecher table of the 2nd-order SPRK method is as follows:(5)In an isotropic lossy material, the Maxwell equations in the time domain can be written as(6)∂E∂t=1ε∇×H-J∂H∂t=-1μ∇×Ewhere H and E denote the magnetic and electric field vectors, respectively, the current density J=σE, and σ, μ, and ε are the conductivity, permeability, and dielectric permittivity, respectively. Letting E=-U and H=∇×A, the generalized Hamiltonian function in lossy media can be expressed as(7)HA,U=∫12μU2+12ε∇×A2-1εJAdVApplying (1), the Maxwell equations can be written in terms of the following canonical equations of the Hamiltonian function:(8)dAdt=dHdU=1μUdUdt=-dHdA=1ε∇2A-σεUwhere ∇2 denotes the Laplacian operator.For the two-dimensional Transverse Magnetic (TM) waves, (8) can be simplified to(9)⁡dAzdt=1μUzdUzdt=1ε∇2Az-σεUzwhere the subscripts of Uz and Az represent the corresponding z-axis parts of the field components.Then,Hx, Hy, and Ez can be derived from Az and Uz as(10)Ez=-UzHx=∂Az∂yHy=-∂Az∂xIn the two-dimensional TM case, the discretization of the SPRK methods is different from the FDTD method [43]. The discrete values of Uz and Az are defined at the same time steps and spatial grid point, as shown in Figure 1, where the time step increment and the spatial interval must be chosen based on the stability requirement to ensure the stability of the SPRK method computation. We use a second-order central difference to approximate the Laplacian operator in this paper.Figure 1
The spatial arrangement of theU and A field components in the 2D TM case.Discretizing (5)–(9), we obtain the following second-order SPRK method iteration scheme:(11)Ai,j1=Ai,jn+Δt2μUi,jn,Ui,jn+1=2ε-Δtσ2ε+ΔtσUi,jn+2Δt2ε+Δtσ∇2Ai,j1,Ai,jn+1=Ai,j1+Δt2μUi,jn+1.where Ui,jn and Ai,jn denote the nodal values of U and A at the nth time step and space point (i,j), and △t represents the time increment.The transmitting boundary condition is adopted to avoid boundary reflections in this paper [44]. The N-order multitransmitting formula (MTF) has the following form:(12)U0n+1=∑j=1N-1j+1CjNUjn+1-j,whereCjn=n!/(n-j)!/j!, Ujn represents the field value at the n-th time step and the spatial computational point -jca△t, and ca expresses the artificial wave velocity. We use a first-order MTF in this paper, and (12) is rewritten as(13)U0n+1=U1n,Taking the correct artificial boundary points as an example, Figure2 shows the typical arrangement of the spatial grid points and computing points. In MTF, the spatial mesh points -j△x generally do not coincide with the spatial computational nodes -jca△, and △x denotes the space interval. A quadratic interpolation scheme is applied to express fields at the computational nodes in terms of those at the spatial mesh points in the domain discretization computation (13).Figure 2
Arrangement of the spatial mesh points and the computational points on the right artificial boundary.The spatial mesh points field values are given by(14)Ui,jn=UnΔt,-jΔx,where the subscripts denote the spatial mesh point.Applying the quadratic interpolation scheme,U1n can be calculated by the relationship:(15)U1n=∑k=13N1,kUi,k-1n,where N1,1=(2-i)(1-i)/2, N1,2=i(2-i), N1,3=i(i-1)/2, i=caΔt/Δx. The points on other boundaries can be calculated in the same method.


### 2.2. Curved Conformal Grid Technology

A conformal grid method based on effective medium parameters is proposed to simulate an arbitrary curved boundary in this paper. A two-dimensional circle is shown as an example to explain the conformal grid method theory in Figure3. Figure 3(a) shows the actual subdivision grid diagram for a circle, the white grid is normal grid points, and the blue grid is actual subdivision grid of a circle. Figure 3(b) shows the subdivision grid diagram for the conventional ladder approximation of a circle and the orange grid is conventional ladder approximation grid of a circle. Figure 3(c) shows the conformal grid subdivision diagram for the circle, and the purple grid is conformal grid. The grid points that are entirely inside the circle or entirely outside of the circle are classified as normal grid points, while the rest are classified as conformal grid points. The effective medium parameters are adopted, to process the conformal grids in this paper.Figure 3
Schematic view for the subdivisions of a circle. (a) An actual model for a circle. (b) The subdivision diagram using the conventional ladder approximation. (c) The conformal grid for a circle.
(a) (b) (c)Next, a conformal cell is extracted from Figure3(c) to illustrate the effective medium parameters for the conformal grid points in the two-dimensional TM wave example. In the symplectic partitioned Runge-Kutta methods, the field components A and U are defined at the same time steps and at the same spatial grid nodes. The calculation for a conformal grid point of the medium in Cartesian coordinates is shown in Figure 4, where F is the sampling point of the field components A and U. Δx and Δy are the width and height of a grid, and Sxy1 and Sxy2 are the area of medium 1 and medium 2, respectively. We assume that the electromagnetic parameters of medium 1 and medium 2 areε1,σ1,μ1 andε2,σ2,μ2, respectively. The field componentsU andA are located at the center of the grid block. The effective values of the dielectric constant, the conductivity, and the permeability coefficient can be obtained by a weighted average over the areas of the grid occupied by the different media.Figure 4
Schematic view for the equivalence of the medium parameters for a conformal grid point. (a) An actual subdivision grid. (b) A conformal grid.The equivalent dielectric constant, conductivity, and permeability coefficient at sampling point F of the field component, as demonstrated in Figure4(b), are given as (16)εzeffF=Sxy1ε1+Sxy2ε2ΔxΔyσzeffF=Sxy1σ1+Sxy2σ2ΔxΔyμzeffF=Sxy1μ1+Sxy2μ2ΔxΔy,where eff denotes the equivalent value of the parameter.Substituting (16) into (11), we obtain the differential iterative equations for the second-order conformal SPRK methods:(17)Ai,j1=Ai,jn+Δt2μeffUi,jn,Ui,jn+1=2εeff-Δtσeff2εeff+ΔtσeffUi,jn+2Δt2εeff+Δtσeff∇2Ai,j1Ai,jn+1=Ai,j1+Δt2μeffUi,jn+1.


### 2.3. Implementation of the Conformal Symplectic Partitioned Runge-Kutta Method with CUDA

Compute unified device architecture (CUDA) programming is well-suited for solving parallel computational problems, where data elements can be divided among multiple threads [45]. The parallel and sequential parts of the CUDA program are executed on the device and the host, respectively. The kernel, as the important component of a CUDA program, executes on a GPU device. CUDA exposes a two-level thread hierarchy to offer developers an optimization method. All threads are divided into blocks, which are grouped into grids [46]. A kernel is executed by being a grid.In this paper, the computation of magnetic and electric field components is considered in a bidimensionalx–y square domain, and the number of threads is decided by the computational domain. As shown in Figure 5, we map each thread to a cell, and each thread block to a set of contiguous cells. In this way, the CUDA grid of thread blocks wraps the complete spatial domain.Figure 5
Mapping between the real domain and the CUDA grid of thread blocks.In the 2D TM case, the CUDA kernel performs the field calculation as depicted in Figure6. There are three CUDA kernels performing the field calculations. The first kernel computesA1, the second updates Un+1, and the third updates An+1. And T is time step; TMax is total time step.Figure 6
A flowchart showing the CUDA implementation of the conformal symplectic partitioned Runge-Kutta methods for determining the field components.In our test, a simple 2D model was used to test the absorption ability of the transmitting boundary condition. We carried out our simulation experiments on a desktop computer with an Intel Core i7-6700K CPU running Windows 7 Pro 64-bit and equipped with an NVIDIA GeForce GTX 1070. For our development environment, we used the CUDA toolkit 7.5, built with Microsoft Visual Studio 2010. The modeled region was a 200 cm × 200 cm square filled with air. A Ricker pulse of central frequency of 1.0 GHz and unity amplitude was used in the center of modeled region. The space increment and time step were chosen as 0.5 cm and 0.01 ns, respectively. The snapshots of the different time steps are shown in Figure7. The diagrams show that the amplitude of U gradually weaken and no reflection appears at the boundary. After the 5 ns time step, the repercussions in modeled region disappear.Figure 7
Distribution of field componentU at different time steps. (a) Snapshots of the U at 2 ns time step. (b) Snapshots of the U at 3 ns time step. (b) Snapshots of the U at 4 ns time step. (d) Snapshots of the U at 5 ns time step.
(a) (b) (c) (d)


## 2.1. Hamiltonian Systems and the Symplectic Partitioned Runge-Kutta Methods [37–39]

A special Hamilton system can be expressed as(1)dpdt=-∂H∂f=pf,gdqdt=∂H∂g=qf,g,where H(p,q) represents the Hamiltonian function, f and g are independent variables, and p and q are dependent variables, such as E and H. The above equations are often solved by the s-stage-partitioned Runge-Kutta methods, which can be specified by the following Butecher table [40]:(2)Applying (2) in (1) yields the Hamiltonian system relations, which are formulated as in the following equations:(3)Pi=pfn+τ∑j=1saijpQjQi=qgn+τ∑j=1sAijqPjfn+1=fn+τ∑i=1sbifQign+1=gn+τ∑i=1sBigPi,i=1,2,…,s,where Qi and Pi are the internal stages corresponding to the variables g and f, and τ denotes the increment of time. If the coefficients of (2) satisfy(4)biAij+Biaji-biBj=0bi=Bii,j=1,…,Lthen the s-stage partitioned Runge-Kutta method is considered symplectic [41, 42]. For example, the Butecher table of the 2nd-order SPRK method is as follows:(5)In an isotropic lossy material, the Maxwell equations in the time domain can be written as(6)∂E∂t=1ε∇×H-J∂H∂t=-1μ∇×Ewhere H and E denote the magnetic and electric field vectors, respectively, the current density J=σE, and σ, μ, and ε are the conductivity, permeability, and dielectric permittivity, respectively. Letting E=-U and H=∇×A, the generalized Hamiltonian function in lossy media can be expressed as(7)HA,U=∫12μU2+12ε∇×A2-1εJAdVApplying (1), the Maxwell equations can be written in terms of the following canonical equations of the Hamiltonian function:(8)dAdt=dHdU=1μUdUdt=-dHdA=1ε∇2A-σεUwhere ∇2 denotes the Laplacian operator.For the two-dimensional Transverse Magnetic (TM) waves, (8) can be simplified to(9)⁡dAzdt=1μUzdUzdt=1ε∇2Az-σεUzwhere the subscripts of Uz and Az represent the corresponding z-axis parts of the field components.Then,Hx, Hy, and Ez can be derived from Az and Uz as(10)Ez=-UzHx=∂Az∂yHy=-∂Az∂xIn the two-dimensional TM case, the discretization of the SPRK methods is different from the FDTD method [43]. The discrete values of Uz and Az are defined at the same time steps and spatial grid point, as shown in Figure 1, where the time step increment and the spatial interval must be chosen based on the stability requirement to ensure the stability of the SPRK method computation. We use a second-order central difference to approximate the Laplacian operator in this paper.Figure 1
The spatial arrangement of theU and A field components in the 2D TM case.Discretizing (5)–(9), we obtain the following second-order SPRK method iteration scheme:(11)Ai,j1=Ai,jn+Δt2μUi,jn,Ui,jn+1=2ε-Δtσ2ε+ΔtσUi,jn+2Δt2ε+Δtσ∇2Ai,j1,Ai,jn+1=Ai,j1+Δt2μUi,jn+1.where Ui,jn and Ai,jn denote the nodal values of U and A at the nth time step and space point (i,j), and △t represents the time increment.The transmitting boundary condition is adopted to avoid boundary reflections in this paper [44]. The N-order multitransmitting formula (MTF) has the following form:(12)U0n+1=∑j=1N-1j+1CjNUjn+1-j,whereCjn=n!/(n-j)!/j!, Ujn represents the field value at the n-th time step and the spatial computational point -jca△t, and ca expresses the artificial wave velocity. We use a first-order MTF in this paper, and (12) is rewritten as(13)U0n+1=U1n,Taking the correct artificial boundary points as an example, Figure2 shows the typical arrangement of the spatial grid points and computing points. In MTF, the spatial mesh points -j△x generally do not coincide with the spatial computational nodes -jca△, and △x denotes the space interval. A quadratic interpolation scheme is applied to express fields at the computational nodes in terms of those at the spatial mesh points in the domain discretization computation (13).Figure 2
Arrangement of the spatial mesh points and the computational points on the right artificial boundary.The spatial mesh points field values are given by(14)Ui,jn=UnΔt,-jΔx,where the subscripts denote the spatial mesh point.Applying the quadratic interpolation scheme,U1n can be calculated by the relationship:(15)U1n=∑k=13N1,kUi,k-1n,where N1,1=(2-i)(1-i)/2, N1,2=i(2-i), N1,3=i(i-1)/2, i=caΔt/Δx. The points on other boundaries can be calculated in the same method.


## 2.2. Curved Conformal Grid Technology

A conformal grid method based on effective medium parameters is proposed to simulate an arbitrary curved boundary in this paper. A two-dimensional circle is shown as an example to explain the conformal grid method theory in Figure3. Figure 3(a) shows the actual subdivision grid diagram for a circle, the white grid is normal grid points, and the blue grid is actual subdivision grid of a circle. Figure 3(b) shows the subdivision grid diagram for the conventional ladder approximation of a circle and the orange grid is conventional ladder approximation grid of a circle. Figure 3(c) shows the conformal grid subdivision diagram for the circle, and the purple grid is conformal grid. The grid points that are entirely inside the circle or entirely outside of the circle are classified as normal grid points, while the rest are classified as conformal grid points. The effective medium parameters are adopted, to process the conformal grids in this paper.Figure 3
Schematic view for the subdivisions of a circle. (a) An actual model for a circle. (b) The subdivision diagram using the conventional ladder approximation. (c) The conformal grid for a circle.
(a) (b) (c)Next, a conformal cell is extracted from Figure3(c) to illustrate the effective medium parameters for the conformal grid points in the two-dimensional TM wave example. In the symplectic partitioned Runge-Kutta methods, the field components A and U are defined at the same time steps and at the same spatial grid nodes. The calculation for a conformal grid point of the medium in Cartesian coordinates is shown in Figure 4, where F is the sampling point of the field components A and U. Δx and Δy are the width and height of a grid, and Sxy1 and Sxy2 are the area of medium 1 and medium 2, respectively. We assume that the electromagnetic parameters of medium 1 and medium 2 areε1,σ1,μ1 andε2,σ2,μ2, respectively. The field componentsU andA are located at the center of the grid block. The effective values of the dielectric constant, the conductivity, and the permeability coefficient can be obtained by a weighted average over the areas of the grid occupied by the different media.Figure 4
Schematic view for the equivalence of the medium parameters for a conformal grid point. (a) An actual subdivision grid. (b) A conformal grid.The equivalent dielectric constant, conductivity, and permeability coefficient at sampling point F of the field component, as demonstrated in Figure4(b), are given as (16)εzeffF=Sxy1ε1+Sxy2ε2ΔxΔyσzeffF=Sxy1σ1+Sxy2σ2ΔxΔyμzeffF=Sxy1μ1+Sxy2μ2ΔxΔy,where eff denotes the equivalent value of the parameter.Substituting (16) into (11), we obtain the differential iterative equations for the second-order conformal SPRK methods:(17)Ai,j1=Ai,jn+Δt2μeffUi,jn,Ui,jn+1=2εeff-Δtσeff2εeff+ΔtσeffUi,jn+2Δt2εeff+Δtσeff∇2Ai,j1Ai,jn+1=Ai,j1+Δt2μeffUi,jn+1.


## 2.3. Implementation of the Conformal Symplectic Partitioned Runge-Kutta Method with CUDA

Compute unified device architecture (CUDA) programming is well-suited for solving parallel computational problems, where data elements can be divided among multiple threads [45]. The parallel and sequential parts of the CUDA program are executed on the device and the host, respectively. The kernel, as the important component of a CUDA program, executes on a GPU device. CUDA exposes a two-level thread hierarchy to offer developers an optimization method. All threads are divided into blocks, which are grouped into grids [46]. A kernel is executed by being a grid.In this paper, the computation of magnetic and electric field components is considered in a bidimensionalx–y square domain, and the number of threads is decided by the computational domain. As shown in Figure 5, we map each thread to a cell, and each thread block to a set of contiguous cells. In this way, the CUDA grid of thread blocks wraps the complete spatial domain.Figure 5
Mapping between the real domain and the CUDA grid of thread blocks.In the 2D TM case, the CUDA kernel performs the field calculation as depicted in Figure6. There are three CUDA kernels performing the field calculations. The first kernel computesA1, the second updates Un+1, and the third updates An+1. And T is time step; TMax is total time step.Figure 6
A flowchart showing the CUDA implementation of the conformal symplectic partitioned Runge-Kutta methods for determining the field components.In our test, a simple 2D model was used to test the absorption ability of the transmitting boundary condition. We carried out our simulation experiments on a desktop computer with an Intel Core i7-6700K CPU running Windows 7 Pro 64-bit and equipped with an NVIDIA GeForce GTX 1070. For our development environment, we used the CUDA toolkit 7.5, built with Microsoft Visual Studio 2010. The modeled region was a 200 cm × 200 cm square filled with air. A Ricker pulse of central frequency of 1.0 GHz and unity amplitude was used in the center of modeled region. The space increment and time step were chosen as 0.5 cm and 0.01 ns, respectively. The snapshots of the different time steps are shown in Figure7. The diagrams show that the amplitude of U gradually weaken and no reflection appears at the boundary. After the 5 ns time step, the repercussions in modeled region disappear.Figure 7
Distribution of field componentU at different time steps. (a) Snapshots of the U at 2 ns time step. (b) Snapshots of the U at 3 ns time step. (b) Snapshots of the U at 4 ns time step. (d) Snapshots of the U at 5 ns time step.
(a) (b) (c) (d)


## 3. Numerical Simulations

### 3.1. Embankment Model

An embankment model was studied, as shown in Figure8, to verify the accuracy and efficiency of the parallel conformal symplectic Euler algorithm. The simulation area was 3.0 m × 2.0 m. The first layer, representing air, had a relative dielectric constant εr=1, and a conductivityσ = 0 mS/m. The second layer modeled asphalt concrete, with a relative dielectric constant εr=6 and a conductivityσ = 1 mS/m. The lowest layer modeled cement-stabilized macadam with εr=30 andσ = 2 mS/m. We assumed that all media were nonmagnetic that their relative permeabilityμ = 1. The time step and spatial increment weredt= 0.01 ns and Δx=Δy=0.5 cm, respectively. The total number of simulated time steps was 2000. A Ricker wavelet of a central frequency of 100 GHz and unity amplitude was adopted to represent the line source excitation in this model, as shown in Figure 9. We placed receivers (Rx) and transmitters (Tx) along the air–earth interface every 0.5 cm. The distance between the receiver and the transmitter was 10 cm.Figure 8
A sketch of the embankment ground-penetrating radar (GPR) model simulated in this study.Figure 9
A Ricker wavelet for a dominant frequency of 1 GHz and unity amplitude.Figure10 shows the single-track radar data of the model at the nodal point (39 cm, 75 cm), obtained with the CUDA-implemented conformal SPRK methods, the CUDA-implemented nonconformal SPRK methods, and the traditional conformal SPRK methods. The traditional conformal SPRK methods required 83.05 s to run, the CUDA-implemented conformal SPRK methods required 6.02 s, and the CUDA-implemented nonconformal SPRK methods required 5.65 s. Results from the CUDA-implemented conformal SPRK methods greatly agree with the traditional conformal SPRK methods, but the proposed CUDA-implemented method reduced the computational effort by more than 92.75%. In comparison to the nonconformal SPRK algorithm, the conformal SPRK algorithm greatly reduced the pseudo-waves generated by the ladder approximation. Figures 11–13 show the wiggle diagrams of the reflection waveforms obtained with the various methods, and red rectangle shows the reflected waves. As shown in Figures 11 and 12, it is obvious that the CUDA-implemented conformal SPRK methods and the traditional conformal SPRK methods have the same accuracy, but the proposed CUDA-implemented method requires approximately one-tenth the time of the traditional method. From Figures 12 and 13, it can be seen that, in comparison to the CUDA-implemented nonconformal SPRK algorithm, the CUDA-implemented conformal SPRK algorithm greatly reduces the pseudo-waves generated by the ladder approximation.Figure 10
A map of the single-track radar data for the embankment model calculated using three methods.Figure 11
The wiggle trace profile of the embankment model simulated by the traditional conformal SPRK methods.Figure 12
The wiggle trace profile of the embankment model simulated by the CUDA-implemented conformal SPRK methods.Figure 13
The wiggle trace profile of the embankment model simulated by the CUDA-implemented nonconformal SPRK methods.


### 3.2. Structural Damage Pavement Model

Cavities are a common quality problem in road engineering. In this case, a three-layer pavement model with cavities was established, and we simulated the GPR profile of this model. As shown in Figure14, the model included an air–earth interface, and a complex three-layered pavement. There was a circular pipe with radii of 5 cm in the first layer, and a 20 cm × 40 cm rectangular cavity in the second layer. The specific model size and dielectric parameters are shown in Figure 14. The permeability of all materials was the free space value, and the excitation source was the same as the first model. The space interval increment and the time step were chosen as Δx=Δy=0.5 cm anddt = 0.01 ns, respectively.Figure 14
A sketch of the structural damage pavement GPR model simulated in this work.We placed receivers and transmitters along the air–earth interface every 1 cm. The distance between the receiver and transmitter was 10 cm. Figures15 and 16 are the plots of the GPR trace profiles that are obtained when simulating this model by using the CUDA-implemented conformal SPRK methods and the CUDA-implemented nonconformal SPRK methods. The yellow circular indicates the position of the pipe, purple rectangle expresses the position of the cavity, and red rectangle shows reflected waves. And From Figures 15 and 16, it can be seen that the boundaries of the pavement model are clearly visible, as are the diffracted waves originating from the structural damage. The results indicate that the CUDA-implemented conformal symplectic partitioned Runge-Kutta algorithm can greatly reduce the pseudo-waves that are generated by the ladder approximation, in comparison to the nonconformal symplectic partitioned Runge-Kutta algorithm.Figure 15
The wiggle trace profile of the three-layer pavement model simulated by the CUDA-implemented conformal SPRK methods.Figure 16
The wiggle trace profile of the three-layer pavement model simulated by the CUDA-implemented nonconformal SPRK methods.


### 3.3. Underground Pipes Model

In this case, we studied the underground pipes model, as shown in Figure17, to verify the accuracy and efficiency between conformal and nonconformal symplectic Euler algorithm and the effect of different dielectric permittivity on reflected wave. This model included an air–earth interface and a sand layer. There were two circular pipes with radii of 5 cm in the sand layer, the A pipe had a relative dielectric constant εr=30 and a conductivityσ = 0 mS/m and B pipe with the relative dielectric constant εr=30 and conductivityσ = 0 mS/m. The excitation source was the same as the first model. The space interval increment and the time step were chosen as Δx=Δy=0.5 cm anddt = 0.01 ns, respectively. The receivers and transmitters are placed along the air–earth interface every 1 cm. The distance between the receiver and transmitter was 10 cm.Figure 17
A sketch of the underground pipes model GPR model simulated in this work.Figures18 and 19 show the wiggle diagrams of the reflection waveforms obtained with the CUDA-implemented conformal and nonconformal SPRK methods. The yellow circular indicates the position of the pipe, orange rectangle expresses A pipe reflected waves, and red rectangle shows B pipe reflected waves. The CUDA-implemented conformal SPRK method has less pseudo-waves than nonconformal method, but the CUDA-implemented conformal method requires more computational time. The CUDA-implemented conformal SPRK method requires 979.056 s and nonconformal method with 761.28 s. And the larger the dielectric permittivity, the more obvious the reflected wave.Figure 18
The wiggle trace profile of the underground pipes model simulated by the CUDA-implemented conformal SPRK methods.Figure 19
The wiggle trace profile of the underground pipes model simulated by the CUDA-implemented nonconformal SPRK methods.


## 3.1. Embankment Model

An embankment model was studied, as shown in Figure8, to verify the accuracy and efficiency of the parallel conformal symplectic Euler algorithm. The simulation area was 3.0 m × 2.0 m. The first layer, representing air, had a relative dielectric constant εr=1, and a conductivityσ = 0 mS/m. The second layer modeled asphalt concrete, with a relative dielectric constant εr=6 and a conductivityσ = 1 mS/m. The lowest layer modeled cement-stabilized macadam with εr=30 andσ = 2 mS/m. We assumed that all media were nonmagnetic that their relative permeabilityμ = 1. The time step and spatial increment weredt= 0.01 ns and Δx=Δy=0.5 cm, respectively. The total number of simulated time steps was 2000. A Ricker wavelet of a central frequency of 100 GHz and unity amplitude was adopted to represent the line source excitation in this model, as shown in Figure 9. We placed receivers (Rx) and transmitters (Tx) along the air–earth interface every 0.5 cm. The distance between the receiver and the transmitter was 10 cm.Figure 8
A sketch of the embankment ground-penetrating radar (GPR) model simulated in this study.Figure 9
A Ricker wavelet for a dominant frequency of 1 GHz and unity amplitude.Figure10 shows the single-track radar data of the model at the nodal point (39 cm, 75 cm), obtained with the CUDA-implemented conformal SPRK methods, the CUDA-implemented nonconformal SPRK methods, and the traditional conformal SPRK methods. The traditional conformal SPRK methods required 83.05 s to run, the CUDA-implemented conformal SPRK methods required 6.02 s, and the CUDA-implemented nonconformal SPRK methods required 5.65 s. Results from the CUDA-implemented conformal SPRK methods greatly agree with the traditional conformal SPRK methods, but the proposed CUDA-implemented method reduced the computational effort by more than 92.75%. In comparison to the nonconformal SPRK algorithm, the conformal SPRK algorithm greatly reduced the pseudo-waves generated by the ladder approximation. Figures 11–13 show the wiggle diagrams of the reflection waveforms obtained with the various methods, and red rectangle shows the reflected waves. As shown in Figures 11 and 12, it is obvious that the CUDA-implemented conformal SPRK methods and the traditional conformal SPRK methods have the same accuracy, but the proposed CUDA-implemented method requires approximately one-tenth the time of the traditional method. From Figures 12 and 13, it can be seen that, in comparison to the CUDA-implemented nonconformal SPRK algorithm, the CUDA-implemented conformal SPRK algorithm greatly reduces the pseudo-waves generated by the ladder approximation.Figure 10
A map of the single-track radar data for the embankment model calculated using three methods.Figure 11
The wiggle trace profile of the embankment model simulated by the traditional conformal SPRK methods.Figure 12
The wiggle trace profile of the embankment model simulated by the CUDA-implemented conformal SPRK methods.Figure 13
The wiggle trace profile of the embankment model simulated by the CUDA-implemented nonconformal SPRK methods.


## 3.2. Structural Damage Pavement Model

Cavities are a common quality problem in road engineering. In this case, a three-layer pavement model with cavities was established, and we simulated the GPR profile of this model. As shown in Figure14, the model included an air–earth interface, and a complex three-layered pavement. There was a circular pipe with radii of 5 cm in the first layer, and a 20 cm × 40 cm rectangular cavity in the second layer. The specific model size and dielectric parameters are shown in Figure 14. The permeability of all materials was the free space value, and the excitation source was the same as the first model. The space interval increment and the time step were chosen as Δx=Δy=0.5 cm anddt = 0.01 ns, respectively.Figure 14
A sketch of the structural damage pavement GPR model simulated in this work.We placed receivers and transmitters along the air–earth interface every 1 cm. The distance between the receiver and transmitter was 10 cm. Figures15 and 16 are the plots of the GPR trace profiles that are obtained when simulating this model by using the CUDA-implemented conformal SPRK methods and the CUDA-implemented nonconformal SPRK methods. The yellow circular indicates the position of the pipe, purple rectangle expresses the position of the cavity, and red rectangle shows reflected waves. And From Figures 15 and 16, it can be seen that the boundaries of the pavement model are clearly visible, as are the diffracted waves originating from the structural damage. The results indicate that the CUDA-implemented conformal symplectic partitioned Runge-Kutta algorithm can greatly reduce the pseudo-waves that are generated by the ladder approximation, in comparison to the nonconformal symplectic partitioned Runge-Kutta algorithm.Figure 15
The wiggle trace profile of the three-layer pavement model simulated by the CUDA-implemented conformal SPRK methods.Figure 16
The wiggle trace profile of the three-layer pavement model simulated by the CUDA-implemented nonconformal SPRK methods.


## 3.3. Underground Pipes Model

In this case, we studied the underground pipes model, as shown in Figure17, to verify the accuracy and efficiency between conformal and nonconformal symplectic Euler algorithm and the effect of different dielectric permittivity on reflected wave. This model included an air–earth interface and a sand layer. There were two circular pipes with radii of 5 cm in the sand layer, the A pipe had a relative dielectric constant εr=30 and a conductivityσ = 0 mS/m and B pipe with the relative dielectric constant εr=30 and conductivityσ = 0 mS/m. The excitation source was the same as the first model. The space interval increment and the time step were chosen as Δx=Δy=0.5 cm anddt = 0.01 ns, respectively. The receivers and transmitters are placed along the air–earth interface every 1 cm. The distance between the receiver and transmitter was 10 cm.Figure 17
A sketch of the underground pipes model GPR model simulated in this work.Figures18 and 19 show the wiggle diagrams of the reflection waveforms obtained with the CUDA-implemented conformal and nonconformal SPRK methods. The yellow circular indicates the position of the pipe, orange rectangle expresses A pipe reflected waves, and red rectangle shows B pipe reflected waves. The CUDA-implemented conformal SPRK method has less pseudo-waves than nonconformal method, but the CUDA-implemented conformal method requires more computational time. The CUDA-implemented conformal SPRK method requires 979.056 s and nonconformal method with 761.28 s. And the larger the dielectric permittivity, the more obvious the reflected wave.Figure 18
The wiggle trace profile of the underground pipes model simulated by the CUDA-implemented conformal SPRK methods.Figure 19
The wiggle trace profile of the underground pipes model simulated by the CUDA-implemented nonconformal SPRK methods.


## 4. Conclusions

This study proposes a precise and efficient forward modeling algorithm of GPR to simulating and analyzing the electromagnetic response of an underground target while using (1) the symplectic partitioned Runge-Kutta algorithm and (2) surface conformal technique, and (3) GPU acceleration technique. The proposed algorithm realizes refined and high-efficient calculation of the electromagnetic response of an arbitrarily shaped underground target and highlights the understanding of the electromagnetic wave propagation laws in the complex underground structure. The reflection waveforms and computational time of the complex underground structure models simulated by CUDA-implemented conformal SPRK methods are compared with those simulated by traditional nonconformal SPRK methods. The CUDA-implemented conformal SPRK algorithm greatly reduces the pseudo-waves and requires approximately one-tenth the time of the traditional method. The performance of the proposed algorithm was validated while using embankment model and structural damage pavement model and underground pipes model numerical simulations. The use of such algorithm is particularly useful because of its accuracy and efficiency and rarity in previously published literature.The authors are now extending the proposed algorithm to three-dimensional GPR forward modeling to improve the processing and interpretation accuracy of the GPR detection data. The proposed algorithm provides an efficient and accurate forward model for the next step of GPR three-dimensional inversion imaging determining accurately and efficiently the distribution of substances in the complex underground structure.



---

*Source: 4025878-2019-06-23.xml*