Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

ICSI: an improved contrast source inversion method for electromagnetic inverse scattering problems

Open Access Open Access

Abstract

An improved contrast source inversion method (ICSI) is proposed to solve the problem of electromagnetic inverse scattering. Specifically, we first resort Fourier bases to represent contrast source, so that we can solve contrast source in the frequency domain, which can accelerate the solution of contrast source. Then, a spatial-frequency domain constraint is designed to ensure that the spanned-domain calculation is synchronous optimal. Afterwards, a multi-round optimization combined with SVD (the singular value decomposition) is developed to alleviate the reliance on initial guesses. Finally, we utilize a frequency-domain filter to eliminate redundant inversion information and narrow the search scope of the solution. Extensive experiments on the synthetic and real data show that ICSI holds a faster computing speed, a better inversion ability of relative permittivities, and a stronger anti-noise ability.

© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

1. Introduction

The electromagnetic inverse scattering problem (EISP) is to research how to estimate an unknown scatter’s attributes, such as geometric attributes (position, shape, size) and physical attributes (refractive index, relative permittivities), according to known incident fields and measured scattered fields. It is widely used in many fields, including medical imaging diagnosis [1], buried target detection [2], remote sensing testing [3], nondestructive testing [4], detection in civil structures [5], etc. The existing solutions to the problem can be divided into physical methods and learning methods. The physical methods mainly construct a class of Lippmann Schwinger equations based on Maxwell’s equations of electromagnetic theory. The physical attributes of the scatterer to be estimated can be obtained by solving the Lippmann Schwinger equation. However, it is challenging in the solving process due to the existence of the multiple-scattering, the limited measured data, and the non-linearity of the equation. To address these challenges, several common methods have been developed, such as contrast source inversion (CSI) [6], subspace optimization method (SOM) [7], contraction integral for inversion (CIE-I) [8], etc. Nevertheless, most of these methods need gradient calculation and iterative update, resulting in low efficiency. To alleviate this problem, the early data-based learning methods [9] are devoted to training artificial neural networks (ANN) as a black box solver. They are simple and have the advantage of real-time in the testing phase, but their inversion performance is affected because of a lack of physical knowledge. Therefore, the model-driven deep learning methods [10] resort physical model to guide the training of ANN, which achieves promising inversion results. As a result, the physical methods still play an irreplaceable role. In this paper, the CSI method, which is one of the common physical methods for reconstructing a scatter’s relative permittivities, is studied and improved.

CSI was first proposed in 1997 [11]. This method only minimizes the numerator of the objective function and updates the contrast based on it, so that the objective function may not be reduced. In 1999, an improved version [6] not only takes into account the contrast of denominator but also adopts the alternative direction method of multiplier (ADMM) [12] and conjugate gradient (CG) to update the contrast source and contrast alternately. However, practical applications often require faster calculation, clearer contours, better inversion, stronger anti-noise ability, and smaller errors. These issues have motivated researchers to constantly explore and refine CSI. [13] solves the contrast source coefficients in the wavelet domain with good sparsity, which is helpful to reduce the number of unknowns and the degree of nonlinearity. [6] presents MR-CSI which uses multiplicative regularization to improve the contours of relative permittivities. In [14], cross-correlated contrast source inversion (CC-CSI) introduces a mismatch measure of data residuals and state residuals to obtain a better inversion result and robustness. CSI and its related methods usually calculate the contrast source in the spatial domain, which has a heavy computational burden [15]. At the same time, redundant information (e.g., measurement noise or multiple-scattering interference [16]) also participates in the reconstruction, which will limit the reconstruction range of relative permittivities and computational efficiency. In addition, most of these methods are optimized by gradient iteration, which depends on initial guess. Although BP method gives an improved initial estimate, it may also lead to poor inversion or even failure.

To overcome the above shortcomings, we propose an improved contrast source inversion method that we call it ICSI. To be specific, inspired by [15,17], we first calculate the Fourier coefficients of contrast source in the frequency domain to accelerate calculation. Then, a spatial-frequency domain constraint is designed in the form of an enhanced Lagrange multiplier to improve the stability of the algorithm and guarantee the synchronicity of spanned-domain calculation. Afterwards, we introduce a frequency-domain filter [18] from the field of image processing to eliminate frequency-domain redundant pieces of information. Finally, a multi-round optimization method combined with SVD is developed to alleviate the reliance on initial guesses. The proposed method is evaluated with synthetic data and real data concurrently. The results show that the inversion ability of relative permittivities of CSI is equal to 2.2 whereas that of ICSI is improved to 2.6 without increasing computational complexity. Furthermore, the capability against AWGN (additive white Gaussian noise) can be increased from 20% to 60%.

The structure of this paper is as follows. Section 2 introduces the formulation and the CSI method. Section 3 introduces the ICSI and algorithm flow. Section 4 gives the results and analysis of the numerical experiments. Section 5 presents conclusions and possibilities for future research.

2. Problem formulation and CSI method

Let us consider a 2-D inhomogeneous domain of interest (DOI), $D\subset R^2$, which is embedded in the background of an unbounded homogeneous free space. There is an unknown scatterer with arbitrary bounded cross-section located in $D$, as shown in Fig. 1. The vacuum permittivity, vacuum permeability, and wave number of background are $\varepsilon _1$, $\mu _1$, and $k_1=\omega \sqrt {\mu _1\varepsilon _1}$. The relative permittivities, permeabilities, and wave number of the unknown scatterer are represented as $\varepsilon _2(\mathbf{r} )$, $\mu _2$($\mu _2=\mu _1$), and $k_2=\omega \sqrt {\mu _2\varepsilon _2(\mathbf{r})}$, $\mathbf{r}\subset D$, where $\omega$ is the angular frequency, and material properties of the unknown scatterer are assumed constant with respect to the ${z}$-axis. Our purpose is to reconstruct the relative permittivities of the unknown scatterer. Such a scenario is illuminated by a set of $N_i$ known incident transverse magnetic (TM) sources whose electric field (incident field) has a time-dependent factor that is equal to exp$\left (-i\omega t \right )$ (omitted hereinafter). The sources are uniform located at $\mathbf{r}_j^i$ with $j=1,2,\ldots,N_i$ in a domain (or on a curve) $S$ outside of and surrounding $D$, where $S$ is the measurement domain. For each source, there are $N_r$ antennas for measurement, and the antennas are uniform located at $\mathbf{r}_q^r$ with $q=1,2,\ldots,N_r$ on the curve of $S$.

 figure: Fig. 1.

Fig. 1. The schematic diagram of two-dimensional EISP.

Download Full Size | PDF

2.1 Formulation

When there is no scatter in free-space, the measurement of the electromagnetic field is the incident field $\mathbf{E^{inc}}(\mathbf{r})$. On the contrary, if a scatter is put in free space, then the electromagnetic field will be affected and the measurement is the total field $\mathbf{E^{tot}}(\mathbf{r})$. Therefore, the scattered field $\mathbf{E^{sca}}(\mathbf{r})$ can be seen as the difference between the total field and the incident field [19]. In other words, the total field can be expressed as

$$\mathbf{E^{tot}}(\mathbf{r})=\mathbf{E^{inc}}(\mathbf{r})+\mathbf{E^{sca}}(\mathbf{r}).$$

Free-space scattering based on Maxwell‘s wave equation is formulated as

