Analysis of GPR Wave Propagation in Complex Underground Structures Using CUDA-Implemented Conformal FDTD Method

Ground penetrating radar (GPR), as a kind of fast, effective, and nondestructive tool, has been widely applied to nondestructive testing of road quality. The finite-difference time-domain method (FDTD) is the common numerical method studying the GPR wave propagation law in layered structure. However, the numerical accuracy and computational efficiency are not high because of the Courant-Friedrichs-Lewy (CFL) stability condition. In order to improve the accuracy and efficiency of FDTD simulation model, a parallel conformal FDTD algorithm based on graphics processor unit (GPU) acceleration technology and surface conformal technique was developed. The numerical simulation results showed that CUDA-implemented conformal FDTD method could greatly reduce computational time and the pseudo-waves generated by the ladder approximation. And the efficiency and accuracy of the proposed method are higher than the traditional FDTD method in simulating GPR wave propagation in two-dimensional (2D) complex underground structures.


Introduction
GPR is a nondestructive way to detect the internal material distribution of underground structure using electromagnetic waves [1]. Its advantages include being fast, continuous, high resolution and no damage. Because of these advantages, it has been widely applied to geological prospecting, civil engineering, highway reconnaissance, and many other areas [2][3][4][5][6]. Numerical simulation of GPR electromagnetic wave propagation in underground structures is helpful to better explain the real GPR profiles. It can provide a way to explore the link between subsurface structures properties and GPR data.
Various approaches have been developed to simulate GPR wave propagation. These include the ray-tracing method [7], the transmission-line-matrix method [8], the finite element method (FEM) [9], the semianalytic mode matching algorithm [10], the pseudo-spectral method [11], the FDTD method [12], the Alternating Direction Implicit-(ADI-) FDTD method [13], and the symplectic method [14]. Among these, the FDTD method proposed by Yee in 1966 [15] has been the most popular simulation tool due to its great stability and convergence. It uses time-step finite-difference to approximate the differential form of Maxwell's equations for Faraday and Ampere's laws and calculates the magnetic and electric field vector components at regularly arranged, distinct points in space. The current state of the art in GPR-oriented FDTD approaches has developed rapidly; for example, a software tool called GprMax has been successfully employed in the detection mechanism of GPR [16], a GPUaccelerated FDTD solver integrating it into open source 2 International Journal of Antennas and Propagation EM simulation software was developed to modeling GPR [17], and a systematic framework for accurate and realistic numerical modeling of GPR for landmine detection has been developed [18]. However, the FDTD method requires large memory and calculation time in simulating GPR wave propagation in 2D complex underground structures due to the limitations of the CFL stability condition [19].
The accuracy and efficiency of forward model are not only related to the performance of the algorithm itself, but also related to computer acceleration technology and modeling method. 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 computation and time also increase correspondingly. For complex problems, the computational resources are difficult to meet. Another approach is to adopt subgrid technology, 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. A last approach is conformal grid technology [20]. This approach can deal with objects with arbitrary surfaces, and the formula deductions and calculations are simple.
In 2006, NVIDIA developed a general-purpose parallel computing architecture called compute unified device architecture (CUDA) that leverages parallel compute engine in NVIDIA graphic processor unit (GPU) to solve many complex computational problems [21]. Due to its excellent programmability and remarkable computational power for parallel executions, CUDA has gained tremendous popularity in the scientific computing society. GPU computing has developed in various fields, such as acoustics [22], biomedical applications [23], plasmonics [24], dispersive media [25], and computational electromagnetism [26][27][28][29][30][31][32][33][34]. FDTD based on the GPU acceleration technique has been applied in the GPR simulation model [35,36], but the pseudo-waves that are generated by the ladder approximation in FDTD modeling method are not considered. This study combines FDTD method and surface conformal technique and GPU acceleration technique proposing a precise and efficient forward modeling algorithm of GPR which can greatly reduce computation time and the pseudo-waves that are generated by the ladder approximation.
In this paper, the simulation model of GPR wave propagation in 2D underground structure was established by CUDA-implemented conformal FDTD method. Line source was chosen to be the incident wave, and second-order Mur absorbing boundary condition was adopted to avoid reflections from the modeled edges. Three numerical examples were provided to verify the high efficiency and accuracy of the CUDA-implemented conformal FDTD method. One obtained the simulated reflection waveform from the twolayer pavement model, the other got the simulated GPR trace profile of the pavement model with structural damage, and the last one demonstrated the GPR simulation section of underground structures with voids.

