1. Introduction
The ability to investigate subsurface scenarios is needed for a variety of applications, including geophysical reconnaissance, archaeology, mine location and detection, environmental monitoring, and ground mapping in civil applications [
1,
2,
3,
4,
5,
6,
7,
8]. In this area, ground penetrating radars (GPRs) are one of the most widely adopted tools. However, GPR data are usually represented through B-scans, which are a time-domain representation of GPR measurements. Even though some post-processing techniques can be applied to these data, such as migration techniques [
1,
9], classical GPRs often require experienced users to interpret data. Moreover, they only provide a qualitative reconstruction of subsurface scenarios, without allowing for a precise dielectric characterization of the inspected targets [
6].
Among the various possible approaches to inspect underground structures, microwave imaging provides the ability to characterize subsoil regions based on scattering measurements collected by antennas placed over the area of interest. These techniques can be used to non-invasively determine the presence and the approximate shape of buried objects (in the case of qualitative methods [
10,
11,
12,
13,
14,
15,
16]) or the distribution of their dielectric properties (in the case of quantitative strategies [
17,
18,
19,
20,
21,
22,
23,
24,
25]).
The scenario in which prospecting is performed makes the imaging problem particularly challenging, especially focusing on quantitative inverse scattering methods. Indeed, in addition to the general issues of the inverse scattering problem, which is inherently ill-posed and nonlinear [
4], there are other issues that are typical of this configuration. First, unlike other applications where it is possible to place antennas around the entire survey area, here, antennas can only be placed on the survey line above the ground or in boreholes, reducing the available information to solve the problem [
26,
27]. In addition, the layering of the environment is another critical element; for instance, probes are positioned at a certain distance from the ground, and the presence of discontinuities in the dielectric properties of the scattering environment further adds a level of complexity due to the insurgence of reflection phenomena [
28].
In order to reconstruct accurate images of objects in the subsurface environments, imaging techniques that integrate electromagnetic models in these settings have been proposed in the literature. Indeed, the ability to accurately model the real environment in which the diagnosis is made is crucial to obtain good results [
29]. Moreover, the exploitation of a suitable model to solve the forward problem enables synthetic data to be obtained and testing methods in a simulated environment. In more detail, two main classes of mathematical models can be distinguished: analytical and numerical ones. The former are characterized by good robustness, but assuming simplified hypotheses [
30,
31,
32,
33]. The second class includes methods based on finite differences in the time domain, the finite element method, and the method of moments [
34,
35,
36]. These classes of approaches, although more computationally expensive, are particularly advantageous when certain geometric configurations and specific soils are present in the domain. In addition, several radar-like techniques [
37,
38,
39,
40] have been recently proposed in this research area to allow soil profiling, so that retrieved soil configuration can be accurately integrated into the electromagnetic model as a priori information [
41]. Finally, hybrid solutions combining analytical and numerical models have also been proposed [
11].
In this paper, a quantitative imaging method based on the inversion of scattering
-parameters combined with an electromagnetic model of finite element (FE) type is proposed and applied to the diagnosis of shallow subsurface structures. In the literature, several approaches have been described to solve the nonlinear inverse scattering problem by deterministic, stochastic, and neural networks techniques [
42,
43,
44,
45,
46,
47]. Among the deterministic approaches, an interesting class is represented by Newton-type schemes [
48,
49,
50,
51]. In the present work, the nonlinear and ill-posed inversion problem is treated using a Newton-type technique formulated in
Lebesgue spaces with variable exponent [
52,
53,
54]. This is an iterative regularization procedure which is able to provide accurate diagnostic results due to an adaptive definition of the exponent function. Moreover, the FE-based electromagnetic model combined with the formulation of
-parameters [
8] allows for a precise modelling of the electromagnetic problem. The approach, which was first developed for stroke imaging [
54], has now been expanded to buried target detection. The main plus of the proposed approach for shallow subsurface inspection is the capability of suitably taking into account the measurements and background configuration inside the inversion procedure. Indeed, in this way, the structure of the measurement system and environments can be considered in the electromagnetic model providing valuable reconstruction results even in this challenging scenario.
The method has been validated by means of numerical data obtained with a two-layer configuration illuminated by open waveguide probes. In detail, the variation in target size, depth, and some other relevant parameters has been evaluated. Finally, a study of the method in the case of non-planar air-soil interface has been conducted; the results of inversion with a priori knowledge of the exact interface have been compared with those achieved by replacing the actual surface profile with a planar one in the embedded model.
The paper is organized as follows: in
Section 2, the inverse scattering formulation along with the finite element approach and the inversion procedure are reported. The results are then presented in
Section 3. Finally, the Conclusions follow in
Section 4.
2. Formulation of the Problem
The shallow subsurface scenario in which detection is performed is shown in
Figure 1. A multistatic and multiview system with
waveguide antennas positioned at height
from air–soil interface
on a measurement line of length
parallel to
-axis is used to illuminate the soil. The air is modelled as vacuum (i.e.,
F/m) and the soil is characterized by a complex dielectric permittivity
.
The antennas are activated one by one, and measurements of transmission
S-parameters are completed between the antenna ports. In addition,
-polarized fields with no components of the electric field parallel to the
plane at an angular frequency
are considered, and the target properties are hypothesized invariant along the
axis. Therefore, a two-dimensional environment is assumed [
55].
The distribution of the dielectric permittivity of the region containing the cross-section of the objects under test is retrieved starting from scattering -parameters between antenna ports and , , where is the -parameter measured in the presence of the target and is the corresponding quantity related to an empty investigation domain (assumed to be available or estimated).
Moreover, the dielectric properties of the investigation domain
can be described by the contrast function
where
is the complex dielectric permittivity of the configuration under test and
is the background value. In a two-dimensional setting, the scattering parameters are related to the contrast function by means of the following data equation [
8]
where
is the incoming wave amplitude on
-th antenna port,
is the
-component of the incident electric field when the
-th antenna acts as a transmitter and
is the
-component of the total electric field when the
-th antenna transmits,
. In detail,
represents the field collected in the presence of the unknown target; therefore, it is a function of the unknown contrast function itself, and this results in a nonlinear relation linking the contrast function with the scattering parameters.
Considering the measurements collected by each antenna for each view, it follows:
This is the nonlinear system to be inverted starting from scattering parameters measurements to retrieve the contrast function.
To handle this problem, an FE approach is applied and integrated inside the inversion performed in variable exponent Lebesgue spaces. In detail, the electromagnetic model expressed through the FE method as applied to the present inverse scattering problem is described in
Section 2.1. Then, in
Section 2.2, the inversion approach is presented.
2.1. FE Approach
A fundamental step for the development of the proposed imaging method is to define a procedure to model the electromagnetic problem, exploiting an FE formulation to compute the electric fields and -parameters inside the inversion procedure.
In more detail, the analyzed measurement system, shown in
Figure 2, is composed of a set of
open waveguides of width
and length
, filled with material of complex dielectric permittivity
. The waveguides are terminated by PEC boundaries on the lateral sides,
, and a waveguide port on the top side,
In each antenna, a reference system
is defined at the waveguide port centered at
where
is the mutual distance between antenna waveguides.
Assuming the
-th waveguide excited with an incoming wave of the
fundamental mode, to compute the electric field, the Helmholtz equation solution is necessary along with boundary conditions in order to take into account the exact structure of the environment and measurement system:
where non-magnetic materials are considered. The following boundary conditions must be imposed: on
, Dirichlet conditions state that
; on waveguide ports, the Dirichlet and Neumann conditions are
and
with
where
, i.e., the
-component of the electric field tangent to the
-th port when the
-th port is in transmitting mode, is:
with
the Kronecker delta,
the amplitude of
-th mode in
th port,
the orthonormal modal function of
modes in the
th waveguide, and
the propagation constant of
modes inside the waveguide [
56].
The first term of Equation (6) represents the incoming
mode feeding the
-th port and the second one represents the outcoming
modes inside the
th waveguide. For the solution of the forward problem, coefficients
are left as unknowns, and the
-parameters are then retrieved as
. Moreover, in the analysed configuration, absorbing boundary conditions (ABC) are required to terminate the simulation domain outside the antenna ports. To this end, the perfectly matched anisotropic absorber (PMA) has been selected [
56,
57].
The simulation domain includes the above-described shallow subsurface scenario comprised air and soil layers and limited by waveguide ports, waveguide PEC walls, and the PMA layer. In order to solve the forward problem with FE formulation, the domain has been partitioned in a mesh of
triangles of dimension
with frontier
Then, the electric field
in each
-th element of the mesh is formulated through first-order basis functions. It is worth highlighting that such a triangular form of subdomains gives the possibility of mapping complex structures in an accurate way, making it suitable to describe the problem at hand. In particular, the first-order triangular basis functions for each
-th triangle,
, satisfy properties
and
, where
are the coordinates of the
-th node of each
-th triangular element and are defined as follows:
The electric field in each
-th triangle can be written as:
By considering the FE approximation of field along with the Helmholtz equation, the electric field is numerically computed following the approach in [
56].
Then, in the soil region, a rectangular investigation domain
of dimension
centered at
has been defined (
Figure 3), which, in its discrete representation, is a subregion of the simulation domain composed of
triangles. Therefore, we can introduce the vector
, which is obtained by approximating the dielectric properties as constant inside each triangle of the domain
, thus, the corresponding contrast vector results
.
In this way, starting from Equation (2), considering the FE formulation, we obtain the following system:
where
is the vector with the nodal values coefficients of the field in each triangle computed by means of the FE method and
is a
matrix of coefficients:
By means of Equation (8), a discrete nonlinear operator linking the contrast vector with the measurements in is defined.
2.2. Inversion Procedure
The imaging problem is solved by inverting Equation (9). In order to address the inversion, an inexact Newton procedure has been exploited where the space
of the unknowns is a variable exponent Lebesgue space
(i.e.,
) and the space
of the data is a space
with constant exponent (i.e.,
) [
53]. In particular, since the discrete case is considered, the exponent function of unknowns’ space
turns out to be a vector
,
being the value of the exponent in each
-th triangle. The constant exponent of data space
is set equal to its spatial average value in
, i.e.,
.
The inexact Newton method solves the inverse problem in a regularized way by minimizing the following residual function:
with
norm in the space
, by moving along a nonstandard gradient direction in the dual space
of
. In detail, the method consists of two nested iteration cycles as summarized in the flow chart in
Figure 4 [
54]. At first, at each
-th iteration of the external cycle, Equation (9) is linearized around its current solution (denoted as
) with null initial value
. Then, a Landweber-like approach in Lebesgue space
is adopted to solve the resulting linearized problem (as shown in the red box of the flowchart). With reference to
Figure 4,
,
, and
represent the duality maps of
, and
, respectively. (
is the dual space of
, with
as the Hölder conjugate of
),
is the Fréchet derivative of
, and
is the step width [
52,
54].
Moreover, the exponent function is adaptively modified during the inversion to achieve accurate reconstruction results. In the first iteration of the external cycle, the exponent vector is set to constant values
since no a priori information is available. Then, at each outer iteration, the current solution
is exploited for updating
, with the
-th component
where
denotes the range of
. This way, values of
close to
are assigned to the portions of
without targets leading to a reduction in ringing effects; values of
close to the maximum are in target regions where a smooth reconstruction is favored.
Finally, the two loops are stopped when the proper convergence criteria are fulfilled. In this case, the number of inner, , and outer iterations, , and the relative variation of the residual between two consequent steps, , are considered as stopping rules.
3. Numerical Validation
In this section, numerical validation of the proposed imaging method is presented. In particular, a shallow subsurface scenario has been simulated by the FE forward solver to produce synthetic data and a reference case is described in
Section 3.1. Then, the behavior of the method was investigated when some variations in the various parameters were applied. In particular, in
Section 3.2 and
Section 3.3, the effects and limitations on the performance of the method versus a variation in target size and target depth have been analysed, respectively. Then, in
Section 3.4, the method has been studied considering different geometries and the number of targets, while in
Section 3.5 the effect of the background uncertainties has been explored. Finally, in
Section 3.6, some numerical tests in the presence of a non-planar air–soil interface are reported; the results achieved with exact interface surface known a priori or replaced with planar surfaces in the inversion procedure have been compared.
3.1. Reference Case
In numerical simulations, synthetic S-parameter data were provided as input to the proposed method to determine the distribution of complex dielectric permittivity within the investigation area. The FE forward solver has been exploited for the generation of numerical data.
Concerning the measurement parameters, a set of
antennas was placed along the measurement line of length
at a distance
mm from eachother and located at
mm over a planar air–soil interface (
Figure 1). The waveguides are filled with a material with a dielectric permittivity
and the dimensions
mm in width and
mm in length. The soil is characterized as dry sand with a permittivity
.
Configurations with and without a target were simulated at frequency
and
mode is considered. The simulation domain is discretized by Gmsh [
58] using the frontal Delaunay algorithm with maximum edge length at the antenna ports and elsewhere with a size of
mm and
mm, respectively. In this way, total and incident
S-parameters have been generated. Total
S-parameters were corrupted with a multiplicative Gaussian noise of 3%.
Once synthetic data have been generated, a rectangular investigation domain of dimension and cm centered in cm is considered. Inside , a coarser discretization of edges mm is considered for the inversion.
The results have been evaluated using the following error metric:
where
is the analysed region inside the investigation domain composed of
elements, with
(i.e., whole domain
or target region
are considered),
is the reconstructed dielectric permittivity in the
-th element, and
is its reference value in the same element.
Initially, the target under consideration has a square cross-section of side , it is centered at and characterized by a dielectric permittivity . The method has been run under the following setting: , , and . Furthermore, solution loops are terminated when a threshold is reached. Concerning the range of variation in the exponent function, a study has been conducted to find its optimum value. Specifically, has been considered with step .
Figure 5 shows the errors in the whole investigation domain,
, and in the target region
. As can be noticed, concerning the error in the investigation domain, an increment versus
can be observed whereas the error inside the targets has a minimum in
. Therefore,
has been selected since it offers the best target reconstruction and high-quality results in the investigation domain. In
Figure 6, the real part of relative dielectric permittivity is reported along with the actual target shape. As can be noticed, the localization of the target is adequate, and its dielectric permittivity has been reconstructed quite well.
3.2. Variation in Target Size
The size of the target has been changed in the first set of numerical simulations. Specifically, the values
cm with a step of 3 cm have been considered for the side length (i.e.,
has been varied between
with
the wavelength in the soil).
Figure 7 shows the reconstructed distributions of the real part of the dielectric permittivity for
cm. In configurations with
cm and
cm [
Figure 7a,b], the reconstruction of the dielectric permittivity is close to the actual value. A good reconstruction is obtained, and the dielectric target is clearly visible and adequately localized. Conversely, considering a smaller target with
cm, the object is barely visible in the reconstruction [
Figure 7c]. In
Figure 8, the reconstruction errors computed in the investigation domain and inside the target region are reported. As can be noticed, the whole domain error
tends to increase with the target size. Instead, the error on the object
has a parabolic-like behavior with a minimum at
cm (
). In particular, very small targets are weak scatterers and become difficult to detect. Indeed, an underestimation of the dielectric properties can be observed by reducing the target’s size until the method is not able to reconstruct the object. On the contrary, slight degradation of the background reconstruction can be seen by increasing the size of the target [
Figure 7a].
In addition, the method was compared with the inversion formulated in the classical Hilbertian space (i.e., ). By comparing the errors achieved with the proposed variable exponent space method, it can be observed that the latter approach allows an accuracy improvement both in the investigation domain and in the target under test. This proves that the inversion performed in variable exponent Lebesgue spaces leads to more accurate reconstructions. Indeed, by defining a proper map of the exponent [Equation (12)], low values are assumed in the background and values close to are in the region where the inhomogeneities are localized. In this way, low values better control sparsity in the background whilst a better estimation of smoothness of the object is reached by assigning relatively high values of exponent in that region.
3.3. Variation in Target Depth
As a further analysis, the depth of the buried target has been varied, i.e.,
cm with a step of 2.5 cm (i.e.,
) has been considered in order to assess the effectiveness and limitations of the method depending on the depth. The reconstructed distributions of the real part of the relative dielectric permittivity,
, are shown in
Figure 9 for the cases
cm,
cm, and
cm. Moreover, in
Figure 10, the scattering
S-parameters simulated when the first antenna acts as the source are compared with those computed inside the inversion procedure from the reconstructed dielectric permittivity for cases
cm and
cm. As can be observed, the reconstructed data match with the simulated ones.
By comparing the results, as the depth increases [
Figure 6 and
Figure 9a,b], the target reconstruction shows a decrease in the accuracy of the dielectric properties up to the point where it is difficult to locate the object [
Figure 9c].
The reconstruction errors in the whole investigation domain and target region are shown in
Figure 11. For low values of depth, the best target and whole domain error are achieved; then, by deepening the buried target, both the error in the target region and in the whole investigation domain has an upward trend. In more detail, beyond the depth of
cm (
), the detection of the object is more difficult, and the quality of the reconstruction deteriorates as evidenced by the error trend, which for
cm shows
.
3.4. Effect of Different Geometries and a Different Number of Targets
In this section, the behavior of the method has been investigated considering different geometries and a different number of targets.
To this end, the target with a square-cross section of side cm analysed in the previous section has been considered and a second target with a circular-cross section has been introduced in the investigation domain. In detail, the second target is “vacuum-filled”, has a radius cm, and is located at .
Figure 12 shows the reconstructed distribution of the real part of the dielectric permittivity together with the actual shapes of the targets. As can be noticed, both targets with different geometries are correctly localized, and although in the case of the object with a circular-cross section a slight underestimation can be observed, a quite good reconstruction is achieved. Moreover, by comparing the reconstruction of the square-cross section target alone [
Figure 7b] and in the presence of a second object, a slight deterioration can be observed, due to the interaction between the objects. This is confirmed by the computed reconstruction errors in the target region and in the whole domain (i.e.,
and
), which are higher than in the single-target case (i.e.,
and
).
3.5. Effect of Uncertainties in Soil Dielectric Properties
In order to assess the robustness of the method versus an inexact knowledge of the dielectric properties of the soil, the background dielectric permittivity used inside the inversion procedure has been set to
. All the other parameters are the same as those described in
Section 3.1. In the first case, an underestimation of the background dielectric permittivity is considered whereas in the second case an overestimation is assumed.
Figure 13 shows the reconstructed distributions of the real part of the dielectric permittivity. The reconstruction errors in the whole investigation domain and in the target region are shown in
Table 1 (for the sake of comparison, we also reported the errors when the exact background properties are known in the inversion analysed in
Section 3.1). The error in the whole investigation domain increases when an inaccurate background permittivity is considered. Moreover, a decrement of the target error can be noticed in the first case because the initial estimation of background permittivity is closer to the properties of the object.
3.6. Variation in Air–Soil Interface Roughness
A further study has been conducted to inspect the robustness of the method in the presence of a non-planar air–soil interface. Moreover, the results have been compared when:
the model endowed in the inversion takes into account the exact interface profile (i.e., assuming that the surface is known a priori);
a planar interface is considered instead (i.e., no a priori information is available).
In the analysed measurement setting, the antennas are located at
cm over the average soil level. To simulate a non-planar interface, in this set of tests, the profile of
is modelled with Catmull–Rom splines [
59]. In particular,
is based on
equally spaced control points located to have a no-planar air-soil interface as shown in
Figure 14. The root mean square height of the interface has been varied between
and
.
A “vacuum-filled” target (i.e., characterized by a permittivity ) centered at with cm has been investigated. The investigation domain is centered at cm. All the other parameters of the forward and inverse procedure are the same as those in the previous sections.
In
Figure 15, the reconstruction errors in the target region and in whole investigation domain are reported.
Figure 16a–c shows the real part of the reconstructed dielectric permittivity when the surface profile is considered as a priori information in the inversion and
Figure 16d–f reports the reconstruction when planar interfaces are adopted in the FE-model inside the inversion procedure.
As can be noticed from the reconstruction of the real part of the permittivity, although the reconstruction is quite successful in both cases, the artifacts in the background appear more pronounced for greater interface roughness when no a priori information about the interface is assumed in the inversion. This is confirmed by the trend of the errors; the error on the object is comparable in the two cases for , while with , the knowledge of the exact interface brings a significant improvement in the reconstruction; the error on the whole domain is better in all cases with a valuable improvement as the roughness increases.
This proves the potentiality of the model embedded in the inversion procedure combined with S-parameters formulation in various operating conditions. Indeed, an accurate structural description of the involved environments can be integrated and taken into account in the inversion procedure, leading to an enhancement in reconstruction results.
4. Conclusions
In this paper, a microwave imaging technique for shallow subsurface prospecting is proposed, whose aim is to provide the distribution of dielectric permittivity from scattering S-parameters measurements. The method combines an FE approach with an inversion procedure in Lebesgue spaces with variable exponent. In this way, the structure of the measurement system together with the environment can be accurately considered together in the electromagnetic model by the FE formulation that is incorporated into the inversion procedure. Furthermore, the adopted inversion technique, which operates in non-Hilbertian Lebesgue spaces, exploits the adaptive update of the exponent function during the inversion procedure leading to good results even in this challenging scenario.
At first, the method is tested by considering a set of waveguide probes to illuminate a “vacuum-filled” target buried in a two-layered configuration with a planar interface. The behavior and limitations of the method such as the size and depth of the target change have been studied. Moreover, the effects of geometries and the number of targets as well as the uncertainties in soil dielectric permittivity have been investigated. In addition, the behavior of the method in the case of the non-planar air–soil interface has been analyzed, as well as the effect of the a priori knowledge of the interface profile inside the inversion procedure. The potentialities of this method for shallow subsurface prospection are shown by the achieved results. Future developments include both the validation with experimentally measured data and the extension to three-dimensional settings.