Complexity
Volume 2019 (2019), Article ID 4025878, 14 pages
https://doi.org/10.1155/2019/4025878

Research Article

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

Hongyuan Fang,1,2 Jianwei Lei,1,2 Man Yang,1,2 and Ziwei Li¹

¹College of Water Conservancy & Environmental Engineering, Zhengzhou University, Zhengzhou 450001, China
²National Local Joint Engineering Laboratory of Major Infrastructure Testing and Rehabilitation Technology, Zhengzhou 450001, China

Correspondence should be addressed to Jianwei Lei; 15333827823@163.com

Received 11 January 2019; Accepted 21 April 2019; Published 23 June 2019

Academic Editor: Mahardhika Pratama

Copyright © 2019 Hongyuan Fang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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: Section 2 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

where represents the Hamiltonian function, and are independent variables, and and are dependent variables, such as and . The above equations are often solved by the -stage-partitioned Runge-Kutta methods, which can be specified by the following Butecher table [40]:

Applying (2) in (1) yields the Hamiltonian system relations, which are formulated as in the following equations:

where and are the internal stages corresponding to the variables and , and denotes the increment of time. If the coefficients of (2) satisfy

then 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:

In an isotropic lossy material, the Maxwell equations in the time domain can be written as

where and denote the magnetic and electric field vectors, respectively, the current density , and , , and are the conductivity, permeability, and dielectric permittivity, respectively. Letting and , the generalized Hamiltonian function in lossy media can be expressed as

Applying (1), the Maxwell equations can be written in terms of the following canonical equations of the Hamiltonian function:

where denotes the Laplacian operator.

For the two-dimensional Transverse Magnetic (TM) waves, (8) can be simplified to

where the subscripts of and represent the corresponding z-axis parts of the field components.

Then, , , and can be derived from and as

In the two-dimensional TM case, the discretization of the SPRK methods is different from the FDTD method [43]. The discrete values of and 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 the and field components in the 2D TM case.

Discretizing (5)–(9), we obtain the following second-order SPRK method iteration scheme:

where and denote the nodal values of and at the th time step and space point , and 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:

where , represents the field value at the -th time step and the spatial computational point , and expresses the artificial wave velocity. We use a first-order MTF in this paper, and (12) is rewritten as

Taking the correct artificial boundary points as an example, Figure 2 shows the typical arrangement of the spatial grid points and computing points. In MTF, the spatial mesh points generally do not coincide with the spatial computational nodes , and 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

where the subscripts denote the spatial mesh point.

Applying the quadratic interpolation scheme, can be calculated by the relationship:

where , , , . 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 Figure 3. 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.

(a)

(b)

(c)

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.

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 symplectic partitioned Runge-Kutta methods, the field components and 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 and . and are the width and height of a grid, and and are the area of medium 1 and medium 2, respectively. We assume that the electromagnetic parameters of medium 1 and medium 2 are ε₁, σ₁, μ₁ and ε₂, σ₂, μ₂, respectively. The field components U and A 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 Figure 4(b), are given as

where denotes the equivalent value of the parameter.

Substituting (16) into (11), we obtain the differential iterative equations for the second-order conformal SPRK methods:

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 bidimensional x–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 Figure 6. There are three CUDA kernels performing the field calculations. The first kernel computes A¹, the second updates , and the third updates . And T is time step; 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 Figure 7. The diagrams show that the amplitude of gradually weaken and no reflection appears at the boundary. After the 5 ns time step, the repercussions in modeled region disappear.

(a)

(b)

(c)

(d)

Figure 7: Distribution of field component at different time steps. (a) Snapshots of the at 2 ns time step. (b) Snapshots of the at 3 ns time step. (b) Snapshots of the at 4 ns time step. (d) Snapshots of the at 5 ns time step.

3. Numerical Simulations

3.1. Embankment Model