FDTD Method
. . e Governing Equations. Among Maxwell's electromagnetic equations, the FDTD method uses two curl equations: Here is the electric field, is the magnetic field, is the dielectric constant, is the conductivity, and is the permeability.
Considering 2D Transverse Magnetic (TM) case, (1) can be written as For 2D TM case, the distribution of the electric and magnetic fields of Yee's cell of the 2D TM wave is shown in Figure 1 [37].
By applying the central finite-difference approximation to the spatial and time derivatives of (2) yields, the 2D TM wave discrete equations can be written as

. . Mur Absorbing Boundary
Conditions. An absorbing boundary condition allows the FDTD to calculate approximately the infinite domain problems with limited computer capacity. Its basic idea is to truncate boundaries in the computing area through special treatment or add the ideal mathematical absorption material to create the incident wave on the interface of the electromagnetic wave reflection phenomenon. The absorption boundary applied in this paper is the second-order Mur absorption boundary [38].
For the TM wave, the distribution map of the field component at the left truncated boundary is shown in Figure 2, and the second-order Mur absorption condition on the left edge can be expressed as Equation (3) is discrete at point (1/2, , ). Using a partial differential in the central difference substitution formula, the resulting difference scheme can be expressed as The points on other boundaries and corners can be calculated using the same method.
. . e Stability Condition. Numerical stability is enforced by the Courant condition, which limits the ratio between the spatial discrete step and time step in the simulation. In 2D grid, this condition is given in the following inequality: where c is the speed of light, Δt is the time step, and Δx and Δy are the spatial step in the x and y direction, respectively.
. . Conformal Grid Technology. We adopted a conformal grid method based on effective medium parameters to simulate an arbitrary curved boundary of subsurface structure. A two-dimensional circle is shown as an example to explain the conformal grid method theory in Figure 3. Figure 3(a) shows the actual subdivision grid diagram for a circle, the white grids are normal grid points and the green 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 blue 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.
Next, a conformal cell is extracted from Figure 3(c) to illustrate the effective medium parameters for the conformal grid points in the two-dimensional TM wave example. In the FDTD methods, the electric field E is located at the midpoint of the conformal cell, and the magnetic field H is located at the edge of the conformal cell. The calculation for a conformal grid point of the medium in Cartesian coordinates is shown in Figure 4, where N is the sampling point of the field components z ; A and B are the sampling points of the field components y ; C and D are the electric field sample points of the field components x . Δx and Δy are the width and height of a grid, and x and y are the length of the different media on a corresponding edge; xy1 and xy2 are the area of media 1 and 2, respectively. We assume that the electromagnetic parameters of media 1 and 2 are 1 , 1 , 1 and 2 , 2 , 2 , respectively.
As the magnetic field node is located at the middle point of the grid edge, the effective values of the magnetic conductivity and permeability at each magnetic field node can be calculated by a weighted average of the length of the different media on a corresponding edge. Similarly, the effective values of the dielectric constant and conductivity  at each electric field node can be calculated by a weighted average of the area occupied by different media on the corresponding surface. The effective magnetic permeability of the A, B, C, and D nodes in Figure 4(b) can be expressed as The effective magnetic permeability of the conformal mesh, respectively, can be written as The effective dielectric constant and the effective conductivity of the N point in Figure 4(b) can be shown as where eff denotes the equivalent value of the parameter. Substituting (7)-(9) into (3), we obtain the differential iterative equations for the conformal FDTD method:

