Abstract
The finite-difference time-domain (FDTD) method adopts the most popular numerical model simulating ground penetrating radar (GPR) wave propagation in an underground structure. However, a staircase approximation method is usually adopted to simulate the curved boundary of an irregular object in the FDTD and symplectic partitioned Runge-Kutta (SPRK) methods. The approximate processing of rectangular mesh parameters will result in calculation errors and virtual surface waves for irregular targets of an underground structure. In this paper, we examine transverse mode (TM) electromagnetic waves with numerical models of electromagnetic wave propagation in geoelectric structures with conformal finite-difference time-domain (CFDTD) method technology in which the effective dielectric parameters are used to accurately simulate the dielectric surface and to absorb waves at the edges of the grid. The third orders of the transmission boundary are used in this paper. Additionally, three complex geocentric models of inclined layered media, spherical media, and three-layered pavement model with structural damages are set up for simulation calculations, then we carry out the actual radar wave detection in a laboratory as the fourth numerical example. Comparison of simulated reflectance waveform of FDTD, symplectic partitioned Runge-Kutta (SPRK), and CFDTD methods shows that at least 50% of the virtual waves can be reduced by using the proposed algorithm. Wiggle diagrams of FDTD and CFDTD methods show that much of the virtual waves have been reduced, and the radar image is clearer than before. This provides a method for the detection of complex geoelectric and layered structures in actual engineering.
1. Introduction
GPR has been widely used in underground engineering because of its complexity and diversity. Through the analysis of computed geoelectrical echo signals, the propagation of GPR signals can be further understood. The method of finite-difference time-domain (FDTD) analysis is the most widely used numerical simulation technique of GPR [1], which has the characteristics of directness and generality, and electromagnetic parameters of the target are reflected in the electromagnetic field of every grid [2]. In order to produce numerical simulation of GPR, mesh generations of complex geoelectric structures were constructed, the extent of which will directly affect the computing accuracy of the calculated results [3].
The reflection waves of radar have been examined and discussed. In a series study, Li et al. [4] have presented a simplified conformal FDTD (CFDTD) scheme to overcome staircase error, and they have compared computed RCS’s with that from the conventional FDTD and other methods and measured values. Chai et al. [5] have proposed a hybrid technique combining the pseudospectral time-domain (PSTD) method with the alternating-direction implicit conformal finite-difference time-domain (ADI-CFDTD) method to solve 3D mixed-scale problems in computational electromagnetics. Fang and Lin [6] have applied the standard FDTD method and the symplectic partitioned Runge-Kutta (SPRK) method to verify the performance of the SPRK algorithms. Dai et al. [7] have introduced element subdivision, interpolation, integration, general synthesis, and GPR finite element wave equation deducing process in details. Lin and Smithe [8] have provided a conformal finite-difference time-domain (CFDTD) method to simulate a gyrotron cavity. Lin et al. [9] have used a 3D conformal finite-difference time-domain particle-in-cell simulation to study the possibility of dynamic after cavity interaction (ACI) in a downscaled gyrotron model. Chen et al. [10] have reported a method for real-time 3D monitoring of thermal therapy through the use of noncontact microwave imaging. Zhang et al. [11] have found that the finite element method (FEM) simulation and measurement of the designed tags show strong resonances in the GPR spectrum. Chen et al. [12] have derived the analytical wavenumber- (-) space domain propagators underlying the fractional Laplacian wave equation of the pseudospectral time-domain (PSTD) method; however, PSTD is limited by the limited memory and limited operation speed of a single machine. Tobon et al. [13] have proposed a new discontinuous Galerkin spectral element time-domain (DG-SETD) method for Maxwell’s equations to analyze three-dimensional (3D) transient electromagnetic phenomena. The SETD method is suitable for the numerical simulation of elastic waves under a complex boundary model, but its meshing is complex [14, 15].
Rectangular grids are typically used for FDTD grid discretization; however, it is impossible to accurately simulate irregularly shaped underground targets [16–20]. There are two methods to improve the accuracy of electromagnetic response calculations for irregularly shaped, subterranean targets. The first method is to increase mesh density, thereby reducing the error of the staircase approximation [21, 22]. Unfortunately, this method increases computation complexity and memory consumption, especially for complex problems. The second method is to use a subgrid algorithm in FDTD [23, 24]. In this calculation method, different grid sizes are adopted for different parts of the structure; the subdivision of grids with irregular shapes is adopted for the target areas, while normal grid sizes are adopted for other areas. By using a subgrid algorithm, the calculation accuracy and stability of this algorithm can be obtained, but in the coarse and fine meshes of the jump interface, it is easy to produce reflection waves, rendering the calculations unstable [25–29].
Due to the discreteness of the FDTD grid, conformal grid technology is a method that does not increase grid density or change the division of the grids. Additionally, it can be applied to a grid containing two or more mediums, as well as be used for the calculation of arbitrary curved surfaces [30]. The iterative formulae of the FDTD method combined with conformal mesh technology for computing the parameters are easily deduced and programmed. This technique has been used in research on the electromagnetic characteristics of the thin absorbing material of stealth fighter jets and to predict the surface current distribution of aircraft carrier flight decks under high-power electromagnetic pulses [31, 32].
In this paper, we examine transverse mode (TM) electromagnetic waves with numerical models of electromagnetic wave propagation in geoelectric structures with CFDTD technology. The third orders of the transmission boundary are used in this paper. By comparing the inclined layered medium with the circular medium models of the GPR profiles among the FDTD, SPRK, and CFDTD methods, we can obtain a reduction of at least 50% of the virtual waves by using the proposed algorithm. The GPR profile of the three-layered pavement model with structural damages is simulated to verify the application of the CFDTD method. To verify the correctness of the CFDTD method, the simulated data are compared with the measured data in the fourth example.
2. CFDTD Method
Maxwell’s curl equations are as follows: where and denote electric field intensity and magnetic flux, respectively. , , , represent dielectric flux density, magnetic flux density, current density, and magnetic flux density, respectively. For TM wave, equations (1) and (2) for 2D can be simplified to where , , , and represent the electrical permittivity, electrical conductivity, magnetic permeability, and magnetic conductivity, respectively.
First of all, we introduce the selection of medium parameters in regular grids. Figure 1(a) shows the integral loop of magnetic flux , and Figure 1(b) shows the integral loop of electric field intensity . It can be seen that magnetic flux is located at the center of the four cell centers in Figure 1(a). The value of the previous time step of the node and the value of the first half of the time step of the four electric field nodes around the node are used to calculate the component of the magnetic field. These four electric field nodes are located at the center of the four planes containing the edge of the cell in which the magnetic field node is located.
(a)
(b)
Then, calculation of the parameters of the magnetic field requires the average of the medium parameters of the four cell centers. The magnetic permeability and magnetic conductivity can be obtained as
It can be seen from Figure 1(b) that the electric field line connects the upper and lower cell centers, so when calculating the parameters of the electric field, the average value of the dielectric parameters of the upper and lower cell centers are required. The parameters of the electric field can be written as
Moreover, we introduce the concept of a conformal mesh. In Figure 2, the discretization technique of the different models of a circular object can be seen. Figure 2(a) shows the actual discreteness of the circular target, and Figure 2(b) shows the traditional FDTD approach, using the staircase approximation method of handling the surface boundary. Figure 2(c) shows the mesh generation model of the CFDTD method. The yellow color in Figure 2(a) represents the actual circular medium, the saddle brown color in Figure 2(b) represents the material parameters used by the traditional FDTD method, the yellow color in Figure 2(c) represents the mesh processed by the traditional FDTD method, and the red color represents the mesh processed by the CFDTD method. The internal and external surface mediums are called the regular grid, and the others are termed the conformal grid. For the conformal grid, the effective medium parameter method has been adopted to calculate it. Figure 3 shows before the effectivity of the medium parameters in a conformal grid in a 3D model. The yellow color in Figure 3 represents the actual circular medium. Figure 4(a) shows before the effectivity of the medium parameters in a conformal grid, and Figure 4(b) shows after the effectivity of the medium parameters in a conformal grid. The yellow color in Figure 4(a) represents the actual circular medium, and the water red color in Figure 4(b) represents the range of parameter filling after treatment by the CFDTD method. The Yee cells of a conformal grid [22, 30] are shown in Figures 3 and 4 in order to illustrate the effective parameters of the conformal grid. In this approach, we suppose that both the electric and magnetic fields inside the distorted cell are located at the same positions as those in the conventional FDTD scheme (Figures 3 and 4), and that Faraday’s law is applied over the entire FDTD cell, rather than only in the distorted part.
(a)
(b)
(c)
Expressed in a rectangular coordinate system, Figure 4 shows the conformal grid calculated for the medium. , , , and are located at the middle of the edge, where and are the sample points of magnetic field and and are the sample points of magnetic field . Point is the sample point for the electric field . The magnetic and electromagnetic parameters of medium 1 and medium 2 are , , , and and , , , and , respectively. and are the space steps of the orthogonal grid, and and are the cell lengths of medium 1 along the and directions, respectively. and are the cell lengths of medium 2 along the and directions, respectively. represents the area of the cell occupied by medium 1, and represents the area of the cell occupied by medium 2. The electric and magnetic fields are sampled at the same position as the conventional grids. The interface of the medium in the cell is approximately a straight line in this method.
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, because the electric field node is located at the center of the grid, 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, as shown in Figure 4. In Figure 4, the effective magnetic conductivity and permeability of the , , , and magnetic field nodes can be expressed as where , , , and represent the magnetic permeability of the four sides of the cell, respectively. , , , and represent the magnetic conductivity of the four sides of the cell, respectively.
As shown in Figures 3 and 4, the components of the magnetic field and are located on the four edges of the cell. We calculate the equivalent parameters of the magnetic field in the conformal cell by using the average of the equivalent parameters of the four edges. With these relationships, the effective magnetic permeability and the effective magnetic conductivity of the conformal mesh, respectively, can be written as
The effective dielectric constant and conductivity of the electric field sampling point can be obtained by taking the weighted average of the area occupied on the corresponding side of medium 1 and medium 2. The component of the electric field is located at the center of the cell. Area weighted average is used to get the equivalent parameters of the electric field. Therefore, the effective dielectric constant and the effective conductivity at sample point of the electric field shown in Figure 4 can be shown to be
After discretization of equations (3), (4), (5), the TM wave of the FDTD formulas can be obtained: where
The substitution of equations (12), (13), (14), (15) into equations (16), (17), (18) result in the CFDTD formulas expressed as where
To avoid numerical instabilities, the time increment must satisfy the Courant condition as follows [22]: where is velocity of light.
3. Transmitting Boundary
Transmitting boundary is a kind of a universally applicable absorption boundary condition based on the general infinite domain calculation model [33–35]. The main point of the method is to represent the one-way wave propagation along the outer normal direction of the truncated boundary of the calculated region as the superposition of a series of nonparallel plane waves [36].
In equation (28), represents any one-way fluctuation, and is the apparent velocity. There are no constraints on and ; they can be arbitrary values, and in general, the values of are different for different values of . Secondly, the concept of artificial wave velocity is proposed, and the absorption boundary condition is directly derived by general expression as shown in equation (23). Therefore, these boundary conditions are independent of the governing equations and applicable to any problem. It can be applied directly to a scalar wave problem, a vector wave problem, an isotropic problem, or an anisotropic problem. In recent years, it has been widely used in various engineering wave numerical simulation problems due to its advantages of easy computer realization and controllable precision.
Figure 5 shows the arrangement of the computational points and the spatial mesh points on the right artificial boundary. The -order multiple transmission formula (MIF) at any point on the right truncation boundary can be expressed as where in and represent the discrete value of field components at time and at space points , and is the artificial wave velocity.
According to the practical application of effectiveness, the absorption effectiveness of the second-order precision transmission boundary is better [37]. In order to improve the accuracy, we use the third-order precision transmission boundary. According to equation (24), the third orders of the transmission boundary can be expressed as
To test the absorbing ability of boundary, a simple 2D model is used. The modeled region is filled with four square meters of medium ( and mS/m), and an excitation source is started at the center of the modeled region. A circular medium ( and mS/m) with a diameter of 0.01 m is located in the middle of the model. A Ricker pulse with a 1 GHz frequency is used, with a time step and space increment chosen to be 0.01 ns and 0.005 m. All materials are assumed to be nonmagnetic () in this paper. Figure 6 shows temporal snapshots of the field component at different time steps.
(a)
(b)
(c)
(d)
It can be observed that the response decays gradually, and no fictitious reflections appear at any boundary. At time step 2.5 ns, the wave has begun to propagate, and at time step 10 ns, there is a smaller reflection at the boundary. At time step 20 ns, the repercussions in the computation domain disappear, thereby displaying the effectiveness of the absorbing boundary.
4. Numerical Examples
4.1. Inclined Layered Model
A two-layered model of the pavement is used to verify the accuracy and efficacy of the proposed algorithm by the implementation of a dipping boundary. The upper layer, representative of asphalt concrete, has and mS/m. The lower layer, representative of cement-stabilized macadam, has and mS/m. An air-earth interface is included in the model at , which is accomplished by adding a thin upper layer with and to the grid. All materials are assumed to be nonmagnetic [6]. The spatial increments and time steps are chosen to be 0.05 m and 0.1 ns, respectively. The simulation time duration is 200 ns, and the Ricker source pulse has a frequency of 100 MHz. We use MATLAB R 2014a to carry out a numerical simulation. Figure 7 shows the schematic diagram of the inclined layered medium. Transmitters and receivers are placed along the air-earth interface every 0.02 m for the reflection survey. The distance between the transmitter and receiver is 0.1 m.
First, we set the transmitter at the point (10 m, 0 m) and the receiver at the point (10.1 m, 0 m). The reflected wave as shown in Figure 8 is obtained by the proposed method and those by the FDTD method and the SPRK scheme. Figure 8 shows the single-trace signal of the layered model during the CFDTD method (black solid line), the SPRK scheme (red dash line), and the FDTD simulation (blue dotted line) for the source located at the point (10 m, 0 m). It can be observed that the FDTD method achieves almost the same level of accuracy as the SPRK scheme [6]. The amplitude of the field component equals -1-6 V/m for the primary peak at 83.1 ns, -1-6 V/m for the second peak at 83.5 ns, and 3-6 V/m for the third peak at 85.9 ns by the CFDTD method. When using the SPRK method, the amplitude of the field component equals -3-6 V/m, -2-6 V/m, and 11-6 V/m at 83.1 ns, 83.5 ns, and 85.9 ns, respectively. Moreover, when using the FDTD method at 83.1 ns, 83.5 ns, and 85.9 ns, the amplitude of the field component equals -2-6 V/m, -1.9-5 V/m, and 9-6 V/m, respectively. By comparing the amplitude of the field component of the proposed algorithm with the SPRK scheme, we can obtain the amplitude of the field component, which will be reduced by 66.7% for the primary peak at 83.1 ns, 50% for the second peak at 83.5 ns, and 72.7% for the third peak at 85.9 ns. In contrast to the proposed method and the FDTD method, the reduction of the virtual wave is 50% for the primary peak at 83.1 ns, 94.8% for the second peak at 83.5 ns, and 66.7% for the third peak at 85.9 ns. In short, it was found that the proposed algorithm can reduce the virtual waves by about 50% when compared with the FDTD method and the SPRK scheme on average.
Figure 9 shows the GPR trace profile of an inclined multilayered medium simulated using the CFDTD method. Figure 10 shows the GPR trace profile of an inclined layered medium simulated using the FDTD method. By comparing Figures 9 and 10, it can be seen that the virtual waves in the circular area marked in red are significantly reduced.
4.2. Circular Cavity Model
As a second example, a GPR-modeled profile of a section of pavement with a circular clearance is simulated. In this case, the subsurface consisted of one separated layer of and mS/m. An air-earth interface is included in the model at m by adding a thin upper layer with and to the grid, as well as two plastic pipes of 50 cm diameter to the layer with a central position fixed at (5 m, 3.5 m) and (15 m, 3.5 m). One of the circular voids is filled with water and has values of and mS/m; the other circular void is filled with air, with values of and mS/m. We assume that all materials are nonmagnetic, with a relative permeability chosen to be 1. The time step employed is 0.1 ns, the space increment is 0.05 m, and the number of iterations is 2000 steps. A Ricker pulse of unity amplitude with a central frequency of 100 MHz is used for excitation, the excitation source moved synchronously from left to right, and its position is one grid from the ground surface. Figure 11 shows the schematic diagram of the circular void model.
Figure 12 shows the single-trace signal of the layered model during the CFDTD method (black solid line), the SPRK scheme (red dash line), and the FDTD simulation (blue dotted line) for the source located at m and m. Figure 13 shows the single-trace signal of the layered model during the CFDTD method (black solid line), the SPRK scheme (red dash line), and the FDTD simulation (blue dotted line) for the source located at m and m. In Figure 12, the whole virtual waves are reduced, and at 101 ns, the virtual reflection wave from the proposed algorithm is 70.4% less than that from the SPRK method and 69.1% less than that from FDTD method. At 161 ns, the virtual wave reduction is 72.6% and 65.6%, respectively. At 179 ns, the virtual wave decreases by 51.2% and 52.3%, respectively. From Figure 13, by comparing the virtual wave produced by the proposed algorithm with the SPRK scheme, it can be found that the virtual wave is reduced by 69.5% for the primary peak at 114 ns, 72.3% for the second peak at 133 ns, and 49.6% for the third peak at 138.5 ns. Moreover, by comparing the virtual wave generated by this algorithm and FDTD algorithm, we can see that the reduction of the virtual wave is 71.4% for the primary peak at 114 ns, 68.1% for the second peak at 133 ns, and 51.2% for the third peak at 138.5 ns. In conclusion, more than 50% of virtual waves can be reduced by using the CFDTD method.
Figure 14 shows the GPR trace profile of a circular void medium model simulated by the CFDTD method, with Figure 15 showing the GPR trace profile of a circular void medium model simulated by the FDTD method. It can be seen from the circular area marked in red in Figures 14 and 15 that most of the virtual waves are basically eliminated.
As can be observed from Figures 14 and 15, excellent agreement is obtained, as two hyperbolic diffraction waves are produced in both models. The diffraction waves of the circular voids filled with water are significantly greater than the diffraction waves of the circular void filled with air. This effect is due to the fact that the relative dielectric constant of water is larger than the dielectric constant of air. The comparison of the GPR trace profile of the CFDTD, SPRK, and FDTD methods has shown that the CFDTD method has eliminated the virtual waves generated by the staircase approximation.
4.3. The Three-Layered Pavement Model
As a third example, the GPR profile of the three-layered pavement model with structural damages is simulated. In this case, the incident wave is a Ricker waveform of unity amplitude and 1 GHz central frequency and the incident wave is considered as a line source. As shown in Figure 16, the pavement structure that suffers structural damages consists of three layers. A 15 mm width crack exists in the first layer. One circular void with a radius of 0.02 m filled with water is located inside the second layer. A rectangular no-dense region with a width of 50 cm and a thickness of 20 cm is located inside the third layer. The electric constants of the materials are also shown in Figure 16; however, the electric constants of the no-dense region need to be specified. We assume that the porosity in this region is 20%, the electric constants of two-tenths of the nodal points in this region are defined as that of the air, and these nodes are randomly distributed. The electric constants of the other nodes are the same as that of the third layer.
The transmitter and receiver move once every 0.02 m to transmit and receive reflected waves. The transmitter is 0.1 m apart from the receiver. The increment of space interval and time step are chosen as 0.005 m and 0.01 ns, respectively. The GPR trace profiles of this model simulated by the FDTD method and the CFDTD scheme are plotted in Figures 17 and 18 accordingly.
As can be observed from these figures, most of the virtual waves can be reduced by using the CFDTD method inside the red box; however, because the void is too small, the effect of conformal treatment on radar reflection is not too obvious.
4.4. The Measurement Model
In order to verify the effect of the inclined layered model and the circular cavity model by the CFDTD method, we carry out the actual radar wave detection as the third example. We dug an underground tunnel that was 1 m long, 1 m wide, and 1 m high in the laboratory. The tunnel was filled with dry sand. A plastic pipe with a diameter of 25.47 mm was buried at a depth of 15.2 cm in the model, and the length of the plastic pipe was 1 m. Figure 19 shows a schematic diagram of the experimental model.
Figure 20 shows the single-trace signal of a circular void medium model at a distance of 0.45 m. The IDS GPR is adopted with a central frequency of 2.5 GHz, 512 sampling points, and a time window of 0.012 ns. The moving step length of the antenna is 0.004 m. Sending and receiving are working at the same time. The measured waveform and simulated waveform at a distance of 0.45 m are compared and shown in Figure 20 (the mesh size is 0.008 m). Figure 21 shows the absolute error of the amplitude of the field component by the CFDTD method, FDTD method, and the actual detection data. The maximum absolute difference between the method is no more than 0.0008 V/m.
Comparing the actual measured data with the ballistic wave diagram of CFDTD calculation method, it can be seen that the simulated data of the CFDTD method is basically consistent with the measured data. The absolute error of the CFDTD method is smaller than the FDTD method.
5. Conclusion
Based on the CFDTD method, we first validated the accuracy of conformal object modeling, and then we established a refined numerical model of complex underground effective structures. Furthermore, effective dielectric parameters were adopted to precisely describe the curved boundary of an irregular object, and three electric models of the inclined layered medium, the circular medium, and the three-layered pavement model were established to simulate the reflection waveform of the GPR detection. Finally, the relevant detected data were compared with the result of numerical results as shown in the fourth example. The comparison clearly shows that the simulated results of the proposed method are fully consistent with the measured values. The absolute error of the CFDTD method is smaller than the FDTD method. Comparing the results of calculations for the FDTD, SPRK, and CFDTD techniques, the virtual waves of the GPR waveform are obviously reduced by at least 50% with the CFDTD method. The CFDTD method is closer to the actual radar detection data, and it can effectively reduce the virtual wave generated by the simulation. The comparison of the CFDTD method with other methods like the pseudospectral time-domain (PSTD) and the spectral element time-domain (SETD) methods will be applied to highlight the advantages of the proposed CFDTD method so as to achieve a tradeoff between accuracy and calculation efficiency of later research.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
This research was funded by the National Key Research and Development Program of China (No. 2017YFC1501204) and the National Natural Science Foundation of China (Nos. 51678536 and 41404096). This project was also supported by the Department of Education’s Production-Study-Research Combined Innovation Funding—“Blue Fire Plan (Huizhou)” (CXZJHZ01742), the Program for Science and Technology Innovation Talents in Universities of Henan Province (Grant No. 19HASTIT043), and the Outstanding Young Talent Research Fund of Zhengzhou University (1621323001). The authors are grateful for these financial supports.