$$\bigtriangledown{\times} \bigtriangledown{\times} \mathbf{E^{tot}}(\mathbf{r})-k_1^2\mathbf{E^{tot}}(\mathbf{r})=k_1^2\mathbf{J}(\mathbf{r}),\;\; \mathbf{r}\in D\cup S$$
where $\mathbf{J}(\mathbf{r})=\mathbf{\chi } \mathbf{E^{tot}}(\mathbf{r})$ is the contrast source, $\mathbf{\chi }$ is the contrast, and $\bigtriangledown$ is the Hamilton operator. The solution of (2) consists of two parts: one is a general solution that is equivalent to the solution of the homogenous Helmholtz equation which corresponds to the right of (2) equals zero, and the other is a special solution of the non-homogeneous Eq. (2). The two parts represent the incident and scattered field respectively. Notice that all possible special solutions of (2) should be zero at infinity from the scatterer to ensure the conservation of energy. By integrating free-space and combining (2) with (1), we can describe (2) with two equations
$$\mathbf{E^{tot}}(\mathbf{r})=\mathbf{E^{inc}}(\mathbf{r})+k_1^2\int_{D}d{\mathbf{r}}'\mathbf{G}(\mathbf{r},{\mathbf{r}}')\mathbf{J}({\mathbf{r}}'),\;\;\mathbf{r}\in D$$
$$\mathbf{E^{sca}}(\mathbf{r})=k_1^2\int_{D}d{\mathbf{r}}'\mathbf{G}(\mathbf{r},{\mathbf{r}}')\mathbf{J}({\mathbf{r}}'),\;\;\mathbf{r}\in S$$
where $\mathbf{G}(\mathbf{r},{\mathbf{r}}')$ is the 2-D free-space Green function. (3) is the Lippmann Schwinger form equation which is commonly known as the state equation that describes the interaction between incident waves and the scatter. (4) is called the data equation and describes the scattered field that is re-radiated by the contrast source in the $S$ domain. However, it is difficult to get the analytical solution of the integral Eq. (3) and (4). So, we exploit the method of moments (MoM) [20] to discretize (3) and (4). Specifically, DOI domain is discretized into $M\times M$ subcells by using the pulse basis function and delta weight function and the center of subcells is located at $\mathbf{r}_n$, with $n=1,2,\ldots,r_{M^2}$. When the subcells are small enough, all the subcells can be approximated as circles of equal area, and the equivalent radius is $a=\sqrt {c/\pi }$ where $c$ is the area of a subcell. Then the relative permittivities of DOI can be regarded as the relative permittivity of each subcell. Thus, (3) and (4) are rewritten as
$$\mathbf{J}=\mathbf{\chi E^{inc}}+\mathbf{\chi G^D{J}}$$
$$\mathbf{{E}^{sca}}=\mathbf{G^S\ J}$$
where $\mathbf{{E}^{inc}}=[\mathbf{E_j^{inc}}(\mathbf{r}_1)$, $\mathbf{E_j^{inc}}(\mathbf{r}_2)$, $\cdots$, $\mathbf{E_j^{inc}}(\mathbf{r}_{M^2})]^T$, $j=1 \cdots N_i$, $\mathbf{r}\in D$, $\mathbf{{E}^{sca}}=[\mathbf{E_j^{sca}}(\mathbf{r}_1)$, $\mathbf{E_j^{sca}}(\mathbf{r}_2)$, $\cdots$, $\mathbf{E_j^{sca}}(\mathbf{r}_{M^2})]^T$, $\mathbf{r}\in S$, $\mathbf{J_j}=[\mathbf{J_j}(\mathbf{r}_1)$, $\mathbf{J_j}(\mathbf{r}_2)$, $\cdots$, $\mathbf{J_j}(\mathbf{r}_{M^2})]^T$, $\mathbf{r}\in D$, and $\mathbf{\chi }(\mathbf{r}_n)=k_1^2[m^2(\mathbf{r}_n)-1]$, $n=1,2,\ldots,r_{M^2}$, $m(\mathbf{r}_n)=k_2(\mathbf{r}_n)/k_1$, $\mathbf{r}\in D$ are the redefined incident field, scattered field, contrast source, and contrast, respectively, and the superscript $T$ denotes transpose operator. $\mathbf{{G}^D}(\cdot )$ and $\mathbf{{G}^S}(\cdot )$ are operators [21]. (5) and (6) are important for the description of EISP. To solve the problem is to solve the unknown $\mathbf{J}$ and $\mathbf{\chi }$ in (5) and (6) first, and then reconstruct the relative permittivities $\varepsilon _2(\mathbf{r} )$ of DOI.

2.2 Classic contrast source inversion (CSI) method

As mentioned in Section 1, CSI is a classical method for the EISP. The objective function of CSI is a linear combination of the data equation and the state equation [6]. It is given by

$$F(\mathbf{J_j},\mathbf{\chi })=\eta^S\Sigma_j \|\mathbf{E_j^{sca}}-\mathbf{G_j^SJ_j}\|_S^2+\eta^D\Sigma_j \|\mathbf{\chi E_j^{inc}}-\mathbf{J_j}+\mathbf{\chi G_j^DJ_j}\|_D^2$$
where $\| \cdot \|_S^2$ and $\| \cdot \|_D^2$ represent the norms of S and DOI domain, $\eta ^S=1/\Sigma _j \|\mathbf{E_j^{sca}} \|_S^2$ and $\eta ^D=1/\Sigma _j \|\mathbf{\chi E_j^{inc}} \|_D^2$ are the normalization terms, $\mathbf{\delta _j^{data}}=\mathbf{E_j^{sca}}-\mathbf{G_j^SJ_j\ }$ is the data residual, and $\mathbf{\delta _j^{state}}=\mathbf{\chi E_j^{inc}\ }-\mathbf{J_j\ }+\mathbf{\chi G_j^DJ_j}$ is the state residual. Although CSI is able to inverse, it updates the contrast source in the spatial domain, so the speed of CSI is slow.

3. Improved contrast source inversion (ICSI) method

In this section, we elaborate on the improvements and algorithm flow of ICSI.

3.1 CSI method based on Fourier transform

Motivated by [15,17], we resort Fourier bases to denote the contrast source and update the Fourier coefficients in the frequency domain, which will bring two advantages. One is that the computational complexity and memory requirements are reduced in frequency domain [15]. Another point is that after the implementation of Fourier transform, we can introduce frequency domain filtering technology to reduce the interference of redundant information. The contrast source is represented as

$$\mathbf{J\ }=vec[\mathrm{IDFT}(\mathbf{\beta } )]=\mathcal{F}^{{-}1}(\mathbf{\beta } )$$
where $\mathbf{\beta }$ is a 2-D discrete Fourier coefficient tensor matrix. The elements in the matrix are the Fourier coefficients of the contrast source. Generally, the elements close to zero correspond to the high-frequency and others correspond to the low-frequency. IDFT is short for inverse discrete Fourier transform. It is denoted as $\mathcal {F}^{-1}=[\mathbf{F_1\ },\mathbf{F_2\ },\ldots,\mathbf{F_q\ },\ldots ]$, $\mathbf{F_q\ }=vec \{ \mathbf{F_{u,v}\ }\}$, where $\mathbf{F_{u,v}\ }$ is a $M\times M$ matrix whose elements are $F_{u,v}(m,n)=exp[2\pi (mu/M+mv/M)]/M$, $m=1,2,\dots,M$, $n=1,2,\dots,M$. Then, the objective function (7) is improved to
$$F(\mathbf{\beta_j\ },\mathbf{\chi } )=\eta^S\Sigma_j \|\mathbf{E_j^{sca}\ }-\mathbf{G_j^S\ }\mathcal{F}^{{-}1}(\mathbf{\beta_j\ })\|_S^2+\eta^D\Sigma_j \|\mathbf{\chi E_j^{inc}}-\mathcal{F}^{{-}1}(\mathbf{\beta_j})+\mathbf{\chi G_j^D\ }\mathcal{F}^{{-}1}(\mathbf{\beta_j\ })\|_D^2.$$

Although Fourier transform improves the calculation efficiency, it ignores "synchronicity" which means concurrently obtaining the optimal solution in the spatial and frequency domain. To be specific, the calculation of the contrast source is carried out in two steps: the first is to calculate $\mathbf{\beta }$ in the frequency domain according to (9), and the second is to calculate $\mathbf{J}$ in the spatial domain by using (8). It is hard to ensure that the final results of $\mathbf{\beta }$ and $\mathbf{J\ }$ are concurrently optimal in both domains after the spanned-domain calculation due to the error produced by the implementation of transform. This problem also happens in [22] and other papers involving spanned-domain calculation. In addition, like other inversion methods [17,23], the inversion method with (9) still has problems such as redundant inversion information and reliance on initial guesses.

3.2 Spatial-frequency domain constraint

To achieve the synchronicity, we design a spatial-frequency domain constraint $\mathcal {F}^{-1}(\mathbf{\beta })-\mathbf{J\ }=0$, and add it to the objective function (9) in the form of an augmented Lagrange multiplier. The final objective function of ICSI is given by

$$\begin{aligned} L_A(\mathbf{\beta_j\ },\mathbf{J_j\ },\mathbf{\chi },\mathbf{\lambda } )=&\eta^S\Sigma_j \left \|\mathbf{E_j^{sca}\ }-\mathbf{G_j^S\ }\mathcal{F}^{{-}1}(\mathbf{\beta_j\ })\right \|_S^2\\ &+\eta^D\Sigma_j \left \| \mathbf{\chi E_j^{inc}\ }-\mathcal{F}^{{-}1}(\mathbf{\beta_j})+\mathbf{\chi G_j^D\ }\mathcal{F}^{{-}1}(\mathbf{\beta_j\ }) \right \|_D^2\\ & -\mathbf{\lambda }^T[\mathcal{F}^{{-}1}(\mathbf{\beta_j\ })-\mathbf{J_j}]+\frac{\mu}{2}\left \| \mathcal{F}^{{-}1}(\mathbf{\beta_j\ })-\mathbf{J_j\ } \right \|^2 ,\\ \end{aligned}$$
where $-\mathbf{\lambda }^T[\mathcal {F}^{-1}(\mathbf{\beta _j\ })-\mathbf{J_j\ }]$ and $\frac {\mu }{2}\left \| \mathcal {F}^{-1}(\mathbf{\beta _j\ })-\mathbf{J_j\ } \right \|^2 \nonumber$ are the multiplier term and the penalty term respectively. Here, $\mathbf{\lambda }$ is a multiplier factor that updates with the iteration of the ICSI algorithm, and $\mu$ is a penalty factor with an empirical assignment range [0.01, 0.1]. The motivation for using the augmented Lagrange multiplier method is as follows. First, it is easy to transform a constrained optimization problem into an unconstrained optimization problem with the method. Second, the multiplier term can measure the overlap of the gradient vectors of the objective function (9) and the spatial-frequency domain constraint. The ICSI objective function will be the minimum when the two gradients overlap. Additionally, the multiplier term can also alleviate the excessive penalty from the penalty term. Last, the penalty term can provide a good correction effect, and make the ICSI algorithm converge faster.

To simplify the calculation, we utilize ADMM [12] to decouple the multi-parameter optimization problem (10) into several single-parameter alternating optimization problems (11). It is similar to the idea of original CSI [11], which updates the contrast source first, and then brings the updated contrast source into the objective function (10) to update the contrast. Corresponding to our problem, we take turns updating the parameters in (11). Specifically, the first step is to update the Fourier coefficient $\mathbf{\beta }$ in the first equation of (11) by conjugate gradient method [6]. The second step is to update the contrast source by bringing the updated $\mathbf{\beta }$ into the second equation of (11). The third step is to update the contrast in the third equation of (11) by the conjugate gradient method and Brent’s method [24,25] under the above updated parameters. Finally, we update $\mathbf{\lambda }$ of the fourth equation of (11) with the updated parameters.

$$\begin{aligned} \mathbf{\beta_{j,n+1}} &= arg\underset{\mathbf{\beta_j}}{min}L_A(\mathbf{\beta_{j,n}\ },\mathbf{J_{j,n}\ },\mathbf{\chi_n\ },\mathbf{\lambda_n\ })\\ \mathbf{J_{j,n+1}} &= arg\underset{\mathbf{J_j}}{min}L_A(\mathbf{\beta_{j,n+1}\ },\mathbf{J_{j,n}},\mathbf{\chi_n\ },\mathbf{\lambda_n\ })\\ \mathbf{\chi_{n+1}} &= arg\underset{\mathbf{\chi}}{min}L_A(\mathbf{\beta_{j,n+1}\ },\mathbf{J_{j,n+1}},\mathbf{\chi_n\ },\mathbf{\lambda_n\ })\\ \mathbf{\lambda_{n+1}} &= \mathbf{\lambda_{n}}-\mu[\mathcal{F}^{{-}1}(\mathbf{\beta_{j,n+1}})-\mathbf{J_{j,n+1}} ] .\\ \end{aligned}$$

3.3 Low-pass filter in the frequency domain

To eliminate the redundant inversion information, we consider using frequency domain low-pass filters that are used in digital image processing [18] to eliminate high-frequency noise and retain low-frequency information. Corresponding to our problem, the Fourier coefficient matrix of contrast source $\mathbf{\beta }$ is solved in frequency domain, and its optimal solution is mostly in the low-frequency range. However, the redundant inversion information generated by measurement noise and multiple-scattering will distribute in the high-frequency range after Fourier transform, which will affect the search for the optimal solution. Therefore, it is reasonable to resort low-pass filters to eliminate the redundant information to reduce the effect. It is worth noting that some details produced by the contour of the scatter also distribute in the high-frequency range after Fourier transform. These details should be retained to provide $\mathbf{\beta }$ with an appropriate search range when we use low-pass filters.

[17] adopts a filter similar to the ideal low pass filter (ILPF) to eliminate the redundant information as shown in Fig. 2. In the centralized spectrum centered on the origin, ILPF determines a circle with its filter radius ($D0$ in Fig. 2, which is also known as the cutoff frequency), which means the low-frequency signals inside the circle can pass through the filter while the high-frequency signals outside the circle can’t. However, as mentioned above, some high-frequency details should be retained and they may distribute near the cutoff frequency. But ILPF completely truncates all high-frequency signals outside the circle, which causes the necessary high-frequency details to be lost. Therefore, we introduce a filter which has a smooth transition near the cutoff frequency instead of a step change, and call it a transition filter(TF), such as Gaussian low pass filter(GLPF), Butterworth low pass filter(BLPF), etc. Here we introduce BLPF in particular, because the performance of BLPF is compatible with ILPF and GLPF. BLPF’s transfer function is $H(u,v)=[1+D(u,v)/D0]^{-2n}$, where $D(u, v)$ is the distance between the spectrum point $(u, v)$ and the center of the spectrum, $n$ is the order. When $n$ is large, BLPF approximates ILPF, on the contrary, it approximates GLPF. The setting of BLPF parameters mainly depends on experience. By observing the centralized spectrum of Fourier coefficients of contrast sources, we suggest that $D0=0.2\pm 0.1$ and $n\in [1,6]$. The radial profile of BLPF and ILPF is shown in Fig. 2.

 figure: Fig. 2.

Fig. 2. The radial profile of BLPF and ILPF. The abscissa represents $D(u,v)$ is the distance between the spectrum point $(u,v)$ and the center of the spectrum, where $D_0$ is related to the cutoff frequency. ${H}$ denotes the operation of low-pass filters.

Download Full Size | PDF

3.4 Multi-round optimization combined with SVD

Whether an inversion method can converge to the global optimal solution is closely related to initial guesses. An inappropriate initial guess may lead to an inaccurate search direction, resulting in convergence to the local optimization solution [14]. To ease the reliance on initial guesses, inspired by SVD [7] and multi-round optimization [17], we propose a multi-round optimization strategy combined with SVD. The research process and implementation details are as follows. First, $\mathbf{G^s}$ is denoted by SVD, $\mathbf{G^s}={\sum }_m\mathbf{\mu _m}{\mathbf{\sigma _m}} \mathbf{\nu _m}^*$, where $\mathbf{\mu _m}$, $\mathbf{\sigma _m}$ and $\mathbf{\nu _m}$ are left singular vector, singular value, and right singular vector, respectively. $*$ denotes Hermite operator. Then the initial guess of the contrast source Fourier coefficients can be represented by the first L singular values from SVD of $\mathbf{G^s}$

$$\mathbf{\beta}=DFT(\underset{m=1}{\overset{L}{\sum}}\frac{\mathbf{\mu_m}^*\mathbf{E^{sca}}}{\mathbf{\sigma _m}} \mathbf{\nu _m}),$$
where DFT is the discrete Fourier transform [7]. In addition, the initial guess of the contrast $\mathbf{\chi }$ is calculated by the backpropagation method (BP) [26] with $\mathbf{\beta }$ in (12), $\chi =\underset {j=1}{\overset {N_i}{\sum }}\mathcal {F}^{-1}(\mathbf{\beta _j\ })[\mathbf{E_{j}^{tot}}]^*/\underset {j=1}{\overset {N_i}{\sum }}| \mathbf{E_{j}^{tot}}|^2$. However, when we conduct experiments, we found that as the relative permittivities increase, the interference of multi-scattering will increase and the contrast error of BP calculation will become larger, all of which will lead to the failure of the inversion. Subsequently, we develop a multi-round optimization strategy with SVD. This strategy means using above mentioned methods for the initial guess calculation of the first round optimization to obtain a rough inversion result, and then taking the updated $\mathbf{\beta }$ and $\mathbf{\chi }$ of the first round as the initial guesses of the second round to obtain a fine profile.

3.5 Algorithm flow of ICSI

The flowchart of ICSI is shown in Fig. 3 and the implementation details are as follows: Initialize the contrast and contrast source according to Section 3.4. Calculate $\mathbf{\delta ^{data}_{j,0}}$ and $\mathbf{\delta ^{state}_{j,0}}$, and let $\mathbf{\lambda _0}=0$. Initialize the gradient and search direction of contrast source as $\mathbf{g^J_{j,0}}=0$ and $\mathbf{\nu ^J_{j,0}}=0$, and do the same for the contrast as $\mathbf{g^\chi _{j,0}}=0$ and $\mathbf{\nu ^\chi _{j,0}}=0$.

 figure: Fig. 3.

Fig. 3. Flowchart of the ICSI procedure.

Download Full Size | PDF

1) Update the Fourier coefficients of contrast source $\mathbf{\beta }$. Firstly,calculate gradient of (10) with respect to $\mathbf{\beta }$:

$$\scalebox{0.93}{$\displaystyle\mathbf{g^{\beta}_{j,n}}=\frac{1}{M^2}\mathcal{F}\left \{ -\eta ^S{ \mathbf{G^S_j}}^H \mathbf{\delta _{n-1}^{data}} -\eta ^D_{n-1}[\mathbf{\delta _{n-1}^{state}}-{ \mathbf{G^D_j}}^H (\mathbf{\chi_{n}\hbox{-}\mathrm{1}}^H \mathbf{\delta _{n-1}^{state}})] +\mu[(\mathcal{F}^{{-}1}(\mathbf{\beta})-\mathbf{J})-\frac{\mathbf{\lambda_n}}{\mathbf{\mu}}]\right \},$}$$
superscript $H$ is a conjugate transpose operator. Secondly, determine the search direction with the Polak-Ribi$\acute {e}$re CG directions [27],
$$\mathbf{\nu ^\beta_{j,n}}=\mathbf{g\ ^\beta_{j,n}\ }+\frac{Re\Sigma_j \left \langle \mathbf{g\ ^\beta_{j,n}},\mathbf{g\ ^\beta_{j,n}}-\mathbf{g\ ^\beta_{j,n-1}} \right\rangle_D}{\Sigma_j\left \| \mathbf{g\ ^\beta_{j,n-1}\ }\right \|_D^2} \mathbf{\nu ^\beta_{j,n-1}}, \quad n \geq 1.$$

Thirdly, calculate the step size

$$\mathbf{\alpha_{j,n}^{\beta}}=\frac{-\operatorname{Re} \Sigma_{j}\left\langle \mathbf{g_{j,n}^{\beta}}, \mathbf{v_{j,n}^{\beta}}\right\rangle_{D}}{\eta^{S} \Sigma_{j}\left\|\mathbf{G_{j}^{S}} \mathcal{F}^{{-}1} \mathbf{v_{j,n}^{\beta}}\right\|_{S}^{2}+\eta_{n-1}^{D} \Sigma_{j}\left\|\mathcal{F}^{{-}1} \mathbf{v_{j,n}^{\beta}}-\mathbf{\chi_{n}\hbox{-}\mathrm{1}\ G_{j}^{D}} \mathcal{F}^{{-}1} \mathbf{v_{j,n}^{\beta}}\right\|_{D}^{2}+\mu\left\|\mathcal{F}^{{-}1} \mathbf{v_{j,n}^{\beta}}\right\|^{2}}.$$

Fourth, update the Fourier coefficients, $\mathbf{\beta _{j,n}}=\mathcal {H}\left (\mathbf{\beta _{j,n}\hbox{-}\mathrm{1}}+\mathbf{\alpha _{j,n}^{\beta }\ v_{j,n}^{\beta }}\right )$.

2) Update contrast source with least square method, $\mathbf{J_{j,n}}=\mathcal {F}^{-1} \mathbf{\beta _{j,n}}-\frac {\mathbf{\lambda _n}}{\mu }$. Then, update the total field $\mathbf{E^{tot}_{j,n}}=\mathbf{E^{inc}_{j}}+\mathbf{G^{D}_{j}J_{j,n}}$, and the residual error of state equation $\mathbf{\delta _{n}^{state\ }}=\mathbf{\chi _{n}\hbox{-}\mathrm{1}\ E_{j,n}^{tot}}-\mathbf{J_{j,n}}$.

