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

Authors: Jianwei Lei; Zibin Wang; Hongyuan Fang; Xin Ding; Xiaowang Zhang; Man Yang; Haitao Wang
Journal: International Journal of Antennas and Propagation
 (2019)

---

## Abstract

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.


---

## Body

## 1. 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–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 GPU-accelerated FDTD solver integrating it into open source 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–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 two-layer 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.


## 2. FDTD Method

### 2.1. The Governing Equations

Among Maxwell’s electromagnetic equations, the FDTD method uses two curl equations:(1)∇×H=ε∂E∂t+σE∇×E=-μ∂H∂tHereE is the electric field, H 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(2)∂Ez∂y=-μ∂Hx∂t∂Ez∂x=μ∂Hy∂t∂Hy∂x-∂Hx∂y=ε∂EZ∂t+σEzFor 2D TM case, the distribution of the electric and magnetic fields of Yee’s cell of the 2D TM wave is shown in Figure1 [37].Figure 1
Schematic diagram of position of the electric and magnetic field vector of the Yee TM wave.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(3)Hxn+1/2i,j+12=Hxn-1/2i,j+12-Δtμ·Ezni,j+1-Ezni,jΔyHyn+1/2i,j+12=Hyn-1/2i,j+12+Δtμ·Ezni,j+1-Ezni,jΔxEzn+1i,j=CA·Ezni,j+CB·Hyn+1/2i+1/2,j-Hyn+1/2i-1/2,jΔx-Hxn+1/2i,j+1/2-Hxn+1/2i,j-1/2Δywhere CA=1-σΔt/2ε/1+σΔt/2ε, CB=Δt/ε/1+σΔt/2ε.


### 2.2. 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 Figure2, and the second-order Mur absorption condition on the left edge can be expressed as(4)∂2Ez∂x-1c·∂Ez∂t-cμ2·∂Hx∂y=0Figure 2
The arrangement of field components at the left truncating boundary.Equation (3) is discrete at point 1/2,j,t. Using a partial differential in the central difference substitution formula, the resulting difference scheme can be expressed as(5)En+10,j=-En-11,j+cΔt-ΔxcΔt+ΔxEn-11,j+En-10,j+2ΔxcΔt+ΔxEn0,j+En1,j+cΔt2Δx2Δy2cΔt+ΔxEn0,j+1-2En0,j+En0,j-1+En1,j+1-2En1,j+En1,j-1The points on other boundaries and corners can be calculated using the same method.


### 2.3. The 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:(6)c·Δt≤1Δx2+Δy2wherecis the speed of light,Δtis the time step, andΔxandΔyare the spatial step in thex andy direction, respectively.


### 2.4. 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 Figure3. 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.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 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 Ez; A and B are the sampling points of the field components Hy; C and D are the electric field sample points of the field components Hx.Δx andΔy are the width and height of a grid, and Lx and Ly are the length of the different media on a corresponding edge; Sxy1 and Sxy2 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.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.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 Figure4(b) can be expressed as(7)μyeffA=Lx1μ1+Lx2μ2ΔyμyeffB=μ2μxeffC=Lx1μ1+Lx2μ2ΔxμxeffD=μ2The effective magnetic permeability of the conformal mesh, respectively, can be written as(8)μeffA=14μyeffA+μyeffB+μxeff+μxeffDThe effective dielectric constant and the effective conductivity of the N point in Figure4(b) can be shown as(9)εzeffN=Sxy1ε1+Sxy2ε2ΔxΔyσzeffN=Sxy1σ1+Sxy2σ2ΔxΔywhereeffdenotes the equivalent value of the parameter.Substituting (7)-(9) into (3), we obtain the differential iterative equations for the conformal FDTD method:(10)Hxn+1/2i,j+1/2=Hxn-1/2i,j+1/2-Δtμeff·Ezni,j+1-Ezni,jΔyHyn+1/2i,j+1/2=Hyn-1/2i,j+1/2+Δtμeff·Ezni,j+1-Ezni,jΔxEzn+1i,j=CA·Ezni,j+CB·Hyn+1/2i+1/2,j-Hyn+1/2i-1/2,jΔx-Hxn+1/2i,j+1/2-Hxn+1/2i,j-1/2Δywhere CA=1-σeffΔt/2εeff/1+σeffΔt/2εeff, CB=Δt/εeff/1+σeffΔt/2εeff.


