# Microwave Imaging of 3D Dielectric Structures by Means of a Newton-CG Method inlp Spaces

Authors: Claudio Estatico; Alessandro Fedeli; Matteo Pastorino; Andrea Randazzo; Emanuele Tavanti
Journal: International Journal of Antennas and Propagation
 (2019)

---

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


---

## Body

## 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 lp space approaches, the lack of a dot product inducing a complete space for p≠2 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 lp spaces, and its convergence and regularization properties have been recently proven [43]. We highlight that the generalization of the CG method to lp spaces is nontrivial, by both a theoretical and a practical point of view. The main difficulty of applying the CG method to lp spaces is that, differing from the Hilbertian l2 case, an explicit formula for computing the step size does not exist. The reason of this difference is that the derivative of the l2 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 lp 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 lp 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 thelp 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 lp 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 inlp 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 volumeVinv of known geometry, as schematically shown in Figure 1. The background medium is characterized by a complex dielectric permittivity ϵb. The inspected scenario is illuminated by a set of known time-harmonic incident electric fields Eivr, where r denotes the position vector and v the index of the considered illumination. The term ejωt, ω being the angular frequency, is omitted in the following. The total electric field Etvr, resulting from the interactions between the incident field and the object, is measured in a given observation domain Dobsv, which surrounds the investigation volume.Figure 1
Schematic representation of the considered imaging configuration.As it is well known, the scattered electric fieldEsvr=Etvr−Eivr can be expressed in terms of the space-dependent dielectric properties in Vinv by means of the following integral relationship (data equation):
(1)Esvr=−kb2∫Vinvcr′Etvr′·G¯¯br,r′dr′,r∈Dobsv,where kb=ωϵbμ0 is the wavenumber in the background medium (assumed to be nonmagnetic, i.e., characterized by magnetic permeability equal to the one of the vacuum, μ0), G¯¯br,r′=I¯¯+kb−2∇∇gbr,r′ is the dyadic Green’s function for free space [45], with gbr,r′=−e−jkbr−r′/4πr−r′ being the three-dimensional scalar Green’s function and I¯¯ being the dyadic identity, and c is the contrast function defined as
(2)cr=ϵrϵb−1.A similar relationship holds true for the total electric field inside the investigation volume (state equation), i.e.,
(3)Eivr=Etvr+kb2∫Vinvcr′Etvr′·G¯¯br,r′dr′,r∈Vinv.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 N nonoverlapping cubic subdomains with center rninv, n=1,…,N, and side d, whereas the field is sampled, for every view, in M points rmobs, m=1,…,M. In this way, the following discrete equations (corresponding to (1) and (3), respectively) are obtained:
(4)e_sv=Gobsvdiag3c_e_tv,e_iv=I−Ginvdiag3c_e_tv,where the involved numerical arrays are defined as follows:
(5)c_=cr1inv,…,crNinvt,(6)e_sv=Es,xvr1obs,…,Es,xvrMobs,Es,yvr1obs,…,Es,yvrMobs,Es,zvr1obs,…,Es,zvrMobst,(7)e_iv=Ei,xvr1inv,…,Ei,xvrNinv,Ei,yvr1inv,…,Ei,yvrNinv,Ei,zvr1inv,…,Ei,zvrNinvt,(8)e_tv=Et,xvr1inv,…,Et,xvrNinv,Et,yvr1inv,…,Et,yvrNinv,Et,zvr1inv,…,Et,zvrNinvt.In (6)–(8) the subscripts x/y/z denote the components, with respect to a rectangular coordinate system, of the corresponding electric field vector. Moreover, in (4), the matrices Gobsv and Ginv 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 I denotes the identity matrix. Finally, the diag3 function is defined as
(9)diag3c_=diagc_000diagc_000diagc_,where 0 is a null matrix of dimension N×N and diagc_ transforms the numerical array c_ in a diagonal matrix whose diagonal elements are the values contained in c_.The discretized data and state equations in (4) are combined together to obtain the following nonlinear relationship:
(10)e_sv=Gobsvdiag3c_I−Ginvdiag3c_−1e_iv.Moreover, in order to exploit the information of the considered multiview setup, the numerical arrays for each illumination are stacked as follows:(11)e_s=e_s1…e_sV=Gobs1diag3c_I−Ginvdiag3c_−1e_i1…GobsVdiag3c_I−Ginvdiag3c_−1e_iV=Fc_,where F is the discrete multiview operator that describes the mapping between the contrast function values c_ and the measured samples of the scattered electric field e_s.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 lp spaces [43]. In particular, the unknown discrete numerical array c_ is assumed to belong to the discrete Banach space C=lpℂN, whereas the data array e_s belongs to the linear space E=lpℂ3MV. 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 c_u (the subscript u denotes the outer iteration index). This loop is initialized with a starting guess c_0, 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., c_0=0_. The linearization is performed by computing a first-order Taylor expansion of the discrete operator F that involves the calculation of the Jacobian matrix Fu′ at c_u, which, for the 3D microwave imaging problem at hand, is defined as [3]
(12)Fu′=Gobs,u1diag3e_t,u1⋮Gobs,uVdiag3e_t,uV,where Gobs,uv and e_t,uv 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 c_u. Such quantities are obtained as [3]
(13)Gobs,uv=GobsvI−diag3c_uGinv−1,e_t,uv=I−Ginvdiag3c_u−1e_iv,and are efficiently computed by using the BiCGSTAB-FFT method [49].Figure 2
Flowchart of the developedlp Newton-CG inversion procedure.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 discretelp 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 1<p<2, to reduce the oversmoothing and ringing effects that are usually associated to the conventional reconstructions in l2 Hilbert spaces [40, 44]. The key point for extending the standard CG inversion scheme to Banach spaces is represented by the duality maps JC, JE, and JC∗ that are used in the update formulas in Figure 2 [51]. In particular, for the considered lp Banach space C, by virtue of the Asplund theorem, the duality map JC is defined as
(14)JCh_=h_p2−ph1p−1signh1⋮hNp−1signhN.Analogous definitions hold for the duality mapJE of data space E and for the one of the dual space of C, JC∗ (where the Hölder conjugate of p, i.e., q=p/p−1, is used instead of p) [51]. Moreover, the parameters αk and βk are computed by using the following relationships:
(15)αk=argminα>0Fu′JC∗JCh_k+αp_k−r_up2,(16)βk=Fu′∗JEFu′h_k−r_u22Fu′∗JEFu′h_k−1−r_u22,where Fu′∗ denotes the Hermitian transpose of the Jacobian matrix Fu′ 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 lp spaces [43]. However, with respect to the Landweber algorithm, the extension to lp settings of the classical CG algorithm poses additional difficulties, since there is not any explicit formula for computing the optimal step size αk (which exists in the standard l2 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ϵ0, magnetic permeability μ0, and wavenumber k0=ωμ0ϵ0). The working angular frequency is ω=2πf, with f=3×108Hz. The target is sequentially illuminated by V=6θ-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 M=82 measurement points uniformly distributed on a sphere of radius RDobs=2λ0, λ0=1m 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 Vinv is a cube of side LVinv=λ0, which has been subdivided into N=8000 cubic voxels inside the inversion procedure. The simulated data have been computed by using a custom method-of-moment numerical code, in which Vinv is partitioned into Nfwd=29791 cubic subdomains. A signal-to-noise ratio SNR=25dB on the total electric field data has been obtained by adding to the simulated values a white Gaussian noise with a zero mean value.Figure 3
Measurement and illumination setup for the numerical validation of the proposed 3D full-vector microwave imaging procedure.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 iterationsUmax=20 and Kmax=10 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 (einv), the average relative errors on the contrast function of the target (et) and of the background (eb), defined as
(17)einv=∑n=1Nc~rninv−crninv2∑n=1Ncrninv2,et/b=1Nt/b∑n∈It/bc~rninv−crninvcrninv+1,where c~rninv is the reconstructed value of the contrast function in the nth subdomain in which the investigation volume Vinv has been partitioned,crninv is the corresponding actual value, It/b are the sets of indexes of the cells inside the target and background regions, respectively, and Nt/b are the corresponding numbers of contained cells. It is interesting to note that the normalized root mean squared error einv 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 p in several scenarios, although it does not have an immediate physical interpretation. Meanwhile, the average relative errors et and eb 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 atrc,1=0.1,0.1,−0.1λ0, characterized by a relative dielectric permittivity ϵr=2.5, and with radius a1∈0.15,0.3λ0. The reconstructed distributions of the relative dielectric permittivity obtained for a1=0.3λ0 with lp space exponents p=popt=1.6 and p=2 (Hilbert space) are reported in Figure 4. The optimal value popt of the parameter p has been selected on the basis of the normalized mean squared error, which is shown in Figure 5 versus p for different values of the sphere radius a1. Considering both reconstructions and errors on the solution, it is clear that the proposed lp space method is able to outperform the classic Hilbert space approach for p∈1,2. Furthermore, the optimal value of p monotonically increases with the target size. For comparison purposes, an inexact Newton approach with Landweber iterations in lp 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 a1=0.3λ0 and for the optimal exponent parameters (popt=1.6 and popt=1.3 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 αk (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).Figure 4
Three-dimensional view of the reconstructed distribution of the dielectric permittivity with (a)p=popt=1.6 and (b) p=2. Single dielectric sphere with a1=0.3λ0.
(a)
(b)Figure 5
Behaviour ofeinv versus the norm parameter p for different values of the radius a1. Single dielectric sphere.Figure 6
Behaviour of the normalized data residuals for the conjugate gradient-based (top,popt=1.6) and Landweber-based (bottom, popt=1.3) methods versus the iteration number. Residuals in the outer inexact Newton (IN) iterations and in the inner conjugate gradient (CG) and Landweber (LW) loops are reported. Single dielectric sphere with a1=0.3λ0.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 atrc,1=0.25,0.25,0λ0, with radius a1=0.125λ0, and relative dielectric permittivity ϵr,1∈2,5. The second target is a z-directed dielectric cylinder with circular cross section, characterized by radius a2=0.125λ0, height h2=0.4λ0, center rc,2=−rc,1, and relative dielectric permittivity ϵr,2=2. The reconstructed distributions of the relative dielectric permittivity for ϵr,1=3 with popt=1.2 and p=2 (Hilbert space) are reported in Figures 7 and 8, respectively. The targets are visible in both cases, but the reconstruction obtained in the lpopt space is more accurate and allows us to appreciate the difference in the dielectric permittivity of the two objects. The behaviour of einv, reported versus the parameter p in Figure 9 for different values of ϵr,1, clearly indicates that the overall best results are achieved with p∈1.2,1.3. The average relative errors on the target and background regions and the numbers of performed Newton iterations Uopt, are reported in Table 1 for p=popt and p=2. 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.Figure 7
Reconstructed distribution of the dielectric permittivity withp=popt=1.2. (a) Three-dimensional view; (b) xz-plane cut (y=−0.25λ0); and (c) yz-plane cut (x=0.25λ0). Dielectric sphere with ϵr,1=3 and cylinder with ϵr,2=2.
(a)
(b)
(c)Figure 8
Reconstructed distribution of the dielectric permittivity withp=2. (a) Three-dimensional view; (b) xz-plane cut (y=−0.25λ0); and (c) yz-plane cut (x=0.25λ0). Dielectric sphere with ϵr,1=3 and cylinder with ϵr,2=2.
(a)
(b)
(c)Figure 9
Behaviour ofeinv versus the parameter p for different values of the relative dielectric permittivity ϵr,1. Dielectric sphere and cylinder.Table 1
Average relative errors (et/b) and number of outer iterations Uopt versus ϵr,1 with p=popt and p=2.
ϵr,1
lpopt
l2 (Hilbert)
popt
et
eb
Uopt
et
eb
Uopt
2
1.2
0.20
0.04
3
0.35
0.07
3
3
1.2
0.28
0.04
3
0.39
0.10
2
4
1.3
0.31
0.05
3
0.40
0.10
2
5
1.2
0.31
0.06
5
0.42
0.11
3However, 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 Figure10. 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.Figure 10
Behaviour of the normalized data residuals for the conjugate gradient-based (top,popt=1.3) and Landweber-based (bottom, popt=1.2) methods versus the iteration number. Residuals in the outer inexact Newton (IN) iterations and in the inner conjugate gradient (CG) and Landweber (LW) loops are reported. Dielectric sphere and cylinder with ϵr,1=4.


### 3.3. Inhomogeneous Targets with Different SNRs

To further validate the approach, an inhomogeneous structure has been considered. It consists of az-directed circular dielectric cylinder with radius a1=0.2λ0, height h1=0.6λ0, and centered at rc=0.15,0.15,0λ0. The lowest half of the cylinder is characterized by a relative dielectric permittivity ϵr,l=2, whereas the upper one has ϵr,u=3. In this case, we also consider a variation in the signal-to-noise ratio SNR∈5,50dB. The normalized mean squared error einv and the values of popt for each SNR are listed in Table 2, where the average relative errors and number of performed outer iterations Uopt are also reported. It can be seen that einv increases, as expected, while Uopt decreases (i.e., less iterations are executed) with a decrease of SNR. The reconstructed distribution of the relative dielectric permittivity for SNR=20dB (popt=1.5) 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 lp space procedure are evident in all the presented cases. Finally, the behaviours of the normalized data residuals are shown in Figure 12 for SNR=30dB and considering the optimal value of the parameter p (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.Table 2
Optimal values ofp, normalized mean squared error (einv), average relative errors (et/b), and number of outer iterations Uopt versus SNR with p=popt and p=2.
SNR (dB)
lpopt
l2 (Hilbert)
popt
einv
et
eb
Uopt
einv
et
eb
Uopt
50
1.8
0.54
0.22
0.09
6
0.54
0.22
0.11
6
40
1.7
0.54
0.23
0.09
7
0.56
0.22
0.13
8
30
1.5
0.53
0.23
0.08
4
0.59
0.22
0.14
4
20
1.5
0.62
0.26
0.13
4
0.71
0.26
0.18
4
10
1.6
0.88
0.41
0.18
3
0.90
0.42
0.20
3
5
1.1
1.02
0.58
0.02
1
1.06
0.37
0.32
3Figure 11
Reconstructed distribution of the dielectric permittivity for the signal-to-noise ratioSNR=20dB. (a) popt=1.5 and (b) p=2. Inhomogeneous cylinder.
(a)
(b)Figure 12
Behaviour of the normalized data residuals for the conjugate gradient-based (top,popt=1.5) and Landweber-based (bottom, popt=1.3) methods versus the iteration number. Residuals in the outer inexact Newton (IN) iterations and in the inner conjugate gradient (CG) and Landweber (LW) loops are reported. Inhomogeneous cylinder with SNR=30dB.


## 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 atrc,1=0.1,0.1,−0.1λ0, characterized by a relative dielectric permittivity ϵr=2.5, and with radius a1∈0.15,0.3λ0. The reconstructed distributions of the relative dielectric permittivity obtained for a1=0.3λ0 with lp space exponents p=popt=1.6 and p=2 (Hilbert space) are reported in Figure 4. The optimal value popt of the parameter p has been selected on the basis of the normalized mean squared error, which is shown in Figure 5 versus p for different values of the sphere radius a1. Considering both reconstructions and errors on the solution, it is clear that the proposed lp space method is able to outperform the classic Hilbert space approach for p∈1,2. Furthermore, the optimal value of p monotonically increases with the target size. For comparison purposes, an inexact Newton approach with Landweber iterations in lp 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 a1=0.3λ0 and for the optimal exponent parameters (popt=1.6 and popt=1.3 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 αk (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).Figure 4
Three-dimensional view of the reconstructed distribution of the dielectric permittivity with (a)p=popt=1.6 and (b) p=2. Single dielectric sphere with a1=0.3λ0.
(a)
(b)Figure 5
Behaviour ofeinv versus the norm parameter p for different values of the radius a1. Single dielectric sphere.Figure 6
Behaviour of the normalized data residuals for the conjugate gradient-based (top,popt=1.6) and Landweber-based (bottom, popt=1.3) methods versus the iteration number. Residuals in the outer inexact Newton (IN) iterations and in the inner conjugate gradient (CG) and Landweber (LW) loops are reported. Single dielectric sphere with a1=0.3λ0.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 atrc,1=0.25,0.25,0λ0, with radius a1=0.125λ0, and relative dielectric permittivity ϵr,1∈2,5. The second target is a z-directed dielectric cylinder with circular cross section, characterized by radius a2=0.125λ0, height h2=0.4λ0, center rc,2=−rc,1, and relative dielectric permittivity ϵr,2=2. The reconstructed distributions of the relative dielectric permittivity for ϵr,1=3 with popt=1.2 and p=2 (Hilbert space) are reported in Figures 7 and 8, respectively. The targets are visible in both cases, but the reconstruction obtained in the lpopt space is more accurate and allows us to appreciate the difference in the dielectric permittivity of the two objects. The behaviour of einv, reported versus the parameter p in Figure 9 for different values of ϵr,1, clearly indicates that the overall best results are achieved with p∈1.2,1.3. The average relative errors on the target and background regions and the numbers of performed Newton iterations Uopt, are reported in Table 1 for p=popt and p=2. 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.Figure 7
Reconstructed distribution of the dielectric permittivity withp=popt=1.2. (a) Three-dimensional view; (b) xz-plane cut (y=−0.25λ0); and (c) yz-plane cut (x=0.25λ0). Dielectric sphere with ϵr,1=3 and cylinder with ϵr,2=2.
(a)
(b)
(c)Figure 8
Reconstructed distribution of the dielectric permittivity withp=2. (a) Three-dimensional view; (b) xz-plane cut (y=−0.25λ0); and (c) yz-plane cut (x=0.25λ0). Dielectric sphere with ϵr,1=3 and cylinder with ϵr,2=2.
(a)
(b)
(c)Figure 9
Behaviour ofeinv versus the parameter p for different values of the relative dielectric permittivity ϵr,1. Dielectric sphere and cylinder.Table 1
Average relative errors (et/b) and number of outer iterations Uopt versus ϵr,1 with p=popt and p=2.
ϵr,1
lpopt
l2 (Hilbert)
popt
et
eb
Uopt
et
eb
Uopt
2
1.2
0.20
0.04
3
0.35
0.07
3
3
1.2
0.28
0.04
3
0.39
0.10
2
4
1.3
0.31
0.05
3
0.40
0.10
2
5
1.2
0.31
0.06
5
0.42
0.11
3However, 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 Figure10. 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.Figure 10
Behaviour of the normalized data residuals for the conjugate gradient-based (top,popt=1.3) and Landweber-based (bottom, popt=1.2) methods versus the iteration number. Residuals in the outer inexact Newton (IN) iterations and in the inner conjugate gradient (CG) and Landweber (LW) loops are reported. Dielectric sphere and cylinder with ϵr,1=4.


## 3.3. Inhomogeneous Targets with Different SNRs

To further validate the approach, an inhomogeneous structure has been considered. It consists of az-directed circular dielectric cylinder with radius a1=0.2λ0, height h1=0.6λ0, and centered at rc=0.15,0.15,0λ0. The lowest half of the cylinder is characterized by a relative dielectric permittivity ϵr,l=2, whereas the upper one has ϵr,u=3. In this case, we also consider a variation in the signal-to-noise ratio SNR∈5,50dB. The normalized mean squared error einv and the values of popt for each SNR are listed in Table 2, where the average relative errors and number of performed outer iterations Uopt are also reported. It can be seen that einv increases, as expected, while Uopt decreases (i.e., less iterations are executed) with a decrease of SNR. The reconstructed distribution of the relative dielectric permittivity for SNR=20dB (popt=1.5) 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 lp space procedure are evident in all the presented cases. Finally, the behaviours of the normalized data residuals are shown in Figure 12 for SNR=30dB and considering the optimal value of the parameter p (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.Table 2
Optimal values ofp, normalized mean squared error (einv), average relative errors (et/b), and number of outer iterations Uopt versus SNR with p=popt and p=2.
SNR (dB)
lpopt
l2 (Hilbert)
popt
einv
et
eb
Uopt
einv
et
eb
Uopt
50
1.8
0.54
0.22
0.09
6
0.54
0.22
0.11
6
40
1.7
0.54
0.23
0.09
7
0.56
0.22
0.13
8
30
1.5
0.53
0.23
0.08
4
0.59
0.22
0.14
4
20
1.5
0.62
0.26
0.13
4
0.71
0.26
0.18
4
10
1.6
0.88
0.41
0.18
3
0.90
0.42
0.20
3
5
1.1
1.02
0.58
0.02
1
1.06
0.37
0.32
3Figure 11
Reconstructed distribution of the dielectric permittivity for the signal-to-noise ratioSNR=20dB. (a) popt=1.5 and (b) p=2. Inhomogeneous cylinder.
(a)
(b)Figure 12
Behaviour of the normalized data residuals for the conjugate gradient-based (top,popt=1.5) and Landweber-based (bottom, popt=1.3) methods versus the iteration number. Residuals in the outer inexact Newton (IN) iterations and in the inner conjugate gradient (CG) and Landweber (LW) loops are reported. Inhomogeneous cylinder with SNR=30dB.


## 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 f=3×109Hz, and the measurement configuration is composed of V=25 sources on a sphere with radius RDobs=1.796m, and M=27 measurement points on a circumference with the same radius. More details about the experimental setup can be found in [54]. The investigation volume Vinv is a cubic region with side LVinv=0.15m centered at the origin, partitioned into N=8000 voxels for the inverse problem solution. The trends of the reconstruction errors versus the parameter p for both the considered targets are shown in Figure 13, whereas the reconstructions for p=popt and p=2 are reported in Figures 14 and 15. These results further confirm the capabilities of the CG-based inexact Newton approach in lp spaces, with 1<p<2, to accurately retrieve a quantitative reconstruction of the scene under test, outperforming the standard method in Hilbert spaces.Figure 13
Behaviour of the reconstruction errors versus the norm parameterp. (a) TwoSpheres and (b) TwoCubes targets.
(a)
(b)Figure 14
Reconstructed distributions of the relative dielectric permittivity with (a)p=popt=1.4 and (b) p=2. TwoSpheres target.
(a)
(b)Figure 15
Reconstructed distributions of the dielectric permittivity with (a)p=popt=1.2 and (b) p=2. TwoCubes target.
(a)
(b)In particular, as shown in Figure13, the background error eb increases as p grows. Such a behaviour can be related to the fact that low values of p 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 p, 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 p 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 popt 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.Table 3
Optimal values ofp, normalized mean squared error (einv), average relative errors (et/b), and number of outer iterations Uopt for the two considered experimental datasets. Comparison with the results obtained with the Landweber-based procedure in [44].
CG
LW
popt
einv
et
eb
Uopt
popt
einv
et
eb
Uopt
TwoSpheres
1.4
0.66
0.31
0.09
6
1.3
0.64
0.33
0.07
12
TwoCubes
1.2
0.74
0.32
0.02
5
1.2
0.70
0.32
0.01
11


## 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 thelp 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 optimallp exponent value popt. 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).



---

*Source: 2841937-2019-05-14.xml*