Full-Wave Modelling of Ground-Penetrating Radars: Antenna Mutual Coupling Phenomena and Sub-Surface Scattering Processes

1. Introduction

Ground-penetrating radar (GPR) technology finds applications in many areas such as geophysical prospecting, archaeology, civil engineering, environmental engineering, and defence applications as a non-invasive sensing tool [3], [6], [18]. One key component in any GPR system is the receiver/transmitter antenna. Desirable features for GPR antennas include efficient radiation of ultra-wideband pulses into the ground, good impedance matching over the operational frequency band, and small size. As the attenuation of radio waves in geophysical media increases with frequency [9], [13], ground-penetrating radars typically operate at frequencies below $1GHz$ [4]. For either impulse [13] or stepped-frequency continuous-wave applications [17], the wider the frequency range, the better the range resolution of the radar. Continuous wave multi-frequency radars are advantageous over impulse radars in coping with dispersion of the medium, the noise level at the receiver end, and the controllability of working frequency. They require, however, mutual coupling between the transmit (Tx) and receive (Rx) antennas, which determines the dynamic range of the system, to be kept as small as possible [12].

In this book chapter, the full-wave analysis of electromagnetic coupling mechanisms between resistively loaded wideband dipole antennas operating in realistic GPR scenarios is carried out. To this end, a locally conformal finite-difference time-domain (FDTD) technique, useful to model electromagnetic structures having complex geometry, is adopted [1], [2]. Such a scheme, necessary to improve the numerical accuracy of the conventional FDTD algorithm [19], [21], by avoiding staircase approximation, is based on the definition of effective material parameters [14], suitable to describe the geometrical and electrical characteristics of the structure under analysis. By doing so, the losses in the soil, as well as the presence of ground-embedded inhomogeneities with arbitrary shape and electrical properties, are properly taken into account. Emphasis is devoted to the investigation of the antenna pair performance for different Tx–Rx separations and elevations over the ground, as well as on scattering from dielectric and metallic pipes buried at different depths and having different geometrical and electrical characteristics. Novelty of the analysis lies in the fact that at the lowest operational frequency both the receive antenna and a pipe are situated in the near-field, whilst at the highest operational frequency only the far field is playing the role. The obtained numerical results provide a physical insight into the underlying mechanisms of subsurface diffraction and antenna mutual coupling processes. This information in turn can be usefully employed to optimize the performance of detection algorithms in terms of clutter rejection.

Finally, a frequency-independent equivalent circuit model of antenna pairs is provided in order to facilitate the design of the RF front-end of ground-penetrating radars by means of suitable software CAD tools. The procedure employed to extract the equivalent circuit is based on a heuristic modification of the Cauer’s network synthesis technique [10] useful to model ohmic and radiation losses. In this way, one can obtain a meaningful description of the natural resonant modes describing the electromagnetic behaviour of antenna pairs for GPR systems.

2. Locally conformal finite-difference time-domain technique

The analysis and design of complex radiating structures requires accurate electromagnetic field prediction models. One such widely used technique is the FDTD algorithm. However, in the conventional formulation proposed by Yee [19], [21], each cell of the computational grid is implicitly supposed to be filled by a homogeneous material. For this reason, the adoption of Cartesian meshes could result in reduced numerical accuracy when structures having curved boundaries have to be modelled. In this case, locally conformal FDTD schemes [1], [2] provide clear advantages over the use of the stair-casing approach or unstructured and stretched space lattices, potentially suffering from significant numerical dispersion and/or instability [19]. Such schemes, necessary to improve the numerical accuracy of the conventional algorithm, are based on the definition of effective material parameters suitable to describe the geometrical and electrical characteristics of the structure under analysis.

In this section, a computationally enhanced formulation of the locally conformal FDTD scheme proposed in [1] is described. To this end, let us consider a three-dimensional domain $D$ filled by a linear, isotropic, non dispersive material, having permittivity$\epsilon (r)$, magnetic permeability$\mu (r)$, and electrical conductivity$\sigma (r)$. In such a domain, a dual-space, non-uniform lattice formed by a primary and secondary mesh is introduced. The primary mesh $M_D$ is composed of space-filling hexahedrons, whose vertices are defined by the Cartesian coordinates:

$$\left\{(x_i,y_j,z_k)|i=\mathrm{0,}\dots ,N_x;j=\mathrm{0,}\dots ,N_y;k=\mathrm{0,}\dots ,N_z\right\}.$$E1

As a consequence, the edge lengths between adjacent vertices in $M_D$ result to be expressed as:

$$\{\begin{array}{ll}\Delta x_i=x_{i+1}-x_i,\hfill & i=\mathrm{0,1,}\dots ,N_x-1,\hfill \\ \Delta y_j=y_{j+1}-y_j,\hfill & j=\mathrm{0,1,}\dots ,N_y-1,\hfill \\ \Delta z_k=z_{k+1}-z_k,\hfill & k=\mathrm{0,1,}\dots ,N_z-1.\hfill \end{array}$$E2

The secondary or dual mesh ${\tilde{M}}_D$ (see Fig. 1) is composed of the closed hexahedrons whose edges penetrate the shared faces of the primary cells and connect the relevant centroids, having coordinates$x_{i+1/2}=x_i+\Delta x_i/2$, $y_{j+1/2}=y_j+\Delta y_j/2$,$z_{k+1/2}=z_k+\Delta z_k/2$. A set of dual edge lengths is then introduced in ${\tilde{M}}_D$ as follows:

$$\{\begin{array}{ll}{\ell }_{x_i}=\left(\Delta x_{i-1}+\Delta x_i\right)/2,\hfill & i=\mathrm{1,2,}\dots ,N_x-1,\hfill \\ {\ell }_{y_j}=\left(\Delta y_{j-1}+\Delta y_j\right)/2,\hfill & j=\mathrm{1,2,}\dots ,N_y-1,\hfill \\ {\ell }_{z_k}=\left(\Delta z_{k-1}+\Delta z_k\right)/2,\hfill & k=\mathrm{1,2,}\dots ,N_z-1.\hfill \end{array}$$E3

As usual, the electric field components are defined along each edge of a primary lattice cell, whereas the magnetic field components are assumed to be located along the edges of the secondary lattice cells. In this formulation, the relationship between $E-$ and $H-$field components is given by Maxwell’s equations expressed in integral form, specifically using Faraday-Neumann’s law and Ampere’s law, respectively. In particular, the enforcement of the Ampere’s law on the generic dual-mesh cell surface ${{\tilde{S}}_x|}_{i+1/\mathrm{2,}j,k}$ having boundary $\partial {{\tilde{S}}_x|}_{i+1/\mathrm{2,}j,k}={{\tilde{C}}_x|}_{i+1/\mathrm{2,}j,k}$ (see Fig. 1) results in the following integral equation:

Figure 1.

Cross-sectional view of the FDTD computational grid in presence of curved boundaries between different dielectric materials.

$${\displaystyle \underset{{{\tilde{C}}_x|}_{i+1/\mathrm{2,}j,k}}{\oint }H(r,t)\cdot dl}={\displaystyle \underset{{{\tilde{S}}_x|}_{i+1/\mathrm{2,}j,k}}{\iint }\sigma (r)E_x(r,t)dS}+\frac{\partial }{\partial t}{\displaystyle \underset{{{\tilde{S}}_x|}_{i+1/\mathrm{2,}j,k}}{\iint }\epsilon (r)E_x(r,t)dS},$$E4

where:

$$\begin{array}{c}{\displaystyle \underset{{{\tilde{C}}_x|}_{i+1/\mathrm{2,}j,k}}{\oint }H\left(r,t\right)\cdot dl}={\ell }_{z_k}\left[H_z\left(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}},z_k,t\right)-H_z\left(x_{i+\frac{1}{2}},y_{j-\frac{1}{2}},z_k,t\right)\right]+\\ -{\ell }_{y_j}\left[H_y\left(x_{i+\frac{1}{2}},y_j,z_{k+\frac{1}{2}},t\right)-H_y\left(x_{i+\frac{1}{2}},y_j,z_{k-\frac{1}{2}},t\right)\right]+o\left({\ell }_{y_j}\right)+o\left({\ell }_{z_k}\right)\end{array}$$E5

as ${\ell }_{y_j}$ and ${\ell }_{z_k}$ tend to zero. Under the assumption that the spatial increments$$\Delta x_i$$, $$\Delta y_j$$, $$\Delta z_k$$of the computational grid are small compared to the minimum working wavelength, the infinitesimal terms of higher order appearing in (5) can be neglected. Furthermore, it should be noticed that the $x-$component of the electric field is continuous along the interfaces crossing ${{\tilde{S}}_x|}_{i+1/\mathrm{2,}j,k}$ so that, under the mentioned hypothesis, the following approximation can be made:

$${\displaystyle \underset{{{\tilde{S}}_x|}_{i+1/\mathrm{2,}j,k}}{\iint }\left\{\begin{array}{c}\epsilon \\ \sigma \end{array}\right\}(r)E_x(r,t)dS}\simeq E_x\left(x_{i+\frac{1}{2}},y_j,z_k,t\right){\displaystyle \underset{{{\tilde{S}}_x|}_{i+1/\mathrm{2,}j,k}}{\iint }\left\{\begin{array}{c}\epsilon \\ \sigma \end{array}\right\}(r)dS}.$$E6

Hence, combining the equations above yields:

$$\begin{array}{c}{{\overline{\epsilon }}_x^{eff}|}_{i+\frac{1}{2},j,k}\frac{\partial }{\partial t}{{ℰ}_x(t)|}_{i+\frac{1}{2},j,k}+{{\overline{\sigma }}_x^{eff}|}_{i+\frac{1}{2},j,k}{{ℰ}_x(t)|}_{i+\frac{1}{2},j,k}=\\ ={{\mathscr{H}}_z(t)|}_{i+\frac{1}{2},j+\frac{1}{2},k}-{{\mathscr{H}}_z(t)|}_{i+\frac{1}{2},j-\frac{1}{2},k}-{{\mathscr{H}}_y(t)|}_{i+\frac{1}{2},j,k+\frac{1}{2}}+{{\mathscr{H}}_y(t)|}_{i+\frac{1}{2},j,k-\frac{1}{2}},\end{array}$$E7

where we have introduced the normalized field quantities:

$${{ℰ}_x(t)|}_{i+\frac{1}{2},j,k}=\Delta x_i{E}_x\left(x_{i+\frac{1}{2}},y_j,z_k,t\right),$$E8

$${{\mathscr{H}}_z(t)|}_{i+\frac{1}{2},j+\frac{1}{2},k}={\ell }_{z_k}{H}_z\left(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}},z_k,t\right),$$E9

$${{\mathscr{H}}_y(t)|}_{i+\frac{1}{2},j,k+\frac{1}{2}}={\ell }_{y_j}{H}_y\left(x_{i+\frac{1}{2}},y_j,z_{k-\frac{1}{2}},t\right),$$E10

and the averaged effective permittivity${{\overline{\epsilon }}_x^{eff}|}_{i+1/\mathrm{2,}j,k}$, and conductivity${{\overline{\sigma }}_x^{eff}|}_{i+1/\mathrm{2,}j,k}$, defined as follows:

$${{\left\{\begin{array}{c}\overline{\epsilon }\\ \overline{\sigma }\end{array}\right\}}_x^{eff}|}_{i+\frac{1}{2},j,k}=\frac{1}{\Delta x_i}{\displaystyle \underset{z_{k-1/2}}{\overset{z_{k+1/2}}{\int }}{\displaystyle \underset{y_{j-1/2}}{\overset{y_{j+1/2}}{\int }}\left\{\begin{array}{c}\epsilon \\ \sigma \end{array}\right\}\left(x_{i+\frac{1}{2}},y,z\right)dydz}}.$$E11

The time derivative in (7) is then evaluated using a central-difference approximation that is second order-accurate if $E-$ and $H-$field components are staggered in time domain [19]. This results in the following explicit time-stepping relation:

$${{ℰ}_x|}_{i+\frac{1}{2},j,k}^{n+1}={{\overline{\alpha }}_x^{(E)}|}_{i+\frac{1}{2},j,k}\text{}{{ℰ}_x|}_{i+\frac{1}{2},j,k}^{n-1}+{{\overline{\beta }}_x^{(E)}|}_{i+\frac{1}{2},j,k}\text{}{{\left(\nabla \times \mathscr{H}\right)}_x|}_{i+\frac{1}{2},j,k}^{n+\frac{1}{2}},$$E12

where:

$${{\left(\nabla \times \mathscr{H}\right)}_x|}_{i+\frac{1}{2},j,k}^{n+\frac{1}{2}}={{\mathscr{H}}_z|}_{i+\frac{1}{2},j+\frac{1}{2},k}^{n+\frac{1}{2}}-{{\mathscr{H}}_z|}_{i+\frac{1}{2},j-\frac{1}{2},k}^{n+\frac{1}{2}}-{{\mathscr{H}}_y|}_{i+\frac{1}{2},j,k+\frac{1}{2}}^{n+\frac{1}{2}}+{{\mathscr{H}}_y|}_{i+\frac{1}{2},j,k-\frac{1}{2}}^{n+\frac{1}{2}}$$E13

denotes the finite-difference expression of the normalized $x-$component of the magnetic field curl at the time step$n+1/2$. In (12), the information regarding the local physical and geometrical properties of the electromagnetic structure under analysis is transferred to the position-dependent coefficients:

$${{\overline{\alpha }}_x^{(E)}|}_{i+\frac{1}{2},j,k}=\frac{1-{{\overline{Q}}_x^{eff}|}_{i+\frac{1}{2},j,k}}{1+{{\overline{Q}}_x^{eff}|}_{i+\frac{1}{2},j,k}},$$E14

$${{\overline{\beta }}_x^{(E)}|}_{i+\frac{1}{2},j,k}=\frac{\Delta t/{{\overline{\epsilon }}_x^{eff}|}_{i+\frac{1}{2},j,k}}{1+{{\overline{Q}}_x^{eff}|}_{i+\frac{1}{2},j,k}},$$E15

with:

$${{\overline{Q}}_x^{eff}|}_{i+\frac{1}{2},j,k}=\frac{{{\overline{\sigma }}_x^{eff}|}_{i+\frac{1}{2},j,k}\Delta t}{2{{\overline{\epsilon }}_x^{eff}|}_{i+\frac{1}{2},j,k}}.$$E16

The update equations of the remaining components of the electric and magnetic field can be easily derived by permuting the spatial indices$$i$$, $$j$$, $$k$$and applying the duality principle in the discrete space.

As it can be readily noticed, the computation of position-dependent coefficients (14)-(16) can be carried out before the FDTD-method time marching starts. As a consequence, unlike in conformal techniques based on stretched space lattices, no additional correction is required in the core of the numerical algorithm. Furthermore, the resulting FDTD update equations (12)-(13) have a very convenient structure, leading to a $14\%$ reduction of the number of floating-point operations needed to determine the unknown field quantities in the generic mesh cell compared to the Yee algorithm [19], [21]. It is also to be pointed out that the proposed scheme has the same numerical stability properties as the conventional FDTD formulation, although it introduces a significant improvement in accuracy over the stair-casing approximation as well as alternative weighted averaging FDTD approaches [7]. In order to assess the effectiveness of the developed technique, several test cases have been considered. For the sake of brevity, only the computation of the fundamental resonant frequency of a metallic cavity loaded with a cylindrical dielectric resonator is presented in the Appendix. The obtained results clearly demonstrate the suitability of the proposed scheme to efficiently handle the problem of modelling antennas for ground-penetrating radar applications, where the accurate characterization of complex metal-dielectric objects having irregular geometry (think about the shape of buried targets and ground-embedded inhomogeneities) is required.

3. The full-wave antenna modelling

It is our intention to focus the attention on the full-wave analysis of a resistively loaded dipole antenna pair located above a ground modelled as a lossy homogeneous half-space having relative permittivity ${\epsilon }_r=6$ and electrical conductivity$\sigma =15mS/m$. The geometry of the structure is depicted in Fig. 2.

Figure 2.

Top and cross-sectional view of a resistively loaded dipole antenna pair located above a lossy homogeneous half space. Structure characteristics:$l_d=40cm$, $D_d=5mm$, $\delta =2.5mm$, ${\epsilon }_r=6$,$\sigma =15mS/m$. The reference system used to express the field quantities is also shown.

The dipoles are denoted as dipole $\#1$ and dipole$\#2$, respectively. In the considered antenna configuration, dipole $\#1$ is driven by a delta-gap voltage source with internal resistance$R_G=R_1=300\Omega$, whereas dipole $\#2$ is closed on a matched load having resistance$R_L=R_2=R_G$. To properly enlarge the antenna bandwidth, thus reducing late-time ringing phenomena, a continuous resistive loading, having Wu-King-like profile [3], [15], [16]:

$${\sigma }_d\left(y\right)={\sigma }_0\left(1-\frac{2\left|y\right|}{l_d}\right),$$E17

has been applied to the flairs of the considered radiators. In (17), ${\sigma }_0$denotes the electrical conductivity value at the input terminals of the antennas ($y=0$), and $l_d=40cm$ is the length of each dipole, assumed to have diameter$D_d=5mm$. In particular, ${\sigma }_0$has been determined by means of a dedicated parametric analysis. In this way, the optimal value${\sigma }_0={\sigma }_{opt}\simeq 200S/m$, resulting in a fractional bandwidth $FBW\simeq 51\%$ centred on the fundamental resonant frequency$f_r\simeq 450MHz$, has been found (see Fig. 3).

Figure 3.

Frequency behaviour of the individual antenna input reflection coefficient for different loading profiles. The antenna is elevated $h_d=3cm$ over the ground.

The FDTD characterization of the structure has been carried out by using a non-uniform computational grid with maximum cell size$\Delta h={\lambda }_{\min }/24=2.5mm$, where ${\lambda }_{\min }\simeq 6cm$ is the wavelength in the ground at the upper $$-10dB$$ cut-off frequency $f_{\max }=1GHz$ of the excitation voltage signal, which is the Gaussian voltage pulse described by the equation:

$$v_G(t)={\widehat{V}}_G\exp \left[-{\left(\frac{t-T_0}{T_G}\right)}^2\right]u(t),$$E18

where$${\widehat{V}}_G=1mV$$, $$T_0=10T_G$$, and:

$$T_G=\frac{\sqrt{\ln 10}}{\pi f_{\max }}\simeq 0.48ns.$$E19

In (18), $$u(t)$$denotes the usual Heaviside function. The source pulse is coupled into the finite-difference equations used to update the time-domain electric field distribution within the feed point of dipole $\#1$ as follows:

$${{ℰ}_y|}_{i_1,j_1+\frac{1}{2},k_1}^{n+1}={{\overline{\alpha }}_y^{(E,G)}|}_{i_1,j_1+\frac{1}{2},k_1}{{ℰ}_y|}_{i_1,j_1+\frac{1}{2},k_1}^{n-1}+{{\overline{\beta }}_y^{(E,G)}|}_{i_1,j_1+\frac{1}{2},k_1}\left[{{\left(\nabla \times \mathscr{H}\right)}_y|}_{i_1,j_1+\frac{1}{2},k_1}^{n+\frac{1}{2}}+{{ℐ}_G|}^{n+\frac{1}{2}}\right],$$E20

where:

$${{\overline{\alpha }}_y^{(E,G)}|}_{i_1,j_1+\frac{1}{2},k_1}=\frac{1-{{\overline{Q}}_y^{(E,G)}|}_{i_1,j_1+\frac{1}{2},k_1}}{1+{{\overline{Q}}_y^{(E,G)}|}_{i_1,j_1+\frac{1}{2},k_1}},$$E21

$${{\overline{\beta }}_y^{(E,G)}|}_{i_1,j_1+\frac{1}{2},k_1}=\frac{\Delta t/{{\overline{\epsilon }}_y^{eff}|}_{i_1,j_1+\frac{1}{2},k_1}}{1+{{\overline{Q}}_y^{(E,G)}|}_{i_1,j_1+\frac{1}{2},k_1}},$$E22

with:

$${{\overline{Q}}_y^{(E,G)}|}_{i_1,j_1+\frac{1}{2},k_1}=\frac{\Delta t}{2{{\overline{\epsilon }}_y^{eff}|}_{i_1,j_1+\frac{1}{2},k_1}{R}_G},$$E23

$$i_1$$, $$j_1$$, $$k_1$$being the spatial indices relevant to the source. In (20), the term:

$${{ℐ}_G|}^{n+\frac{1}{2}}={\frac{v_G(t)}{R_G}|}_{t=(n+1/2)\Delta t}$$E24

denotes the time-discretised nominal current delivered by the generator. Similarly, the electric field distribution within the feed point of dipole $\#2$ is updated by using the finite-difference equation:

$${{ℰ}_y|}_{i_2,j_2+\frac{1}{2},k_2}^{n+1}={{\overline{\alpha }}_y^{(E,L)}|}_{i_2,j_2+\frac{1}{2},k_2}{{ℰ}_y|}_{i_2,j_2+\frac{1}{2},k_2}^{n-1}+{{\overline{\beta }}_y^{(E,L)}|}_{i_2,j_2+\frac{1}{2},k_2}{{\left(\nabla \times \mathscr{H}\right)}_y|}_{i_2,j_2+\frac{1}{2},k_2}^{n+\frac{1}{2}},$$E25

where:

$${{\overline{\alpha }}_y^{(E,L)}|}_{i_2,j_2+\frac{1}{2},k_2}=\frac{1-{{\overline{Q}}_y^{(E,L)}|}_{i_2,j_2+\frac{1}{2},k_2}}{1+{{\overline{Q}}_y^{(E,L)}|}_{i_2,j_2+\frac{1}{2},k_2}},$$E26

$${{\overline{\beta }}_y^{(E,L)}|}_{i_2,j_2+\frac{1}{2},k_2}=\frac{\Delta t/{{\overline{\epsilon }}_y^{eff}|}_{i_2,j_2+\frac{1}{2},k_2}}{1+{{\overline{Q}}_y^{(E,L)}|}_{i_2,j_2+\frac{1}{2},k_2}},$$E27

$${{\overline{Q}}_y^{(E,L)}|}_{i_2,j_2+\frac{1}{2},k_2}=\frac{\Delta t}{2{{\overline{\epsilon }}_y^{eff}|}_{i_2,j_2+\frac{1}{2},k_2}{R}_L},$$E28

and with$$i_2$$, $$j_2$$, $$k_2$$denoting the spatial indices relevant to the load. The total voltage and current signals excited at the input terminals of the antenna pair, regarded as a two-port microwave network, are readily computed as:

$$v_m(t)=-{\displaystyle {\int }_{C_{V_m}}E(r,t)\cdot dl}\simeq -{{ℰ}_y(t)|}_{i_m,j_m+\frac{1}{2},k_m},$$E29

$$i_m(t)={\displaystyle {\oint }_{C_{I_m}}H\left(r,t\right)\cdot dl}\simeq {{\mathscr{H}}_x|}_{i_m,j_m+\frac{1}{2},k_m+\frac{1}{2}}^{n+\frac{1}{2}}-{{\mathscr{H}}_x|}_{i_m,j_m+\frac{1}{2},k_m-\frac{1}{2}}^{n+\frac{1}{2}}-{{\mathscr{H}}_z|}_{i_m+\frac{1}{2},j_m+\frac{1}{2},k_m}^{n+\frac{1}{2}}+{{\mathscr{H}}_z|}_{i_m-\frac{1}{2},j_m+\frac{1}{2},k_m}^{n+\frac{1}{2}},$$E30

where $C_{V_m}$ is an open contour extending along the delta gap, and $C_{I_m}$ a closed contour path wrapping around the driving point of dipole $\#m$ ($m=\mathrm{1,2}$). Under the mentioned assumptions, the normalized incident and reflected waves are evaluated as:

$$a_1(f)=\frac{1}{2}\left[\frac{V_1(f)}{\sqrt{R_0}}+I_1(f)\sqrt{R_0}\right],$$E31

$$b_1(f)=\frac{1}{2}\left[\frac{V_1(f)}{\sqrt{R_0}}-I_1(f)\sqrt{R_0}\right],$$E32

$$a_2(f)=0,$$E33

$$b_2(f)=\frac{V_2(f)}{\sqrt{R_0}},$$E34

where $R_0=R_G$ denotes the reference resistance and$$V_m(f)=F\left[v_m(t)\right]$$, $$I_m(f)=F\left[i_m(t)\right]$$, $$F[\cdot ]$$being the usual Fourier transform operator. Therefore, the scattering parameters of the structure can be easily determined as:

$$S_{11}(f)=\frac{b_1(f)}{a_1(f)},$$E35

$$S_{21}(f)=\frac{b_2(f)}{a_1(f)}.$$E36

As it appears from Fig. 4a, the return-loss level is slightly affected by the Tx–Rx antenna separation that, on the other hand, is primarily responsible for the parasitic coupling level between the radiating elements. The impact of the antenna elevation above the ground has been also analyzed (see Fig. 4b). It is worth noting that, as $$h_d$$ decreases, the fundamental resonant frequency of the dipole is shifted down because of the proximity effect of the soil. On the other hand, the ground influence on the $$S_{21}$$ parameter is remarkable only at high frequencies, where the coupling level between the two radiating elements tends to decrease as the dipoles approach the air-ground interface.

In the performed numerical computations, a ten-cell uniaxial perfectly matched layer (UPML) absorbing boundary condition for lossy media [19] has been used at the outer FDTD mesh boundary to simulate the extension of the space lattice to infinity. As outlined in [19], the UPML is indeed perfectly matched to the inhomogeneous medium formed by the upper air region and the lossy material modelling the soil. So, no spurious numerical reflections take place at the air-ground interface. In particular, a quartic polynomial grading of the UPML conductivity profile has been selected in order to have a nominal reflection error$$R_{PML}\simeq e^{-16}$$.

Figure 4.

Frequency behaviour of the scattering parameters of the dipole pair for different Tx–Rx separations (a) and antenna elevations (b) over the ground modelled as a lossy homogeneous half-space with electrical properties ${\epsilon }_r=6$ and$\sigma =15mS/m$.

4. The radar detection of buried pipes

In this section, emphasis is devoted to the analysis of the dipole antenna pair located above a lossy homogeneous/inhomogeneous material half space where an infinitely long pipe is buried (see Fig. 5). In such configuration, the transmit element of the radar unit emits an electromagnetic pulse that propagates into the ground, where it interacts with the target, modelled as a $y-$directed circular cylinder having diameter$D_p=30cm$, buried at a depth$h_p=40cm$. This interaction results in a diffracted electromagnetic field which is measured by the receive element of the radar. By changing the location of the radar on the soil interface and recording the output of the receive antenna as function of time (or frequency) and radar location, one obtains the scattering data, which can be processed to get an image of the subsurface.

Figure 5.

Top and cross-sectional view of a resistively loaded dipole antenna pair located above a lossy inhomogeneous ground, where an infinitely long pipe is buried. Structure characteristics:$l_d=40cm$, $D_d=5mm$, $\delta =2.5mm$, $s_d=20cm$, $h_d=3cm$, $D_p=30cm$,$h_p=40cm$.

The parasitic coupling level between transmit and receive antennas is a critical parameter in the design of ground-penetrating radars and satisfactory levels are usually achieved by empirical design methods. Anyway, the prediction of coupling levels already at the design stage enhances structure reliability, while also improving design cycle. To this end, the locally conformal FDTD model presented in Section 2 has been usefully adopted. In this way, as it can be noticed in Fig. 6, it has been found that the antenna return-loss parameter $$S_{11}$$ is not strongly affected by the buried target which, conversely, has a significant impact on the frequency behaviour of the coupling level between the radiating elements, due to the electromagnetic field interference processes occurring at the receiver end.

Figure 6.

Frequency behaviour of the scattering parameters of the dipole antenna pair for different electrical properties of the buried pipe. Structure characteristics:$l_d=40cm$, $D_d=5mm$, $\delta =2.5mm$, $s_d=20cm$, $h_d=3cm$, $D_p=30cm$,$h_p=40cm$.

4.1. Analysis of sub-surface scattering processes

The considered antenna pair has been used to analyze the sub-surface scattering phenomena arising from the field interaction with a PVC pipe buried at a depth$h_p=30cm$. As it can be easily inferred, the intensity of the radio signal diffracted by the pipe and measured at the input terminals of the receive antenna strongly depends on the electrical and geometrical properties of the target. In particular, the peak-to-peak level of the signal increases as the diameter, and hence the radar cross section of the pipe becomes larger (see Fig. 7). Another noticeable phenomenon is the sub-surface excitation of creeping waves. Such waves, propagating along the pipe surface, give rise to late-time pulse contributions in the received radio signal which can be clearly noticed in Fig. 7 as the radius of the object increases. Furthermore, it is worth noting that the strength of the creeping wave contribution tends to reduce with the pipe size because of the field attenuation relevant to the extra-path length.

Figure 7.

Time-domain behaviour of the radio signal diffracted by a PVC pipe buried at a depth $h_p=30cm$ in a lossy ground with electrical properties ${\epsilon }_r=6$ and$\sigma =15mS/m$. The effect of the sub-surface excitation of creeping waves can be noticed.

4.2. Impact of ground-embedded inhomogeneities

An invariable feature of real-life soils is heterogeneity. Without taking into account the inhomogeneities altering the idealized nature of the considered ground model, it becomes a futile effort to design a complex GPR system that will perform well over a real-life soil. To overcome this limitation, a realistic ground model has been developed by simulating small ellipsoidal scatterers embedded in the soil (see Fig. 5). The size, location and electrical properties of these inhomogeneities are determined randomly within preset limits. In particular, the maximum dimensions of the scatterers are$a_x=10cm$, $a_y=10cm$, and $a_z=2cm$ in$x$, $y$, and $z$ coordinate direction, respectively. In addition, all inhomogeneities have randomly selected values of the relative permittivity according to the following Gaussian probability distribution:

$${\phi }_{{\mu }_{\epsilon },{\sigma }_{\epsilon }^2}\left({\epsilon }_r\right)=\frac{1}{\sqrt{2\pi }{\sigma }_{\epsilon }}{e}^{-\frac{{\left({\epsilon }_r-{\mu }_{\epsilon }\right)}^2}{2{\sigma }_{\epsilon }^2}},$$E37

with mean${\mu }_{\epsilon }=13.7$, and standard deviation ${\sigma }_{\epsilon }=4.2$ (see Fig. 8). As it appears from Fig. 9, the ground-embedded inhomogeneities have a considerable impact on the coupling coefficient $S_{21}$ of the dipole pair especially at high frequencies ($f>0.6GHz$) where the radiated field assumes a quasi-optical behaviour and the diffraction processes arising from the interaction with the inhomogeneities tend to be significant, and could mask the detection of the buried target (see Fig. 10).

Figure 8.

Gaussian probability distribution of the relative permittivity of ground-embedded inhomogeneities.

Figure 9.

Effect of the ground inhomogeneities on the frequency behaviour of the scattering parameters of the dipole pair. Structure characteristics:$l_d=40cm$, $D_d=5mm$, $\delta =2.5mm$, $s_d=20cm$, $h_d=3cm$, $D_p=30cm$,$h_p=40cm$.

Figure 10.

Spatial distribution of the electric field excited along the vertical cut-plane of the antenna pair elevated over a lossy inhomogeneous ground, where an infinitely long metallic pipe is buried. Structure characteristics:$l_d=40cm$, $D_d=5mm$, $\delta =2.5mm$, $s_d=20cm$, $h_d=3cm$, $D_p=30cm$,$h_p=40cm$. Working frequency: $f=500MHz$(a), $f=1GHz$(b).

As outlined in [11], the rigorous analysis of the aforementioned subsurface scattering phenomena is very important in order to asses the suitability of detection and imaging algorithms for GPR applications in realistic scenarios.

5. Frequency-independent equivalent circuit model

Typically, electromagnetic field solvers and measurement systems, such as network analyzers, generate $S-$parameter representations of microwave components and antennas. However, circuit simulators, such as SPICE, require conventional equivalent circuits with lumped frequency-independent parameters that can be conveniently modelled within CAD software. In view of this, a technique for extracting the equivalent circuit from a given $S-$parameter representation is highly desirable [20].

The proposed procedure for the equivalent circuit extraction is based on a heuristic modification of the Cauer’s network synthesis technique [10]. Resistive elements are introduced to model metal, dielectric, and radiation losses. The scattering matrix $\underset{¯}{S}$ of the antenna pair is modelled by means of an equivalent network composed of shunt $\pi -$networks as shown in Fig. 11, where ideal transformers are to be included in order to assure the physical realizability of the network. In particular, RLC $\pi -$subnetworks are required to properly model the high-frequency resonant phenomena taking place into the structure, whereas the capacitive $\pi -$subnetwork is needed to describe the circuital behaviour of the antenna pair in the quasi-static regime. A simple analysis of the network shown in Fig. 11 allows determining the admittance matrix ${\underset{¯}{Y}}^{eq}$ of the equivalent circuit. Then, the corresponding scattering matrix ${\underset{¯}{S}}^{eq}$ can be easily computed as follows:

$${\underset{¯}{S}}^{eq}=\left(\underset{¯}{I}-{\underset{¯}{\widehat{Y}}}^{eq}\right)\cdot {\left(\underset{¯}{I}+{\underset{¯}{\widehat{Y}}}^{eq}\right)}^{-1},$$E38

where $\underset{¯}{I}$ is the identity matrix, and:

$${\underset{¯}{\widehat{Y}}}^{eq}=R_0{\underset{¯}{Y}}^{eq},$$E39

$R_0$denoting the selected reference resistance at the input ports of the transmit and receive antennas. The synthesis procedure is based on an iterative non-linear fitting procedure [8] aimed to minimize, across a specific frequency band$B$, the mean-square error functional:

$$e_B\left(R_{ij}^{(p)},L_{ij}^{(p)},C_{ij}^{(p)}\right)=\sqrt{{\displaystyle \underset{B}{\int }{\Vert {\underset{¯}{S}}^{FDTD}\left(f\right)-{\underset{¯}{S}}^{eq}\left(f,R_{ij}^{(p)},L_{ij}^{(p)},C_{ij}^{(p)}\right)\Vert }^{2}df}},$$E40

${\underset{¯}{S}}^{FDTD}$being the frequency-dependent scattering matrix of the structure computed by means of the locally conformal FDTD technique. In this way, it has been possible to evaluate numerically the circuital parameters, listed in Table 1, of the radar unit located above the lossy homogeneous material half space with relative permittivity ${\epsilon }_r=6$ and electrical conductivity$\sigma =15mS/m$. As shown in Fig. 12, a few resonant $\pi -$subnetworks are necessary to describe adequately the antenna return-loss and mutual coupling coefficient in the frequency band of interest.

Figure 11.

Frequency-independent equivalent circuit model. Quasi-static (a), and $p-$th high-frequency (b) $\pi -$subnetwork. For symmetrical structures the transformation ratio $n^{(p)}$ is equal to $\pm 1$ for each subnetwork.

Figure 12.

Magnitude of the scattering parameters of the dipole pair versus frequency. A good agreement between the FDTD numerical results and those obtained by means of the equivalent circuit model is observed. Structure characteristics:$l_d=40cm$, $D_d=5mm$, $\delta =2.5mm$, $s_d=20cm$,$h_d=3cm$.

Table 1.

Circuital parameters relevant to the equivalent circuit of the antenna pair. Structure characteristics:$l_d=40cm$, $D_d=5mm$, $\delta =2.5mm$, $s_d=20cm$,$h_d=3cm$.

6. Conclusion

The full-wave analysis of electromagnetic sensing of buried pipes with GPR in realistic scenarios has been carried out. An enhanced locally conformal FDTD technique, useful to accurately model complex electromagnetic structures as well as ground-embedded inhomogeneities with arbitrary shape and material parameters, has been adopted. By using this scheme, an extensive parametric analysis of the antenna scattering parameters and radiated near-field spatial distribution has been performed for different Tx–Rx antenna separations and elevations over the ground, taking into account the presence of buried metallic and dielectric targets, as well as soil-embedded ellipsoidal inhomogeneities with arbitrary size, location and electrical properties. The obtained numerical results provide a physical insight into the underlying mechanisms of subsurface scattering and antenna mutual coupling processes. Finally, a frequency-independent equivalent circuit, useful to be employed in CAD tools, has been derived from the antenna scattering parameters, showing that including the effect of just a few resonant modes yields high numerical accuracy.

7. Appendix

In order to validate the accuracy of the proposed locally conformal FDTD scheme a number of test cases have been considered. Here the results obtained for the computation of the fundamental resonant frequency of a dielectric resonator enclosed in a metallic cavity are presented. The structure under consideration (see Fig. 13a) has been already analyzed in [5]. It consists of a perfectly conducting metallic cavity of dimensions $a=b=50mm$ and$c=30mm$, loaded with a cylindrical dielectric (ceramic) puck having diameter$D=36mm$,

Figure 14.

Geometry of a dielectric loaded rectangular cavity (a), and behaviour of the relevant fundamental resonant frequency $f_r$ as function of the FDTD mesh size $\Delta h$ (b). Shown is the confidence region where the relative error $e_r$ with respect to the reference resonant frequency ${\nu }_r=1.625GHz$ [5] is smaller than$0.1\%$. Structure characteristics:$a=b=50mm$, $c=30mm$, $D=36mm$, $t=16mm$,$h=7mm$. Relative permittivity of the dielectric puck:${\epsilon }_r=37$.

height $t=16mm$ and relative dielectric constant${\epsilon }_r=37$. The puck is suspended at a distance of $h=7mm$ from the bottom of the cavity. Since the dielectric permittivity of the resonator is rather high, the effect of the orthogonal Cartesian mesh being not conform to the resonator shape is expected to be noticeable. Here the structure is analyzed by means of a standard FDTD scheme featuring the traditional staircase approximation of the resonator’s contour, and by means of the weighted averaging approach proposed in [7], and the locally conformal FDTD technique detailed in Section III. The numerical results obtained from these FDTD schemes are compared against the ones reported in [5] resulting from the use of a commercial Transmission Line Matrix (TLM) method-based solver. To this end, a cubic FDTD mesh having fixed spatial increment $\Delta h$ has been adopted to analyze the structure. As it appears in Fig. 13b, this example clearly demonstrates the suitability of the proposed approach to efficiently handle complex metal-dielectric structures with curved boundaries. The proposed locally FDTD scheme introduces a significant improvement in accuracy over the stair-casing approximation, converging very quickly to the reference value. Such feature is thus of crucial importance to optimize the design of antennas for ground-penetrating radar applications.