## 2.1. The Governing Equations

Among Maxwell’s electromagnetic equations, the FDTD method uses two curl equations:(1)∇×H=ε∂E∂t+σE∇×E=-μ∂H∂tHereE is the electric field, H 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(2)∂Ez∂y=-μ∂Hx∂t∂Ez∂x=μ∂Hy∂t∂Hy∂x-∂Hx∂y=ε∂EZ∂t+σEzFor 2D TM case, the distribution of the electric and magnetic fields of Yee’s cell of the 2D TM wave is shown in Figure1 [37].Figure 1
Schematic diagram of position of the electric and magnetic field vector of the Yee TM wave.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(3)Hxn+1/2i,j+12=Hxn-1/2i,j+12-Δtμ·Ezni,j+1-Ezni,jΔyHyn+1/2i,j+12=Hyn-1/2i,j+12+Δtμ·Ezni,j+1-Ezni,jΔxEzn+1i,j=CA·Ezni,j+CB·Hyn+1/2i+1/2,j-Hyn+1/2i-1/2,jΔx-Hxn+1/2i,j+1/2-Hxn+1/2i,j-1/2Δywhere CA=1-σΔt/2ε/1+σΔt/2ε, CB=Δt/ε/1+σΔt/2ε.


## 2.2. 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 Figure2, and the second-order Mur absorption condition on the left edge can be expressed as(4)∂2Ez∂x-1c·∂Ez∂t-cμ2·∂Hx∂y=0Figure 2
The arrangement of field components at the left truncating boundary.Equation (3) is discrete at point 1/2,j,t. Using a partial differential in the central difference substitution formula, the resulting difference scheme can be expressed as(5)En+10,j=-En-11,j+cΔt-ΔxcΔt+ΔxEn-11,j+En-10,j+2ΔxcΔt+ΔxEn0,j+En1,j+cΔt2Δx2Δy2cΔt+ΔxEn0,j+1-2En0,j+En0,j-1+En1,j+1-2En1,j+En1,j-1The points on other boundaries and corners can be calculated using the same method.


## 2.3. The 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:(6)c·Δt≤1Δx2+Δy2wherecis the speed of light,Δtis the time step, andΔxandΔyare the spatial step in thex andy direction, respectively.


## 2.4. 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 Figure3. 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.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 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 Ez; A and B are the sampling points of the field components Hy; C and D are the electric field sample points of the field components Hx.Δx andΔy are the width and height of a grid, and Lx and Ly are the length of the different media on a corresponding edge; Sxy1 and Sxy2 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.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.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 Figure4(b) can be expressed as(7)μyeffA=Lx1μ1+Lx2μ2ΔyμyeffB=μ2μxeffC=Lx1μ1+Lx2μ2ΔxμxeffD=μ2The effective magnetic permeability of the conformal mesh, respectively, can be written as(8)μeffA=14μyeffA+μyeffB+μxeff+μxeffDThe effective dielectric constant and the effective conductivity of the N point in Figure4(b) can be shown as(9)εzeffN=Sxy1ε1+Sxy2ε2ΔxΔyσzeffN=Sxy1σ1+Sxy2σ2ΔxΔywhereeffdenotes the equivalent value of the parameter.Substituting (7)-(9) into (3), we obtain the differential iterative equations for the conformal FDTD method:(10)Hxn+1/2i,j+1/2=Hxn-1/2i,j+1/2-Δtμeff·Ezni,j+1-Ezni,jΔyHyn+1/2i,j+1/2=Hyn-1/2i,j+1/2+Δtμeff·Ezni,j+1-Ezni,jΔxEzn+1i,j=CA·Ezni,j+CB·Hyn+1/2i+1/2,j-Hyn+1/2i-1/2,jΔx-Hxn+1/2i,j+1/2-Hxn+1/2i,j-1/2Δywhere CA=1-σeffΔt/2εeff/1+σeffΔt/2εeff, CB=Δt/εeff/1+σeffΔt/2εeff.


## 3. The Simulation of GPR Waves on GPU

### 3.1. 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 andH 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.Figure 5
Mapping between real domain and CUDA grid of thread blocks.In 2D case, the flowchart of CUDA kernel performing the field calculation is depicted in Figure6 [39]. There are two CUDA kernels performing the field calculations. The first kernel computesE, and the second updatesH. T is time step; TMax is total time step.Figure 6
The flowchart of the CUDA-implemented FDTD method.


### 3.2. The 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 reflectionAccording 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 Figure7.Figure 7
The forward modeling procedure of GPR electromagnetic wave in the pavement structures.


### 3.3. The 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 sourceE(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 anddt = 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 CUDA-implemented FDTD method, we calculated the 600-time-stepEz 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.Figure 8
Field value of the observation point with two methods.As shown in Figure9, it was the 600-time-step Ez 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.Figure 9
Variation diagram of theEz field value. (a) Traditional FDTD method. (b) CUDA-implemented FDTD method.
(a) (b)


## 3.1. 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 andH 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.Figure 5
Mapping between real domain and CUDA grid of thread blocks.In 2D case, the flowchart of CUDA kernel performing the field calculation is depicted in Figure6 [39]. There are two CUDA kernels performing the field calculations. The first kernel computesE, and the second updatesH. T is time step; TMax is total time step.Figure 6
The flowchart of the CUDA-implemented FDTD method.


## 3.2. The 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 reflectionAccording 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 Figure7.Figure 7
The forward modeling procedure of GPR electromagnetic wave in the pavement structures.


## 3.3. The 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 sourceE(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 anddt = 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 CUDA-implemented FDTD method, we calculated the 600-time-stepEz 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.Figure 8
Field value of the observation point with two methods.As shown in Figure9, it was the 600-time-step Ez 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.Figure 9
Variation diagram of theEz field value. (a) Traditional FDTD method. (b) CUDA-implemented FDTD method.
(a) (b)


## 4. Numerical Simulation

Firstly, we verified the high efficiency of the CUDA-implemented FDTD method by simulating a two-layer pavement model. As shown in Figure10, the top layer of the model represented air with εr=1,σ =0; the second layer represented asphalt concrete with εr=6,σ =1 mS / m; the third layer was the base with εr=30 andσ =2 mS / m. We assumed that all materials are nonmagnetic. The time step and spatial increment were set equal todt = 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.Figure 10
Computational domain of the three-layered model.Figure 11
Map of Ricker waveform.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 Figure12. 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.Figure 12
Map of the single-track radar data calculated using two methods.As shown in Figure12, it can prove that CUDA-implemented FDTD method has the same accuracy as the traditional FDTD method. The speedup ratios achieved are shown in Table 1. When the number of units is small, the 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.Table 1
Speedup ratio of the two methods.
FDTD Cells CPU Times (s) GPU Times (s) Speedup Ratio 128 × 128 1.063 0.398 2.67 256 × 256 4.290 0.743 5.77 512 × 512 18.981 1.879 10.11 1024 × 1024 72.834 6.773 10.75 1536 × 1536 158.545 15.016 10.56 2048 × 2048 299.919 26.276 11.41 2560 ×2560 475.560 41.643 11.58 3072 × 3072 701.439 59.439 11.81 4096 × 4096 1357.060 105.684 13.02GPR profiles of a structurally damaged pavement were simulated in second example. As shown in Figure13, the first layer of the model was an air layer, and the bottom was three-layer 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 asdt= 0.01 ns andΔx = Δy =dl = 0.5 cm, respectively. GPR tracking profiles of the model simulated by CUDA-implemented FDTD method and traditional FDTD method were plotted in Figures 14 and 15.Figure 13
Schematic diagram of a two-layer pavement model with no-dense region.Figure 14
GPR trace profile with traditional FDTD method.Figure 15
GPR trace profile with CUDA-implemented FDTD method.As shown in Figures14 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 Figure16, 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 asdt = 0.01 ns andΔx = Δy =dl = 0.5 cm, respectively.Figure 16
Sketch of the complex underground structure model.We placed receivers and transmitters along the air-earth interface every 1 cm. The distance between the receiver and transmitter was 10 cm. Figures17 and 18 were GPR simulation section using CUDA-implemented FDTD method 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 pseudo-waves than the nonconformal modeling method. The CUDA-implemented 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.Figure 17
GPR simulation section with CUDA-implemented FDTD method.Figure 18
GPR simulation section with traditional FDTD method.


## 5. 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.



---

*Source: 5043028-2019-08-19.xml*