Microwave-Based Subsurface Characterization through a Combined Finite Element and Variable Exponent Spaces Technique

A microwave characterization technique to inspect subsurface scenarios is proposed and numerically assessed in this paper. The approach is based on a combination of finite element electromagnetic modeling and an inversion procedure in Lebesgue spaces with variable exponents. The former allows for description of the measurement system and subsurface scenario with high accuracy, while the latter exploits the adaptive definition of exponent function to achieve improved results in the regularized solution of the inverse scattering problem. The method has been assessed with numerical simulations regarding two-layered environments with both planar and non-planar air–soil interfaces. The results show the capabilities of the method of detecting buried objects in different operative conditions.


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

Formulation of the Problem
The shallow subsurface scenario in which detection is performed is shown in Figure 1. A multistatic and multiview system with W waveguide antennas positioned at height h from air-soil interface Λ on a measurement line of length l m parallel to x-axis is used to illuminate the soil. The air is modelled as vacuum (i.e., ε 0 8.85 × 10 −12 F/m) and the soil is characterized by a complex dielectric permittivity ε b .
illuminate the soil. The air is modelled as vacuum (i.e., ≃ 8.85 10 F/m) and the soil is characterized by a complex dielectric permittivity .
The antennas are activated one by one, and measurements of transmission Sparameters 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 crosssection 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] − 2 ( , ) ( , ) ( , ) = 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 thecomponent 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: The antennas are activated one by one, and measurements of transmission S-parameters are completed between the antenna ports. In addition, z-polarized fields with no components of the electric field parallel to the xy plane at an angular frequency ω are considered, and the target properties are hypothesized invariant along the z axis. Therefore, a twodimensional environment is assumed [55].
The distribution of the dielectric permittivity of the region R containing the crosssection of the objects under test is retrieved starting from scattering S-parameters between antenna ports i and w, δŜ iw = S iw tot − S iw inc , where S iw tot is the S-parameter measured in the presence of the target and S iw inc is the corresponding quantity related to an empty investigation domain (assumed to be available or estimated).
Moreover, the dielectric properties of the investigation domain R can be described by the contrast function where ε(x, y) is the complex dielectric permittivity of the configuration under test and ε b 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 c m is the incoming wave amplitude on m-th antenna port, E i inc is the z-component of the incident electric field when the i-th antenna acts as a transmitter and E w tot is the z-component of the total electric field when the w-th antenna transmits, i = w. In detail, E w tot 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.

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 S-parameters inside the inversion procedure.
In more detail, the analyzed measurement system, shown in Figure 2, is composed of a set of W open waveguides of width a and length b, filled with material of complex dielectric permittivity ε wg . The waveguides are terminated by PEC boundaries on the lateral sides, Π PEC , and a waveguide port on the top side, Π w . In each antenna, a reference system x (w) , y (w) is defined at the waveguide port centered at where s is the mutual distance between antenna waveguides. 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 , / . 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 Assuming the i-th waveguide excited with an incoming wave of the TE 10 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 Π PEC , Dirichlet conditions state that E i inc/tot Π PEC = 0; on waveguide ports, wg,inc/tot , i.e., the z-component of the electric field tangent to the w-th port when the i-th port is in transmitting mode, is: with δ iw the Kronecker delta, A (i,w) m,inc/tot the amplitude of m-th mode in wth port, e (w) m the orthonormal modal function of TE m0 modes in the wth waveguide, and β m the propagation constant of TE m0 modes inside the waveguide [56].
The first term of Equation (6) represents the incoming TE 10 mode feeding the i-th port and the second one represents the outcoming TE m0 modes inside the wth waveguide. For the solution of the forward problem, coefficients A (i,w) 1,inc/tot are left as unknowns, and the S-parameters are then retrieved as S iw inc/tot = A (i,w) 1,inc/tot . 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 N triangles of dimension Σ n with frontier Π (n) . Then, the electric field E i inc/tot (x, y) in each n-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 n-th are the coordinates of the u-th node of each n-th triangular element and are defined as follows: The electric field in each n-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 R of dimension l R × h R centered at (x R , y R ) has been defined ( Figure 3), which, in its discrete representation, is a subregion of the simulation domain composed of N R triangles. Therefore, we can introduce the vector ε ∼ = ε (1) , . . . , ε (N R ) T , which is obtained by approximating the dielectric properties as constant inside each triangle of the domain R, thus, the corresponding contrast vector results τ ∼ = τ (1) , . . . , τ (N R ) T . In this way, starting from Equation (2), considering the FE formulation, we obtain the following system: T is the vector with the nodal values coefficients of the field in each triangle computed by means of the FE method and [T (n) ] is a 3 × 3 matrix of coefficients: By means of Equation (8), a discrete nonlinear operator D(τ) linking the contrast vector with the measurements in δŜ is defined.  In this way, starting from Equation (2), considering the FE formulation, we obtain the following system:

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 X of the unknowns is a variable exponent Lebesgue space L p(·) (i.e., τ ∈ X ⊆ L p(·) ) and the space Y of the data is a space L p a with constant exponent (i.e., δŜ ∈ Y ⊆ L p a ) [53]. In particular, since the discrete case is considered, the exponent function of unknowns' space X turns out to be a vector p = [p (1) , . . . ,p (N R ) T , p (n) being the value of the exponent in each n-th triangle. The constant exponent of data space Y is set equal to its spatial average value in R, i.e., p a = 1 The inexact Newton method solves the inverse problem in a regularized way by minimizing the following residual function: with · Y norm in the space Y, by moving along a nonstandard gradient direction in the dual space X * of X. In detail, the method consists of two nested iteration cycles as summarized in the flow chart in Figure 4 [54]. At first, at each k-th iteration of the external cycle, Equation (9) is linearized around its current solution (denoted as τ k ) with null initial value τ 0 = 0. Then, a Landweber-like approach in Lebesgue space L p(·) is adopted to solve the resulting linearized problem (as shown in the red box of the flowchart). With reference to Figure 4, J X , J X * , and J Y represent the duality maps of X, X * , and Y, respectively. (X * = L p(·) * is the dual space of X, with p(·) * as the Hölder conjugate of p(·)), D k is the Fréchet derivative of D, and α = D k −2 is the step width [52,54]. 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 = 0. 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 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 p 0 = [p start , . . . ,p start ] T since no a priori information is available. Then, at each outer iteration, the current solution τ k is exploited for updating p, with the l-th component where [p min , p MAX ] denotes the range of p. This way, values of p close to p min are assigned to the portions of R without targets leading to a reduction in ringing effects; values of p 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, I max , and outer iterations, K max , and the relative variation of the residual between two consequent steps, ∆ R , are considered as stopping rules.

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 Sections 3.2 and 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.

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 W = 10 antennas was placed along the measurement line of length l = 1.095 m at a distance s = 55 mm from eachother and located at h = 50 mm over a planar air-soil interface ( Figure 1). The waveguides are filled with a material with a dielectric permittivity ε wg = 25ε 0 and the dimensions a = 60 mm in width and b = 80 mm in length. The soil is characterized as dry sand with a permittivity Configurations with and without a target were simulated at frequency f = 550 MHz and M = 1 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 s wg = 2 mm and s q = 4 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 R of dimension l R = 60 cm and h R = 50 cm centered in (x R , y R ) = (0, −30) cm is considered. Inside R, a coarser discretization of edges s inv = 7 mm is considered for the inversion.
The results have been evaluated using the following error metric: where Reg is the analysed region inside the investigation domain composed of N Reg elements, with Reg = {R, tar} (i.e., whole domain R or target region tar are considered), ε (n) is the reconstructed dielectric permittivity in the n-th element, and ε (n) is its reference value in the same element. Initially, the target under consideration has a square cross-section of side l t = 15 cm, it is centered at (x t , y t ) = (−15, −15) cm and characterized by a dielectric permittivity ε t = ε 0 . The method has been run under the following setting: K max = I max = 100, p start = p min , and p MAX = 2. Furthermore, solution loops are terminated when a threshold ∆ R = 1% is reached. Concerning the range of variation in the exponent function, a study has been conducted to find its optimum value. Specifically, p min = [1.2, 1.9] has been considered with step 0.1. Figure 5 shows the errors in the whole investigation domain, e R , and in the target region e tar . As can be noticed, concerning the error in the investigation domain, an increment versus p min can be observed whereas the error inside the targets has a minimum in p min = 1.3. Therefore, p min = 1.3 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.
Sensors 2023, 23, x FOR PEER REVIEW 9 of 19 tivity = . The method has been run under the following setting: = = 100, = , and = 2. Furthermore, solution loops are terminated when a threshold Δ ℛ = 1% is reached. Concerning the range of variation in the exponent function, a study has been conducted to find its optimum value. Specifically, = [1.2,1.9] has been considered with step 0.1. 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 = 1.3. Therefore, = 1.3 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.

Variation in Target Size
The size of the target has been changed in the first set of numerical simulations. Specifically, the values ∈ [3,18] cm with a step of 3 cm have been considered for the side length (i.e., has been varied between [0.095 , 0.572 ] with the wavelength in the soil). Figure 7 shows the reconstructed distributions of the real part of the dielectric permittivity for = {3,12,18} cm. In configurations with = 18 cm and = 12 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 = 3 cm, the object is barely visible in the reconstruction [ Figure 7c]. In Figure 8, the reconstruction errors , and = 2. Furthermore, solution loops are terminated when a threshold Δ ℛ = 1% is reached. Concerning the range of variation in the exponent function, a study has been conducted to find its optimum value. Specifically, = [1.2,1.9] has been considered with step 0.1. 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 = 1.3. Therefore, = 1.3 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.

Variation in Target Size
The size of the target has been changed in the first set of numerical simulations. Specifically, the values ∈ [3,18] cm with a step of 3 cm have been considered for the side length (i.e., has been varied between [0.095 , 0.572 ] with the wavelength in the soil). Figure 7 shows the reconstructed distributions of the real part of the dielectric permittivity for = {3,12,18} cm. In configurations with = 18 cm and = 12 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 = 3 cm, the object is barely visible in the reconstruction [ Figure 7c]. In Figure 8, the reconstruction errors

Variation in Target Size
The size of the target has been changed in the first set of numerical simulations. Specifically, the values l t ∈ [3,18] cm with a step of 3 cm have been considered for the side length (i.e., l t has been varied between [0.095λ b , 0.572λ b ] with λ b the wavelength in the soil). Figure 7 shows the reconstructed distributions of the real part of the dielectric permittivity for l t = {3, 12, 18} cm. In configurations with l t = 18 cm and l t = 12 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 l t = 3 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 e R tends to increase with the target size. Instead, the error on the object e tar has a parabolic-like behavior with a minimum at l t = 9 cm (0.286λ b ). 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]. 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 = 9 cm (0.286 ). 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].   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 = 9 cm (0.286 ). 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., X ⊆ L 2 , Y ⊆ L 2 ). 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 p MAX 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.

Variation in Target Depth
As a further analysis, the depth of the buried target has been varied, i.e., y t ∈ [−15, −35] cm with a step of 2.5 cm (i.e., y t ∈ [−0.477λ b , −1.112λ b ]) 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, Re{ε r }, are shown in Figure 9 for the cases y t = −32.5 cm, y t = −25 cm, and y t = −20 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 y t = −25 cm and y t = −20 cm. As can be observed, the reconstructed data match with the simulated ones. 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.

Variation in Target Depth
As a further analysis, the depth of the buried target has been varied, i.e., ∈ [−15, −35] cm with a step of 2.5 cm (i.e., ∈ [−0.477 , −1.112 ]) 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 = −32.5 cm, = −25 cm, and = −20 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 = −25 cm and = −20 cm. As can be observed, the reconstructed data match with the simulated ones.  By comparing the results, as the depth increases [ Figures 6 and 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 y t = −25 cm (y t = −0.794λ b ), the detection of the object is more difficult, and the quality of the reconstruction deteriorates as evidenced by the error trend, which for y t < −25 cm shows e R > 0.6. By comparing the results, as the depth increases [ Figures 6 and 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 = −25 cm ( = −0.794 ), the detection of the object is more difficult, and the quality of the reconstruction deteriorates as evidenced by the error trend, which for −25 cm shows 0.6.

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 = 12 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 "vacuumfilled", has a radius = 6 cm, and is located at ( , ) = (17, −20) cm. 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., =  By comparing the results, as the depth increases [ Figures 6 and 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 = −25 cm ( = −0.794 ), the detection of the object is more difficult, and the quality of the reconstruction deteriorates as evidenced by the error trend, which for −25 cm shows 0.6.

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 = 12 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 "vacuumfilled", has a radius = 6 cm, and is located at ( , ) = (17, −20) cm. 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., = Figure 11. Reconstruction errors in the whole domain and inside the target versus the target size, y t .

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 l t = 12 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 r = 6 cm, and is located at (x t2 , y t2 ) = (17, −20) cm. 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., e tar = 0.553 and e R = 0.074), which are higher than in the single-target case (i.e., e tar = 0.521 and e R = 0.035).

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 ε * r,b = ε r,b ± 0.5. 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.

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 , * = , 0.5. 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.

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 , * = , 0.5. 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.

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 h = 7 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 C = 7 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 h rms = λ 10 and h rms = λ 30 . 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 ℎ = 7 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 = 7 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 ℎ = 10 and ℎ = 30 .
A "vacuum-filled" target (i.e., characterized by a permittivity = ) centered at ( , ) = ( − 20, − 20) cm with = 12 cm has been investigated. The investigation domain is centered at ( , ) = (0, −32) 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.  A "vacuum-filled" target (i.e., characterized by a permittivity ε t = ε 0 ) centered at (x t , y t ) = (−20, −20) cm with l t = 12 cm has been investigated. The investigation domain is centered at (x R , y R ) = (0, −32) 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 h rms ≤ λ 20 , while with h rms = λ 10 , 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. 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 ℎ 20 , while with ℎ = 10 , 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.

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. 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.

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.

Data Availability Statement:
The numerical data presented in this study are available from the corresponding authors on request.