An embankment model was studied, as shown in Figure 8, 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 , and a conductivity σ = 0 mS/m. The second layer modeled asphalt concrete, with a relative dielectric constant and a conductivity σ = 1 mS/m. The lowest layer modeled cement-stabilized macadam with and σ = 2 mS/m. We assumed that all media were nonmagnetic that their relative permeability μ = 1. The time step and spatial increment were dt = 0.01 ns and 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.

Figure 10 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 Figure 14, 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 cm and dt = 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. Figures 15 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 Figure 17, 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 and a conductivity σ = 0 mS/m and B pipe with the relative dielectric constant 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 cm and dt = 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.

Figures 18 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.

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 is no conflict of interest regarding the publication of this paper.

Acknowledgments

This work was supported by the National Key Research and Development Program of China (No. 2017YFC1501204), the National Natural Science Foundation of China (No. 51678536, 41404096), the Program for Science and Technology Innovation Talents in Universities of Henan Province (Grant No. 19HASTIT043), the Outstanding Young Talent Research Fund of Zhengzhou University (1621323001), and the Program for Innovative Research Team (in Science and Technology) in University of Henan Province (18IRTSTHN007). The authors would like to thank for these financial supports.

References

D. J. Daniels, Ground Penetrating Radar, The Intuition of Electrical Engineering, London, UK, 2nd edition, 2004.

K. R. Maser, T. J. Holland, R. Roberts, and J. Popovics, “NDE methods for quality assurance of new pavement thickness,” International Journal of Pavement Engineering, vol. 7, no. 1, pp. 1–10, 2006.

B. Park, J. Kim, J. Lee, M.-S. Kang, and Y.-K. An, “Underground object classification for urban roads using instantaneous phase analysis of ground-penetrating radar (GPR) data,” Remote Sensing, vol. 10, no. 9, 2018.

O. Shamir, N. Goldshleger, U. Basson, and M. Reshef, “Laboratory measurements of subsurface spatial moisture content by ground-penetrating radar (GPR) diffraction and reflection imaging of agricultural soils,” Remote Sensing, vol. 10, no. 10, article 1667, 2018.

L. Yi, L. Zou, and M. Sato, “Practical approach for high-resolution airport pavement inspection with the yakumo multistatic array ground-penetrating radar system,” Sensors, vol. 18, no. 8, 2018.

A. Benedetto, F. Tosti, L. Bianchini Ciampoli, and F. D’Amico, “An overview of ground-penetrating radar signal processing techniques for road inspections,” Signal Processing, vol. 132, pp. 201–209, 2017.

I. Prokopovich, A. Popov, L. Pajewski, and M. Marciniak, “Application of coupled-wave Wentzel-Kramers-Brillouin approximation to ground penetrating radar,” Remote Sensing, vol. 10, no. 1, 2018.

J. Zhang, S. Ye, L. Yi et al., “A hybrid method applied to improve the efficiency of full-waveform inversion for pavement characterization,” Sensors, vol. 18, no. 9, article 2916, 2018.

L. Nuzzo, G. Leucci, S. Negri, M. T. Carrozzo, and T. Quarta, “Application of 3D visualization techniques in the analysis of GPR data for archaeology,” Annals of Geophysics, vol. 45, no. 2, pp. 321–338, 2002.

R. Persico, Introduction to Ground Penetrating Radar: Inverse Scattering and Data Processing, Wiley-IEEE Press, Hoboken, NJ, USA, 2014.

X. Hu, L. Ye, S. Pang, and J. Shan, “Semi-global filtering of airborne LiDAR data for fast extraction of digital terrain models,” Remote Sensing, vol. 7, no. 8, pp. 10996–11015, 2015.

L. Zhang, Z. Zeng, J. Li et al., “Parameter estimation of lunar regolith from lunar penetrating radar data,” Sensors, vol. 18, no. 9, 2018.

