Hostname: page-component-745bb68f8f-l4dxg Total loading time: 0 Render date: 2025-02-11T15:18:03.603Z Has data issue: false hasContentIssue false

Alternating direction method of multiplier for solving electromagnetic inverse scattering problems

Published online by Cambridge University Press:  11 March 2020

Jian Liu
Affiliation:
School of Information Engineering, Nanchang University, 999 Xuefu Avenue, Honggutan New District, Nanchang, Jiangxi, China
Huilin Zhou*
Affiliation:
School of Information Engineering, Nanchang University, 999 Xuefu Avenue, Honggutan New District, Nanchang, Jiangxi, China
Liangbing Chen
Affiliation:
School of Information Engineering, Nanchang University, 999 Xuefu Avenue, Honggutan New District, Nanchang, Jiangxi, China
Tao Ouyang
Affiliation:
School of Information Engineering, Nanchang University, 999 Xuefu Avenue, Honggutan New District, Nanchang, Jiangxi, China
*
Author for correspondence: Huilin Zhou, E-mail: zhouhuilin@ncu.edu.cn
Rights & Permissions [Opens in a new window]

Abstract

In this paper, a novel alternating direction method of multiplier (ADMM) is proposed to solve the inverse scattering problems. The proposed method is suitable for a wide range of applications with electromagnetic detection. In order to solve the internal ill-posed problem of the integral equation and make the reconstructed images more closer to the ground truth, first, the augmented Lagrangian method is introduced to transform the complex constrained optimization problem into the extremum problem of unconstrained cost function. Therefore, two artificial regularization factors of the cost function are optimized. Then, this proposed method decomposes the unconstrained global problem in the inversion process into three linear sub-problem forms of contrast source function, contrast function, and dual variables. And the form of the updated algebra for each sub-problem is not complicated. By cross-iterating and updating contrast source function, contrast function, and dual variables, the global minimization of the cost function can be accurately found. Finally, the proposed method is compared with the existing well-known iterative method for solving the inverse scattering problem. Both the numerical and experimental tests verify the validity and accuracy of the proposed ADMM.

Type
Radar
Copyright
Copyright © Cambridge University Press and the European Microwave Association 2020

Introduction

Electromagnetic inverse scattering imaging technique hascontactless and non-destructive characteristics when measuring geometric shapes and permittivities of unknown objects. Therefore, it plays an important role in the fields of biomedical imaging, microwave imaging, remote sensing, non-destructive evaluation,security checks, and so on [Reference Song, Li, Yang, Xu and Abubakar1Reference Zhuge and Yarovoy5]. Specifically, scattered field data and other prior information tested outside the domain of interest (DoI) are used to reconstruct unknown objects in the DoI. Since electromagnetic inverse scattering problem (EISP) is inherently nonlinear and ill-posed, it is challenging to obtain a solution of the global optimization in an efficient and steady approach. It is partly ascribable to the amount of unknowns of the electric field integral equation. Various methods have been proposed to deal with these difficulties. Generally speaking, these algorithms are divided into two categories; namely, the linear algorithms and the nonlinear algorithms. In these leading edge algorithms, the first Born approximations and the Rytov approximations belong to the linear algorithms [Reference Poli, Oliveri and Massa6,Reference Barriĺĺre, Idier, Goussard and Laurin7]. Contrast source-type inversion (CSI) method [Reference van den Berg and Kleinman8], subspace optimization method (SOM) [Reference Chen9], and variational Born iterative method [Reference Liang, Qiu, Han, Zhu, Liu, Liu, Liu, Fang and Liu10] belong to the nonlinear algorithms.

These advanced algorithms can effectively solve most inverse problems, but there are still some subtle defects. In linear algorithms, multiple scattering effects are not considered. The scattering field outside the DoI is simply ignored, and the total field is usually replaced by the incident field. Once the EISP is linearized, researchers can effectively solve it by using the preferred reconstruction algorithms. However, these algorithms are usually sparse and iterative [Reference Sung and Dasari11,Reference Kim, Zhou, Mir, Babacan, Carney, Goddard and Popescu12], and are only applied to weak scatterers. Compared with the linear algorithms, the quality of the spatial distribution of the reconstructed permittivity obtained by the nonlinear method is excellent. The nonlinear algorithm considers multiple electromagnetic scattering effects in the EISP and can make more accurate reconstruction of scatterers by increasing computational complexity of the reconstruction. In addition, nonlinear optimization algorithms take a long time to reconstruct scatterers. Therefore, it is not suitable for instantaneous applications or very complex problems. Moreover, it still needs to effectively combine a priori knowledge, starting from the initial guess of the spatial distribution of the permittivity in the EISP.

Optimization problems in imaging are often nontrivial to solve. When solving the EISP, the nonlinearity and the ill-posedness of its integral equation will lead to the complexity of the solving process [Reference Abdullah and Louis13]. To this end, we propose an eligible alternating direction method of multiplier (ADMM) to solve the EISP. This method has a unique advantage. It is different from the classical optimizations and its iteration assures the convergence of the reconstructed images in an iterative process. Aiming at the application of the ADMM in the EISP, the ADMM is suitable for dealing with large-scale data due to its significant low computational complexity. In the framework of the ADMM, the complex cost function in the EISP is transformed into an unconstrained problem. Then, by employing the decomposition-coordination method, the unconstrained global problem is decomposed into three sub-problems that are easier to solve. The specific expression form of each sub-problem is iteratively updated by different solving methods. By implementing the cross-iterating method, the global minimization of the cost function can be accurately found. This method provides relevant theoretical guidance and practical tests for the EISP of large-scale data. In addition, itis aware that there also exists another work that employs the ADMM to the EISP [Reference Li, Wang, Ding, Liu, Xia and Cui14]. Particularly, Li et al.regarded the equation of the external scattering field as a constraint in the ADMM framework and selected the equation of the internal total field as the objective function. By contrast, in our work, the proposed method regards the equation of the internal total field as a constraint in the ADMM framework and choose the equation of the external scattering field as the objective function. As the external scattering field is often corrupted by large noise in [Reference Chen9,Reference Liu15Reference Su, Xu, Lu and Jin17], thus our formulation can be more suitable to describe this circumstance. Furthermore, in order to pursue better robust recovery, more regularizations such as L1-norm are considered in our work.

