Efficacious GPR Implementations of Z-Transform-Based Hybrid LOD-FDTD with Subgridding Scheme: Theoretical Formalism and Numerical Study

Abstract

Ground penetrating radar (GPR) forward modeling is one of the core geophysical research topics and also the primary task of simulating ground penetrating radar system. It is a process of simulating the propagation laws and characteristics of electromagnetic waves in simulated space when the distribution of internal parameters in the exploration region is known. And the finite-difference-time-domain (FDTD) method has the characteristics of simulating the space-time transient evolution of electromagnetic wave, whose numerical method is simple and easy to program, so it has become one of the most extensively utilized methods in GPR forward modeling. It is generally known that the conventional FDTD approach requires finer uniform Yee cell all the time to produce satisfactory accuracies from numerical simulations of the GPR. However, the smaller temporal incremental has to be adopted due to the lower spatial incremental, which would dramatically weaken the advantage of the FDTD method. To solve this issue, the subgridding-technique-based hybrid local-one-dimensional FDTD (LOD-FDTD) is applied in this work to modeling the classical GPR scenarios. In this method, the unconditional-stable LOD-FDTD is employed in the fine-grid domain, while the traditional FDTD is used in the coarse-grid domain, which could avoid the oversampling problem in the local domain if the uniform fine-grid scheme is adopted. Meanwhile due to the unconditional stability of the LOD-FDTD, the larger time step, derived from the coarse grid which satisfies the Courant-Friedrichs-Lewy (CFL) stability condition, could be utilized in the whole domain so that the long-time interpolation process could be circumvented. Additionally, the proposed approach could be arbitrarily adjusted by means of different ratio of both coarse- and fine-grid, and hence it holds much higher generality. As compared with the auxiliary differential equation (ADE) technique, the Z-transform method is integrated into FDTD methods for modeling multi-pole Debye-based dispersive media in this method, resulting in more direct numerical implementations and fewer computing steps. Finally, three different classical GPR problems are carried out to validate accuracies and efficiencies of the proposed method.

1. Introduction

Ground penetrating radar (GPR), consisting of antenna, control unit, display unit, storage unit and collateral equipments, is a non-destructive measurement technology in geophysical exploration methods [1,2,3]. It applies electromagnetic (EM) waves to locate targets or interfaces concealed in visually opaque objects or underground media. The geological radar transmits high-frequency EM waves through the transmitting antenna, and receives reflected EM waves by the receiving antenna. With the help of the reflected EM wave data returned by the receiving antenna, the specific parameter characteristics could be extracted, and then combined with certain data processing and analysis means, therefore, the spatial position and distribution of different media can be inferred [4,5,6]. When the GPR method is used for detection, it mainly depends on the comparison of the returned data to infer and analyze the attributes of the detected target, and the reference data for comparison mainly derives from field sampling and numerical simulation. The former relies on experiments to finish data collection, which increases the complexity to implement. Numerical simulation could easily obtain comparative data through establishing relevant environmental models, including soil, rock, groundwater distribution and abnormal targets with high dielectric properties [7,8,9,10,11]. This method can save the cost of human, material and financial resources in the process of experiment and system improvement, and accelerate the progress of theoretical and experimental research. This numerical simulation method could be denoted as the GPR forward simulation, which is one of the main research contents of geophysics [12,13]. In GPR forward modeling, the finite-difference time-domain (FDTD) is a commonly employed numerical method at present [14,15,16], since the FDTD method has the characteristics of simple concept, flexibility and strong versatility [17,18]. Simultaneously, the detection scene composed of the detected target and the background soil always keep the characteristics of non-uniformity, electrical characteristics, dispersion characteristics and anisotropy [19,20,21,22]. Fortunately, the EM characteristics of these structures could be easily simulated by the Yee-cell-based FDTD method which has been applied to the research of GPR system for a long time.

Giannopoulos et al. developed the GPRMax, a kind of GPR forward modeling software based on the FDTD method [23]; Literature [24] studied the forward simulation and GPR imaging in tunnel lining leakage detection and grouting evaluation; and Shaari studied the polarity of GPR antenna targets [25]; with although the FDTD method has numerous advantages in the GPR forward simulation, for the EM model containing complex high dielectric constant, medium with dispersion characteristics and small electrical structure, the FDTD method is required to adopt a finer grid to well simulate the propagation characteristics of EM wave in the region with drastic changes in amplitude, phase, and other information, and ensures that the numerical method has low numerical dispersion error and anisotropy error [26,27]. Applying the finer-grid generation method would greatly increase the memory occupation of the computer. Moreover, the FDTD method uses the explicit difference strategy, and its time step is limited by the spatial cell size, namely it needs to meet the Courant-Friedrichs-Lewy (CFL) stability condition [28,29], which makes the GPR simulation system with fine-grid discretization have poor efficiency, and hence significantly affects the computational advantage of the FDTD method.

With the aim of addressing this problem, to the best of our knowledge, the subgridding technology by adopting multiscale meshing scheme is a quite effective means [30,31]. The strategy we apply is that the fine-grid discretization is adopted in the domain where the amplitude of EM field changes violently due to the non-uniformity of medium properties over a small scale, and coarse-grid discretization is employed in other domains, which could circumvent the oversampling of the whole computational domain, and hence could considerably improve the computational efficiency and reduce the memory. The essence of subgridding technology is the temporal and spatial information exchange of EM field components between coarse and fine grids [32]. As a result of different cell sizes for coarse and fine grids, each subdomain oughts to meet its corresponding CFL stability condition. In order to obtain stable FDTD solutions, the computational schemes of early subgridding technology mainly fall into the following two categories: (1) the whole computational domain, including coarse- and fine-grid domains, adopts a unified time step that meets the CFL stability condition of fine grid for numerical calculation [33]. (2) The coarse- and fine-grid domains are numerically calculated using time steps that meet their respective CFL conditions [34,35]. The former method could avoid interpolation process in time between coarse and fine grids, and only requires spatial interpolation, which is easier to obtain numerically. However, on account of the small time step used in the entire computational domain, this method exhibits poor efficiency. Compared with the first method, the second one is more efficient, but it requires double interpolation operation in both time and space, hence the numerical implementation turns out to be more sophisticated.

