Next Article in Journal
A Second Harmonic Wave Angle Sensor with a Collimated Beam of Femtosecond Laser
Previous Article in Journal
Numerical Simulation of the Enrichment of Chemotactic Bacteria in Oil-Water Two-Phase Transfer Fields of Heterogeneous Porous Media
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Modeling GPR Wave Propagation in Complex Underground Structures Using Conformal ADI-FDTD Algorithm

1
School of Water Conservancy Engineering, Zhengzhou University, Zhengzhou 450001, China
2
School of Mechanical Electronic & Information Engineering, China University of Mining & Technology, Beijing 100083, China
*
Author to whom correspondence should be addressed.
Submission received: 19 April 2022 / Revised: 17 May 2022 / Accepted: 19 May 2022 / Published: 21 May 2022
(This article belongs to the Topic Pipeline and Underground Space Technology)

Abstract

:
Ground Penetrating Radar (GPR) is a shallow geophysical method for detecting and locating subsurface targets. The GPR image echo characteristics of complex underground structures can be obtained by carrying out GPR forward modeling research. The traditional finite-difference time-domain (FDTD) method has low efficiency and accuracy. The alternating direction implicit FDTD (ADI-FDTD) algorithm surmounts the stability limitations of the traditional FDTD method, making it possible to select a larger time step for higher computational efficiency. For circular underground structures, a pseudowave produced by the ladder approximation method can be corrected using the surface conformal technique. This paper proposes a high-efficiency and high-accuracy GPR forward modeling method that combines the ADI-FDTD algorithm and surface conformal technology. The performance of the conformal ADI-FDTD algorithm is verified by a simple two-layer model. Based on the proposed algorithm, the GPR image features of three complex underground structure models are obtained. Finally, a field experiment is used to support the accuracy and usefulness of the conformal ADI-FDTD algorithm. The numerical simulation results and experimental results show that the conformal ADI-FDTD algorithm reduces the pseudodiffraction wave caused by the ladder approximation method and can significantly improve the computing efficiency for complex underground structure models.

1. Introduction

With rapid urban development, the scale of underground pipeline networks has gradually expanded [1,2]. Pipelines with different functions and purposes are arranged alternately, and the arrangement of the underground pipeline network is very complicated. As detailed distribution information for the entire underground pipeline network cannot be accurately obtained, accidents (e.g., the collapse of old pipelines or the cutting-off of existing pipelines) often occur [3,4,5,6]. In addition, road collapse accidents due to pipeline leaks occur frequently. The quick and accurate determination of the type, size, and location of underground pipelines and hidden diseases around underground pipelines has become necessary in urban road construction and maintenance projects.
GPR is widely regarded as one of the most powerful and useful geophysical methods [7]. GPR has the advantage of non-destructive detection and is fast, lightweight, and easy to operate with strong anti-interference and high resolution, etc. compared with other Non-Destructive Testing (NDT) methods such as seismic, sonic, temperature, resistivity, logging, and magnetic exploration methods [8,9]. Information about the location, size, and dielectric properties of complex underground structures can be obtained by inversion analysis of the GPR echo signals [10,11,12,13]. Through accurate and efficient GPR forward modeling of complex underground structure models, the propagation law of GPR electromagnetic waves in underground structures can be obtained [14,15,16]. This lays the theoretical foundation for inversion analysis and helps to improve the interpretation and processing accuracy of GPR-measured data. At present, the most commonly used GPR forward simulation techniques include the finite element method (FEM) [17,18], the ray tracing method (RTM) [19,20], the method of moment (MOM) [21], the finite difference time domain method (FDTD) [22,23,24], the pseudospectral time domain method (PSTD) [25], and the symplectic algorithm [26,27,28], among others. Although research into GPR forward simulation has achieved fruitful results, these algorithms still have some limitations in terms of computational accuracy and efficiency. For example, the FEM may appear as a “pseudosolving” phenomenon during calculation; the RTM cannot consider the dynamic characteristics of the electromagnetic wave of the ground penetrating radar, while the time step of the symplectic algorithm needs to satisfy the Courant–Friedrichs–Lewy (CFL) stability condition, and the computational efficiency is limited. Therefore, it is necessary to propose a GPR forward modeling method with high computational accuracy and fast computational efficiency.
The ADI-FDTD algorithm overcomes the CFL stability condition and can employ a larger time step to improve the computing simulation efficiency [29,30,31]. The efficiency and accuracy of GPR forward simulation are also related to the modeling method. Some methods are available to improve the accuracy of GPR simulation calculations for circular structures, such as sub-grid technology and conformal grid technology. The sub-grid technique has been widely used to improve the calculation accuracy, but it generates pseudowaves at the interface of the fine and coarse grids [32,33]. The derivations and calculations for conformal grid technology are simple, making it very suitable for circular underground pipes [34,35,36].
In this paper, the efficient and accurate GPR forward models are established to simulate GPR electromagnetic wave propagation in underground structures, employing the ADI-FDTD algorithm and surface conformal technology. We obtained the reflection profile images of multi-pipe models and complex underground pipe models with hidden diseases. By analyzing the profile images, the spatial propagation characteristics of radar waves can be more clearly understood, and the interpretation accuracy of the data can be improved. A field experiment proved the correctness and effectiveness of the proposed algorithm in actual detection. The numerical simulation results and field experiment indicate that the conformal ADI-FDTD method requires less computation time and can greatly reduce the pseudowave effect compared with the traditional FDTD.

2. Methodology

2.1. ADI-FDTD Algorithm

For transverse magnetic (TM) waves in two-dimensional lossy media, the Maxwell equations [37] are expressed as
μ H x t = E z y ,
μ H y t = E z x ,
ε E z t + σ E z = H y x H x y ,
where Hx and Hy are the magnetic field intensities in the x and y directions; Ez is the electric field intensity in the z direction; σ is the electrical conductivity; ε is the dielectric constant; μ is the magnetic permeability.
In the ADI-FDTD algorithm, we divide a time step into two sub-steps and discretize them in different ways. In the sub-step nn + 1/2, the x-direction derivative adopts the implicit difference format, and the y-direction derivative adopts the explicit difference format. In the sub-step of n + 1/2→n + 1, the x-direction derivative adopts the explicit difference format, and the y-direction derivative adopts the implicit difference format.
In the nn + 1/2 time step, the difference scheme of Equation (1) is as follows:
H x n + 1 2 ( i , j + 1 2 ) = H x n ( i , j + 1 2 ) Δ t 2 μ · 1 Δ y · [ E z n ( i , j + 1 ) E z n ( i , j ) ] .
Similarly, the difference discretization scheme of Equation (2) is
H y n + 1 2 ( i + 1 2 , j ) = H y n ( i + 1 2 , j ) + Δ t 2 μ · 1 Δ x · [ E z n + 1 2 ( i + 1 , j ) E z n + 1 2 ( i , j ) ] .
For Equation (3), the Ez component takes the mean value of time nΔt and time (n + 1/2)Δt, where the difference discretization scheme of Equation (3) is
E z n + 1 2 ( i , j ) = C A · E z n ( i , j ) + C B · { 1 Δ x · [ H y n + 1 2 ( i + 1 2 , j ) H y n + 1 2 ( i 1 2 , j ) ] 1 Δ y · [ H x n ( i , j + 1 2 ) H x n ( i , j 1 2 ) ] }
where
C A = 4 ε σ Δ t 4 ε + σ Δ t       C B = 2 Δ t 4 ε + σ Δ t .
Equations (4)–(6) are the time domain advancing formulas of the electromagnetic field from step nn + 1/2. Equation (4) is called the display format because the right side of the formula only involves the field value at time n. However, both sides of Equations (5) and (6) include the simultaneous field value, which is called the implicit format, which cannot be used to calculate the display time advance. To solve this problem, we substitute Equation (5) into Equation (6) and eliminate H y n + 1 / 2 , and we can derive
Δ t 2 μ Δ x 2 C B · E z n + 1 2 ( i 1 , j ) + [ 1 + Δ t μ Δ x 2 · C B ] · E z n + 1 2 ( i , j ) Δ t 2 μ Δ x 2 · C B · E z n + 1 2 ( i + 1 , j ) = C A · E z n ( i , j ) + C B Δ x · [ H y n ( i + 1 2 , j ) H y n ( i 1 2 , j ) ] C B Δ y · [ H x n ( i , j + 1 2 ) H x n ( i , j 1 2 ) ] ,
where
a i = Δ t 2 μ Δ x 2 · C B             b i = 1 + Δ t μ Δ x 2 · C B             c i = Δ t 2 μ Δ x 2 · C B d i = C A · E z n ( i , j ) + C B Δ x · [ H y n ( i + 1 2 , j ) H y n ( i 1 2 , j ) ] C B Δ y · [ H x n ( i , j + 1 2 ) H x n ( i , j 1 2 ) ] .
Then, Equation (8) can be expressed as
a i E z n + 1 2 ( i 1 , j ) + b i E z n + 1 2 ( i , j ) + c i E z n + 1 2 ( i + 1 , j ) = d i ,
and Equation (9) is written in matrix form as
A X = Y ,
where
X = [ E z n + 1 2 ( 1 , j ) E z n + 1 2 ( i , j ) E z n + 1 2 ( i m a x , j ) ] , Y = [ d 1 d i d i m a x ] , A = [ b 1 c 1 0 0 0 a 2 b 2 c 2 0 0 0 0 0 0 a i m a x 1 b i m a x 1 c i m a x 1 0 0 0 a i m a x b i m a x ] .
Here, the matrix A is a tridiagonal strip matrix, and the first-row elements (b1 and c1) and the last-row elements (aimax and bimax) are given by absorption boundary conditions.
In the n + 1/2→n + 1 time step, the difference discretization schemes of Equations (1)–(3) are
H x n + 1 ( i , j + 1 2 ) = H x n + 1 2 ( i , j + 1 2 ) Δ t 2 μ · 1 Δ y · [ E z n + 1 ( i , j + 1 ) E z n + 1 ( i , j ) ] ,
H y n + 1 ( i + 1 2 , j ) = H y n + 1 2 ( i + 1 2 , j ) + Δ t 2 μ · 1 Δ x · [ E z n + 1 2 ( i + 1 , j ) E z n + 1 2 ( i , j ) ] ,
E z n + 1 ( i , j ) = C A · E z n + 1 2 ( i , j ) + C B · { 1 Δ x · [ H y n + 1 2 ( i + 1 2 , j ) H y n + 1 2 ( i 1 2 , j ) ] 1 Δ y · [ H x n + 1 ( i , j + 1 2 ) H x n + 1 ( i , j 1 2 ) ] } .
The time domain advance of the electromagnetic field in sub-step n + 1/2→n + 1 can be calculated by Equations (13)–(15), as Equations (13) and (15) contain the simultaneous field value (i.e., implicit format) and cannot be used to calculate the display time advance. From Equations (13) and (15), we can derive
Δ t 2 μ Δ y 2 · C B · E z n + 1 ( i , j 1 ) + [ 1 + Δ t μ Δ y 2 · C B ] · E z n + 1 ( i , j ) Δ t 2 μ Δ y 2 · C B · E z n + 1 ( i , j + 1 ) = C A · E z n + 1 2 ( i , j ) C B Δ y · [ H x n + 1 2 ( i , j + 1 2 ) H x n + 1 2 ( i , j 1 2 ) ] + C B Δ x · [ H y n + 1 2 ( i + 1 2 , j ) H y n + 1 2 ( i 1 2 , j ) ] ,
where
a j = Δ t 2 μ Δ y 2 · C B           b j = 1 + Δ t μ Δ y 2 · C B           c j = Δ t 2 μ Δ y 2 · C B   d j = C A · E z n + 1 2 ( i , j ) C B Δ y · [ H x n + 1 2 ( i , j + 1 2 ) H x n + 1 2 ( i , j 1 2 ) ] + C B Δ x · [ H y n + 1 2 ( i + 1 2 , j ) H y n + 1 2 ( i 1 2 , j ) ] .
Then, Equation (16) can be rewritten as
a j E z n + 1 ( i , j 1 ) + b j E z n + 1 ( i , j ) + c j E z n + 1 ( i , j + 1 ) = d j ,
and Equation (17) can be written in the following matrix form:
A X = Y ,
where
X = [ E z n + 1 ( i , 1 ) E z n + 1 ( i , j ) E z n + 1 ( i , j m a x ) ] , Y = [ d 1 d j d j m a x ] , A = [ b 1 c 1 0 0 0 a 2 b 2 c 2 0 0 0 0 0 0 a j m a x 1 b j m a x 1 c j m a x 1 0 0 0 a j m a x b j m a x ] .
Here, matrix A is a tridiagonal strip matrix, and the first-row elements (b1 and c1) and the last-row elements (ajmax and bjmax) are given by absorption boundary conditions.
In conclusion, the steps for calculating the 2D TM wave by ADI-FDTD are as follows:
  • In the nn + 1/2 time step:
(1) Calculate Hxn+1/2 by Equation (4);
(2) Calculate Ezn+1/2 by Equation (10);
(3) Calculate Hyn+1/2 by Equation (5).
  • In the n + 1/2→n + 1 time step:
(1) Calculate Hyn+1 by Equation (14);
(2) Calculate Ezn+1 by Equation (18);
(3) Calculate Hxn+1 by Equation (13).

2.2. Surface Conformal Technology

