ABSTRACT

As a highly efficient absorbing boundary condition, Perfectly Matched Layer (PML) has been widely used in Finite Difference Time Domain (FDTD) simulation of Ground Penetrating Radar (GPR) based on the first order electromagnetic wave equation. However, the PML boundary condition is difficult to apply in GPR Finite Element Time Domain (FETD) simulation based on the second order electromagnetic wave equation. This paper developed a non-split perfectly matched layer (NPML) boundary condition for GPR FETD simulation based on the second order electromagnetic wave equation. Taking two-dimensional TM wave equation as an example, the second order frequency domain equation of GPR was derived according to the definition of complex extending coordinate transformation. Then it transformed into time domain by means of auxiliary differential equation method, and its FETD equation is derived based on Galerkin method. On this basis, a GPR FETD forward program based on NPML boundary condition is developed. The merits of NPML boundary condition are certified by compared with wave field snapshots, signal and reflection errors of homogeneous medium model with split and non-split PML boundary conditions. The comparison demonstrated that the NPML algorithm can reduce memory occupation and improve calculation efficiency. Furthermore, numerical simulation of a complex model verifies the good absorption effects of the NPML boundary condition in complex structures.

Keywords:

Non-Split Perfectly Matched Layer (NPML), Ground Penetrating Radar (GPR), Second Order Wave Equation, Finite Element Time Domain (FETD)

1. Introduction

Numerical modeling of ground penetrating radar (GPR) is an important means to study high-frequency electromagnetic wave detection, which plays a critical role in theoretical research of electromagnetic wave propagation in underground structures and guiding for processing and interpretation of actual data [1] . In the field of GPR forward modeling, finite-difference time-domain (FDTD) method is very popular owing to its advantages, including direct calculation in time-domain, easy programming, saving memory space and calculation time [2] [3] [4] [5] . However, it is difficult for FDTD to solve the problems such as undulating surface, free boundary and serious numerical dispersion. Compared with FDTD, pseudo spectral time-domain (PSTD) method is a global calculation method characterized by its high calculation accuracy and weak numerical dispersion [6] [7] [8] . However, PSTD would result in numerical instabilities when dealing with the boundary problems of strongly inhomogeneous medium, and it is also difficult to solve the problems of undulating surface and free boundary condition [9] . In recent years, finite element time-domain (FETD) method has been introduced into the research of GPR numerical modeling as it can well represent the geometric models of complex medium and structures and meet free boundary condition naturally. Shen et al. [10] first derived the finite element equation of GPR wave in low-loss medium and carried out the forward modeling of step model. Di and Wang [11] derived the GPR finite element wave equation with an attenuation term and achieved the FDTD forward modeling of GPR for complex medium. Di and Wang [12] carried out the forward modeling and migration imaging of GPR in complex medium using FETD method. Xi et al. [13] constructed the FETD with 20-node isoparametric elements to simulate Pulse-type GPR; Lu et al. [14] applied FETD to the forward modeling of GPR in Debye dispersive medium. Dai et al. [15] applied the FETD based on biquadratic interpolation and triangulation to GPR forward modeling, achieving higher simulation accuracy than FETD based on linear interpolation and rectangular partition. Feng et al. [5] proposed a hybrid boundary condition suitable for FETD by analyzing the respective advantages of transmitting and Sarma boundary conditions, improving the accuracy of GPR forward modeling. Varela et al. [16] used the FETD algorithm to simulate and analyze the electromagnetic response to the detection of GPR on concrete bridge. Feng et al. [17] applied the wavelet interpolation basis function to the FETD algorithm and carried out the GPR FETD forward modeling, demonstrating the advantages of the algorithm with numerical examples. Zarei et al. [18] applied the orthogonal interpolation basis function to the FETD algorithm, showing that FETD is easier to eliminate the influence of numerical dispersion than FDTD with numerical examples. Feng et al. [19] [20] derived the FETD forward algorithm and FETD/FDTD hybrid algorithm with irregular quadrilateral meshes, achieving high-precision forward modeling of GPR in complex geologic bodies.

When using a computer to carry out GPR numerical simulation, due to limited memory space, mesh space is always truncated somewhere to form a finite region. Then strong non-physical electromagnetic reflection interference waves will be generated at the truncations of the mesh space. Therefore, the truncation boundaries must be properly treated to eliminate or weaken such spurious boundary reflections. To this end, predecessors have done a lot of research work and developed many kind of boundary conditions, such as Sarma boundary condition [5] [21] , Mur boundary condition [22] [23] [24] , paraxial approximation boundary condition [25] [26] , super absorbing boundary condition [27] [28] , perfectly matched layer (PML) boundary condition [29] , etc. Among them, the PML boundary condition has the optimal absorption effect and been widely used in GPR forward modeling with FDTD method [30] - [36] . However, the theory of PML proposed by Berenger does not follow from Maxwell's equations and its physical mechanism is fuzzy. In addition, the calculation of electromagnetic field splitting based on PML boundary condition increases the computational memory and difficulties of numerical implementation [34] [37] . Moreover, the PML boundary condition only has a good effect on absorbing traveling wave but has a poor effect on absorbing low-frequency wave, grazing angle wave and evanescent wave with small incident angles. In order to improve the absorption effects of PML boundary condition and reduce the difficulties in numerical implementation, scholars have done a lot of researches on PML absorbing boundary condition and developed uniaxial anisotropic PML (UPML) [14] [38] [39] [40] [41] [42] , convolution PML (CPML) [2] [43] [44] [45] and complex frequency shift PML (CFS-PML) (Complex Frequency Shifted PML, CFS-PML) [34] [43] [46] [47] [48] .

It should be noted that the aforementioned PMLs and improved boundary conditions mostly are proposed for electromagnetic field modeling based on FDTD method. FDTD is a numerical calculation method based on first-order electromagnetic wave equations while FETD is a method based on second-order equations. Therefore, the PML boundary conditions widely used in FDTD cannot be directly applied to FETD. To solve this problem, Komatitsch and Tromp [49] firstly derived the split PML absorbing boundary condition for the second-order elastic wave equations by referring the idea of constructing PML boundary condition for the first-order elastic wave equations in the elastic wave forward modeling based on FETD. Numerical examples showed that good absorption effects can be achieved at truncated artificial boundaries by incorporating the PML boundary condition into the FETD equations by using a variational formulation. Liu et al. [50] [51] and Liu et al. [52] applied the PML boundary condition to the elastic wave FETD forward modeling based on triangular mesh. However, a third-order time derivative term is included in the PML boundary condition equations and the elastic wave field needs to be split, greatly increasing the computational time and memory space. Basu et al. [53] [54] proposed a non-split PML boundary condition suitable for second-order elastic wave FETD forward modeling with high calculation accuracy and great improvement in calculation efficiency. Currently, the PML boundary condition based on second-order wave equations has been widely used in FETD numerical modeling of acoustic wave [55] [56] [57] , elastic wave [50] [52] [58] [59] and surface wave [60] , and has been continuously developed [61] [62] .

Based on the researches of Komatitsch, Tromp, Basu, Matzen, Liu, etc., and considering the FETD’s advantages of easiness to deal with irregular meshes and free surface boundary conditions as well as flexible representation of irregular boundaries and complex geometries, we propose an efficient non-split PML boundary condition by using complex stretching coordinate transformation and auxiliary variables. Furthermore, we give the computation format of FETD in time domain by using Galerkin approximation technique and develop the FETD forward modeling algorithm of second-order GPR wave equation based on the non-split PML boundary condition. The numerical examples show that, compared with the split PML boundary condition, the non-split PML boundary condition has made up the deficiency of splitting in electromagnetic wave field, ensuring the calculation accuracy and improving the calculation efficiency of full-wave-fields numerical modeling.

2. Second-Order Electromagnetic Wave Equations in Complex Stretching Coordinate System

According to the electromagnetic wave theory, the propagation of GPR electromagnetic wave in underground medium satisfies the Maxwell’s equations [5] . Considering the strike direction of the two-dimensional geo-electric model is along the z-axis, the second-order time-domain electromagnetic wave equation can be expressed as follows:

$\mu \epsilon \frac{{\partial }^2{E}_z}{\partial t^2}+\mu \sigma \frac{\partial E_z}{\partial t}=\frac{{\partial }^2{E}_z}{\partial x^2}+\frac{{\partial }^2{E}_z}{\partial y^2}.$ (1)

where, $E_z$ is the electric field strength (V/m) in the z direction; $H_x,H_y$ respectively are the magnetic field strength (A/m) in the x, y direction; t is the time (s), $\mu ,\sigma ,\epsilon$ respectively are the permeability (H/m), electrical conductivity (S/m) and dielectric constant (F/m) of the medium.

According to the PML boundary condition theory [63] , by introducing the complex stretching coordinate system into the frequency domain, the PML boundary condition can be expressed as follows:

$\tilde{p}=p-\frac{i}{\omega }{\displaystyle {\int }_0^p{d}_p\left(s\right)\text{d}s},\left(p=x,p=y\right).$ (2)

where, $\tilde{p}$ is the complex coordinate, $i=\sqrt{-1}$ is the imaginary unit, $d_p\left(p\right)$ is the boundary attenuation coefficient, which is a real function attenuating with the coordinate p, and $\omega$ is the angular frequency.

From Equation (2), we can respectively obtain the first and second order partial derivatives relation for p and $\tilde{p}$ as follows:

$\frac{\partial }{\partial \tilde{p}}=\frac{i\omega }{i\omega +d_p}\frac{\partial }{\partial p},\frac{{\partial }^2}{\partial {\tilde{p}}^2}={\left(\frac{i\omega }{i\omega +d_p}\right)}^2\frac{{\partial }^2}{\partial p^2}-\frac{{\left(i\omega \right)}^2{d^{\prime }}_p}{i\omega +d_p}\frac{\partial }{\partial p}.$ (3)

Transforming Equation (1) to the frequency domain, we can obtain the following equation:

$\frac{{\partial }^2{\tilde{E}}_z}{\partial x^2}+\frac{{\partial }^2{\tilde{E}}_z}{\partial y^2}=\left(-{\omega }^2\mu \epsilon +i\omega \mu \sigma \right){\tilde{E}}_z.$ (4)

where ${\tilde{E}}_z$ is a Fourier transform of $E_z$ with respect to time. Substituting Equation (3) into Equation (4), we can derive the second-order electromagnetic wave equation in the complex stretching coordinate system as follows:

$\frac{i\omega }{i\omega +d_x}\frac{\partial }{\partial x}\left(\frac{i\omega }{i\omega +d_x}\frac{\partial {\tilde{E}}_z}{\partial x}\right)+\frac{i\omega }{i\omega +d_y}\frac{\partial }{\partial y}\left(\frac{i\omega }{i\omega +d_y}\frac{\partial {\tilde{E}}_z}{\partial y}\right)=\left(-{\omega }^2\mu \epsilon +i\omega \mu \sigma \right){\tilde{E}}_z$ (5)

Deriving Equation (5) we can obtain:

$\begin{array}{l}{\left(\frac{i\omega }{i\omega +d_x}\right)}^2\frac{{\partial }^2{\tilde{E}}_z}{\partial x^2}+\frac{{\omega }^2{d^{\prime }}_x}{{\left(i\omega +d_x\right)}^3}\frac{\partial {\tilde{E}}_z}{\partial x}+{\left(\frac{i\omega }{i\omega +d_y}\right)}^2\frac{{\partial }^2{\tilde{E}}_z}{\partial y^2}+\frac{{\omega }^2{d^{\prime }}_y}{{\left(i\omega +d_y\right)}^3}\frac{\partial {\tilde{E}}_z}{\partial y}\\ =\left(-{\omega }^2\mu \epsilon +i\omega \mu \sigma \right){\tilde{E}}_z\end{array}$ (6)

3. PML Boundary Condition

3.1. Traditional Split PML Boundary Condition

In the second-order electromagnetic wave Equation (6), the electromagnetic field is split into four terms under the traditional split PML boundary condition:

${\tilde{E}}_z={\tilde{E}}_{z,1}+{\tilde{E}}_{z,2}+{\tilde{E}}_{z,3}+{\tilde{E}}_{z,4},$ (7)

The split results corresponding to Equation (6) can be expressed as:

$\{\begin{array}{l}\left(-{\omega }^2\mu \epsilon +i\omega \mu \sigma \right){\tilde{E}}_{z,1}={\left(\frac{i\omega }{i\omega +d_x}\right)}^2\frac{{\partial }^2{\tilde{E}}_z}{\partial x^2},\\ \left(-{\omega }^2\mu \epsilon +i\omega \mu \sigma \right){\tilde{E}}_{z,2}=\frac{{\omega }^2{d}^{\prime }\left(x\right)}{{\left(i\omega +d_x\right)}^3}\frac{\partial {\tilde{E}}_z}{\partial x},\\ \left(-{\omega }^2\mu \epsilon +i\omega \mu \sigma \right){\tilde{E}}_{z,3}={\left(\frac{i\omega }{i\omega +d_y}\right)}^2\frac{{\partial }^2{\tilde{E}}_z}{\partial y^2},\\ \left(-{\omega }^2\mu \epsilon +i\omega \mu \sigma \right){\tilde{E}}_{z,4}+\frac{{\omega }^2{d^{\prime }}_y}{{\left(i\omega +d_y\right)}^3}\frac{\partial {\tilde{E}}_z}{\partial y}.\end{array}$ (8)

Performing an inverse Fourier transform on both side of Equation (8) with respect to $\omega$ . Then the time domain wave equation satisfied in the PML region is as follows:

$\{\begin{array}{l}\left(\mu \epsilon +\mu \sigma {\displaystyle {\int }_0^t}\right){\left({\partial }_t+d_x\right)}^2{E}_{z,1}=\frac{{\partial }^2{E}_z}{\partial x^2}\\ \left(\mu \epsilon +\mu \sigma {\displaystyle {\int }_0^t}\right){\left({\partial }_t+d_x\right)}^3{E}_{z,2}={d^{\prime }}_x\frac{\partial E_z}{\partial x}\\ \left(\mu \epsilon +\mu \sigma {\displaystyle {\int }_0^t}\right){\left({\partial }_t+d_x\right)}^2{E}_{z,3}=\frac{{\partial }^2{E}_z}{\partial y^2}\\ \left(\mu \epsilon +\mu \sigma {\displaystyle {\int }_0^t}\right){\left({\partial }_t+d_x\right)}^3{E}_{z,4}={d^{\prime }}_y\frac{\partial E_z}{\partial y}\end{array}$ (9)

Equation (9) is the second-order GPR wave equation based on the split PML boundary condition. Where, the second and fourth equations need to calculate the third derivative of the electric field ${\tilde{E}}_z$ with respect to time and first integral, which will take more computation time. To avoid the calculation of the third derivative of ${\tilde{E}}_z$ with respect to time, the intermediate variables ${\tilde{P}}_x,{\tilde{P}}_y$ are introduced. Letting $\left({\partial }_t+d_x\right){\tilde{E}}_{z,2}={\tilde{P}}_x,\left({\partial }_t+d_y\right){\tilde{E}}_{z,4}={\tilde{P}}_y$ , the second and fourth equations in Equation (9) can be rewritten as:

$\{\begin{array}{l}\left(\mu \epsilon +\mu \sigma {\displaystyle {\int }_0^t}\right){\left({\partial }_t+d_x\right)}^2{P}_x={d^{\prime }}_x\frac{\partial E_z}{\partial x}\\ \left(\mu \epsilon +\mu \sigma {\displaystyle {\int }_0^t}\right){\left({\partial }_t+d_y\right)}^2{P}_y={d^{\prime }}_y\frac{\partial E_z}{\partial y}\end{array}$ (10)

3.2. Improved Split PML Boundary Condition

To reduce the terms of split in electromagnetic field, the first two terms and the last two terms of Equation (6) respectively are combined into one term [50] . Meanwhile the intermediate variables ${\tilde{P}}_x,{\tilde{P}}_y$ are introduced. Then the split result in Equation (6) can be expressed as：

$\{\begin{array}{l}\left(\mu \epsilon +\mu \sigma /i\omega \right){\left(i\omega +d_x\right)}^2{\tilde{E}}_{1,z}=\frac{{\partial }^2{\tilde{E}}_z}{\partial x^2}-{d^{\prime }}_x{\tilde{P}}_x\\ \left(\mu \epsilon +\mu \sigma /i\omega \right){\left(i\omega +d_y\right)}^2{\tilde{E}}_{2,z}=\frac{{\partial }^2{\tilde{E}}_z}{\partial y^2}-{d^{\prime }}_y{\tilde{P}}_y\\ \left(i\omega +d\left(x\right)\right){\tilde{P}}_x=\frac{\partial {\tilde{E}}_z}{\partial x}\\ \left(i\omega +d\left(y\right)\right){\tilde{P}}_y=\frac{\partial {\tilde{E}}_z}{\partial y}\end{array}$ (11)

According to the Fourier transform theory, the time domain signal $U\left(t\right)$ and its frequency signal $\tilde{U}\left(\omega \right)$ satisfied that $U\left(t\right)=i\omega \cdot \tilde{U}\left(\omega \right)$ , applied it into Equation (11), we can obtain Equation (12) as follows:

$\{\begin{array}{l}\left(\mu \epsilon +\mu \sigma {\displaystyle {\int }_0^t}\right){\left({\partial }_t+d_x\right)}^2{E}_{1,z}=\frac{{\partial }^2{E}_z}{\partial x^2}-{d^{\prime }}_x{P}_x\\ \left(\mu \epsilon +\mu \sigma {\displaystyle {\int }_0^t}\right){\left({\partial }_t+d_y\right)}^2{E}_{2,z}=\frac{{\partial }^2{E}_z}{\partial y^2}-{d^{\prime }}_y{P}_y\\ \left({\partial }_t+d_x\right)P_x=\frac{\partial E_z}{\partial x}\\ \left({\partial }_t+d_y\right)P_y=\frac{\partial E_z}{\partial y}\end{array}$ (12)

The comparison between Equation (9) and Equation (12) shows that the improved split PML boundary condition has the advantages over the traditional split PML boundary condition: 1) avoiding the calculation of third derivative with respect to time by splitting the displacement into two terms in the absorbing layer; 2) reducing calculation amount and improving calculation speed since the electromagnetic field is only split into two terms.

3.3. Non-Split PML Boundary Condition

To avoid electromagnetic field splitting in Equation (5), the variables ${\beta }_x,{\beta }_y$ are introduced [59] [62] and defined as:

${\beta }_x=1+d_x/i\omega ,{\beta }_y=1+d_y/i\omega .$ (13)

The following identities can be obtained by deriving Equation (13):

$\{\begin{array}{l}{\left(i\omega \right)}^2{\beta }_x{\beta }_y={\left(i\omega \right)}^2+i\omega \left(d_x+d_y\right)+d_x{d}_y\\ {\beta }_y/{\beta }_x=1+\left(d_y-d_x\right)/\left(i\omega +d_x\right)\\ {\beta }_x/{\beta }_y=1+\left(d_x-d_y\right)/\left(i\omega +d_y\right)\end{array}$ (14)

Multiplying both sides of Equation (5) by ${\beta }_x{\beta }_y$ and substituting Equation (14) into it, we can obtain:

$\begin{array}{l}\left({\left(i\omega \right)}^2+i\omega \left(d_x+d_y\right)+d_x{d}_y\right)\mu \epsilon {\tilde{E}}_z+\left(i\omega +\left(d_x+d_y\right)+\frac{d_x{d}_y}{i\omega }\right)\mu \sigma {\tilde{E}}_z\\ =\frac{{\partial }^2{\tilde{E}}_z}{\partial x^2}+\frac{\partial }{\partial x}\left(\frac{d_y-d_x}{i\omega +d_x}\frac{\partial {\tilde{E}}_z}{\partial x}\right)+\frac{{\partial }^2{\tilde{E}}_z}{\partial y^2}+\frac{\partial }{\partial y}\left(\frac{d_x-d_y}{i\omega +d_y}\frac{\partial {\tilde{E}}_z}{\partial y}\right)\end{array}$ (15)

Introducing the intermediate variables $P_x,P_y,Q$ , defined as:

${\tilde{P}}_x=\frac{d_y-d_x}{i\omega +d_x}\frac{\partial {\tilde{E}}_z}{\partial x},P_y=\frac{d_x-d_y}{i\omega +d_y}\frac{\partial {\tilde{E}}_z}{\partial y},Q=\frac{d_x{d}_y}{i\omega }{\tilde{E}}_z.$ (16)

Substituting Equation (16) into Equation (15) and deriving it, we can obtain:

$\begin{array}{l}\left[{\left(i\omega \right)}^2+i\omega \left(d_x+d_y\right)+d_x{d}_y\right]\mu \epsilon {\tilde{E}}_z+\left[i\omega +\left(d_x+d_y\right)\right]\mu \sigma {\tilde{E}}_z+\mu \sigma \tilde{Q}\\ =\frac{{\partial }^2{\tilde{E}}_z}{\partial x^2}+\frac{\partial {\tilde{P}}_x}{\partial x}+\frac{{\partial }^2{\tilde{E}}_z}{\partial y^2}+\frac{\partial {\tilde{P}}_y}{\partial y}\end{array}$ (17)

Performing an inverse Fourier transform on Equation (16) and Equation (17) to time domain, we can obtain:

$\{\begin{array}{l}\mu \epsilon \frac{{\partial }^2{E}_z}{\partial t^2}+\left[\left(d_x+d_y\right)\mu \epsilon +\mu \sigma \right]\frac{\partial E_z}{\partial t}+\left[d_x{d}_y\mu \epsilon +\left(d_x+d_y\right)\mu \sigma \right]E_z\\ +\mu \sigma Q=\frac{{\partial }^2{E}_z}{\partial x^2}+\frac{\partial P_x}{\partial x}+\frac{{\partial }^2{E}_z}{\partial y^2}+\frac{\partial P_y}{\partial y},\\ \frac{\partial P_x}{\partial t}+d_x{P}_x=\left(d_y-d_x\right)\frac{\partial E_z}{\partial x},\\ \frac{\partial P_y}{\partial t}+d_y{P}_y=\left(d_x-d_y\right)\frac{\partial E_z}{\partial y},\\ \frac{\partial Q}{\partial t}=d_x{d}_y{E}_z.\end{array}$ (18)

Equation (18) is the second-order GPR wave equation based on the non-split PML boundary condition. Comparison among Equations (9), (12) and (18) shows that the non-split PML boundary condition neither needs split the electromagnetic field, nor needs the calculation of third derivative with respect to time, which means less calculation amount and higher calculation speed.

4. Finite Element Time Domain GPR Equation Based on Non-Split PML Boundary Condition

According to the principle of the Galerkin method [51] , the weak form of Equation (18) is:

$\{\begin{array}{l}{\displaystyle {\int }_{\Omega }\frac{{\partial }^2{E}_z}{\partial t^2}\phi \text{d}\Omega }+{\displaystyle {\int }_{\Omega }\left[\left(d_x+d_y\right)+\frac{\sigma }{\epsilon }\right]\frac{\partial E_z}{\partial t}\phi \text{d}\Omega }+{\displaystyle {\int }_{\Omega }\left[d_x{d}_y+\left(d_x+d_y\right)\frac{\sigma }{\epsilon }\right]E_z\phi \text{d}\Omega }\\ +{\displaystyle {\int }_{\Omega }\frac{\sigma }{\epsilon }Q\phi \text{d}\Omega }-{\displaystyle {\int }_{\Gamma }\frac{1}{\mu \epsilon }\left(\frac{\partial E_z}{\partial x}{n}_x+\frac{\partial E_z}{\partial y}{n}_y\right)\phi \text{d}\Gamma }+{\displaystyle {\int }_{\Omega }\frac{1}{\mu \epsilon }\left(\frac{\partial E_z}{\partial x}\frac{\partial \phi }{\partial x}+\frac{\partial E_z}{\partial y}\frac{\partial \phi }{\partial y}\right)\text{d}\Omega }\\ =-{\displaystyle {\int }_{\Omega }\frac{1}{\mu \epsilon }\frac{\partial P_x}{\partial x}\text{d}\Omega }-{\displaystyle {\int }_{\Omega }\frac{1}{\mu \epsilon }\frac{\partial P_y}{\partial y}\text{d}\Omega }\\ {\displaystyle {\int }_{\Omega }\frac{\partial P_x}{\partial t}{\phi }_x\text{d}\Omega }+{\displaystyle {\int }_{\Omega }{d}_x{P}_x\text{d}\Omega }={\displaystyle {\int }_{\Omega }\left(d_y-d_x\right)\frac{\partial E_z}{\partial x}\text{d}\Omega }\\ {\displaystyle {\int }_{\Omega }\frac{\partial P_y}{\partial t}{\phi }_y\text{d}\Omega }+{\displaystyle {\int }_{\Omega }{d}_y{P}_y\text{d}\Omega }={\displaystyle {\int }_{\Omega }\left(d_x-d_y\right)\frac{\partial E_z}{\partial y}\text{d}\Omega }\\ {\displaystyle {\int }_{\Omega }\frac{\partial Q}{\partial t}\phi \text{d}\Omega }={\displaystyle {\int }_{\Omega }{d}_x{d}_y{E}_z\text{d}\Omega }\end{array}$ (19)

where, $\Gamma ,\Omega$ respectively are the calculation region and its boundary. $n={\left(n_x,n_y\right)}^{\text{T}}$ is the outer normal vector of the boundary. $\phi ={\left({\phi }_x,{\phi }_y\right)}^{\text{T}}$ is the virtual displacement vector. Dirichlet boundary condition with a displacement of 0 is used for the PML boundary and free boundary condition is used for the upper boundary. The finite element equation corresponding to Equation (19) can be derived with the Galerkin finite element method:

$\{\begin{array}{l}M{\ddot{E}}_z+M_1{\dot{E}}_z+M_2{E}_z+M_{3}Q=K{E}_z+K_x{P}_x+K_y{P}_y\\ M{\dot{P}}_x+C_x{P}_x=C_{xx}{E}_z\\ M{\dot{P}}_y+C_y{P}_y=C_{yy}{E}_z\\ M\dot{Q}=C_{xy}{E}_z\end{array}$ (20)

Equation (20) is the FETD GPR wave equation based on the non-split PML boundary condition. Where, the coefficient matrixes are expressed as follows:

$\{\begin{array}{l}M={\displaystyle {\int }_{\Omega }N{N}^{\text{T}}\text{d}\Omega },M_1={\displaystyle {\int }_{\Omega }\left(d_x+d_y+\sigma /\epsilon \right)N{N}^{\text{T}}\text{d}\Omega },\\ M_2={\displaystyle {\int }_{\Omega }\left(d_x{d}_y+\left(d_x+d_y\right)\sigma /\epsilon \right)N{N}^{\text{T}}\text{d}\Omega },\\ M_3=\left(\sigma /\epsilon \right){\displaystyle {\int }_{\Omega }N{N}^{\text{T}}\text{d}\Omega },K=\frac{1}{\mu \epsilon }{\displaystyle {\int }_{\Omega }\left(\frac{\partial N}{\partial x}\frac{N^{\text{T}}}{\partial x}+\frac{\partial N}{\partial y}\frac{N^{\text{T}}}{\partial y}\right)\text{d}\Omega },\\ K_x=\frac{1}{\mu \epsilon }{\displaystyle {\int }_{\Omega }\left(\frac{\partial N}{\partial x}\frac{N^{\text{T}}}{\partial x}\right)\text{d}\Omega },K_y=\frac{1}{\mu \epsilon }{\displaystyle {\int }_{\Omega }\left(\frac{\partial N}{\partial y}\frac{N^{\text{T}}}{\partial y}\right)\text{d}\Omega },\\ C_x={\displaystyle {\int }_{\Omega }{d}_{x}N{N}^{\text{T}}\text{d}\Omega },C_y={\displaystyle {\int }_{\Omega }{d}_{y}N{N}^{\text{T}}\text{d}\Omega },\\ C_{xx}={\displaystyle {\int }_{\Omega }\left(d_y-d_x\right)\frac{\partial N}{\partial x}{N}^{\text{T}}\text{d}\Omega },\\ C_{yy}={\displaystyle {\int }_{\Omega }\left(d_x-d_y\right)\frac{\partial N}{\partial y}{N}^{\text{T}}\text{d}\Omega },\\ C_{xy}={\displaystyle {\int }_{\Omega }{d}_x{d}_{y}N{N}^{\text{T}}\text{d}\Omega },\end{array}$ (21)

where, $N$ is the shape function of the calculation region, $M$ is the mass matrix, $C$ is the damping matrix, $S$ is the load vector; $\ddot{U}$ and $\dot{U}$ indicates the first and second derivative of the electric field $E_z$ with respect to time, respectively.

The general form of the equations in Equation (21) is:

$M\ddot{U}+C\dot{U}+KU=S$ (22)

To solve Equation (22), we use the Newmark difference algorithm which is a time integration algorithm with energy conservation. Using the algorithm, we can update the point $\left(U^{n+1},{\dot{U}}^{n+1}\right)$ based on the point $\left(U^n,{\dot{U}}^n\right)$ . And the iterative formula of the Newmark time-stepping algorithm for Equation (22) is as follows:

$\{\begin{array}{l}{U}^{n+1}=U^n+\Delta t{\dot{U}}^n+\Delta t^2\left[\left(0.5-\beta \right)a^n+\beta a^{n+1}\right],\\ {\dot{U}}^{n+1}={\dot{U}}^n+\Delta t\left[\left(1-\gamma \right)a^n+\gamma a^{n+1}\right],\\ a^{n+1}=M^{-1}\left(S-C{\dot{U}}^n-K{U}^n\right).\end{array}$ (23)

The above matrix formed by finite elements is a large sparse matrix which requires a large memory space for storing all elements. In this paper, the stiffness matrix adopts compressed storage row (CSR) format and only needs to store its non-zero elements, which can greatly reduce the storage space needed [51] . In addition, a lumped mass matrix technique is adopted to avoid the inverse operation of the matrix $M$ and improve the calculation efficiency [52] .

5. Calculation Examples

5.1. Analysis of Absorption Effect and Calculation Efficiency of PML

To verify the advantages of the non-split PML boundary condition proposed in this paper, a homogeneous medium model of 2.2 m × 2.2 m is established with a relative dielectric constant of 5.0 and electrical resistivity of 0.001 Ω·m. For simplicity, the model space is discretized using structured 4-node quadrilateral meshes with a total of 220 × 220 meshes with a size of 0.01 m × 0.01 m. The PML absorbing boundary layer around the model is 0.2 m thick, occupying 20 meshes as a unit as shown in Figure 1. The forward modeling calculation takes 30 ns. To ensure stable update of display time, the time step is set to be 0.01 ns. A zero-phase Ricker wavelet with a center frequency of 500 MHz is used as the excitation source and located in the upper left corner of the model (0.3 m, 0.5 m), and the receiver is located at the center of the model (2.3 m, 2.2 m).

The model is performed with FETD forward modeling under PML-free, split PML and non-split PML boundary conditions respectively, and the snapshots of wave fields at different times are obtained as shown in Figures 2-4. Figure 2(a) and Figure 2(b) respectively are the snapshots of GPR wavefields at time of 4 ns and 8.5 ns under PML-free boundary condition. Figure 2(a) shows strong reflections from the upper and left boundaries when propagating at time of 4 ns. At

8.5 ns, the reflections from the upper and left boundaries continue to spread to the middle of the model. Meanwhile, the wave fronts have propagated to the right and bottom boundaries and been reflected. The above non-physical reflected waves seriously interfere with the propagation of GPR waves in the model space.

Figure 3(a) and Figure 3(b), Figure 4(a) and Figure 4(b) are the snapshots of GPR wavefields under the split and non-split PML boundary conditions, respectively. The figures show that there is almost no difference between the split and non-split PML boundary conditions. They are both ideal to absorb waves of different incident angles without obvious spurious reflections near the boundaries. When the GPR waves propagate into the PML layers, there are no visible reflections as shown in the wave field snapshots at 8.5 ns.

Figure 5 shows the variation curves of electric field energy with time in the calculation region under three kinds of boundary conditions. And from 0 ns to 6 ns, the variations of electric field energy with time are identical because the electromagnetic waves only propagate in the calculation region. After 6 ns, when no boundary condition is imposed, electric field energy always fluctuates around 110 J and propagates with time, which will seriously affect the energy distribution in the calculation region. When the PML boundary condition is imposed,

electric field energy gradually propagates into the PML layers after 6 ns and attenuates rapidly to zero after 18 ns. It is clear that both of the PML boundary conditions have good absorption effects.

Figure 6 shows the comparison among the GPR waveforms at the receiver respectively calculated under the split and non-split PML boundary conditions and their reference solution. The reference solution used in this paper is calculated using the same numerical method in a large enough region (to ensure that spurious reflections do not reach the receiver during the modeling time) where no boundary conditions are imposed, and the relative positions of the fluctuation source and the receiver remain unchanged. It can be seen from the figure that calculation results under the two PML boundary conditions are consistent with the reference solution, which means that both of them are relatively accurate for simulating electromagnetic waves. Meanwhile, it also proves that the PML boundary condition built in this paper has good absorption effects without generating spurious reflections from the artificial truncated boundaries.

In order to better analyze the absorption effects of the two PML boundary conditions, we use the formula $Erro{r}_{db}=20\times \log \left(\left|E_S-E_{ref}\right|/\left|E_{ref\max }\right|\right)$ to calculate the reflection errors at the receiver under the two different boundary conditions, as shown in Figure 7. Compared with the split PML boundary condition, the reflection errors under the non-split PML boundary condition are significantly reduced by 20 dB on average.

For better analysis of the calculation efficiency of the split and non-split PML boundary conditions, we program the FETD GPR forward modeling under the two boundary conditions using the Matlab platform and test them on the same computer (Lenovo Think Centre M8300t). The statistics of CPU time consuming and memory occupation under the two boundary conditions is shown as Table 1. Compared with the split PML boundary condition, the calculation speed under the non-split PML boundary condition are doubled and the memory usage is also reduced by about 1/4 times. Thus in conclusion, the non-split PML boundary condition can ensure calculation accuracy, improve calculation efficiency and reduce memory space, which means a significant advantage in rapid and high-precision FETD modeling of complex GPR models.

5.2. Forward Modeling Examples

To verify the absorption effects of the non-split PML boundary condition in the FETD forward modeling of complex GPR model, a complex geo-electric model is established as shown in Figure 8. The model is divided into two layers by an undulating interface, including an upper layer medium of concrete with relative dielectric constant of 5.0 and electrical conductivity of 0.001 S∙m⁻¹ and a lower layer medium of soil with relative dielectric constant of 10.0 and electrical conductivity of 0.002 S∙m⁻¹. In the lower layer medium, three circular anomalies with a radius of 0.1m are buried at different depths. From left to right, their relative dielectric constants and center coordinates are 20, 15, 81 and (1.4 m, 1.0 m), (2.4 m, 1.1 m), (3.4 m, 1.2 m), respectively. The model is discretized by quadrilateral meshes with a total of 480 × 180 meshes and a distance of 0.01 m between each other. In addition, there are 20 PML boundary layers on the outer boundary of the calculation region. At last, a zero-phase Ricker wavelet with a center frequency of 500 MHz is given as an excitation source on the surface, and FETD method based on the non-split PML boundary condition is used for forward modeling with a time step of 0.02 ns and a window length of 40 ns.

Figure 9(a) shows the profile of GPR forward modeling based on the FETD algorithm under the non-split PML boundary condition. In the profile, there are no obvious reflections from the artificial truncation boundaries, indicating that the non-split PML boundary condition has a good effect on absorbing strong reflections from the artificial truncation boundaries. In the vicinity of 10 ns, from left to right, there are three obvious diffraction peaks, which correspond to the two lowest points and one highest point of the undulating surface respectively. The undulation shape of the undulating surface is clearly discernible in the profile, and is consistent with the real shape. From 20 ns to 30 ns, the three circular anomalies in the lower layer medium present obvious hyperbolic diffraction waves which greatly extend on both sides. In addition, the multiple diffractions have occurred due to a large difference of the dielectric constant between

the water pipe medium and the background medium. Based on wide angle method, the profile of GPR forward modeling is shown as Figure 9(b). Compared with Figure 9(a), the shape of the undulating surface is not as obvious as that obtained by profiling method. Moreover, due to the influence of the undulating interface, the occurrence time, shape and vertex position of the hyperbolic reflections from the three circular anomalies in the lower layer medium are different from those obtained by profiling method.

To better understand the absorption effects of non-split PML boundary condition on the strong reflections from model boundaries, the excitation source is positioned in the center of the surface and the wave field snapshots of different times are obtained through the FETD forward modeling, as shown in Figure 10. In the wave field snapshots, the thicknesses of the PMLs on the left, right and lower boundaries are preserved. In Figure 10(a), the electromagnetic waves spread out in semi-concentric circles and just meet the undulating interface. In Figure 10(b), the electromagnetic waves are reflected from the undulating interface and the reflections propagate upward. In Figure 10(c), the wave fronts continue to propagate forward meeting the undulating interface in a gradually increasing range; meanwhile, the reflections from the undulating interface propagate to the upper interface and are fully absorbed after entering the PML layers. In Figure 10(d), due to the undulation of the interface, the wave fronts are irregular and reflections are generated from the circular anomalies at the center when the wave fronts meet them. In Figure 10(e) and Figure 10(f), the

wave fronts propagate to the three circular anomalies, from which the reflected waves spread and propagate; meanwhile, the reflected waves from the undulating interface all have entered the PML layers and have been completely absorbed. In Figure 10(g), the wave fronts meet the truncated boundaries of the model without spurious reflections; meanwhile, the diffraction waves caused by the three circular anomalies are clearly visible. In Figure 10(h), the wave fronts enter the PML layers without spurious reflections, which indicates that the non-split PML boundary condition used in this paper has good absorption effects and can be applied to the FETD forward modeling of complex GPR models.

6. Conclusion

In this paper, based on the second order electromagnetic wave equation, an efficient non-split PML boundary condition is proposed by using the complex stretching coordinates and introducing the auxiliary functions. Furthermore, the FETD GPR wave equation under the non-split PML boundary condition is derived based on the Galerkin method. The forward modeling examples show that spurious reflections from the boundaries can be well absorbed under the non-split PML boundary condition. Compared with the split PML boundary condition, the non-split PML boundary condition can reduce memory occupation and improve calculation efficiency without loss of forward modeling accuracy. Therefore, based on the convenient and efficient non-split PML boundary condition and combined with the FETD method, which can represent complex geometric models easily by using irregular quadrilateral and triangular meshes, we can perform rapid and highly accurate FETD forward modeling of complex GPR models.

Acknowledgements

This research is financially supported by the National Natural Science Foundation of China (41574078, 41604102 and 41604039) and supported by the National Natural Science Foundation of Guangxi (2016GXNSFBA380082, 2016GXNSFBA380215, 2018GXNSFAA138059).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

Cite this paper

Zhang, Z., Wang, H.H., Wang, M.L., Guo, X. and Guo, G.H. (2019) Non-Split PML Boundary Condition for Finite Element Time-Domain Modeling of Ground Penetrating Radar. Journal of Applied Mathematics and Physics, 7, 1077-1096. https://doi.org/10.4236/jamp.2019.75073

References