Contraction Integral Equation for Three-Dimensional Electromagnetic Inverse Scattering Problems

Abstract

Inverse scattering problems (ISPs) stand at the center of many important imaging applications, such as geophysical explorations, industrial non-destructive testing, bio-medical imaging, etc. Recently, a new type of contraction integral equation for inversion (CIE-I) has been proposed to tackle the two-dimensional electromagnetic ISPs, in which the usually employed Lippmann–Schwinger integral equation (LSIE) is transformed into a new form with a modified medium contrast via a contraction mapping. With the CIE-I, the multiple scattering effects, i.e., the physical reason for the nonlinearity in the ISPs, is substantially suppressed in estimating the modified contrast, without compromising physical modeling. In this paper, we firstly propose to implement this new CIE-I for the three-dimensional ISPs. With the help of the FFT type twofold subspace-based optimization method (TSOM), when handling the highly nonlinear problems with strong scatterers, those with higher contrast and/or larger dimensions (in terms of wavelengths), the performance of the inversions with CIE-I is much better than the ones with the LSIE, wherein inversions usually converge to local minima that may be far away from the solution. In addition, when handling the moderate scatterers (those the LSIE modeling can still handle), the convergence speed of the proposed method with CIE-I is much faster than the one with the LSIE. Secondly, we propose to relax the contraction mapping condition, i.e., different contraction mappings are used in updating contrast sources and contrast, and we find that the convergence can be further accelerated. Several numerical tests illustrate the aforementioned interests.

1. Introduction

Inverse scattering problems (ISPs) in electromagnetics and acoustics are of great interest in industries due to important imaging applications in various areas, such as geophysical survey, non-destructive testing, ground-penetrating radar, bio-medical imaging, etc., as solving the ISPs provides rich information about the unknown targets, such as locations, shapes, and material distributions within some structures [1,2,3,4]. For instance, microwave imaging has been used to inspect the abnormalities in human bodies like bleeding in the head and tumours in breasts. On the other hand, they are also quite important due to the representative difficulties in solving a large group of inverse problems concerning waves and fields, and thus researchers in mathematics, physics, and engineering societies have devoted great efforts to improving efficiencies and accuracies of the numerical solvers [5]. As shown in Figure 1, solving ISPs is to determine the unknown scatterers within a domain, when the domain is illuminated by several different incidences and we can measure the scattered fields outside the domain (usually contaminated by noise) for each incidence. It is well known that the main difficulties in such ISPs are their two intrinsic properties, i.e., ill-posedness and nonlinearity [1]. As most of the practical problems are three-dimensional (3-D) ones, which require much more computational resources than those in two-dimensional cases, they are often more difficult to solve due to the data deficiency (measurement aperture usually covers a small solid angle), and therefore also stronger ill-posedness and nonlinearity. Great efforts have been paid to tackle these demanding problems, for instance [6,7,8,9].

Methods for solving ISPs can be catalogued into two types of optimization approaches: deterministic and stochastic. The deterministic type of the inversion methods has been developed for decades, such as the contrast source inversion (CSI) method [6,10,11], the Born iterative method and distorted Born iterative method [12,13], the level set method [14,15], the subspace-based optimization method (SOM) [8,16,17,18], etc. The second type, the stochastic type of the inversion methods, usually employs a group of initial guesses and uses the stochastic optimization scheme to minimize the objective function, such as the genetic algorithm and the evolutionary optimization. The stochastic type methods increase the possibilities of finding the global minimum rather than being trapped in a local minimum as the deterministic optimization techniques [19,20]. Both techniques have also been applied to solve 3-D inverse scattering problems [6,21,22,23,24]. Other than the quantitative methods mentioned above, some qualitative methods, such as the linear sampling method [25] and other methods, like [26], are proposed to retrieve the geometric supports of the unknown scatterers.

To tackle the nonlinearity, the major difficulty in efficiently solving the ISPs, different types of integral equations have been proposed, such as in [27,28,29]. In particular, motivated by the contraction integral equation (CIE) in solving the direct scattering problems with highly conductive background media [30], a new type of contraction integral equation for inversion (CIE-I) has been proposed [29] to tackle the two-dimensional highly nonlinear ISPs by transforming the usually employed Lippmann–Schwinger integral equation (LSIE) into a new form with a modified medium contrast via a contraction mapping. With the CIE-I, the global multiple scattering effects (MSE) are substantially suppressed in estimating the modified contrast in the CSI type methods. With the FFT type twofold SOM regularization scheme [8], the inversion solver with CIE-I is capable of effectively alleviating the nonlinearity of the two-dimensional ISPs, especially those highly nonlinear ones with strong scatterers with large contrast and/or large dimensions (in terms of wavelengths). In this paper, this new CIE-I will be implemented in tackling the 3-D ISPs.

In summary, the contributions of the paper include:

• The CIE-I is firstly implemented to tackle the computationally costly full 3-D ISPs, to address the highly nonlinear 3-D ISPs, and to accelerate the convergence of the inversions.
• A relaxed type of inversion scheme based on CIE-I is proposed, with different auxiliary parameters $\beta$ (the parameter in CIE-I to control the portion of MSE in estimating the contrast) in updating the contrast sources and in updating the contrast. This means to further accelerate the convergence of the inversions.
• Several numerical tests are provided with details, for the sake of further algorithmic studies.

