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

Authors: Claudio Estatico; Matteo Pastorino; Andrea Randazzo
Journal: International Journal of Antennas and Propagation
 (2013)

---

## Abstract

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.


---

## Body

## 1. 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–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–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 [27–49].More recently, the enormous improvements in computer powers have allowed reconsidering the 3D formulation. In addition, new efficient solving methods have been developed in the past years. Most of them allow a more or less straightforward extension for the inspection of 3D targets. Some examples of these approaches are reported in [50–69] and in the references therein.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–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 ill-posedness 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 Section2. Section 3 reports some initial numerical results aimed at validating the inversion procedure.  Finally, conclusions are drawn in Section 4.


## 2. Mathematical Formulation

Let us consider the configuration of Figure1. An unknown dielectric target is located inside a predefined volume (in the following denoted as investigation domain, Vinv). The object is illuminated by one or more incident electric fields Einc(r) (generated by proper antennas located outside the investigation area). The total electric field Etot(r) produced by the interaction between the target and the illuminating radiation is collected by means of RX antennas located in a measurement domain Vobs. 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.Figure 1
Imaging configuration.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, aejωt time dependence is assumed and the corresponding terms are omitted in the following.


### 2.1. Electromagnetic Model

The scattered electric fieldEscatt(r)=Etot(r)-Einc(r) is related to the dielectric properties of the investigated area Vinv by means of the following two integral equations [1], usually referred to as data and state equations,
(1)Escatt(r)=Gdata(cEtot)(r),r∈Vobs,Einc(r)=Etot(r)-Gstate(cEtot)(r),r∈Vinv,
where c(r)=ϵr(r)-1 is the contrast function (with ϵr being the space dependent relative complex dielectric permittivity of the investigation domain) and the linear operators Gdata and Gstate are defined as
(2)Gdataf(r)=-k02∫Vinvf(r)·G-0(r,r′)dr′,r∈Vobs,Gstatef(r)=-k02∫Vinvf(r)·G-0(r,r′)dr′,r∈Vinv
with k0=ωϵ0μ0 being the free-space wavenumber and
(3)G-0(r,r′)=-14π(I-+∇∇k02)e-jk|r-r′||r-r′|
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(4)Escatt(r)=Gdata(C((I-GstateC)-1Einc))(r),
where r∈Vobs and C is the operator such that C(f)(r)=c(r)f(r), r∈Vinv. By defining the nonlinear operator F(c)=Gdata(C((I-GstateC)-1Einc)) that maps the contrast function c with the scattered electric field Escatt, (4) can be written in compact form as follows:
(5)F(c)(r)=Escatt(r),r∈Vobs.The nonlinear equation (5) models the full inverse scattering problem: given the scattered electric field Escatt (i.e., the data), measured in the measurement domain Vobs, find the contrast function c (i.e., the unknown) in the investigation domain Vinv such that  F(c)=Escatt.


### 2.2. Regularized Inversion Algorithm

The integral equation (5) for the computation of the contrast function c turns out to be ill posed, and its solution requires a regularization algorithm. To this end, the inner/outer regularizing scheme developed in [70–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 settingk=0 (k denotes the outer iteration index) and by choosing a starting guess c0 (if no a priori information is available, just void domain can be used, i.e., c0=0).(2)
Linearize (5) in order to obtain the following linear equation (“outer loop”)
(6)Fck′hk(r)=Ek(r),r∈Vobs,whereFck′ is the Frechét derivative of F at ck and Ek(r)=Escatt(r)-F(ck)(r),r∈Vobs.(3)
Compute a regularized solution of the linear equation (6) with respect to the unknown hk(r), r∈Vinv, by using the following truncated Landweber algorithm (“inner loop”, where l denotes the corresponding inner iteration index)
(7)hk,0(r)=0,hk,l+1(r)=hk,l(r)-βFck′*(Fck′hk,l(r)-Ek(r)),whereβ=∥Fck′*Fck′∥-1 is a constant values. Here Fck′* denotes the adjoint operator of Fck′. The inner iterations are stopped when a maximum number of iterations L is reached.(4)
Update the current solution as(8)ck+1(r)=ck(r)+hk,L(r).(5)
Setk=k+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 K is reached.We briefly recall that the number of inner stepL 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 Ek), 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 operatorF is given by
(9)Fc′h(r)=Gdatac(hEtotc)(r),r∈Vobs,
where Etotc(r)=((I-GstateC)-1Einc)(r), r∈Vinv, and the operator Gdatac is given by
(10)Gdatacf(r)=-k02∫Vinvf(r)·G-c(r,r′)dr′,r∈Vobs
with G-c being the dyadic Green’s function for an inhomogeneous background characterized by a contrast function c.


### 2.3. 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 Vinv is discretized into N cubic voxels of side l and centers rninv, n=1,…,N, and the observation domain is composed by M measurement points located at positions rmmeas,m=1,…,M.Consequently, the discrete version of (4) is
(11)escatt=GdataD3(c)(I-GstateD3(c))-1einc,
where c=[c(r1inv)⋯c(rNinv)]t is an array containing the values of the contrast function inside the voxels,
(12)etot/scatt/inc=[etot/scatt/incxetot/scatt/incyetot/scatt/incz]
with etot/incx/y/z=[Etot/incx/y/z(r1inv)⋯Etot/incx/y/z(rNinv)]t and escattx/y/z=[Escattx/y/z(r1meas)⋯Escattx/y/z(rMmeas)]t (the superscript denotes the components of the field vector), and D3(·) is the discrete diagonal operator defined as
(13)D3(c)=[c1⋱cN0c1⋱cN0c1⋱cN].In (11) the matrices Gdata and Gstate 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.


## 2.1. Electromagnetic Model

The scattered electric fieldEscatt(r)=Etot(r)-Einc(r) is related to the dielectric properties of the investigated area Vinv by means of the following two integral equations [1], usually referred to as data and state equations,
(1)Escatt(r)=Gdata(cEtot)(r),r∈Vobs,Einc(r)=Etot(r)-Gstate(cEtot)(r),r∈Vinv,
where c(r)=ϵr(r)-1 is the contrast function (with ϵr being the space dependent relative complex dielectric permittivity of the investigation domain) and the linear operators Gdata and Gstate are defined as
(2)Gdataf(r)=-k02∫Vinvf(r)·G-0(r,r′)dr′,r∈Vobs,Gstatef(r)=-k02∫Vinvf(r)·G-0(r,r′)dr′,r∈Vinv
with k0=ωϵ0μ0 being the free-space wavenumber and
(3)G-0(r,r′)=-14π(I-+∇∇k02)e-jk|r-r′||r-r′|
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(4)Escatt(r)=Gdata(C((I-GstateC)-1Einc))(r),
where r∈Vobs and C is the operator such that C(f)(r)=c(r)f(r), r∈Vinv. By defining the nonlinear operator F(c)=Gdata(C((I-GstateC)-1Einc)) that maps the contrast function c with the scattered electric field Escatt, (4) can be written in compact form as follows:
(5)F(c)(r)=Escatt(r),r∈Vobs.The nonlinear equation (5) models the full inverse scattering problem: given the scattered electric field Escatt (i.e., the data), measured in the measurement domain Vobs, find the contrast function c (i.e., the unknown) in the investigation domain Vinv such that  F(c)=Escatt.


## 2.2. Regularized Inversion Algorithm

The integral equation (5) for the computation of the contrast function c turns out to be ill posed, and its solution requires a regularization algorithm. To this end, the inner/outer regularizing scheme developed in [70–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 settingk=0 (k denotes the outer iteration index) and by choosing a starting guess c0 (if no a priori information is available, just void domain can be used, i.e., c0=0).(2)
Linearize (5) in order to obtain the following linear equation (“outer loop”)
(6)Fck′hk(r)=Ek(r),r∈Vobs,whereFck′ is the Frechét derivative of F at ck and Ek(r)=Escatt(r)-F(ck)(r),r∈Vobs.(3)
Compute a regularized solution of the linear equation (6) with respect to the unknown hk(r), r∈Vinv, by using the following truncated Landweber algorithm (“inner loop”, where l denotes the corresponding inner iteration index)
(7)hk,0(r)=0,hk,l+1(r)=hk,l(r)-βFck′*(Fck′hk,l(r)-Ek(r)),whereβ=∥Fck′*Fck′∥-1 is a constant values. Here Fck′* denotes the adjoint operator of Fck′. The inner iterations are stopped when a maximum number of iterations L is reached.(4)
Update the current solution as(8)ck+1(r)=ck(r)+hk,L(r).(5)
Setk=k+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 K is reached.We briefly recall that the number of inner stepL 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 Ek), 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 operatorF is given by
(9)Fc′h(r)=Gdatac(hEtotc)(r),r∈Vobs,
where Etotc(r)=((I-GstateC)-1Einc)(r), r∈Vinv, and the operator Gdatac is given by
(10)Gdatacf(r)=-k02∫Vinvf(r)·G-c(r,r′)dr′,r∈Vobs
with G-c being the dyadic Green’s function for an inhomogeneous background characterized by a contrast function c.


## 2.3. 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 Vinv is discretized into N cubic voxels of side l and centers rninv, n=1,…,N, and the observation domain is composed by M measurement points located at positions rmmeas,m=1,…,M.Consequently, the discrete version of (4) is
(11)escatt=GdataD3(c)(I-GstateD3(c))-1einc,
where c=[c(r1inv)⋯c(rNinv)]t is an array containing the values of the contrast function inside the voxels,
(12)etot/scatt/inc=[etot/scatt/incxetot/scatt/incyetot/scatt/incz]
with etot/incx/y/z=[Etot/incx/y/z(r1inv)⋯Etot/incx/y/z(rNinv)]t and escattx/y/z=[Escattx/y/z(r1meas)⋯Escattx/y/z(rMmeas)]t (the superscript denotes the components of the field vector), and D3(·) is the discrete diagonal operator defined as
(13)D3(c)=[c1⋱cN0c1⋱cN0c1⋱cN].In (11) the matrices Gdata and Gstate 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.


## 3. 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 fromNS uniformly distributed directions, is used in order to increase the available information. The scattered field is collected into NM measurement points uniformly distributed on a sphere of radius RM 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 inrc=(0.1λ0,0.1λ0,0.1λ0), whereas its radius and relative dielectric permittivity are equal to l=0.5λ0 and ϵr=2.0, respectively. The investigation area is a cubic volume of side D=λ0, which has been partitioned into N=20×20×20=8000 subdomains. The number of views is NS=6, whereas the number of measurement points is NM=82 and the radius of the measurement sphere is RM=λ0. The maximum number of allowed outer iterations has been set equal to K=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:
(14)edata=∑s=1S∥escatt-F(c)∥2∑s=1S∥escatt∥2,eRMSE=∥c-cactual∥∥cactual∥.The first error measure,edata, can be used in real setting for implementing the discrepancy principle or the GCV, while the second one, eRMSE, can be used only for numerical testing, since it requires the knowledge of the target cactual.Table1 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 L allows to obtain slightly better results. On the contrary, when noise is present, low values of L 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.Table 1
Error measures for different values of the signal-to-noise ratio and of the number of inner iterations. Homogeneous dielectric cube.
SNR
L
e
data
e
RMSE
K
opt
∞
1
0.088
0.57
20
2
0.06
0.54
20
3
0.046
0.52
20
4
0.038
0.52
20
5
0.032
0.51
20
10
0.027
0.51
20
30
1
0.48
0.64
20
2
0.49
0.65
12
3
0.48
0.65
9
4
0.48
0.66
8
5
0.48
0.67
8
10
0.48
0.72
6
20
1
0.62
0.66
16
2
0.62
0.67
11
3
0.61
0.68
8
4
0.62
0.68
6
5
0.62
0.69
6
10
0.61
0.74
6
10
1
0.87
0.73
7
2
0.89
0.75
5
3
0.88
0.76
4
4
0.89
0.88
2
5
0.88
0.81
3
10
0.91
0.91
2Moreover, Figure2 shows the behavior of the error measures versus the number of outer iterations and for different values of L. The semiconvergence [79] effect is clearly visible in Figure 2(b), confirming that lower values of L provide a stronger regularization, leading however to a higher number of iterations needed to reach convergence.Error measures versus the number of outer iterations for different values of the number of inner iterations. Homogeneous dielectric cube.SNR=20 dB.
(a)
(b)An example of the reconstructed distribution of the relative dielectric permittivity, concerning the case in whichSNR=20 dB and L=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.Reconstructed distribution of the relative dielectric permittivity. (a) 3D view; (b)xy plane (z=0.1λ0); (c) xz plane (y=0.1λ0); (d) yz plane (x=0.1λ0). Homogeneous dielectric cube. SNR=20 dB. L=3.
(a)
(b)
(c)
(d)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 Table2 (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).Table 2
Error measures for different values of the number of views. Homogeneous dielectric cube.
N
S
e
data
e
RMSE
K
opt
6
0.52
0.66
5
12
0.54
0.57
9
18
0.56
0.55
9
29
0.54
0.53
12Table 3
Error measures for different values of the number of measurement points. Homogeneous dielectric cube.
N
M
e
data
e
RMSE
K
opt
82
0.53
0.67
6
122
0.48
0.64
8
170
0.46
0.61
12
212
0.43
0.58
13Finally, 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 Table4. 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.Table 4
Error measures for different values of the relative dielectric permittivity. Homogeneous dielectric cube.
ε
r
e
data
e
RMSE
K
opt
1.5
0.78
0.7
3
1.75
0.63
0.66
6
2
0.53
0.64
10
2.25
0.44
0.61
14
2.5
0.38
0.59
16In 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:l=0.9λ0, center rc=(0,0,0), and ϵr=1.5. The inclusion has a spherical shape with radius rs=0.25λ0 and center rs=(0,0,0). The investigation area is a cubic volume of side D=1.5λ0, which has been partitioned into N=28×28×28=21952 subdomains. The maximum number of allowed outer iterations has been set equal to K=20, whereas the number of inner iterations is L=2. The number of views is NS=6, whereas NM=101 measurement points uniformly distributed on a sphere of radius RM=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.Reconstructed distribution of the relative dielectric permittivity. (a) 3D view; (b)xy plane (z=0); (c) xz plane (y=0); (d) yz plane (x=0). Dielectric cube with a void inclusion. SNR=25 dB. L=2.
(a)
(b)
(c)
(d)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 byl=0.7λ0, center rc=(0.4λ0,0.4λ0,0.4λ0), and ϵr=1.5. The second one is a sphere with radius rs=0.4λ0, centered located at rs=(0.5λ0,-0.5λ0,-0.5λ0), and relative dielectric permittivity ϵr=1.5. The investigation area is a cubic volume of side D=2λ0, which has been partitioned into N=30×30×30=27000 subdomains. NS=6 views have been used and, for every view, the field is collected in NM=101 points on a sphere of radius RM=3λ0. Similarly to the previous case, the maximum number of allowed outer iterations has been set equal to K=20 and the number of inner iterations is L=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 yz plane considered in Figure 5(d) contains the both targets.Reconstructed distribution of the relative dielectric permittivity. (a) 3D view; (b)xy plane (z=0.4λ0); (c) xz plane (y=0.4λ0); (d) yz plane (x=0.4λ0). Two separate objects. SNR=25 dB. L=2.
(a)
(b)
(c)
(d)


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



---

*Source: 407607-2013-08-18.xml*