The Simulation of GPR Waves on GPU
. . CUDA-Implemented FDTD Method. CUDA programming is well suited to solve problems that can be expressed as parallel computing, which achieves application acceleration performance by mapping data elements to multiple threads. The kernel is a function that is executed on the device, which is the most important component of the CUDA program. The kernel is always represented as a sequential program on the host and executed on the device. CUDA exposes a two-level thread hierarchy that provides an optimization method for developers. All threads are divided into blocks, all of which are grouped into a single grid. The kernel executes in the grid of thread blocks. E and H field components calculate in a bidimensional xy square domain. The number of threads is decided by the calculation domain, which is the number of Yee units, because we map threads to units. As shown in Figure 5, the strategy adopted is mapping each thread to Yee's cell and each thread block to a group of cells. Each thread deals with one  Yee cell during the parallel computing. In this way, the entire calculation domain is wrapped by the CUDA grid. In 2D case, the flowchart of CUDA kernel performing the field calculation is depicted in Figure 6 [39]. There are two CUDA kernels performing the field calculations. The first kernel computes E, and the second updates H. T is time step; T Max is total time step.
. . e Forward Modeling. The steps of the forward modeling of the propagation of GPR in the layer pavement system are as follows: (1) Input basic parameters of pavement structure and build pavement structure model (2) Obtain and input the source of excitation; make numerical difference a discrete source and make it a point corresponding to the time step (3) Calculate the propagation characteristics of GPR electromagnetic waves in the pavement structure with the finite-difference principle of FDTD, including reflection, loss characteristics, analysis, and calculation of the general field (4) Calculate the reflection field (5) Output the signal of radar reflection According to the above steps, a forward model of the GPR propagation in the pavement structure is established. The block diagram of the program is shown in Figure 7. . . e Accuracy of the Proposed Algorithm. The following numerical simulations were performed on a desktop computer with Intel Core i7-6700K CPU and NVIDIA Geforce GTX 1070 Graphics card. The desktop operating system is Windows 7 Pro 64-bit. The development tool of CUDA program is the CUDA Toolkit 7.5 built with Microsoft Visual Studio 2010. In the following, we verified the accuracy of the absorption boundary and algorithm by a simple example.
In the 2D FDTD region in free space, the source of excitation was a sine point source E(t) = sin(2 ⋅ t ⋅ dt ⋅ m) (m= 1.5e13), and Mur boundary was used at the regional truncation boundary. In order to ensure the stability and convergence of FDTD method, the spatial discrete grid and the time step were chosen as Δx = Δy = 0.25 cm and dt = 0.005 ns, respectively. The point source was initiated in the center of the 50 cm × 50 cm square free space by using the 2D TM wave as an example. A CUDA calculation model was established for each thread to solve a Yee cellular (magnetic) field. There were 16 × 16 threads in each block, and the simulated area was divided into 200 × 200 grids. The point source and observation point were located at (25 cm, 25 cm) and (25 cm, 30 cm), respectively. Using traditional FDTD method and CUDAimplemented FDTD method, we calculated the 600-timestep z field change with the same parameter calculation. As shown in Figure 8, the simulation results of traditional FDTD method and CUDA-implemented FDTD method were with excellent agreement. As shown in Figure 9, it was the 600-time-step z field change diagram calculated by the CUDA-implemented FDTD method and traditional FDTD method. It can be seen that GPR electromagnetic wave reached the boundary without causing reflection due to the setting of absorbing boundary. It also simulated intuitively the GPR electromagnetic wave propagation in free space and proved the accuracy of the Mur absorbing boundary.

