Bistatic Landmine and IED Detection Combining Vehicle and Drone Mounted GPR Sensors

Abstract

This work proposes a novel Ground Penetrating Radar (GPR) system to detect landmines and Improvised Explosive Devices (IEDs). The system, which was numerically evaluated, is composed of a transmitter placed on a vehicle and looking forward and a receiver mounted on a drone and looking downwards. This combination offers both a good penetration and a high resolution, enabling the detection of non-metallic targets and mitigating the clutter at the air–soil interface. First, a fast ray tracing simulator was developed to find proper configurations of the system. Then, these configurations were validated using a full wave simulator, considering a flat and a rough surface. All simulations were post-processed using a fast and accurate Synthetic Aperture Radar (SAR) algorithm that takes into account the constitutive parameters of the soil. The SAR images for all configurations were compared, concluding that the proposed contribution greatly improves the target detection and the surface clutter reduction over conventional forward-looking GPR systems.

1. Introduction

The non-invasive detection of hidden or buried objects has attracted an increasing interest due to its practical applicability in several fields such as civil engineering (structural and road inspection), security and defense (landmine detection), and archeology, among others [1]. These techniques are able to detect the concealed objects without physically interacting with them or the surrounding medium. Furthermore, they can even be used to image the inspected area. Electromagnetic induction, thermal imaging, nuclear quadrupole resonance or Ground Penetrating Radar (GPR) are some well-known examples of non-invasive techniques.

Among these techniques, GPR has been widely used for subsurface imaging applications [2]. It is based on transmitting an electromagnetic wave and detecting the scattered waves at the air–soil interface and from the buried targets, providing a radar image of the underground. One of its main advantages is that it can detect both metallic and dielectric targets. However, this technique is quite sensitive to the soil heterogeneity, the soil surface roughness and the possible low contrast between the soil and a non-metallic target [3]. As a result, it requires careful configuration and advanced signal processing techniques to overcome these issues and improve the detectability of the system.

GPR systems can be classified using different criteria. According to the distance between the antennas and the soil, they can be classified as ground-coupled or air-launched systems. The former usually allow a better penetration into the soil and are less affected by the reflections at the air–soil interface (provided the antennas are well-matched to the soil impedance). However, they need to be in contact with the soil, which also slows down the scanning speed and should be avoided when searching for dangerous targets such as landmines and Improvised Explosive Devices (IEDs). The latter avoid the interaction with the soil, but the strong obscuring clutter (due to impedance mismatch at the rough air–soil interface) greatly compromises the detection of buried targets. GPR systems can also be classified as Forward-Looking GPR (FLGPR) [4] and Down-Looking GPR (DLGPR) [5]. In vehicle mounted FLGPR systems, the antennas look ahead of a vehicle, with an angle of incidence that helps to maximize TM (Transverse Magnetic) waves penetration into the soil and/or to minimize reflections from the air–soil interface backscattered to the receiver. However, they have lower resolution (being difficult to distinguish whether the targets are over or under the surface) and sensitivity (since much of a flat-topped target’s reflections are in the forward opposite direction from the transmitter). Concerning DLGPR systems, the antennas are perpendicular to the soil surface, which yields higher resolution at the expense of stronger clutter.

In landmine detection, the scanning system must keep a safety distance from the inspected area in order to avoid the threat of explosion, which thus strongly favors FLGPR. To address this issue, a GPR system on board an Unmanned Aerial Vehicle (UAV) or a drone has been recently presented for subsurface imaging [6,7]. This system provides high resolution subsurface images since it allows the coherent combination of measurements using a Synthetic Aperture Radar (SAR) algorithm. However, the strong clutter at the air–soil interface clearly degrades the detection capabilities, especially when the contrast between the target and the soil is low. It must also be noticed that, although bistatic SAR systems have gained an increasing interest in the last years [8], most GPR systems (both DLGPR and FLGPR) adopt a monostatic or a quasi-monostatic configuration.

