Bistatic landmine and IED detection combining vehicle and drone mounted GPR sensors

remote sensing Letter Bistatic Landmine and IED Detection Combining Vehicle and Drone Mounted GPR Sensors Maria Garcia-Fernandez 1,* , Ann Morgenthaler 2, Yuri Alvarez-Lopez 1, Fernando Las Heras 1 and Carey Rappaport 2 1Group of Signal Theory and Communications, University of Oviedo, 33203 Gijón, Spain; [email protected] (Y.A.-L.); [email protected] (F.L.H.) 2Department of Electrical and Computer Engineering, Northeastern University, Boston, MA 02115, USA; [email protected] (A.M.); [email protected] (C.R.) *Correspondence: gar[email protected] Received: 31 July 2019; Accepted: 29 September 2019; Published: 2 October 2019   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. Keywords: landmine detection; Improvised Explosive Device (IED); Ground Penetrating Radar (GPR); drone; bistatic radar 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 Remote Sens. 2019,11, 2299; doi:10.3390/rs11192299 www.mdpi.com/journal/remotesensing Remote Sens. 2019,11, 2299 2 of 14 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. TX RX Target Figure 1. Scheme of the novel GPR system combining FLGPR and DLGPR. 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 rt ( t= 1, ..., T ) and, for each transmitter, there are R receivers (RX) at positions rt r (where Remote Sens. 2019,11, 2299 3 of 14 subindex r= 1, ..., R denotes the receiver and superindex t the transmitter). The soil is characterized by its relative permittivity ers and its conductivity σs . The target is assumed to have a circular shape, with radius δtg and centered at coordinates ( xtg , ytg ). It is characterized by its constitutive parameters ertg and σtg . The simulation is performed at N frequencies, assuming either TE (Transverse Electric) or TM (Transverse Magnetic) polarization. ... ... ... ... ... R rx T tx Target Soil Figure 2. 2D GPR modeling scenario. 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 Einc is the incident field amplitude; At and Ar are used to take into account the transmitter and receiver antenna beamwidth; Gin and Gout are the in-plane and out-of-plane geometrical spreading factors [ 10 ]; Γ and τ denote the reflection and transmission coefficients; Rm , αm and β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 ( Rm ), which is then multiplied by the propagation constants so as to perform a multifrequency simulation. Esoil =EincAtAr GinGout Γair–soil exp(−jβ0R0)(1) Etarget =Einc AtAr GinGout τair–soilΓsoil-targetτsoil-air exp(−αsRs)exp(−j(β0R0+βsRs)) (2) For these ray-tracing simulations, each transmitter and receiver is characterized by its angle of incidence with respect to the soil ( θit and θit r , respectively) and by its 3-dB beamwidth. The terms At and Ar model the antenna assuming a cosq pattern (where q is calculated according to the antenna beamwidth). Remote Sens. 2019,11, 2299 4 of 14 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 ( Γ and τ ), 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 (xtg −δtg , ytg) and (xtg +δtg , ytg) ). 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 Etarget . A detailed flowchart of the implemented approach is shown in Figure 3. Calculate specular ray at air-soil interface Estimate refraction points at the air-soil interface for points (xtg-δtg, ytg) and (xtg+δtg, ytg) Launch Nr rays between these points Find intersection with soil surface (i.e. refraction point) Calculate transmitted ray (from air to soil) Intersection with target? Calculate reflected ray (from target to soil) Intersection with soil surface? Calculate transmitted ray (from soil to air) Calculate distance from the closest ray dr dr < th For each ray For each pair TX-RX For each RX For each TX This ray hits the RX Compute Etarget Compute Esoil Yes Yes Yes Figure 3. Flowchart of the ray-tracing implementation. Remote Sens. 2019,11, 2299 5 of 14 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 r0 of the investigation domain is given by Equation (3), where Rt r is the path length between the t-th transmitter (located at rt), the point where the reflectivity is calculated r0and the r-th receiver rt r. ρ(r0) = N ∑ n=1 T ∑ t=1 R ∑ r=1 E(fn,rt,rt r)exp(+jβ0Rt r)(3) Assuming free-space propagation, Rt r would be equal to krt−r0k+krt r−r0k . 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 √ersd 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 r0 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, Rt r is modified so as to consider the permittivity of the soil [ 16 , 17 ]. Thus, Rt r is given by Equation (4), where ns=√εrs−1−√εrs and the other parameters are defined according to the scheme shown in Figure 4. Rt r=2dpεrs−1+dt(dt−dnscos(2φt)) dt+dnssin2(2φt)+dr(dr−dnscos(2φr)) dr+dnssin2(2φr)(4) Remote Sens. 2019,11, 2299 6 of 14 Figure 4. Scheme for estimating the path length (the green dashed line represents the true ray path). 3. Results 3.1. Scenario Configuration The scenario simulated with these methods consists of a low moisture sandy soil (with ers= 2.5 and σs= 0.0125 S/m) where a target of 2 cm radius is buried at 25 cm depth ( xtg = 0 m and ytg =− 0.25 m). Both metallic and dielectric targets are considered. If the target is dielectric, it is modeled as trinitrotoluene (TNT) with ertg = 2.9 and σ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 ◦ 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 ◦ (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 ◦ 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◦for the TX and 0◦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. Remote Sens. 2019,11, 2299 7 of 14 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 Figures 5–7for 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 Ez (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). Frequency [GHz] 3.5 4 4.5 5 5.5 Real Part -1 -0.5 0 0.5 1 RT FDFD (a) Frequency [GHz] 3.5 4 4.5 5 5.5 Imaginary Part -1 -0.5 0 0.5 1 RT FDFD (b) Figure 5. Normalized scattered field at the first transmitter-receiver position ( x=− 0.8 m): multimonostatic scenario with a metallic target. Real part (a) and imaginary part (b). Frequency [GHz] 3.5 4 4.5 5 5.5 Real Part -1 -0.5 0 0.5 1 RT FDFD (a) Frequency [GHz] 3.5 4 4.5 5 5.5 Imaginary Part -1 -0.5 0 0.5 1 RT FDFD (b) Figure 6. Normalized scattered field in the middle of the observation domain ( x= 0 m): multimonostatic scenario with a metallic target. Real part (a) and imaginary part (b). Remote Sens. 2019,11, 2299 8 of 14 (a) (b) Figure 7. B-Scan comparison from RT (a) and FDFD (b) simulations: multimonostatic scenario with a metallic 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. (a) (b) Figure 8. SAR image from RT ( a ) and FDFD ( b ) simulations: multimonostatic scenario with a metallic target. 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. Remote Sens. 2019,11, 2299 9 of 14 (a) (b) Figure 9. SAR image from RT ( a ) and FDFD ( b ) simulations: multimonostatic scenario with a dielectric target. 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 Figures 10 and 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. (a) (b) Figure 10. SAR image from RT ( a ) and FDFD ( b ) simulations: multistatic scenario with a dielectric target.