The surface conformal technique was used to simulate the mesh generation of circular underground structures. Figure 1 shows schematic diagrams of the grid division of the circular structure model constructed by different methods. Figure 1a shows the actual subdivision grid of the circular structure, where the white squares are ordinary and the red squares are the actual subdivisions. Figure 1b shows the conventional ladder approximation subdivision grid for a circular structure, where the green squares are the subdivisions obtained by the step approximation method. The subdivision grid of a circular structure based on surface conformal technology is shown in Figure 1c, where the orange squares are conformal areas and the blue squares are non-conformal areas.
A conformal grid point in Figure 1c is taken as an example in order to demonstrate the selection of equivalent medium parameters for the conformal grid region in 2D TM waves. As shown in Figure 2, where F is the sampling point of the electric field Ez; A, B and C, D are sampling points of the magnetic field Hx and Hy, respectively; Δy and Δx are the height and width of the grid, respectively; Sxy1 and Sxy2 are the areas of media 1 and 2 in the conformal grid, respectively; Lx1 and Lx2 are the lengths occupied by media 1 and 2 on the edge of magnetic field node B, respectively; and Ly1 and Ly2 are the lengths occupied by media 1 and 2 on the edge of magnetic field node D, respectively.
Considering the two-dimensional TM wave, for the conformal grid of an ideal conductor (metal conductor), the recursive formula of a magnetic field node remains unchanged and is still calculated according to conventional ADI-FDTD. For the recursive formula of the electric field node, we make the following changes:
In the nn + 1/2 time step, Equation (6) is changed to
E z n + 1 2 ( i , j ) = C A · E z n ( i , j ) + C B · { l y S · [ H y n + 1 2 ( i + 1 2 , j ) H y n + 1 2 ( i 1 2 , j ) ] l x S · [ H x n ( i , j + 1 2 ) H x n ( i , j 1 2 ) ] } .
where lx and ly represent the length of the corresponding edge of the electric field node outside the conductor, and S represents the area of the cell outside the conductor.
We substitute H y n + 1 / 2 in Equation (5) into Equation (21) and obtain the update equation about E z n + 1 / 2
a i E z n + 1 2 ( i 1 , j ) + b i E z n + 1 2 ( i , j ) + c i E z n + 1 2 ( i + 1 , j ) = d i ,
where
a i = Δ t 2 μ Δ x · l y S · C B           b i = 1 + Δ t μ Δ x · l y S · C B           c i = Δ t 2 μ Δ x · l y S · C B d i = C A · E z n ( i , j ) + C B · l y S · [ H y n ( i + 1 2 , j ) H y n ( i 1 2 , j ) ] C B · l x S · [ H x n ( i , j + 1 2 ) H x n ( i , j 1 2 ) ] .
In the n + 1/2→n + 1 time step, the difference discretization schemes of Equation (15) are
E z n + 1 ( i , j ) = C A · E z n + 1 2 ( i , j ) + C B · { l y S · [ H y n + 1 2 ( i + 1 2 , j ) H y n + 1 2 ( i 1 2 , j ) ] l x S · [ H x n + 1 ( i , j + 1 2 ) H x n + 1 ( i , j 1 2 ) ] } ,
Similarly, we can obtain the iterative format of E z n + 1
a j E z n + 1 ( i , j 1 ) + b j E z n + 1 ( i , j ) + c j E z n + 1 ( i , j + 1 ) = d j ,
where
a j = Δ t 2 μ Δ y · l x S · C B             b j = 1 + Δ t μ Δ y · l x S · C B             c j = Δ t 2 μ Δ y · l x S · C B   d j = C A · E z n + 1 2 ( i , j ) C B · l x S · [ H x n + 1 2 ( i , j + 1 2 ) H x n + 1 2 ( i , j 1 2 ) ] + C B · l y S · [ H y n + 1 2 ( i + 1 2 , j ) H y n + 1 2 ( i 1 2 , j ) ] .
For media (non-metallic conductor) surface conformal meshes, the permittivity, permeability, and conductivity of dielectrics 1 and 2 are assumed to be ε1, μ1, and σ1 and ε2, μ2, and σ2, respectively. A weighted average of the electromagnetic coefficients, according to the sizes of the conformal grids occupied by media 1 and 2, respectively, provides the following equivalent medium parameters for the dielectric constant, permeability, and conductivity:
ε e f f = L x 1 ε 1 + L x 2 ε 2 δ ,   σ e f f = L x 1 σ 1 + L x 2 σ 2 δ ,   μ e f f = S x y 1 μ 1 + S x y 2 μ 2 δ 2 .
In the nn + 1/2 time step, the difference discretization schemes of the electric field node and the magnetic field node are written as
H x n + 1 2 ( i , j + 1 2 ) = H x n ( i , j + 1 2 ) Δ t 2 μ e f f · 1 Δ y · [ E z n ( i , j + 1 ) E z n ( i , j ) ] ,
H y n + 1 2 ( i + 1 2 , j ) = H y n ( i + 1 2 , j ) + Δ t 2 μ e f f · 1 Δ x · [ E z n + 1 2 ( i + 1 , j ) E z n + 1 2 ( i , j ) ] ,
a i E z n + 1 2 ( i 1 , j ) + b i E z n + 1 2 ( i , j ) + c i E z n + 1 2 ( i + 1 , j ) = d i ,
where
C A = 4 ε e f f σ e f f Δ t 4 ε e f f + σ e f f Δ t             C B = 2 Δ t 4 ε e f f + σ e f f Δ t a i = Δ t 2 μ e f f Δ x 2 · C B             b i = 1 + Δ t μ e f f Δ x 2 · C B             c i = Δ t 2 μ e f f Δ x 2 · C B d i = C A · E z n ( i , j ) + C B Δ x · [ H y n ( i + 1 2 , j ) H y n ( i 1 2 , j ) ] C B Δ y · [ H x n ( i , j + 1 2 ) H x n ( i , j 1 2 ) ] .
In the n + 1/2→n + 1 time step, the difference discretization schemes of the electric field node and the magnetic field node are
H x n + 1 ( i , j + 1 2 ) = H x n + 1 2 ( i , j + 1 2 ) Δ t 2 μ e f f · 1 Δ y · [ E z n + 1 ( i , j + 1 ) E z n + 1 ( i , j ) ] ,
H y n + 1 ( i + 1 2 , j ) = H y n + 1 2 ( i + 1 2 , j ) + Δ t 2 μ e f f · 1 Δ x · [ E z n + 1 2 ( i + 1 , j ) E z n + 1 2 ( i , j ) ] ,
a j E z n + 1 ( i , j 1 ) + b j E z n + 1 ( i , j ) + c j E z n + 1 ( i , j + 1 ) = d j ,
where
C A = 4 ε e f f σ e f f Δ t 4 ε e f f + σ e f f Δ t             C B = 2 Δ t 4 ε e f f + σ e f f Δ t a j = Δ t 2 μ e f f Δ y 2 · C B             b j = 1 + Δ t μ e f f Δ y 2 · C B             c j = Δ t 2 μ e f f Δ y 2 · C B   d j = C A · E z n + 1 2 ( i , j ) C B Δ y · [ H x n + 1 2 ( i , j + 1 2 ) H x n + 1 2 ( i , j 1 2 ) ] + C B Δ x · [ H y n + 1 2 ( i + 1 2 , j ) H y n + 1 2 ( i 1 2 , j ) ] .

2.3. UPML Absorbing Boundary Condition