It would be desirable to combine the advantages of FLGPR and DLGPR systems in order to obtain both good penetration into the soil and high resolution. This article is devoted to analyzing this novel GPR configuration. As shown in Figure 1, a transmitter is placed on a vehicle with the antenna looking ahead and a receiver is placed on a UAV with the antenna pointing straight down to the soil surface. A fast ray-tracing method has been developed to find feasible configurations of the system. Then, the resulting configurations were accurately analyzed with a Finite-Difference Frequency-Domain (FDFD) method. These configurations as well as the multimonostatic configuration (corresponding to just a DLGPR system on board a moving UAV) were compared by post-processing the simulated scattered field with a SAR algorithm.

2. Methodology

2.1. Scenario

Ray-tracing and FDFD methods are used to simulate a 2D GPR scenario, such as the one shown in Figure 2, with a target buried in the soil. In this scenario, there are T transmitters (TX) placed at positions ${\mathit{r}}_t$ ($t=1,\dots ,T$) and, for each transmitter, there are R receivers (RX) at positions ${\mathit{r}}_r^t$ (where subindex $r=1,\dots ,R$ denotes the receiver and superindex t the transmitter). The soil is characterized by its relative permittivity ${{ϵ}_r}_s$ and its conductivity ${\sigma }_s$. The target is assumed to have a circular shape, with radius ${\delta }_{tg}$ and centered at coordinates ($x_{tg},y_{tg}$). It is characterized by its constitutive parameters ${{ϵ}_r}_{tg}$ and ${\sigma }_{tg}$. The simulation is performed at N frequencies, assuming either TE (Transverse Electric) or TM (Transverse Magnetic) polarization.

2.2. Ray-Tracing

Ray-tracing (RT) is a geometrical optics method that models propagation by following straight rays [9]. Although it is less accurate than conventional full-wave methods, it requires a much lower computational effort. Thus, it is useful for fast modeling of large scenarios at several frequencies and for several transmitter and receiver positions.

2.2.1. Field Computation

The contribution to the electric field of a ray impinging a given receiver in air is calculated according to Equation (1) or Equation (2), depending on whether the ray comes from the reflection at the soil interface or from the target. In these equations $E_{\mathrm{inc}}$ is the incident field amplitude; $A_t$ and $A_r$ are used to take into account the transmitter and receiver antenna beamwidth; $G_{\mathrm{in}}$ and $G_{\mathrm{out}}$ are the in-plane and out-of-plane geometrical spreading factors [10]; $\mathrm{\Gamma }$ and $\tau$ denote the reflection and transmission coefficients; $R_m$, ${\alpha }_m$ and ${\beta }_m$ are the total ray path-length and the attenuation and phase constants in medium m (where $m=0$ is air and $m=s$ is soil). The ray-tracing method implemented in this contribution calculates the path length in each medium ($R_m$), which is then multiplied by the propagation constants so as to perform a multifrequency simulation.

$$E_{\mathrm{soil}}=\frac{E_{\mathrm{inc}}{A}_t{A}_r}{G_{\mathrm{in}}{G}_{\mathrm{out}}}{\mathrm{\Gamma }}_{\mathrm{air}‒\mathrm{soil}}exp(-j{\beta }_0{R}_0)$$

(1)

$$\begin{array}{cc}\hfill E_{\mathrm{target}}=& \frac{E_{\mathrm{inc}}{A}_t{A}_r}{G_{\mathrm{in}}{G}_{\mathrm{out}}}{\tau }_{\mathrm{air}‒\mathrm{soil}}{\mathrm{\Gamma }}_{\mathrm{soil}‒\mathrm{target}}{\tau }_{\mathrm{soil}‒\mathrm{air}}exp(-{\alpha }_s{R}_s)exp(-j({\beta }_0{R}_0+{\beta }_s{R}_s))\hfill \end{array}$$

(2)

For these ray-tracing simulations, each transmitter and receiver is characterized by its angle of incidence with respect to the soil (${{\theta }_i}_t$ and ${{\theta }_i}_r^t$, respectively) and by its 3-dB beamwidth. The terms $A_t$ and $A_r$ model the antenna assuming a $co{s}^q$ pattern (where q is calculated according to the antenna beamwidth).

Assuming a moderately lossy soil, multiple reflections are not considered and the reflection and transmission angles are calculated using Snell’s law without taking into account the conductivity of the soil. These angles are then used to compute the reflection and transmission coefficients ($\mathrm{\Gamma }$ and $\tau$), which do incorporate the conductivity of the soil.

2.2.2. Implementation

Usually, many rays are launched from each transmitter for proper illumination of the scenario [11]. However, to reduce the computational time required for multiple rays, only the rays that come from the specular reflection at the air–soil interface and the rays coming from the target are used. The former can be calculated directly using simple geometrical relations. For the latter, we estimate the angles of the incident rays that hit the left and the right sides of the target (i.e., the points $(x_{tg}-{\delta }_{tg},y_{tg})$ and $(x_{tg}+{\delta }_{tg},y_{tg})$). These estimations require first calculating the refraction points at the air–soil interface for those points on the target. Then, many rays are launched between the computed angles. If one of these rays (after its reflection at the target) is closer than a given threshold ($th$) to a receiver, it is assumed that ray hits that receiver and it is used to compute $E_{target}$. A detailed flowchart of the implemented approach is shown in Figure 3.

2.3. FDFD

The 2D Finite Difference Frequency Domain (FDFD) algorithm is well suited to nearfield analysis of dielectric or metal targets from about 0.1 to 30 wavelengths in size placed in lossy, rough dielectric backgrounds [12]. This algorithm simulates nearfield scattering from objects that are of electrical sizes that are particularly difficult to model (less than about 30 wavelengths), filling a desirable niche between geometric optics methods (high frequency or electrically large scatterers) [13] and Born approximation methods (low frequency or electrically small scatterers) [14]. The scattering objects and backgrounds can both be lossy dielectrics of any contrast and any loss tangent. The 2DFDFD algorithm subdivides space into uniform Yee cells and applies simple finite differences to describe the 2D partial differential Helmholtz wave equation for which one dimension, typically z, is invariant. Termination of the space is done by using a perfectly matched layer (PML) to minimize scattering from the computational boundaries [15]. Compared with 3DFDFD methods which require iterative (and slow) GMRES (Generalized Minimal Residual Method) or LGMRES (“loose” GMRES) solvers, the simpler 2DFDFD algorithm uses direct matrix inversion for both the TM and TE subclasses of problems. Computational time is relatively fast and complex geometries are easy to model.

2.4. Inversion

To compare the results for the different configurations, the simulated field is represented in the time-domain (B-scan) and post-processed with a SAR algorithm. SAR reflectivity at point ${\mathit{r}}^{\prime }$ of the investigation domain is given by Equation (3), where $R_r^t$ is the path length between the t-th transmitter (located at ${\mathit{r}}_t$), the point where the reflectivity is calculated ${\mathit{r}}^{\prime }$ and the r-th receiver ${\mathit{r}}_r^t$.

$$\rho \left({\mathit{r}}^{\prime }\right)={\displaystyle \sum }_{n=1}^N{\displaystyle \sum }_{t=1}^T{\displaystyle \sum }_{r=1}^{R}E(f_n,{\mathit{r}}_t,{\mathit{r}}_r^t)exp(+j{\beta }_0{R}_r^t)$$

(3)

Assuming free-space propagation, $R_r^t$ would be equal to $\parallel {\mathit{r}}_t-{\mathit{r}}^{\prime }\parallel +\parallel {\mathit{r}}_r^t-{\mathit{r}}^{\prime }\parallel$. This assumption provides good results when the incidence angle is close to normal incidence and the permittivity and conductivity of the soil are low. When these conditions are fulfilled, it is possible to detect the object in the SAR image at approximately $\sqrt{{{ϵ}_r}_s}d$ depth (being d the true depth of the buried target). To obtain better results and to detect the object at its real depth, the constitutive parameters of the soil must be taken into account. The common approach consists of calculating the refraction point at the air–soil interface (for each point ${\mathit{r}}^{\prime }$ in the investigation domain, and each combination of transmitter and receiver positions). This requires solving a fourth-order equation derived for Snell’s Law. However, instead of calculating the refraction point, $R_r^t$ is modified so as to consider the permittivity of the soil [16,17]. Thus, $R_r^t$ is given by Equation (4), where $n_s=\sqrt{{{\epsilon }_r}_s-1}-\sqrt{{{\epsilon }_r}_s}$ and the other parameters are defined according to the scheme shown in Figure 4.

$$\begin{array}{c}\hfill R_r^t=2d\sqrt{{{\epsilon }_r}_s-1}+\frac{d_t(d_t-d{n}_{s}cos\left(2{\varphi }_t\right))}{d_t+d{n}_s{sin}^2\left(2{\varphi }_t\right)}+\frac{d_r(d_r-d{n}_{s}cos\left(2{\varphi }_r\right))}{d_r+d{n}_s{sin}^2\left(2{\varphi }_r\right)}\end{array}$$

(4)

3. Results

3.1. Scenario Configuration

The scenario simulated with these methods consists of a low moisture sandy soil (with ${ϵ}_{r_s}=2.5$ and ${\sigma }_s=0.0125$ S/m) where a target of 2 cm radius is buried at 25 cm depth ($x_{tg}=0$ m and $y_{tg}=-0.25$ m). Both metallic and dielectric targets are considered. If the target is dielectric, it is modeled as trinitrotoluene (TNT) with ${ϵ}_{r_{tg}}=2.9$ and ${\sigma }_{tg}=0$ S/m. Simulations are performed between $3.5$ and $5.5$ GHz at 10 MHz steps considering TE polarization. Although these frequencies are higher than those commonly used in GPR, they have been chosen so that the radar could be light enough to be mounted on board a UAV (as it has been already proved in the prototype shown in [6]). In the ray-tracing simulation, the antenna beamwidth is ${30}^{\circ }$ and the results are contaminated with white Gaussian noise, resulting in a signal to noise ratio of 30 dB. It must be noted that the direct signal between the TX and RX has not been included in the simulations, since it is expected to be removed from the received signal in a real implementation of the system (thanks to the fact that it will arrive earlier than the signal coming from the soil reflection and it will likely to be stronger).

Three different configurations have been simulated:

• Multimonostatic, where the TX–RX (drone mounted transceiver) is placed at 65 different positions between down-track positions $x=-0.8$ m and $x=0.8$ m at $y=1$ m height. The angle of incidence is $0^{\circ }$ (i.e., the antennas are aligned perpendicular to the soil surface, with main beam pointing straight down).
• Multistatic, where the TX is placed at a fixed position (on a vehicle, at down-track position $x=-20$ m and height $y=2.5$ m) with main beam pointing at an angle of incidence of ${83}^{\circ }$ with the nominal ground surface, and the drone-mounted RX is looking downward and is moved to the same positions as in the multimonostatic case.
• Multibistatic, where the vehicle-mounted TX is placed at $y=2.5$ m height and is moved between down-track positions $x=-20.8$ m and $x=-19.2$ m and the drone-mounted RX is moved between the same positions as in the multimonostatic case. Thus, both TX and RX are moved coherently. The angles of incidence are ${83}^{\circ }$ for the TX and $0^{\circ }$ for the RX.

The positions of the TX–RX in the multimonostatic configuration were set according to those already used in previous experimental work. The positions of the TX–RX and the angle of incidence in the multistatic and multibistatic configurations were found using ray tracing simulations, so as to be able to detect dielectric targets. The performance of all configurations were then verified with FDFD.

3.2. Initial Comparison: Scattered Field and B-Scan

Before applying the inversion algorithm, the simulated scattered fields obtained with each method were compared, in both the frequency and time domains. This comparison is shown in Figure 5, Figure 6 and Figure 7 for the multimonostatic scenario with a metallic target buried in the soil. The normalized scattered field in the frequency domain is shown for two observation domain positions: $x=-0.8$ m (Figure 5) and $x=0$ m (Figure 6). The inverse Fourier Transform is used to compute the scattered field in the time domain (B-scan), as shown in Figure 7. The agreement between the scattered field simulated with RT and TE FDFD modeling $E_z$ (or TM${}_z$, relative to z-axis) is good. The main difference is that in FDFD the amplitude at the air–soil interface is larger than the amplitude at the target, whereas in RT both amplitudes are similar. This might be due to the fact that RT only considers the specular reflection at the air–soil interface. This fact also explains that the scattered fields in the frequency domain are more similar at $x=-0.8$ m (left side of the observation domain) than at $x=0$ m (center of the observation domain, exactly over the target).

3.3. SAR Image Comparison

The final goal was to compare the SAR images for each configuration (multimonostatic, multistatic and multibistatic) to determine the best configuration. First, the SAR images were obtained from the RT simulations and then the results were verified with the FDFD simulations.

3.3.1. Multimonostatic Simulations

The SAR image of the multimonostatic scenario with a buried metallic target is shown in Figure 8. Both the interface and the object are clearly detected. The reflectivity at the interface is larger in the FDFD simulation, as expected from the previous discussion.

When the buried target is dielectric (TNT), it is hardly detected in the SAR image (as shown in Figure 9). Thus, in accordance with the initial hypothesis, the antenna configuration must be improved to be able to detect non-metallic targets.

3.3.2. Multistatic and Multibistatic Simulations

Since the goal is to detect non-metallic targets, the multistatic and multibistatic simulations comparison was performed when the buried target is dielectric. SAR images are shown in Figure 10 and Figure 11 for the multistatic and multibistatic scenarios, respectively. In RT simulations, the specular reflections from the soil surface do not reach the receiver. Therefore, in FDFD simulations, the known flat ground background is removed, thus showing only the target-scattered response. The results are almost the same for both configurations, where the dielectric object is clearly distinguishable. There is also a good agreement between the RT and FDFD simulations, thus it can be concluded that RT is a useful tool for designing new GPR configurations.

3.3.3. Effect of Inversion Path Length

As aforementioned, the permittivity of the soil must be known or estimated in order to obtain an accurate SAR image. If the permittivity is not taken into account, the object is detected deeper and bigger than expected, as shown for the multimonostatic case in Figure 12a. Furthermore, if the permittivity or the angle of incidence is high, the resulting SAR image is considerably distorted, as shown for the multibistatic scenario in Figure 12b.

3.4. Computational Performance

The computational time required by the forward wave modeling (RT and FDFD) for each scenario is shown in

Table 1. Comparison of computational time.

Method	Multimonostatic	Multistatic	Multibistatic
RT	6.1 s	0.5 s	9.2 s
FDFD	28.5 h	32.2 min	32.7 h

. It was measured running the codes in a conventional laptop with 16 GB of RAM and Intel Core i7-6700HQ. As expected, RT is much faster than FDFD. It must be noticed that the pathlength is computed only once for each ray, since it is independent of the frequency. Furthermore, RT can be parallelized, which would reduce the computational time even more.

Regarding the inversion step, the SAR algorithm takes $5.1$ s to run. It must be notice that it can also be easily parallelized to reduce the computational time.

3.5. Effect of Rough Surface

One of the main difficulties for detecting buried objects with GPR is the strong random clutter produced by the reflection at the rough air–soil interface (due to impedance mismatch). If the interface is rough, this reflection cannot be easily removed and it clearly worsens the detection capability of the system. Therefore, it was necessary to analyze the performance of the different configurations with a rough surface. This analysis was performed with FDFD, considering a rough surface with an average height of 2 cm and a correlation length of 60 cm.