3) Update the contrast $\mathbf{\chi }$. Firstly, calculate gradient of (10) with respect to $\mathbf{\chi }$,

$$\mathbf{g_{j,n}^{\chi}}=\eta_{n-1}^{D} \frac{\sum_{j} \mathbf{E_{j,n}^{tot}\ \delta_{n}^{state\ }}}{\sum_{j}\left|\mathbf{E_{j,n}^{tot}}\right|^{2}}.$$

Secondly, determine the search direction with the Polak-Ribi$\acute {e}$re CG directions,

$$\mathbf{v_{j,n}^{\chi}}=\mathbf{g_{j,n}^{\chi}}+\frac{\operatorname{Re} \Sigma_{j}\left\langle \mathbf{g_{j,n}^{\chi}}, \mathbf{g_{j,n}^{\chi}}-\mathbf{g_{j,n}\hbox{-}\mathrm{1}^{\chi}}\right\rangle_{D}}{\Sigma_j\left\|\mathbf{g_{j,n}\hbox{-}\mathrm{1}^{\chi}}\right\|_{D}^{2}} \mathbf{v_{j,n}\hbox{-}\mathrm{1}^{\chi}}, \quad n \geq 1$$

Thirdly, we calculate the step size by Brent’s method [25], which is an efficient way to find the minimum value of a single-variable quadratic function.

$$\alpha^\chi_{n}=Brent[L_{A}(\alpha^\chi_{n-1} )] = Brent[\frac{\Sigma_{j}\left\|\left(\mathbf{\chi_{n-1}}+\alpha_{ n-1}^{\chi} \mathbf{v_{j,\ n}^{\chi}}\right) \mathbf{E_{j,\ n}^{tot}}-\mathbf{J_{j,\ n}}\right\|_{D}^{2}}{\Sigma_{j}\left\|\left(\mathbf{\chi_{n-1}}+\alpha_{ n-1}^{\chi} \mathbf{v_{j,\ n}^{\chi}}\right) \mathbf{E_{j}^{inc}}\right\|_{D}^{2}}].$$

Fourth, update the contrast, $\mathbf{\chi _{n}\ }=\mathbf{\chi _{n}\hbox{-}\mathrm{1}}-\alpha _{ n}^{\chi } \mathbf{v_{j,n}^{\chi }}$.

4) Update the panel factor $\mathbf{\lambda _n}$ as shown in Eq. (11)

5) Quit the iteration when the termination conditions(the maximum number of iterations or threshold of relative error) are matched, or return to step 1).

4. Numerical results

In this section, synthetic and real data are used to assess the performance of ICSI. We consider two synthetic data: lossless homogeneous scatterers and lossy non-homogeneous scatterers. The real data are the FoamDielExt and FoamTwinDiel data from Institute Fresnel [28]. In addition to the regular inversion, anti-noise, and time consumption experiments, we also conduct several discussion experiments to evaluate the effects of different filters or the total variational regularization. To numerically appraise our method, the reconstruction relative error is defined as $err=\left \|\mathbf{\chi _t}-\mathbf{\chi _r} \right \|_2/\left \| \mathbf{\chi _t} \right \|_2$ [21], where, $\mathbf{\chi _r}$ and $\mathbf{\chi _t}$ are the reconstructed contrast and the true contrast.

4.1 Experiments on lossless homogeneous scatters

Using "Austria" [21] as a lossless uniform scatterer is a well-known challenge in the field of inverse scattering. The "Austria" consists of two small discs and a large ring. The discs have a radius of 0.2 m and are centered at $(0.3, 0.6)$ m and $(-0.3, 0.6)$ m respectively. The ring has an inner radius of 0.3 m and an exterior radius of 0.6 m.

The experimental setup is as follows. The relative permittivity of free space $\varepsilon _1$ is 1. The 16 sources ($N_i=16$) and 32 receivers ($N_r=32$) are equally placed on a circle of radius 3 m centered at (0,0) m. The frequency of the incident wave is 400MHz. DOI is a square test domain of $2\times 2$ m$^2$ and is discretized into $64 \times 64 (M=64)$ grids.

Next, the experimental results and analysis of inversion ability, anti-noise ability, filter selection, the effect of total variational regularization, and time consumption will be presented. For simplification, the four improvements of ICSI, i.e., Fourier transform, spatial-frequency domain constraint, Low-pass filter, and multi-round optimization combined with SVD, are abbreviated as FFT, SFDC, LPF, and MROSVD, respectively. To verify the effectiveness of each improvement in the following experiments, we denote our ICSI method without a Low-pass filter as ICSI-w/oLPF (CSI+FFT+SFDC+MROSVD), our ICSI method containing only FFT as CSI+FFT, and our ICSI method containing FFT and SFDC as CSI+FFT+SFDC.

4.1.1 Inversion ability

Fig. 4 and 5 summarize the comparison results between ICSI and CSI methods with regard to inversion ability. From the figures, we have the following observations: 1) The effect of CSI+FFT (Fig. 4(c)) is not as good as that of CSI (Fig. 4(b)). As described in Section 3.1, FFT can speed up the computation but does not always satisfy the synchronicity after the cross-domain computation. 2) The result shown in Fig. 4(d) obviously indicates that the constraint SFDC can help improve the inversion ability. 3) Fig. 4(e) shows a better result, which proves that MROSVD works. 4) When the relative permittivity exceeds 2.2, the relative permittivity of CSI inversion is not only far greater than the value we set, but also seriously distorted, as shown in Fig. 5(a) and (d). However, our ICSI can still roughly reconstruct the profile, as shown in Fig. 5(b), (c) and (f). 5) The results of Fig. 4(e), Fig. 5(b) and Fig. 5(e) mean that the interference from multiple scattering increases as the relative permittivity increases, which is discussed in Section 3.3. Fig. 4(f), Fig. 5(c) and Fig. 5(f) show LPF can indeed help ICSI reduce this interference to improve the inversion ability.

 figure: Fig. 4.