In order to circumvent interpolation processing in time, and at the same time, the whole computational domain could adopt a larger time step, some hybrid subgridding technologies have been proposed within the fine-grid domain, weakly conditionally-stable and unconditionally-stable FDTD methods, such as the hybrid-implicit-explicit (HIE) method [36], the alternating-direction-implicit FDTD (ADI-FDTD) method [37,38], the Crank-Nicolson FDTD (CN-FDTD) method [39], and the spatial filtering FDTD method [40] with all of them are applied to calculate the fine-grid domain. Because these methods have the characteristics of expandable stability conditions or unconditional stability, the fine-grid region could adopt the time step that meets the coarse-grid CFL condition so that the whole calculation could choose a larger time step for numerical calculations, which improves the computational efficiency, avoids the interpolation processing in time, and reduces complexity of the numerical implementation. However, the hybrid HIE subgridding method illustrated in [36] did not analyze the treatment of dispersive media. Researchers [37,38] only gave the application of ADI-FDTD in two-dimensional (2D) subgridding technology, and then the computational efficiency of ADI-FDTD was not high. Next, the hybrid subgridding technology introduced in [39] adopted the CN-FDTD method, which also had the problem of high computational complexity. The hybrid subgridding technology presented in [40] employed the spatial filtering FDTD method. Although the spatial filtering FDTD method retained the explicit iterative characteristics of the traditional FDTD method, its computational accuracy was exceedingly affected by the dielectric properties of the materials in the computational region. As an improvement of ADI-FDTD method, the local-one-dimensional FDTD (LOD-FDTD) method has simpler implementation and more efficient computational performance [41,42]. Literature [43] introduced the hybrid LOD-FDTD subgridding technology for simulating periodic metal nanoparticle arrays, but only introduced the calculation scheme in the 2D case. Concurrently, most of the above work only studied the interpolation calculation strategy with odd coarse- and fine-grid proportion, and did not propose a general subgridding technology scheme that could be divided into arbitrary proportions.

In view of this, we propose a three-dimensional (3D) hybrid LOD-FDTD method with subgridding technology for numerical GPR simulations. The LOD-FDTD method is employed for numerical calculation in fine-grid domain, and the traditional FDTD method is applied to calculation in coarse-grid region. Meanwhile, we introduce a more general hybrid multiscale grid processing technology, which could realize arbitrary-scale fine-grid generation in local coarse-grid region and perceive high versatility. Moreover, aiming at the dispersion properties of soil and detection targets, the LOD-FDTD method based on the Z-transform technique (ZT-LOD-FDTD) is proposed to deal with the multi-level Debye dispersive media. Compared with the ADE method [44], the Z-transform method [45] has no further use for retaining abundant values of the previous moments, occupies less memory, and hence its numerical implementation is more direct. Furthermore, the convolutional perfectly matched layer (CPML) developed based on the Z-transform method is utilized to truncate the soil dispersion domain [46].

The rest of this paper is organized as follows: in the first part of Section 2, for the multi-level Debye dispersion model, the numerical iterative formulations of the ZT-LOD-FDTD method in the frequency domain are derived. In the second part of the Section 2, the spatial modeling technology that could realize arbitrary-scale generation of coarse- and fine-grid is proposed based on the LOD-FDTD method. In the Section 3, the correctness and effectiveness of the proposed method are verified through the numerical simulation of several common classical GPR scenarios. Finally, we draw the conclusion in the Section 4.

2. Materials and Methods

2.1. LOD-FDTD Method Based on the Z-Transform Technique

The splitting-based LOD scheme is applied to the Maxwell’s equations [41,42], and the formulations could be written in the form of matrix as

$$\left(I_6-\frac{\Delta t}{2}{A}_1\right){\mathit{u}}^{n+\frac{1}{2}}=\left(I_6+\frac{\Delta t}{2}{A}_1\right){\mathit{u}}^n,$$

(1)

$$\left(I_6-\frac{\Delta t}{2}{A}_2\right){\mathit{u}}^{n+1}=\left(I_6+\frac{\Delta t}{2}{A}_2\right){\mathit{u}}^{n+\frac{1}{2}},$$

(2)

where

$$\mathit{u}={\left[D_x,D_y,D_z,H_x,H_y,H_z\right]}^T,$$

(3)

$$A_1=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& \frac{\partial }{\partial y}\\ 0& 0& 0& \frac{\partial }{\partial z}& 0& 0\\ 0& 0& 0& 0& \frac{\partial }{\partial x}& 0\\ 0& \frac{1}{\mu }\frac{\partial }{\partial z}& 0& 0& 0& 0\\ 0& 0& \frac{1}{\mu }\frac{\partial }{\partial x}& 0& 0& 0\\ \frac{1}{\mu }\frac{\partial }{\partial \mathrm{y}}& 0& 0& 0& 0& 0\end{array}\right]A_2=\left[\begin{array}{cccccc}0& 0& 0& 0& -\frac{\partial }{\partial z}& 0\\ 0& 0& 0& 0& 0& -\frac{\partial }{\partial x}\\ 0& 0& 0& -\frac{\partial }{\partial y}& 0& 0\\ 0& 0& \frac{-1}{\mu }\frac{\partial }{\partial y}& 0& 0& 0\\ \frac{-1}{\mu }\frac{\partial }{\partial z}& 0& 0& 0& 0& 0\\ 0& \frac{-1}{\mu }\frac{\partial }{\partial x}& 0& 0& 0& 0\end{array}\right],$$

(4)

At the same time, the constitutive relation between the D and E is defined as

$$\mathit{D}\left(\omega \right)={\epsilon }_0{\epsilon }_r\left(\omega \right)\mathit{E}\left(\omega \right),$$

(5)

For the multipole Debye dispersive media, the constitutive relation of Equation (5) could be expressed as

$$\mathit{D}(\omega )={\epsilon }_0({\epsilon }_{\infty }+\frac{\sigma }{j\omega {\epsilon }_0}+{\displaystyle \sum _{p=1}^{N_p}}\frac{\Delta {\epsilon }_p}{1+j\omega {\tau }_p})+\mathit{E}(\omega )$$

(6)

where ε₀ is denoted as the dielectric constant in the vacuum, σ is the electrical conductivity, ∆ε_p = ε_s,p − ε_∞,p is the difference between the relative dielectric constant in the zero frequency and that in the infinite frequency. τ_p is the pole relaxation time, N_p is the pole number in the sensitivity response.

Transforming the Equation (6) from the frequency domain to the Z domain, we have

$$\mathit{D}(z)={\epsilon }_0{\epsilon }_{\infty }\mathit{E}(z)+{\displaystyle \sum _{p=1}^{N_p}{\mathcal{F}}_{\mathrm{p}}(z)}+\mathcal{R}(z),$$

(7)

where ${\mathcal{F}}_p$ and $\mathcal{R}$ are the auxiliary variables

$${\mathbf{ℱ}}_p(z)=\frac{{\epsilon }_0\Delta \epsilon \Delta t/(2{\tau }_p)}{1-e^{-\Delta t/(2{\tau }_p)}{z}^{-1/2}}\mathit{E}(z),$$

(8)

$$\mathbf{ℛ}(z)=\frac{\sigma \Delta t/2}{1-z^{-1/2}}\mathit{E}(z),$$

(9)

We convert Equations (7)–(9) from the Z domain to the time domain, namely replacing $\mathit{P}(z)$ and $z^{-1/2}\mathit{P}(z)$ with ${\mathit{P}}^n$ and ${\mathit{P}}^{n-1/2}$, ($\mathit{P}=\left(\mathit{D},\mathit{E},{\mathit{ℱ}}_p,\mathit{ℛ}\right)$), and then incorporate them into the Equations (1) and (2) so as to obtain the iterative formulations of the ZT-LOD-FDTD of the sub time step from n to n+1/2:

$${\epsilon }_0{\epsilon }_{\infty }{\mathit{E}}^{n+\frac{1}{2}}-\frac{\Delta t^2}{4{\mu }_0}AA{\mathit{E}}^{n+\frac{1}{2}}+{\displaystyle \mathit{\sum }_{p=1}^{N_p}{\mathit{ℱ}}_p^{n+\frac{1}{2}}}+{\mathit{ℛ}}^{n+\frac{1}{2}}={\epsilon }_0{\epsilon }_{\infty }{\mathit{E}}^n+{\displaystyle \sum _{p=1}^{N_p}{\mathit{ℱ}}_p^n}+{\mathit{ℛ}}^n+\Delta tB{\mathit{H}}^n+\frac{\Delta t^2}{4{\mu }_0}AA{\mathit{E}}^n,$$

(10)

$${\mathit{H}}^{n+\frac{1}{2}}={\mathit{H}}^n+\frac{1}{2}\frac{\Delta t}{{\mu }_0}\left(A{\mathit{E}}^n-A{\mathit{E}}^{n+\frac{1}{2}}\right),$$

(11)

For the sub time step from n + 1/2 to n + 1:

$${\epsilon }_0{\epsilon }_{\infty }{\mathit{E}}^{n+1}-\frac{\Delta t^2}{4{\mu }_0}BB{\mathit{E}}^{n+1}+{\displaystyle \sum _{p=1}^{N_p}{\mathit{ℱ}}_p^{n+1}}+{\mathcal{R}}^{n+1}={\epsilon }_0{\epsilon }_{\infty }{\mathit{E}}^{n+\frac{1}{2}}+{\displaystyle \sum _{p=1}^{N_p}{\mathcal{F}}_p^{\mathrm{n}+\frac{1}{2}}}+{\mathcal{R}}^{\mathrm{n}+\frac{1}{2}}-\Delta tA{\mathit{H}}^{n+\frac{1}{2}}+\frac{\Delta t^2}{4{\mu }_0}BB{\mathit{E}}^{n+\frac{1}{2}},$$

(12)

$${\mathit{H}}^{n+1}={\mathit{H}}^{n+\frac{1}{2}}-\frac{1}{2}\frac{\Delta t}{{\mu }_0}\left(B{\mathit{E}}^{n+1}-B{\mathit{E}}^{n+\frac{1}{2}}\right),$$

(13)

where

$${\mathit{ℱ}}_p^{n+\frac{1}{2}}=e^{-\frac{\Delta t}{2{\tau }_p}}{\mathit{ℱ}}_p^n+\frac{\Delta \epsilon \Delta t}{2{\tau }_p}{\mathit{E}}^{n+\frac{1}{2}},$$

(14)

$${\mathit{ℛ}}^{n+\frac{1}{2}}={\mathit{ℛ}}^n+\frac{\sigma \Delta t}{2}{\mathit{E}}^{n+\frac{1}{2}},$$

(15)

$$A=\left(\begin{array}{ccc}0& \frac{\partial }{\partial z}& 0\\ 0& 0& \frac{\partial }{\partial x}\\ \frac{\partial }{\partial y}& 0& 0\end{array}\right)\hspace{1em}B=\left(\begin{array}{ccc}0& 0& \frac{\partial }{\partial y}\\ \frac{\partial }{\partial z}& 0& 0\\ 0& \frac{\partial }{\partial x}& 0\end{array}\right),$$

(16)

Afterward, Equations (14) and (15) are incorporated into Equation (10), the iterative formulation of E in the moment of n + 1/2 could be obtained. Taking x direction as an example, we have