The remainder of this paper is organized as follows. In Section “Inverse scatteringproblem,” the general model and basic theory of the EISP are introduced. In Section “Method,” the formulas and procedure of the ADMM are presented to solve the EISP. In Section “Experimental validation,” numerical and experimental examples demonstrate the effectiveness and accuracy of the ADMM. And conclusions are drawn in Section “Conclusion.”

Notation: $\overline {\overline \chi }$ and $\overline \omega$ denote the matrix and vector, respectively. Then, A T denotes the transpose of a matrix A. | · |F denotes the Frobenius norm of a matrix. In addition, in the reconstructed images, the units of the horizontal and vertical coordinates are meters. Finally, ψ( · ) denotes contraction operator.

Inverse scattering problem

The configuration of Fig. 1 is the same as the configuration of [Reference Abubakar and van den Berg18]. A two-dimensional(2-D) transverse-magnetic (TM) case (${\hat z}$ is the longitudinal direction) is considered in an EISP. As shown in Fig. 1, in free-space background, nonmagnetic scatterers are located in a DoI. The sources and receivers are equally placed outside the DoI, and their relevant position vectors are recorded as ${{\bf r}}_ {{\bf j}}\; \lpar {j = 1\comma\; 2\comma\; \ldots \comma\; {N_j}} \rpar$ and ${{\bf r}}_ {{\bf s}}\; \lpar {s = 1\comma\; 2\comma\; \ldots \comma\; {N_s}} \rpar$, respectively. A total number of N j line sources illuminate harmonic electromagnetic waves through the unknown scatterers. Then, the scattered electric field generated by each incidence is measured by an array of N s antennas, so the size of the obtained total data of the scattered field is N jN s. And we discretize the DoI into the total number of M small subunits. The unknown scatterers are located in the DoI where the background medium is evenly distributed, so the relevant permittivity of the background is ɛ0 and the relevant permeability is μ0.

Fig. 1. Test system of a 2-D inverse scattering problem, where unknown scatterers are located in a DoI. The sources and receivers are evenly distributed on a circle, respectively.

The total field in the DoI can be expressed as:

(1)$$E^{tot}_j\lpar {\bf r}\rpar =E^{inc}_j\lpar {\bf r}\rpar +k_b^2\int_D {{{G}_D}\lpar {\bf r}\comma\; {\bf r}^{\prime}\rpar }\chi\lpar {\bf r}^{\prime}\rpar E^{tot}_j\lpar {\bf r}^{\prime}\rpar d{\bf r}^{\prime}\comma\; $$

where $E^{inc}_j\lpar {\bf r}\rpar$ represents the incident field emitted by the jth antenna. The scattered field outside the DoI satisfies the following equation:

(2)$$E^{sca}_j\left({{{\bf r}_{\bf s}}} \right)= k_b^2\int_D {{G}_S}\left({{{\bf r}_{\bf s}}\comma\; {\bf r}^{\prime}} \right)\chi\left({{\bf r}^{\prime}} \right)E^{tot}_j\left({{\bf r}^{\prime}} \right)d{\bf r}^{\prime}\comma\; $$

where $G\lpar {\bf r}\comma\; {\bf r}^{\prime}\rpar = \lpar {1}/{{4j}}\rpar H_0^{\lpar 2\rpar }\lpar {k_b}\left \vert {{\bf r} - {\bf r}^{\prime}}\right \vert \rpar$ is a 2-D scalar Green's function and $H_0^{\lpar 2\rpar }\lpar{\cdot}\rpar$ is the zeroth order of Hankel function of the second kind. ${k_b} = \omega \sqrt {{\varepsilon _0}{\mu _0}} = 2\pi /{\lambda _0}$ represents the wave number of the background medium. ω is the angular frequency of the incident electromagnetic wave, and λ0 represents the wavelength. The contrast χ(r′) of the object can be expressed as χ(r′) = (ɛ(r′) − ɛ0)/ɛ0 + i(σ(r′)/ωɛ0). In order to process data conveniently, equations (1) and (2) are written in the form of matrixes hereinafter.

Method

The ADMM is a simple and effective algorithm for solving the distributed concave optimization problem. Therefore, it is used to solve imaging reconstruction for the sparse model [Reference Ji and Ye19]. The proposed algorithm mainly solves two types of problems. One is the constraint problem, and the other is the optimization problem. The ADMM transforms the global problem of a complex cost function into an unconstrained problem. The unconstrained problem is decomposed into several sub-problems by the decomposition-coordination method. Moreover, the analytical formulas of the transformed local sub-problems are relatively simple, so it is easy to find the iterative updating formula, and there is no need to constrain and converge for each local sub-problem. The ADMM can be used to effectively solve the L1-norm regularization problem [Reference Deng and Yin20]. Therefore, the ADMM has a large advantage in inversion speed and reconstruction accuracy because of its unique characteristics.

Therefore, the expression of the cost function in the EISP can be written as:

(3)$$\eqalign{L\left({\overline{\overline \chi } \comma\; {{\overline \omega }_j}\comma\; j = {{1\comma\; 2\comma\; }}\ldots} \right)& = {\eta _r}\sum\limits_{\,j = 1}^{{N_j}} {\left\Vert {\overline E^{sca}_j - {{\overline{ \overline G} }_S}{{\overline \omega }_j}} \right\Vert _2^2} \cr & \quad + {\eta _D}\sum\limits_{\,j = 1}^{{N_j}} {\left\Vert {{{\overline \omega }_j} - \overline{\overline \chi } \lpar \overline E^{inc}_j + {{\overline{\overline G} }_D}{{\overline \omega }_j}\rpar } \right\Vert _2^2}\comma\; }$$

