Next Article in Journal
Three-Dimensional Structure and Transport Properties of Dust Aerosols in Central Asia—New Insights from CALIOP Observations, 2007–2022
Next Article in Special Issue
Range Limitations in Microwave Quantum Radar
Previous Article in Journal
Hyperspectral Image Classification Based on Double-Branch Multi-Scale Dual-Attention Network
Previous Article in Special Issue
Approximate Evaluation of the Resolution in Near Field Remote Sensing
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Quantitative Inversion of Multiantenna Ground-Penetrating Radar Data with Modeling Error Correction Based on Long Short-Term Memory Cells

Department of Electrical, Electronic, Telecommunications Engineering and Naval Architecture (DITEN), University of Genoa, 16145 Genoa, Italy
*
Author to whom correspondence should be addressed.
Remote Sens. 2024, 16(12), 2050; https://doi.org/10.3390/rs16122050
Submission received: 8 May 2024 / Revised: 31 May 2024 / Accepted: 4 June 2024 / Published: 7 June 2024
(This article belongs to the Special Issue Microwave Tomography: Advancements and Applications)

Abstract

:
Quantitative inversion of GPR data opens the door to precise characterization of underground environments. However, in order to make the inverse scattering problem solution easier from a computational viewpoint, simplifying assumptions are often applied, i.e., two-dimensional approximations or the consideration of idealized field probes and electromagnetic sources. These assumptions usually produce modeling errors, which can degrade the dielectric reconstruction results considerably. In this article, a processing step based on long short-term memory cells is proposed for the first time to correct the modeling error in a multiantenna GPR setting. In particular, time-domain GPR data are fed into a neural network trained with couples of finite-difference time-domain simulations, where a set of sample targets are simulated in both realistic and idealized configurations. Once trained, the neural network outputs an approximation of multiantenna GPR data as they are collected by an ideal two-dimensional measurement setup. The inversion of the processed data is then accomplished by means of a regularizing Newton-based nonlinear scheme with variable exponent Lebesgue space formulation. A numerical study has been conducted to assess the capabilities of the proposed inversion methodology. The results indicate the possibility of effectively compensating for modeling error in the considered test cases.

1. Introduction

Ground-penetrating radar (GPR) has been revealed to be a very useful tool in many different fields, ranging from civil and environmental engineering [1] to geophysics [2], and archeology [3] to forensic science [4]. Its use has the advantage of being fully nondestructive and safe for operators, and its ability to provide relevant information about different kinds of buried targets has been proved in several contexts.
However, there is a well-known issue connected to the GPR data interpretation [5]. Customarily, GPR data are provided to the user in the form of B-scans, i.e., a time-domain representation of the signals measured by GPR antennas in different spatial locations. In such a representation, buried targets cause the presence of reflection features that, typically, have a hyperbolic shape [6].
The problem of data interpretation is further complicated in the case of multiantenna GPR systems, i.e., composed of a set of antennas that alternately operate in transmitting and receiving mode (for example, one of them acts in turn as the source, and all the remaining ones receive the scattered fields). On one hand, multiantenna GPRs are able to provide much more information about buried targets, even with the device operating in a single position (since, for each location of the source, scattered field responses can be collected at various points). On the other hand, the obtained received signals versus time are more complex and they are difficult to interpret directly.
These issues motivate the use of suitable processing techniques to extract target information from GPR data [7]. Several GPR processing techniques have been developed in the past few decades, and they can be roughly subdivided into two main classes: qualitative approaches (i.e., aimed at obtaining only some information about targets, such as their position, shape, and size) [8,9,10], and quantitative methods (i.e., whose goal is a complete electromagnetic characterization of the targets of interest) [11,12,13].
The second class of methods, in particular, will be considered in the present work. Their objective is very interesting because of the possibility of providing key information about the target materials, thus providing a precise characterization of the underground environment [14]. Nevertheless, quantitative inversion is very challenging. This is because of (a) the nonlinear nature of the link between the electromagnetic field scattered by the buried objects and their properties, and (b) several complications arising from the measurement setup.
To address point (a), many techniques have been proposed in the scientific literature to solve the nonlinear inverse scattering problem [15,16,17,18,19,20]. Most of them, due to the high computational cost of solving this problem, are based on two-dimensional (2D) formulations [21,22,23,24]. They have the benefit of reducing the electromagnetic problem to a 2D scalar one, saving a significant amount of the required computational resources, but at the same time, they are based on significant assumptions on the electromagnetic model.
Another key aspect of GPR data inversion is point (b). Indeed, the measurement configurations typically adopted in GPR data collection only allow for measurements above the region under testing and, in particular, above the air–soil interface, or in borehole configurations, with a consequent strong reduction in the available information compared to standard microwave tomography [25]. Moreover, antennas generally operate in near-field conditions, and this leads to significant mutual interactions among them and strong interactions with targets [26]. Neglecting these effects, i.e., directly treating 3D data gathered with real antennas as 2D canonical data, may lead to inaccurate inversion results due to the disagreement between the actual and modeled configuration (i.e., the so-called modeling error).
To address this kind of problem, adequate techniques must be used. For example, some frequency-domain calibration schemes have been proposed in [27,28,29]. Recently, neural networks are also emerging to perform this task [30,31,32]. In particular, [31] introduced time-domain preprocessing of quasi-monostatic data acquired in free space configurations using long short-term memory (LSTM) cells [33]. LSTM cells are a specific kind of recurrent neural network that are attracting researchers, even in the field of electromagnetics, due to their impressive ability to process time series of data [12,34,35]. In the proposed approach, an LSTM cell, combined with a fully connected neural network, processes the time-domain data measured from the receiving antenna to perform a modeling error correction, i.e., to estimate the corresponding measurements as they were given by a 2D canonical configuration with ideal probes and sources.
In this work, the scheme initially proposed in [31] is extended to deal with a multiantenna GPR configuration for the detection of buried objects. Dealing with such configurations implies a change in both the inputs and outputs of the proposed neural network, which are now arrays of measurements containing the scattered-field responses for all the possible transmitting/receiving antenna pairs, as well as in the network training.
The proposed technique is evaluated in connection with a multifrequency Newton-based nonlinear imaging method [36]. The inversion algorithm has a variable exponent L p space formulation [37,38] to exploit the well-known benefits of performing regularization outside the conventional context of Hilbert spaces [39,40]. The inversion procedure is formulated under 2D canonical settings. Its input are the scattered-field data corrected by the LSTM-based method, and its output is a 2D quantitative image of the dielectric properties of the targets under testing.
Both the modeling error correction approach and the capability of obtaining a dielectric reconstruction with the proposed techniques have been assessed in a simulated scenario involving 3D configurations with cylindrical targets, where the forward electromagnetic problem is solved using the finite-difference time-domain (FDTD) method.
The paper is organized as follows. Section 2 provides the formulation of the LSTM-based GPR data processing method and outlines the quantitative inversion approach. Section 3 describes the obtained results. The discussion is reported in Section 4.

2. Formulation of the Method

The configuration of the measurement setup and a basic scheme of the proposed processing approach are shown in Figure 1. A half-space configuration is considered, where the ground ( y < 0 ) is modeled as a homogeneous material characterized by a dielectric permittivity ϵ g and electric conductivity σ g , and air ( y > 0 ) is characterized by dielectric permittivity ϵ 0 equal to that of a vacuum. The magnetic permeability is assumed to be equal to the vacuum value μ 0 everywhere. Data are collected by a multiantenna GPR system composed of a set of S antennas placed at a height h above ground level. Such antennas allow for the collection of multiview and multistatic data, i.e., one of the antennas, in turn, acts as a source, and the remaining K = S 1 elements are employed in receiving mode as measurement probes. In particular, K scalar measurements in the time domain are gathered for each position of the source.

2.1. Modeling Error Correction Approach Based on LSTM Cells

Let us denote the scalar response measured in the time domain by the k th antenna of the s th view ( k = 1 ,   ,   K ; s = 1 ,   , S ) as u ~ s c s , k ( t ) . In our case, we consider
u ~ s c s , k t = e ~ z , s c ( s )   r m s ( s , k ) , t ,
i.e., the measured scalar responses are the z components of the electric field in the k -th measurement point of the s -th view r m s ( s , k ) . The set of all points r m s ( s , k ) ,   k = 1 ,   ,   K constitutes the measurement domain for the s -th view, D m e a s ( s ) .
An array that groups all the available scalar measurements at the antenna ports for all the considered views can be defined as
U ~ s c t = u ~ s c 1,1 t ,   , u ~ s c 1 , S 1 t , , u ~ s c S , 1 t ,   , u ~ s c S , S 1 t   T .
The size of this array is equal to M = S S 1 , that is, the total number of multiantenna measurements for all the possible pairs of transmitting and receiving elements.
The array U ~ s c t , which is a collection of all the multistatic and multiview measurements of the scattered field at the time t , represents the data used by the proposed method. In the practical implementation, the time-domain data are collected in a time window with t 0 ,   T in a set of discrete time steps separated by an interval Δ t . Therefore, the used data have the form of a numeric matrix (i.e., the multistatic–multiview GPR B-scan), which is processed line-by-line starting from t = 0 . For the sake of simplicity, in the following formulation, we describe the processing approach focusing on a generic time instant t .
Let us consider an LSTM cell with U ~ s c t as the input (Figure 2). This cell is described by two recurrent vectors called cell state  c U ~ s c , t R H and hidden state  h U ~ s c , t R H , where H is the state size. The vectors c and h are used to “selectively remember” the information collected from the previous time history of the input signal and come from the previously processed time step.
The internal structure of the LSTM is composed of four operators, namely F , G , I , and O , which are conventionally called forget gate, cell candidate, input gate, and output gate, respectively. Each of these operators sums, with proper weights, the input array at time t , U ~ s c t , with the hidden state taken from the previous time step, h U ~ s c , t Δ t and adds a bias vector. The output is then passed to a gate activation function (sigmoid or hyperbolic tangent). In particular:
  • The function of the forget gate is defined as
F U ~ s c , t = S ( W F U ~ s c t + R F h U ~ s c , t Δ t + b F ) ,
where W F R H × M , R F R H × H , and b F R H are real-valued parameters, and S is a sigmoid gate activation function, i.e., considering a vector a = a 1 , ,   a H T it is defined as
S a = 1 + e a 1 1 , 1 + e a 1 2 , ,   1 + e a H 1 T .
  • The cell candidate G is described as
G U ~ s c , t = T W G U ~ s c t + R G h U ~ s c , t Δ t + b G ,
where
T a = tanh a 1 , tanh a 2 , ,   tanh a S T
is the function used for gate activation, and related parameters are W G R H × M , R G R H × H , b G R S .
  • The input gate I is defined as
I U ~ s c , t = S W I U ~ s c t + R I h U ~ s c , t Δ t + b I ,
with W I R H × M , R I R H × H and b I R H .
  • The function of the output gate O is
O U ~ s c , t = S W O U ~ s c t + R O h U ~ s c , t Δ t + b O ,
with W O R H × M , R O R H × H , and b O R H .
Based on the output of these four gates, the cell state c U ~ s c , t and hidden state h U ~ s c , t are, respectively, calculated as
c U ~ s c , t = F U ~ s c , t     c U ~ s c , t Δ t + I U ~ s c , t     G U ~ s c , t ,
h U ~ s c , t = O U ~ s c , t     T c U ~ s c , t ,
where indicates the Hadamard product. The hidden state h U ~ s c , t is then processed by a fully connected neural network composed of three layers, where the first two layers are composed of L neurons, and the last layer has M neurons (since the final output should have the same size as the input array, i.e., the total number of multiantenna measurements). In more detail, the first layers provide the following outputs:
N 1 U ~ s c , t = T W N 11 T h U ~ s c , t + b N 11 ,   ,   W N 1 L T h U ~ s c , t + b N 1 L T
N 2 U ~ s c , t = T W N 21 T N 1 U ~ s c , t + b N 21 ,   ,   W N 2 L T N 1 U ~ s c , t + b N 2 L T
where W N 11 ,   ,   W N 1 L R H ,   W N 21 ,   ,   W N 2 L R L ,   b N 11 ,   ,   b N 1 L R , b N 21 ,   ,   b N 2 L R . The final output is given by
N U ~ s c , t = W N 31 T N 2 U ~ s c , t + b N 31 ,   ,   W N 3 M T N 2 U ~ s c , t + b N 3 M T
with W N 31 ,   ,   W N 3 M R L ,   b N 31 ,   ,   b N 3 M R .
The key point is that N U ~ s c , t is used to approximate the set of multiantenna measurements that would be obtained if both sources and probes are ideal, in 2D settings. In other words, by suitable training of the network, we enforce that
U s c t = N U ~ s c , t
where
U s c t = u s c 1,1 t ,   , u s c 1 , S 1 t , , u s c S , 1 t ,   , u s c S , S 1 t   T .
is an array of the same size as U ~ s c that contains the corresponding set of measurements of the scattered electric field u s c s , k t = e z , s c ( s )   r m s ( s , k ) , t as obtained by ideal non-perturbing field probes, in the presence of infinite line current sources, and in a 2D configuration, i.e., with transverse magnetic (TM) z -polarized fields and targets approximated as cylinders with uniform properties on the z axis.
With this definition of inputs ( U ~ s c ) and outputs ( U s c ), the neural-network-based operator N defined by Equations (3)–(14) is trained to be able to correct the modeling error, extracting from non-ideal measurements an approximation of the measured values obtained in canonical 2D conditions. As in [31], the neural network coefficients inside the definition of the operator N are found by a training procedure performed with pairs of time-domain electromagnetic simulations, where the same set of configurations is simulated in both 3D realistic and 2D canonical settings.

2.2. Nonlinear Inverse Scattering Method

The data obtained by the proposed approach are then processed by a tomographic nonlinear inverse scattering method that aims at reconstructing the dielectric properties of targets contained in a given 2D subsurface region of interest ( D s u b , see Figure 1). In more detail, a Fast Fourier Transform (FFT) algorithm is applied to the obtained approximated values of e z , s c ( s )   r m s ( s , k ) , t to extract the frequency-domain scattered field at frequency f f 1 , ,   f F (i.e., F frequencies). These values, indicated as E z , s c ( s )   r m s ( s , k ) , f (an e j 2 π f t time dependance is assumed and dropped), are used as input data by the inverse scattering algorithm, which assumes that incident electromagnetic fields are TM- z and that buried objects are infinite cylinders along the z axis, characterized by dielectric properties ϵ r t and σ ( r t ) , with r t = x x ^ + y y ^ D s u b . With these assumptions, the z component of the scattered electric field at the measurement points can be connected to the dielectric properties of the buried targets by the following relationship (data equation):
E z , s c ( s )   r m s ( s , k ) , f = k g 2 D s u b T f x r t E z , t o t s r t , f g s r m s ( s , k ) , r t , f d r t ,   s = 1 ,   ,   S k = 1 ,   ,   K
where
x r t = ϵ r t ϵ g a ϵ σ r t σ g a σ T ,   r t D s u b
represents the unknown of the inverse problem, T f = a ϵ ϵ 0 / ϵ g , a σ j 2 π f ϵ g 1 , a ϵ , a σ are normalization coefficients [37], k g 2 = 4 π 2 f 2 μ 0 ϵ g , g s is the half-space Green’s function [41], and finally E z , t o t s r t , f is the z component of the total electric field inside the subsurface investigation domain. The latter quantity is also unknown and depends on x , making the problem nonlinear. Therefore, Equation (16) needs to be combined with the state equation
E z , t o t s r t , f = E z , i n c s r t , f k g 2 D s u b T f x r t E z , t o t s r t , f g s r t , r t , f d r t ,   r t D s u b
where E z , i n c s r t , f = j 2 π f μ 0 I g s r t , r t x s , f is the z component of the incident electric field in the s -th view due to a line current source located at r t x s .
Considering all the available measurement points and frequencies, the inverse problem is then formulated as a system of equations [36] given by
E z , s c ( 1 )   r m s 1,1 , f 1 E z , s c ( S )   r m s S , K , f F = k g 2 D s u b T f 1 x r t E z , t o t 1 r t , f 1 g s r m s ( s , k ) , r t , f 1 d r t k g 2 D s u b T f F x r t E z , t o t S r t , f F g s r m s ( S , K ) , r t , f F d r t .
This system can be written more compactly as
y = F x ,
where y Y is the array that contains all the measured scattered fields, x X is the unknown, and F : X Y is a nonlinear operator.
The solution of the inverse scattering problem described in Equation (20) is carried out by a Newton-based inversion algorithm [24]. In particular, the algorithm is iterative and composed of an external linearization loop that, in each j -th step, calculates the Fréchet derivative of the operator F around a current tentative solution x j ( F x j ), leading to the linearized equation
y F x j F x j Δ x j = 0 .
The unknown of the linearized problem at the j -th step is then Δ x j , which is found by a truncated Landweber-like regularization approach [37], formulated by considering X as a variable exponent L p j r t space. The exponent p j r t is an adaptive function computed as p 1 r t = 2 Δ p in the first step ( Δ p is a parameter describing the exponent variation) and p j r t = 2 + Δ p x j 1 r t / max r t D sub x j 1 r t 1 , r t D s u b in the following ones. Conversely, Y is assumed as an L p a space, p a , being equal to the mean value of p j . After the solution of the linearized equation, the global solution is updated as x j + 1 = x j + Δ x j , and the process is iterated until convergence criteria are met.

3. Results

The proposed technique has been tested within a simulated environment, where the scattered-field data have been calculated with the FDTD method implemented by the gprMax open-source software (version: 3.1.7) [42].
First of all, as mentioned in Section 2.1, the modeling error correction network has been trained through pairs of time-domain simulations performed in 3D realistic settings (to obtain U ~ s c t ), and in the corresponding 2D canonical conditions (to obtain U s c t ).
An example of the configurations adopted in these paired FDTD simulations is given in Figure 3. In both situations, the soil (i.e., the lower half space with y < 0 ) has been characterized as a homogeneous material with randomly chosen dielectric properties ϵ g 3 ,   5 ϵ 0 and σ g 5 × 10 3 ,   1.5 × 10 2   S / m . Random circular cylinders placed inside the ground with z -directed axes have been considered, characterized by dielectric properties ϵ 1 ,   9 ϵ 0 and σ 0 ,   3 × 10 2 S / m , centered at coordinates x c 0.625 ,   0.625   m and y c 0.675 ,   0.02   m , with diameter d 0.06 ,   0.32   m .
Scattered field data are collected by a set of S = 10 antennas above the soil with uniform spacing on an x -directed measurement line of length L = 1.1   m at height h = 0.05   m over the soil level. In more detail, in the s -th view, the transmitting antenna is centered at
r t x s = L s 1 S 1 d s , h , 0
with d s = 0.55   m , and the K = 9 receiving antennas are located at
r m s s , k = L k 1 S 1 d s , h , 0 ,   k = 1 ,   ,   s 1 L k S 1 d s , h , 0 ,   k = s ,   ,   K .
In the 3D simulations, antennas have a bowtie shape with size l = 0.05   m , flare angle θ = 90 ° , are parallel to the ground and polarized on the z axis. The FDTD simulation domain has x , y , and z sides of length D x = 1.5   m , D y = 1   m , and D z = 0.5   m , respectively. It has been discretized into 6 × 10 6 cubic cells ( 300 × 200 × 100 ) with sides 5   mm in length.
The vertical slice of the 3D simulation domain located on the x y plane at z = 0 has also been the subject of 2D simulations, which have been carried out with a rectangular domain of sides D x = 1.5   m and D y = 1   m discretized into 6 × 10 4 square cells ( 300 × 200 ) with sides 5   mm in length. In the 2D case, transmitting antennas have been represented by infinite z -directed line current sources, and the scattered field has been sampled by ideal nonperturbing field probes.
In both 2D and 3D conditions, the simulation domains are terminated by a 10-cell wide perfectly matched layer. A time window of length T = 15   ns has been simulated, with a time step of 9.63   ps . Source antennas have been excited by a Ricker waveform with a central frequency of 1   GHz .
Simulations have been carried out on a Lubuntu-based personal computer equipped with an 8-core Intel® (Santa Clara, CA, USA) Core™ i7-11700F central processing unit (CPU) with 2.50 GHz base clock frequency, 16 GB of random-access memory (RAM), and an NVIDIA (Santa Clara, CA, USA) GeForce RTX 3070 graphical processing unit (GPU) with 8 GB of dedicated RAM (CUDA version: 12.2). A single 3D case (i.e., S = 10 simulations with different positions of the source) requires 50 s to run and 485 MB RAM on the host plus 617 MB on the GPU, whereas 4 s, 189 MB RAM, and 271 MB GPU RAM are required to run a 2D case.
A total number of 400 cases have been simulated in 3D and 2D settings (i.e., 8000 total FDTD simulations, 72,000 A-scans), where all the random variables have been defined with uniform distributions in the given ranges. Among them, 397 cases have been used inside the training set, and the remaining ones have been considered for validation (corresponding to 60 different FDTD simulations, 540 A-scans). The network has been trained by the ADAM algorithm [43] for 100 epochs, with mean squared error loss function, a mini-batch size of 50, and a learning rate equal to 10 2 .

3.1. Assessment of the Modeling Error Correction Approach

The neural network used to perform modeling error correction has been assessed versus some of the main hyperparameters, i.e., the LSTM cell state size (indicated in Section 2.1 as the variable H ) and the number of neurons present in each fully connected (FC) layer (i.e., the variable L ). In particular, the values H 100 ,   200 ,   300 and L = 25 ,   50 ,   100 ,   200 have been considered.
The accuracy of the reconstruction of the canonical scattered-field data starting from 3D simulations has been evaluated by means of the relative root mean squared error
e U = 0 T U s c t U s c a c t t d t 0 T U s c a c t t d t
where U s c t is the reconstructed scattered-field measurement array obtained from the neural network, and U s c a c t t contains the corresponding actual values from 2D simulations. The obtained average of e U in the validation set versus LSTM state size and number of FC neurons of the two hidden layers is reported in Figure 4. As can be seen, for all values of L , the error always decreases by increasing the state size up to H = 300 . Conversely, if the FC neurons are increased, an initial decrease in e U is observed, followed by a subsequent increase. The lowest error in the reconstruction of scattered field data has been obtained with LSTM state size of H = 300 and with L = 100 neurons in the two hidden layers. Therefore, these values have been selected for the following analysis.
In order to assess the performance of the proposed approach, the behavior of the processing network has been tested in another configuration not included in the training and validation set. In this case, soil is characterized by ϵ g = 4 ϵ 0 and σ g = 10 2   S / m , and a circular dielectric cylinder with ϵ = 8.5 ϵ 0 , σ = 2 × 10 2   S / m , diameter d = 0.12   m , and centered at x c = 0.3   m , y c = 0.15   m has been considered. Also in this case, both 3D and 2D simulations of the same target have been performed for comparison purposes.
The multiantenna GPR B-scan resulting from 3D simulations (i.e., U ~ s c t ) is shown in Figure 5a (measurement point index reported in abscissa is computed as k + s 1 K ). The corresponding quantity obtained from 2D simulations ( U s c a c t t ) and the result of the LSTM-based network applied to 3D data ( U s c t ) are shown in Figure 5b and Figure 5c, respectively. The difference between the 2D canonical data and corrected 3D data, U s c t U s c a c t t , is reported in Figure 5d. Moreover, an example of an A-scan is reported in Figure 5e to better compare the results. In general, the proposed method allows for a good estimation of the scattered field in canonical conditions, although it is affected by some noise when the amplitude of the scattered field is very low. As can be seen in Figure 5d, the estimation of the first reflection is always better than that of the second one, which presents some differences compared to the canonical data. In addition, the reconstruction of the second reflection seems slightly less accurate than the first one. In this case, the value of the relative root mean squared error on the scattered field is e U = 0.66 .

3.2. Reconstruction of Buried Dielectric Targets

The B-scans presented in Figure 5 have also been exploited to obtain a quantitative reconstruction of the dielectric properties of the buried target by adopting the inversion method outlined in Section 2.2.
To this end, scattered-field data have been extracted at F = 12 frequencies equally spaced between 100   MHz and 600   MHz . Data have been corrupted with an additive white Gaussian noise with mean value equal to zero and a signal-to-noise ratio equal to 20   dB . The investigation domain D s u b is a rectangular region on the x y plane between x m i n = 0.5   m , y m i n = 0.55   m , and x m a x = 0.5   m , y m a x = 0.05   m . The region has been discretized into 800 square cells ( 40 × 20 ) with sides 2.5   cm in length.
The inversion method has been executed with exponent function variation Δ p = 0.5 and the following stopping criteria: 20 maximum Newton iterations, 10 maximum Landweber iterations, minimum relative variation in the data residual between two successive steps equal to 10 2 .
The reconstructed distributions of the relative dielectric permittivity obtained by applying the inversion method to the raw 3D data, to the canonical 2D data, and to the 3D data corrected by the proposed approach are shown in Figure 6. In particular, the results achieved by the proposed variable exponent Lebesgue space inversion method are compared with those obtained with a nonlinear Hilbert space inversion algorithm [44]. It is evident, as expected, that both inversion methods applied to the raw 3D data do not lead to satisfactory results (Figure 6b,c). More correct estimations of the target properties are obtained using the 2D data (Figure 6d,e). Compared to the Hilbert space method (Figure 6e), the proposed variable exponent approach allows for a significantly better characterization of the buried cylinder (Figure 6d). However, the 2D data used to produce these reconstructions come from a simplified solution of the forward problem, which closely matches the one adopted for the inversion, but is also quite far from realistic settings. Figure 6f,g show the reconstructions obtained by processing the corrected 3D data. Clearly, with the variable exponent inversion method (Figure 6f), the dielectric permittivity of the buried target, although slightly overestimated, is retrieved in a more accurate way with respect to the Hilbert space result (Figure 6g). Overall, the quality of the reconstruction is not far from the one obtained with canonical 2D data. However, some artifacts appear. In particular, two regions with reconstructed dielectric permittivity lower than the background appear above ( y 0.05   m ) and below the target ( y 0.25   m ). Moreover, a barely visible ghost target is present close to the point x 0.05   m , y 0.15   m , with an increased value of permittivity compared to the background. Some minor ringing effects can be seen surrounding these regions. The insurgence of such artifacts can be mainly ascribed to the differences between the actual 2D data and their approximation obtained with the proposed LSTM-based network, visible in Figure 5d. In more detail, these artifacts are related to the imperfections in the retrieval of the second reflection features (left part of Figure 5d) and to the presence of some noise where scattered-field values are relatively low (right part of Figure 5d).
A quantitative assessment of the reconstruction performance has been obtained by computing the average relative reconstruction error on the contrast function in the target region ( e T ) and in the background ( e B ) as
e T = 1 A D T D T χ f 1 r t χ f 1 a c t r t χ f 1 a c t r t + 1 d r t
e B = 1 A D B D B χ f 1 r t χ f 1 a c t r t χ f 1 a c t r t + 1 d r t
where the integrals are computed in the region occupied by the target ( D T ) and in the background ( D B = D s u b \ D T ), A D T and A D B indicate the corresponding areas, χ f 1 r t = T f 1 x r t is the reconstructed value of the contrast function at the first frequency f 1 (considered as reference), and χ f 1 a c t r t is its actual value.
The computed values of e T and e B obtained with raw 3D data, canonical 2D data and 3D data corrected by the proposed approach are reported in Table 1. Moreover, a comparison has been carried out between the variable exponent Lebesgue space inversion algorithm and the Hilbert space method.
As confirmed from the observation of the reconstructed dielectric permittivity in Figure 6, when raw 3D data are considered, the error on the background is minimal but the error on the target is large (i.e., the target is not detected correctly). In the other two cases (i.e., reconstruction with 2D canonical data and with corrected 3D data) the error on the target is lower and quite similar, although slightly higher in the last case. Analyzing the reported values, it also appears that the variable exponent Lebesgue space inversion algorithm generally provides lower reconstruction errors, especially in the target region.

4. Discussion

In this paper, a methodology to correct the modeling error by exploiting a neural-network-based approach has been considered. The network operates in the time domain by exploiting LSTM cells. This kind of processing method, previously evaluated in simpler configurations, has been applied here for the first time to process scattered-field data acquired by a multiantenna GPR measurement setup.
The neural network has been trained with pairs of 3D and 2D numerical simulations, where the 3D ones involved a set of bowtie antennas placed above the ground, which are replaced by ideal line current sources and ideal field probes in 2D settings. Testing has also been carried out in a simulated environment, after an evaluation of the behavior of the processing network by varying some of its main hyperparameters.
The results are quite promising and show that the proposed architecture is able to extract from 3D realistic data an acceptable estimation of the scattered field, as measured in ideal canonic conditions. The processed data, if provided in input to a quantitative inverse-scattering procedure, give rise to a dielectric reconstruction that is not far from the one achieved with ideal 2D data, and much better than the one obtained by directly processing the raw 3D data. In particular, a nonlinear Newton-based inversion algorithm formulated in the framework of variable exponent Lebesgue spaces has been considered here, but the proposed processing method, in principle, can be applied to every inversion procedure that aims at reconstructing the properties of unknown buried targets starting from measurements of the scattered field or related scalar responses (e.g., S-parameters) or vector quantities.
Future research directions may include the test of the proposed approach in different configurations and in conjunction with other inverse scattering approaches. Moreover, the assessment of the technique in the presence of real experimental data will be considered. Finally, a very important step will be the validation of the proposed method in practical scenarios (e.g., detection of buried objects in various subsurface configurations, road assessment, analysis of civil structures).

Author Contributions

Conceptualization, A.F., V.S., and A.R.; methodology, A.F., V.S., and A.R.; software, A.F., V.S., and A.R.; validation, A.F., V.S., and A.R.; writing—original draft preparation, A.F., V.S., and A.R.; writing—review and editing, A.F., V.S., and A.R. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors on request.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Benedetto, A.; Pajewski, L. Civil Engineering Applications of Ground Penetrating Radar; Springer: Cham, Switzerland, 2015; ISBN 978-3-319-04813-0. [Google Scholar]
  2. Busch, S.; Weihermüller, L.; Huisman, J.A.; Steelman, C.M.; Endres, A.L.; Vereecken, H.; Van, D.K. Coupled Hydrogeophysical Inversion of Time-Lapse Surface GPR Data to Estimate Hydraulic Properties of a Layered Subsurface. Water Resour. Res. 2013, 49, 8480–8494. [Google Scholar] [CrossRef]
  3. Goodman, D.; Piro, S. GPR Remote Sensing in Archaeology; Geotechnologies and the Environment; Springer: Berlin, Germany, 2013; ISBN 978-3-642-31856-6. [Google Scholar]
  4. Almeida, E.R.; Porsani, J.L.; Catapano, I.; Gennarelli, G.; Soldovieri, F. Microwave Tomography-Enhanced GPR in Forensic Surveys: The Case Study of a Tropical Environment. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2016, 9, 115–124. [Google Scholar] [CrossRef]
  5. Persico, R.; Ludeno, G.; Soldovieri, F.; De Coster, A.; Lambot, S. Improvement of Ground Penetrating Radar (GPR) Data Interpretability by an Enhanced Inverse Scattering Strategy. Surv. Geophys. 2018, 39, 1069–1079. [Google Scholar] [CrossRef]
  6. Persico, R. Introduction to Ground Penetrating Radar: Inverse Scattering and Data Processing; Wiley: Hoboken, NJ, USA, 2014; ISBN 978-1-118-30500-3. [Google Scholar]
  7. Catapano, I.; Randazzo, A.; Slob, E.; Solimene, R. GPR Imaging via Qualitative and Quantitative Approaches. In Civil Engineering Applications of Ground Penetrating Radar; Benedetto, A., Pajewski, L., Eds.; Springer International Publishing: Cham, Switzerland, 2015; pp. 239–280. ISBN 978-3-319-04812-3. [Google Scholar]
  8. Özdemir, C.; Demirci, Ş.; Yiğit, E.; Yilmaz, B. A Review on Migration Methods in B-Scan Ground Penetrating Radar Imaging. Available online: https://www.hindawi.com/journals/mpe/2014/280738/ (accessed on 4 January 2018).
  9. Ludeno, G.; Gennarelli, G.; Lambot, S.; Soldovieri, F.; Catapano, I. A Comparison of Linear Inverse Scattering Models for Contactless GPR Imaging. IEEE Trans. Geosci. Remote Sens. 2020, 58, 7305–7316. [Google Scholar] [CrossRef]
  10. Ambrosanio, M.; Bevacqua, M.T.; Isernia, T.; Pascazio, V. Performance Analysis of Tomographic Methods against Experimental Contactless Multistatic Ground Penetrating Radar. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2021, 14, 1171–1183. [Google Scholar] [CrossRef]
  11. Li, L.; Wang, L.G.; Teixeira, F.L.; Liu, C.; Nehorai, A.; Cui, T.J. DeepNIS: Deep Neural Network for Nonlinear Electromagnetic Inverse Scattering. IEEE Trans. Antennas Propag. 2019, 67, 1819–1825. [Google Scholar] [CrossRef]
  12. Wang, J.; Liu, H.; Jiang, P.; Wang, Z.; Sui, Q.; Zhang, F. GPRI2Net: A Deep-Neural-Network-Based Ground Penetrating Radar Data Inversion and Object Identification Framework for Consecutive and Long Survey Lines. IEEE Trans. Geosci. Remote Sens. 2022, 60, 1–20. [Google Scholar] [CrossRef]
  13. Schenone, V.; Fedeli, A.; Estatico, C.; Pastorino, M.; Randazzo, A. Microwave Imaging of Mixed Metallic-Dielectric Configurations via a Finite Element-Based Variable-Exponent Approach. URSI Radio Sci. Lett. 2021, 3, 1–5. [Google Scholar] [CrossRef] [PubMed]
  14. Di Donato, L.; Crocco, L. Model-Based Quantitative Cross-Borehole GPR Imaging via Virtual Experiments. IEEE Trans. Geosci. Remote Sens. 2015, 53, 4178–4185. [Google Scholar] [CrossRef]
  15. Abubakar, A.; Habashy, T.M.; Van den Berg, P.M. Nonlinear Inversion of Multi-Frequency Microwave Fresnel Data Using the Multiplicative Regularized Contrast Source Inversion. Prog. Electromagn. Res. 2006, 62, 193–201. [Google Scholar] [CrossRef]
  16. Gurbuz, T.U.; Aslanyurek, B.; Karabulut, E.P.; Akduman, I. An Efficient Nonlinear Imaging Approach for Dielectric Objects Buried under a Rough Surface. IEEE Trans. Geosci. Remote Sens. 2014, 52, 3013–3022. [Google Scholar] [CrossRef]
  17. Desmal, A.; Bağcı, H. A Preconditioned Inexact Newton Method for Nonlinear Sparse Electromagnetic Imaging. IEEE Geosci. Remote Sens. Lett. 2015, 12, 532–536. [Google Scholar] [CrossRef]
  18. Wei, Z.; Chen, X. Deep-Learning Schemes for Full-Wave Nonlinear Inverse Scattering Problems. IEEE Trans. Geosci. Remote Sens. 2019, 57, 1849–1860. [Google Scholar] [CrossRef]
  19. Schenone, V.; Estatico, C.; Gragnani, G.L.; Pastorino, M.; Randazzo, A.; Fedeli, A. Microwave-Based Subsurface Characterization through a Combined Finite Element and Variable Exponent Spaces Technique. Sensors 2023, 23, 167. [Google Scholar] [CrossRef] [PubMed]
  20. Estatico, C.; Schenone, V.; Fedeli, A.; Randazzo, A. Application of a Mild Data-Driven Technique to Lippmann-Schwinger Inverse Scattering in Variable-Exponent Lebesgue Spaces for Microwave Imaging. Inverse Probl. 2024, 40, 065007. [Google Scholar] [CrossRef]
  21. Salucci, M.; Oliveri, G.; Randazzo, A.; Pastorino, M.; Massa, A. Electromagnetic Subsurface Prospecting by a Fully Nonlinear Multifocusing Inexact Newton Method. J. Opt. Soc. Am. A 2014, 31, 2618. [Google Scholar] [CrossRef] [PubMed]
  22. Schenone, V.; Fedeli, A.; Estatico, C.; Pastorino, M.; Randazzo, A. Experimental Assessment of a Novel Hybrid Scheme for Quantitative GPR Imaging. IEEE Geosci. Remote Sens. Lett. 2022, 19, 1–5. [Google Scholar] [CrossRef]
  23. Takahashi, S.; Suzuki, K.; Hanabusa, T.; Kidera, S. Microwave Subsurface Imaging Method by Incorporating Radar and Tomographic Approaches. IEEE Trans. Antennas Propag. 2022, 70, 11009–11023. [Google Scholar] [CrossRef]
  24. Schenone, V.; Estatico, C.; Pastorino, M.; Randazzo, A.; Fedeli, A. Electromagnetic Imaging in Stratified Media by Means of a Finite-Element Variable-Exponent Inversion Approach. URSI Radio Sci. Lett. 2022, 4, 1–4. [Google Scholar] [CrossRef]
  25. Bucci, O.M.; Crocco, L.; Isernia, T.; Pascazio, V. Subsurface Inverse Scattering Problems: Quantifying, Qualifying, and Achieving the Available Information. IEEE Trans. Geosci. Remote Sens. 2001, 39, 2527–2538. [Google Scholar] [CrossRef]
  26. André, F.; Lambot, S. Intrinsic Modeling of Near-Field Electromagnetic Induction Antennas for Layered Medium Characterization. IEEE Trans. Geosci. Remote Sens. 2014, 52, 7457–7469. [Google Scholar] [CrossRef]
  27. Lambot, S.; Slob, E.C.; van den Bosch, I.; Stockbroeckx, B.; Vanclooster, M. Modeling of Ground-Penetrating Radar for Accurate Characterization of Subsurface Electric Properties. IEEE Trans. Geosci. Remote Sens. 2004, 42, 2555–2568. [Google Scholar] [CrossRef]
  28. Lambot, S.; André, F. Full-Wave Modeling of near-Field Radar Data for Planar Layered Media Reconstruction. IEEE Trans. Geosci. Remote Sens. 2014, 52, 2295–2303. [Google Scholar] [CrossRef]
  29. Ostadrahimi, M.; Mojabi, P.; Gilmore, C.; Zakaria, A.; Noghanian, S.; Pistorius, S.; LoVetri, J. Analysis of Incident Field Modeling and Incident/Scattered Field Calibration Techniques in Microwave Tomography. IEEE Antennas Wirel. Propag. Lett. 2011, 10, 900–903. [Google Scholar] [CrossRef]
  30. Zhang, J.W.; Ye, S.B.; Liu, H.; Yi, L.; Fang, G.Y. Filtering out Antenna Effects from GPR Data by an RBF Neural Network. IEEE Geosci. Remote Sens. Lett. 2019, 16, 1378–1382. [Google Scholar] [CrossRef]
  31. Fedeli, A. Microwave Tomography with LSTM-Based Processing of the Scattered Field. IEEE Open J. Antennas Propag. 2021, 2, 213–223. [Google Scholar] [CrossRef]
  32. Hanabusa, T.; Morooka, T.; Kidera, S. Deep-Learning-Based Calibration in Contrast Source Inversion Based Microwave Subsurface Imaging. IEEE Geosci. Remote Sens. Lett. 2022, 19, 1–5. [Google Scholar] [CrossRef]
  33. Greff, K.; Srivastava, R.K.; Koutník, J.; Steunebrink, B.R.; Schmidhuber, J. LSTM: A Search Space Odyssey. IEEE Trans. Neural Netw. Learn. Syst. 2017, 28, 2222–2232. [Google Scholar] [CrossRef] [PubMed]
  34. Lei, W.; Luo, J.; Hou, F.; Xu, L.; Wang, R.; Jiang, X. Underground Cylindrical Objects Detection and Diameter Identification in GPR B-Scans via the CNN-LSTM Framework. Electronics 2020, 9, 1804. [Google Scholar] [CrossRef]
  35. Noakoasteen, O.; Wang, S.; Peng, Z.; Christodoulou, C. Physics-Informed Deep Neural Networks for Transient Electromagnetic Analysis. IEEE Open J. Antennas Propag. 2020, 1, 404–412. [Google Scholar] [CrossRef]
  36. Fedeli, A.; Estatico, C.; Randazzo, A.; Pastorino, M. Multifrequency Microwave Tomography in Lebesgue Spaces with Nonconstant Exponents. URSI Radio Sci. Lett. 2020, 2, 1–4. [Google Scholar] [CrossRef] [PubMed]
  37. Estatico, C.; Fedeli, A.; Pastorino, M.; Randazzo, A. Quantitative Microwave Imaging Method in Lebesgue Spaces with Nonconstant Exponents. IEEE Trans. Antennas Propag. 2018, 66, 7282–7294. [Google Scholar] [CrossRef]
  38. Fedeli, A.; Schenone, V.; Randazzo, A.; Pastorino, M.; Henriksson, T.; Semenov, S. Nonlinear S-Parameters Inversion for Stroke Imaging. IEEE Trans. Microw. Theory Tech. 2021, 69, 1760–1771. [Google Scholar] [CrossRef]
  39. Schuster, T.; Kaltenbacher, B.; Hofmann, B.; Kazimierski, K.S. Regularization Methods in Banach Spaces; De Gruyter: Berlin, Germany, 2012; ISBN 978-3-11-025524-9. [Google Scholar]
  40. Schenone, V.; Fedeli, A.; Estatico, C.; Pastorino, M.; Randazzo, A. Detection of Failures in Antenna Arrays through a Lebesgue-Space Approach. IEEE Open J. Antennas Propag. 2022, 3, 652–662. [Google Scholar] [CrossRef]
  41. Stinson, D.C. Intermediate Mathematics of Electromagnetics; Prentice-Hall: Englewood Cliffs, NJ, USA, 1976; ISBN 0-13-470633-1. [Google Scholar]
  42. Warren, C.; Giannopoulos, A.; Giannakis, I. gprMax: Open Source Software to Simulate Electromagnetic Wave Propagation for Ground Penetrating Radar. Comput. Phys. Commun. 2016, 209, 163–170. [Google Scholar] [CrossRef]
  43. Kingma, D.P.; Ba, J. ADAM: A Method for Stochastic Optimization. arXiv 2017, arXiv:1412.6980v9. [Google Scholar]
  44. Estatico, C.; Fedeli, A.; Pastorino, M.; Randazzo, A. Microwave Imaging by Means of Lebesgue-Space Inversion: An Overview. Electronics 2019, 8, 945. [Google Scholar] [CrossRef]