The SAR images for the multimonostatic and multibistatic configurations with a buried dielectric target are shown in Figure 13. In the multimonostatic configuration, the target is hardly detected (similarly as for the flat surface). In the multibistatic configuration, the target is still detected, although there are some stronger reflections at the air–soil interface, which worsens the quality of the SAR image. The inversion results for the multistatic configuration are similar to the multibistatic arrangement. The target detection is improved relative to the mutimonostatic configuration because, in the multistatic and multibistatic configurations, rays specularly reflected from the ground surface do not reach the receiver. Using a forward-stationed drone receiver also offers advantages of higher signal strength and closer angular proximity to the maximum target-scattered response, relative to a vehicle-mounted receiver [16]. Thus, it can be concluded that the proposed multibistatic configuration allows to detect dielectric targets even under rough surfaces.

3.6. Effect of Polarization

All previous simulations were performed considering TE polarization. Since the reflection coefficient is smaller for TM polarization (which provides better penetration into the soil), an improvement in the results could be expected with TM. However, when the surface is rough, there is not only better penetration, but also slightly higher clutter levels with TM. As a result, the contrast between the target and the soil remains almost the same, as shown in Figure 14 for the multibistatic configuration. Nonetheless, there is less clutter below the soil surface, which helps to facilitate the detection.

4. Analysis

As mentioned above, the target detection and surface clutter reduction in the proposed distributed GPR configuration are better than conventional GPR architectures. The previous analysis was performed for a small generic target (both dielectric and metallic), whose behavior is similar to a point source.

To further analyze this system for the proposed application (landmine and IED detection), a target with characteristics similar to a plastic PMN-type landmine was considered. In particular, the target, which is fully composed by TNT, has a cross-section of 10 cm × 4 cm and it has been buried at 25 cm depth under the same rough surface used previously. This is one of the most challenging scenarios for this application mainly due to the major part of the reflected signal from the target being specularly-reflected away from the receiver, the low-contrast between the target and the soil, and the clutter due to the soil surface roughness.

The SAR images obtained from FDFD simulations for the configurations previously compared (multimonostatic DLGPR, multistatic and multibistatic FLGPR–DLGPR) as well as the multimonostatic FLGPR configuration were compared in this scenario. In the multimonostatic FLGPR configuration (shown in this section for comparison purposes), the TX–RX is moved between down-track positions $x=-20.8$ m and $x=-19.2$ m at $y=2.5$ m height, with main beam pointing at an angle of ${83}^{\circ }$. The results for all these configurations are shown in Figure 15. As explained above, the multistatic and multibistatic FLGPR–DLGPR configurations show similar results (Figure 15c,d, respectively). In the multimonostatic configurations (Figure 15a,b), the target cannot be detected at all. Furthermore, in the multimonostatic FLGPR case (Figure 15b), it can also be seen that there is poor resolution in vertical position, which implies that it is not possible to distinguish even whether the target is above or below the surface (in case it was detected) [16]. For this configuration, the bright high-intensity bands shown in the SAR image are generated entirely from the rough surface scattering, and their extended lengths are due to poor synthetic focusing in the down range direction. Their change of direction is due to the fact that the refraction at the soil surface is taken into account in the inversion. However, in the proposed multibistatic FLGPR–DLGPR configuration (Figure 15d), two reflections coming from the target (at $-15$ dB mainly from its left and right sides) can be clearly distinguished, thus allowing the detection of the target.

5. Conclusions

A novel GPR system architecture was designed using a fast ray tracing algorithm. The proposed architecture was validated using FDFD, even taking into account the surface roughness. This new system exploits the advantages of FLGPR (with the transmitter looking ahead of a vehicle) and DLGPR (with the receiver on a drone), mainly in terms of less clutter and good penetration and resolution. Compared to a multimonostic FLGPR (on a vehicle at a stand-off distance), it provides higher signal strength and resolution, since the receiver is closer to the investigation domain. Compared to a multimonostatic DLGPR (on board a UAV), it helps to mitigate the strong clutter from the air–soil interface. For the analyzed aperture size and TX–RX distance, it was also shown that the results with the transmitter stationary (multistatic) or moving synchronously with the receiver (multibistatic) are almost the same.