$$\begin{array}{l}a{E}_x^{n+\frac{1}{2}}(i+\frac{1}{2},j-1,k)+b{E}_x^{n+\frac{1}{2}}(i+\frac{1}{2},j,k)+c{E}_x^{n+\frac{1}{2}}(i+\frac{1}{2},j+1,k)\\ =w_1{E}_x^n(i+\frac{1}{2},j,k)+{\displaystyle \sum _{p=1}^{N_p}{w}_{2p}{\mathcal{F}}_{xp}^n(i+\frac{1}{2},j,k)}\\ +w_3\left(E_x^n(i+\frac{1}{2},j-1,k)+E_x^n(i+\frac{1}{2},j+1,k)\right)\\ +\frac{\Delta t}{{\epsilon }_0\Delta y}\left(H_z^n(i+\frac{1}{2},j+\frac{1}{2},k)-H_z^n(i+\frac{1}{2},j-\frac{1}{2},k)\right)\end{array}$$

(17)

where $a=c=-\frac{\Delta t^2}{4{\mu }_0{\epsilon }_0\Delta y^2}$, $b={\epsilon }_{\infty }+{\displaystyle \sum _{p=1}^{N_p}\frac{\Delta {\epsilon }_p\Delta t}{2{\tau }_p}}+\frac{\Delta t^2}{2{\mu }_0{\epsilon }_0\Delta y^2}+\frac{\sigma }{2\Delta t}$, $w_1={\epsilon }_{\infty }-\frac{\Delta t^2}{2{\mu }_0{\epsilon }_0\Delta y^2}$, $w_{2p}=\frac{e^{-\Delta t/2{\tau }_p}}{{\epsilon }_0}$, $w_3=\frac{\Delta t^2}{4{\mu }_0{\epsilon }_0\Delta y^2}$.

In the same way, the iterative formulation of E_x in the moment of n+1 could be also obtained

$$\begin{array}{l}a{E}_x^{n+1}(i+\frac{1}{2},j,k-1)+b{E}_x^{n+1}(i+\frac{1}{2},j,k)+c{E}_x^{n+1}(i+\frac{1}{2},j,k+1)\\ =w_1{E}_x^{n+\frac{1}{2}}(i+\frac{1}{2},j,k)+{\displaystyle \sum _{p=1}^{N_p}{w}_{2p}{\mathcal{F}}_{xp}^{n+\frac{1}{2}}(i+\frac{1}{2},j,k)}\\ +w_3\left(E_x^{n+\frac{1}{2}}(i+\frac{1}{2},j,k-1)+E_x^{n+\frac{1}{2}}(i+\frac{1}{2},j,k+1)\right)\\ -\frac{\Delta t}{{\epsilon }_0\Delta z}\left(H_y^{n+\frac{1}{2}}(i+\frac{1}{2},j,k+\frac{1}{2})-H_y^{n+\frac{1}{2}}(i+\frac{1}{2},j,k-\frac{1}{2})\right)\end{array}$$

(18)

where $a=c=-\frac{\Delta t^2}{4{\mu }_0{\epsilon }_0\Delta z^2}$, $b={\epsilon }_{\infty }+{\displaystyle \sum _{p=1}^{N_p}\frac{\Delta {\epsilon }_p\Delta t}{2{\tau }_p}}+\frac{\Delta t^2}{2{\mu }_0{\epsilon }_0\Delta z^2}+\frac{\sigma }{2\Delta t}$, $w_1={\epsilon }_{\infty }-\frac{\Delta t^2}{2{\mu }_0{\epsilon }_0\Delta z^2}$, $w_{2p}=\frac{e^{-\Delta t/2{\tau }_p}}{{\epsilon }_0}$, $w_3=\frac{\Delta t^2}{4{\mu }_0{\epsilon }_0\Delta z^2}$.

Since Equation (15) is incorporated into Equations (10) and (12), there is no auxiliary variable $\mathcal{R}$ in Equations (17) and (18). In the calculation process, only the auxiliary variable ${\mathcal{F}}_{xp}$ needs to be updated iteratively (Equation (14)), which need the value at the previous time and the value of the electric field at the current time.

2.2. Hybrid LOD-FDTD with the Subgridding Scheme

In this work, the hybrid computational scheme of multiscale grid is proposed by combining the traditional FDTD method with the LOD-FDTD method. The EM components in the coarse-grid region are numerically calculated through the traditional FDTD method, while the EM components in the fine-grid domain are numerically solved by the LOD-FDTD method (such as Equations (17) and (18)). Depending on the unconditional stability of the LOD-FDTD method, the fine-grid region could adopt the time step that meets the CFL stability conditions of the coarse-grid time step, therefore, the coarse-grid region only is required to be interpolated in space, not in time. The purpose of spatial interpolation is that the electric field value e of fine grid on the interface of coarse and fine grids cannot be solved iteratively by the LOD-FDTD method. As shown in Figure 1, the electric field value E of coarse grid on the interface needs to be interpolated to the electric field value e of fine grid on the interface so that the internal EM field components of fine grid domain could be solved collaboratively. In 3D computational space, there are six interfaces in the coarse and fine grid domains. For the convenience of description, we take the electric field value e_x on the coarse and fine grid interface seen in Figure 1 as an example to introduce the relevant interpolation processing. Meanwhile, the similar method could be applied to other interfaces.

As shown in Figure 2, it is the diagram of electric field distribution between coarse and fine grids on the interface $j_F=j_A$($j_f=1$) of Figure 1, the subscript F is denoted as the coarse-grid coordinate, the f is defined as the fine-grid coordinate. Taking the interpolation processing of fine-grid electric field from $e_x(i_f+\frac{1}{2}+s_1,1,1)$ to $e_x(i_f+\frac{1}{2}+s_1,1,R\times (k_A-k_B))$ in the box of Figure 2 as an example, where $s_1$ is the positive integer, $1\le s_1\le R$, R is the ratio between coarse and fine grids, $i_F$ is expressed as the grid coordinate ($i_A\le i_F<i_B$) of the coarse grid along x direction on the interface. i_f is defined as the grid coordinate (i_f = R × (i_F − i_A). (i_A, i_B)) of the fine grid along x direction on the interface, (j_A, j_B), (k_A, k_B) are respectively the front and rear interfaces in the x, y, and z directions of coarse-grid domain. The spatial interpolation formulations shown below could be utilized to obtain the fine-grid electric fields from $e_x(i_f+\frac{1}{2}+s_1,1,1)$ to $e_x(i_f+\frac{1}{2}+s_1,1,R\times (k_A-k_B))$ in the box.

