Abstract
An increasing number of practical applications of three-dimensional microwave imaging require accurate and efficient inversion techniques. In this context, a full-wave 3D inverse-scattering method, aimed at characterizing dielectric targets, is described in this paper. In particular, the inversion approach has a Newton-based structure, in which the internal linear solver is a conjugate gradient-like algorithm in spaces. The presented results, which include the inversion of both numerical and experimental scattered-field data obtained in the presence of homogeneous and inhomogeneous targets, validate the reconstruction capabilities of the proposed technique.
1. Introduction
The possibility of obtaining a point-by-point characterization of dielectric structures from noninvasive measurements of the scattered electromagnetic field has attracted the attention of scientists and researchers for decades [1–10]. Basically, the quantitative retrieval of the dielectric properties of a structure under test from scattered-field data is associated with the solution of an inverse problem [3], which is nonlinear (if no simplifying model approximations are made) and typically ill-posed. Several two- and three-dimensional approaches have been devised in the past years for the numerical solution of such a problem. Under the two-dimensional approximation, which assumes a particular form of the electromagnetic field (e.g., transverse-magnetic) and that the targets are invariant along one axis, the problem at hand can be significantly simplified [3]. Based on this 2D assumption, many inverse scattering techniques have been proposed [11–16] and subsequently exploited in experimental systems [17–19]. Among them, the family of the so-called Newton-based deterministic methods appears really valuable in the reconstruction process, providing accurate characterization results. The common aspect of such iterative techniques is that the original nonlinear problem is firstly linearized around the current solution estimate and then approximately solved by means of a linear regularization method. For their attractive features, Newton-based techniques have been chosen and specialized for different applications, ranging from the detection of buried objects [20, 21] and the imaging of civil structures and wood [19] to the biomedical uses for breast imaging [22, 23] and brain stroke detection [24].
Clearly, the development of both measurement systems and reconstruction techniques useful to address the full three-dimensional electromagnetic problem is more challenging. First of all, it is necessary to develop effective measurement systems able to acquire the needed field samples. In this framework, different setups have been proposed in the literature [25–28]. Secondly, the development of 3D imaging techniques should be guided by a responsible use of computational resources, a problem which unfortunately does not disappear even with the enormous increase in the computer power nowadays available. However, thanks to the relevant advances carried out by several research groups all around the world, the practical 3D imaging of complex structures does not seem an unreachable goal anymore [29–33].
Along with satisfying and encouraging results, the development of microwave imaging methods opens a series of theoretical and computational issues that have been only partially addressed. In particular, from a theoretical perspective, a critical point of the most deterministic approaches is the need for a priori information, which is sometimes crucial for a successful inversion with complex configurations [34, 35]. Another key issue is the accuracy of the provided solution, frequently affected by ringing phenomena and oversmoothing effects when the reconstruction procedure is formulated in the standard mathematical framework, that is, in the framework of Hilbert spaces. In order to overcome the latter problem, compressive sensing and sparsity-promoting techniques [36–39], as well as Banach space formulations, have been proposed, the last with both Landweber- [40, 41] and conjugate gradient-based inner solvers [42]. As regards Banach space methods and in particular space approaches, the lack of a dot product inducing a complete space for does not allow to use spectral tools such as the singular value decomposition, and therefore, much more involving convex analysis tools need to be used to characterize the inversion procedures. Anyway, despite such difficulties, the conjugate gradient (CG) method has been extensively analysed in the framework of minimization into spaces, and its convergence and regularization properties have been recently proven [43]. We highlight that the generalization of the CG method to spaces is nontrivial, by both a theoretical and a practical point of view. The main difficulty of applying the CG method to spaces is that, differing from the Hilbertian case, an explicit formula for computing the step size does not exist. The reason of this difference is that the derivative of the cost function for the computation of its stationary points leads to a linear operator, so that we can explicitly compute its solutions. Conversely, the derivative of the cost function leads to a nonlinear operator, and an explicit and closed-form formula for its solution does not exist. On this ground, the computation of the optimal step size unfortunately requires a one-dimensional iterative minimization procedure in the case. On the other hand, the CG method usually provides a faster convergence speed, especially when compared with the basic Landweber one.
Despite the increased complexity of the space approaches, the obtained results in both 2D [41, 42] and 3D configurations [44] are promising. So far, only the Landweber solver has been extensively validated for the use inside space inexact Newton techniques, but its convergence speed sometimes appears very low, leading to a considerable number of required iterations of the method. Such a behaviour may be problematic, especially in the 3D case where the computational cost for each iteration is quite significant. Consequently, in this paper, the Newton conjugate gradient approach, which has been found to provide a higher convergence rate in 2D microwave imaging applications [42], is extended for the first time to deal with three-dimensional configurations.
This paper has the following structure: the next section provides the basic mathematical background useful for understanding the proposed Newton-CG inversion method in spaces. The subsequent sections report a simulated and experimental validation campaign aimed at testing the inversion strategy in different operating conditions. Finally, conclusions are outlined, as well as some future goals.
2. Mathematical Formulation
In the considered microwave imaging problem, an unknown target is located in an investigation volume of known geometry, as schematically shown in Figure 1. The background medium is characterized by a complex dielectric permittivity . The inspected scenario is illuminated by a set of known time-harmonic incident electric fields , where denotes the position vector and the index of the considered illumination. The term , being the angular frequency, is omitted in the following. The total electric field , resulting from the interactions between the incident field and the object, is measured in a given observation domain , which surrounds the investigation volume.
As it is well known, the scattered electric field can be expressed in terms of the space-dependent dielectric properties in by means of the following integral relationship (data equation): where is the wavenumber in the background medium (assumed to be nonmagnetic, i.e., characterized by magnetic permeability equal to the one of the vacuum, ), is the dyadic Green’s function for free space [45], with being the three-dimensional scalar Green’s function and being the dyadic identity, and is the contrast function defined as
A similar relationship holds true for the total electric field inside the investigation volume (state equation), i.e.,
In order to numerically solve the scattering problem, equations (1) and (3) are discretized by using the method of moments with pulse basis functions and Dirac’s delta testing functions [46]. In particular, the investigation domain is subdivided into nonoverlapping cubic subdomains with center , , and side , whereas the field is sampled, for every view, in points , . In this way, the following discrete equations (corresponding to (1) and (3), respectively) are obtained: where the involved numerical arrays are defined as follows:
In (6)–(8) the subscripts // denote the components, with respect to a rectangular coordinate system, of the corresponding electric field vector. Moreover, in (4), the matrices and contain the contributions of the integrals of the Green’s dyadic function over the subdomains [3], which are analytically computed by using the relationship in [47], whereas denotes the identity matrix. Finally, the function is defined as where is a null matrix of dimension and transforms the numerical array in a diagonal matrix whose diagonal elements are the values contained in .
The discretized data and state equations in (4) are combined together to obtain the following nonlinear relationship:
Moreover, in order to exploit the information of the considered multiview setup, the numerical arrays for each illumination are stacked as follows: where is the discrete multiview operator that describes the mapping between the contrast function values and the measured samples of the scattered electric field .
In this paper, (11) is inverted by using a Newton scheme [40, 48] with a CG inner loop performing a regularization in the framework of the discrete spaces [43]. In particular, the unknown discrete numerical array is assumed to belong to the discrete Banach space , whereas the data array belongs to the linear space . The developed inversion procedure is schematized in Figure 2. Fundamentally, the algorithm is based on two nested loops. The outer one iteratively performs a linearization of (11) around the current estimate of the discrete distribution of dielectric properties described by the array (the subscript denotes the outer iteration index). This loop is initialized with a starting guess , which may incorporate the eventually available a priori information about the target. If such information is not available, a void volume is assumed, i.e., The linearization is performed by computing a first-order Taylor expansion of the discrete operator that involves the calculation of the Jacobian matrix at , which, for the 3D microwave imaging problem at hand, is defined as [3] where and are matrices and arrays containing the contributions of the integrals of the inhomogeneous dyadic Green’s function and of the internal total electric field due to the dielectric distribution described by . Such quantities are obtained as [3] and are efficiently computed by using the BiCGSTAB-FFT method [49].
The inner loop is used to solve the linearized equation, and it is based on a CG-like scheme performing a regularization in the framework of the discrete Banach spaces [50]. Such inversion procedure has been firstly proposed for microwave imaging of 2D structures in [42], and it is extended for the first time to 3D settings in this paper. The main advantage of the regularization in the more general Banach spaces is that they allow, for , to reduce the oversmoothing and ringing effects that are usually associated to the conventional reconstructions in Hilbert spaces [40, 44]. The key point for extending the standard CG inversion scheme to Banach spaces is represented by the duality maps , , and that are used in the update formulas in Figure 2 [51]. In particular, for the considered Banach space , by virtue of the Asplund theorem, the duality map is defined as
Analogous definitions hold for the duality map of data space and for the one of the dual space of , (where the Hölder conjugate of , i.e., , is used instead of ) [51]. Moreover, the parameters and are computed by using the following relationships: where denotes the Hermitian transpose of the Jacobian matrix and the minimization problem in (15) is solved by using the secant method [52].
It is important to remark that the main difference with respect to the approach presented in [44] is related to the adoption of the CG inner solver in each Newton linearization step instead of the basic Landweber algorithm. Such a modification allows to exploit the faster convergence and good regularization properties of the CG method, which still hold in spaces [43]. However, with respect to the Landweber algorithm, the extension to settings of the classical CG algorithm poses additional difficulties, since there is not any explicit formula for computing the optimal step size (which exists in the standard case). On this ground, we choose to apply a 1D iterative minimization approach based on the secant method to solve the optimization problem in (15).
3. Numerical Results
An initial validation campaign of the proposed full-wave inversion scheme has been carried out with simulated data. In the considered configuration, the background medium is vacuum (characterized by a dielectric permittivity , magnetic permeability , and wavenumber ). The working angular frequency is , with . The target is sequentially illuminated by -polarized plane waves with unit amplitude, whose directions of propagation are shown in Figure 3 (red arrows). For each view, the total electric field is collected in measurement points uniformly distributed on a sphere of radius , being the wavelength in a free space (blue dots in Figure 3). The number of measurement points has been chosen slightly larger than the minimum value needed to correctly sample the field, as suggested by the degrees of freedom theory [53]. Since the computational requirements (memory and computational time) increase almost linearly with the number of views, a limited number of illumination directions have been instead adopted. However, this value has been empirically found to be able to provide good results for the considered configurations. The investigation volume is a cube of side , which has been subdivided into cubic voxels inside the inversion procedure. The simulated data have been computed by using a custom method-of-moment numerical code, in which is partitioned into cubic subdomains. A signal-to-noise ratio on the total electric field data has been obtained by adding to the simulated values a white Gaussian noise with a zero mean value.
The iterative loops of the inversion procedure have been terminated when the relative variation of the data residual is below the predefined thresholds 0.01 (for the outer loop) and 0.05 (for the inner loop) or when the maximum numbers of iterations and are reached in the outer and inner loops, respectively. The accuracy of the dielectric reconstructions has been evaluated by using the normalized root mean squared error on the solution (), the average relative errors on the contrast function of the target () and of the background (), defined as where is the reconstructed value of the contrast function in the th subdomain in which the investigation volume has been partitioned, is the corresponding actual value, are the sets of indexes of the cells inside the target and background regions, respectively, and are the corresponding numbers of contained cells. It is interesting to note that the normalized root mean squared error represents a good compromise between the errors on the background and on the target. Such a trade-off is useful for characterizing the optimal value of the norm parameter in several scenarios, although it does not have an immediate physical interpretation. Meanwhile, the average relative errors and allow to obtain a more practical measure of the error on the two subparts of the domain occupied by the target and background.
3.1. Single Dielectric Sphere
First of all, the reconstruction capabilities of the proposed approach have been assessed by considering a single dielectric sphere centered at , characterized by a relative dielectric permittivity , and with radius . The reconstructed distributions of the relative dielectric permittivity obtained for with space exponents and (Hilbert space) are reported in Figure 4. The optimal value of the parameter has been selected on the basis of the normalized mean squared error, which is shown in Figure 5 versus for different values of the sphere radius . Considering both reconstructions and errors on the solution, it is clear that the proposed space method is able to outperform the classic Hilbert space approach for . Furthermore, the optimal value of monotonically increases with the target size. For comparison purposes, an inexact Newton approach with Landweber iterations in spaces has been also applied to the present case, with the same stopping criteria. The behaviours of the normalized data residuals in the inner and outer loops versus the iteration number are reported in Figure 6 for both methods, with and for the optimal exponent parameters ( and for the CG- and Landweber-based methods, respectively). A significantly faster minimization of the residual, associated with a lower number of both inner and outer iterations required to meet the convergence criteria, is observed with the proposed CG-based approach. Concerning the computational times, each CG iteration took in average about 4.5 s, whereas the Landweber one about 0.5 s (on a PC equipped with an Intel® Core™ i5 CPU and 8 GB of RAM). The significantly larger time needed by the CG method is mainly related to the solution of the minimization problem in (15), which is solved by a secant method in the present implementation. However, more advanced line-search procedures can be adopted in order to reduce the computational burden associated to the calculation of the step size (such as the derivative-free one used in [43]). Nevertheless, the overall Newton scheme combined with the CG solver requires a lower number of outer steps (5 in this case) to reach the final solution than a corresponding method with the Landweber inner loop (for which 9 iterations are needed).
(a)
(b)
The outer steps (which are the same in both methods) are by far the most computational expensive, since at every iteration it is necessary to update the Fréchet derivative and to compute the current estimate of the scattered field. In the considered case, such operations took in average about 200 s. Consequently, the overall computational time of the Newton scheme with the Landweber solver is significantly larger than the one obtained by using the CG method.
3.2. Separate Scatterers with Different Properties
We consider now the reconstruction of two separate dielectric scatterers. The first target is a dielectric sphere centered at , with radius , and relative dielectric permittivity . The second target is a -directed dielectric cylinder with circular cross section, characterized by radius , height , center , and relative dielectric permittivity . The reconstructed distributions of the relative dielectric permittivity for with and (Hilbert space) are reported in Figures 7 and 8, respectively. The targets are visible in both cases, but the reconstruction obtained in the space is more accurate and allows us to appreciate the difference in the dielectric permittivity of the two objects. The behaviour of , reported versus the parameter in Figure 9 for different values of , clearly indicates that the overall best results are achieved with . The average relative errors on the target and background regions and the numbers of performed Newton iterations , are reported in Table 1 for and . It is worth noting that an increase in the dielectric permittivity of the sphere corresponds to a slight decrease in the reconstruction accuracy, accompanied with an increase in the number of required Newton iterations.
(a)
(b)
(c)
(a)
(b)
(c)
However, in all cases, the proposed method gives the best results. Similar to the previous section, an example of the behaviours of the normalized data residuals is reported in Figure 10. As can be seen, in this configuration, the number of inner iterations is comparable, although the CG method still requires a slightly lower number of steps in the first outer iteration.
3.3. Inhomogeneous Targets with Different SNRs
To further validate the approach, an inhomogeneous structure has been considered. It consists of a -directed circular dielectric cylinder with radius , height , and centered at . The lowest half of the cylinder is characterized by a relative dielectric permittivity , whereas the upper one has . In this case, we also consider a variation in the signal-to-noise ratio . The normalized mean squared error and the values of for each SNR are listed in Table 2, where the average relative errors and number of performed outer iterations are also reported. It can be seen that increases, as expected, while decreases (i.e., less iterations are executed) with a decrease of SNR. The reconstructed distribution of the relative dielectric permittivity for () is shown in Figure 11. In the same figure, the corresponding reconstruction obtained in Hilbert spaces is reported, too. As highlighted by the errors in Table 2, the advantages of using the space procedure are evident in all the presented cases. Finally, the behaviours of the normalized data residuals are shown in Figure 12 for and considering the optimal value of the parameter (1.5 and 1.3 for the CG- and LW-based approaches, respectively). In this case, too, the inner CG solver requires less iterations, and also the overall Newton scheme stops earlier.
(a)
(b)
4. Experimental Results
Finally, the proposed approach has been tested against experimental data. The dataset provided by the Institut Fresnel [54] has been used, since it represents a standard and well-accepted benchmark for testing inverse-scattering procedures. In particular, the TwoSpheres and TwoCubes targets have been considered. In this case, the operating frequency is , and the measurement configuration is composed of sources on a sphere with radius , and measurement points on a circumference with the same radius. More details about the experimental setup can be found in [54]. The investigation volume is a cubic region with side centered at the origin, partitioned into voxels for the inverse problem solution. The trends of the reconstruction errors versus the parameter for both the considered targets are shown in Figure 13, whereas the reconstructions for and are reported in Figures 14 and 15. These results further confirm the capabilities of the CG-based inexact Newton approach in spaces, with , to accurately retrieve a quantitative reconstruction of the scene under test, outperforming the standard method in Hilbert spaces.
(a)
(b)
(a)
(b)
(a)
(b)
In particular, as shown in Figure 13, the background error increases as grows. Such a behaviour can be related to the fact that low values of usually enhance the sparsity of the solution, by also reducing the ringing effects and artefacts on the background (for which the values of the contrast function are zero) [40, 44]. On the contrary, the object is usually better retrieved with higher values of , although increasing too much, such a parameter produces an oversmoothing of the solution and thus a worse reconstruction of localized objects. Consequently, the optimal value of is a compromise between these two opposite behaviours, and in several cases, it depends upon the target size.
Finally, it is worth noticing that the reconstruction results compare very well with the ones obtained by using the Landweber-based procedure in [44]. In fact, the optimal value of the norm parameter is about the same (considering the same stopping criteria adopted for the CG procedure). Moreover, the shape, positions, and dielectric properties of the targets are correctly identified with the new approach, too. The reconstruction errors for the optimal value of the norm parameter are reported in Table 3. For comparison purposes, the corresponding values obtained by using the Landweber-based inversion procedure in [44] are also provided. As can be seen, the two approaches give rise to comparable results, although the normalized root mean square errors for the CG-based approach are slightly larger. Such a difference can be ascribed to the small variations in the background error, which is however related to a large part of the investigation domain. Nevertheless, the relative errors on the reconstruction of the objects are equal or lower for the CG-based approach. Moreover, as observed with simulated data, a lower number of outer iterations are required to satisfy the adopted stopping criteria.
5. Conclusions
A microwave imaging approach for reconstructing the three-dimensional distribution of the dielectric properties of unknown targets has been presented in this paper. The developed approach is based on the use of a Newton scheme, in which at each step the linearized problem is solved by using a conjugate gradient-like algorithm performing a regularization in the framework of the spaces. Such an approach, initially developed by the present authors for solving two-dimensional inverse scattering problems, has been extended for the first time to three-dimensional settings by taking into account the full-vector nature of the measured samples of the fields. The reported numerical and experimental results show that the proposed inversion procedure is able to effectively reconstruct the considered targets, by also providing less oversmoothing and ringing effects compared to Hilbert space regularization schemes. Moreover, with respect to previous approaches in which the inner loop was solved by using a Landweber-type algorithm, the new method presents a faster convergence.
However, a relevant remark should be devoted to the selection of the optimal exponent value . From the observed results, the choice of this parameter is critical for obtaining reliable results, and the optimal value depends on the problem and target configuration. No unequivocal rules give general criteria for this selection, and therefore, the best values have been selected here a posteriori, knowing the actual configuration and evaluating an error metric. It has been recently found that Lebesgue space approaches with nonconstant exponents may be attractive in dealing with such an issue [55]. However, so far, only 2D formulations have been proposed. Therefore, the use of these nonconventional spaces will be hopefully explored even in 3D settings as a further development, based on the findings reported in the present work.
Moreover, further developments will be also devoted to the validation of the approach in more complex configurations, where the wanted targets are embedded in different background media or hosting objects (e.g., in biomedical applications as well as in ground penetrating radar imaging).
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.