where ${\eta _r} = {1}/{{\sum \nolimits _{j = 1}^{{N_j}} {\left \Vert {\overline E^{sca}_j} \right \Vert _2^2} }}$, and ${\eta _D} = {1}/{{\sum \nolimits _{j = 1}^{{N_j}} {\left \Vert {\overline E^{inc}_j} \right \Vert _2^2} }}$.

Equation (3) is transformed into a form with constrained optimization problems, which can be expressed as:

(4)$$\eqalign{& \arg \mathop {\min \limits_{{\overline{\overline \chi}}, {{\overline \omega }_j}}} \; {L_{uc}} \left({\overline{\overline \chi }}, {{\overline \omega }_j} E^{tot}_j, \; j = 1, \; 2, \; \ldots \right) \cr & {\rm s.t.} \; {\overline \omega_j} = \overline{\overline \chi } \left({\overline E^{inc}_j + {{\overline{\overline G} }_D}{{\overline \omega }_j}} \right), j = 1,\; 2,\; \ldots,} $$

where the expression of data fidelity item L uc in equation (4) is

(5)$${L_{uc}}\left({\overline{\overline \chi } \comma\; {{\overline \omega }_j}\comma\; \overline E^{tot}_j\comma\; j = 1\comma\; 2\comma\; \ldots}\right)={1\over 2}\sum\limits_{\,j = 1}^{{N_j}} {\left\Vert {\overline E^{sca}_j - {{\overline{\overline G} }_S}{{\overline \omega }_j}} \right\Vert _2^2}.$$

We extend the traditional real-valued augmented Lagrangian algorithm to the complex-valued domain. By its augmented Lagrangian form, equation (4) is transformed into a non-constrained problem, so that the augmented Lagrangian optimization problem at the complex-valued domain can be transformed into the following form:

(6)$$\mathop{\max}\limits_{\overline y}\mathop {\min }\limits_{{\overline{\overline \chi } \comma {{\overline \omega }_j}\comma \overline E^{tot}_j}}{L_{AL}}\left({\overline{\overline \chi } \comma\; {{\overline \omega }_j}\comma\; \overline E ^{tot}_j\comma\; j = 1\comma\; 2\comma\; \ldots} \right).$$

The augmented Lagrangian function L AL at the complex-valued domain in equation (6) is expressed as:

(7)$$\eqalign{{L_{AL}}\left({\overline{\overline \chi } \comma\; {{\overline \omega }_j}\comma\; \overline E^{tot}_j\comma\; \overline y} \right)& = {L_{uc}} + {\sum\limits_{\,j = 1}^{{N_j}} {\left\langle {\overline y \comma\; {{\overline \omega }_j} - \overline{\overline \chi } \left({\overline E^{inc}_j + {{\overline{\overline G} }_D}{{\overline \omega }_j}} \right)} \right\rangle } } \cr & \quad + {\rho \over 2}\sum\limits_{\,j = 1}^{{N_j}} {\left\Vert {{{\overline \omega }_j} - \overline{\overline \chi } \left({\overline E^{inc}_j + {{\overline{\overline G} }_D}{{\overline \omega }_j}} \right)} \right\Vert _2^2} {{ }}.}$$

Due to the nonlinearity and the ill-posedness of the electromagnetic integral equation, the stability and the uniqueness of the unknown contrast to be reconstructed cannot be guaranteed. In this case, in the process of the reconstruction of the scatterer, the range of the solution is expanded to some extent, so computational cost increases. The expression of equation (7) is combined with the L1-norm regularization, so the concrete expression of the cost function of the ADMM can be written as:

(8)$$\eqalign{{L_\rho }\left({\overline{\overline \chi } \comma\; {{\overline \omega }_j}\comma\; \overline E^{tot}_j\comma\; {{\overline y }}} \right)& = {1\over 2}\sum\limits_{\,j = 1}^{{N_j}} {\left\Vert {\overline E^{sca}_j - {{\overline{\overline G} }_S}{{\overline \omega }_j}} \right\Vert _2^2}\cr & \quad + \lambda {\left\Vert {\overline{\overline \chi } } \right\Vert _1}{{ + }} {\sum\limits_{\,j = 1}^{{N_j}} {\left\langle {\overline y \comma\; {{\overline \omega }_j} - \overline{\overline \chi } \left({\overline E^{inc}_j + {{\overline G }_D}{{\overline \omega }_j}} \right)} \right\rangle } } \cr & \quad + {\rho \over 2}\sum\limits_{\,j = 1}^{{N_j}} {\left\Vert {{{\overline \omega }_j} - \overline{\overline \chi } \left({\overline E^{inc}_j + {{\overline{\overline G} }_D}{{\overline \omega }_j}} \right)} \right\Vert _2^2} \comma\; }$$

where the second term on the right side ofequation (8) is the L1-norm regularization term, and λ represents the regularization coefficient. ${\overline y}$ denotes the dual value of the complex-valued Lagrangian (Lagrangian multiplier) for the balance of the data fidelity item and the second penalty. ρ is a parameter of the penalty term, and ρ > 0.

When using the ADMM framework to optimize the EISP, the solving process mainly includes three layers of variables. The observed variable of the first layer needs to obtain the scattering field data $\overline E_j^{sca}$ outside the DoI. The variable of the second layer is the contrast source ${\overline \omega }_j$ that is introduced to simplify the number of unknown variables. The variable of the third layer is the contrast $\overline {\overline \chi }$ to be reconstructed. n represents the number of iterations. When processing in a single layer, the variables are independent of each other.

Cost function of equation (8) is decomposed into the following three iterative equations:

(9)$${\overline \omega _{n + 1}} =\arg\mathop {\min }\limits_{\overline \omega}L\rho \left({{\overline \omega}\comma\; {{\overline {\overline\chi} }_n}\comma\; {{\overline y }_n}} \right)\comma\; $$
(10)$${\overline {\overline \chi } _{n + 1}} =\arg\mathop {\min }\limits_{\overline {\overline\chi }}L\rho \left({{{\overline \omega_{n+1} }}}\comma\; {\overline{\overline \chi } }\comma\; {{\overline y }_n} \right)\comma\; $$
(11)$${\overline y _{n + 1}} =\arg\mathop {\min }\limits_{\overline y}L\rho \left({{{\overline \omega_{n+1} }}}\comma\; {{\overline{ \overline \chi }_{n + 1}}}\comma\; {\overline y } \right).$$

According to equations (9), (10), and (11), the following three sub-problems are handled separately:

(12)$$\eqalign{\overline \omega _{n + 1}^{\;j} & = \arg\mathop {\min }\limits_{\overline \omega} {1\over 2}\sum\limits_{\,j = 1}^{{N_j}} {\left\Vert {\overline E^{sca}_j - {{\overline{\overline G} }_{S}}\overline \omega _n^{\;j}} \right\Vert _2^2} \cr & \quad +\sum\limits_{\,j = 1}^{{N_j}}{\left\langle {\overline y _n^{\;j}\comma\; \overline \omega _n^{\;j}-\overline{\overline \chi }\overline E _{\,j\comma n}^{tot} } \right\rangle } + {\rho \over 2}\sum\limits_{\,j = 1}^{{N_j}} {\left\Vert {\overline \omega _n^{\;j} - \overline{\overline \chi }\overline E_{\,j\comma n}^{tot} } \right\Vert _2^2} \comma\; }$$
(13)$$\eqalign{{\overline{\overline \chi } _{n + 1}} & = \arg\mathop {\min }\limits_{\overline {\overline\chi }} {1\over 2}\sum\limits_{\,j = 1}^{{N_j}} {\left\Vert {\overline \omega _{n + 1}^{\;j} - {{\overline{\overline \chi } }_n}\overline E_{\,j\comma n+1}^{tot} } \right\Vert _2^2} \cr & \quad + \sum\limits_{\,j = 1}^{{N_j}} {\left\langle {{{\overline y }_n}\comma\; \overline \omega _{n + 1}^{\;j} - {{\overline{\overline \chi } }_n}\overline E_{\,j\comma n+1}^{tot} } \right\rangle } {\rm{ + }}\lambda {\left\Vert {\overline{\overline \chi } } \right\Vert _1}\comma\; }$$
(14)$${\overline y _{n + 1}}=\arg\mathop {\min }\limits_{\overline y}\sum\limits_{\,j = 1}^{{N_j}}{\left\langle {{{\overline y }_n}\comma\; \overline \omega _{n + 1}^{\;j}-{{\overline{\overline \chi } }_{n + 1}}\overline E_{\,j\comma n+1}^{tot} } \right\rangle }\comma\; $$

where $\overline E_{j\comma n}^{tot}= \overline E^{inc}_j + {\overline {\overline G }_D}\overline \omega _n^{\;j}$.

For the sake of simplicity, the above three sub-problems are respectively simplified. The scaled dual variable is introduced, and $\overline u = {{\overline {{y}} }}/{\rho }$ is defined, thereby obtaining:

(15)$$\eqalign{\overline \omega _{n + 1}^{\;j} & = \arg\mathop {\min }\limits_{\overline \omega} {1\over 2}\sum\limits_{\,j = 1}^{{N_j}} {\left\Vert {\overline E^{sca}_j - {{\overline G }_S}\overline \omega _n^{\;j}} \right\Vert _2^2} \cr & \quad + {\rho \over 2}\sum\limits_{\,j = 1}^{{N_j}} {\left\Vert {{{\overline{\overline \chi } }_n}\overline E_{\,j\comma n}^{tot} - \overline \omega _n^{\;j} + {{\overline u }_n}} \right\Vert _2^2}\comma\; }$$
(16)$$\eqalign{\overline{\overline \chi } _{n + 1}& = \arg\mathop {\min }\limits_{\overline {\overline\chi }} {\rho \over 2}\sum\limits_{\,j = 1}^{{N_j}} {\left\Vert {{{\overline{\overline \chi } }_n}\overline E_{\,j\comma n+1}^{tot} - \overline \omega _{n + 1}^{\;j} + {{\overline u }_n}} \right\Vert _2^2}\cr & \quad + \lambda {\left\Vert {{{\overline{\overline \chi } }_n}} \right\Vert _1}\comma\; }$$
(17)$${\overline u _{n + 1}} = {\overline u _n} + \left({{{\overline{\overline \chi } }_{n + 1}}\overline E_{\,j\comma n+1}^{tot} - \overline \omega _{n + 1}^{\;j}} \right).$$

According to equations (15), (16), and (17), the contrast source function $\overline \omega$ is solved by the least squares method. Then, the gradient of $\overline \omega$ is solved via equation (15) and the gradient is zero. The collation can be obtained, and the updated equation of the contrast source function $\overline \omega ^{\;j}$ is

(18)$$\overline \omega _{n + 1}^{\;j} =PQ\comma\; $$

where $P={\lpar {\overline {\overline G} _S^T{{\overline {\overline G} }_S} + \rho {{\lpar {\overline {\overline I} - \overline {\overline \chi } {{\overline {\overline G} }_D}} \rpar }^T}\lpar {\overline {\overline I} - \overline {\overline \chi } {{\overline {\overline G} }_D}} \rpar } \rpar ^{-1}}$, and $Q=( \overline {\overline G} _S^T\overline E^{sca}_j + \rho {{\left( {\overline{\overline I} - \overline{\overline \chi } {{\overline{\overline G} }_D}} \right)}^T}\left( {\overline{\overline \chi } \overline E^{inc}_j - {{\overline u }_n}} \right) )$. $\overline {\overline I}$ is the identity matrix. According to the form of the shrinkage threshold, the update of the contrast function $\overline {\overline \chi }$ in equation (16) can be solved. According to the essence of the proposed algorithm, in order to simplify equation, the above equation (16) can be recorded as:

(19)$$\eqalign{F\left({\overline{\overline \chi } } \right)& = \arg\mathop {\min }\limits_{\overline {\overline\chi }}{\rho \over 2}\left\Vert {\overline{\overline \chi }\overline E_{\,j\comma n+1}^{tot}\! -\! \left({\overline \omega _{n + 1}^{\;j}\! -\! {{\overline u }_n}} \right)} \right\Vert _2^2 \quad +\; \lambda {\left\Vert {\overline{\overline \chi } } \right\Vert _1}.}$$

For the solution of equation (19), it is generally solved by the iterative soft thresholding algorithm [Reference Wang, Aboutanios and Amin21]. Using this algorithm, the updating algebra of the contrast function $\overline {\overline \chi }$ in the iterative process of inversion is written as:

(20)$${\overline{\overline \chi } _{n + 1}} = \rho {\psi _\lambda }\left({{{\overline{\overline \chi } }_n} + \overline E_{\,j\comma n+1}^{tot}\left({{{\overline{\overline \chi } }_n}\overline E_{\,j\comma n+1}^{tot} - \overline b } \right)} \right)\comma\; $$

where ${\psi _\lambda }\lpar {{{\overline {\overline \chi } }_i}} \rpar = {{\rm sgn}} \lpar {{{\overline {\overline \chi } }_i}} \rpar {\lpar {\left \vert {{{\overline {\overline \chi } }_i}} \right \vert - \lambda } \rpar _{\rm {\ +\ }}}$. And $\overline b = \overline \omega _{n + 1}^{\;j} - {\overline u _n}$. In each updating iteration, once the solution is obtained, the soft threshold needs to be shrunk.

The gradient descent method is used to solve the scaled dual variable. Therefore, the iterative process of the ADMM in the EISP can be obtained by combining equations (17), (18), and (20).

The procedures of the detailed iteration are stated as follows.

  • Step 1: Set n = 0. Calculate the relevant parameters, and let $\overline E^{tot}_j = \overline E^{inc}_j$ at initialization. Calculate $\overline {\overline \omega }_0$ by the backpropagation (BP) method [Reference van den Berg and Kleinman8]. Then, $\overline {\overline \chi }_0$ can be obtained.

  • Step 2: n = n + 1. $\overline \omega _{n - 1}^{\;j}$ and ${\overline {\overline \chi } _{n - 1}}$ are known at the next iteration. Therefore, the data item error is ${\overline r _n}\;{{ = }}\;\overline E^{sca}_j - {\overline {\overline G} _S}\overline \omega _{n - 1}^{\;j}$, and the status item error is ${\overline s _n} = \overline \omega _{n - 1}^{\;j} - {\overline {\overline \chi } _{n - 1}}\lpar {\overline E^{inc}_j + {{\overline {\overline G} }_D}\overline \omega _{n - 1}^{\;j}} \rpar$.

  • Step 3: Use equation (18) to calculate the contrast source function ${\overline \omega ^{\;j}}$. Then, update the total field $\overline E^{tot}_j = \overline E^{inc}_j + {\overline {\overline G} _D}{\overline \omega ^{\;j}}$.

  • Step 4: Similar to the update of the contrast source function ${\overline \omega ^{\;j}}$, the contrast function $\overline {\overline \chi }$ is updated according to equation (20).

  • Step 5: Update the Lagrange multiplier (dual variable) $\overline y$ using equation (17).

  • Step 6: If the cost function ${L_\rho }$ in equation (8) at this time meets the termination condition, the iteration is terminated. Otherwise, return to step 2 again for the next iteration.

As can be seen from the basic iterative process of the ADMM, the ADMM has the characteristics of distributed optimization. Compared with the CSI and the SOM, two artificial regularization factors are avoided. The ADMM is suitable for solving large-scale sparse matrices, and has the characteristics of fast convergence and high quality of inversion reconstruction.

Experimental validation

This section considers some reconstructed results constructed using numerical and experimental data at a single frequency to evaluate the performance of the proposed ADMM. In these tests, we use handwritten words in EMNIST database [Reference Cohen, Afshar, Tapson and van Schaik22], which includes some commonly numbers, English letters, and Greek letters. However, due to the irregularities and the singular characteristics of these handwritten words, it is very difficult to restore the original shape and the permittivity. In addition, we compare the ADMM with the well-known CSI and the SOM in this field. Finally, the ADMM is applied to the experimental data to verify its high efficiency.

Numerical results

In the forward model, the handwritten words in EMNIST database were used as scatterers for 2-D modeling. These handwritten images are assigned to a DoI with the size of 2 × 2 m2. What's more, the DoI is divided into 112 × 112 pixels. In the inversion process, the DoI is divided into 64 × 64 pixels in order to avoid the inverse crisis. In total, 16 line sources and 32 line receivers are equally placed on a circle of radius of 3 m, wherein the center of the circle is at the origin. In addition, the operating frequency of is 400 MHz. These scatterers have a permittivity of 2. The background permittivity is 1. Therefore, method of moment is used to obtain the corresponding scattering field matrix $\overline {\overline K}$ whose size is N s × N j. On this basis, additive white Gaussian noise matrix $\overline {\overline n}$ is added to $\overline {\overline K}$. Therefore, $\overline {\overline K} + \overline {\overline n}$ is the final scattering field matrix as a known condition to reconstruct scatterers for the EISP. The noise level is described as $\lpar {\left \Vert {\overline {\overline n} } \right \Vert _F}/{\left \Vert {\overline {\overline K} } \right \Vert _F}\rpar \times 100\percnt$. A priori information is that the scatterers are lossless and have nonnegative contrast [Reference van den Berg and Kleinman8]. Finally, the spatial distribution of the permittivity of the scatterer is reconstructed by using the CSI, the SOM, and the ADMM, respectively.

