Microwave Imaging of Three-Dimensional Targets by Means of an Inexact-Newton-Based Inversion Algorithm

A microwave imaging method previously developed for tomographic inspection of dielectric targets is extended to three-dimensional objects. The approach is based on the full vector equations of the electromagnetic inverse scattering problem. The ill-posedness of the problem is faced by the application of an inexact-Newton method. Preliminary reconstruction results are reported.


Introduction
Microwave imaging is a technique aimed at inspecting unknown targets by using interrogating electromagnetic waves at microwave frequencies.As it is well known, the image formation is based on the solution of an electromagnetic inverse scattering problem.One of the main features of microwave imaging is the ability of providing the distributions of the dielectric parameters (e.g., dielectric permittivity and electric conductivity), which cannot be directly retrieved by using other diagnostic techniques [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18].From a theoretical point of view, this achievement can be obtained even in the case of three-dimensional (3D) configurations.As a matter of fact, some initial works in this area were just oriented to 3D imaging [19][20][21][22][23][24][25][26].Unfortunately, the computational resources needed to face a full vector 3D inverse scattering problem has been prohibitive for a long time.Consequently, the research community mainly focused his efforts on two-dimensional (2D) settings.In particular, the presence of cylindrical targets (uniform and homogeneous along an infinite direction) has been usually assumed, leading to the well known microwave tomography .
In this context, the present paper introduces for the first time an extension to 3D scatterers of a diagnostic tomographic approach recently developed by the authors [70][71][72][73][74].It is based on a two-step inexact-Newton method for the solution of two electric field integral equations (EFIEs), namely the data equation and the state equation, which represent the physical model of the whole scattering process.By a theoretical point of view, the approach is quite a classic one.Anyway, the particular inexact-Newton iterative solving procedure has been found to be quite robust against the illposedness of the considered inverse scattering problem both by assuming synthetic or experimental data.
In the 3D case under investigation, the full vector equations, written in terms of the Green's dyadic function for free space, are discretized and combined in order to obtain a single functional equation, which is iteratively solved by Gauss-Newton linearization (outer step) and successive regularization by truncated Landweber method (inner step).
The paper is organized as follows.The mathematical formulation of the developed approach is discussed in Section 2. Section 3 reports some initial numerical results aimed at validating the inversion procedure.Finally, conclusions are drawn in Section 4.

Mathematical Formulation
Let us consider the configuration of Figure 1.An unknown dielectric target is located inside a predefined volume (in the following denoted as investigation domain,  inv ).The object is illuminated by one or more incident electric fields E inc (r) (generated by proper antennas located outside the investigation area).The total electric field E tot (r) produced by the interaction between the target and the illuminating radiation is collected by means of RX antennas located in a measurement domain  obs .As we will better analyze in Section 3, we mention that in our numerical simulation the investigation domain is a cube, while the TX and RX antennas are uniformly distributed on an outside sphere of predefined radius.
For sake of simplicity, a single-view case is described in the following.The extension of the mathematical formulation to the multiview case is however straightforward.Moreover, a   time dependence is assumed and the corresponding terms are omitted in the following.

Electromagnetic Model.
The scattered electric field E scatt (r) = E tot (r) − E inc (r) is related to the dielectric properties of the investigated area  inv by means of the following two integral equations [1], usually referred to as data and state equations, where (r) =   (r) − 1 is the contrast function (with   being the space dependent relative complex dielectric permittivity of the investigation domain) and the linear operators G data and G state are defined as with  0 =  √  0  0 being the free-space wavenumber and the free-space dyadic Green's function.
The data and state equations can be merged in order to obtain a single nonlinear operator equation relating the contrast function and the scattered electric field as follows where r ∈  obs and C is the operator such that C(f)(r) = (r)f(r), r ∈  inv .By defining the nonlinear operator F() = G data (C((I − G state C) −1 E inc )) that maps the contrast function  with the scattered electric field E scatt , (4) can be written in compact form as follows: The nonlinear equation ( 5) models the full inverse scattering problem: given the scattered electric field E scatt (i.e., the data), measured in the measurement domain  obs , find the contrast function  (i.e., the unknown) in the investigation domain  inv such that F() = E scatt .

Regularized Inversion
Algorithm.The integral equation ( 5) for the computation of the contrast function  turns out to be ill posed, and its solution requires a regularization algorithm.To this end, the inner/outer regularizing scheme developed in [70][71][72] for the 2D case is here generalized to a 3D configuration.The inner/outer scheme can be summarized as follows.Any outer step is a basic Gauss-Newton linearization of the nonlinear equation ( 5).Such a linearized equation is then solved by means of the truncated Landweber iterative method, which represents the inner step.
In particular, the algorithm works as follows.
(1) Initialize the outer loop by setting  = 0 ( denotes the outer iteration index) and by choosing a starting guess  0 (if no a priori information is available, just void domain can be used, i.e.,  0 = 0).
(2) Linearize ( 5) in order to obtain the following linear equation ("outer loop") where (3) Compute a regularized solution of the linear equation ( 6) with respect to the unknown ℎ  (r), r ∈  inv , by using the following truncated Landweber algorithm ("inner loop", where  denotes the corresponding inner iteration index) where denotes the adjoint operator of F    .The inner iterations are stopped when a maximum number of iterations  is reached.
International Journal of Antennas and Propagation 3 (4) Update the current solution as (5) Set  =  + 1 and iterate from step (2) until a predefined stopping rule (such as discrepancy principle or GCV) is satisfied or a maximum number of outer iterations  is reached.
We briefly recall that the number of inner step  acts as the regularization parameter of the truncated Landweber algorithm for the solution of any Newton equation (6).Basically, during the first iterations the Landweber method is able to filter out the components usually corrupted by noise (i.e., the highest frequencies Fourier components of the data E  ), so that an early stop of the inner iterations gives rise to a regularization effect [71].
Similarly to the two-dimensional case, the Frechét derivative of the operator F is given by where

and the operator G 𝑐
data is given by (10) with G  being the dyadic Green's function for an inhomogeneous background characterized by a contrast function .

Discretization of the Scattering Equations.
The data and state scattering equations (1) are discretized by using pulse basis functions and Dirac's delta weighting functions [75].In particular, the investigation domain  inv is discretized into  cubic voxels of side  and centers r inv  ,  = 1, . . ., , and the observation domain is composed by  measurement points located at positions r meas  ,  = 1, . . ., .Consequently, the discrete version of ( 4) is  (the superscript denotes the components of the field vector), and  3 (⋅) is the discrete diagonal operator defined as In (11) the matrices G data and G state are the discrete counterparts of the linear operators defined in (2) and they contain the integrals of the components of the corresponding dyadic Green's function over any voxels' volume [76].The discrete version of the Frechét derivative used in ( 6) is obtained in a similar way.
It is worth noting that, in the numerical implementation of the approach, a BiCGStab-FFT approach [77,78] has been used to speed up the computation of the total internal field and of the inhomogeneous dyadic Green's function.

Numerical Results
The developed approach has been tested by means of several numerical simulations.In all cases, a plane wave illumination is assumed.Moreover, a multiview configuration, in which the object is illuminated by subsequent different waves impinging from   uniformly distributed directions, is used in order to increase the available information.The scattered field is collected into   measurement points uniformly distributed on a sphere of radius   for any impinging wave.The synthetic data used for the inversion have been obtained by using a numerical simulator based on the Method of Moments [75].A finer discretization than that used in the inversion code is employed in order to avoid inverse crimes.Moreover, the computed data have been corrupted with a Gaussian noise with zero-mean value and variance corresponding to a signal-to-noise ratio SNR.
As a first case, a single homogeneous cubic target is considered.The center of the object is located in r  = (0.1 0 , 0.1 0 , 0.1 0 ), whereas its radius and relative dielectric permittivity are equal to  = 0.5 0 and   = 2.0, respectively.The investigation area is a cubic volume of side  =  0 , which has been partitioned into  = 20 × 20 × 20 = 8000 subdomains.The number of views is   = 6, whereas the number of measurement points is   = 82 and the radius of the measurement sphere is   =  0 .The maximum number of allowed outer iterations has been set equal to  = 20.The performance of the approach have been firstly analyzed with respect to the signal-to-noise ratio and to the number of iterations of the inner loop.The following two error measures are used: The first error measure,  data , can be used in real setting for implementing the discrepancy principle or the GCV, while the second one,  RMSE , can be used only for numerical testing, since it requires the knowledge of the target c actual .
Table 1 reports the errors obtained for the considered cases and the number of outer iterations needed to reach the best reconstruction.As expected, with noiseless data the inversion algorithm is always able to converge to a correct solution and a higher value of  allows to obtain slightly better results.On the contrary, when noise is present, low values of  provide a stronger regularization, thus allowing for a better reconstruction.It is worth noting that in these cases, a higher number of iterations are needed to reach the best solution, leading to a higher computational time.Moreover, Figure 2 shows the behavior of the error measures versus the number of outer iterations and for different values of .The semiconvergence [79] effect is clearly visible in Figure 2(b), confirming that lower values of  provide a stronger regularization, leading however to a higher number of iterations needed to reach convergence.
An example of the reconstructed distribution of the relative dielectric permittivity, concerning the case in which SNR = 20 dB and  = 3, is reported in Figure 3.As can be seen, the target is correctly identified.As expected, the shape is smoothed due to the regularization effect of the inversion algorithm.
The effects of the number of views and of the measurement point number on the quality of the reconstruction have also been evaluated.The obtained results are summarized in Table 2 (errors versus number of views) and Table 3 (errors versus number of measurement points).As expected, as the number of views or the number of measurement points increases, the reconstruction error decreases.However, the computational time increases as well (since the dimensions of the matrices are higher and more iterations are needed).
Finally, the performance of the inversion approach with respect to the value of the relative dielectric permittivity of the target has been evaluated, too.The obtained results are reported in Table 4.As can be seen, in all the considered cases the developed algorithm provides comparable reconstruction results.However, as expected, for the higher values of dielectric permittivity, a higher number of iterations is required for obtaining a good solution.This is basically due to the effects of the Newton linearization, which is less accurate for strong scatterers (i.e., large values dielectric permittivity), with respect to the original nonlinear operator.
In the second case, a nonhomogeneous structure has been considered.The target is a cube with a void inclusion.The parameters of the cube are as follows:  = 0.9 0 , center r  = (0, 0, 0), and   = 1.5.The inclusion has a spherical shape with radius   = 0.25 0 and center r  = (0, 0, 0).The investigation area is a cubic volume of side  = 1.5 0 , which has been partitioned into  = 28 × 28 × 28 = 21952 subdomains.The maximum number of allowed outer iterations has been set equal to  = 20, whereas the number of inner iterations is  = 2.The number of views is   = 6, whereas   = 101 measurement points uniformly distributed on a sphere of radius   = 1.5 0 has been considered.The distribution of the reconstructed dielectric permittivity is shown in Figure 4.As can be seen, the reconstruction algorithm provides a good  reconstruction of the target, providing a good restoration of the void inclusion too.In the third case, a more complex configuration with two different separate objects inside the investigation area is assumed.The first one is a cube characterized by  = 0.7 0 , center r  = (0.4 0 , 0.4 0 , 0.4 0 ), and   = 1.5.The second one is a sphere with radius   = 0.4 0 , centered located at r  = (0.5 0 , −0.5 0 , −0.5 0 ), and relative dielectric permittivity   = 1.5.The investigation area is a cubic volume of side  = 2 0 , which has been partitioned into  = 30 × 30 × 30 = 27000 subdomains.  = 6 views have been used and, for every view, the field is collected in   = 101 points on a sphere of radius   = √ 3 0 .Similarly to the previous case, the maximum number of allowed outer iterations has been set equal to  = 20 and the number of inner iterations is  = 2.The reconstructed distribution of the relative dielectric permittivity is reported in Figure 5.As can be seen from this figure, in this case, too, the two targets are correctly reconstructed.It is worth noting that the two cuts in Figures 5(b) and 5(c) refer to two planes passing only through the first object.Consequently, only the cubic target is visible in these two figures.On the contrary, the  plane considered in Figure 5(d) contains the both targets.

Conclusions
In this paper, an algorithm for three-dimensional nondestructive diagnostics through microwave inverse scattering has been analyzed.The dielectric properties of an object (i.e., its contrast function) has to be restored starting from the scattered electric field generated by the interaction with a known electric incident field.The developed algorithm is based on the computation of a regularized solution of the nonlinear equation that relates the scattered electric field in an external observation domain (i.e., the data) with the contrast function in an inaccessible investigation domain (i.e., the unknown).In particular, the approach is based on an outer-inner iterative scheme, where the outer iteration is International Journal of Antennas and Propagation the basic Newton linearization and the inner iteration is the truncated Landweber method.Several preliminary synthetic numerical tests have been performed in order to validate the developed method.Both homogenous and nonhomogeneous targets have been considered, and in all cases the inversion scheme provided good restorations also for noisy data.A natural subsequent step that will be pursued in the future is the evaluation of the performance of the method against real data.
where c = [(r inv 1 ) ⋅ ⋅ ⋅ (r inv  )]  is an array containing the values of the contrast function inside the voxels, e tot/scatt/inc = [ [

Figure 2 :Figure 3 :
Figure 2: Error measures versus the number of outer iterations for different values of the number of inner iterations.Homogeneous dielectric cube.SNR = 20 dB.

Table 1 :
Error measures for different values of the signal-to-noise ratio and of the number of inner iterations.Homogeneous dielectric cube.

Table 2 :
Error measures for different values of the number of views.Homogeneous dielectric cube.

Table 3 :
Error measures for different values of the number of measurement points.Homogeneous dielectric cube.

Table 4 :
Error measures for different values of the relative dielectric permittivity.Homogeneous dielectric cube.