Figure 1. Configuration of the GPR multiantenna measurement setup (considering, as an example, the s -th view) and basic scheme of the proposed processing approach.
Figure 1. Configuration of the GPR multiantenna measurement setup (considering, as an example, the s -th view) and basic scheme of the proposed processing approach.
Remotesensing 16 02050 g001
Figure 2. Proposed modeling error correction strategy based on the use of an LSTM cell.
Figure 2. Proposed modeling error correction strategy based on the use of an LSTM cell.
Remotesensing 16 02050 g002
Figure 3. Example of paired FDTD simulations adopted to train the neural network for modeling error correction: (a) 3D simulation with bowtie antennas; (b) 2D simulation of the same case with line-current sources and ideal field probes.
Figure 3. Example of paired FDTD simulations adopted to train the neural network for modeling error correction: (a) 3D simulation with bowtie antennas; (b) 2D simulation of the same case with line-current sources and ideal field probes.
Remotesensing 16 02050 g003
Figure 4. Average relative root mean squared error on the reconstruction of the scattered electric field by the proposed neural network ( e U ) versus the number of neurons of FC layers ( L ) and the LSTM cell state size ( H ).
Figure 4. Average relative root mean squared error on the reconstruction of the scattered electric field by the proposed neural network ( e U ) versus the number of neurons of FC layers ( L ) and the LSTM cell state size ( H ).
Remotesensing 16 02050 g004
Figure 5. Assessment of the processing of scattered field by the proposed method: (a) B-scan of the raw 3D data, U ~ s c t ; (b) B-scan of the 2D canonical data, U s c a c t t ; (c) B-scan of the 3D data corrected by the LSTM-based method, U s c t ; (d) difference between 2D canonical data and corrected 3D data, U s c t U s c a c t t ; (e) example of A-scan ( s = 1 , k = 3 ).
Figure 5. Assessment of the processing of scattered field by the proposed method: (a) B-scan of the raw 3D data, U ~ s c t ; (b) B-scan of the 2D canonical data, U s c a c t t ; (c) B-scan of the 3D data corrected by the LSTM-based method, U s c t ; (d) difference between 2D canonical data and corrected 3D data, U s c t U s c a c t t ; (e) example of A-scan ( s = 1 , k = 3 ).
Remotesensing 16 02050 g005
Figure 6. Distributions of the relative dielectric permittivity inside the subsurface investigation domain: (a) actual values; reconstruction using raw 3D data: (b) variable exponent Lebesgue space inversion algorithm, and (c) Hilbert space inversion algorithm; reconstruction using 2D canonical data: (d) variable exponent Lebesgue space inversion algorithm, and (e) Hilbert space inversion algorithm; reconstruction using 3D data corrected by the proposed approach: (f) variable exponent Lebesgue space inversion algorithm, and (g) Hilbert space inversion algorithm. White circle represents the true profile of the target.
Figure 6. Distributions of the relative dielectric permittivity inside the subsurface investigation domain: (a) actual values; reconstruction using raw 3D data: (b) variable exponent Lebesgue space inversion algorithm, and (c) Hilbert space inversion algorithm; reconstruction using 2D canonical data: (d) variable exponent Lebesgue space inversion algorithm, and (e) Hilbert space inversion algorithm; reconstruction using 3D data corrected by the proposed approach: (f) variable exponent Lebesgue space inversion algorithm, and (g) Hilbert space inversion algorithm. White circle represents the true profile of the target.
Remotesensing 16 02050 g006
Table 1. Average relative reconstruction errors in target and background regions obtained using 3D data, 2D data, and 3D data corrected by the proposed approach inside the quantitative inversion algorithm. Comparison between inversion algorithm formulated in variable exponent Lebesgue spaces and Hilbert space algorithm.
Table 1. Average relative reconstruction errors in target and background regions obtained using 3D data, 2D data, and 3D data corrected by the proposed approach inside the quantitative inversion algorithm. Comparison between inversion algorithm formulated in variable exponent Lebesgue spaces and Hilbert space algorithm.
Inversion AlgorithmCaseTarget Error, e T Background Error, e B
Variable exponent Lebesgue spaceRaw 3D data0.513470.01615
2D data0.263520.05091
Corrected 3D data0.299920.03832
Hilbert spaceRaw 3D data0.525410.02063
2D data0.307970.05646
Corrected 3D data0.421720.03629
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Fedeli, A.; Schenone, V.; Randazzo, A. Quantitative Inversion of Multiantenna Ground-Penetrating Radar Data with Modeling Error Correction Based on Long Short-Term Memory Cells. Remote Sens. 2024, 16, 2050. https://doi.org/10.3390/rs16122050

AMA Style

Fedeli A, Schenone V, Randazzo A. Quantitative Inversion of Multiantenna Ground-Penetrating Radar Data with Modeling Error Correction Based on Long Short-Term Memory Cells. Remote Sensing. 2024; 16(12):2050. https://doi.org/10.3390/rs16122050

Chicago/Turabian Style

Fedeli, Alessandro, Valentina Schenone, and Andrea Randazzo. 2024. "Quantitative Inversion of Multiantenna Ground-Penetrating Radar Data with Modeling Error Correction Based on Long Short-Term Memory Cells" Remote Sensing 16, no. 12: 2050. https://doi.org/10.3390/rs16122050

APA Style

Fedeli, A., Schenone, V., & Randazzo, A. (2024). Quantitative Inversion of Multiantenna Ground-Penetrating Radar Data with Modeling Error Correction Based on Long Short-Term Memory Cells. Remote Sensing, 16(12), 2050. https://doi.org/10.3390/rs16122050

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop