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
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 (
) is modeled as a homogeneous material characterized by a dielectric permittivity
and electric conductivity
, and air (
) is characterized by dielectric permittivity
equal to that of a vacuum. The magnetic permeability is assumed to be equal to the vacuum value
everywhere. Data are collected by a multiantenna GPR system composed of a set of
antennas placed at a height
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
elements are employed in receiving mode as measurement probes. In particular,
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
th antenna of the
th view (
) as
. In our case, we consider
i.e., the measured scalar responses are the
components of the electric field in the
-th measurement point of the
-th view
. The set of all points
constitutes the
measurement domain for the
-th view,
.
An array that groups all the available scalar measurements at the antenna ports for all the considered views can be defined as
The size of this array is equal to , that is, the total number of multiantenna measurements for all the possible pairs of transmitting and receiving elements.
The array , which is a collection of all the multistatic and multiview measurements of the scattered field at the time , represents the data used by the proposed method. In the practical implementation, the time-domain data are collected in a time window with in a set of discrete time steps separated by an interval . 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 . For the sake of simplicity, in the following formulation, we describe the processing approach focusing on a generic time instant .
Let us consider an LSTM cell with
as the input (
Figure 2). This cell is described by two recurrent vectors called
cell state and
hidden state , where
is the state size. The vectors
and
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 , , , and , 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 , , with the hidden state taken from the previous time step, and adds a bias vector. The output is then passed to a gate activation function (sigmoid or hyperbolic tangent). In particular:
where
,
, and
are real-valued parameters, and
is a sigmoid gate activation function, i.e., considering a vector
it is defined as
where
is the function used for gate activation, and related parameters are
,
,
.
with
,
and
.
with
,
, and
.
Based on the output of these four gates, the cell state
and hidden state
are, respectively, calculated as
where
indicates the Hadamard product. The hidden state
is then processed by a fully connected neural network composed of three layers, where the first two layers are composed of
neurons, and the last layer has
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:
where
,
,
,
. The final output is given by
with
,
.
The key point is that
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
where
is an array of the same size as
that contains the corresponding set of measurements of the scattered electric field
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)
-polarized fields and targets approximated as cylinders with uniform properties on the
axis.
With this definition of inputs (
) and outputs (
), the neural-network-based operator
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
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 (
, see
Figure 1). In more detail, a Fast Fourier Transform (FFT) algorithm is applied to the obtained approximated values of
to extract the frequency-domain scattered field at frequency
(i.e.,
frequencies). These values, indicated as
(an
time dependance is assumed and dropped), are used as input data by the inverse scattering algorithm, which assumes that incident electromagnetic fields are TM-
and that buried objects are infinite cylinders along the
axis, characterized by dielectric properties
and
, with
. With these assumptions, the
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):
where
represents the unknown of the inverse problem,
,
,
are normalization coefficients [
37],
,
is the half-space Green’s function [
41], and finally
is the
component of the total electric field inside the subsurface investigation domain. The latter quantity is also unknown and depends on
, making the problem nonlinear. Therefore, Equation (16) needs to be combined with the
state equation
where
is the
component of the incident electric field in the
-th view due to a line current source located at
.
Considering all the available measurement points and frequencies, the inverse problem is then formulated as a system of equations [
36] given by
This system can be written more compactly as
where
is the array that contains all the measured scattered fields,
is the unknown, and
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
-th step, calculates the Fréchet derivative of the operator
around a current tentative solution
(
), leading to the linearized equation
The unknown of the linearized problem at the
-th step is then
, which is found by a truncated Landweber-like regularization approach [
37], formulated by considering
as a variable exponent
space. The exponent
is an adaptive function computed as
in the first step (
is a parameter describing the exponent variation) and
in the following ones. Conversely,
is assumed as an
space,
, being equal to the mean value of
. After the solution of the linearized equation, the global solution is updated as
, 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
), and in the corresponding 2D canonical conditions (to obtain
).
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
) has been characterized as a homogeneous material with randomly chosen dielectric properties
and
. Random circular cylinders placed inside the ground with
-directed axes have been considered, characterized by dielectric properties
and
centered at coordinates
and
, with diameter
.
Scattered field data are collected by a set of
antennas above the soil with uniform spacing on an
-directed measurement line of length
at height
over the soil level. In more detail, in the
-th view, the transmitting antenna is centered at
with
, and the
receiving antennas are located at
In the 3D simulations, antennas have a bowtie shape with size , flare angle , are parallel to the ground and polarized on the axis. The FDTD simulation domain has , , and sides of length , , and , respectively. It has been discretized into cubic cells () with sides in length.
The vertical slice of the 3D simulation domain located on the plane at has also been the subject of 2D simulations, which have been carried out with a rectangular domain of sides and discretized into square cells () with sides in length. In the 2D case, transmitting antennas have been represented by infinite -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 has been simulated, with a time step of . Source antennas have been excited by a Ricker waveform with a central frequency of .
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., 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
.
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
) and the number of neurons present in each fully connected (FC) layer (i.e., the variable
). In particular, the values
and
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
where
is the reconstructed scattered-field measurement array obtained from the neural network, and
contains the corresponding actual values from 2D simulations. The obtained average of
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
, the error always decreases by increasing the state size up to
. Conversely, if the FC neurons are increased, an initial decrease in
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
and with
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 and , and a circular dielectric cylinder with , , diameter , and centered at , 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.,
) is shown in
Figure 5a (measurement point index reported in abscissa is computed as
). The corresponding quantity obtained from 2D simulations (
) and the result of the LSTM-based network applied to 3D data (
) are shown in
Figure 5b and
Figure 5c, respectively. The difference between the 2D canonical data and corrected 3D data,
, 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
.
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 frequencies equally spaced between and . Data have been corrupted with an additive white Gaussian noise with mean value equal to zero and a signal-to-noise ratio equal to . The investigation domain is a rectangular region on the plane between , , and , . The region has been discretized into square cells () with sides in length.
The inversion method has been executed with exponent function variation 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 .
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 (
) and below the target (
). Moreover, a barely visible ghost target is present close to the point
,
, 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 (
) and in the background (
) as
where the integrals are computed in the region occupied by the target (
) and in the background (
),
and
indicate the corresponding areas,
is the reconstructed value of the contrast function at the first frequency
(considered as reference), and
is its actual value.
The computed values of
and
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).