Fig. 4. Comprisions of inversion ability when the relative permittivity is 2.2. (a) True profile. (b) CSI(err=0.1485). (c)CSI+FFT(err=0.1728). (d) CSI+FFT+SFDC(err=0.1412). (e) ICSI-w/oLPF(err=0.1362). (f) ICSI(err=0.1226).

Download Full Size | PDF

 figure: Fig. 5.

Fig. 5. Reconstruction results of CSI and ICSI when relative permittivity(RP) is 2.4 and 2.6. (a) CSI (BP=2.4, err=0.4966). (b) ICSI-w/oLPF (RP=2.4, err=0.1513). (c) ICSI (RP=2.4, err=0.1498). (d) CSI (RP=2.6, err=0.6128). (e) ICSI-w/oLPF (RP=2.6, err=0.3418). (f) ICSI (RP=2.6, err=0.1891).

Download Full Size | PDF

4.1.2 Anti-noise ability

This subsection investigates the anti-noise ability of our proposed ICSI method and compares it with CSI method. Fig. 6 and 7 report the results when the relative permittivity is 2.0. According to the results, we can notice that: 1) When the addition of AWGN is 10%, the anti-noise ability of ICSI (Fig. 6(c)) is much better than that of CSI (Fig. 6(a)). And even without a filter, our ICSI-w/oLPF (Fig. 6(b)) is also better than CSI. 2) As AWGN increases to 20%, the anti-noise ability of CSI (Fig. 6(d)) gets worse, but our methods (Fig. 6(e) and Fig. 6(f)) still have a promising performance. 3) When we continue to add AWGN, we find that CSI is gradually distorted, as shown in Fig. 7(d), Fig. 7(e) and Fig. 7(f). Especially when AWGN=60%, our method can still roughly reconstruct the outline of the scatterer, as shown in Fig. 7(c), but CSI is distorted as shown in Fig. 7(f). From the above experimental results, the increase of noise leads to the increase of reconstruction error. This indicates that noise has a significant impact on the reconstruction quality, which requires that the model should have higher robustness. Another discovery provides a rationale for the selection of BLPF parameters. With the increase of noise, the cut-off frequency $D0$ and order $n$ of the filter should gradually decrease.

 figure: Fig. 6.

Fig. 6. Comparisons of the anti-noise ability when the AWGN levels are 10% and 20%. (a) CSI(AWGN = 10%, err=0.1793). (b) ICSI-w/oLPF(AWGN = 10%, err=0.1495). (c) ICSI(AWGN = 10%, D0=0.3, n=5, err=0.1125). (d) CSI(AWGN = 20%, err=0.1885). (e) ICSI-w/oLPF(AWGN = 20%, err=0.1654). (f) ICSI(AWGN = 20%, D0=0.3, n=5, err=0.1267).

Download Full Size | PDF

 figure: Fig. 7.

Fig. 7. Reconstruction results of high AWGN levels. (a) ICSI with AWGN=40%(D0=0.2, n=4, err=0.1493). (b) ICSI with AWGN=50%(D0=0.2, n=4, err=0.1603). (c) ICSI with AWGN=60%(D0=0.2, n=3, err=0.1804). (d) CSI with AWGN=40%(err=0.1902). (e) CSI with AWGN=50%(err=0.1922). (f) CSI with AWGN=60%(err=0.1985).

Download Full Size | PDF

4.1.3 Comparison between different filters

This subsection studies the effects of different filters, such as the ideal low-pass filter (ILPF) and Butterworth low-pass filter (BLPF) . In the following experiments, the relative permittivity is 2.2, and the AWGN = 20%. From the presented results in Fig. 8, we observe that: 1) The effect of ILPF (Fig. 8(b)) is better than ICSI-w/oLPF (Fig. 8(a)). Since ILPF can filter out redundant information in high-frequency, as mentioned in Section 3.3. 2) By using BLPF, we get a clearer and more complete contour than that of ILPF, as shown in Fig. 8(c). This is because BILP belongs to the TF. TF has a transition at the cut-off frequency, it can preserve some useful details around the cut-off frequency, as described in Section 3.3. However, ILPF would discard these useful details. In addition, BLPF has the maximum flat amplitude frequency response. The flatness ensures that the original signal would not be attenuated by the filtering.

 figure: Fig. 8.

Fig. 8. Performance of different filters. (a) The performance of ICSI-w/oLPF, which is without filters(err=0.1654)). (b) The performance of ICSI using ILPF(D0=0.3, err=0.1351). (c) The performance of ICSI using BLPF(D0=0.3, n=4, err=0.1267)).

Download Full Size | PDF

4.1.4 Effect of the total variational regularization on ICSI

Similar to [6,29], we intend to study the effect of the total variational (TV) regularization in this subsection. A $L^1$-norm TV regularization was added to (10), then the objective function of a TV-enhanced ICSI is rewritten as

$$\begin{aligned} L_A(\mathbf{\beta_j\ },\mathbf{J_j\ },\mathbf{\chi },\mathbf{\lambda } )&= \eta^S\Sigma_j \left \|\mathbf{E_j^{sca}\ }-\mathbf{G_j^S\ }\mathcal{F}^{{-}1}(\mathbf{\beta_j\ })\right \|_S^2\\ &+\eta^D\Sigma_j \left \| \mathbf{\chi E_j^{inc}\ }-\mathcal{F}^{{-}1}(\mathbf{\beta_j})+\mathbf{\chi G_j^D\ }\mathcal{F}^{{-}1}(\mathbf{\beta_j\ }) \right \|_D^2\\ & -\mathbf{\lambda_1\ }^T[\mathcal{F}^{{-}1}(\mathbf{\beta_j\ })-\mathbf{J_j}]+\frac{\mu}{2}\left \| \mathcal{F}^{{-}1}(\mathbf{\beta_j\ })-\mathbf{J_j\ } \right \|^2\\ & +\left \| \mathbf{w} \right \|_1 -\mathbf{\lambda_{2}}^{T}\left(\mathbf{w}-D \mathbf{\chi_{j}}\right)+\frac{\mu_{2}}{2}\left\|\mathbf{w\ }-D \mathbf{\chi_{j}}\right\|_{2}^{2} .\\ \end{aligned}$$

From the results in Fig. 9, the TV-enhanced ICSI brings obvious improvement, compared to Fig. 5(c) and (d). Especially when the relative permittivity is 2.4, Fig. 9(a) exhibits a continuous and undistorted contour, and the inversion value of relative permittivity is very close to the preset true value. This is because the $L^1$-norm TV regularization can conduce to the contrast piecewise constant. In other words, it is a priori information that can alleviate the pathology of the mathematical problems in EISP.

 figure: Fig. 9.

Fig. 9. Results of the TV-enhanced ICSI. (a) The relative permittivity is 2.4(err=0.1266). (b) The relative permittivity is 2.6(err=0.1834).

Download Full Size | PDF

4.1.5 Comparison of time consumption

In this subsection, we discuss the time consumption comparison between our ICSI method and the CSI method. The comparison is conducted on a computer with AMD R7 CPU, 16G RAM, and MATLAB 2017a. The relative permittivity is 2.2.

The results in Table 1 show that our proposed ICSI method consumes less computational time than the CSI method. It demonstrates that the spanned-domain calculation based on FFT indeed makes ICSI computationally efficient.

Tables Icon

Table 1. Comparison of time consumption

4.1.6 Comparison of initialization

In this part, we compare the BP initialization method with our initialization method (Section 3.4). To be specific, we apply these two initialization methods to CSI and ICSI inversion methods respectively, and the results are shown in Fig. 10, where the relative permittivity is set to 2.2. From the figures, it can be seen that: 1) Compared with the original combination of CSI and BP(Fig. 10(b)), our initialization method works well with CSI, and helps CSI to obtain better reconstruction results and smaller relative errors, as shown in Fig. 10(a). 2) When these two initialization methods are applied to ICSI, our initialization method is also better,as shown in Fig. 10(c) and Fig. 10(d). It is obviously that our initialization method has better generality and more suitable for estimating the initial guess.

 figure: Fig. 10.

