Nonlinear Electromagnetic Inverse Scattering Imaging Based on IN-LSQR

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.


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 twoorder 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][15][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.

Nonlinear Electromagnetic Inverse Scattering Problem.
The imaging geometry model of the inverse scattering problem is shown in Figure 1.Let D represent the target domain with unknown target scatterers.Let ε r r be the unknown permittivity of target scatterer, and in the background medium, permittivity ε r = ε 0 .Let Ω represent the observation domain surrounded by N r receiving antennas.The magnetic permeability of both the observation domain and the background medium is μ 0 .N T transmitting antennas and N r receiving antennas are deployed with an equal space in the circular trajectories T and R.
The incident electric field generated by the ith source is E inc i r .Upon illumination by E inc i r , the equivalent electric current density J i r is induced on the target domain D, and J i r generates scattered filed E sca i r .J i r = χ r E tot i r , here note that the contrast function of scatterer is χ r = ε r r /ε 0 − 1 and the total electric field is is the equation of field state, for r ∈ D, given as 2 International Journal of Antennas and Propagation Here, G D r, r = H 2 0 k 0 r − r /4j, r, r ∈ D represents the 2D Green function of target domain D. The wavenumber is k 0 = ω ε 0 μ 0 = 2π/λ 0 , where the angular frequency is denoted by ω and the wavelength is represented by λ 0 .
is the second kind Hankel function of zero order.G R r q , r = H 2 0 k 0 r q − r /4j, r q ∈ Ω, 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 where z r = χ r , J 1 r , … , J N T r T and y r = 0, … , 0, 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 L z = y, where the inexact Newton method is given by Here, k is the number of outer inexact Newton iterations.L z k represents the first-order multivariate Fréchet derivative of L z evaluated at z = z k .The solution z r is updated by 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 T , according to (1), (2), and (3), let i J i represent the derivatives of L d i χ, J i and L r i J i with respect to J i and χ, respectively.These operators are expressed as follows: Remarkably, L d i χ, J i only depends on J i and χ, and L r i J i only depends on J i .Thus, 5), (6), and ( 7 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 z k .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.

Inexact
Newton Iterations with LSQR.The internal linear equation in ( 4) is a large-scale sparse system.The solution Δz k 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 International Journal of Antennas and Propagation LSQR is an iterative regularization method, and the number of iterations N in is the regularization parameter.Implementation of this algorithm is carried out in two steps.
First, Lanczos bidiagonalization of L z k , namely, L z k = UB i V T , obtains the orthogonal matrix Lanczos iteration begins at the orthogonal basis vectors u i , v i and the positive constant α i , β i .
The iteration i of this method is given as where 16) can also be expressed as Input: number of transmitting antennas N T ; number of receiving antennas N r ; the number of square cells N d ; scattered field data measured E sca Output: solution of (3) z r Initialization: initial solution z 0 r = 0 Outer inexact Newton iterations for k = 1, 2, 3, … Calculate L z k on the basis of ( 8) and ( 9), and calculate L z k on the basis of ( 10), (11), and (12) and stored in sparse compressed form.1.1 Internal initialization  (18), where ρ i+1 and ϕ i+1 are the initial value of the next iteration.where e T i+1 = 0, 0, … , 1 T is the i + 1 row of the ndimensional identity matrix.Denoting the solution we need as Δz

Update Δz
The residual vector evaluates at t i+1 = β 1 e 1 − B i k i , so one can obtain the residual as 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 min is a bidiagonal matrix, it can be easily solved by the QR decomposition method, and then one can obtain Δz i k according to Δz The internal LSQR iteration obtains the solution Δz k , and the external update z k+1 = z k + Δz k is applied to get the solution to the nonlinear inverse scattering problem.
Since the number of grid cells holds that N d ≫ N r and N d ≫ N T , compared with the Landweber iterative method, the LSQR does not require the singular decomposition of Fréchet derivative matrix with computational complexity of O N d 3 .The LSQR algorithm mainly calculates matrix-vector multiplication, L z k ν or L z k T u.Thus, the computational complexity can be reduced to O N d 2 .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 z k .Consequently, the computational efficiency is significantly improved by use of sparse characteristics in L z k .

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 Figure 2. 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.
The large-scale sparse characteristics of Fréchet derivative matrix was described above.Table 1 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.
To enable comparison of the inversion imaging quality of the proposed IN-LSQR method, Figure 3 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 E sca was synthetically generated with 20 dB noise.The inner iterations      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 error err: 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.E sca in the experiments with dielTM and twodielTM was generated with 20 dB,   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 N in = 20 for both IN-LSQR and IN-LW.
Figure 5(a) shows the actual profile of FoamTwinDiel.

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 illconditioned 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.

Figure 1 :
Figure 1: Geometric model of inverse scattering imaging.
), the target domain D is divided into N d square grid cells.The centres of the square cells and the locations of the receivers are expressed by r d n , n = 1, … , N d and r r m , m = 1, … , N r , respectively.To promote clarity, The centres of the square cells are represented by r d p , p = 1, … , N d .The L z k and nonzero ele-

1 Algorithm 1 :
and repeat steps 1.2-1.4.Otherwise, terminate the internal iteration, achieving Δz k , the solution to(13). 2 Update the solution to (3) using z k+1 = z k + Δz k 3 When the change of z r 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. Algorithm flow of IN-LSQR. 4 International Journal of Antennas and Propagation

Figure 4 :
Figure 4: Relationship between contrast error and iterations: (a) err k computed by IN-LSQR and IN-LW different noise levels in dielTM; (b) err k computed by IN-LSQR and IN-LW for different noise levels in twodielTM.

6
International Journal of Antennas and Propagation N in = 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 IN-LSQR method, where the location, number of scatterers, and the contrast profile are clearly identified.

Figure 5 :
Figure 5: FoamTwinDiel: (a) actual profile of FoamTwinDiel; (b) err k computed by IN-LSQR and IN-LW for different noise levels; (c) err k 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.

Figure 5 (
b) shows plots of err k 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 err k computed by IN-LSQR and CSI for different noise levels.As can be seen, err k 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.
… , v i and bidiagonal matrix B i .1.3 Calculate the intermediate variable of the QR decomposition process in

Table 1 :
Proportion of nonzero elements in matrix L z k .Number of cells in target domain (N d ) Size of L z k Number of nonzero elements in L z k

Table 2 :
Comparison of IN-LW and IN-LSQR methods in terms of computation time and contrast error at k = 1, 5, 10, 15.