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.

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 [1416]. 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 Figure 1. Let D represent the target domain with unknown target scatterers. Let be the unknown permittivity of target scatterer, and in the background medium, permittivity . Let represent the observation domain surrounded by receiving antennas. The magnetic permeability of both the observation domain and the background medium is μ0. transmitting antennas and receiving antennas are deployed with an equal space in the circular trajectories T and R.

The incident electric field generated by the th source is . Upon illumination by , the equivalent electric current density is induced on the target domain D, and generates scattered filed . , here note that the contrast function of scatterer is and the total electric field is for [6]. is the equation of field state, for , given as

is the data equation. For the scattered field measured by the receiving antennas is given by

Here, represents the 2D Green function of target domain D. The wavenumber is , where the angular frequency is denoted by ω and the wavelength is represented by . is the second kind Hankel function of zero order. 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 can be obtained by solving where and , .

2.2. Inexact Newton Formulation

Equation (3) describes the ill-posedness and nonlinearity of electromagnetic inverse imaging, and the optimal solution of 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 , where the inexact Newton method is given by

Here, is the number of outer inexact Newton iterations. represents the first-order multivariate Fréchet derivative of evaluated at . The solution is updated by , where is obtained by solving the inner ill-conditioned linear equation.

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 , according to (1), (2), and (3), let , , and represent the derivatives of and with respect to and , respectively. These operators are expressed as follows:

Remarkably, only depends on and , and only depends on . Thus, , , and for . is a sparse matrix with nonzero elements , , .

To discretize (5), (6), and (7), the target domain D is divided into square grid cells. The centres of the square cells and the locations of the receivers are expressed by and , respectively. To promote clarity, The centres of the square cells are represented by . The and nonzero elements of are where represents the diagonalization operation, and is a unit matrix. Equations (8), (9), (10), (11), and (12) describe the large-scale sparse structure of the Fréchet derivative matrix . 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 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

LSQR is an iterative regularization method, and the number of iterations is the regularization parameter. Implementation of this algorithm is carried out in two steps.

First, Lanczos bidiagonalization of , namely, , obtains the orthogonal matrix , and bidiagonal matrix . The form of is shown as

Lanczos iteration begins at the orthogonal basis vectors and the positive constant .

The iteration of this method is given as where , , . Equation (16) can also be expressed as where is the row of the n-dimensional identity matrix. Denoting the solution we need as , the residual vector in iteration is . The residual vector evaluates at , 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 . Since is a bidiagonal matrix, it can be easily solved by the QR decomposition method, and then one can obtain according to . The proposed IN-LSQR algorithm is shown in Algorithm 1.

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)
Initialization: initial solution
Outer inexact Newton iterations for
Calculate on the basis of (8) and (9), and calculate on the basis of (10), (11), and (12) and stored in sparse compressed form.
1.1 Internal initialization , , , , , . is a nonnegative real number
 Internal iterations for
1.2 Lanczos bidiagonalization of on the basis of (8), gives orthogonal matrix , and bidiagonal matrix .
1.3 Calculate the intermediate variable of the QR decomposition process in (18),
                      
                      
                      
                      
                     
                     
                      
where and are the initial value of the next iteration.
1.4 Update and , according to
                     
                    
1.5 If , then and repeat steps 1.2–1.4. Otherwise, terminate the internal iteration, achieving , the solution to (13).
2 Update the solution to (3) using
3 When the change of in two consecutive iterations is less than the termination condition, terminate the outer inexact Newton iteration, else and return to step 1.

The internal LSQR iteration obtains the solution , and the external update is applied to get the solution to the nonlinear inverse scattering problem.

Since the number of grid cells holds that and , 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 . The LSQR algorithm mainly calculates matrix-vector multiplication, or . Thus, the computational complexity can be reduced to . 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 . Consequently, the computational efficiency is significantly improved by use of sparse characteristics in .

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

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 Esca was synthetically generated with 20 dB noise. The inner iterations for both IN-LSQR and IN-LW, the outer iteration . 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.

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 : where represents the actual contrast in the target domain, and is the reconstructed contrast result at every inexact Newton iteration . in the experiments with dielTM and twodielTM was generated with 20 dB, 10 dB, and 5 dB Gaussian noise. Here, the internal iteration for both IN-LSQR and IN-LW. Figure 4 shows 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 in the experiments of dielTM and twodielTM with 20 dB noise. The computer used had an Intel [email protected] GHZ CPU, 64 GB RAM and MATLAB 2014b. According to Figure 4 and Table 2, computed by IN-LSQR converges faster than the IN-LW, and 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.

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 with diameter being 31 mm, where one of the cylinders is embedded in a larger cylinder with 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 for both IN-LSQR and IN-LW.

Figure 5(a) shows the actual profile of FoamTwinDiel. Figure 5(b) shows plots of 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 computed by IN-LSQR and CSI for different noise levels. As can be seen, 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 , respectively. The profile result reconstructed by IN-LSQR, similar to CSI, is clearly more accurate than that of IN-LW.

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.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (Grant nos. 61561034, 61661028, and 41505015), the Natural Science Foundation of Jiangxi Province, China (Grant no. 2015BAB207001), and the Key Technology Research and Development Program of Jiangxi Province, China (Grant no. 20151BBE50090).