GPR electromagnetic wave propagation in underground structures is an open domain problem. Thus, it is necessary to set reasonable absorbing boundary conditions at the truncated boundary of the computational region. In this paper, the uniaxial anisotropic absorption layer (UPML) with easy programming and simple iterative formula is used as the absorption boundary [38]. The Maxwell curl equation in the UPML medium can be written as:
× H = j ω ε S ¯ · E × E = j ω μ S ¯ · H
where S ¯ has the characteristics of a uniaxial anisotropic medium, expressed as:
S ¯ = [ s y s z / s x 0 0 0 s z s x / s y 0 0 0 s x s y / s z ] .
Here, si = κi + σi/jωε0 (i = x, y, z) is the diagonal tensor of a uniaxially anisotropic medium (for a 2D GPR wave, sz is 1), σi is the attenuation factor of the UPML region, and κi is used to absorb the modulated wave that reaches the UPML layer.
For the 2D TM wave, Ex = 0, Ey = 0, and Hz = 0, and the ADI-FDTD formula under the UPML boundary conditions is shown below.
When Dz = εsyEz, the time advance formula of Hx, HyDz is as follows:
H y x H x y = κ x D z t + σ x ε 0 D z ,
and the time advance formula of DzEz is:
D z t = ε 1 κ y E z t + ε 1 σ y ε 0 E z ,
When Bx = μ1Hx/sx and By = μ1Hy/sy, the time advance formulae for EzBx, By are as follows:
E z y = κ y B x t σ y ε 0 B x ,
E z x = κ x B y t + σ x ε 0 B y ,
and the time advance formulae for Bx, ByHx and Hy are as follows:
κ x B x t + σ x ε 0 B x = μ 1 H x t ,
κ y B y t + σ y ε 0 B y = μ 1 H y t ,
setting
a x = 2 κ x / Δ t + σ x / 2 ε 0             a y = 2 κ y / Δ t + σ y / 2 ε 0             e r = 2 / ( ε 1 Δ t ) b x = 2 κ x / Δ t σ x / 2 ε 0             b y = 2 κ y / Δ t σ y / 2 ε 0             h r = Δ t / ( 2 μ 1 )
In sub-step nn + 1/2, discretize equation (38) as
D z n + 1 2 ( i , j ) = b x a x · D z n ( i , j ) + 1 a x · { 1 Δ x · [ H y n + 1 2 ( i + 1 2 , j ) H y n + 1 2 ( i 1 2 , j ) ] 1 Δ y · [ H x n ( i , j + 1 2 ) H x n ( i , j 1 2 ) ] } .
The discretization scheme of Equation (39) is
E z n + 1 2 ( i , j ) = b y a y · E z n ( i , j ) + e r a y · [ D z n + 1 2 ( i , j ) D z n ( i , j ) ] ,
where
D z n + 1 2 ( i , j ) = D z n ( i , j ) + a y e r · E z n + 1 2 ( i , j ) b y e r · E z n ( i , j ) .
The discretization schemes of Equations (40)–(43), respectively, are
B x n + 1 2 ( i , j + 1 2 ) = b y a y · B x n ( i , j + 1 2 ) 1 a y · 1 Δ y · [ E z n ( i , j + 1 ) E z n ( i , j ) ] ,
B y n + 1 2 ( i + 1 2 , j ) = b x a x · B y n ( i + 1 2 , j ) + 1 a x · 1 Δ x · [ E z n + 1 2 ( i + 1 , j ) E z n + 1 2 ( i , j ) ] ,
H x n + 1 2 ( i , j + 1 2 ) = H x n ( i , j + 1 2 ) + h r · [ a x · B x n + 1 2 ( i , j + 1 2 ) b x · B x n ( i , j + 1 2 ) ] ,
H y n + 1 2 ( i + 1 2 , j ) = H y n ( i + 1 2 , j ) + h r · [ a y · B y n + 1 2 ( i + 1 2 , j ) b y · B y n ( i + 1 2 , j ) ] ,
In Equations (45)–(51), except for Equation (48), both sides contain the field component at time (n + 1/2). By eliminating the simultaneous components from Equations (45)–(51), except for Ezn+1/2, we obtain
a i E z n + 1 2 ( i 1 , j ) + b i E z n + 1 2 ( i , j ) + c i E z n + 1 2 ( i + 1 , j ) = d i ,
where
a i = h r ( Δ x ) 2 · a y a x           b i = a x · a y e r + h r ( Δ x ) 2 · 2 · a y a x           c i = h r ( Δ x ) 2 · [ a y ( i + 1 2 , j ) / a x ( i + 1 2 , j ) ] d i = a x · b y e r E z n ( i , j ) + { 1 Δ x · [ H y n ( i + 1 2 , j ) H y n ( i 1 2 , j ) ] 1 Δ y · [ H x n ( i , j + 1 2 ) H x n ( i , j 1 2 ) ] } + ( b x a x ) · D z n ( i , j ) + h r Δ x · [ a y · b x a x b y ] · B y n ( i + 1 2 , j ) [ a y · b x a x b y ] · B y n ( i 1 2 , j ) .
In the sub-step of n + 1/2→n + 1, the difference scheme of Equation (38) is the same as above, which can be expressed as
D z n + 1 ( i , j ) = b x a x · D z n + 1 2 ( i , j ) + 1 a x · { 1 Δ x · [ H y n + 1 2 ( i + 1 2 , j ) H y n + 1 2 ( i 1 2 , j ) ] 1 Δ y · [ H x n + 1 ( i , j + 1 2 ) H x n + 1 ( i , j 1 2 ) ] } .
Equation (39) can be discretized as
E z n + 1 ( i , j ) = b y a y · E z n + 1 2 ( i , j ) + e r a y · [ D z n + 1 ( i , j ) D z n + 1 2 ( i , j ) ] ,
where
D z n + 1 ( i , j ) = D z n + 1 2 ( i , j ) + a y e r · E z n + 1 ( i , j ) b y e r · E z n + 1 2 ( i , j ) .
Equations (40)–(43), respectively, can be discretized as
B x n + 1 ( i , j + 1 2 ) = b y a y · B x n + 1 2 ( i , j + 1 2 ) 1 a y · 1 Δ y · [ E z n + 1 ( i , j + 1 ) E z n + 1 ( i , j ) ] ,
B y n + 1 ( i + 1 2 , j ) = b x a x · B y n + 1 2 ( i + 1 2 , j ) + 1 a x · 1 Δ x · [ E z n + 1 2 ( i + 1 , j ) E z n + 1 2 ( i , j ) ] ,
H x n + 1 ( i , j + 1 2 ) = H x n + 1 2 ( i , j + 1 2 ) + h r · [ a x · B x n + 1 ( i , j + 1 2 ) b x · B x n + 1 2 ( i , j + 1 2 ) ] ,
H y n + 1 ( i + 1 2 , j ) = H y n + 1 2 ( i + 1 2 , j ) + h r · [ a y · B y n + 1 ( i + 1 2 , j ) b y · B y n + 1 2 ( i + 1 2 , j ) ] .
Eliminating the simultaneous components of Equations (54)–(60), the time domain advance calculation is obtained as follows:
a j E z n + 1 ( i , j 1 ) + b j E z n + 1 ( i , j ) + c j E z n + 1 ( i , j + 1 ) = d j ,
where
a i = h r ( Δ y ) 2 · a x a y           b i = a x · a y e r + h r ( Δ y ) 2 · 2 a x a y           c i = h r ( Δ y ) 2 · a x a y d i = a x · b y e r E z n + 1 2 ( i , j ) + { 1 Δ x · [ H y n + 1 2 ( i + 1 2 , j ) H y n + 1 2 ( i 1 2 , j ) ] 1 Δ y · [ H x n + 1 2 ( i , j + 1 2 ) H x n + 1 2 ( i , j 1 2 ) ] } + ( b x a x ) · D z n + 1 2 ( i , j ) h r Δ y · [ ( a x · b y a y b x ) · B x n + 1 2 ( i , j + 1 2 ) ( a x · b y a y b x ) · B x n + 1 2 ( i , j 1 2 ) ] .
The steps of ADI-FDTD for the UPML boundary of 2D TM waves are as follows:
  • In the time step of nn + 1/2:
(1) Calculate Ezn+1/2 by Equation (52);
(2) Calculate Dzn+1/2 by Equation (47);
(3) Calculate Bxn+1/2 and Byn+1/2 by Equations (48) and (49);
(4) Calculate Hxn+1/2 and Hyn+1/2 by Equations (50) and (51).
  • In the time step of n+1/2→n+1:
(1) Calculate Ezn+1 by Equation (61);
(2) Calculate Dzn+1 by Equation (56),
(3) Calculate Bxn+1 and Byn+1 by Equations (57) and (58);
(4) Calculate Hxn+1 and Hyn+1 by Equations (59) and (60).
We used a simple air model to verify the absorption effect of the UPML boundary conditions. Visual Studio 2010 was used as a development tool, and the CPU was an Intel Core i5-4200H with an NVIDIA GeForce GTX 950m. The computing environment above supported all the computing processes in this paper. The simulation area was a 2.0 m × 2.0 m rectangular area; the spatial step was 0.005 m, and the time step was Δt = 0.01 ns. A Ricker wavelet was added at the positive center of the model’s area (see Figure 3). Figure 4 shows the snapshots of the Ez field at different moments. It can be seen that the UPML absorption boundary condition had no obvious reflections, and the absorption effect was good.

3. Numerical Simulations

3.1. Two-Layer Medium Model with Circular Cavity

As shown in Figure 5, the model area was set to 2.0 m × 2.0 m. The upper layer of the model was a 0.2 m air layer, and the lower layer was a 1.8 m clay layer. The main frequency of the excitation source was a 1 GHz Ricker wavelet, which was located at a 0.2 m depth in the model, and the transmitter and receiver were 0.1 m apart. The relative dielectric constant of the clay layer was 12, and the conductivity was 2 mS/m. A circular cavity with a diameter of 0.1 m was set in the clay layer. The relative dielectric constant and conductivity were 30 and 0 mS/m, respectively. The relative magnetic conductivity of all the material was 1. We set the spatial step size to 0.005 m, the time step size to 0.01 ns, and the thickness of the UPML to 0.01 m.
The time step of the conventional FDTD algorithm needs to satisfy the CFL stability condition. However, the ADI-FDTD algorithm is unconditionally stable and can save computation time by choosing a larger time step. Figure 6 shows a comparison of the Ez field distributions of the conformal ADI-FDTD with different CFLN, the non-conformal FDTD, and the non-conformal ADI-FDTD. The CFLN is the ratio of the time step of the ADI-FDTD algorithm to that of the FDTD algorithm. It can be seen that the results of the non-conformal FDTD, the non-conformal ADI FDTD algorithm, and the conformal ADI FDTD algorithm have good consistency in different time steps. Table 1 provides the calculation time of the different algorithms when simulating the circular cavity model, respectively. When CFLN is 5, compared with the non-conformal FDTD algorithm, the time step of the conformal ADI-FDTD algorithm is increased to 5 times, and the number of iterations is reduced to 1/5. At this point, the computation time of the conformal ADI-FDTD algorithm is 961 s and that of the non-conformal FDTD algorithm is 1849 s, saving about 48% of the time. Figure 7 shows the GPR profiles obtained by the non-conformal FDTD method, non-conformal ADI-FDTD method, and the conformal ADI-FDTD method. We indicate multiple reflections within the cavity with red boxes. As can be seen from Figure 7, the images of the non-conformal FDTD method and the non-conformal ADI-FDTD method are consistent. The conformal ADI-FDTD method has fewer multiple reflections than the non-conformal ADI-FDTD method. At the same time, the conformal ADI-FDTD method is closer to the conformal FDTD method than the non-conformal ADI-FDTD method. This indicates that the conformal ADI-FDTD method can effectively reduce the false reflection waves caused by conventional ladder approximation when simulating an underground circular cavity. Table 1 and Figure 7 show that the conformal ADI-FDTD algorithm can greatly improve the computational efficiency and can also effectively improve the computational accuracy.

3.2. Single-Pipe Model with Cavity Disease

This model simulated two types of the cavity disease around the pipeline, which further verifies the accuracy of the conformal algorithm. The model had two layers with an air layer on top and a clay layer on the bottom. As shown in Figure 8, there were some irregular cavities above the concrete pipe with an outer diameter of 0.6 m and an inner diameter of 0.48 m. Figure 9 shows a concrete pipe with an outer diameter of 1.2 m and an inner diameter of 0.96 m with some irregular cavities underneath. The spatial steps and time steps were 0.005 m and 0.01 ns, respectively, and the simulation iterations were 5000. Table 2 shows the conductivity and relative permittivity of different materials in the model. The excitation source of the model was the same as above.
By analyzing Figure 10, the non-conformal model profile only shows one distinct reflected wave at the bottom of the pipe. The profile of the conformal model shows that there are two reflected waves on the upper and lower sides of the pipe, and the reflected waves of the irregular cavity disease above the pipe are incredibly apparent. As shown in Figure 11, both the non-conformal and conformal model profiles clearly show the reflected waves at the top and bottom of the pipe, but the non-conformal model profile only shows the reflected waves from the cavity damage on the left side, while the conformal model profile completely shows the reflected waves distributed on both sides of the cavity damage. As can be seen in Figure 10 and Figure 11, the circular pipe grid points are more accurately processed using the conformal grid technique, and the conformal ADI FDTD algorithm is more accurate.

3.3. Complex Multi-Pipe Model

This model used the conformal ADI-FDTD algorithm to simulate the complex multi-tube model. Figure 12 and Figure 13 show the schematic diagrams of complex multi-pipe models. Figure 12 shows an underground concrete multi-pipe model, with pipe diameters of 0.03 m, 0.06 m, and 0.1 m. Figure 13 shows an underground multi-pipe model with pipe diameters of 0.4 m. The pipe materials were concrete, metal, and PVC plastic. The space step and time step were set as 0.005 m and 0.01 ns, and simulation iterations were 5000. The model had two layers: the upper layer was an air layer, and the lower layer was a clay layer. Table 2 provides the relative permittivity and conductivity. The excitation source and GPR system were the same as for the single-pipe model.
Figure 14 and Figure 15 show the GPR B-scan images obtained by the underground multi-pipe model using the conformal ADI-FDTD algorithm. Figure 14, G1, G2, and G3 indicate the intersection points of the hyperbolic asymptotes at the tops of concrete pipes with diameters of 0.3 m, 0.6 m, and 1.0 m, respectively. Reflected waves are evident at the top and bottom of the pipeline; the more reflective the diameter, the more pronounced the reflected waves are. There are multiple reflected and stray waves at the bottom of the profile image. In Figure 15, the PVC plastic pipe on the right side of the profile has weak reflected waves, and the concrete pipeline located on the left side has distinct hyperbolic reflected waves. The reflected waves are evident at the top of the middle metal tube, while not at the bottom. In addition, some interfering stray waves and multiple reflected waves can be observed at the bottom of the image.
As shown in Figure 14 and Figure 15, for an underground multi-pipe model with the same depth and materials but different diameters, an increase in diameter causes the intersection position of the hyperbolic asymptote at the top of the pipe to become higher, and the reflection waves at the bottom and top of the pipe will be more obvious. As the thickness of the pipe wall is related to the pipe’s diameter, two obvious hyperbolas with larger diameters can be seen at the bottom and top of the concrete pipe. For underground multi-pipe models with the same diameter and composed of different materials, the closer the dielectric properties of the pipe material are to those of the surrounding clay, the weaker its response to the electromagnetic waves and the less pronounced the diffraction hyperbolic features are at the top and bottom. Due to the strong response of electromagnetic waves to metals, the diffraction hyperbolic characteristics at the top of the metal pipe were very obvious, while there were no hyperbolic characteristics at the bottom of the pipe and fewer multiple reflections occurred inside.