When $s_1\le \frac{R}{2}$:

$$\left[\begin{array}{c}{e}_x(i_f+\frac{1}{2}+s_1,1,1)\\ e_x(i_f+\frac{1}{2}+s_1,1,2)\\ ⋮\\ e_x(i_f+\frac{1}{2}+s_1,1,R(k_B-k_A))\end{array}\right]=M_1\left[\begin{array}{c}{E}_x(i_F+\frac{1}{2},j_A,k_A)\\ E_x(i_F+\frac{1}{2},j_A,k_A+1)\\ ⋮\\ E_x(i_F+\frac{1}{2},j_A,k_B)\end{array}\right]+M_2\left[\begin{array}{c}{E}_x(i_F-\frac{1}{2},j_A,k_A)\\ E_x(i_F-\frac{1}{2},j_A,k_A+1)\\ ⋮\\ E_x(i_F-\frac{1}{2},j_A,k_B)\end{array}\right],$$

(19)

When $\frac{R}{2}<s_1\le R$:

$$\left[\begin{array}{c}{e}_x(i_f+\frac{1}{2}+s_1,1,1)\\ e_x(i_f+\frac{1}{2}+s_1,1,2)\\ ⋮\\ e_x(i_f+\frac{1}{2}+s_1,1,R(k_B-k_A))\end{array}\right]=M_3\left[\begin{array}{c}{E}_x(i_F+\frac{1}{2},j_A,k_A)\\ E_x(i_F+\frac{1}{2},j_A,k_A+1)\\ ⋮\\ E_x(i_F+\frac{1}{2},j_A,k_B)\end{array}\right]+M_4\left[\begin{array}{c}{E}_x(i_F+\frac{3}{2},j_A,k_A)\\ E_x(i_F+\frac{3}{2},j_A,k_A+1)\\ ⋮\\ E_x(i_F+\frac{3}{2},j_A,k_B)\end{array}\right],$$

(20)

where

$$M_n=\left[\begin{array}{cccccc}{P}_n& Q_n& 0& \cdots & 0& 0\\ 0& P_n& Q_n& \cdots & 0& 0\\ ⋮& ⋮& \ddots & \ddots & ⋮& ⋮\\ 0& 0& 0& \cdots & P_n& Q_n\end{array}\right],$$

(21)

$$P_1=\left[\begin{array}{c}\frac{(R+2s_1-1)R}{2R^2}\\ \frac{(R+2s_1-1)(R-1)}{2R^2}\\ ⋮\\ \frac{(R+2s_1-1)}{2R^2}\end{array}\right]Q_1=\left[\begin{array}{c}0\\ \frac{(R+2s_1-1)}{2R^2}\\ ⋮\\ \frac{(R+2s_1-1)(R-1)}{2R^2}\end{array}\right],$$

(22)

$$P_2=\left[\begin{array}{c}\frac{(R-2s_1+1)R}{2R^2}\\ \frac{(R-2s_1+1)(R-1)}{2R^2}\\ ⋮\\ \frac{(R-2s_1+1)}{2R^2}\end{array}\right]Q_2=\left[\begin{array}{c}0\\ \frac{(R-2s_1+1)}{2R^2}\\ ⋮\\ \frac{(R-2s_1+1)(R-1)}{2R^2}\end{array}\right],$$

(23)

$$P_3=\left[\begin{array}{c}\frac{(3R-2s_1+1)R}{2R^2}\\ \frac{(3R-2s_1+1)(R-1)}{2R^2}\\ ⋮\\ \frac{(3R-2s_1+1)}{2R^2}\end{array}\right]Q_3=\left[\begin{array}{c}0\\ \frac{(3R-2s_1+1)}{2R^2}\\ ⋮\\ \frac{(3R-2s_1+1)(R-1)}{2R^2}\end{array}\right],$$

(24)

$$P_4=\left[\begin{array}{c}\frac{(2s_1-R-1)R}{2R^2}\\ \frac{(2s_1-R-1)(R-1)}{2R^2}\\ ⋮\\ \frac{(2s_1-R-1)}{2R^2}\end{array}\right]Q_4=\left[\begin{array}{c}0\\ \frac{(2s_1-R-1)}{2R^2}\\ ⋮\\ \frac{(2s_1-R-1)(R-1)}{2R^2}\end{array}\right],$$

(25)

It should be noted in Figure 2 that the ratio of the coarse- and fine-grid generation is the odd. Similarly, the case of the even ratio of coarse- and fine-grid generation could be obtained in the same way.

After adopting the spatial interpolation processing, the iterative formulations of electric field $e_x(i_f+\frac{1}{2}+s_1,2,k_f+s_2)$ on the surface j_f = 2 of fine-grid domain need to be modified, taking the equation (Equation (17)) in the moment of n + 1/2 as an example, we have

When $s_1\le \frac{R}{2}$:

$$\begin{array}{l}b{e}_x^{n+\frac{1}{2}}(i_f+\frac{1}{2}+s_1,2,k_f+s_2)+c{e}_x^{n+\frac{1}{2}}(i_f+\frac{1}{2}+s_1,3,k_f+s_2)\\ =w_1{e}_x^n(i_f+\frac{1}{2}+s_1,2,k_f+s_2)+{\displaystyle \sum _{p=1}^{N_p}{w}_{2p}{\mathcal{F}}_{xp}^n(i_f+\frac{1}{2}+s_1,2,k_f+s_2)}\\ +w_3\left(e_x^n(i_f+\frac{1}{2}+s_1,1,k_f+s_2)+e_x^n(i_f+\frac{1}{2}+s_1,3,k_f+s_2)\right)\\ +\frac{\Delta t}{{\epsilon }_0\Delta y}\left(h_z^n(i_f+\frac{1}{2}+s_1,\frac{3}{2},k+s_2)-h_z^n(i_f+\frac{1}{2}+s_1,\frac{1}{2},k_f+s_2)\right)\\ -a\left[\begin{array}{l}\frac{(R+2s_1-1)(R-s_2+1)}{2R^2}{E}_x^n(i_F+\frac{1}{2},j_A,k_F)+\frac{(R-2s_1+1)(R-s_2+1)}{2R^2}{E}_x^n(i_F-\frac{1}{2},j_A,k_F)\\ +\frac{(R+2s_1-1)(s_2-1)}{2R^2}{E}_x^n(i_F+\frac{1}{2},j_A,k_F+1)+\frac{(R-2s_1+1)(s_2-1)}{2R^2}{E}_x^n(i_F-\frac{1}{2},j_A,k_F+1)\end{array}\right]\end{array}$$

(26)

when $\frac{R}{2}<s_1\le R$:

$$\begin{array}{l}b{e}_x^{n+\frac{1}{2}}(i_f+\frac{1}{2}+s_1,2,k_f+s_2)+c{e}_x^{n+\frac{1}{2}}(i_f+\frac{1}{2}+s_1,3,k_f+s_2)\\ =w_1{e}_x^n(i_f+\frac{1}{2}+s_1,2,k_f+s_2)+{\displaystyle \sum _{p=1}^{N_p}{w}_{2p}{\mathcal{F}}_{xp}^n(i_f+\frac{1}{2}+s_1,2,k_f+s_2)}\\ +w_3\left(e_x^n(i_f+\frac{1}{2}+s_1,1,k_f+s_2)+e_x^n(i_f+\frac{1}{2}+s_1,3,k_f+s_2)\right)\\ +\frac{\Delta t}{{\epsilon }_0\Delta y}\left(h_z^n(i_f+\frac{1}{2}+s_1,\frac{3}{2},k_f+s_2)-h_z^n(i_f+\frac{1}{2}+s_1,\frac{1}{2},k_f+s_2)\right)\\ -a\left[\begin{array}{l}\frac{(3R-2s_1+1)(R-s_2+1)}{2R^2}{E}_x^n(i_F+\frac{1}{2},j_A,k_F)+\frac{(2s_1-R-1)(R-s_2+1)}{2R^2}{E}_x^n(i_F+\frac{3}{2},j_A,k_F)\\ +\frac{(3R-2s_1+1)(s_2-1)}{2R^2}{E}_x^n(i_F+\frac{1}{2},j_A,k_F+1)+\frac{(2s_1-R-1)(s_2-1)}{2R^2}{E}_x^n(i_F+\frac{3}{2},j_A,k_F+1)\end{array}\right]\end{array}$$

(27)

where $e_x^n(i_f+\frac{1}{2}+s_1,1,k_f+s_2)$ is derived from Equations (19) and (20). s₂ is defined as the positive integer, $1\le s_2\le R$, $k_F$ is expressed as the grid coordinate ($k_A\le k_F<k_B$) along z direction on the interface of coarse grid. k_f is the grid coordinate ($k_f=R\times (k_F-k_A)$) along z direction on the interface of fine grid.

After finishing the iteration of EM field components in the fine-grid region, The EM field components of coarse-grid domain in the fine-grid domain need to be interpolated in order that the EM information exchange between coarse and fine grids could be achieved. Taking x direction as an instance, we have

If R is the odd:

$$E_x^{n+\frac{1}{2}}\left(i_F+\frac{1}{2},j_F,k_F\right)=e_x^{n+\frac{1}{2}}\left(i_f+\frac{R}{2},j_f+1,k_f+1\right),$$

(28)

$$H_x^{n+1}\left(i_F,j_F+\frac{1}{2},k_F+\frac{1}{2}\right)=h_x^{n+1}\left(i_f+1,j_f+\frac{R}{2},k_f+\frac{R}{2}\right),$$

(29)

If R is the even:

$$E_x^{n+\frac{1}{2}}\left(i_F+\frac{1}{2},j_F,k_F\right)=\left[e_x^{n+\frac{1}{2}}\left(i_f+\frac{R+1}{2},j_f+1,k_f+1\right)+e_x^{n+\frac{1}{2}}\left(i_f+\frac{R+3}{2},j_f+1,k_f+1\right)\right]/2,$$

(30)

$$\begin{array}{ll}{H}_x^{n+1}\left(i_F,j_F+\frac{1}{2},k_F+\frac{1}{2}\right)& =[h_x^{n+1}\left(i_f+1,j_f+\frac{R+1}{2},k_f+\frac{R+1}{2}\right)\\ & +h_x^{n+1}\left(i_f+1,j_f+\frac{R+3}{2},k_f+\frac{R+1}{2}\right)\\ & +h_x^{n+1}\left(i_f+1,j_f+\frac{R+1}{2},k_f+\frac{R+3}{2}\right)\\ & +h_x^{n+1}\left(i_f+1,j_f+\frac{R+3}{2},k_f+\frac{R+3}{2}\right)]/4\end{array}$$

(31)

Finally, we present computational steps of the hybrid subgridding LOD-FDTD method (as shown in Figure 3), for instance, (1) Employing the FDTD method to solve the EM field values (E and H) in the coarse-grid domain; (2) the interpolation is carried out on the interface of fine-grid domain by using Equations (19) and (20) so as to obtain the electric field value in the fine-grid domain; (3) the electric field ex in the fine-grid domain could be solved by applying Equations (17), (18), (26) and (27) (a similar way could be carried out to achieve ey and ez in the fine-grid domain); (4) the magnetic field value h in the fine-grid domain could be solved by means of Equation (11) and Equation (13); (5) finally the EM field values in the coarse-grid domain are corrected by adopting Equations (28)–(31).

3. Results and Discussion

3.1. Rectangular Waveguide

In order to demonstrate the accuracy of proposed ZT-LOD-FDTD method, a simulation model which has analytical solution is investigated [47]. As shown in Figure 4, a 3D rectangular waveguide model with a dispersive dielectric slab is considered, and the property of the dispersive dielectric slab can be described by the one-pole lossy Debye model. A plane wave with Gaussian pulse is performed and propagate along y-direction, and the parameters in the Gaussian function are set as Ts = 0.1 ns, t₀ = 5 Ts. The convolutional PML (CPML) is loaded at both ends of the propagation direction, and the periodic boundary condition (PBC) is set at the x and z directions. The grid size is chosen as 0.1 mm in x, y and z direction, and the thickness of the slab is set as 0.2 mm, and the parameters of the soil ε_∞ = 29.9, ∆ε₁ = 18, τ₁ = 43.6 fs, and σ = 0.54 S/m.

Figure 5 shows the transmission coefficients obtained by the ZT-LOD-FDTD method (CFLN = 1, 5, 10). For comparison, the result obtained by the analytical method is used as the benchmark. It is clearly shown that the results obtained by the ZT-LOD-FDTD method with different time steps are agreed with that of the analytical results, which testify the ZT-LOD-FDTD approach as an efficient technique to simulate Debye dispersive media.

With the purpose of validating the correctness and effectiveness of the proposed hybrid LOD-FDTD with subgridding scheme, three classical GPR scenarios are simulated by the proposed method. The simulation platform of the numerical examples in this work is: 32.0 GB RAM, Intel(R) Core(TM) i7-11700 CPU @ 3.80GHz. For truncating the physical domain, the CPML implementation based on Z-transform method is adopted. The underground detection targets comprise two EM structures: dielectric sphere and dielectric cylinder. In the free space, the transmitting and receiving antennas are respectively used as the excitation source point and the detection point, which are employed to transmit and receive detection signals. The excitation source in the time domain is defined as

$$pulse(t)=\{\begin{array}{l}-\frac{2\pi }{\tau }{\displaystyle \sum _{n=1}^3{a}_{n}n\sin (\frac{2\pi nt}{\tau })\begin{array}{cc}& 0<t<\tau \end{array}}\\ \begin{array}{ccc}\begin{array}{cc}& \end{array}& & \end{array}0\begin{array}{ccc}\begin{array}{cc}\begin{array}{cc}& \end{array}& \end{array}& & \mathrm{other}\end{array}\end{array},$$

(32)

where the centre frequency is 300 MHz, τ = 1.7/f_c, the coefficient a₁ = −0.493, a₂ = 0.144, a₃ = −0.01.

3.2. Dielectric Sphere

Considering the actual situation, the density, humidity and other properties of soil would vary due to its depth. As shown in Figure 6, the simulation domain in this example includes air, soil layer 1 and soil layer 2.

The detection target is a medium sphere whose upper and lower parts are located in these two soil media respectively. The dispersion property of soil is described by adopting the multi-level Debye model, and the parameters of the first-layer soil are ε_∞ = 3.0, ∆ε₁ = 0.10, ∆ε₂ = 0.25, τ₁ = 2.78 ns, τ₂ = 0.113 ns, and σ = 0.415 × 10⁻³ S/m. ones of second-layer soil are ε_∞ = 4.5, ∆ε₁ = 2.10, ∆ε₂ = 0.70, τ₁ = 4.08 ns, τ₂ = 0.261 ns, and σ = 1.11 × 10⁻³ S/m. The whole computational domain in size is 1.8 m by 1.8 m by 2.4 m in the x, y, and z directions. The spatial incremental size in the coarse-grid domain is set as Δx = Δy = Δz = 0.03 m, and the temporal incremental size in the coarse-grid domain is Δt = 0.0625 ns. The whole simulation domain occupies 60 × 60 × 80 Yee cells. The transmitting and receiving antennas are placed in the free space, and are 0.18 m far away from the ground, and the distance between transmitting and receiving antennas (Tx-Rx) is 0.45 m in the y direction. The subgridding domain in size is 0.3 m by 0.3 m by 0.3 m, the detection target is a dielectric sphere buried 0.3 m below the ground, with the relative dielectric constant of 30 and a radius of 0.12 m. To validate the correctness and effectiveness of the proposed LOD-FDTD method with subgridding scheme, the computational result of uniform FDTD method based on the fine grid (Δx_f = Δy_f = Δz_f = 0.01 m) is viewed as the reference solution.

It shows the results of the hybrid LOD-FDTD with subgridding scheme by using two ratios of grid generation with 1:3 and 1:4 respectively, and the result of the uniform fine-grid FDTD method are illustrated in Figure 7. It could be found from the agreement of waveforms that the subgridding method proposed in this work is suitable for odd and even discretizations.

Besides, for quantitatively analysing the numerical error, Equation (33) is adopted to calculate the relative error of the hybrid subgridding LOD-FDTD method. As reflected in Figure 8, the relative errors generated by odd or even times are relatively small, and meanwhile as the meshing ratio increases, the numerical error decreases, which further verifies the high versatility of the proposed method. As a result, it could be observed that using a larger time step in the fine-grid domain would not have a great influence on the numerical dispersion difference.

$$E_r=\frac{\left|E{z}_{\mathrm{fine}}-E{z}_{\mathrm{subgrid}}\right|}{\max \left|E{z}_{\mathrm{fine}}\right|}\times 100\%,$$

(33)

where Ez is obtained at the Rx grid location.

The received signal when the distance between the receiving antenna and the transmitting antenna changes from 0 m to 0.39 m are seen in Figure 9. Through the normalized time-domain waveform, it could be found that the amplitude of EM radiation scattered by the target and received by the receiving end would increase with the increase of Tx-Rx distance.

Ultimately,

Table 1. Memory Cost (mb) and Execution Times (s) for Fine-Grid FDTD and Hybrid Subgrid Method.

Method	Memory Cost (MB)	Execution Time (s)
Fine-grid FDTD	798.2	829.5
subgrid method (R = 3)	142.4	62.0
subgrid method (R = 4)	240.4	88.6

shows the memory occupation and CPU time of the two methods. It could be seen from Table 1 that the proposed hybrid subgridding LOD-FDTD method could substantially reduce the computational memory and CPU time.

3.3. Multi-Target Detection with Dielectric Cylinder and Dielectric Sphere

Considering the diversity of underground targets, the GPR scenario with multiple detection targets is simulated in the second case. In this example, both the sphere and the cylinder are placed in the same soil media.

As plotted in Figure 10, the entire domain presents the same characteristics with the Case A, the CPML we adopt is 10 layer thickness. The parameter of soil layer ε_∞ = 3.2, ∆ε₁ = 0.75, ∆ε₂ = 0.3, τ₁ = 2.71 ns, τ₂ = 0.108 ns, σ = 1.11 × 10⁻³ S/m. These two subgridding domains in size are 0.3 m by 0.3 m by 0.3 m with these two detected targets are buried 0.3 m away from the ground surface. The distance between the centers of these two targets is 0.15 m in the x direction. The dielectric constant of sphere is ${\epsilon }_{\mathrm{r}}=30$, its radius is 0.12 m. Meanwhile, the dielectric constant of cylinder is ${\epsilon }_{\mathrm{r}}=50$, its radius and length are respectively 0.09 m and 0.24 m. To validate the correctness and effectiveness of the proposed hybrid subgridding LOD-FDTD method, the result obtained by using uniform FDTD method with the fine-grid discretization is considered as the reference solution.

As reflected in Figure 11, the time-domain waveforms are shown at the receiving antenna calculated by the hybrid LOD-FDTD method with different meshing ratios and the traditional FDTD method with a uniform fine-grid discretization. Good consistency of the numerical results between them is observed. Furthermore, Figure 12 reflects the calculation error between the hybrid LOD-FDTD and the reference solution in the multi-target detecting scenario. As a result, the aforementioned result verifies that the recommended approach is correct. As illustrated in Figure 13, the transient changes of the electric field at different moments are extracted through the electric field snapshot. The transition from the coarse- to the fine-grid domain could be smoothed by the electric field, which reflects the stability and correctness of the coarse- and fine-grid field value exchange of the proposed hybrid subgridding LOD-FDTD method.

Table 2. Memory Cost (MB) and Execution Time (s) Consumption for Fine-Grid FDTD and Hybrid Subgridding Methods for Case 3.3.

Method	Memory Cost (MB)	Execution Time (s)
Fine-grid FDTD	1866.6	1182.4
subgrid method (R = 3)	188.8	75.4
subgrid method (R = 4)	390.0	129.9

shows the memory occupation and CPU time of the two methods. It could be seen from Table 2 that the proposed hybrid subgridding LOD-FDTD method could greatly reduce the computational memory and CPU time.

3.4. Dispersive Sphere

In the practical GPR detection, we often encounter the detection target whose dielectric coefficient is related to the frequency. Therefore, this example would simulate this kind of GPR scenario, replacing the dielectric sphere in Case A with the dispersive dielectric sphere in the same computational domain, and only considering monolayer soil. The dispersive dielectric sphere is simulated by the multi-level Debye model, its parameters are ε_∞ = 168.0151, ∆ε₁ = 34.2007, ∆ε₂ = 12.0813, τ₁ = 13.56 ms, τ₂ = 0.108 ms.

As seen in Figure 14, the time-domain waveforms of the electric field of the receiving antenna are respectively calculated by employing the hybrid subgridding LOD-FDTD method with the coarse-to-fine-grid ratio of 1:3 and 1:4 and the uniform fine-grid FDTD method. Taking the calculation result of uniform fine-grid FDTD as the reference solution, as shown in Figure 15, the errors of the two meshing methods are quantitatively analyzed. Besides, as we could observe that the time-domain waveforms are in the good consistency and the relative error is small, which verifies the correctness of the proposed method. Finally, as shown in

Table 3. Memory Cost (MB) and Execution Times (s) for Hybrid Subgridding and Fine-Grid FDTD Methods for Case 3.4.

Method	Memory Cost (MB)	Execution Time (s)
Fine-grid FDTD	610.5	798.2
subgrid method (R = 3)	59.2	142.2
subgrid method (R = 4)	86.7	240.2

, these three methods we proposed could effectively decrease the computational memory and CPU time.

4. Conclusions and Future Prospects

The explicit characteristic of regular FDTD method makes it difficult to trade off the accuracy and efficiency in the GPR scene simulation. In the GPR system, the multiscale characteristics of the EM structure and the amplitude of the EM field of the detection target change dramatically. With a view to ensuring the accuracy of numerical calculation, in general a smaller spatial incremental is adopted for grid generation, which greatly reduces the advantage of FDTD method in the GPR application. Hence the hybrid subgridding LOD-FDTD method is proposed in this work. The unconditionally stable LOD-FDTD method is adopted in the fine-grid domain, and the conventional FDTD method is used in other coarse-grid domains. This multiscale hybrid subgridding method could effectively simulate 3D GPR problems. Generally speaking, the main contributions of this work are as follows: (1) the LOD-FDTD method based on Z-transform technique is proposed to simulate multi-level Debye dispersive media, which is more direct in numerical realization than the convolutional method and auxiliary differential equation (ADE) method; (2) A hybrid subgridding technique combining the LOD-FDTD method and the traditional FDTD method is proposed. The unconditional stability of LOD-FDTD method enables the whole calculation domain to employ a uniform coarse-grid time step for numerical calculation. Therefore, only spatial interpolation processing is required, which reduces the application difficulty of subgridding method; (3) the subgridding method proposed in this work could divide the coarse grid into fine grids at any scale in the computational domain, which shows high versatility. Finally, three GPR scenarios are simulated to verify the correctness of the proposed method. Through numerical simulation, it can be found that the GPR system simulated by the proposed hybrid subgrid LOD-FDTD technology is about 5–10 times of the traditional FDTD method. If the GPR system needs more refined grid division in local areas of the simulation space, the computational advantages will be more obvious. The multi-scale subgrid technology based on unconditionally stable FDTD method is a research hotspot in the field of GPR forward simulation. The method proposed in this paper has a certain guiding role in the development of this field, and also provides research ideas for the development of other efficient GPR forward simulation methods.