Fig. 10. Comparison of two initialization methods. (a) CSI uses our method as the initial guess for inversion (L=20, err=0.1396). (b) CSI uses BP as the initial guess for inversion (err=0.1485). (c) ICSI uses our method as the initial guess for inversion (L=20, err=0.1226). (d)ICSI uses BP as the initial guess for inversion (err=0.1402).

Download Full Size | PDF

4.2 Experiments on lossy non-homogeneous scatterers

In this section, we conduct experiments on lossy non-homogeneous scatterers to verify the reconstruction effect of ICSI and compare it with CSI. The true profile consists of two discs and a coated rectangle. The discs have a radius of 0.3 m, and their centers are (-0.4, 0.6) m and (0.4, 0.6) m respectively. The relative permittivity of the left disc is 2.5 and that of the right disc is 3.0. The center of the coated rectangle is located at (0,-0.3) m, the sides of the inner rectangle are 0.6 m and 1 m, and those of the outer are 1 m and 1.6 m. The relative permittivity of the inner rectangle is 2.0+i0.5 and that of the outer is 1.5+i0.2. The remainder of the experimental setup is the same as Section 4.1. Fig. 11 and Fig. 12 show the results of ICSI and CSI respectively, where the inversion results for the real part of the relative permittivity are in the first row and the results for the corresponding imaginary part are in the second row. It can be seen from the results that 1) When there is no noise, both ICSI and CSI can reconstruct the profile, as shown in Fig. 11(b), (f), and Fig. 12(a), (d), but the effect and error of ICSI are better than CSI. 2) When AWGN = 30%, and even when AWGN is increased to 50%, ICSI can still roughly reconstruct the profil, as shown in Fig. 11(c), (d), (g), (h), but the results of CSI become distorted, as shown in Fig. 12(b), (c), (e), (f).

 figure: Fig. 11.

Fig. 11. Reconstruction of lossy non-homogeneous scatterers by ICSI. (a) Real part of true profile. (b) Real part inversion(non-noise, err=0.1498). (c) Real part inversion(AWGN=30%, err=0.1572). (d) Real part inversion(AWGN=50%, err=0.1762). (e) Imaginary part of true profile. (f) Imaginary part inversion(non-noise, err=0.1498). (g) Imaginary part inversion(AWGN=30%, err=0.1572). (h) Imaginary part inversion(AWGN=50%, err=0.1762).

Download Full Size | PDF

 figure: Fig. 12.

Fig. 12. Reconstruction of lossy non-homogeneous scatterers by CSI. (a) Real part inversion(non-noise, err=0.1672). (b) Real part inversion(AWGN=30%, err=0.2035). (c) Real part inversion(AWGN=50%, err=0.2803). (d) Imaginary part inversion(non-noise, err=0.1672). (e) Imaginary part inversion(AWGN=30%, err=0.2035). (f) Imaginary part inversion(AWGN=50%, err=0.2803).

Download Full Size | PDF

4.3 Experiments on real data

In this section, we consider the real data carried out by Institute Fresnel in the year 2005. We choose “FoamDielExt" and “FoamTwinDiel" to further validate ICSI in the TM cases. There are two types of cylinders, a big cylinder with a diameter of 0.08 m and the relative permittivity is $1.45\pm 0.15$, and a small cylinder (Babylon) with a diameter of 0.031 m and the relative permittivity is $3\pm 0.3$. The "FoamDielExt" consist of one larger circular dielectric cylinder with a smaller adjacent one outside, as shown in Fig. 13(a). And it is measured by 8 sources and 241 receiving antennas at different frequencies. The distance from the source and receiver to the center of the DOI is 1.67 m. The “FoamTwinDiel" consists of one larger circular dielectric cylinder with a smaller one embedded inside and a smaller adjacent one outside, as shown in Fig. 14(a). Other settings are the same as FoamDielExt. The DOI is 0.15m $\times$ 0.15m and is discretized into 64 $\times$ 64 grids. To calibrate the real data, multiply the data by a single complex coefficient derived only from the ratio of the measured incident field to the receiver’s simulated incidence field located opposite the source [28]. As can be seen from the results in Fig. 13 and 14, in the test of real data, our method still presents promising performance as usual.

 figure: Fig. 13.

Fig. 13. Reconstructed results against FoamDielExt data set at 3GHz and 5GHz. (a) True profile of FoamDielExt. (b) Reconstruction at 3GHz(err=0.1764). (c) Reconstruction at 5GHz(err=0.2269).

Download Full Size | PDF

 figure: Fig. 14.

Fig. 14. Reconstructed results against FoamTwinDiel data set at 3GHz. (a) True profile of FoamTwinDiel. (b) Reconstruction at 3GHz(err=0.1658).

Download Full Size | PDF

5. Conclusion

An improved contrast source inversion (ICSI) method is proposed to solve the electromagnetic inverse scattering problems based on CSI. The improvements include: using Fourier bases to denote the contrast source to accelerate calculation in the frequency domain; designing a spatial-frequency domain constraint to make the spanned-domain calculation synchronously optimal; exploiting a frequency domain filter to eliminate the redundant inversion information and narrow the search for optimal solutions; developing a multi-round optimization combined with SVD to alleviate reliance on initial guesses. By processing the lossless homogeneous scatterers, we demonstrate the advantage of the proposed improvements case by case and compare them with CSI. In addition, we also evaluate ICSI in the experiments with lossy non-homogeneous scatterers and the real data. In the future, we will try the combination of ANN and physical methods to solve more complex electromagnetic inverse scattering problems.

Funding

National Natural Science Foundation of China (62041302).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. P. Colli Franzone, L. Pavarino, G. Savarè, and A. Quarteroni, Complex Systems in Biomedicine (Springer, 2006).

2. O. M. Bucci, L. Crocco, T. Isernia, and V. Pascazio, “Subsurface inverse scattering problems: quantifying, qualifying, and achieving the available information,” IEEE Trans. Geosci. Remote Sensing 39(11), 2527–2538 (2001). [CrossRef]  

3. H. Kagiwada, R. Kalaba, S. Timko, and S. Ueno, “Associate memories for system identification: Inverse problems in remote sensing,” Math. Comput. Model. 14, 200–202 (1990). [CrossRef]  

4. R. Zoughi, Microwave non-destructive testing and evaluation principles, vol. 4 (Springer Science & Business Media, 2000).

5. A. Randazzo and C. Estatico, “A regularisation scheme for electromagnetic inverse problems: application to crack detection in civil structures,” Nondestr. Test. Eval. 27(3), 189–197 (2012). [CrossRef]  

6. P. M. van den Berg, A. Van Broekhoven, and A. Abubakar, “Extended contrast source inversion,” Inverse Probl. 15(5), 1325–1344 (1999). [CrossRef]  

7. X. Chen, “Subspace-based optimization method for solving inverse-scattering problems,” IEEE Trans. Geosci. Remote Sensing 48(1), 42–49 (2010). [CrossRef]  

8. Y. Zhong, M. Lambert, D. Lesselier, and X. Chen, “A new integral equation method to solve highly nonlinear inverse scattering problems,” IEEE Trans. Antennas Propag. 64(5), 1788–1799 (2016). [CrossRef]  

9. Z. Wei and X. Chen, “Deep-learning schemes for full-wave nonlinear inverse scattering problems,” IEEE Trans. Geosci. Remote Sensing 57(4), 1849–1860 (2019). [CrossRef]  

10. J. Liu, H. Zhou, T. Ouyang, Q. Liu, and Y. Wang, “Physical model-inspired deep unrolling network for solving nonlinear inverse scattering problems,” IEEE Trans. Antennas Propag. 70(2), 1236–1249 (2022). [CrossRef]  

11. P. M. Van Den Berg and R. E. Kleinman, “A contrast source inversion method,” Inverse Probl. 13(6), 1607–1620 (1997). [CrossRef]  

12. S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Foundations Trends Mach. learning 3(1), 1–122 (2010). [CrossRef]  

13. O. Semerci, M. Li, A. Abubakar, and T. Habashy, “A contrast source inversion method using truncated wavelet representations,” Proceedings of the 2012 IEEE International Symposium on Antennas and Propagation, (IEEE, 20121–2).

14. S. Sun, B. J. Kooij, T. Jin, and A. G. Yarovoy, “Cross-correlated contrast source inversion,” IEEE Trans. Antennas Propag. 65(5), 2592–2603 (2017). [CrossRef]  

15. R. Pierri, F. Soldovieri, A. Liseno, and F. De Blasio, “Dielectric profiles reconstruction via the quadratic approach in 2-d geometry from multifrequency and multifrequency/multiview data,” IEEE transactions on geoscience and remote sensing 40(12), 2709–2718 (2002). [CrossRef]  

16. X. Chen, Computational methods for electromagnetic inverse scattering (John Wiley & Sons, 2018).

17. Y. Zhong and X. Chen, “An fft twofold subspace-based optimization method for solving electromagnetic inverse scattering problems,” IEEE Trans. Antennas Propag. 59(3), 914–927 (2011). [CrossRef]  

18. R. C. Gonzalez and R. E. Woods, Digital Image Processing (3rd Edition) (Person education, 2007).

19. M. I. Mishchenko, L. D. Travis, and A. A. Lacis, Scattering, absorption, and emission of light by small particles (Cambridge university press, 2002).

20. A. Lakhtakia, “Strong and weak forms of the method of moments and the coupled dipole method for scattering of time-harmonic electromagnetic fields,” Int. J. Mod. Phys. C 03(03), 583–603 (1992). [CrossRef]  

21. X. Chen, Computational Methods for Electromagnetic Inverse Scattering (Computational Methods for Electromagnetic Inverse Scattering, 2018).

22. T. Yin, Z. Wei, and X. Chen, “Wavelet transform subspace-based optimization method for inverse scattering,” IEEE J. on Multiscale Multiphysics Comput. Tech. 3, 176–184 (2018). [CrossRef]  

23. K. Xu, L. Zhang, and Z. Wei, “Fourier bases-expansion contraction integral equation for inversion highly nonlinear inverse scattering problem,” IEEE Trans. Microwave Theory Tech. 68(6), 2206–2214 (2020). [CrossRef]  

24. R. P. Brent, Algorithms for minimization without derivatives (Courier Corporation, 2013).

25. G. E. Forsythe, “Computer methods for mathematical computations,” Prentice-Hall series in automatic computation259, (1977).

26. K. Belkebir, P. C. Chaumet, and A. Sentenac, “Superresolution in total internal reflection tomography,” J. Opt. Soc. Am. A 22(9), 1889–1897 (2005). [CrossRef]  

27. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical recipes 3rd edition: The art of scientific computing (Cambridge university press, 2007).

28. J.-M. Geffrin, P. Sabouroux, and C. Eyraud, “Free space experimental scattering database continuation: experimental set-up and measurement precision,” Inverse Probl. 21(6), S117–S130 (2005). [CrossRef]  

29. J. Yang, Y. Zhang, and W. Yin, “A fast alternating direction method for tvl1-l2 signal reconstruction from partial fourier data,” IEEE J. Sel. Top. Signal Process 4(2), 288–297 (2010). [CrossRef]  

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (14)

Fig. 1.
Fig. 1. The schematic diagram of two-dimensional EISP.
Fig. 2.
Fig. 2. The radial profile of BLPF and ILPF. The abscissa represents $D(u,v)$ is the distance between the spectrum point $(u,v)$ and the center of the spectrum, where $D_0$ is related to the cutoff frequency. ${H}$ denotes the operation of low-pass filters.
Fig. 3.
Fig. 3. Flowchart of the ICSI procedure.
Fig. 4.
Fig. 4. Comprisions of inversion ability when the relative permittivity is 2.2. (a) True profile. (b) CSI(err=0.1485). (c)CSI+FFT(err=0.1728). (d) CSI+FFT+SFDC(err=0.1412). (e) ICSI-w/oLPF(err=0.1362). (f) ICSI(err=0.1226).
Fig. 5.
Fig. 5. Reconstruction results of CSI and ICSI when relative permittivity(RP) is 2.4 and 2.6. (a) CSI (BP=2.4, err=0.4966). (b) ICSI-w/oLPF (RP=2.4, err=0.1513). (c) ICSI (RP=2.4, err=0.1498). (d) CSI (RP=2.6, err=0.6128). (e) ICSI-w/oLPF (RP=2.6, err=0.3418). (f) ICSI (RP=2.6, err=0.1891).
Fig. 6.
Fig. 6. Comparisons of the anti-noise ability when the AWGN levels are 10% and 20%. (a) CSI(AWGN = 10%, err=0.1793). (b) ICSI-w/oLPF(AWGN = 10%, err=0.1495). (c) ICSI(AWGN = 10%, D0=0.3, n=5, err=0.1125). (d) CSI(AWGN = 20%, err=0.1885). (e) ICSI-w/oLPF(AWGN = 20%, err=0.1654). (f) ICSI(AWGN = 20%, D0=0.3, n=5, err=0.1267).
Fig. 7.
Fig. 7. Reconstruction results of high AWGN levels. (a) ICSI with AWGN=40%(D0=0.2, n=4, err=0.1493). (b) ICSI with AWGN=50%(D0=0.2, n=4, err=0.1603). (c) ICSI with AWGN=60%(D0=0.2, n=3, err=0.1804). (d) CSI with AWGN=40%(err=0.1902). (e) CSI with AWGN=50%(err=0.1922). (f) CSI with AWGN=60%(err=0.1985).
Fig. 8.
Fig. 8. Performance of different filters. (a) The performance of ICSI-w/oLPF, which is without filters(err=0.1654)). (b) The performance of ICSI using ILPF(D0=0.3, err=0.1351). (c) The performance of ICSI using BLPF(D0=0.3, n=4, err=0.1267)).
Fig. 9.
Fig. 9. Results of the TV-enhanced ICSI. (a) The relative permittivity is 2.4(err=0.1266). (b) The relative permittivity is 2.6(err=0.1834).
Fig. 10.
Fig. 10. Comparison of two initialization methods. (a) CSI uses our method as the initial guess for inversion (L=20, err=0.1396). (b) CSI uses BP as the initial guess for inversion (err=0.1485). (c) ICSI uses our method as the initial guess for inversion (L=20, err=0.1226). (d)ICSI uses BP as the initial guess for inversion (err=0.1402).
Fig. 11.
Fig. 11. Reconstruction of lossy non-homogeneous scatterers by ICSI. (a) Real part of true profile. (b) Real part inversion(non-noise, err=0.1498). (c) Real part inversion(AWGN=30%, err=0.1572). (d) Real part inversion(AWGN=50%, err=0.1762). (e) Imaginary part of true profile. (f) Imaginary part inversion(non-noise, err=0.1498). (g) Imaginary part inversion(AWGN=30%, err=0.1572). (h) Imaginary part inversion(AWGN=50%, err=0.1762).
Fig. 12.
Fig. 12. Reconstruction of lossy non-homogeneous scatterers by CSI. (a) Real part inversion(non-noise, err=0.1672). (b) Real part inversion(AWGN=30%, err=0.2035). (c) Real part inversion(AWGN=50%, err=0.2803). (d) Imaginary part inversion(non-noise, err=0.1672). (e) Imaginary part inversion(AWGN=30%, err=0.2035). (f) Imaginary part inversion(AWGN=50%, err=0.2803).
Fig. 13.
Fig. 13. Reconstructed results against FoamDielExt data set at 3GHz and 5GHz. (a) True profile of FoamDielExt. (b) Reconstruction at 3GHz(err=0.1764). (c) Reconstruction at 5GHz(err=0.2269).
Fig. 14.
Fig. 14. Reconstructed results against FoamTwinDiel data set at 3GHz. (a) True profile of FoamTwinDiel. (b) Reconstruction at 3GHz(err=0.1658).

Tables (1)

Tables Icon

Table 1. Comparison of time consumption

Equations (19)

Equations on this page are rendered with MathJax. Learn more.

$$\mathbf{E^{tot}}(\mathbf{r})=\mathbf{E^{inc}}(\mathbf{r})+\mathbf{E^{sca}}(\mathbf{r}).$$
$$\bigtriangledown{\times} \bigtriangledown{\times} \mathbf{E^{tot}}(\mathbf{r})-k_1^2\mathbf{E^{tot}}(\mathbf{r})=k_1^2\mathbf{J}(\mathbf{r}),\;\; \mathbf{r}\in D\cup S$$
$$\mathbf{E^{tot}}(\mathbf{r})=\mathbf{E^{inc}}(\mathbf{r})+k_1^2\int_{D}d{\mathbf{r}}'\mathbf{G}(\mathbf{r},{\mathbf{r}}')\mathbf{J}({\mathbf{r}}'),\;\;\mathbf{r}\in D$$
$$\mathbf{E^{sca}}(\mathbf{r})=k_1^2\int_{D}d{\mathbf{r}}'\mathbf{G}(\mathbf{r},{\mathbf{r}}')\mathbf{J}({\mathbf{r}}'),\;\;\mathbf{r}\in S$$
$$\mathbf{J}=\mathbf{\chi E^{inc}}+\mathbf{\chi G^D{J}}$$
$$\mathbf{{E}^{sca}}=\mathbf{G^S\ J}$$
$$F(\mathbf{J_j},\mathbf{\chi })=\eta^S\Sigma_j \|\mathbf{E_j^{sca}}-\mathbf{G_j^SJ_j}\|_S^2+\eta^D\Sigma_j \|\mathbf{\chi E_j^{inc}}-\mathbf{J_j}+\mathbf{\chi G_j^DJ_j}\|_D^2$$
$$\mathbf{J\ }=vec[\mathrm{IDFT}(\mathbf{\beta } )]=\mathcal{F}^{{-}1}(\mathbf{\beta } )$$
$$F(\mathbf{\beta_j\ },\mathbf{\chi } )=\eta^S\Sigma_j \|\mathbf{E_j^{sca}\ }-\mathbf{G_j^S\ }\mathcal{F}^{{-}1}(\mathbf{\beta_j\ })\|_S^2+\eta^D\Sigma_j \|\mathbf{\chi E_j^{inc}}-\mathcal{F}^{{-}1}(\mathbf{\beta_j})+\mathbf{\chi G_j^D\ }\mathcal{F}^{{-}1}(\mathbf{\beta_j\ })\|_D^2.$$
$$\begin{aligned} L_A(\mathbf{\beta_j\ },\mathbf{J_j\ },\mathbf{\chi },\mathbf{\lambda } )=&\eta^S\Sigma_j \left \|\mathbf{E_j^{sca}\ }-\mathbf{G_j^S\ }\mathcal{F}^{{-}1}(\mathbf{\beta_j\ })\right \|_S^2\\ &+\eta^D\Sigma_j \left \| \mathbf{\chi E_j^{inc}\ }-\mathcal{F}^{{-}1}(\mathbf{\beta_j})+\mathbf{\chi G_j^D\ }\mathcal{F}^{{-}1}(\mathbf{\beta_j\ }) \right \|_D^2\\ & -\mathbf{\lambda }^T[\mathcal{F}^{{-}1}(\mathbf{\beta_j\ })-\mathbf{J_j}]+\frac{\mu}{2}\left \| \mathcal{F}^{{-}1}(\mathbf{\beta_j\ })-\mathbf{J_j\ } \right \|^2 ,\\ \end{aligned}$$
$$\begin{aligned} \mathbf{\beta_{j,n+1}} &= arg\underset{\mathbf{\beta_j}}{min}L_A(\mathbf{\beta_{j,n}\ },\mathbf{J_{j,n}\ },\mathbf{\chi_n\ },\mathbf{\lambda_n\ })\\ \mathbf{J_{j,n+1}} &= arg\underset{\mathbf{J_j}}{min}L_A(\mathbf{\beta_{j,n+1}\ },\mathbf{J_{j,n}},\mathbf{\chi_n\ },\mathbf{\lambda_n\ })\\ \mathbf{\chi_{n+1}} &= arg\underset{\mathbf{\chi}}{min}L_A(\mathbf{\beta_{j,n+1}\ },\mathbf{J_{j,n+1}},\mathbf{\chi_n\ },\mathbf{\lambda_n\ })\\ \mathbf{\lambda_{n+1}} &= \mathbf{\lambda_{n}}-\mu[\mathcal{F}^{{-}1}(\mathbf{\beta_{j,n+1}})-\mathbf{J_{j,n+1}} ] .\\ \end{aligned}$$
$$\mathbf{\beta}=DFT(\underset{m=1}{\overset{L}{\sum}}\frac{\mathbf{\mu_m}^*\mathbf{E^{sca}}}{\mathbf{\sigma _m}} \mathbf{\nu _m}),$$
$$\scalebox{0.93}{$\displaystyle\mathbf{g^{\beta}_{j,n}}=\frac{1}{M^2}\mathcal{F}\left \{ -\eta ^S{ \mathbf{G^S_j}}^H \mathbf{\delta _{n-1}^{data}} -\eta ^D_{n-1}[\mathbf{\delta _{n-1}^{state}}-{ \mathbf{G^D_j}}^H (\mathbf{\chi_{n}\hbox{-}\mathrm{1}}^H \mathbf{\delta _{n-1}^{state}})] +\mu[(\mathcal{F}^{{-}1}(\mathbf{\beta})-\mathbf{J})-\frac{\mathbf{\lambda_n}}{\mathbf{\mu}}]\right \},$}$$
$$\mathbf{\nu ^\beta_{j,n}}=\mathbf{g\ ^\beta_{j,n}\ }+\frac{Re\Sigma_j \left \langle \mathbf{g\ ^\beta_{j,n}},\mathbf{g\ ^\beta_{j,n}}-\mathbf{g\ ^\beta_{j,n-1}} \right\rangle_D}{\Sigma_j\left \| \mathbf{g\ ^\beta_{j,n-1}\ }\right \|_D^2} \mathbf{\nu ^\beta_{j,n-1}}, \quad n \geq 1.$$
$$\mathbf{\alpha_{j,n}^{\beta}}=\frac{-\operatorname{Re} \Sigma_{j}\left\langle \mathbf{g_{j,n}^{\beta}}, \mathbf{v_{j,n}^{\beta}}\right\rangle_{D}}{\eta^{S} \Sigma_{j}\left\|\mathbf{G_{j}^{S}} \mathcal{F}^{{-}1} \mathbf{v_{j,n}^{\beta}}\right\|_{S}^{2}+\eta_{n-1}^{D} \Sigma_{j}\left\|\mathcal{F}^{{-}1} \mathbf{v_{j,n}^{\beta}}-\mathbf{\chi_{n}\hbox{-}\mathrm{1}\ G_{j}^{D}} \mathcal{F}^{{-}1} \mathbf{v_{j,n}^{\beta}}\right\|_{D}^{2}+\mu\left\|\mathcal{F}^{{-}1} \mathbf{v_{j,n}^{\beta}}\right\|^{2}}.$$
$$\mathbf{g_{j,n}^{\chi}}=\eta_{n-1}^{D} \frac{\sum_{j} \mathbf{E_{j,n}^{tot}\ \delta_{n}^{state\ }}}{\sum_{j}\left|\mathbf{E_{j,n}^{tot}}\right|^{2}}.$$
$$\mathbf{v_{j,n}^{\chi}}=\mathbf{g_{j,n}^{\chi}}+\frac{\operatorname{Re} \Sigma_{j}\left\langle \mathbf{g_{j,n}^{\chi}}, \mathbf{g_{j,n}^{\chi}}-\mathbf{g_{j,n}\hbox{-}\mathrm{1}^{\chi}}\right\rangle_{D}}{\Sigma_j\left\|\mathbf{g_{j,n}\hbox{-}\mathrm{1}^{\chi}}\right\|_{D}^{2}} \mathbf{v_{j,n}\hbox{-}\mathrm{1}^{\chi}}, \quad n \geq 1$$
$$\alpha^\chi_{n}=Brent[L_{A}(\alpha^\chi_{n-1} )] = Brent[\frac{\Sigma_{j}\left\|\left(\mathbf{\chi_{n-1}}+\alpha_{ n-1}^{\chi} \mathbf{v_{j,\ n}^{\chi}}\right) \mathbf{E_{j,\ n}^{tot}}-\mathbf{J_{j,\ n}}\right\|_{D}^{2}}{\Sigma_{j}\left\|\left(\mathbf{\chi_{n-1}}+\alpha_{ n-1}^{\chi} \mathbf{v_{j,\ n}^{\chi}}\right) \mathbf{E_{j}^{inc}}\right\|_{D}^{2}}].$$
$$\begin{aligned} L_A(\mathbf{\beta_j\ },\mathbf{J_j\ },\mathbf{\chi },\mathbf{\lambda } )&= \eta^S\Sigma_j \left \|\mathbf{E_j^{sca}\ }-\mathbf{G_j^S\ }\mathcal{F}^{{-}1}(\mathbf{\beta_j\ })\right \|_S^2\\ &+\eta^D\Sigma_j \left \| \mathbf{\chi E_j^{inc}\ }-\mathcal{F}^{{-}1}(\mathbf{\beta_j})+\mathbf{\chi G_j^D\ }\mathcal{F}^{{-}1}(\mathbf{\beta_j\ }) \right \|_D^2\\ & -\mathbf{\lambda_1\ }^T[\mathcal{F}^{{-}1}(\mathbf{\beta_j\ })-\mathbf{J_j}]+\frac{\mu}{2}\left \| \mathcal{F}^{{-}1}(\mathbf{\beta_j\ })-\mathbf{J_j\ } \right \|^2\\ & +\left \| \mathbf{w} \right \|_1 -\mathbf{\lambda_{2}}^{T}\left(\mathbf{w}-D \mathbf{\chi_{j}}\right)+\frac{\mu_{2}}{2}\left\|\mathbf{w\ }-D \mathbf{\chi_{j}}\right\|_{2}^{2} .\\ \end{aligned}$$
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.