3.4. Complex Underground Structures

We used the conformal ADI-FDTD algorithm to simulate complex underground structure models. Figure 16 shows schematic diagrams of rectangular and circular void models. There was a metal pipe located 0.5 m below the surface with an outside diameter of 0.3 m and an inside diameter of 0.2 m. There was a rectangular cavity of a size of 0.1 m on the left side of the metal pipe. On the other side, there was a circular cavity of a diameter of 0.1 m. The space step and time step were set as 0.005 m and 0.01 ns, and simulation iterations were 5000. The model had two layers: the upper layer was an air layer, and the lower layer was a clay layer. Table 2 provides the relative permittivity and conductivity. The excitation source and GPR system were the same as for the previous model.
Figure 17 shows the GPR B-scan image obtained by the underground cavity model using the conformal ADI-FDTD algorithm. By analyzing Figure 17, the metal pipe produces strong reflected and diffracted waves. It is due to the strong response of the electromagnetic waves to metals. The upper interface of the rectangular cavity is still a flat interface in the radar forward image. Meanwhile, the lower interface is similar to the upper interface but retains weak reflected wave energy. The circular cavity produces a clear reflected wave, and the diffracted wave is weaker due to its smaller size. Moreover, some mutual interference clutters and multiple reflection waves can be observed at the bottom of the image. By analyzing the GPR scans, the characteristics of the model can be inferred more comprehensively and accurately, helping to interpret the GPR information.

4. Field Experiment

The correctness and effectiveness of the conformal ADI-FDTD algorithm were verified by a GPR detection experiment of rectangular voids. As shown in Figure 18, a carton with a length of 0.5 m, a width of 0.3 m, and a height of 0.16 m was buried in the soil, and the top surface of the carton was 0.16m from the ground. In this experiment, GSSI SIR-4000 road radar was used for detection, and the center frequency was 400 MHz.
The experimental model was simplified to a two-dimensional case, as shown in Figure 19. It was assumed that the relative permittivity of the soil layer was 6.8 and the conductivity was 3 mS/m. The time step was set to 0.001 ns; the space step was set to 0.005 m, and the simulation iterations were taken as 2000. Figure 20 was a single-channel wave comparison diagram of the numerical simulation results of the conformal ADI-FDTD algorithm and the measured results.
It can be seen from Figure 20 that the simulated waveform obtained by the algorithm in this paper is in good agreement with the measured waveform in terms of amplitude and delay, with differences only in details. This is mainly due to the fact that we considered the soil layer as a homogeneous and isotropic medium, which led to the difference between the assumed dielectric constant of the material and the real dielectric constant of the material during the simulation. This example verifies the applicability and effectiveness of the proposed algorithm for practical engineering inspection.

5. Conclusions

In this study, an efficient and accurate GPR forward model based on the ADI-FDTD method and surface conformal technology was established. The GPR image characteristics of the different shaped cavity diseases and multiple underground pipeline structures with different materials and diameters were obtained. The model provides a basis for further processing and interpretation of the measured data of GPR electromagnetic waves. The numerical simulation results demonstrate that the conformal ADI-FDTD method can greatly reduce the errors caused by the ladder approximation method, which is traditionally employed to divide a circular pipe grid. Moreover, when the time step is increased to 5 times, the conformal ADI-FDTD algorithm can save 48% of the computation time compared to the traditional FDTD method. The correctness and effectiveness of the conformal ADI-FDTD algorithm was verified by a field experiment. In future research, we intend to conduct GPR inversion analyses of the underground pipe structures. Through a 2D GPR image inversion analysis of the underground pipe structures, we will determine the positions and sizes of cavity diseases in underground pipes.

Author Contributions

Methodology, Y.L. and F.W.; software, Y.L. and N.W.; validation, J.L.; writing—original draft preparation, Y.L.; writing—review and editing, J.L. and N.W.; supervision, F.W. and C.L.; investigation, C.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by [the National Key Research and Development Program of China] grant number [2017YFC1501204], [the National Natural Science Foundation of China] grant number [11771407], [the Program for Science and Technology Innovation Talents in Universities of Henan Province] grant number [19HASTIT043], and [the Outstanding Young Talent Research Fund of Zhengzhou University] grant number [1621323001]. And The APC was funded by [the Program for Innovative Research Team at the University of Henan Province] grant number [18IRTSTHN007].

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

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. 11771407), the Program for Science and Technology Innovation Talents in Universities of Henan Province (No. 19HASTIT043), the Outstanding Young Talent Research Fund of Zhengzhou University (1621323001), and the Program for Innovative Research Team at the University of Henan Province (18IRTSTHN007). The authors would like to thank for these financial supports.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Zhang, Z.; Xu, H.; Li, X. Application education about urban underground pipeline management of point-source information system. In Proceedings of the 2010 International Conference on Computer and Communication Technologies in Agriculture Engineering, Chengdu, China, 12–13 June 2010; pp. 238–241. [Google Scholar]
  2. Xu, H.; Zhang, Q.; Li, X. Based on 3D technology application research about urban underground pipeline management information system. Appl. Mech. Mater. 2010, 29–32, 2227–2232. [Google Scholar] [CrossRef]
  3. Pang, L.; Lan, G.; Tao, Y. Design and Implementation of 3D Urban Underground Pipe Network System. In Proceedings of the IOP Conference Series: Earth and Environmental Science, Nanjing, China, 13–14 April 2019; Volume 283, p. 012013. [Google Scholar]
  4. Zhu, W.; Sun, P.; Ma, Y.; Song, R.; He, S.; Liu, K.H. Method of Underground Pipeline Monitoring Point Based on Cavity Detection. In Proceedings of the Information Computing and Applications, Chengde, China, 14–16 September 2012; Volume 7473, pp. 75–81. [Google Scholar]
  5. Zhao, S.; Al-Qadi, I. Pavement drainage pipe condition assessment by GPR image reconstruction using FDTD modeling. Constr. Build. Mater. 2017, 154, 1283–1293. [Google Scholar] [CrossRef]
  6. Zhang, P.; Muhammat, N. A New Method on Probing and Interpreting Underground Pipeline’s Diameter by GPR. In Tunneling and Underground Construction; American Society of Civil Engineers: Reston, VA, USA, 2014; pp. 855–861. [Google Scholar]
  7. Ma, F. Recent Advances in the GPR Detection of Grouting Defects behind Shield Tunnel Segments. Remote Sens. 2021, 13, 4596. [Google Scholar]
  8. Bai, X.; An, W.; Wang, B.; Jiang, J.; Zhang, Y.; Zhang, J. Automatic Identification of Underground Pipeline Based on Ground Penetrating Radar. In Proceedings of the International Conference on Wireless and Satellite Systems, Harbin, China, 12–13 January 2019. [Google Scholar]
  9. Yang, H.W.; Yang, Z.K.; Pei, Y.K. Ground-penetrating radar for soil and underground pipelines using the forward modeling simulation method. Optik 2014, 125, 7075–7079. [Google Scholar] [CrossRef]
  10. Keskinen, J.; Looms, M.C.; Klotzsche, A.; Nielsen, L. Practical data acquisition strategy for time-lapse experiments using crosshole GPR and full-waveform inversion. J. Appl. Geophys. 2021, 191, 104362. [Google Scholar] [CrossRef]
  11. Kaplanvural, S.; Peken, E.; Zkap, K. 1D waveform inversion of GPR trace by particle swarm optimization. J. Appl. Geophys. 2020, 181, 104157. [Google Scholar] [CrossRef]
  12. Zhang, F.; Liu, B.; Wang, J.; Li, Y.; Nie, L.; Wang, Z.; Zhang, C. Two-dimensional Time-domain Full Waveform Inversion of On-ground Common-offset GPR Data Based on Integral Preprocessing. J. Environ. Eng. Geophys. 2020, 25, 369–380. [Google Scholar] [CrossRef]
  13. Motevalli, Z.; Zakeri, B. Time-Domain Spectral Inversion Method for Characterization of Subsurface Layers in Ground-Penetrating-Radar (GPR) Applications. Appl. Comput. Electromagn. Soc. J. 2019, 34, 93–99. [Google Scholar]
  14. Yang, F.; Qiao, X.; Zhang, Y.; Xu, X. Prediction method of underground pipeline based on hyperbolic asymptote of GPR image. In Proceedings of the 15th International Conference on Ground Penetrating Radar, Brussels, Belgium, 30 June–4 July 2014. [Google Scholar]
  15. Caselle, C.; Bonetto, S.; Comina, C.; Stocco, S. GPR surveys for the prevention of karst risk in underground gypsum quarries. Tunn. Undergr. Space Technol. 2019, 95, 10313. [Google Scholar] [CrossRef]
  16. Park, B.; Kim, J.; Lee, J.; Kang, M.-S.; An, Y.-K. Underground Object Classification for Urban Roads Using Instantaneous Phase Analysis of Ground-Penetrating Radar (GPR) Data. Remote Sens. 2018, 10, 1417. [Google Scholar] [CrossRef] [Green Version]
  17. Xu, X.; Brekke, C.; Doulgeris, A.P.; Melandsø, F. Numerical Analysis of Microwave Scattering from Layered Sea Ice Based on the Finite Element Method. Remote Sens. 2018, 10, 1332. [Google Scholar] [CrossRef] [Green Version]
  18. Zhang, W.Y.; Hao, T.; Chang, Y.; Zhao, Y. Time-frequency analysis of enhanced GPR detection of RF tagged buried plastic pipes. NDT E Int. 2017, 92, 88–96. [Google Scholar] [CrossRef]
  19. Langan, R.T. Seismic tomography: The accurate and efficient tracing of rays through heterogeneous media. Seg Tech. Program Expand. Abstr. 1984, 3, 856. [Google Scholar]
  20. Goodman, D. Ground-penetrating radar simulation in engineering and archaeology. Geophysics 1994, 59, 224–232. [Google Scholar] [CrossRef]
  21. Yu, J.; Sheng, X. Calculation of electromagnetic scattering from subsurface three-dimensional targets by moment method. J. Electron. Inf. Technol. 2006, 28, 950–954. [Google Scholar]
  22. Yee, K.S. Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media. IEEE Trans. Antennas Propag. 1966, 14, 302–307. [Google Scholar]
  23. Li, Y.; Zhao, Z.; Xu, W.; Liu, Z.; Wang, X. An effective FDTD model for GPR to detect the material of hard objects buried in tillage soil layer. Soil Tillage Res. 2019, 195, 104353. [Google Scholar] [CrossRef]
  24. Alsharahi, G.; Faize, A.; Mostapha, A.M.M.; Driouach, A. 2D FDTD Simulation to Study Response of GPR Signals in Homogeneous and Inhomogeneous Mediums. Int. J. Commun. Antenna Propag. 2016, 6, 153. [Google Scholar] [CrossRef]
  25. Os, A.; Jlv, A.; Hv, B. Two-step perfectly matched layer for arbitrary-order pseudo-spectral analytical time-domain methods. Comput. Phys. Commun. 2019, 235, 102–110. [Google Scholar]
  26. Yang, M.; Fang, H.; Wang, F.; Jia, H.; Lei, J.; Zhang, D. The Three Dimension First-Order Symplectic Partitioned Runge-Kutta Scheme Simulation for GPR Wave Propagation in Pavement Structure. IEEE Access 2019, 7, 151705–151712. [Google Scholar] [CrossRef]
  27. Fang, H.; Gao, L. Symplectic partitioned Runge-Kutta methods for two-dimensional numerical model of ground penetrating radar. Comput. Geosci. 2012, 49, 323–329. [Google Scholar] [CrossRef]
  28. Lei, J.; Wang, Z.; Fang, H.; Ding, X.; Zhang, X.; Yang, M.; Wang, H. Analysis of GPR Wave Propagation in Complex Underground Structures Using CUDA-Implemented Conformal FDTD Method. Int. J. Antennas Propag. 2019, 2019, 5043028. [Google Scholar] [CrossRef]
  29. Yang, Y.; Chen, R.S.; Ye, Z.B.; Wang, Z.B. Analysis of Planar Antennas Using Unconditionally Stable Three-Dimensional ADI-FDTD Method. In Proceedings of the 2005 IEEE Antennas and Propagation Society International Symposium, Washington, DC, USA, 3–8 July 2005; Volume 1B, pp. 170–173. [Google Scholar]
  30. Perhirin, S.; Auffret, Y. A low consumption electronic system developed for a 10km long all-optical extension dedicated to sea floor observatories using power-over-fiber technology and SPI protocol. Microw. Opt. Technol. Lett. 2013, 55, 2562–2568. [Google Scholar]
  31. Yang, Y.; Chen, R.S.; Tang, W.C.; Sha, K.; Yung, E.K.N. Analysis of planar circuits using an unconditionally stable 3D ADI-FDTD method. Microw. Opt. Technol. Lett. 2005, 46, 175–179. [Google Scholar] [CrossRef]
  32. Ahmed, I.; Chen, Z. A hybrid ADI-FDTD subgridding scheme for efficient electromagnetic computation. Int. J. Numer. Model. Electron. Netw. Devices Fields 2004, 17, 237–249. [Google Scholar] [CrossRef]
  33. Zivanovic, S.S.; Yee, K.S.; Mei, K.K. A subgridding method for the time-domain finite-difference method to solve Maxwell’s equations. IEEE Trans. Microw. Theory Tech. 2002, 39, 471–479. [Google Scholar] [CrossRef]
  34. Yu, W.; Mittra, R. A conformal FDTD algorithm for modeling perfectly conducting objects with curve-shaped surfaces and edges. Microw. Opt. Technol. Lett. 2015, 27, 136–138. [Google Scholar] [CrossRef]
  35. Cabello, M.R.; Angulo, L.D.; Alvarez, J.; Bretones, A.R.; Gutierrez, G.G.; Garcia, S.G. A New Efficient and Stable 3D Conformal FDTD. IEEE Microw. Wirel. Compon. Lett. 2016, 26, 553–555. [Google Scholar] [CrossRef] [Green Version]
  36. Zagorodnov, I.; Schuhmann, R.; Weiland, T. Conformal FDTD-methods to avoid time step reduction with and without cell enlargement. J. Comput. Phys. 2007, 225, 1493–1507. [Google Scholar] [CrossRef]
  37. Namiki, T. A new FDTD algorithm based on alternating-direction implicit method. Microw. Theory Tech. IEEE Trans. 1999, 47, 2003–2007. [Google Scholar] [CrossRef]
  38. Feng, D.S.; Dai, Q.W. GPR numerical simulation of full wave field based on UPML boundary condition of ADI-FDTD. NDT E Int. 2011, 44, 495–504. [Google Scholar]