Numerical Simulation
Firstly, we verified the high efficiency of the CUDAimplemented FDTD method by simulating a two-layer pavement model. As shown in Figure 10, the top layer of the model represented air with = 1, = 0; the second layer represented asphalt concrete with = 6, = 1 mS / m; the third layer was the base with = 30 and = 2 mS / m. We assumed that all materials are nonmagnetic. The time step and spatial increment were set equal to dt = 0.005 ns and Δx = Δy = dl = 0.25 cm, respectively. We used a Ricker wavelet with a unit amplitude and a center frequency of 1.0 GHz (Figure 11) to represent the line source of the GPR transmitter.
The transmitter was at the point (15 cm, 100 cm) and the receiver was at the point (15 cm, 110 cm), respectively. The simulation results of traditional FDTD method and CUDA-implemented FDTD method after 1000 time steps were shown in Figure 12. Excellent agreement was reached. It was proved that traditional FDTD method and CUDA-implemented FDTD method have the same accuracy. Moreover, the CUDA-implemented FDTD method consumed 26.276 s, while the traditional FDTD method increased to 299.919 s. The former saved 91.24% of the time. The result demonstrated that CUDA-implemented FDTD method greatly reduced computation time.
As shown in Figure 12, it can prove that CUDAimplemented FDTD method has the same accuracy as the traditional FDTD method. The speedup ratios achieved are shown in Table 1 acceleration ratio is also small, so the acceleration performance obtained is not significant. When the number of units is getting larger, better acceleration performance will be obtained. In particular, when the FDTD cells are up to 4096×4096, we can obtain more than 13 times speedup ratio. GPR profiles of a structurally damaged pavement were simulated in second example. As shown in Figure 13, the first layer of the model was an air layer, and the bottom was threelayer pavement with structural damage. There was a 10 mm wide crack in first layer. There were a circular no-dense region with radius of 5 cm and a 30 cm × 40 cm rectangular cavity in second layer. In this example, assuming that the porosity of the region was 10%, all dielectric parameters of the material were shown in Figure 13. Receivers and transmitters were placed along the air-earth interface every 1 cm to measure the reflection waveform of the GPR wave in the pavement structure. The distance between receiver and transmitter was 5 cm. The incident wave in first example was used to represent the line source of the transmitter in this model. The increments of the time step and spatial interval were selected as dt = 0.01 ns and Δx = Δy = dl = 0.5 cm, respectively. GPR tracking profiles of the model simulated by CUDAimplemented FDTD method and traditional FDTD method were plotted in Figures 14 and 15.   As shown in Figures 14 and 15, a good agreement has been reached. However, the traditional FDTD method required 2.7357 h; the CUDA-implemented FDTD method required 0.1771 h. The results showed that the reflection waveform of CUDA-implemented FDTD method was consistent with the traditional FDTD method, but the CUDA-implemented method reduced the computational time by more than 93.5%.
In this example, a complex underground structure with two cavities was established. We verified the high efficiency and accuracy of CUDA-implemented conformal FDTD method by using conformal and nonconformal modeling method for two cavities, respectively. As shown in Figure 16, this model included an air-earth interface and a soil layer with two cavities with a radius of 5 cm. The model size and dielectric parameters were shown in Figure 16. The excitation source was the same as the second model. The time step and spatial interval increment were selected as dt = 0.01 ns and Δx = Δy = dl = 0.5 cm, respectively.
We placed receivers and transmitters along the air-earth interface every 1 cm. The distance between the receiver and transmitter was 10 cm. Figures 17 and 18 were GPR simulation section using CUDA-implemented FDTD method International Journal of Antennas and Propagation and traditional FDTD method. The right cavity used conformal modeling method, the left one with nonconformal modeling method. And from Figures 17 and 18, it can be seen that the conformal modeling method has fewer pseudowaves than the nonconformal modeling method. The CUDAimplemented FDTD method required 967.23 s, and traditional FDTD method required 11606.76 s. The results indicate that CUDA-implemented conformal FDTD algorithm can greatly reduce computational time and the pseudo-waves that are generated by the ladder approximation.

Conclusions
This paper presented a CUDA-implemented conformal FDTD method for simulating GPR wave propagation in 2D low-loss medium. The high efficiency and accuracy of the method were verified by three numerical examples. The results showed that CUDA-implemented conformal FDTD algorithm can greatly reduce computation time and the pseudo-waves that were generated by the ladder approximation. In addition, we obtained the GPR simulation profile of complex underground structure by forward modeling using the CUDA-implemented conformal FDTD method. By analyzing the simulated profiles of GPR, we can deepen our understanding of the propagation of GPR waves in underground space and better interpret the real GPR profiles. This paper focuses on 2D isotropic conditions. However, this method can be extended to three-dimensional (3D) or anisotropic situations. Through further research, we will develop a CUDA-implemented conformal FDTD to simulate the propagation of GPR waves in anisotropic media or 3D structures.

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.