Microwave Imaging by Means of Lebesgue-Space Inversion: An Overview

: An overview of the recent advancements in the development of microwave imaging procedures based on the exploitation of the regularization theory in Lebesgue spaces is reported in this paper. Such inversion schemes have been found to provide accurate results in several microwave imaging scenarios, thanks to the di ﬀ erent geometrical properties that Lebesgue spaces can exhibit with respect to the more classical Hilbert ones. Moreover, the recent extension to the more general case of variable-exponent Lebesgue spaces is also addressed. Experimental results involving reference data are shown for supporting the theoretical description of the approaches.


Introduction
Microwave imaging is attracting an ever-growing interest in several applicative areas, ranging from non-destructive evaluation of civil and geophysical structures to biomedical diagnostics [1][2][3][4][5][6][7][8][9][10][11]. Such an interest is due mainly to the attracting features of microwave techniques. Actually, they are in principle able to provide a reconstruction of the physical and geometrical properties of the targets under test, e.g., the distributions of the dielectric properties, which cannot be directly obtained by other more consolidated systems. Moreover, microwaves are non-ionizing radiations and, consequently, since the power required for the inspection is usually quite low, they are safe for the operators and the inspected bodies. Finally, from a technological point of view, the imaging systems share many components with other apparatuses, such as the ones used in telecommunications, making it possible to use relatively cheap and off-the-shelf components.
Taking a mathematical point of view, it is required to "invert" the relationship that maps the dielectric properties of the inspected region with the measured electric field. Such a relationship is, in general, nonlinear and strongly ill-posed. Consequently, proper inversion procedures need to be devised, to consider both these theoretical problems. To this end, several approaches have been proposed in the last years [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26]. They have been investigated in the context of Hilbert spaces, in most cases. Furthermore, in this case it is possible to use the spectral theory to analyze the convergence and regularization properties of the solving method. However, despite this potential advantage that highly simplifies the mathematical study, regularization methods in Hilbert spaces in several cases lead to smoothed solutions, due to the filtering effects produced by regularization. This kind of behavior may be problematic in numerous imaging applications, especially when small targets needs to be retrieved. Compressive sensing strategies, which allow to retrieve sparse solutions, have been proposed recently to mitigate this drawback. More recently, regularization methods developed in the framework of the regularization theory in the more general Banach spaces also have been introduced in the mathematical literature [27][28][29] and investigated for microwave imaging applications [30][31][32][33].

Overview of the Inverse-Scattering Problem Formulation
Let us consider the imaging configuration sketched in Figure 1. An unknown target, characterized by a complex relative dielectric permittivity * ( ) = ( ) − ( )/ , being the position vector, the real part of the relative dielectric permittivity, and the (equivalent) electric conductivity, is located in a known investigation domain Ω. Briefly, the single-view case is described in the following, being the extension to multi-view multi-illumination configurations straightforward. The target is illuminated by an incident electric field , and the total electric field resulting from the interactions between the objects and the interrogating radiation is collected in a set of points , , = 1, … , , located in a predefined observation domain Γ. A timeharmonic dependence is assumed and the exp( ) term is omitted in the following. As it is well known, the relationship between the scattered electric field ( ) = ( ) − ( ) and the dielectric properties of the inspected region can be written in terms of the following integral equation: where ( ) = * ( )/ − 1, ∈ Ω, is the contrast function, is the relative dielectric permittivity of the background medium, which is assumed to be constant, = is the squared wavenumber in the background medium, and is the corresponding dyadic Green's function. It is worth As it is well known, the relationship between the scattered electric field E scatt (r) = E tot (r) − E inc (r) and the dielectric properties of the inspected region can be written in terms of the following integral equation: where c(r) = * r (r)/ b − 1, r ∈ Ω, is the contrast function, b is the relative dielectric permittivity of the background medium, which is assumed to be constant, k 2 b = ω 2 µ b is the squared wavenumber in the background medium, and = G is the corresponding dyadic Green's function. It is worth remarking that (1) is nonlinear with respect to the contrast function c, since the total electric field E r inside Ω itself depends on the dielectric properties of the inspected region and it is given by a relationship similar to (1) [7]. For sake of simplicity, a tomographic configuration is considered from now on. Particularly, it is assumed that the targets have cylindrical symmetry, i.e., c(r) = c(r t ), r t being the position vector in the transverse plane, and the incident electromagnetic field is transverse magnetic, i.e., E inc (r) = E inc (r t )ẑ. Under such hypotheses, only the z-component of the total electric field is different from zero and the scattering problem reduces to a two-dimensional and scalar one, i.e., where g is the two-dimensional Green's function for the background configuration. It is worth noting that E tot (r t ) is still unknown inside the two-dimensional cross-section Ω t and, thus, a second equation is needed, i.e.: Combining the data Equation (2) and the state Equation (3), the following nonlinear functional equation is obtained [7]: where I denotes the identity operator. Equation (4) can be formally written in compact form as: where x ∈ X is the unknown function to be retrieved (the contrast function inside the investigation domain Ω t ), y ∈ Y represents the data of the problem (the scattered electric field in the observation domain Γ), and F : X → Y is the nonlinear operator defined by (4).

Newton-Type Methods in Banach Spaces
The nonlinear functional Equation (5) is often solved by performing an iterative minimization of a residual (or cost) functional Φ : X → R , such as: where F : X → Y is a (Fréchet) differentiable operator between the two functional spaces X and Y. The basic minimization approach is the one-step gradient method, which is based on iterated linearizations of the nonlinear functional F .

Classical Mathematical Foundations in Hilbert Spaces
Taking the simplest case, that is, in the Hilbertian (i.e., Euclidean) spaces of the square integrable functions X = L 2 (Ω t ) and Y = L 2 (Γ), starting from a (suitable) initial guess x 0 ∈ X, the gradient method provides a sequence (x n ) n∈N ⊂ X that approximates a solution of (5) and is computed by the following one-step iterative scheme: x n+1 = x n − τ n p n , where τ n > 0 is a suitable step size and p n is the gradient of the residual functional Φ at x n , that is: Recalling that ∇Φ(x n ) is the direction of steepest increase of the residual functional Φ in a neighborhood of x n , the iterative step (7) gives rise to a reduction of Φ, provided that a suitable step size τ n > 0 is used. Moreover, since in any Hilbert space ∇ 1 2 u 2 Y = u, by using the chain rule for derivatives of composite functions, p n can be computed as follows: where F * x n : Y → X denotes the adjoint operator of the Fréchet derivative at point x n , indicated as F x n . The method (7) consequently reads as follows: where the step size has to be chosen to induce a decrease of the residual functional Φ, that is, to allow that Φ(x n+1 ) < Φ(x n ) [29]. When F : X → Y is linear, the Fréchet derivative coincides with F itself, that is, F x n ≡ F , and hence (10) reduces to the Landweber method when a constant step-size τ n ≡ τ ∈ 0, 2/ F 2 2 is fixed, or to the Steepest Descent method when the Cauchy optimal step size τ n = is used [35]. We remark that, by different choices of the ascent direction p n , we obtain other iterative minimization schemes, such us the conjugate gradient method [35], which is generally more powerful than the basic and simplest ones we briefly review here.
The one-step method (10) for the nonlinear Equation (5) is intrinsically related to the local first order linearization F x n of the nonlinear operator F , and it can be viewed as the simplest application of an (inexact) Newton scheme for non-linear equations. Indeed, let us consider the Taylor expansion of center x n and increment h n , that is F (x n + h n ) = F (x n ) + F x n h n + O h n 2 . Using the first order approximation, rather than solving the nonlinear functional Equation (5), we can consider the associated linear equation F (x n ) + F x n h n = y w.r.t. the unknown h n as follows: Using the Newton method, the least squares solution of the linear equation (11) gives rise to a new element: which, in general, is a better approximation of the solution x of the nonlinear Equation (5), and the procedure is then iterated according to a certain stopping rule. It is interesting to notice that, in real applications, the linear Equation (11) of the n-th Newton method is solved by means of iterative minimization methods. This way, the full solution algorithm is made of two nested iterative schemes, that is, the whole algorithm involves an outer-inner iterative procedure, since each outer Newton step is solved by means of a sequence of inner minimization steps. As a basic example, in the following we explicitly consider a Newton method with inner Landweber iterations, described in Algorithm 1.

Algorithm 1.
Two level (outer-inner iterations) inexact Newton method for the nonlinear Equation (5) (I) Let X and Y be two Hilbert spaces, and X 0 ∈ Y be an initial guess (X 0 = 0 is used when no a-priori information is available). Set the initial outer iteration index to n = 0. (II) OUTER STEP: Linearize (5) by means of the Fréchet derivative F x n at point x n and consider the associated linearized system (11), that is F x n h n = Y − F (x n ), with respect to the unknown h n . (III) INNER STEP: Find a (regularized) solution of the linear Equation (11) by means of an iterative minimization, with respect to h, of the n-th residual 1 2 F x n h − (y − F (x n )) 2 Y . Specifically, let h n,0 = 0 ∈ X be the inner initial guess. Then, for k = 0, 1, 2 . . ., compute: until a certain stopping rule is satisfied (e.g., a maximum number of inner iterations K LW is reached or the norm of the functional 1 2 F x n h n,k+1 − (y − F (x n )) 2 Y falls below a specified threshold). The obtained regularized solution of the n-th linear system (11) is denoted as h n,k .
(IV) Update the current (outer step) solution by setting: (V) IF a predefined stopping rule (e.g., based on the discrepancy principle [36]) on the outer iteration x n+1 is satisfied THEN return x n+1 (and STOP); ELSE continue with the subsequent outer iteration, by setting n = n + 1 and going to step II.
The outer-inner inexact Newton scheme (13)- (14) for the nonlinear functional equation (5) in the classical Hilbert space setting is a regularization algorithm that has been widely applied to nonlinear inverse scattering problems [37,38]. We recall that the Newton scheme is called "inexact" because each linear system is not solved exactly, but its solution is just iteratively approximated. It is interesting to notice that the outer-inner inexact Newton scheme (13)- (14) is an extension of the (single-step) Landweber iterative method (10). Indeed, if the inner iterations of (13) are always stopped at the very first iteration, that is, if the number of inner iteration is always fixed to k max ≡ 1 (which means that the linearized equation (11) is solved with a first, and very low, level of accuracy), then the two-steps scheme (13)- (14) coincides with the one-step scheme (10). Hence, it is quite evident that the method (13)- (14) overcomes (10) in both speed and quality, because at each Newton iteration (11), the associated linear equation is solved with a higher level of accuracy.

Extension to Banach Spaces
The extension of the Newton method to Banach spaces is not straightforward. Generally, given a linear operator L : X → Y between two Banach spaces X and Y, its adjoint operator L * acts between the associated dual spaces Y * and X * , that is, L * : Y * → X * . Generally speaking, we recall that the dual space B * of a Banach space B is the space of all the linear functionals from B to the real values, that is, B * = {b * : B → R, linear}. When the Banach space is a Hilbert space, by virtue of the Riesz representation theorem [27], given a linear functional b * ∈ B * , we can identify uniquely b * with the unique vector b ∈ B such that b * (z) = b, z , ∀z ∈ B, being ·, · the scalar product of the Hilbert space B. Hence, in this case, the dual space B * is isometrically isomorph to the space B. This is the reason why the adjoint operator L * in Hilbert spaces is always represented as L * : Y → X (rather than the formal correct definition L * : Y * → X * ), since both the isometric isomorphisms I Y : Y → Y * and I X : X * → X are implicitly applied.
Let us come back to the Landweber iterations (10) or (13), which are therein defined onto Hilbert spaces X and Y, so that F * x n : Y * → X * can be identified with F * x n : Y → X . Regarding this, (10) is well defined, since the operator F * x n : Y → X acts rightly on the residual F (x n ) − y ∈ Y. The same holds true for (13), since F x n h n,k − (y − F (x n )) ∈ Y. Moreover, the subtraction in (10) is performed between the two operands x n and τ n F * x n (F (x n ) − y), both belonging to X. The same applies for (13), where h n,k and τ n F * x n F x n h n,k − (y − F (x n )) belong to X. However, to extend the iterations (10) and (13) to Banach space setting, their forms need to be modified, since the terms F * x n (F (x n ) − y) of (10) as well as F * x n F x n h n,k − (y − F (x n )) of (13) are no longer correct, because now the adjoint operator The key tools for the generalization to Banach spaces are the so-called duality mappings [27]. Usually, a duality map is a special function that associates an element of a Banach space B with an element of its dual B * , and it is useful when B is not isomorph to B * . The duality map has an illustrative meaning in the context of minimization of convex functionals, as explained by the Asplund Theorem [27]. To this aim, given a convex functional f : B → R , we recall that the subdifferential of f is the multi-valued operator ∂ f : where, for any s ∈ B and s * ∈ B * , we have used the so-called pairing notation s * (s) = s * , s = s, s * . The subdifferential extends the concept of gradient to general Banach spaces. Indeed, if the convex functional is differentiable, then is b * ∈ B * is unique and can be identified as its gradient, since it holds is a duality map of B, which is then defined as follows: We recall that, if the Banach space is a Hilbert space, then J B 2 is simply identified with the identity operator, since Generally, this is not true in Banach spaces. However, thanks to the Asplund Theorem, from (16) and by using again the chain rule for the derivatives of composite functions as in (9), the subgradient of the residual functional 1 r F (x n ) − y r Y at point x n in Banach spaces can be computed explicitly as: We can notice that, differing from the well-known least square term F * x n (F (x n ) − y) of (9) in Hilbert space, we have now the term F * x n J Y r (F (x n ) − y), which is well defined in Banach spaces since F * x n : Y * → X * and J Y r (F (x n ) − y) ∈ Y * . Anyway, to obtain a well-defined generalization of the iterative step (10) in Banach spaces, we have to consider also that the addendum F * x n J Y r (F (x n ) − y) now belongs to X * , that is, the dual space of X. Hence, the sum has to be computed into the space X * as follows: where now x * n = J X r (x n ) ∈ X * . Subsequently, it is necessary to return to the original space X, i.e.: During this iteration, it is required that the space X is reflexive, that is, X * * is isomorph to X, so that J X * r * x * n+1 ⊆ X. Additionally, we point out that, in general, any duality map is a multi-valued map, and an arbitrary choice of a single element is implicitly assumed in both J Y r and J X * r * . Anyway, in our application, the Banach spaces are always L p (Ω t ) with 1 < p < +∞, and any duality map is always single-valued in these functional spaces [27].
The same arguments of (18) and (19) apply to the inexact Newton approach too, leading to a new version of the inner iterations (13) in Banach spaces, involving the duality maps J Y r and J X * r * . Specifically, Step III of Algorithm 1 now reads: with respect to h, of the n-th residual 1 2 F x n h − (y − F (x n )) 2 Y . Specifically, let h n,0 = 0 ∈ X be the inner initial guess. Then, for k = 0, 1, 2, . . ., compute: until a certain stopping rule is satisfied (e.g., a maximum number of inner iterations K LW is reached or the norm of the functional 1 2 F x n h n,k+1 − (y − F (x n )) 2 Y falls below a specified threshold). The obtained regularized solution of the n-th linear system (11) is denoted as h n,k .

The Role of the Exponent Parameter p in the L p Lebesgue Spaces Solution
The Landweber inner method (20) for the linearized outer system (11) is conceived in the framework of the regularization theory in Banach spaces. Considering our applicative case, X and Y are the Banach spaces of the Lebesgue p-summable functions, L p (Ω t ), with 1 < p < +∞. Any L p space with 1 < p < ∞ is reflexive, smooth and strictly convex [27], so that: consequently no other local minima arise in the context of L p regularization, and differentiable; (ii) the duality maps J Y r and J X * r * are always single-valued.
Based on this, recalling that the L p -norm is defined as: the duality map J L p r , which is the differential of 1 r u r L p , for 1 < p, r < ∞, can be explicitly computed as [28]: where the sign function is defined as sign(u) = u/|u| when u 0 and zero otherwise. If u ∈ L p (Ω t ) then J L p r (u) ∈ (L p ) * (Ω t ), which is isometrically isomorph to L p * (Ω t ), and the function on the right-hand side of (22) effectively belongs to L p * (Ω t ) [28], being p * the Hölder conjugate of p, that is, 1 p + 1 p * = 1. Moreover, if the Banach space is the Hilbert space L 2 , then J L 2 2 reduces to the identity operator, that is, J L 2 2 (u) = u, ∀u ∈ L 2 (Ω t ), as expected due to the isometric isomorphism between (L 2 ) * and L 2 .
Differing from the power parameter r, which acts merely as a scaling factor, the parameter p > 1 of the chosen L p space has a crucial meaning. Indeed in (22), the value of p gives rise to different amplifications of the small and large components of its argument u. As an example, let us consider a value 1 < p < 2: Then |u l | p−1 |u l | for small |u l | < 1, and |u l | p−1 |u l | for large |u l | > 1, which means that the duality map J L p r of (22) emphasizes the small components and reduces the large ones (and obviously the behavior is the opposite for p > 2). As a general comment, solving the functional Equation (5) in the framework of a L p space gives rise, just heuristically, to a Tikhonov-like regularization algorithm R α by considering α = p − 1, where α > 0 is the regularization parameter. This means that we obtain low regularization (i.e., low filtering and low smoothness), for small p > 1 close to 1, corresponding to a regularization parameter α close to 0, and high regularization (i.e., high filtering and high smoothness), for large p 1, corresponding to a regularization parameter α 0. Indeed, the numerical examples in the next Subsection briefly will show that, for p ≥ 2, some oversmoothing effects appear and discontinuities between different scattering media are not well reconstructed, as generally happens with too large Tikhonov regularization parameters α 0. Quite the opposite, with smaller p 1, the restoration of the discontinuities is more accurate, although some instability and noise amplification may arise, as usually encountered with a too small choice of the Tikhonov regularization parameter α 0.

A Reconstruction Example
An example of reconstruction obtained by applying the fixed-exponent Lebesgue space inversion procedure presented above is reported here. In particular, the FoamDielExtTM target from the reference measured data provided by the Institut Fresnel is considered [39]. Such a target is composed by two adjacent cylinders: The first one has center in (0, 0) cm, radius 4 cm, and relative dielectric permittivity 1.45, whereas the second one is centered in (−5.55, 0) cm, has radius 1.55 cm, and relative dielectric permittivity 3 (in both cases, the electric conductivity is negligible). The object is illuminated by horn antennas located in S = 8 positions uniformly spaced on a circumference of radius 1.67 m, and, for each view, the scattered field is collected in 241 points uniformly spaced on an arc of 270 • on the same circumference. Data acquired at different frequencies are also available in the range [2][3][4][5][6][7][8][9][10] GHz, with 1 GHz step. The details of the measurement setup can be found in [40]. During the inversion procedure, the assumed investigation region is a square domain of side 20 cm, which has been discretized into 63 × 63 square subdomains. The maximum number of outer and inner iterations have been fixed to N IN = 50 and K LW = 10, respectively. Moreover, to test the approach in optimal conditions, the loops are stopped when the variations of the normalized root mean square error NRMSE = c − c / c , c being the actual value of the contrast function, fall below 0.5%.
The reconstructed distribution of the relative dielectric permittivity for some of the available working frequencies are shown in Figure 2. The results obtained with the optimal value of the norm parameter p are reported in the first row, whereas the corresponding Hilbert-space reconstructions are provided for comparison in the second row. As can be seen, the Lebesgue-space method is able to reconstruct correctly the inspected scenario in all cases, providing a quite accurate reconstruction of the cylinders' cross section, both in terms of dimensions and dielectric permittivity. Using the Hilbert-space procedure it still is possible to identify both targets, but a significant over-smoothing effect is present (which prevents a good reconstruction of the small cylinder) and the ringing in the background is higher. This basic numerical example confirms the behavior of the inversion procedure in L p Lebesgue spaces, with respect to different choices of the exponent p parameter, as briefly discussed at the end of the previous Subsection. Indeed, in Figure 2d-f, where p = 2, oversmooting effects are quite evident, whilst in Figure 2a-c, where p = 1.2 and 1.3, the reconstructions are less oversmoothed (but some numerical instabilities may occur, especially for smaller p, not shown here). These two kinds of results usually are associated with too high and too low regularization in classical Tikhonov Hilbertian approaches, respectively.
Such considerations are also supported by the values of the mean relative errors reported in Table 1, defined as: where * r denotes the reconstructed dielectric permittivity, and Ω obj , Ω bg are the subdomains occupied by the target and by the background, respectively. discussed at the end of the previous Subsection. Indeed, in Figure 2d,e,f, where = 2, oversmooting effects are quite evident, whilst in Figure 2a,b,c, where = 1.2 and 1.3, the reconstructions are less oversmoothed (but some numerical instabilities may occur, especially for smaller , not shown here). These two kinds of results usually are associated with too high and too low regularization in classical Tikhonov Hilbertian approaches, respectively. Such considerations are also supported by the values of the mean relative errors reported in Table 1, defined as: where ̃ * denotes the reconstructed dielectric permittivity, and Ω , Ω are the subdomains occupied by the target and by the background, respectively.  Table 1 also reports the computational data for the considered cases (on a computer equipped with an Intel i5-8265u CPU and 16 GB of RAM). Specifically, the numbers of performed outer iterations, N * in , and the corresponding average computational time per iteration, t m , are provided. As expected, the time needed for performing a single outer iteration is similar between the optimal norm and the Hilbert-space approaches. However, in all cases, less iterations are needed when considering the optimal Lebesgue-space procedure, allowing a faster reconstruction.
It is worth remarking that the over-smoothing effect in the Hilbert-space solution can be reduced by varying the regularization parameter, which in the present approach is represented by the number of iterations. Table 2 reports the reconstruction errors obtained with different values of the parameters K LW in the case f = 2 GHz (the threshold on the NRMSE is not set). As can be seen, the object error decreases when higher values of K LW are used (i.e., a lower regularization is performed), becoming comparable with the ones provided by the optimal value of the norm parameter (in particular, the peak value of the dielectric permittivity is closer). However, the background error increases significantly, producing a greater overall reconstruction error. 0.14 0.15 0.14 100 0.14 0.14 0.14 As an example, Figure 3 shows the mean relative reconstruction errors versus the norm parameter p for the case of the lowest frequency (i.e., f = 2 GHz). The overall reconstruction error e inv presents a minimum value corresponding to the optimal norm parameter p opt . Following that, the error increases monotonically with p. Moreover, as expected, the background error is always increasing, since low values of p produce sparser solutions. Concerning the object error, in this case a minimum is present at p opt , and after an initial increase it becomes almost constant. Similar trends can be observed for the other frequencies.
To summarize, Figure 4 shows the behavior of the residual functional ( Figure 4a) and of the normalized root mean square error (Figure 4b) versus the outer iteration number for the case f = 2 GHz. As can be seen, for all values of p the algorithm converges after few iterations (between 5 and 8).

Multifrequency Lebesgue-Space Inversion
To increase the available information, multi-frequency data can be exploited. It is assumed to have at our disposal the scattered-field data collected for a set of values of the angular frequencies , = 1, … , . Subsequently, a subscript is added to the frequency-dependent functions and operators to specify at which frequency they refer. However, since the contrast function depends upon the frequency, it is necessary to modify the problem formulation [31,41]. To explain, let us assume that the dielectric permittivity and the electric conductivity do not depend on the frequency (i.e., dispersion is neglected). The contrast function ( ) in this case can be rewritten as: where is a reference frequency and the new unknown, which does not depend on the frequency, is given by . The inverse scattering problem at a single frequency can be then described by the following equation: where ℱ is given by (4).

Multifrequency Lebesgue-Space Inversion
To increase the available information, multi-frequency data can be exploited. It is assumed to have at our disposal the scattered-field data collected for a set of values of the angular frequencies , = 1, … , . Subsequently, a subscript is added to the frequency-dependent functions and operators to specify at which frequency they refer. However, since the contrast function depends upon the frequency, it is necessary to modify the problem formulation [31,41]. To explain, let us assume that the dielectric permittivity and the electric conductivity do not depend on the frequency (i.e., dispersion is neglected). The contrast function ( ) in this case can be rewritten as: where is a reference frequency and the new unknown, which does not depend on the frequency, is given by . The inverse scattering problem at a single frequency can be then described by the following equation: where ℱ is given by (4).

Multifrequency Lebesgue-Space Inversion
To increase the available information, multi-frequency data can be exploited. It is assumed to have at our disposal the scattered-field data collected for a set of F values of the angular frequencies ω f , f = 1, . . . , F. Subsequently, a subscript ω is added to the frequency-dependent functions and operators to specify at which frequency they refer. However, since the contrast function depends upon the frequency, it is necessary to modify the problem formulation [31,41]. To explain, let us assume that the dielectric permittivity and the electric conductivity do not depend on the frequency (i.e., dispersion is neglected). The contrast function c ω (r t ) in this case can be rewritten as: where ω 0 is a reference frequency and the new unknown, which does not depend on the frequency, is given by x MF . The inverse scattering problem at a single frequency ω can be then described by the following equation: where F ω is given by (4). It is worth noting that now x MF is a real-valued function, whereas the electric fields are complex-valued. To work with all real-valued functions, the data are split into their real ( ) and imaginary ( ) parts. Applying (25) to all the available F frequencies and stacking all the resulting functional equations, the multi-frequency operator equation to be inverted can be then formally written as: Concerning the Fréchet derivative of the new operator, needed in the outer linearization step, it can be written in terms of the derivative of the complex single-frequency operator as: It is worth remarking that the multifrequency inversion technique is based on the use of the Lebesgue-space theory discussed in Section 3. Consequently, the same value of the norm parameter p is used for all the elements of the data vector y MF (i.e., the scattered fields at the different frequencies). A possible future extension of the approach would be to use different values of p for the different frequencies (i.e., for different parts of the data vector) by exploiting the variable-exponent Lebesgue space theory discussed in Section 5.

A Reconstruction Example
To show the effectiveness of the multi-frequency approach combined with the Lebesgue-space inversion procedure, the FoamDielExtTM target considered in Section 3.4 has been used again as a reference target. As previously described, the Fresnel data are available in the range 2-10 GHz. Consequently, the F frequencies considered in the inversion are f i = 2 + (i − 1) GHz, i = 1, . . . , F. Figure 5 shows the results obtained by considering F = 2, F = 4, and F = 7. As can be seen, increasing the number of processed data allows one to improve the reconstruction quality with respect to the single-frequency inversion procedure, especially concerning the edges of the two cylinders. This improvement is also confirmed by the reconstruction errors reported in Table 3. As expected, all the reconstruction errors decrease when F increases, although for F higher than about 5 a slower improvement can be observed. Even in this case, the use of a norm with an exponent lower than 2 always produces better reconstructions than the corresponding standard Hilbert-space approach. However, when a high number of frequencies are used, the advantages of using lower exponent parameters become less significant (although always present). It is worth remarking that, in this case, the Hilbert-space approach could provide results comparable to the optimal Lebesgue-space procedure by considering a higher number of frequencies. As an example, the reconstruction errors obtained with the Hilbert-space approach with 7 frequencies are comparable to the ones of the optimal Lebesgue-space procedure with 4 frequencies. However, the computational cost is significantly higher (the needed time per iteration is doubled and the number of iterations is three times higher). Concerning the computational times, it should also be noted that in the multifrequency case they are higher than the single-frequency case and increase with the number of processed frequencies. Such an increase is due to the higher dimensions of the involved matrices. Moreover, in the current implementation, several sub-matrix and sub-vector extractions are performed for building the Fréchet derivative in (27)-(28), leading to an additional computational burden.

The Novel Variable-Exponent Approach
To summarize, the inverse scattering results obtained by using the inexact Newton method as a regularization algorithm for the nonlinear equation (5) in some Lebesgue spaces show less presence of over-smoothing in the reconstructed properties of the targets. Naturally, this is an interesting property for all the applications in which an accurate dielectric reconstruction is required. Despite these advantages, a critical issue should be tackled, related to the selection of the value of that gives rise to the best results. This parameter is essential, since it allows the definition of the space in which the algorithm operates and, thus, the corresponding norm. Confirmed by heuristic analyses, it turns out that the background of a dielectric reconstruction obtained by an -space inversion with low values of the exponent > 1 is cleaner from ringing effects. Additionally, discontinuities in the dielectric properties of the targets are retrieved with a higher level of accuracy. Conversely, values of 2 are more suitable for reconstructing the internal part of scatterers with homogeneous dielectric characteristics. Taking the point of view of applications, a sort of compromise between these trends has to be found. Some indications have been derived from numerical and experimental analyses, in which this parameter has been studied with regard to different target typologies, their dielectric properties, size, and the amount of data signal-to-noise ratio [30,41]. Despite this, with the tools described so far, only an a-posteriori selection of the optimal has been done. Clearly, such an optimal choice is feasible only in controlled conditions, when the actual solution is known.
To overcome such a limitation, an innovative approach working in Lebesgue spaces (⋅) (i.e., where the exponent is a non-constant function) has been recently devised in [42]. The key point of this method is the possibility to have a variable parameter (⋅) in the investigation domain, which is not limited anymore to be a unique constant value. The considered functional spaces are different from the Lebesgue spaces of -integrable functions (in which is constant) and the power turns out to be a function of the position, in other words. The result is the freedom of defining a

The Novel Variable-Exponent Approach
To summarize, the inverse scattering results obtained by using the inexact Newton method as a regularization algorithm for the nonlinear Equation (5) in some L p Lebesgue spaces show less presence of over-smoothing in the reconstructed properties of the targets. Naturally, this is an interesting property for all the applications in which an accurate dielectric reconstruction is required. Despite these advantages, a critical issue should be tackled, related to the selection of the value of p that gives rise to the best results. This parameter is essential, since it allows the definition of the L p space in which the algorithm operates and, thus, the corresponding norm. Confirmed by heuristic analyses, it turns out that the background of a dielectric reconstruction obtained by an L p -space inversion with low values of the exponent p > 1 is cleaner from ringing effects. Additionally, discontinuities in the dielectric properties of the targets are retrieved with a higher level of accuracy. Conversely, values of p ≈ 2 are more suitable for reconstructing the internal part of scatterers with homogeneous dielectric characteristics. Taking the point of view of applications, a sort of compromise between these trends has to be found. Some indications have been derived from numerical and experimental analyses, in which this parameter has been studied with regard to different target typologies, their dielectric properties, size, and the amount of data signal-to-noise ratio [30,41]. Despite this, with the tools described so far, only an a-posteriori selection of the optimal p has been done. Clearly, such an optimal choice is feasible only in controlled conditions, when the actual solution is known.
To overcome such a limitation, an innovative approach working in Lebesgue spaces L p(·) (i.e., where the exponent is a non-constant function) has been recently devised in [42]. The key point of this method is the possibility to have a variable parameter p(·) in the investigation domain, which is not limited anymore to be a unique constant value. The considered functional spaces are different from the Lebesgue spaces of p-integrable functions L p (in which p is constant) and the power p turns out to be a function of the position, in other words. The result is the freedom of defining a different value of p for each point inside the region under inspection: This way, the sparsity of the background can be enforced with low values of p, whereas values of p closer to 2 can be used to achieve a more accurate reconstruction of the properties inside the scatterers.
The price to pay is a more complex mathematical structure of the Lebesgue spaces with variable exponents L p(·) compared to the Lebesgue spaces L p (i.e., with constant exponent). An example is the computation of the L p(·) -space norm. First, to calculate the norm, the evaluation of the modulus ρ p(·) is necessary, i.e.: which generalizes the argument of the (constant-value) p-root of (21) to an exponent function 1 < p(r t ) < +∞, r t ∈ Ω t [42]. The p-root is required in the definition (21) for guaranteeing that the homogeneity property of any norm holds, that is αv L p = |α| v L p for any real value α. Moreover, the p-root of the constant case cannot be directly extended to the non-constant case, because there is not a unique value to use. The L p(·) -norm, which extends the L p -space definition (21) and also is known as a Luxemburg norm [42], is found, thus, by solving a one-dimensional minimization problem, that is: When the exponent p(r t ) is a constant value (denoted asp), it clearly results that: where the identity ρ p(·) v γ = 1 γ p ρp(v) holds only in this constant case (otherwise it cannot be formally written).
Considering the above concepts, the extension of the Landweber-type iterations in (20) to the variable-exponent Lebesgue space X = L p(·) relies on the calculation of the L p(·) -space duality maps and has the following form: where J L p(·) r and J L q(·) r * are the duality maps of Y = L p(·) and L q(·) , and the function q(·) represents the Hölder conjugate of p(·), defined as 1 p(r t ) + 1 q(r t ) = 1. This iterative algorithm, although formally similar to (20), has relevant differences inside the expression of the duality map J L p(·) r (v), which is now defined as: This expression has been found by exploiting the Fréchet derivative of the norm in L p(·) spaces [43] and needs to evaluate the Luxemburg norm (30) of v. Taking a purely theoretical point of view, it is important to note that, differing from the conventional constant exponent case, the duality map J L q(·) r * is only an approximation of the respective duality map of the dual space of L p(·) which is used in the iterative scheme.

Strategies for Choosing the Variable Exponent Function
The possibility of arbitrarily selecting the exponent p(·), which is a function with values 1 < p(·) < +∞, can be exploited to apply different exponent parameters to different regions of the investigation domain. However, the choice of such a function p(·) has not been discussed yet in this survey. The basic idea behind our proposal has been explained in the previous subsection: Exponents p close to 1 are useful in the background (to improve the sparsity and reduce the ringing phenomena), while exponents p close to 2 are useful inside the scattering objects (to estimate the true values of their dielectric distribution). Although background and location of the objects are not known as prior information, but rather their retrieval is the aim of the proposed algorithms, an estimate can be obtained by the first iterations of any classical methodology. This way, we can consider the first reconstruction (even if really poor) as basic information for constructing the exponent function p(·) for the space of unknown X. Concerning the data space Y, the norm parameter p is kept constant and equal to the average value of p(r t ), r t ∈ Ω t .
Based on this paradigm, we have developed an adaptive and automatic iterative procedure for defining the nonconstant exponent p(·). An improved accuracy, as well as a greater stability compared to reconstruction algorithms in constant-exponent L p spaces, characterize the novel procedure. To detail more, a map of the values of p(r t ), for r t ∈ Ω t , is updated at each Newton iteration on the basis of the retrieved values of the contrast function inside the investigation domain Ω t . The exponent function is chosen such that p min ≤ p(r t ) ≤ p max (with p min > 1 and p max > p min ), where p min and p max are two values fixed before starting the inversion process, which identify the range of variation. The update equation of the exponent p(r t ) is given by: (34) in which n denotes the current inexact-Newton iteration. The mapping function Υ(·) has values from [0, 1] to [0, 1] and presents a monotonic increase. The application of (34) results in a non-linear scaling of the values of the contrast function magnitude (divided by its maximum value) in the given interval range for the exponents, [p min , p max ]. Since the background region (i.e., the set of points of Ω t where no targets are present) is characterized by low values of |x n |, the corresponding values of p(r t ) given by (34) will approach the minimum of the range p min . In this way, the region outside targets would present a reduction of background ringing and artifacts. Conversely, in the regions internal to the scattering targets (where |x n | is close to its maximum) p(r t ) will approach p max , providing a good retrieval of the properties of large regions without overestimation peaks. Equation (34) defines p(r t ) by exploiting the reconstructed contrast function, which is typically available after the first outer step of the inexact-Newton loop. To start the inexact-Newton iterations, different strategies can be followed based on the availability of a-priori information. When an initial guess is available before executing the inversion (that is, we know a starting approximate solution x 0 0), (34) can still be applied considering such an estimated value of the unknown. Alternatively, when this kind of information is not available, a constant function p(·) = p start can be considered, obtaining the same initial iteration as the fixed-exponent approach (20).

A Reconstruction Example
To give an example, the FoamDielExtTM target considered in Section 3.4 has been reconstructed by using the variable-exponent approach. The parameters of the outer and inner loops are the same as in the constant exponent case. The power function p(·) has been defined with p min = 1.4 and p max = 2. Five different typologies of mapping functions are experimented:

5.
Sine → Υ(u) = 1 2 sin πu − π 2 + 1 2 The reconstruction results are summarized in Table 4, which reports the reconstruction errors obtained with the variable-exponent procedure when considering the data at 2, 3 and 4 GHz. Regarding this configuration, the different mapping functions provide comparable errors. Moreover, the obtained results are similar to the ones provided by the optimal norm parameter in the fixed-exponent algorithm. The corresponding computational times are reported in Table 5. As can be seen, the time needed for performing a single outer iteration is only slightly higher than the corresponding time for the fixed exponent case (see Table 1). Such a small increase is due to the increased complexity in the computation of the Luxemburg norm. However, the number of needed iterations is still comparable to the fixed-exponent case, thus the overall computational burden is only slightly greater. It is worth remarking that, although the reconstruction quality (both in terms of computational cost and errors) is comparable, the main advantages of the proposed variable L p approach are related to the fact that it does not require an a-priori selection of the norm parameter p, whose optimal value can be selected only a-posteriori and with some knowledge of the expected results. Actually, the user is only required to select a range of values for this parameter, whose choice is far less critical. Figure 6 shows the reconstructed distribution of the relative dielectric permittivity for a linear mapping function.
Observing the results, we notice a correct retrieval of the properties of both targets. Additionally, the L p(·) -space reconstructions look very similar to the ones obtained when considering the fixed-exponent method with an optimal choice of the exponent p. This is further confirmed by the cuts of the relative dielectric permittivity along the x axis reported in Figure 7.

Conclusions
An overview of the recent advancements in the development of microwave imaging methods based on inversion procedures in Lebesgue spaces has been reported. First, an introductory description of the inverse-scattering problem formulation has been provided. Subsequently, the application of Newton-based techniques combined with regularization strategies in Hilbert and fixed-exponent Lebesgue spaces has been reviewed, including multi-frequency approaches. Next, the more recent variable-exponent space methods have been detailed. The adaptive procedure used to define, iteration by iteration, the power function of the variable-exponent spaces removes the need for choosing a fixed space exponent value. Previously, the choice of the "best" fixed exponent was possible only with an a-posteriori error analysis, based on known configurations. Therefore, the applicability of this kind of inversion methods in real environments is now significantly increased, in particular when no a-priori information about the target is available. The considered inversion techniques have been analyzed and compared by considering experimental data, confirming their effectiveness also under realistic operating conditions. Author Contributions: Conceptualization, C.E., A.F., M.P., and A.R.; methodology, C.E., A.F., M.P., and A.R.; formal analysis, C.E., A.F., M.P., and A.R.; investigation, C.E., A.F., M.P., and A.R.; validation, C.E., A.F., M.P., and A.R.

Conflicts of Interest:
The authors declare no conflict of interest.

Conclusions
An overview of the recent advancements in the development of microwave imaging methods based on inversion procedures in Lebesgue spaces has been reported. First, an introductory description of the inverse-scattering problem formulation has been provided. Subsequently, the application of Newton-based techniques combined with regularization strategies in Hilbert and fixed-exponent Lebesgue spaces has been reviewed, including multi-frequency approaches. Next, the more recent variable-exponent space methods have been detailed. The adaptive procedure used to define, iteration by iteration, the power function of the variable-exponent spaces removes the need for choosing a fixed L p space exponent value. Previously, the choice of the "best" fixed exponent was possible only with an a-posteriori error analysis, based on known configurations. Therefore, the applicability of this kind of inversion methods in real environments is now significantly increased, in particular when no a-priori information about the target is available. The considered inversion techniques have been analyzed and compared by considering experimental data, confirming their effectiveness also under realistic operating conditions.

Conflicts of Interest:
The authors declare no conflict of interest.