Figure 1. Grid dissection of the circular structure: (a) Actual grid dissection of a circular structure; (b) Grid dissection of circular structure by ladder approximation method; (c) Dissection of circular structure by conformal grid technology.
Figure 1. Grid dissection of the circular structure: (a) Actual grid dissection of a circular structure; (b) Grid dissection of circular structure by ladder approximation method; (c) Dissection of circular structure by conformal grid technology.
Applsci 12 05219 g001
Figure 2. Equivalent schematic of the media parameters at the conformal grid points: (a) Actual grids; (b) Conformal grids.
Figure 2. Equivalent schematic of the media parameters at the conformal grid points: (a) Actual grids; (b) Conformal grids.
Applsci 12 05219 g002
Figure 3. The waveform of the Ricker wavelet.
Figure 3. The waveform of the Ricker wavelet.
Applsci 12 05219 g003
Figure 4. Snapshots of the wave field of Ez at different times: (a) At 2 ns; (b) At 3 ns; (c) At 4 ns; (d) At 5 ns.
Figure 4. Snapshots of the wave field of Ez at different times: (a) At 2 ns; (b) At 3 ns; (c) At 4 ns; (d) At 5 ns.
Applsci 12 05219 g004
Figure 5. Schematic diagram of the circular cavity model.
Figure 5. Schematic diagram of the circular cavity model.
Applsci 12 05219 g005
Figure 6. Comparison of Ez field distribution between ADI−FDTD and FDTD simulation.
Figure 6. Comparison of Ez field distribution between ADI−FDTD and FDTD simulation.
Applsci 12 05219 g006
Figure 7. GPR B-scan images of the two-layer cavity model: (a) Obtained by conformal ADI-FDTD method; (b) Obtained by conformal FDTD method; (c) Obtained by non-conformal ADI-FDTD method; (d) Obtained by non-conformal FDTD method.
Figure 7. GPR B-scan images of the two-layer cavity model: (a) Obtained by conformal ADI-FDTD method; (b) Obtained by conformal FDTD method; (c) Obtained by non-conformal ADI-FDTD method; (d) Obtained by non-conformal FDTD method.
Applsci 12 05219 g007aApplsci 12 05219 g007b
Figure 8. Schematic diagram of model with cavity disease on the upper part of pipeline.
Figure 8. Schematic diagram of model with cavity disease on the upper part of pipeline.
Applsci 12 05219 g008
Figure 9. Schematic diagram of model with cavity disease on the lower part of pipeline.
Figure 9. Schematic diagram of model with cavity disease on the lower part of pipeline.
Applsci 12 05219 g009
Figure 10. GPR B-scan image of a disease model with an irregular cavity above pipeline: (a) Obtained by non-conformal method; (b) Obtained by the conformal method.
Figure 10. GPR B-scan image of a disease model with an irregular cavity above pipeline: (a) Obtained by non-conformal method; (b) Obtained by the conformal method.
Applsci 12 05219 g010
Figure 11. GPR B-scan image of a disease model with an irregular cavity below pipeline: (a) Obtained by non-conformal method; (b) Obtained by the conformal method.
Figure 11. GPR B-scan image of a disease model with an irregular cavity below pipeline: (a) Obtained by non-conformal method; (b) Obtained by the conformal method.
Applsci 12 05219 g011
Figure 12. Underground multi-pipe model with different pipe diameters.
Figure 12. Underground multi-pipe model with different pipe diameters.
Applsci 12 05219 g012
Figure 13. Underground multi-pipe model with different pipe materials.
Figure 13. Underground multi-pipe model with different pipe materials.
Applsci 12 05219 g013
Figure 14. Simulation results of underground multi−pipe model with different pipe diameters.
Figure 14. Simulation results of underground multi−pipe model with different pipe diameters.
Applsci 12 05219 g014
Figure 15. Simulation results of underground multi−pipe model with different pipe materials.
Figure 15. Simulation results of underground multi−pipe model with different pipe materials.
Applsci 12 05219 g015
Figure 16. Schematic diagram of complex underground cavity model.
Figure 16. Schematic diagram of complex underground cavity model.
Applsci 12 05219 g016
Figure 17. GPR B−scan images of underground cavity mode.
Figure 17. GPR B−scan images of underground cavity mode.
Applsci 12 05219 g017
Figure 18. Road radar rectangular cavity experimental diagram.
Figure 18. Road radar rectangular cavity experimental diagram.
Applsci 12 05219 g018
Figure 19. Two-dimensional experimental model.
Figure 19. Two-dimensional experimental model.
Applsci 12 05219 g019
Figure 20. Single−channel wave comparison diagram of simulated and measured results.
Figure 20. Single−channel wave comparison diagram of simulated and measured results.
Applsci 12 05219 g020
Table 1. Computational resource usage of different methods.
Table 1. Computational resource usage of different methods.
MethodMemory (Mb)Δt (ns)IterationsTime (s)
Non-conformal FDTD10.160.0150001849
Conformal FDTD10.200.0150001897
Non-conformal ADI-FDTD22.390.0150004765
Conformal ADI-FDTD-CFLN = 122.480.0150004812
Conformal ADI-FDTD-CFLN = 222.480.0225002403
Conformal ADI-FDTD-CFLN = 322.480.0316671608
Conformal ADI-FDTD-CFLN = 422.480.0412501196
Conformal ADI-FDTD-CFLN = 522.480.051000961
Table 2. Conductivity and dielectric constant of different media.
Table 2. Conductivity and dielectric constant of different media.
Materialεrσ (S/m)
Air10
Clay120.002
PVC plastic30
Concrete60.001
Metal11.0 × 106
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Li, Y.; Wang, N.; Lei, J.; Wang, F.; Li, C. Modeling GPR Wave Propagation in Complex Underground Structures Using Conformal ADI-FDTD Algorithm. Appl. Sci. 2022, 12, 5219. https://0-doi-org.brum.beds.ac.uk/10.3390/app12105219

AMA Style

Li Y, Wang N, Lei J, Wang F, Li C. Modeling GPR Wave Propagation in Complex Underground Structures Using Conformal ADI-FDTD Algorithm. Applied Sciences. 2022; 12(10):5219. https://0-doi-org.brum.beds.ac.uk/10.3390/app12105219

Chicago/Turabian Style

Li, Yinping, Niannian Wang, Jianwei Lei, Fuming Wang, and Ce Li. 2022. "Modeling GPR Wave Propagation in Complex Underground Structures Using Conformal ADI-FDTD Algorithm" Applied Sciences 12, no. 10: 5219. https://0-doi-org.brum.beds.ac.uk/10.3390/app12105219

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop