# Nonlinear Electromagnetic Inverse Scattering Imaging Based on IN-LSQR

Authors: Huilin Zhou; Youwen Liu; Yuhao Wang; Liangbing Chen; Rongxing Duan
Journal: International Journal of Antennas and Propagation
 (2018)

---

## Abstract

A nonlinear inversion scheme is proposed for electromagnetic inverse scattering imaging. It exploits inexact Newton (IN) and least square QR factorization (LSQR) methods to tackle the nonlinearity and ill-posedness of the electromagnetic inverse scattering problem. A nonlinear model of the inverse scattering in functional form is developed. At every IN iteration, the sparse storage method is adopted to solve the storage and computational bottleneck of Fréchet derivative matrix, a large-scale sparse Jacobian matrix. Moreover, to address the slow convergence problem encountered in the inexact Newton solution via Landweber iterations, an LSQR algorithm is proposed for obtaining a better solution of the internal large-scale sparse linear equations in the IN step. Numerical results demonstrate the applicability of the proposed IN-LSQR method to quantitative inversion of scatterer electric performance parameters. Moreover, compared with the inexact Newton method based on Landweber iterations, the proposed method significantly improves the convergence rate with less computational and storage cost.


---

## Body

## 1. Introduction

Electromagnetic inverse scattering imaging is a typical inverse problem. It uses the scattering data and the wave propagation model to reconstruct the shape or the electrical parameter distribution of the unknown scatterer. The solution to this inverse scattering problem has been widely studied due to its widespread applications in various fields, including nondestructive detection [1], geophysics [2], ground penetrating radar imaging [3], and medical microwave imaging [4]. Among them, the inherent ill-posedness and nonlinearity of the inverse scattering problem are the major challenges in the inversion solution.Inversion methods can be classified into linear and nonlinear methods. For weak scattering targets, the linear method uses the Born approximation or the Rytov approximation model [5] to approximate the nonlinear relationship between the scattering field and the scatter electrical performance parameters. It results in an inaccurate reconstruction for strong scattering targets. In order to address this issue, it is necessary to directly tackle the nonlinear relationship. In this case, the nonlinear methods have been proposed. Nonlinear methods can be divided into the deterministic inversion methods and the stochastic inversion methods. Deterministic inversion methods, such as contrast source inversion (CSI) [6] and subspace optimization method (SOM) [7], characterize the nonlinear model by constructing the objective function in the form of an iterative optimization and can reconstruct the target’s electrical performance parameters well. However, they easily fall into the locally optimal solution. Stochastic inversion methods include genetic algorithms [8] and particle swarm optimization method [9], which use the concepts of biology or solid annealing to search the optimal solution of objective function or fitness function at random, without providing a priori knowledge of the target region and the number of target scatterers, however, the computational cost is enormous.As a nonlinear inversion method, inexact Newton method has been extensively utilized in microwave imaging because of its quadratic convergence rate and efficiency in reconstructing domains with strong scattering targets. In this method, the outer loop locally linearizes the nonlinear model about the scattering field and total field integral equation by the inexact Newton frame as primary step and then solves the obtained linear equations via the Landweber iteration or iterative shrinkage threshold method (the inner-loop calculation) [10, 11]. Unlike other Newton inversion methods, such as the Gauss Newton inversion method [12], the inexact Newton method does not involve the storage and computation of the large-scale Hessian matrix.Various researchers have proposed using the inexact Newton method for microwave imaging. For example, Estatico et al. [13] proposed an inversion method in which two-order Born approximation is used to construct the nonlinear model, and the ill-posedness of the internal linear equation is circumvented via the threshold Landweber iteration. However, the storage and computational cost which result from Fréchet derivative matrix in the inexact Newton method are primary issues that have to be taken into account. In [14], the Fréchet derivative matrix is described as a large-scale sparse Jacobian matrix. To this end, sparse matrix triple compression storage is applied as an efficient form for the calculation of inner large-scale sparse linear equation. In addition, the point that the Landweber iteration is actually the steepest descent method with a fixed step size and has the problem of slow convergence is demonstrated in [14–16]. To this end, considering the sparse storage form and the least square QR decomposition algorithm (LSQR) having a stable convergence property and superior convergence rate for the ill-conditioned linear equation [17], in this paper, the IN-LSQR algorithm is proposed to solve the nonlinear inverse scattering problem. It exploits the inexact Newton and the large-scale sparse characteristics of Fréchet derivative matrix.Starting from the nonlinear model of inverse scattering in functional form, this paper presents experimental results that demonstrate the efficiency of the proposed IN-LSQR method. It is shown that the convergence rate is significantly improved and the computational and storage cost is reduced for the inversion process.


## 2. Formulation

### 2.1. Nonlinear Electromagnetic Inverse Scattering Problem

The imaging geometry model of the inverse scattering problem is shown in Figure1. Let D represent the target domain with unknown target scatterers. Let εrr be the unknown permittivity of target scatterer, and in the background medium, permittivity εr=ε0. Let Ω represent the observation domain surrounded by Nr receiving antennas. The magnetic permeability of both the observation domain and the background medium is μ0. NT transmitting antennas and Nr receiving antennas are deployed with an equal space in the circular trajectories T and R.Figure 1
Geometric model of inverse scattering imaging.The incident electric field generated by theith source is Eiincr. Upon illumination by Eiincr, the equivalent electric current density Jir is induced on the target domain D, and Jir generates scattered filed Eiscar. Jir=χrEitotr, here note that the contrast function of scatterer is χr=εrr/ε0−1 and the total electric field is Eitotr=Eiincr+Eiscar for r∈D [6]. Lid. is the equation of field state, for r∈D, given as
(1)Lidχ,Ji=Jir−χrEiincr−k02χr∫DJir'GDr,r'dr'=0,r,r'∈D.Lir. is the data equation. For the scattered field Eiscar measured by the receiving antennas is given by
(2)LirJi=Eiscar=k02∫DJir'GRrq,r'dr',rq∈Ω,r'∈D.Here,GDr,r'=H02k0r−r'/4j,r,r'∈D represents the 2D Green function of target domain D. The wavenumber is k0=ωε0μ0=2π/λ0, where the angular frequency is denoted by ω and the wavelength is represented by λ0. H02 is the second kind Hankel function of zero order. GRrq,r'=H02k0rq−r'/4j,rq∈Ω,r'∈D is Green function of target domain to positions of receiving antennas.In order to solve the nonlinear inverse scattering problem given by (1) and (2), based on the idea of [14, 15], we cascade (1) and (2) in a functional form and χr can be obtained by solving
(3)Lz−y=0,where zr=χr,J1r,…,JNTrT and yr=0,…,0,E1scar,E2scar,…,ENTscarT, Lz=L1dχ,J1,…,LNTdχ,J1,L1rJ1,…,LNTrJNTT.


### 2.2. Inexact Newton Formulation

Equation (3) describes the ill-posedness and nonlinearity of electromagnetic inverse imaging, and the optimal solution of χr is obtained from the optimization scheme. The inexact Newton scheme is an efficient method to solve the nonlinear inverse problem in electromagnetic imaging. It does this by locally linearizing the nonlinear problem using the inexact Newton frame as the primary step and then solving the obtained linear systems (the inner-loop calculation) [18]. The nonlinear relation model presented as (3) can be described by Lz=y, where the inexact Newton method is given by
(4)L'zkΔzk=y−Lzk,zk+1=zk+Δzk.Here,k is the number of outer inexact Newton iterations. L'zk represents the first-order multivariate Fréchet derivative of Lz evaluated at z=zk. The solution zr is updated by zk+1=zk+Δzk, where Δzk is obtained by solving the inner ill-conditioned linear equationL'zkΔzk=y−Lzk.The Fréchet derivative matrix in (4) is a Jacobian matrix with large dimensionality. Therefore, the storage and computation cost derived from the Fréchet derivative matrix is the key problem in the inversion scheme. Since L'zk=∂Lzk/∂χ,∂Lzk/∂J1,…,∂Lzk/∂JNTT, according to (1), (2), and (3), let ∂JiLidχ,Ji, ∂χLidχ,Ji, and ∂JiLirJi represent the derivatives of Lidχ,Ji and LirJi with respect to Ji and χ, respectively. These operators are expressed as follows:
(5)∂JiLidχ,Ji=ΔJir−k02χr∫DΔJir'GDr,r'dr',r,r'∈D,(6)∂χLidχ,Ji=Δχr−Eiincr−k02∫DJir'GDr,r'dr',r,r'∈D,(7)∂JiLirJi=k02∫DΔJir'GRrq,r'dr',rq∈Ω,r'∈D.Remarkably,Lidχ,Ji only depends on Ji and χ, and LirJi only depends on Ji. Thus, ∂JmLidχ,Ji=0, ∂JmLirJi=0, and ∂χLirJi=0 for i≠m. L'zk is a sparse matrix with nonzero elements L'zki,1=∂χLidχ,Ji∣zk, L'zki,i+1=∂JiLidχ,Ji∣zk, L'zki+NT,i+1=∂JiLirJi∣zk.To discretize (5), (6), and (7), the target domain D is divided into Nd square grid cells. The centres of the square cells and the locations of the receivers are expressed by rnd,n=1,…,Nd and rmr,m=1,…,Nr, respectively. To promote clarity, The centres of the square cells are represented by rpd,p=1,…,Nd. The Lzk and nonzero elements of L'z∣zk are
(8)Lzki−1Nd+p=Jik−DχkEiinc−DχkGDJikp,(9)LzkNTNd+i−1Nr+m=GRJikm,(10)L'zzki−1Nd+n,n=D−Eiinc−GDJikn,n,(11)L'zzki−1Nd+p,iNd+n=I−DχkGDp,n,(12)L'zzkNTNd+i−1Nr+m,iNd+n=GRm,n.where D. represents the diagonalization operation, and I is a unit matrix. Equations (8), (9), (10), (11), and (12) describe the large-scale sparse structure of the Fréchet derivative matrix L'zk. In this paper, the sparse matrix triple compression storage is used for the Fréchet derivative matrix, and only nonzero elements in (10), (11), and (12) and their row and column indexes are stored, solving the storage cost deriving from the Fréchet derivative matrix in the framework of the inexact Newton method.


### 2.3. Inexact Newton Iterations with LSQR

The internal linear equation in (4) is a large-scale sparse system. The solution Δzk can be obtained by the LSQR algorithm. Thus, the problem of solving the above internal linear equation is transformed into an optimization problem, given as
(13)Δzk=minΔzkL'zkΔzk−y−Lzk2.LSQR is an iterative regularization method, and the number of iterationsNin is the regularization parameter. Implementation of this algorithm is carried out in two steps.First, Lanczos bidiagonalization ofL'zk, namely, L'zk=UBiVT, obtains the orthogonal matrix Ui=u1,u2,…,ui, Vi=v1,v2,…,vi and bidiagonal matrix Bi. The form of Bi is shown as
(14)Bi=α1β2α2β3⋱⋱αiβi+1.Lanczos iteration begins at the orthogonal basis vectorsui,vi and the positive constant αi,βi.
(15)β1μ1=y−Lzkβ1=y−Lzkα1ν1=L'zkTμ1α1=L'zkTu1.The iterationi of this method is given as
(16)βi+1=L'zkνi−αiμi,αi+1=L'zkTμi+1−βi+1νi,βi+1μi+1=L'zkνi−αiμi,αi+1νi+1=L'zkTμi+1−βi+1νi,i=1,2,…,Nin,where αi≥0, βi≥0, ui=vi=1. Equation (16) can also be expressed as
(17)Ui+1β1e1=y−Lzk,L'zkVi=Ui+1Bi,L'zkTUi+1=ViBiT+αi+1vi+1ei+1T,where ei+1T=0,0,…,1T is the i+1 row of the n-dimensional identity matrix. Denoting the solution we need as Δzki=Viki, the residual vector in iteration i is ri=y−Lzk−L'zkΔzki. The residual vector evaluates at ti+1=β1e1−Biki, so one can obtain the residual as
(18)ri=y−Lzk−L'zkΔzki=Ui+1β1e1−L'zkViki=Ui+1β1e1−Ui+1Biki=Ui+1ti+1.Second, solve the least square problem after bidiagonalization. From the residual described in (18), the solution of the internal sparse linear equation in (4) is transformed into the least square solution of minkiβ1e1−Biki2. Since Bi is a bidiagonal matrix, it can be easily solved by the QR decomposition method, and then one can obtain Δzki according to Δzki=Viki. The proposed IN-LSQR algorithm is shown in Algorithm 1.Algorithm 1: Algorithm flow of IN-LSQR.
Input: number of transmitting antennas NT; number of receiving antennas Nr; the number of square cells Nd; scattered field data measured Esca
Output: solution of (3) zr
Initialization: initial solution z0r=0
Outer inexact Newton iterations fork=1,2,3,…
CalculateLzk on the basis of (8) and (9), and calculate L'zk on the basis of (10), (11), and (12) and stored in sparse compressed form.
1.1 Internal initializationβ1μ1=y−Lzk, α1ν1=L'zkTu1, ω1=ν1, Δzk0=0, ϕ¯1=β1=y−Lzk, ρ¯1=α1=L'zkTμ1. ϕ¯1,ρ¯1,α1,β1 is a nonnegative real number
Internal iterations fori=1,2,3,…,Nin
1.2 Lanczos bidiagonalization ofL'zk on the basis of (8), L'zk=UBiVT gives orthogonal matrix Ui=u1,u2,…,ui, Vi=v1,v2,…,vi and bidiagonal matrix Bi.
1.3 Calculate the intermediate variable of the QR decomposition process in (18),
ρi+1=ρ¯2i+β2i+11/2,
ci=ρ¯iρi,
si=βi+1ρi,
θi+1=siαi+1,
ϕi=ciϕ¯i,
ρ¯i+1=−ciαi+1¯,
ϕ¯i+1=siϕ¯i,
whereρ¯i+1 and ϕ¯i+1 are the initial value of the next iteration.
1.4 UpdateΔzki and ω, according to
Δzki=Δzki−1+ϕiρiωi,
ωi+1=νi+1−θi+1ρiωi.
1.5 Ifi<Nin, then i=i+1 and repeat steps 1.2–1.4. Otherwise, terminate the internal iteration, achieving Δzk, the solution to (13).
2 Update the solution to (3) using zk+1=zk+Δzk
3 When the change of zr in two consecutive iterations is less than the termination condition, terminate the outer inexact Newton iteration, else k=k+1 and return to step 1.The internal LSQR iteration obtains the solutionΔzk, and the external update zk+1=zk+Δzk is applied to get the solution to the nonlinear inverse scattering problem.Since the number of grid cells holds thatNd≫Nr and Nd≫NT, compared with the Landweber iterative method, the LSQR method does not require the singular value decomposition of Fréchet derivative matrix with computational complexity of ONd3. The LSQR algorithm mainly calculates matrix-vector multiplication, L'zkν or L'zkTu. Thus, the computational complexity can be reduced to ONd2. Furthermore, the proposed IN-LSQR method uses the sparse matrix triple compression form to store Fréchet derivative matrix, reducing the storage cost of algorithm in practical situation. The large-scale sparse structure of Fréchet derivative matrix does not change in the iterative process, and matrix-vector multiplication in IN-LSQR only involves the nonzero elements of L'zk. Consequently, the computational efficiency is significantly improved by use of sparse characteristics in L'zk.


## 2.1. Nonlinear Electromagnetic Inverse Scattering Problem

The imaging geometry model of the inverse scattering problem is shown in Figure1. Let D represent the target domain with unknown target scatterers. Let εrr be the unknown permittivity of target scatterer, and in the background medium, permittivity εr=ε0. Let Ω represent the observation domain surrounded by Nr receiving antennas. The magnetic permeability of both the observation domain and the background medium is μ0. NT transmitting antennas and Nr receiving antennas are deployed with an equal space in the circular trajectories T and R.Figure 1
Geometric model of inverse scattering imaging.The incident electric field generated by theith source is Eiincr. Upon illumination by Eiincr, the equivalent electric current density Jir is induced on the target domain D, and Jir generates scattered filed Eiscar. Jir=χrEitotr, here note that the contrast function of scatterer is χr=εrr/ε0−1 and the total electric field is Eitotr=Eiincr+Eiscar for r∈D [6]. Lid. is the equation of field state, for r∈D, given as
(1)Lidχ,Ji=Jir−χrEiincr−k02χr∫DJir'GDr,r'dr'=0,r,r'∈D.Lir. is the data equation. For the scattered field Eiscar measured by the receiving antennas is given by
(2)LirJi=Eiscar=k02∫DJir'GRrq,r'dr',rq∈Ω,r'∈D.Here,GDr,r'=H02k0r−r'/4j,r,r'∈D represents the 2D Green function of target domain D. The wavenumber is k0=ωε0μ0=2π/λ0, where the angular frequency is denoted by ω and the wavelength is represented by λ0. H02 is the second kind Hankel function of zero order. GRrq,r'=H02k0rq−r'/4j,rq∈Ω,r'∈D is Green function of target domain to positions of receiving antennas.In order to solve the nonlinear inverse scattering problem given by (1) and (2), based on the idea of [14, 15], we cascade (1) and (2) in a functional form and χr can be obtained by solving
(3)Lz−y=0,where zr=χr,J1r,…,JNTrT and yr=0,…,0,E1scar,E2scar,…,ENTscarT, Lz=L1dχ,J1,…,LNTdχ,J1,L1rJ1,…,LNTrJNTT.


## 2.2. Inexact Newton Formulation

Equation (3) describes the ill-posedness and nonlinearity of electromagnetic inverse imaging, and the optimal solution of χr is obtained from the optimization scheme. The inexact Newton scheme is an efficient method to solve the nonlinear inverse problem in electromagnetic imaging. It does this by locally linearizing the nonlinear problem using the inexact Newton frame as the primary step and then solving the obtained linear systems (the inner-loop calculation) [18]. The nonlinear relation model presented as (3) can be described by Lz=y, where the inexact Newton method is given by
(4)L'zkΔzk=y−Lzk,zk+1=zk+Δzk.Here,k is the number of outer inexact Newton iterations. L'zk represents the first-order multivariate Fréchet derivative of Lz evaluated at z=zk. The solution zr is updated by zk+1=zk+Δzk, where Δzk is obtained by solving the inner ill-conditioned linear equationL'zkΔzk=y−Lzk.The Fréchet derivative matrix in (4) is a Jacobian matrix with large dimensionality. Therefore, the storage and computation cost derived from the Fréchet derivative matrix is the key problem in the inversion scheme. Since L'zk=∂Lzk/∂χ,∂Lzk/∂J1,…,∂Lzk/∂JNTT, according to (1), (2), and (3), let ∂JiLidχ,Ji, ∂χLidχ,Ji, and ∂JiLirJi represent the derivatives of Lidχ,Ji and LirJi with respect to Ji and χ, respectively. These operators are expressed as follows:
(5)∂JiLidχ,Ji=ΔJir−k02χr∫DΔJir'GDr,r'dr',r,r'∈D,(6)∂χLidχ,Ji=Δχr−Eiincr−k02∫DJir'GDr,r'dr',r,r'∈D,(7)∂JiLirJi=k02∫DΔJir'GRrq,r'dr',rq∈Ω,r'∈D.Remarkably,Lidχ,Ji only depends on Ji and χ, and LirJi only depends on Ji. Thus, ∂JmLidχ,Ji=0, ∂JmLirJi=0, and ∂χLirJi=0 for i≠m. L'zk is a sparse matrix with nonzero elements L'zki,1=∂χLidχ,Ji∣zk, L'zki,i+1=∂JiLidχ,Ji∣zk, L'zki+NT,i+1=∂JiLirJi∣zk.To discretize (5), (6), and (7), the target domain D is divided into Nd square grid cells. The centres of the square cells and the locations of the receivers are expressed by rnd,n=1,…,Nd and rmr,m=1,…,Nr, respectively. To promote clarity, The centres of the square cells are represented by rpd,p=1,…,Nd. The Lzk and nonzero elements of L'z∣zk are
(8)Lzki−1Nd+p=Jik−DχkEiinc−DχkGDJikp,(9)LzkNTNd+i−1Nr+m=GRJikm,(10)L'zzki−1Nd+n,n=D−Eiinc−GDJikn,n,(11)L'zzki−1Nd+p,iNd+n=I−DχkGDp,n,(12)L'zzkNTNd+i−1Nr+m,iNd+n=GRm,n.where D. represents the diagonalization operation, and I is a unit matrix. Equations (8), (9), (10), (11), and (12) describe the large-scale sparse structure of the Fréchet derivative matrix L'zk. In this paper, the sparse matrix triple compression storage is used for the Fréchet derivative matrix, and only nonzero elements in (10), (11), and (12) and their row and column indexes are stored, solving the storage cost deriving from the Fréchet derivative matrix in the framework of the inexact Newton method.


## 2.3. Inexact Newton Iterations with LSQR

The internal linear equation in (4) is a large-scale sparse system. The solution Δzk can be obtained by the LSQR algorithm. Thus, the problem of solving the above internal linear equation is transformed into an optimization problem, given as
(13)Δzk=minΔzkL'zkΔzk−y−Lzk2.LSQR is an iterative regularization method, and the number of iterationsNin is the regularization parameter. Implementation of this algorithm is carried out in two steps.First, Lanczos bidiagonalization ofL'zk, namely, L'zk=UBiVT, obtains the orthogonal matrix Ui=u1,u2,…,ui, Vi=v1,v2,…,vi and bidiagonal matrix Bi. The form of Bi is shown as
(14)Bi=α1β2α2β3⋱⋱αiβi+1.Lanczos iteration begins at the orthogonal basis vectorsui,vi and the positive constant αi,βi.
(15)β1μ1=y−Lzkβ1=y−Lzkα1ν1=L'zkTμ1α1=L'zkTu1.The iterationi of this method is given as
(16)βi+1=L'zkνi−αiμi,αi+1=L'zkTμi+1−βi+1νi,βi+1μi+1=L'zkνi−αiμi,αi+1νi+1=L'zkTμi+1−βi+1νi,i=1,2,…,Nin,where αi≥0, βi≥0, ui=vi=1. Equation (16) can also be expressed as
(17)Ui+1β1e1=y−Lzk,L'zkVi=Ui+1Bi,L'zkTUi+1=ViBiT+αi+1vi+1ei+1T,where ei+1T=0,0,…,1T is the i+1 row of the n-dimensional identity matrix. Denoting the solution we need as Δzki=Viki, the residual vector in iteration i is ri=y−Lzk−L'zkΔzki. The residual vector evaluates at ti+1=β1e1−Biki, so one can obtain the residual as
(18)ri=y−Lzk−L'zkΔzki=Ui+1β1e1−L'zkViki=Ui+1β1e1−Ui+1Biki=Ui+1ti+1.Second, solve the least square problem after bidiagonalization. From the residual described in (18), the solution of the internal sparse linear equation in (4) is transformed into the least square solution of minkiβ1e1−Biki2. Since Bi is a bidiagonal matrix, it can be easily solved by the QR decomposition method, and then one can obtain Δzki according to Δzki=Viki. The proposed IN-LSQR algorithm is shown in Algorithm 1.Algorithm 1: Algorithm flow of IN-LSQR.
Input: number of transmitting antennas NT; number of receiving antennas Nr; the number of square cells Nd; scattered field data measured Esca
Output: solution of (3) zr
Initialization: initial solution z0r=0
Outer inexact Newton iterations fork=1,2,3,…
CalculateLzk on the basis of (8) and (9), and calculate L'zk on the basis of (10), (11), and (12) and stored in sparse compressed form.
1.1 Internal initializationβ1μ1=y−Lzk, α1ν1=L'zkTu1, ω1=ν1, Δzk0=0, ϕ¯1=β1=y−Lzk, ρ¯1=α1=L'zkTμ1. ϕ¯1,ρ¯1,α1,β1 is a nonnegative real number
Internal iterations fori=1,2,3,…,Nin
1.2 Lanczos bidiagonalization ofL'zk on the basis of (8), L'zk=UBiVT gives orthogonal matrix Ui=u1,u2,…,ui, Vi=v1,v2,…,vi and bidiagonal matrix Bi.
1.3 Calculate the intermediate variable of the QR decomposition process in (18),
ρi+1=ρ¯2i+β2i+11/2,
ci=ρ¯iρi,
si=βi+1ρi,
θi+1=siαi+1,
ϕi=ciϕ¯i,
ρ¯i+1=−ciαi+1¯,
ϕ¯i+1=siϕ¯i,
whereρ¯i+1 and ϕ¯i+1 are the initial value of the next iteration.
1.4 UpdateΔzki and ω, according to
Δzki=Δzki−1+ϕiρiωi,
ωi+1=νi+1−θi+1ρiωi.
1.5 Ifi<Nin, then i=i+1 and repeat steps 1.2–1.4. Otherwise, terminate the internal iteration, achieving Δzk, the solution to (13).
2 Update the solution to (3) using zk+1=zk+Δzk
3 When the change of zr in two consecutive iterations is less than the termination condition, terminate the outer inexact Newton iteration, else k=k+1 and return to step 1.The internal LSQR iteration obtains the solutionΔzk, and the external update zk+1=zk+Δzk is applied to get the solution to the nonlinear inverse scattering problem.Since the number of grid cells holds thatNd≫Nr and Nd≫NT, compared with the Landweber iterative method, the LSQR method does not require the singular value decomposition of Fréchet derivative matrix with computational complexity of ONd3. The LSQR algorithm mainly calculates matrix-vector multiplication, L'zkν or L'zkTu. Thus, the computational complexity can be reduced to ONd2. Furthermore, the proposed IN-LSQR method uses the sparse matrix triple compression form to store Fréchet derivative matrix, reducing the storage cost of algorithm in practical situation. The large-scale sparse structure of Fréchet derivative matrix does not change in the iterative process, and matrix-vector multiplication in IN-LSQR only involves the nonzero elements of L'zk. Consequently, the computational efficiency is significantly improved by use of sparse characteristics in L'zk.


## 3. Experimental Results

In this section, the efficiency of our proposed IN-LSQR method is demonstrated via three sets of scattered field measurements from Fresnel laboratory. The experimental data can be found in the files “dielTM_dec8f,” “twodielTM_8f,” and “FoamTwinDielTM” [19, 20].


### 3.1. dielTM and twodielTM

The actual profile in the experiment is shown in Figure2. Figure 2(a) shows the scene of dielTM. The single cylindrical scatterer is located on the right side and from the center to the origin is 30 mm. Figure 2(b) shows the scene of twodielTM, in which two cylindrical scatterers are located on both sides of the origin, and the distance from the center of scatterer to the origin is 45 mm. In dielTM and twodielTM, the scatterers have a radius of 15 mm and a dielectric permittivity of 3 ± 0.3. The transmitter-receiver configuration is described in [19]. Free space is assumed as the background medium; hence, the contrast of scatterer χ=εr−ε0/ε0=εr−1. In this paper, we adopt the measured scattered field data at frequency f = 2GHz. λ is the wavelength of the electromagnetic wave, and the target domain of size λ×λ is divided into square cells of size λ/40×λ/40.Figure 2
Actual profile of the experiments: (a) dielTM; (b) twodielTM.
(a)
(b)The large-scale sparse characteristics of Fréchet derivative matrix was described above. Table1 shows the scale and nonzero elements of Fréchet derivative matrix, where the scattered field of dielTM has a signal to noise ratio of 20 dB. Thus, from Table 1, the sparse matrix triple compression storage can work as expected, avoiding the storage cost of massive zero elements.Table 1
Proportion of nonzero elements in matrixL'zk.
Number of cells in target domain (Nd)
Size ofL'zk
Number of nonzero elements inL'zk
Proportion of nonzero elements
20 × 20
16164 × 14800
6,480,000
0.02709
30 × 30
34164 × 33300
30,780,000
0.02706
40 × 40
59364 × 59200
95,040,000
0.02704To enable comparison of the inversion imaging quality of the proposed IN-LSQR method, Figure3 plots the contrast profile for the experiments involving dielTM and twodielTM, which were recovered by two methods: (1) the inexact Newton method based on Landweber iteration (inexact Newton-Landweber, IN-LW) and (2) the proposed IN-LSQR method. Here, the measured scattered field Esca was synthetically generated with 20 dB noise. The inner iterations Nin=40 for both IN-LSQR and IN-LW, the outer iteration k=40. Figures 3(a) and 3(c) show the results recovered by IN-LW; Figures 3(b) and 3(d) show the results recovered by IN-LSQR. As shown in Figure 3, the LN-LSQR method proposed in this paper effectively reconstructed the contrast profile of the scatterers in the experiments. Compared with the IN-LW method, as can be seen in Figures 3(c) and 3(d), the IN-LSQR method has better results than the IN-LD method in the contrast of scatterers. The results demonstrate the efficiency and accuracy of the proposed IN-LSQR method, where the location, number of scatterers, and the contrast profile are clearly identified.Figure 3
The reconstructed contrast profile in dielTM and twodielTM: (a) dielTM using IN-LW; (b) dielTM using IN-LSQR; (c) twodielTM using IN-LW; (d) twodielTM using IN-LSQR.
(a)
(b)
(c)
(d)It is worth discussing the convergence rate of the methods in the presence of noise. Here, the quality of reconstructed profile of the contrast function is evaluated in terms of the relative norm errorerr:
(19)err=χk,com−χref2χref2,where χref represents the actual contrast in the target domain, and χk,com is the reconstructed contrast result χr at every inexact Newton iteration k. Esca in the experiments with dielTM and twodielTM was generated with 20 dB, 10 dB, and 5 dB Gaussian noise. Here, the internal iteration Nin=40 for both IN-LSQR and IN-LW. Figure 4 shows errk computed by the IN-LSQR and the IN-LW for three levels of noise. Moreover, Table 2 lists the computation time and contrast error of IN-LW and IN-LSQR methods at k=1,5,10,15 in the experiments of dielTM and twodielTM with 20 dB noise. The computer used had an Intel E5-2620@2.1 GHZ CPU, 64 GB RAM and MATLAB 2014b. According to Figure 4 and Table 2, errk computed by IN-LSQR converges faster than the IN-LW, and errk in IN-LSQR converges to a smaller value than that of IN-LW at different noise levels. It means the IN-LSQR can get reconstruction result with less number of iterations in different levels of noise. Combining the results of Table 2, the actual computation time of IN-LSQR is less than that of IN-LW.Figure 4
Relationship between contrast error and iterations: (a)errk computed by IN-LSQR and IN-LW for different noise levels in dielTM; (b) errk computed by IN-LSQR and IN-LW for different noise levels in twodielTM.
(a)
(b)Table 2
Comparison of IN-LW and IN-LSQR methods in terms of computation time and contrast error atk=1,5,10,15.
Experiment
Method
Iterationk=1
Iterationk=5
Iterationk=10
Iterationk=15
Time (s)
err
Time (s)
err
Time (s)
err
Time (s)
err
dielTM in Figure2(a)
IN-LW
427.950
0.971
3.725 × 103
0.893
9.784 × 103
0.823
1.098 × 104
0.768
IN-LSQR
35.236
0.898
164.792
0.668
319.208
0.609
470.326
0.570
twodielTM in Figure2(b)
IN-LW
445.85
0.962
3.824 × 103
0.876
9.809 × 103
0.801
1.177 × 104
0.758
IN-LSQR
35.558
0.915
169.672
0.693
330.833
0.644
488.199
0.611


### 3.2. FoamTwinDiel

The proposed IN-LSQR was also tested using the experimental data of FoamTwinDiel [20]. The FoamTwinDiel profile includes two smaller cylinders of permittivity εr=3±0.3 with diameter being 31 mm, where one of the cylinders is embedded in a larger cylinder with εr=1.45±0.15 and the diameter is 80 mm. This data set comprises measurements from 18 transmitters and 241 receivers. The frequency, size of target domain, and sample cells are the same as those in dielTM and twodielTM. The noise level in the scattered field data changes from 20 dB, 10 dB, and 5 dB. Here, the internal iteration is Nin=20 for both IN-LSQR and IN-LW.Figure5(a) shows the actual profile of FoamTwinDiel. Figure 5(b) shows plots of errk computed by IN-LSQR and IN-LW for noise levels of 20 dB, 10 dB, and 5 dB. This figure illustrates that IN-LSQR converges faster and is more immune to noise than IN-LW. Moreover, Figure 5(c) plots errk computed by IN-LSQR and CSI for different noise levels. As can be seen, errk in CSI converges to a smaller value in the presence of 20 dB noise; however, IN-LSQR converges faster and performs better in 10 dB and 5 dB noise. Figures 5(d)–5(f) show plots of the contrast profile recovered by IN-LSQR, IN-LW, and CSI with 20 dB noise at iteration k=50, respectively. The profile result reconstructed by IN-LSQR, similar to CSI, is clearly more accurate than that of IN-LW.Figure 5
FoamTwinDiel: (a) actual profile of FoamTwinDiel; (b)errk computed by IN-LSQR and IN-LW for different noise levels; (c) errk computed by IN-LSQR and CSI for different noise levels; (d) contrast profile recovered by IN-LSQR for 20 dB noise; (e) contrast profile recovered by IN-LW for 20 dB noise; (f) contrast profile recovered by CSI for 20 dB noise.
(a)
(b)
(c)
(d)
(e)
(f)


## 3.1. dielTM and twodielTM

The actual profile in the experiment is shown in Figure2. Figure 2(a) shows the scene of dielTM. The single cylindrical scatterer is located on the right side and from the center to the origin is 30 mm. Figure 2(b) shows the scene of twodielTM, in which two cylindrical scatterers are located on both sides of the origin, and the distance from the center of scatterer to the origin is 45 mm. In dielTM and twodielTM, the scatterers have a radius of 15 mm and a dielectric permittivity of 3 ± 0.3. The transmitter-receiver configuration is described in [19]. Free space is assumed as the background medium; hence, the contrast of scatterer χ=εr−ε0/ε0=εr−1. In this paper, we adopt the measured scattered field data at frequency f = 2GHz. λ is the wavelength of the electromagnetic wave, and the target domain of size λ×λ is divided into square cells of size λ/40×λ/40.Figure 2
Actual profile of the experiments: (a) dielTM; (b) twodielTM.
(a)
(b)The large-scale sparse characteristics of Fréchet derivative matrix was described above. Table1 shows the scale and nonzero elements of Fréchet derivative matrix, where the scattered field of dielTM has a signal to noise ratio of 20 dB. Thus, from Table 1, the sparse matrix triple compression storage can work as expected, avoiding the storage cost of massive zero elements.Table 1
Proportion of nonzero elements in matrixL'zk.
Number of cells in target domain (Nd)
Size ofL'zk
Number of nonzero elements inL'zk
Proportion of nonzero elements
20 × 20
16164 × 14800
6,480,000
0.02709
30 × 30
34164 × 33300
30,780,000
0.02706
40 × 40
59364 × 59200
95,040,000
0.02704To enable comparison of the inversion imaging quality of the proposed IN-LSQR method, Figure3 plots the contrast profile for the experiments involving dielTM and twodielTM, which were recovered by two methods: (1) the inexact Newton method based on Landweber iteration (inexact Newton-Landweber, IN-LW) and (2) the proposed IN-LSQR method. Here, the measured scattered field Esca was synthetically generated with 20 dB noise. The inner iterations Nin=40 for both IN-LSQR and IN-LW, the outer iteration k=40. Figures 3(a) and 3(c) show the results recovered by IN-LW; Figures 3(b) and 3(d) show the results recovered by IN-LSQR. As shown in Figure 3, the LN-LSQR method proposed in this paper effectively reconstructed the contrast profile of the scatterers in the experiments. Compared with the IN-LW method, as can be seen in Figures 3(c) and 3(d), the IN-LSQR method has better results than the IN-LD method in the contrast of scatterers. The results demonstrate the efficiency and accuracy of the proposed IN-LSQR method, where the location, number of scatterers, and the contrast profile are clearly identified.Figure 3
The reconstructed contrast profile in dielTM and twodielTM: (a) dielTM using IN-LW; (b) dielTM using IN-LSQR; (c) twodielTM using IN-LW; (d) twodielTM using IN-LSQR.
(a)
(b)
(c)
(d)It is worth discussing the convergence rate of the methods in the presence of noise. Here, the quality of reconstructed profile of the contrast function is evaluated in terms of the relative norm errorerr:
(19)err=χk,com−χref2χref2,where χref represents the actual contrast in the target domain, and χk,com is the reconstructed contrast result χr at every inexact Newton iteration k. Esca in the experiments with dielTM and twodielTM was generated with 20 dB, 10 dB, and 5 dB Gaussian noise. Here, the internal iteration Nin=40 for both IN-LSQR and IN-LW. Figure 4 shows errk computed by the IN-LSQR and the IN-LW for three levels of noise. Moreover, Table 2 lists the computation time and contrast error of IN-LW and IN-LSQR methods at k=1,5,10,15 in the experiments of dielTM and twodielTM with 20 dB noise. The computer used had an Intel E5-2620@2.1 GHZ CPU, 64 GB RAM and MATLAB 2014b. According to Figure 4 and Table 2, errk computed by IN-LSQR converges faster than the IN-LW, and errk in IN-LSQR converges to a smaller value than that of IN-LW at different noise levels. It means the IN-LSQR can get reconstruction result with less number of iterations in different levels of noise. Combining the results of Table 2, the actual computation time of IN-LSQR is less than that of IN-LW.Figure 4
Relationship between contrast error and iterations: (a)errk computed by IN-LSQR and IN-LW for different noise levels in dielTM; (b) errk computed by IN-LSQR and IN-LW for different noise levels in twodielTM.
(a)
(b)Table 2
Comparison of IN-LW and IN-LSQR methods in terms of computation time and contrast error atk=1,5,10,15.
Experiment
Method
Iterationk=1
Iterationk=5
Iterationk=10
Iterationk=15
Time (s)
err
Time (s)
err
Time (s)
err
Time (s)
err
dielTM in Figure2(a)
IN-LW
427.950
0.971
3.725 × 103
0.893
9.784 × 103
0.823
1.098 × 104
0.768
IN-LSQR
35.236
0.898
164.792
0.668
319.208
0.609
470.326
0.570
twodielTM in Figure2(b)
IN-LW
445.85
0.962
3.824 × 103
0.876
9.809 × 103
0.801
1.177 × 104
0.758
IN-LSQR
35.558
0.915
169.672
0.693
330.833
0.644
488.199
0.611


## 3.2. FoamTwinDiel

The proposed IN-LSQR was also tested using the experimental data of FoamTwinDiel [20]. The FoamTwinDiel profile includes two smaller cylinders of permittivity εr=3±0.3 with diameter being 31 mm, where one of the cylinders is embedded in a larger cylinder with εr=1.45±0.15 and the diameter is 80 mm. This data set comprises measurements from 18 transmitters and 241 receivers. The frequency, size of target domain, and sample cells are the same as those in dielTM and twodielTM. The noise level in the scattered field data changes from 20 dB, 10 dB, and 5 dB. Here, the internal iteration is Nin=20 for both IN-LSQR and IN-LW.Figure5(a) shows the actual profile of FoamTwinDiel. Figure 5(b) shows plots of errk computed by IN-LSQR and IN-LW for noise levels of 20 dB, 10 dB, and 5 dB. This figure illustrates that IN-LSQR converges faster and is more immune to noise than IN-LW. Moreover, Figure 5(c) plots errk computed by IN-LSQR and CSI for different noise levels. As can be seen, errk in CSI converges to a smaller value in the presence of 20 dB noise; however, IN-LSQR converges faster and performs better in 10 dB and 5 dB noise. Figures 5(d)–5(f) show plots of the contrast profile recovered by IN-LSQR, IN-LW, and CSI with 20 dB noise at iteration k=50, respectively. The profile result reconstructed by IN-LSQR, similar to CSI, is clearly more accurate than that of IN-LW.Figure 5
FoamTwinDiel: (a) actual profile of FoamTwinDiel; (b)errk computed by IN-LSQR and IN-LW for different noise levels; (c) errk computed by IN-LSQR and CSI for different noise levels; (d) contrast profile recovered by IN-LSQR for 20 dB noise; (e) contrast profile recovered by IN-LW for 20 dB noise; (f) contrast profile recovered by CSI for 20 dB noise.
(a)
(b)
(c)
(d)
(e)
(f)


## 4. Conclusions

In this paper, considering the nonlinear relation of electromagnetic inverse scattering imaging, an IN-LSQR algorithm was proposed. The proposed algorithm utilizes sparse matrix triple compression storage for the large-scale sparse Fréchet derivative matrix and LSQR for stable solution of the ill-conditioned linear equation. Numerical results demonstrated the efficiency of IN-LSQR in reconstructing contrast profile of target scatterers with different noise levels. Further, the proposed IN-LSQR improves the convergence rate with lower computational and storage cost.



---

*Source: 2794646-2018-08-02.xml*