After the Introduction, the proposed 3-D inversion method is detailed in Section 2. In Section 3, three numerical examples are given to validate the proposed method. Finally, conclusions are drawn.

2. Inversion with CIE-I

In this paper, the domains of interest (DoI) are chosen to be rectangular cuboid in order to implement the conjugate gradient fast Fourier transform (CG-FFT) scheme and apply the Fourier basis in TSOM. For the convenience of reading, in this paper, we denote the one-dimensional tensor as $\overline{a}$, two-dimensional tensor as $\overline{\overline{a}}$, three-dimensional tensor as $\widehat{a}$, and four-dimensional tensor as $\widehat{\widehat{a}}$. Unless otherwise specified, the subscript of the tensors denotes the index of the element, such as ${\overline{\overline{a}}}_{m,n}$ denotes the element in $\overline{\overline{a}}$ with index $\{m,n\}$. We use bold symbols to denote vectorial physical quantities, such as the positions $\mathit{r}$ and the electric fields $\mathit{E}$ in 3-D cases.

2.1. 3-D Modeling

In 3-D cases, there are $N_{\mathrm{i}}$ incident waves from different angles impinging onto the rectangular cuboid DoI $\mathcal{D}$ ($\mathcal{D}\subset {\mathcal{R}}^3$, the background 3-D homogeneous medium with permittivity ${ϵ}_0$ and permeability ${\mu }_0$), where nonmagnetic scatterers are located, and these incident waves are expressed as ${\mathit{E}}_l^{\mathrm{inc}}(\mathit{r})$, $l=1,2,\cdots ,N_{\mathrm{i}}$, $\mathit{r}\in \mathcal{D}$. For each incidence, the scattered fields are collected by $N_{\mathrm{r}}$ receivers located at ${\mathit{r}}_j^{\prime }$, $j=1,2,\cdots ,N_{\mathrm{r}}$. With all such information, including every incident field inside the domain of interest and the corresponding scattered fields at the positions of all detectors, we aim at determining the dielectric profile $ϵ(\mathit{r})$, $\mathit{r}\in \mathcal{D}$.

The scattering models are governed by the following electric field volume integral equation based on the well-known Lippmann–Schwinger integral equation (LSIE). For the $l^{\mathrm{th}}$ incidence, the field equation in the domain $\mathcal{D}$ is expressed as

$${\mathit{I}}_l(\mathit{r})=\chi (\mathit{r}){\mathit{E}}_l^{\mathrm{inc}}(\mathit{r})+\chi (\mathit{r})({\mathit{G}}_{\mathrm{D}}^{3\mathrm{D}}{\mathit{I}}_l)(\mathit{r}),$$

(1)

where $\chi (\mathit{r})=\left({ϵ}_{\mathrm{r}}(\mathit{r})-1\right)$ is the contrast, ${ϵ}_{\mathrm{r}}(\mathit{r})$, ${\mathit{I}}_l(\mathit{r})=\chi (\mathit{r}){\mathit{E}}_l^{\mathrm{tot}}(\mathit{r})$ and ${\mathit{E}}_l^{\mathrm{inc}}(\mathit{r})$ are the relative permittivity, the contrast source and the incident electric field at $\mathit{r}$, respectively, and $({\mathit{G}}_{\mathrm{D}}^{3\mathrm{D}}{\mathit{I}}_l)(\mathit{r})$ is an integral operator with the dyadic Green’s function as the integral kernel, which can be written as [8,31],

$$({\mathit{G}}_{\mathrm{D}}^{3\mathrm{D}}{\mathit{I}}_l)(\mathit{r}):=k_0^2{\int }_{\mathcal{D}}\left(\overline{\overline{I}}+\frac{\nabla \nabla }{k_0^2}\right)·g(\mathit{r},{\mathit{r}}^{\prime }){\mathit{I}}_l({\mathit{r}}^{\prime })d{\mathit{r}}^{\prime },$$

(2)

with $k_0=\omega \sqrt{{ϵ}_0{\mu }_0}$ as the wave number of the background medium, and $g(\mathit{r},{\mathit{r}}^{\prime })$ as the 3-D Green’s function for the background homogeneous medium. Notice that the ${\mathit{I}}_l(\mathit{r})$ is the contrast source in this paper, while, in [8], it is the physically induced current that includes a multiplicative factor $-i\omega {ϵ}_0$ compared to the contrast source. Otherwise, as we know, the dyadic Green’s function is composed of nine scalar elements, which represent the mapping from the three components of the contrast source to the three components of the scattered fields inside the DoI.

The nonlinearity of ISPs comes from the MSE, as shown in Equation (1). To alleviate the nonlinearity of the model, in [29], by some mathematical manipulation on Equation (1), another new-type integral equation, which is denoted as CIE-I herein, can be obtained

$$\beta (\mathit{r}){\mathit{I}}_l(\mathit{r})=R(\mathit{r})\beta (\mathit{r}){\mathit{I}}_l(\mathit{r})+R(\mathit{r})\left[{\mathit{E}}_l^{\mathrm{inc}}(\mathit{r})+({\mathit{G}}_{\mathrm{D}}^{3\mathrm{D}}{\mathit{I}}_l)(\mathit{r})\right],$$

(3)

where $R(\mathit{r})=\beta (\mathit{r})\chi (\mathit{r})/\left[\beta (\mathit{r})\chi (\mathit{r})+1\right]$ is the modified contrast function, and $\beta (\mathit{r})$ is a chosen auxiliary parameter to control the portion of the MSE in estimating the contrast [29], which can be a constant or a variable at the different position in the DoI. For the convenience of discretizing the equation, we rewrite the dyadic Green’s function as

$$(G_{\mathrm{D};uv}^{3\mathrm{D}}{I}_{l;v})(\mathit{r}):=\left\{\begin{array}{cc}{k}_0^2(1+{\displaystyle \frac{1}{k_0^2}}{\displaystyle \frac{{\partial }^2}{\partial {\iota }_u^2}}){\int }_{\mathcal{D}}g(\mathit{r},{\mathit{r}}^{\prime })I_{l;v}({\mathit{r}}^{\prime })d{\mathit{r}}^{\prime },\hfill & \mathrm{if}u=v,\hfill \\ {\displaystyle \frac{{\partial }^2}{\partial {\iota }_u\partial {\iota }_v}}{\int }_{\mathcal{D}}g(\mathit{r},{\mathit{r}}^{\prime })I_{l;v}({\mathit{r}}^{\prime })d{\mathit{r}}^{\prime },\hfill & \mathrm{if}u\ne v,\hfill \end{array}\right.$$

(4)

where $u,v=1,2,3$ represent the x-, y- and z-components of a vector, respectively, and ${\iota }_1=x$, ${\iota }_2=y$, and ${\iota }_3=z$.

Following the conventions in [8], we can rewrite the vectorial Equation (3) into three coupled scalar equations in the discrete forms. Thus, we first discretize the rectangular domain of interest into small cuboid subdomains, whose dimensions are much smaller than the wavelength and the center of which are located at ${\mathit{r}}_{m,n,p}$, with m, n, p integers and $m\in [1,M_1]$, $n\in [1,M_2]$ and $p\in [1,M_3]$. Here, $M_1$, $M_2$, and $M_3$ are the total number of subdomains along x-, y-, and z-directions, respectively, and we let $M=M_1\times M_2\times M_3$ be the total number of the subdomains. With such discretization, we have

$${\widehat{\beta }}_{m,n,p}{\widehat{I}}_{l;u;m,n,p}={\widehat{R}}_{m,n,p}{\widehat{\beta }}_{m,n,p}{\widehat{I}}_{l;u;m,n,p}+{\widehat{R}}_{m,n,p}\left[{\widehat{E}}_{l;u;m,n,p}^{\mathrm{inc}}+\sum _{v=1}^3{\widehat{G}}_{\mathrm{D};uv;m,n,p}^{3\mathrm{D}}({\widehat{I}}_{l;v})\right],$$

(5)

where ${\widehat{R}}_{m,n,p}={\widehat{\beta }}_{m,n,p}{\widehat{\chi }}_{m,n,p}/\left[{\widehat{\beta }}_{m,n,p}{\widehat{\chi }}_{m,n,p}+1\right]$, ${\widehat{\chi }}_{m,n,p}$ is the contrast at ${\mathit{r}}_{m,n,p}$, whereas ${\widehat{E}}_{l;u;m,n,p}^{\mathrm{inc}}$ and ${\widehat{I}}_{l;u;m,n,p}$ are the incident electric field and the induced current at ${\mathit{r}}_{m,n,p}$, respectively. Subscript $u=1,2,3$ denotes the x, y, and z components of a vector. Note that the convolution-type operators ${\widehat{G}}_{\mathrm{D};uv}^{3\mathrm{D}}$ are obtained via the Equation (4). For further details of the discretization of Equation (1) and the finite difference scheme to generate $({\mathit{G}}_{\mathrm{D}}^{3\mathrm{D}}{\mathit{I}}_l)(\mathit{r})$, please refer to the Appendix of [6].

Similarly, the integral operator relating the contrast sources and the scattered fields could also be expressed as the summation of the contribution from all the subdomains,

$${\overline{E}}_l^{\mathrm{sca}}={\overline{\overline{G}}}_{\mathrm{S}}^{3\mathrm{D}}·{\overline{I}}_l,$$

(6)

where ${\overline{E}}_l^{\mathrm{sca}}={\left[{{\overline{E}}_{l;1}^{\mathrm{sca}}}^{\mathrm{T}},{{\overline{E}}_{l;2}^{\mathrm{sca}}}^{\mathrm{T}},{{\overline{E}}_{l;3}^{\mathrm{sca}}}^{\mathrm{T}}\right]}^{\mathrm{T}}$ is a $3N_{\mathrm{r}}$ dimensional vector with ${\overline{E}}_{l;\kappa }^{\mathrm{sca}}={\left[E_{l;\kappa ;1}^{\mathrm{sca}},E_{l;\kappa ;2}^{\mathrm{sca}},\cdots ,E_{l;\kappa ;N_{\mathrm{r}}}^{\mathrm{sca}}\right]}^{\mathrm{T}}$ ($\kappa =1,2,3$ denotes the x, y, z component of the corresponding vector, respectively), ${\overline{I}}_l$ is a $3M$ dimensional vector obtained by ${\overline{I}}_l=\mathrm{vec}\left\{{\widehat{\widehat{I}}}_l\right\}$. In a 3-D scenario, the vectorization operation $\mathrm{vec}\left\{·\right\}$ is defined to vectorize a four-dimensional tensor into a vector, i.e., if ${\overline{I}}_l=\mathrm{vec}\left\{{\widehat{\widehat{I}}}_l\right\}$, we have ${\overline{I}}_{l;\vartheta }={\widehat{\widehat{I}}}_{l;m,n,p,\kappa }$ with $\vartheta =(\kappa -1)\times M+(p-1)\times (M_1\times M_2)+(n-1)\times M_1+m$. In Equation (6), the scattering operator is defined as

$${\overline{\overline{G}}}_{\mathrm{S}}^{3\mathrm{D}}=\left[\begin{array}{ccc}{\overline{\overline{G}}}_{\mathrm{S};11}& {\overline{\overline{G}}}_{\mathrm{S};12}& {\overline{\overline{G}}}_{\mathrm{S};13}\\ {\overline{\overline{G}}}_{\mathrm{S};21}& {\overline{\overline{G}}}_{\mathrm{S};22}& {\overline{\overline{G}}}_{\mathrm{S};23}\\ {\overline{\overline{G}}}_{\mathrm{S};31}& {\overline{\overline{G}}}_{\mathrm{S};32}& {\overline{\overline{G}}}_{\mathrm{S};33}\end{array}\right],$$

(7)

a $3N_{\mathrm{r}}\times 3M$ matrix, with ${\overline{\overline{G}}}_{\mathrm{S};uv}$, a $N_{\mathrm{r}}\times M$ matrix, the mapping from the v component of the induced current to the u component of scattered fields (the subscripts $u,v=1,2,3$ are not indexed for tensor elements). The explicit expression of ${\overline{\overline{G}}}_{\mathrm{S};uv}$ is

$$\begin{array}{ll}{\overline{\overline{G}}}_{\mathrm{S};uv}^{3\mathrm{D}}(a,b)=& \left\{\left(k_0^2+\frac{i{k}_0}{R_{a,b}}-\frac{1}{R_{a,b}^2}\right)\delta (u-v)+\right.\hfill \\ & \left.\left[{({\mathit{r}}_a^{\prime })}_u-{({\mathit{r}}_{m,n,p})}_u\right]\left[{({\mathit{r}}_a^{\prime })}_v-{({\mathit{r}}_{m,n,p})}_v\right]\left(-\frac{k_0^2}{R_{a,b}^2}-\frac{3i{k}_0}{R_{a,b}^3}+\frac{3}{R_{a,b}^4}\right)\right\}g({\mathit{r}}_a^{\prime },{\mathit{r}}_{m,n,p}),\hfill \end{array}$$

(8)

where $\delta (y)$ is 1 when $y=0$ and is 0 otherwise, $R_{a,b}=\left|{\mathit{r}}_a^{\prime }-{\mathit{r}}_{m,n,p}\right|$ in which $b=(p-1)\times (M_1\times M_2)+(n-1)\times M_1+m$ with $a=1,2,\cdots ,N_{\mathrm{r}}$, $b=1,2,\cdots ,M$, $m=1,2,\cdots ,M_1$, $n=1,2,\cdots ,M_2$, and $p=1,2,\cdots ,M_3$. Here, ${(\mathit{r})}_u$ denotes the u component of $\mathit{r}$.

2.2. Objective Function for Inversions

In this subsection, we build the objective function used in the proposed inversion method, in which we will use the FFT-TSOM [8] as the regularization to stabilize the inversion, as done in [29]. Emphasize that the twofold subspace constraints have different regularization effects. The first fold, the original SOM, balances the two mismatches in the objective function, whereas the second fold, the TSOM, is the key to stabilize the inversions with CIE-I. Details are given in below.

For the original SOM part, as given in [16], with the spectral information of ${\overline{\overline{G}}}_{\mathrm{S}}^{3\mathrm{D}}$ (the singular value decomposition–SVD of ${\overline{\overline{G}}}_{\mathrm{S}}^{3\mathrm{D}}$ tells ${\overline{\overline{G}}}_{\mathrm{S}}^{3\mathrm{D}}·{\overline{v}}_j^{\mathrm{S}}={\sigma }_j^{\mathrm{S}}{\overline{u}}_j^{\mathrm{S}}$ and the complexity of thin SVD of ${\overline{\overline{G}}}_{\mathrm{S}}^{3\mathrm{D}}$ is $O(27N_{\mathrm{r}}^{2}M)$, assuming that the singular values ${\sigma }_j^{\mathrm{S}}$ is a non-increasing sequence), the contrast sources can be decomposed into two parts, deterministic part of the contrast sources (DPCS) and ambiguous part of the contrast sources (APCS), the former being obtained as

$${\overline{I}}_l^{\mathrm{d}}=\sum _{j=1}^L\frac{{\overline{u}}_j^{\mathrm{S}*}·{\overline{E}}_l^{\mathrm{sca}}}{{\sigma }_j^{\mathrm{S}}}{\overline{v}}_j^{\mathrm{S}}={\overline{\overline{V}}}_{\mathrm{S}}^{+}·{\overline{\alpha }}_l^{+},$$

(9)

where ${\overline{\overline{V}}}_{\mathrm{S}}^{+}=\left[{\overline{v}}_1^{\mathrm{S}},{\overline{v}}_2^{\mathrm{S}},\cdots ,{\overline{v}}_L^{\mathrm{S}}\right]$, ${\overline{\alpha }}_l^{+}={\left[{\alpha }_{l;1}^{+},{\alpha }_{l;2}^{+},\cdots ,{\alpha }_{l;L}^{+}\right]}^{\mathrm{T}}$ with ${\alpha }_{l;j}^{+}=({\overline{u}}_j^{\mathrm{S}*}·{\overline{E}}_l^{\mathrm{sca}})/{\sigma }_j^{\mathrm{S}}$, $j=1,2,\cdots ,L$, and the superscript ∗ denotes the Hermitian operation while superscript + refers to the dominant current subspace, the subspace corresponding to the dominant singular values. The value of L is chosen according to the noise level [16]. Later, we will see how this might work after introducing the objective function for the inversion.

For the second fold subspace constraint, according to the FFT-TSOM [8], the APCS can be written as

$${\overline{I}}_l^{\mathrm{a}}({\widehat{\widehat{\gamma }}}_l)={\overline{I}}_l^{\mathrm{tmp}}-{\overline{\overline{V}}}_{\mathrm{S}}^{+}·\left({\overline{\overline{V}}}_{\mathrm{S}}^{+*}·{\overline{I}}_l^{\mathrm{tmp}}\right),$$

(10)

where ${\widehat{\widehat{\gamma }}}_l=\left[{\widehat{\gamma }}_{x;l};{\widehat{\gamma }}_{y;l};{\widehat{\gamma }}_{z;l}\right]$, ${\overline{I}}_l^{\mathrm{tmp}}=\left[{\overline{I}}_{x;l}^{\mathrm{tmp}};{\overline{I}}_{y;l}^{\mathrm{tmp}};{\overline{I}}_{z;l}^{\mathrm{tmp}}\right]$ with

$${\overline{I}}_{u;l}^{\mathrm{tmp}}=\mathrm{vec}\left\{\mathrm{IDFT}\left\{{\widehat{\gamma }}_{u;l}\right\}\right\}$$

(11)

and $u=x$, y, and z, $\mathrm{IDFT}$ is the inverse discrete Fourier transform operator, the $\mathrm{vec}\{·\}$ is the vectorization operator. Note that the inverse discrete Fourier transform (IDFT) is performed by the 3-D FFT algorithm, the computational complexity of which is $O(M{\log }_{2}M)$, with $M=M_1\times M_2\times M_3$.

By using the low-frequency Fourier components, we are able to constrain the APCS within a low-dimensional subspace. The reason [17] is to use only the contrast sources components in this subspace that is influential to the scattered fields within the DoI, such that the stability of the inversions is substantially increased. To reduce the computational costs, such a subspace can be approximately spanned by low-frequency Fourier bases [8]. Here, if we use all Fourier bases, the construction of the APCS becomes the one in the original SOM. Otherwise, we can use the low-frequency components, and set the coefficients for those high-frequency components as zeros. This can be achieved via using a mask with eight corners being 1, with size $M_F$ ($1\le M_F\le M_1/2,M_2/2,M_3/2$), and other positions being 0. The details can be found in [8].

Having expressed the contrast sources as aforementioned, it is straightforward to define the objective function. Firstly, it is natural to give the mismatch of the scattered fields by

$${\Delta }_l^{\mathrm{fie}}({\widehat{\widehat{\gamma }}}_l)={∥{\overline{\overline{G}}}_{\mathrm{S}}·{\overline{I}}_l^{\mathrm{a}}+{\overline{\overline{G}}}_{\mathrm{S}}·{\overline{I}}_l^{\mathrm{d}}-{\overline{E}}_l^{\mathrm{sca}}∥}^2,$$

(12)

where ${\overline{I}}_l^{\mathrm{d}}$ and ${\overline{I}}_l^{\mathrm{a}}$ are as in Equations (9) and (10), respectively, and $∥·∥$ denotes the $L^2$ norm of a tensor. The current equation in Equation (5) is another key equation to satisfy. Using the APCS construction Equation (10), we define an operator as

$${\left({\mathcal{L}}^{3\mathrm{D}}\left({\widehat{\widehat{\gamma }}}_l\right)\right)}_{\kappa ;m,n,p}=\beta {\widehat{I}}_{l;\kappa ;m,n,p}^{\mathrm{a}}-\beta {\widehat{R}}_{m,n,p}{\widehat{I}}_{l;\kappa ;m,n,p}^{\mathrm{a}}-{\widehat{R}}_{m,n,p}\left[\sum _{v=1}^3{\widehat{G}}_{\mathrm{D};\kappa v;m,n,p}^{3\mathrm{D}}({\widehat{I}}_{l;v}^{\mathrm{a}})\right].$$

(13)

With this definition, we could write the mismatch of Equation (5) as

$${\Delta }_l^{\mathrm{cur}}({\widehat{\widehat{\gamma }}}_l,\widehat{R})={∥{\mathcal{L}}^{3\mathrm{D}}({\widehat{\widehat{\gamma }}}_l)-{\widehat{\widehat{\mathrm{\Gamma }}}}_l^{3\mathrm{D}}∥}^2,$$

(14)

where ${\widehat{\widehat{\mathrm{\Gamma }}}}_{l;\kappa ;m,n,p}^{3\mathrm{D}}=-\beta {\widehat{I}}_{l;\kappa ;m,n,p}^{\mathrm{d}}+\beta {\widehat{R}}_{m,n,p}{\widehat{I}}_{l;\kappa ;m,n,p}^{\mathrm{d}}+{\widehat{R}}_{m,n,p}\left[\left\{{\sum }_{v=1}^3{\widehat{G}}_{\mathrm{D};\kappa v;m,n,p}^{3\mathrm{D}}({\widehat{I}}_{l;v}^d)\right\}+{\widehat{\widehat{E}}}_{l;\kappa ;m,n,p}^{\mathrm{inc}}\right]$. Then, the objective function can be given as

$$f({\widehat{\widehat{\gamma }}}_1,{\widehat{\widehat{\gamma }}}_2,\cdots ,{\widehat{\widehat{\gamma }}}_{N_{\mathrm{i}}},\widehat{R})=\sum _{l=1}^{N_{\mathrm{i}}}\left({\Delta }_l^{\mathrm{fie}}/{∥{\overline{E}}_l^{\mathrm{sca}}∥}^2+{\Delta }_l^{\mathrm{cur}}/{∥{\overline{\overline{E}}}_l^{\mathrm{inc}}∥}^2\right).$$

(15)

The inversion is to minimize this objective function. As in [8], the conjugate gradient (CG) type algorithm that is used in contrast source inversion (CSI) method is adopted to minimize this nonlinear problem by alternatively updating the ${\widehat{\widehat{\gamma }}}_l$ and $\widehat{R}$ at every iteration of the optimization.

2.3. Sketch of the Inversion Method

Following the inversion method in [8,29], we summarize it as follows:

• Set the background medium and null APCS as the initial guesses and choose an L value such that the corresponding first L singular values of ${\overline{\overline{G}}}_{\mathrm{S}}^{3\mathrm{D}}$ are larger than the noise level (assuming that the noise is a white Gaussian one).
• Set proper values for $\beta$ in CIE-I modeling and proper value for $M_F$ to control the number of Fourier bases being used.
• Carry out the CG type optimization algorithm to alternatively update the two types of variables, where the APCS is updated with a one-step Polak–Ribière CG scheme and the contrast is updated with the least squares method.
• Stop the optimization if a termination condition is met, which can be a maximum number of iterations or a pre-defined relative change of APCS coefficients.
• If the maximum number of rounds of inversion is met, go to Step 6. Otherwise, the obtained contrast and APCS will be used as the initial guesses for the next round of optimization with smaller $\beta$ and larger $M_F$ in Step 3.
• Output the obtained contrast.

As mentioned in [29], the first round inversion with a large $\beta$ and with a low-dimensional subspace enables a fast convergence to a meaningful coarse result that could be used as an initial guess in the second round. By increasing the dimension of the APCS subspace and including more MSE in estimating the contrast, the inversion gives us a result with better resolution. Since in the second round a (supposedly) good initial guess is given, the convergence is also very fast. For some difficult problems, the inversion could be carried out with multi-round optimizations by gradually increasing the dimensions of the APCS subspace in each round while decreasing the values of $\beta$ in estimating the contrast.

2.4. Updating Contrast and Contrast Sources with Different $\beta$

As mentioned above, there are two updates at each iteration for the two types of unknowns. In [29], the same $\beta$ is used in each round of optimization for both updates. In the first round optimization in inversions, to tackle the nonlinearity, a large $\beta$ is needed in updating the modified contrast, and therefore we use the same $\beta$ for the update of APCS. In the subsequent rounds (instead of the first round) of optimization, for the sake of stability, a small $\beta$ is used in estimating the modified contrast to cope with a high-dimensional APCS subspace, as discussed in [29]. However, there is not such a need to use a small $\beta$ in updating the APCS. Consequently, in the second round of optimization, we can use a CIE-I model with larger value ${\beta }_1$ for the update of the APCS and another with smaller value ${\beta }_2$ for the update of the contrast at the same iteration. The same applies to the third round optimization, if there is such a need. The purpose of doing so is to further accelerate the convergence of these computationally burdensome 3-D inversions after the first round optimization.

As for the value range for ${\beta }_2$, we carry out a similar calculation as in [29], which reflects the norm of the scattering operator that maps the contrast sources to the scattered fields within DoI. The results shown in Figure 2 indicate that, if we want to suppress the MSE in estimating the contrast when using the least squares method, we might need to choose a value for ${\beta }_2$ that is larger than 3.5, and this is the guideline for the first round optimization. For the second or subsequent rounds, as good initial guesses are provided, we can include more MSE to retrieve the fine features of the unknown scatterers.

3. Numerical Simulations

In this section, we will test the proposed inversion methods in three examples. We will use the same physical setup of sources and receivers for the 3-D example in [8], as shown in Figure 3. The DoI are all the same, a box with size $3\lambda \times 3\lambda \times 3\lambda$. The box is illuminated by 60 electric dipole antennas at 300 MHz ($\lambda =1$ m in air), located at three circles (with 20 dipole antennas evenly distributed on each) with the same radius 3 m. The three circles are in $x-y$, $y-z$ and $x-z$ planes, and their centers are at (0.2, 0, –0.1), (0.1, 0, –0.15), and (–0.05, 0.1, 0), respectively. The direction of the electric dipole sources in the $y-z$ plane are in the x-direction, while those in the $x-z$ and $x-y$ planes are in the y- and z-directions, respectively. Scattered fields are measured by 60 receivers, located at the same positions as the 60 dipole sources. We collect all three components of the fields. Consequently, we have a $60\times 180$ synthetic data matrix. In all three tests, the synthetic data are calculated with a $60\times 60\times 60$ mesh, while, in the inversions, a $30\times 30\times 30$ mesh is used. All synthetic data are contaminated with additive Gaussian white noise with level of $10\%$.

The termination condition in each round of optimization is to reach a pre-defined relative changing rate of the APCS coefficients, given as

$${\delta }^{3\mathrm{D}}=\sqrt{\frac{1}{N_{\mathrm{i}}}{\sum }_{l=1}^{N_{\mathrm{i}}}\frac{\parallel {\widehat{\widehat{\gamma }}}_l^{(h)}-{\widehat{\widehat{\gamma }}}_l^{(h-1)}{\parallel }^2}{\parallel {\widehat{\widehat{\gamma }}}_l^{(h-1)}{\parallel }^2}},$$

(16)

where ${\widehat{\widehat{\gamma }}}_l^{(h)}$ is the APCS Fourier coefficient for the $l^{\mathrm{th}}$ incidence at $h^{\mathrm{th}}$ iteration. We set ${\delta }^{3\mathrm{D}}<{10}^{-3}$ as the termination condition.

The reconstruction results are quantitatively evaluated by the following error:

$$Err=\sqrt{\frac{1}{M}{\sum }_{m,n,p}\frac{|{\widehat{ϵ}}_{\mathrm{r};m,n,p}^{\mathrm{est}}-{\widehat{ϵ}}_{\mathrm{r};m,n,p}^{\mathrm{tr}}{|}^2}{|{\widehat{ϵ}}_{\mathrm{r};m,n,p}^{\mathrm{tr}}{|}^2}}.$$

(17)

When running on a workstation with eight threads and 32 Gb RAM, every iteration of the proposed inversion method costs about 25 seconds CPU time in MATLAB, including one update on the contrast sources and one update on the contrast. Due to the independence between different incidences, we use 8 MATLAB workers to update the contrast sources in parallel.

3.1. Example 1

We reuse the 3-D test in [8], where a coated cube is employed. The cube is with an outer layer (${ϵ}_{\mathrm{r}2}=1.5+i0.3$) and an inner layer (${ϵ}_{\mathrm{r}1}=2+i0.8$), where the outer edge length is $b=2\lambda$ and inner edge length is $a=1\lambda$, as shown in Figure 4a. In [8], we observe that the inversion with LSIE and FFT-TSOM is able to satisfactorily retrieve the cube while the LSIE with SOM fails. Here, we illustrate how the proposed inversion method with CIE-I and FFT-TSOM performs. Similarly, we carry out two rounds of optimization for this inversion as done in [8] with $M_F=6$ and 10, but with CIE-I, we use $\beta =6$ in the first round and $\beta =1$ for the second round. The reconstruction results are shown in Figure 5, with $\Re \left\{{\widehat{ϵ}}_{\mathrm{r}}^{\mathrm{est}}\right\}$ as the real part of the reconstructed relative permittivity and $\Im \left\{{\widehat{ϵ}}_{\mathrm{r}}^{\mathrm{est}}\right\}$ as the imaginary part, where we see after two rounds of optimization that the coated cube is successfully found.

Now, we compare the convergence speeds of inversions with LSIE and CIE-I. We plot the reconstruction errors for both cases, as shown in Figure 6, in which the inversion with CIE-I appears to be converging much faster than the one with LSIE, with the same reconstruction quality. This confirms again the speeding convergence of inversions in two-dimensional cases shown in [29]. We would like to emphasize that this fast convergence property is extremely important in handling the 3-D ISPs, since the computational costs in 3-D problems are much higher than in two-dimensional cases.

With this example, we further test the relaxed scheme with different $\beta$ for updating the APCS and the contrast. The termination condition remains the same, ${\delta }^{3\mathrm{D}}<{10}^{-3}$. We apply this scheme in the second round optimization, with ${\beta }_1=6$ for the APCS update and ${\beta }_2=1$ for the contrast update. For comparison, the recorded reconstruction error at each iteration is plotted in Figure 6 as well, which shows that the convergence is further accelerated compared with the case of using ${\beta }_1=1$.

3.2. Example 2

In this test, we use a profile with four cubes with increasing contrasts. Having edge length $0.8\lambda$, they are centered at $(-0.7,0.7,-0.7)\lambda$, $(0.7,-0.7,-0.7)\lambda$, $(-0.7,-0.7,0.7)\lambda$, and $(0.7,0.7,0.7)\lambda$, with ${ϵ}_{\mathrm{r}}=2,3,4$ and 5, respectively, as shown in Figure 7.

Two inversions are carried out with LSIE and CIE-I, each with two rounds ($M_F=6$ and 10). For CIE-I, the choices of $\beta$ are 6 and 1 in the two rounds. The reconstruction results and errors are shown in Figure 8, Figure 9 and Figure 10. From these results, we clearly observe that the inversion with LSIE is only able to find the bottom two cubes with lower contrasts, whereas the inversion with CIE-I succeeds in finding all four of four cubes. Reconstruction errors displayed in Figure 9 confirm it, inversion with LSIE diverging while converging with CIE-I.

In addition, we test the inversion scheme with different $\beta$ for updating the APCS and for updating the contrast. Firstly, we choose ${\beta }_1=6$ for the update of APCS and ${\beta }_2=1$ for the update of the contrast in the second round of optimization with $M_F=10$, which however immaturely converges. This might be due to the reason of using a large ${\beta }_1$. Then, we choose ${\beta }_1=3$, and the optimization converges to a good solution with a convergence speed faster than when using ${\beta }_1=1$.

From this example, the inversion with CIE-I is able to handle the highly nonlinear problems consisting of strong scatterers, those with either large contrasts or large dimensions (in terms of wavelength), while the inversion with LSIE may fail.

3.3. Example 3

We will use another highly nonlinear example to confirm the large difference between the performances of inversion methods with CIE-I and LSIE. In this example, we use a profile that is similar to the famous two-dimensional “Austria” that has an annular and two separated disks. Here, we use a coated cube with a hollow inside and two rods outside the cube, as shown in Figure 11, all of which are with ${ϵ}_{\mathrm{r}}=2.5$.

We still carry out inversions with LSIE and CIE-I with two rounds of optimization, the first with $M_F=6$ and the second with $M_F=10$. For the CIE-I, ${\beta }_1={\beta }_2=6$ is used in the first round, and ${\beta }_1={\beta }_2=1$ in the second round. The reconstructions and errors are shown in Figure 12, Figure 13a and Figure 14. We see that the inversion with LSIE fails to find the profile while the one with CIE-I succeeds. This can also be seen from the errors of the reconstructions in Figure 13a, where the error by the inversion with LSIE diverges.

In the reconstruction results in Figure 14, we see that there is an artefact within the coated cube, where the medium is supposed to be air. This artefact appears in the reconstructed result after the first round optimization with $M_F=6$, meaning that the inversion, though capable of retrieving most of the main features of the scatterers, is still not stable enough. Therefore, we carry out another inversion with three rounds of optimization, with $M_F=4$, 8, and 10 in each round, and $\beta =6$, 1, and $0.5$, respectively. The final reconstruction result is given in Figure 15, and we see that there is no more such an artefact.

At last, we also carry out the inversions with different $\beta$ values for the updates of the APCS and the contrast, and the results are shown in Figure 13. We see that with larger ${\beta }_1$ value for the update of the APCS in the second round optimization, the convergence speed is indeed faster.

From this example, we again confirm the great resolvability of the inversion method with CIE-I against the high nonlinearity of the ISPs.

4. Conclusions

In this paper, we propose an inversion method with CIE-I to tackle the 3-D electromagnetic ISPs, especially for highly nonlinear problems with strong scatterers, those with high contrasts or electrically large dimensions. With the CIE-I modeling, the multiple scattering effects in estimating the contrast can be suppressed such that the nonlinearity of inversion can be effectively alleviated. Together with the FFT-TSOM regularization scheme, this largely increases the ability of retrieving the profile of strong scatterers, and effectively accelerates the convergence when handling the moderate scatterers. In addition, by relaxing the contraction mapping, i.e., different contraction mappings are used in updating the contrast sources and in updating the contrast, the convergence of the inversions can be further accelerated. Through numerical examples, we clearly show that the inversions with CIE-I outperform the ones with LSIE, in terms of the resolvability against the nonlinearity and the convergence speed. As shown in [32], the resolvability of the inversion solver with CIE-I against nonlinearity can be further improved in the two-dimensional case if proper regularization techniques are implemented, which should be expected in the 3-D case as well.