G. Ludeno, L. Capozzoli, E. Rizzo, F. Soldovieri, and I. Catapano, “A microwave tomography strategy for underwater imaging via ground penetrating radar,” Remote Sensing, vol. 10, no. 9, 2018.

F. Soldovieri, O. Lopera, and S. Lambot, “Combination of advanced inversion techniques for an accurate target localization via GPR for demining applications,” IEEE Transactions on Geoscience and Remote Sensing, vol. 49, no. 1, pp. 451–461, 2011.

W. Bi, Y. Zhao, C. An, and S. Hu, “Clutter elimination and random-noise denoising of GPR signals using an SVD method based on the Hankel matrix in the local frequency domain,” Sensors, vol. 18, no. 10, 2018.

X. Xu, C. Brekke, A. P. Doulgeris, and F. Melandso, “Numerical analysis of microwave scattering from layered sea ice based on the finite element method,” Remote Sensing, vol. 10, no. 9, 2018.

A. Villarino, B. Riveiro, D. Gonzalez-Aguilera, and L. J. Sánchez-Aparicio, “The integration of geotechnologies in the evaluation of a wine cellar structure through the finite element method,” Remote Sensing, vol. 6, no. 11, pp. 11107–11126, 2014.

Y. Huang, J. Zhang, and Q. H. Liu, “Three-dimensional GPR ray tracing based on wavefront expansion with irregular cells,” IEEE Transactions on Geoscience and Remote Sensing, vol. 49, no. 2, pp. 679–687, 2011.

M. Van Leeuwen, N. C. Coops, and T. Andrew Black, “Using stochastic ray tracing to simulate a dense time series of gross primary productivity,” Remote Sensing, vol. 7, no. 12, pp. 17272–17290, 2015.

G.-X. Fan, Q. H. Liu, and J. S. Hesthaven, “Multidomain pseudospectral time-domain simulations of scattering by objects buried in lossy media,” IEEE Transactions on Geoscience and Remote Sensing, vol. 40, no. 6, pp. 1366–1373, 2002.

D. Vande Ginste, H. Rogier, F. Olyslager, and D. De Zutter, “A fast multipole method for layered media based on the application of perfectly matched layers-the 2-D case,” IEEE Transactions on Antennas and Propagation, vol. 52, no. 10, pp. 2631–2640, 2004.

J. Irving and R. Knight, “Numerical modeling of ground-penetrating radar in 2-D using MATLAB,” Computers Geosicences, vol. 32, no. 9, pp. 1247–1258, 2006.

T. M. Millington and N. J. Cassidy, “Optimising GPR modelling: a practical, multi-threaded approach to 3D FDTD numerical modelling,” Computers & Geosciences, vol. 36, no. 9, pp. 1135–1144, 2010.

D. Feng, X. Wang, and B. Zhang, “Specific evaluation of tunnel lining multi-defects by all-refined GPR simulation method using hybrid algorithm of FETD and FDTD,” Construction and Building Materials, vol. 185, pp. 220–229, 2018.

D.-S. Feng and Q.-W. Dai, “GPR numerical simulation of full wave field based on UPML boundary condition of ADI-FDTD,” NDT & E International, vol. 44, no. 6, pp. 495–504, 2011.

H. Fang, G. Lin, and R. Zhang, “The first-order symplectic euler method for simulation of GPR wave propagation in pavement structure,” IEEE Transactions on Geoscience and Remote Sensing, vol. 51, no. 1, pp. 93–98, 2013.

C. Xiang, F. Kong, K. Li, and M. Liu, “A high-order symplectic FDTD scheme for the maxwell-schrodinger system,” IEEE Journal of Quantum Electronics, vol. 54, no. 1, 2018.

X.-K. Wei, X. Zhang, N. Diamanti, W. Shao, and C. D. Sarris, “Subgridded FDTD modeling of ground penetrating radar scenarios beyond the courant stability limit,” IEEE Transactions on Geoscience and Remote Sensing, vol. 55, no. 12, pp. 7189–7198, 2017.