In order to evaluate the quality of these algorithms, the mean square error of the reconstructed permittivity is defined as:

(21)$$MSE = \sqrt {{1\over N}\sum\limits_{{m}}^M {{{{{\left\vert {\varepsilon _m^{inv} - \varepsilon _m^{gt}} \right\vert }^2}}\over {{{\left\vert {\varepsilon _m^{gt}} \right\vert }^2}}}} }\comma\; $$

where M is the number of grids in the whole DoI. $\varepsilon _m^{inv}$ denotes the permittivity of each iteration and $\varepsilon _m^{gt}$ represents the ground truth of scatterers in the DoI.

In the first example, there are four representative profiles in Fig. 2 to be reconstructed, including the English letter G, the number 8, the Greek letter ɛ and the famous “Austria ring” profile composed of two ear-shaped disks and one ring. In particular, we often use “Austria ring” as a benchmark to evaluate the performance of electromagnetic inverse scattering algorithms. Obviously, these irregular handwritten words are more challenging than regular shapes. Three different methods are used for each profile to reconstruct the spatial distribution of its permittivity with noise-free scattered field data. The permittivities of reconstructed images in the ADMM are distributed in the range of 1.7–2.3, and the shape also has less difference from the ground truth. Compared with the CSI and the SOM, the ADMM has achieved satisfactory results. Some scatterers in the CSI and the SOM are somewhat distorted and inaccurate, while the ADMM does not. Especially, the lack of information appears at some corners of the image in (b) and (c) of Fig. 2. Therefore, it can be observed that the proposed ADMM has obvious advantages over the CSI and the SOM. These numerical results show that the proposed ADMM can effectively solve the inverse scattering problem.

Fig. 2. Tests on “Austria ring” profile and some profiles in EMNIST dataset with noise-free scattered field data. Here, the ground truth profiles in EMNIST dataset are displayed in (a), (e), and (i), respectively. The ground truth of “Austria ring”in (m). Reconstructed images (b), (f), (j), and (n) by the CSI; (c), (g), (k), and (o) by the SOM; (d), (h), (l), and (p) by the ADMM, respectively.

In the second example, we test the more challenging “Austria ring” profile with different signal to noise ratios (SNRs) in Fig. 3. Under the condition of 15% Gaussian noises (SNR = 12 dB), the spatial distribution of the permittivity of the scatterer can be well reconstructed in the CSI, the SOM, and the ADMM, respectively. However, as the noise level increases, the reconstructed images of three algorithms show distortion of shape and inaccuracy of permittivity, and the original profiles cannot be completely restored. At 50% Gaussian noises (SNR = 6 dB), it can be seen from (g) and (h) that the CSI and the SOM have great defects for the shape reconstruction of two disks of “Austria ring.” Although the defects are also obvious, the position of disks and ring is correctly reconstructed (i). It can be seen that the noise immunity of the ADMM is better than that of the CSI and the SOM at low SNRs.

Fig. 3. Testing the noise immunity of the CSI, the SOM, and the ADMM in “Austria ring”. Top images, middle images, and bottom images are reconstructed from the scattered field data with 25% (SNR = 12 dB), 35% (SNR = 9 dB), and 50% (SNR = 6 dB) Gaussian noises, respectively. Left images, middle images, and right images are reconstructed by the CSI, the SOM, and the ADMM, respectively.

As presented in Fig. 4, the convergent speed of MSE curve is still fast under the noisy condition. Due to the use of the BP algorithm, the difference of initial MSE is very small. And under 50 iterations, all three algorithms can meet the convergent condition. At SNR = 12 dB, the convergent speed of the ADMM is much faster than the convergent speed of the SOM and the CSI, so the imaging speed is faster. At different SNRs, the convergent speed of the ADMM is still very fast. At low SNRs, the convergent speed of the ADMM is reduced. Therefore, the ADMM converges faster than other algorithms.

Fig. 4. Comparison of convergent trajectories in the first 50 iterations for different algorithms and different SNRs.

In the third example, in order to illustrate the insensitivity of the present method to large noise corruption, reconstruction comparison between IS-ADMM [Reference Li, Wang, Ding, Liu, Xia and Cui14] and the proposed ADMM is conducted. In the experiment, the scatterer in Fig. 5(a) contains a square with a side length of 0.6 m and a circle with a radius of 0.3 m. Their relative permittivities are 2. At the same time, three noises with different SNRs are introduced into the scatterer, i.e. no noise, SNR = 12 dB, and SNR = 9 dB. The constructions from the three observed measurements are depicted in Fig. 5. Generally, it can be seen that the scatterer can be roughly reconstructed using two different methods. Obviously, the reconstructions with the proposed ADMM are clearer under the same SNR, indicating that the proposed ADMM exhibits better robustness. Since the proposed ADMM additionally adds the L1-norm as a constraint to the solution in the cost function, the reconstruction with the ADMM performs more smoothly at the edges and has a more uniform distribution.

Fig. 5. Comparisons of the reconstructions by IS-ADMM in [Reference Li, Wang, Ding, Liu, Xia and Cui14], and the proposed ADMM under different noises. (a) Ground truth. (b) No noise with IS-ADMM. (c) SNR = 12 dB with IS-ADMM. (d) SNR = 9 dB with IS-ADMM. (e) No noise with ADMM. (f) SNR = 12 dB with ADMM. (g) SNR = 9 dB with ADMM.

Experimental data

In order to better verify the proposed ADMM, experimental data from Institute Fresnel [Reference Geffrin, Sabouroux and Eyraud23] is used to reconstruct the image. A “FoamTwinDiel” profile with TM case is applied in the ADMM. There are 241 line receivers and 18 line sources. As can be seen from Fig. 6(a), there are two same plastic cylinders (berylon) with a diameter of 31 mm. Its permittivity is ɛ = 3 ± 0.3. A foam cylinder (SAITEC SBF 300) has a diameter of 80 mm with the permittivity ɛ = 1.45 ± 0.15. One plastic cylinder is embedded in the foam cylinder, and the other plastic cylinder is tangent to the foam cylinder. The ± sign indicates that the permittivity of the experimental measurement is inaccurate. Because the operating frequency is changed to 5 GHz, the DoIis alsochanged to 0.15 × 0.15 m2. At this moment, MSE becomes the mean square error of the permittivity of the current iteration and the previous iteration.

