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 1 GHz[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 Dfilled by a linear, isotropic, non dispersive material, having permittivityε(r), magnetic permeabilityμ(r), and electrical conductivityσ(r). In such a domain, a dual-space, non-uniform lattice formed by a primary and secondary mesh is introduced. The primary mesh MDis composed of space-filling hexahedrons, whose vertices are defined by the Cartesian coordinates:

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

The secondary or dual mesh 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 coordinatesxi+1/2=xi+Δ xi/2, yj+1/2=yj+Δ yj/2,zk+1/2=zk+Δ zk/2. A set of dual edge lengths is then introduced in M˜Das follows:

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 S˜x|i+1/2,j,khaving boundary ∂S˜x|i+1/2,j,k=C˜x|i+1/2,j,k(see Fig. 1) results in the following integral equation:

where:

as ℓy jand ℓz ktend to zero. Under the assumption that the spatial incrementsΔ xi, Δ yj, Δ zkof 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 S˜x|i+1/2,j,kso that, under the mentioned hypothesis, the following approximation can be made:

Hence, combining the equations above yields:

where we have introduced the normalized field quantities:

and the averaged effective permittivityε¯xeff|i+1/2,j,k, and conductivityσ¯xeff|i+1/2,j,k, defined as follows:

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:

where:

denotes the finite-difference expression of the normalized x−component of the magnetic field curl at the time stepn+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:

with:

The update equations of the remaining components of the electric and magnetic field can be easily derived by permuting the spatial indicesi, j, kand 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 εr=6and electrical conductivityσ=15 mS/m. The geometry of the structure is depicted in Fig. 2.

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

has been applied to the flairs of the considered radiators. In (17), σ0denotes the electrical conductivity value at the input terminals of the antennas (y=0), and ld=40 cmis the length of each dipole, assumed to have diameterDd=5 mm. In particular, σ0has been determined by means of a dedicated parametric analysis. In this way, the optimal valueσ0=σopt≃200 S/m, resulting in a fractional bandwidth FBW≃51%centred on the fundamental resonant frequencyfr≃450 MHz, has been found (see Fig. 3).

The FDTD characterization of the structure has been carried out by using a non-uniform computational grid with maximum cell sizeΔh=λmin/24=2.5 mm, where λmin≃6 cmis the wavelength in the ground at the upper −10dBcut-off frequency fmax=1 GHzof the excitation voltage signal, which is the Gaussian voltage pulse described by the equation:

whereV^G=1 mV, T0=10 TG, and:

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 #1as follows:

where:

with:

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

where:

and withi2, j2, k2denoting 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:

where CVmis an open contour extending along the delta gap, and CIma closed contour path wrapping around the driving point of dipole #m(m=1,2). Under the mentioned assumptions, the normalized incident and reflected waves are evaluated as:

where R0=RGdenotes the reference resistance andVm(f)=F[vm(t)], Im(f)=F[im(t)], F[ ⋅ ]being the usual Fourier transform operator. Therefore, the scattering parameters of the structure can be easily determined as:

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 hddecreases, 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 S21parameter 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 errorRPML≃e−16.

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 diameterDp=30 cm, buried at a depthhp=40 cm. 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.

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 S11is 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.

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 areax=10 cm, ay=10 cm, and az=2 cminx, y, and zcoordinate direction, respectively. In addition, all inhomogeneities have randomly selected values of the relative permittivity according to the following Gaussian probability distribution:

φμε,σε2(εr)=12π σεe−(εr−με)22σε2 ,E37

with meanμε=13.7, and standard deviation σε=4.2(see Fig. 8). As it appears from Fig. 9, the ground-embedded inhomogeneities have a considerable impact on the coupling coefficient S21of the dipole pair especially at high frequencies (f>0.6 GHz) 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).

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 S¯of the antenna pair is modelled by means of an equivalent network composed of shunt π−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 π−subnetworks are required to properly model the high-frequency resonant phenomena taking place into the structure, whereas the capacitive π−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 Y¯eqof the equivalent circuit. Then, the corresponding scattering matrix S¯eqcan be easily computed as follows:

where I¯is the identity matrix, and:

R0denoting 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 bandB, the mean-square error functional:

S¯FDTDbeing 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 εr=6and electrical conductivityσ=15 mS/m. As shown in Fig. 12, a few resonant π−subnetworks are necessary to describe adequately the antenna return-loss and mutual coupling coefficient in the frequency band of interest.

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=50 mmandc=30 mm, loaded with a cylindrical dielectric (ceramic) puck having diameterD=36 mm,

height t=16 mmand relative dielectric constantεr=37. The puck is suspended at a distance of h=7 mmfrom 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 Δhhas 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.