K. Hong, Y. Gao, A. Ullah et al., “Multi-scale CFD modeling of gas-solid bubbling fluidization accounting for sub-grid information,” Advanced Powder Technology, vol. 29, no. 3, pp. 488–498, 2018.

X.-J. Hu, D.-B. Ge, and B. Wei, “Study on MCFDTD for 3-D coated targets by using effective parameters,” Systems Engineering and Electronics, vol. 28, no. 11, pp. 1652–1667, 2006.

Y. Huang, Z. Chen, B. Wu et al., “Estimating roof solar energy potential in the downtown area using a GPU-accelerated solar radiation model and airborne LiDAR data,” Remote Sensing, vol. 7, no. 12, pp. 17212–17233, 2015.

L. Linghu, J. Wu, Z. Wu, and X. Wang, “Parallel computation of EM backscattering from large three-dimensional sea surface with CUDA,” Sensors, vol. 18, no. 11, article 3656, 2018.

E. Martel, R. Lazcano, J. López et al., “Implementation of the principal component analysis onto high-performance computer facilities for hyperspectral dimensionality reduction: results and comparisons,” Remote Sensing, vol. 10, no. 6, 2018.

J. J. Lopez, D. Carnicero, N. Ferrando, and J. Escolano, “Parallelization of the finite-difference time-domain method for room acoustics modelling based on CUDA,” Mathematical and Computer Modelling, vol. 57, no. 7, pp. 1822–1831, 2013.

E. Torti, G. Florimbi, F. Castelli et al., “Parallel K-means clustering for brain cancer detection using hyperspectral images,” Electronics, vol. 7, no. 11, article 283, 2018.

C.-C. Kung, “Study on consulting air combat simulation of cluster UAV based on mixed parallel computing framework of graphics processing unit,” Electronics, vol. 7, no. 9, 2018.

T. Hirono, W. W. Lui, and K. Yokoyama, “Time-domain simulation of electromagnetic field using a symplectic integrator,” IEEE Microwave and Guided Wave Letters, vol. 7, no. 9, pp. 279–281, 1997.

Z.-X. Huang and X.-L. Wu, “Symplectic partitioned Runge-Kutta scheme for Maxwell’s equations,” International Journal of Quantum Chemistry, vol. 106, no. 4, pp. 839–842, 2006.

L.-L. Jiang, J.-F. Mao, and X.-L. Wu, “Symplectic finite-difference time-domain method for Maxwell equations,” IEEE Transactions on Magnetics, vol. 42, no. 8, pp. 1991–1995, 2006.

J. C. Butcher, Numerical Methods for Ordinary Differential Equations, John Wiley & Sons, Chichester, UK, 2nd edition, 2008.

T. Hirono, W. W. Lui, K. Yokoyama, and S. Seki, “Stability and numerical dispersion of symplectic fourth-order time-domain schemes for optical field simulation,” Journal of Lightwave Technology, vol. 16, no. 10, pp. 1915–1920, 1998.

I. Saitoh and N. Takahashi, “Stability of symplectic finite-difference time-domain methods,” IEEE Transactions on Magnetics, vol. 38, no. 2, pp. 665–668, 2002.

H. Y. Fang and G. Lin, “Symplectic partitioned Runge-Kutta methods for two-dimensional numerical model of ground penetrating radar,” Computers & Geosciences, vol. 49, pp. 323–329, 2012.

Z. P. Liao, “Extrapolation non-reflecting boundary conditions,” Wave Motion, vol. 24, no. 2, pp. 117–138, 1996.

NVIDIA Corporation, NVIDIA CUDA C Programming Guide, version 9.1 edition, 2017.

H. Jiang, S. Chen, D. Li, C. Wang, and J. Yang, “Papaya tree detection with UAV images using a GPU-accelerated scale-space filtering method,” Remote Sensing, vol. 9, no. 7, article 721, 2017.