Fig. 6. “FoamTwinDiel” profiles. (a) The original shape. Reconstructed images of “FoamTwinDiel” profile in (b) the CSI, (c) the SOM, and (d) the ADMM, respectively.

As presented in Fig. 6, in the reconstructed scatterers, the reconstructions of the foam cylinder and the external plastic cylinder have achieved satisfactory results in the CSI, the SOM, and the ADMM. Unfortunately, images in the CSI and the SOM have obvious distortions. The plastic cylinder embedded in the foam cylinder apparently fails to reconstruct while the ADMM can do it successfully. The spatial distribution of the permittivity of the reconstructed “FoamTwinDiel” profile in the ADMM is given in Fig. 6(d). It can be seen that the permittivities of plastic cylinders are in the range of 2.5–3.5 and the permittivities of the foam cylinder are in the range of 1.2–1.7 with the ADMM in Fig. 6(d). The spatial position is basically accurate. The ADMM is significantly better than the CSI and the SOM in reconstruction. The experimental results show that the proposed ADMM can also show an excellent performance for the experimental data. This shows the high efficiency and the accurateness of the ADMM.

Maximum 50 iterations are set in the inversion process, and all three algorithms can reach convergent conditions. As shown in Fig. 7, the CSI, the SOM, and the ADMM can converge quickly. But, it is clear that MSE of the ADMM drops faster, which means that the reconstructed images are closer to the ground truth more quickly. Consequently, the ADMM is more efficient.

Fig. 7. Comparison of convergent trajectories in the first 50 iterations for different algorithms with “FoamTwinDiel” profile.

Conclusion

In this paper, we propose an ADMM to solve the EISP. During the inversion process, the augmented Lagrangian method is introduced to transform the complex nonlinear cost function into an unconstrained global problem. Two artificial regularization factors in the CSI and the SOM are optimized. By the decomposition-coordination method, this unconstrained problem is decomposed into three solved sub-problems easily. The global optimization process of the cost function is achieved by alternately updating the contrast source function, the contrast function, and dual variables. In addition, we found that ρ = 0.2 is the best choice according to the trial and error method. Moreover, we prove that the anti-noise performance of the ADMM is better than other algorithms in the numerical results. Not only is the reconstruction of the scatterer highly efficient in numerical experiments, but it can also effectively reconstruct complex profiles in practical applications.

Although the proposed method can outperform others in the performance of reconstructing scatterers, it still has some shortcomings. The selection of some parameters (such as ρ) has subjective dependence and non-generality. In future work, we will consider the deep unrolling network. The knowledge in the imaging field and a priori information are merged into the deep network to adaptively acquire some parameters in the ADMM.

Jian Liu was born in 1995 in Jiangxi, China. He received his bachelor's degree in 2018 at the Department of Electronic Information, School of Information Engineering, Nanchang University. He is currently pursuing his master's degree in electromagnetic field at Nanchang University. He is mainly engaged in the research of inverse scattering imaging methods and radar signal processing.

Huilin Zhou was born in Jiangxi, China, in 1979. He received his Ph.D. degree in space physics from Wuhan University, Wuhan, China, in 2006. His research interests include ultra-wideband radar imaging and radar signal processing. He is currently the director of the Department of Electronic Information Engineering of Nanchang University, a peer reviewer of the National Natural Science Foundation of China, and a member of the China Electronics Society. He is an anonymous reviewer of international SCI journals such as Advanced Robot Systems, Electromagnetic Waves and Applications, and Applied Computing Electromagnetics.

Liangbing Chen received his B.E. and Ph.D. degrees in electrical engineering from the Huazhong University of Science and Technology, Wuhan, China, in 2005 and 2010, respectively. Presently, he is a researcher in the Department of Electronic Information Engineering, Nanchang University, Nanchang, China. His research interests include microwave remote sensing, array signal processing, and passive millimeter-wave imaging.

Tao Ouyang is currently pursuing his master's degree in electromagnetic field at Nanchang University, Jiangxi, China. His research interests include inverse scattering and radar signal processing.

References

1.Song, X, Li, M, Yang, F, Xu, S and Abubakar, A (2019) Study on joint inversion algorithm of acoustic and electromagnetic data in biomedical imaging. IEEE Journal on Multiscale and Multiphysics Computational Techniques, 4, 211.CrossRefGoogle Scholar
2.Palmeri, R, Bevacqua, MT, Crocco, L, Isernia, T and Donato, LD (2017) Microwave imaging via distorted iterated virtual experiments. IEEE Transactions on Antennas and Propagation, 65, 829838.CrossRefGoogle Scholar
3.Chen, G, Stang, J, Haynes, M, Leuthardt, E and Moghaddam, M (2018) Real-time three-dimensional microwave monitoring of interstitial thermal therapy. IEEE Transactions on Biomedical Engineering, 65, 528538.CrossRefGoogle ScholarPubMed
4.Ireland, D, Bialkowski, K and Abbosh, A (2013) Microwave imaging for brain stroke detection using Born iterative method. IET Microwaves, Antennas & Propagation, 7, 909915.CrossRefGoogle Scholar
5.Zhuge, X and Yarovoy, AG (2011) A sparse aperture MIMO-SAR-based UWB imaging system for concealed weapon detection. IEEE Transactions on Geoscience and Remote Sensing, 49, 509518.CrossRefGoogle Scholar
6.Poli, L, Oliveri, G and Massa, A (2012) Microwave imaging within the first-order Born approximation by means of the contrast-field Bayesian compressive sensing. IEEE Transactions on Antennas and Propagation, 60, 28652879.CrossRefGoogle Scholar
7.Barriĺĺre, P, Idier, J, Goussard, Y and Laurin, J (2010) Fast solutions of the 2D inverse scattering problem based on a TSVD approximation of the internal field for the forward model. IEEE Transactions on Antennas and Propagation, 58, 40154024.CrossRefGoogle Scholar
8.van den Berg, PM and Kleinman, RE (1997) A contrast source inversion method. Inverse Problems, 13, 16071620.CrossRefGoogle Scholar
9.Chen, X (2010) Subspace-based optimization method for solving inverse-scattering problems. IEEE Transactions on Geoscience and Remote Sensing, 48, 4249.CrossRefGoogle Scholar
10.Liang, B, Qiu, C, Han, F, Zhu, C, Liu, N, Liu, H, Liu, F, Fang, G and Liu, Q (2018) A new inversion method based on distorted Born iterative method for grounded electrical source airborne transient electromagnetics. IEEE Transactions on Geoscience and Remote Sensing, 56, 877887.CrossRefGoogle Scholar
11.Sung, Y and Dasari, RR (2011) Deterministic regularization of three-dimensional optical diffraction tomography. Journal of the Optical Society of America, 28, 15541561.CrossRefGoogle ScholarPubMed
12.Kim, T, Zhou, R, Mir, M, Babacan, SD, Carney, PS, Goddard, LL and Popescu, G (2014) White-light diffraction tomography of unlabelled live cells. Nature Photonics, 8, 256263.CrossRefGoogle Scholar
13.Abdullah, H and Louis, AK (1999) The approximate inverse for solving an inverse scattering problem for acoustic waves in an inhomogeneous medium. Inverse Problems, 15, 12131229.CrossRefGoogle Scholar
14.Li, L, Wang, LG, Ding, J, Liu, PK, Xia, MY and Cui, TJ (2017) A probabilistic model for the nonlinear electromagnetic inverse scattering: TM case. IEEE Transactions on Antennas and Propagation, 65, 59845991.CrossRefGoogle Scholar
15.Liu, Z (2019) A new scheme based on Born iterative method for solving inverse scattering problems with noise disturbance. IEEE Geoscience and Remote Sensing Letters, 16, 10211025.CrossRefGoogle Scholar
16.Liu, Z and Nie, Z (2019) Subspace-based variational Born iterative method for solving inverse scattering problems. IEEE Geoscience and Remote Sensing Letters, 16, 10171020.CrossRefGoogle Scholar
17.Su, H, Xu, F, Lu, S and Jin, Y (2016) Iterative ADMM for inverse FE-BI problem: a potential solution to radio tomography of asteroids. IEEE Geoscience and Remote Sensing Letters, 54, 52265238.CrossRefGoogle Scholar
18.Abubakar, A and van den Berg, PM (2002) The contrast source inversion method for location and shape reconstructions. Inverse Problems, 18, 495510.CrossRefGoogle Scholar
19.Ji, S and Ye, J (2009) An accelerated gradient method for trace norm minimization. The 26th Annual International Conference on Machine Learning, Canada, June 2009. Montreal, Canada: International Machine Learning Society (IMLS).CrossRefGoogle Scholar
20.Deng, W and Yin, W (2016) On the global and linear convergence of the generalized alternating direction method of multipliers. Journal of Scientific Computing, 66, 889916.CrossRefGoogle Scholar
21.Wang, X, Aboutanios, E and Amin, MG (2014) Thinned array beampattern synthesis by iterative soft-thresholding-based optimization algorithms. IEEE Transactions on Antennas and Propagation, 62, 61026113.CrossRefGoogle Scholar
22.Cohen, G, Afshar, S, Tapson, J and van Schaik, A (2017) EMNIST: An extension of MNIST to handwritten letters. Arxiv:1702.05373.Google Scholar
23.Geffrin, JM, Sabouroux, P and Eyraud, C (2005) Free space experimental scattering database continuation: experimental set-up and measurement precision. Inverse Problems, 21, S117S130.CrossRefGoogle Scholar
Figure 0

Fig. 1. Test system of a 2-D inverse scattering problem, where unknown scatterers are located in a DoI. The sources and receivers are evenly distributed on a circle, respectively.

Figure 1

Fig. 2. Tests on “Austria ring” profile and some profiles in EMNIST dataset with noise-free scattered field data. Here, the ground truth profiles in EMNIST dataset are displayed in (a), (e), and (i), respectively. The ground truth of “Austria ring”in (m). Reconstructed images (b), (f), (j), and (n) by the CSI; (c), (g), (k), and (o) by the SOM; (d), (h), (l), and (p) by the ADMM, respectively.

Figure 2

Fig. 3. Testing the noise immunity of the CSI, the SOM, and the ADMM in “Austria ring”. Top images, middle images, and bottom images are reconstructed from the scattered field data with 25% (SNR = 12 dB), 35% (SNR = 9 dB), and 50% (SNR = 6 dB) Gaussian noises, respectively. Left images, middle images, and right images are reconstructed by the CSI, the SOM, and the ADMM, respectively.

Figure 3

Fig. 4. Comparison of convergent trajectories in the first 50 iterations for different algorithms and different SNRs.

Figure 4

Fig. 5. Comparisons of the reconstructions by IS-ADMM in [14], and the proposed ADMM under different noises. (a) Ground truth. (b) No noise with IS-ADMM. (c) SNR = 12 dB with IS-ADMM. (d) SNR = 9 dB with IS-ADMM. (e) No noise with ADMM. (f) SNR = 12 dB with ADMM. (g) SNR = 9 dB with ADMM.

Figure 5

Fig. 6. “FoamTwinDiel” profiles. (a) The original shape. Reconstructed images of “FoamTwinDiel” profile in (b) the CSI, (c) the SOM, and (d) the ADMM, respectively.

Figure 6

Fig. 7. Comparison of convergent trajectories in the first 50 iterations for different algorithms with “FoamTwinDiel” profile.