Inversion of Multi-frequency Data with the Cross-Correlated Contrast Source Inversion Method

TL;DR

This paper introduces a multi-frequency extension of the cross-correlated contrast source inversion (CC-CSI) method for electromagnetic inverse scattering, demonstrating improved performance over existing methods through numerical and experimental results.

Contribution

The paper proposes a novel multi-frequency CC-CSI method that enhances inversion accuracy for complex scatterers using electromagnetic data.

Findings

Multi-frequency CC-CSI outperforms MR-CSI in complex scatterer inversion.

Numerical and experimental results validate the effectiveness of the proposed method.

The method is applicable to both TM and TE polarizations.

Abstract

Cross-correlated contrast source inversion (CC-CSI) is a non-linear iterative inversion method that is proposed recently for solving the inverse scattering problems. In CC-CSI, a cross-correlated error is constructed and introduced to the cost functional, which improves the inversion ability when compared to the classical design of the cost functional by exploiting the mismatch between the data error and state error. In this paper, the multi-frequency inversion for electromagnetic waves is considered and a multi-frequency version of CC-CSI is proposed. Numerical and experimental inversion results of both transverse magnetic (TM) and transverse electric (TE) polarization demonstrate that, when multi-frequency data are available, CC-CSI still outperforms the multiplicative-regularized CSI method (MR-CSI) in the inversion of more complicated scatterers.

Figures

Equations

$${\mathbf{y}}_p$$

$${\mathbf{j}}_p$$

$${\mathbf{\rho }}_{p,i,\ell -1/2}$$

$${\mathbf{\gamma }}_{p,i,\ell -1/2}$$

$${\mathbf{\xi }}_{p,i,\ell -1/2}={\mathbf{y}}_{p,i}-{\mathbf{\Phi }}_{p,i}\left({\mathbf{\chi }}_{i,\ell -1}{\mathbf{e}}_{p,i}^{\text{inc}}+{\mathbf{\chi }}_{i,\ell -1}{\mathbf{A}}_i^{-1}{\mathbf{j}}_{p,i,\ell -1}\right),$$

$$\begin{array}{rl}{\mathcal{C}}_{\text{MF-CC-CSI},\ell -1/2}=& \sum _{i=1}^I{\eta }_i^{\mathcal{S}}\sum _{p=1}^P{\Vert {\mathbf{\rho }}_{p,i,\ell -1/2}\Vert }_{\mathcal{S}}^2+\sum _{i=1}^I{\eta }_{i,\ell -1}^{\mathcal{D}}\sum _{p=1}^P{\Vert {\mathbf{\gamma }}_{p,i,\ell -1/2}\Vert }_{\mathcal{D}}^2+\\ & \sum _{i=1}^I{\eta }_i^{\mathcal{S}}\sum _{p=1}^P{\Vert {\mathbf{\xi }}_{p,i,\ell -1/2}\Vert }_{\mathcal{S}}^2,\end{array}$$

$${\eta }_i^{\mathcal{S}}=\sum _{p=1}^P{\Vert {\mathbf{y}}_{p,i}\Vert }_{\mathcal{S}}^2,\text{and}{\eta }_{i,\ell -1}^{\mathcal{D}}=\sum _{p=1}^P{\Vert {\mathbf{\chi }}_{i,\ell -1}{\mathbf{e}}_{p,i}^{\text{inc}}\Vert }_{\mathcal{D}}^2,$$

$$\begin{array}{rl}{\mathbf{g}}_{p,i}=& -2{\eta }_i^{\mathcal{S}}{\mathbf{\Phi }}_{p,i}^H{\mathbf{\rho }}_{p,i,\ell -1}+2{\eta }_{i,\ell -1}^{\mathcal{D}}{\left({\mathbf{\chi }}_{i,\ell -1}{\mathbf{A}}_i^{-1}-\mathbf{I}\right)}^H{\mathbf{\gamma }}_{p,i,\ell -1}-\\ & 2{\eta }_i^{\mathcal{S}}{\left({\mathbf{\Phi }}_{p,i}{\mathbf{\chi }}_{i,\ell -1}{\mathbf{A}}_i^{-1}\right)}^H{\mathbf{\xi }}_{p,i,\ell -1}.\end{array}$$

$${\mathbf{j}}_{p,i,\ell }={\mathbf{j}}_{p,i,\ell -1}+{\alpha }_{p,i,\ell }{\mathbf{\nu }}_{p,i,\ell },$$

$${\mathbf{\nu }}_{p,i,\ell }=\{\begin{array}{ll}0,& \ell =0,\\ {\mathbf{g}}_{p,i,\ell }+\frac{{\sum }_{p^{\prime }}{⟨{\mathbf{g}}_{p^{\prime },i,\ell },{\mathbf{g}}_{p^{\prime },i,\ell }-{\mathbf{g}}_{p^{\prime },i,\ell -1}⟩}_{\mathcal{D}}}{{\sum }_{p^{\prime }}{\Vert {\mathbf{g}}_{p^{\prime },i,\ell -1}\Vert }_{\mathcal{D}}^2}{\mathbf{\nu }}_{p,i,\ell -1},& \ell \ge 1.\end{array}$$

$${\mathbf{g}}_{p,i,\ell }={{\mathbf{g}}_{p,i}|}_{{\mathbf{j}}_{p,i}={\mathbf{j}}_{p,i,\ell -1}}.$$

$${{\mathcal{C}}_{\text{MF-CC-CSI},\ell -1/2}|}_{{\mathbf{j}}_{p,i}={\mathbf{j}}_{p,i,\ell -1}+{\alpha }_{p,i}{\mathbf{\nu }}_{p,i,\ell }}.$$

$${\mathbf{e}}_{p,i,\ell }^{\text{tot}}={\mathbf{e}}_{p,i,\ell -1}^{\text{tot}}+{\alpha }_{p,i,\ell }{\mathbf{A}}_i^{-1}{\mathbf{\nu }}_{p,i,\ell }.$$

$${\mathbf{\chi }}_{1,\ell -1}=\Re \left\{{\mathbf{\chi }}_{i,\ell -1}\right\}+\frac{{\omega }_1}{{\omega }_i}\Im \left\{{\mathbf{\chi }}_{i,\ell -1}\right\},i=1,2,3,\cdots ,I,$$

$${\mathbf{\chi }}_{\ell }={\mathbf{\chi }}_{\ell -1}+{\beta }_{\ell }{\mathbf{\nu }}_{\ell }^{\chi },$$

$${\mathbf{\nu }}_{\ell }^{\chi }=\{\begin{array}{ll}0,& n=0,\\ {\mathbf{g}}_{\ell }^{\chi }+\frac{{⟨{\mathbf{g}}_{\ell }^{\chi },{\mathbf{g}}_{\ell }^{\chi }-{\mathbf{g}}_{\ell -1}^{\chi }⟩}_{\mathcal{D}}}{{\Vert {\mathbf{g}}_{\ell -1}^{\chi }\Vert }_{\mathcal{D}}^2}{\mathbf{\nu }}_{\ell -1}^{\chi },& n\ge 1,\end{array}$$

$${\mathcal{C}}_{\text{MF-CC-CSI},\ell }=\sum _{i=1}^I{\eta }_{i,\ell -1}^{\mathcal{D}}\sum _{p=1}^P{\Vert {\mathbf{\gamma }}_{p,i,\ell }\Vert }_{\mathcal{D}}^2+\sum _{i=1}^I{\eta }_i^{\mathcal{S}}\sum _{p=1}^P{\Vert {\mathbf{\xi }}_{p,i,\ell }\Vert }_{\mathcal{S}}^2.$$

$${\mathbf{\gamma }}_{p,i,\ell }={\mathbf{\chi }}_{i,\ell -1}{\mathbf{e}}_{p,i}^{\text{inc}}+{\mathbf{\chi }}_{i,\ell -1}{\mathbf{A}}_i^{-1}{\mathbf{j}}_{p,i,\ell }-{\mathbf{j}}_{p,i,\ell },$$

$${\mathbf{\xi }}_{p,i,\ell }={\mathbf{y}}_{p,i}-{\mathbf{\Phi }}_{p,i}\left({\mathbf{\chi }}_{i,\ell -1}{\mathbf{e}}_{p,i}^{\text{inc}}+{\mathbf{\chi }}_{i,\ell -1}{\mathbf{A}}_i^{-1}{\mathbf{j}}_{p,i,\ell }\right).$$

$${\mathbf{g}}_{\ell }^{\chi }=\frac{2\Re \left\{{\sum }_{i=1}^I{\mathbf{g}}_{i,\ell }^{\chi }\right\}}{{\sum }_{i=1}^I{\sum }_{p=1}^P{\mathbf{e}}_{p,i,\ell }^{\text{tot}}\overline{{\mathbf{e}}_{p,i,\ell }^{\text{tot}}}}+\text{i}\frac{2\Im \left\{{\sum }_{i=1}^I\frac{{\omega }_1}{{\omega }_i}{\mathbf{g}}_{i,\ell }^{\chi }\right\}}{{\sum }_{i=1}^I{\left(\frac{{\omega }_1}{{\omega }_i}\right)}^2{\sum }_{p=1}^P{\mathbf{e}}_{p,i,\ell }^{\text{tot}}\overline{{\mathbf{e}}_{p,i,\ell }^{\text{tot}}}},$$

$${\mathbf{g}}_{i,\ell }^{\chi }={\eta }_{i,\ell -1}^{\mathcal{D}}\sum _{p=1}^P\overline{{\mathbf{e}}_{p,i,\ell }^{\text{tot}}}{\mathbf{\gamma }}_{p,i,\ell }-{\eta }_i^{\mathcal{S}}\sum _{p=1}^P\overline{{\mathbf{e}}_{p,i,\ell }^{\text{tot}}}{\mathbf{\Phi }}_{p,i}^H{\mathbf{\xi }}_{p,i,\ell }.$$

$$\sum _{i=1}^I\frac{{\sum }_{p=1}^P{\Vert \left({\mathbf{\chi }}_{i,\ell -1}+\beta {\mathbf{\chi }}_{i,\ell }\right){\mathbf{e}}_{p,i,\ell }^{\text{tot}}-{\mathbf{j}}_{p,i,\ell }\Vert }_{\mathcal{D}}^2}{{\sum }_{p=1}^P{\Vert \left({\mathbf{\chi }}_{i,\ell -1}+\beta {\mathbf{\chi }}_{i,\ell }\right){\mathbf{e}}_{p,i}^{\text{inc}}\Vert }_{\mathcal{D}}^2}+\sum _{i=1}^I{\eta }_i^{\mathcal{S}}\sum _{p=1}^P{\Vert {\mathbf{y}}_{p,i}-{\mathbf{\Phi }}_{p,i}\left({\mathbf{\chi }}_{i,\ell -1}+\beta {\mathbf{\chi }}_{i,\ell }\right){\mathbf{e}}_{p,i,\ell }^{\text{tot}}\Vert }_{\mathcal{S}}^2.$$

$${\mathbf{j}}_{p,i,0}=\frac{{\Vert {\mathbf{\Phi }}_{p,i}^H{\mathbf{y}}_{p,i}\Vert }_{\mathcal{D}}^2}{{\Vert {\mathbf{\Phi }}_{p,i}{\mathbf{\Phi }}_{p,i}^H{\mathbf{y}}_{p,i}\Vert }_{\mathcal{S}}^2}{\mathbf{\Phi }}_{p,i}^H{\mathbf{y}}_{p,i},$$

$${\mathbf{e}}_{p,i,0}^{\text{tot}}={\mathbf{e}}_{p,i}^{\text{inc}}+{\mathbf{A}}_i^{-1}{\mathbf{j}}_{p,i,0}.$$

$${\mathbf{\chi }}_{1,0}=\frac{\Re \left\{\sum _{i=1}^I\sum _{p=1}^P{\mathbf{j}}_{p,i,0}\overline{{\mathbf{e}}_{p,i,0}^{\text{tot}}}\right\}}{\sum _{i=1}^I\sum _{p=1}^P{\mathbf{e}}_{p,i,0}^{\text{tot}}\overline{{\mathbf{e}}_{p,i,0}^{\text{tot}}}}+\text{i}\frac{\Im \left\{\sum _{i=1}^I\frac{{\omega }_1}{{\omega }_i}\sum _{p=1}^P{\mathbf{j}}_{p,i,0}\overline{{\mathbf{e}}_{p,i,0}^{\text{tot}}}\right\}}{\sum _{i=1}^I{\left(\frac{{\omega }_1}{{\omega }_i}\right)}^2\sum _{p=1}^P{\mathbf{e}}_{p,i,0}^{\text{tot}}\overline{{\mathbf{e}}_{p,i,0}^{\text{tot}}}}.$$

$$err:={\Vert \hat{\mathbf{\chi }}-\mathbf{\chi }\Vert }_2/{\Vert \mathbf{\chi }\Vert }_2,$$

$$\mathbf{x}:=\left\{{\mathbf{\chi }}_i,{\mathbf{e}}_{p,i}^{\text{tot}},p=1,2,3,\cdots ,P,i=1,2,3,\cdots ,I\right\}.$$

$$\mathbf{x}\left({\beta }_1,{\beta }_2\right):={\beta }_2\left[\left({\beta }_1+1\right){\mathbf{x}}_{\text{MF-CC-CSI}}^{\text{opt}}-{\beta }_1{\mathbf{x}}^{\text{act}}\right]-\left({\beta }_2-1\right){\beta }_1{\mathbf{x}}_{\text{MF-MR-CSI}}^{\text{opt}},$$

$$\begin{array}{rl}& {{\mathcal{C}}_{\text{MF-CC-CSI},\ell -1/2}|}_{{\mathbf{j}}_{p,i}={\mathbf{j}}_{p,i,\ell -1}+{\alpha }_{p,i}{\mathbf{\nu }}_{p,i,\ell }}=\sum _{i=1}^I{\eta }_i^{\mathcal{S}}\sum _{p=1}^P{\Vert {\mathbf{\rho }}_{p,i,\ell -1/2}-{\alpha }_{p,i}{\mathbf{\Phi }}_i{\mathbf{\nu }}_{p,i,\ell }\Vert }^2+\\ & \sum _{i=1}^I{\eta }_{i,\ell -1}^{\mathcal{D}}\sum _{p=1}^P{\Vert {\mathbf{\gamma }}_{p,i,\ell -1/2}+{\alpha }_{p,i}\left({\mathbf{\chi }}_{\ell -1}{\mathcal{M}}_{\mathcal{D}}{\mathbf{A}}_i^{-1}-\mathbf{I}\right){\mathbf{\nu }}_{p,i,\ell }\Vert }^2+\\ & {\sum _{i=1}^I{\eta }_i^{\mathcal{S}}\sum _{p=1}^P{\Vert {\mathbf{\xi }}_{p,i,\ell -1/2}-{\alpha }_{p,i}{\mathbf{\Phi }}_i{\mathbf{\chi }}_{\ell -1}{\mathcal{M}}_{\mathcal{D}}{\mathbf{A}}_i^{-1}{\mathbf{\nu }}_{p,i,\ell }\Vert }^2|}_{{\mathbf{j}}_{p,i}={\mathbf{j}}_{p,i,\ell -1}+{\alpha }_{p,i}{\mathbf{\nu }}_{p,i,\ell }}.\end{array}$$

$${{\mathcal{C}}_{\text{MF-CC-CSI},\ell -1/2}|}_{{\mathbf{j}}_{p,i}={\mathbf{j}}_{p,i,\ell -1}+{\alpha }_{p,i}{\mathbf{\nu }}_{p,i,\ell }}=\sum _{j=0}^2\sum _{i=1}^I\left(a_{p,i,j}+b_{p,i,j}+c_{p,i,j}\right){\alpha }_{p,i}^j.$$

$$\begin{array}{rl}{\alpha }_{p,i,\ell }& =\underset{{\alpha }_{p,i}}{\text{arg}\max }\left\{{{\mathcal{C}}_{\text{MF-CC-CSI},\ell -1/2}|}_{{\mathbf{j}}_{p,i}={\mathbf{j}}_{p,i,\ell -1}+{\alpha }_{p,i}{\mathbf{\nu }}_{p,i,\ell }}\right\}\\ & =-\frac{1}{2}\frac{a_{p,i,1}+b_{p,i,1}+c_{p,i,1}}{a_{p,i,2}+b_{p,i,2}+c_{p,i,2}}.\end{array}$$

Full text

1 Introduction

Determining the inhomogeneity in a certain region of medium by probing the scattered fields (electromagnetic and acoustic) is a common problem arising from many different fields in science and engineering such as remote sensing, biomedical imaging, geophysical exploration, nondestructive testing, etc. The information of interest could be the morphological information or the values of the inhomogeneities. In this paper, efforts have been made to study the inversion of the dielectric parameters, which is referred to as the inverse medium problem for electromagnetic waves [1]. The research output is of course applicable for achieving the support of the objects, and it also applies to the inversion for acoustic waves.

Since the inverse scattering problem we will be discussing is in the resonance range, i.e., the wavelength is comparable to the dimension of the objects, it turns out to be inherently nonlinear. Severe ill-posedness also accompanies this inverse scattering problem, which comes from both the inherent challenges of the physical problem itself and the incompleteness and/or inaccuracy of the measurement data domain. There are at least two main macro-classes of fully-nonlinear iterative inversion approaches: stochastic methods and deterministic ones. In cases where the dimension of the solution space is not huge, stochastic global optimization methods [2, 3, 4, 5] are good candidates to search for the global optimal solution. Among the deterministic approaches, one is to use the linear sampling method together with a knowledge of the first transmission eigenvalue [6] or with several so-called “virtual experiments” [7]. A more accurate alternative is to recover the dielectric parameters and the total fields iteratively using an optimization method [8, 9, 10, 11, 12]. The iterative algorithms have been further improved over the recent decades by using regularization constraint [13, 14, 15], multi-scaling technique [16], wavelet transformation [17], etc. As an example, sparsity constraints can be appropriately exploited and compressive sensing-based techniques are proved to be effective in both qualitative imaging [18, 19, 20] and quantitative ones [21, 22, 23, 15]. However, it is worth noting that the cost functional remains unchanged in the aforementioned research work, which consists of two error terms: the data error and the state error. This classical construction of the cost functional has been changed recently by the work of [24], in which a cross-correlated error was proposed as a measure of the mismatch between the two error terms. Subsequently, the cross-correlated error was introduced into the cost function, which leads to a novel inversion method, referred to as the cross-correlated contrast source inversion (CC-CSI) method. It has been demonstrated that CC-CSI shows better inversion performance than the classical contrast source inversion (CSI) method [11] and multiplicative regularized CSI (MR-CSI) method [13], especially for inverting more complicated objects of higher contrast values. However, CC-CSI was only tested with single frequency in the work of [24], and it is still open whether it shows superiority when multi-frequency data is available, because the reliability of MR-CSI can be improved indeed by exploitation of multiple frequencies [25].

In this paper, we consider the inversion of multi-frequency electromagnetic data with CC-CSI. A multi-frequency version of CC-CSI (MF-CC-CSI) is proposed, which is able to process the multi-frequency data simultaneously. To validate the advantage of the proposed method, comparison has been made to the multi-frequency version of MR-CSI (MF-MR-CSI) proposed in [25]. Numerical results show that, compared to MF-MR-CSI, a more stabilized solution can be obtained by MF-CC-CSI. The $2$-nd Fresnel datasets in the year of 2005 [26] have been selected as the experimental data for further validation. In the remainder of this paper, both transverse magnetic (TM) and transverse electric (TE) polarizations are considered. Problem statement and formulation of the proposed MF-CC-CSI method are introduced in Section 2, the numerical simulation is presented in Section 3, and the experimental data inversion is given in Section 4. The main body of this paper is finalized with conclusion in Section 5.

2 Problem Statement and Formulation

2.1 Problem Statement

The multi-frequency inversion problem discussed in this paper is assumed to be non-dispersive, i.e., the contrast to be inverted is independent of frequency. Let us consider the canonical 2-D inverse scattering problem in a known background $\mathcal{D}$. The measurement domain $\mathcal{S}$ contains the sources denoted by the subscript $p\in\{1,2,3...,P\}$ and the receivers denoted by the subscript $q\in\{1,2,3,...,Q\}$. Equal subscript means the same position. We use a right-handed coordinate system in which the unit vector in the invariant direction points out of the paper. The time factor of $\exp(\text{i}\omega t)$ is considered in this paper, where $\text{i}^{2}=-1$.

Now the data equation and the state equation can be formulated based on finite difference frequency domain (FDFD) scheme as (see [24])

[TABLE]

respectively, with $p=1,2,\dots,P$, where, $\bm{y}_{p}$ represents the measurement data; $\bm{\Phi}_{p}:=\mathcal{M}_{\mathcal{S},p}\bm{A}^{-1}$ is the measurement matrix; $\mathcal{M}_{\mathcal{S},p}$ is an operator that interpolates field values defined at the finite-difference grid points to the appropriate receiver positions; $\bm{A}$ is the FDFD stiffness matrix; $\bm{e}^{\text{inc}}_{p}$ and $\bm{e}_{p}$ are the incident and total electric fields in the form of a column vector; $\bm{\chi}$ is the complex contrast consisting of the contrast permittivity, $\Delta\bm{\varepsilon}$, and the contrast conductivity, $\Delta\bm{\sigma}$, i.e., $\bm{\chi}:=\Delta\bm{\varepsilon}-\text{i}\Delta\bm{\sigma}/\omega$; $\bm{j}_{p}:=\bm{\chi}\bm{e}_{p}$ is defined as the component-wise multiplication of the contrast and the total fields, which is referred to as the contrast sources. In the remainder of this paper, $\omega^{2}$ is incorporated into $\bm{\Phi}_{p}$ and $\bm{A}$ for the sake of conciseness. Both of the matrices contain the background information. The inverse problem is to reconstruct the contrast $\bm{\chi}$ from the incomplete and/or inaccurate measurement data, $\bm{y}_{p}$.

2.2 Multi-frequency CC-CSI

2.2.1 Modified Cost Functional, $\mathcal{C}_{\text{MF-CC-CSI},\ell-1/2}$

Consider the inversion of multi-frequency data, the subscript $i$ is used to represent the $i$-th frequency. The data error equation and the state error equation for the update of the contrast sources are defined respectively as follows

[TABLE]

with $p=1,2,3,\cdots,P$, $i=1,2,3,\cdots,I$. Here, $(\ell-1/2)$ means the update of the contrast sources taking place after the $(\ell-1)$-th iteration and before the $\ell$-th update of the contrast. Equation (2a) is in the measurement domain and Equation (2b) is in the field domain. The latter is used to monitor the behavior of the solution in the field domain and check if it satisfies Maxwell’s equations. Since Equation (2b) is always not perfectly satisfied, the “$=$” is supposed to be a “$\approx$”. Note that the solution is monitored only in the field domain, the mismatch in the field domain should also be monitored back in the measurement domain. Otherwise, the design of the cost functional is logically not complete. In order to fill this gap, we define a new equation in the measurement domain as follows

[TABLE]

with $p=1,2,3,\cdots,P$, $i=1,2,3,\cdots,I$. We refer to this equation as multi-frequency cross-correlated error equation. Consequently, the cost functional for the update of the contrast sources is defined as follows

[TABLE]

where, $\eta_{i}^{\mathcal{S}}$ and $\eta_{i,\ell-1}^{\mathcal{D}}$ are defined as

[TABLE]

respectively.

2.2.2 Updating the Contrast Sources

The gradient (Fréchet derivative) of the modified cost functional with respect to the contrast source $\bm{j}_{p,i}$ is

[TABLE]

Now suppose $\bm{j}_{p,i,\ell-1}$ and $\bm{\chi}_{i,\ell-1}$ are known, then we update $\bm{j}_{p,i,\ell-1}$ by

[TABLE]

where $\alpha_{p,i,\ell}$ is constant and the update direction $\bm{\nu}_{p,i,\ell}$ is chosen as the Polak-Ribière conjugate gradient directions, which is given by

[TABLE]

where

[TABLE]

The step size $\alpha_{p,i,\ell}$ is the minimizer of the cost functional

[TABLE]

See Appendix A for the derivation of $\alpha_{p,i,\ell}$. Following the update of the contrast sources, the total fields are updated by

[TABLE]

Consequently, the state error and the cross-correlated error are supposed to be updated as well: $\bm{\gamma}_{p,i,\ell-1}\rightarrow\bm{\gamma}_{p,i,\ell}$, $\bm{\xi}_{p,i,\ell-1}\rightarrow\bm{\xi}_{p,i,\ell}$.

2.2.3 Modified Cost Functional, $\mathcal{C}_{\text{MF-CC-CSI},\ell}$, and Updating the Contrast

Suppose $\bm{\chi}_{i,\ell-1}$ is known and consider the relation

[TABLE]

where, $\Re$ and $\Im$ respectively represent the real parts and imaginary parts of complex numbers, we therefore define $\bm{\chi}_{\ell}=\bm{\chi}_{1,\ell}$. Once the contrast source $\bm{j}_{p,i,\ell}$ is determined, we update the contrast by

[TABLE]

where, $\beta_{\ell}$ is the step size, and the update direction, $\bm{\nu}^{\chi}_{\ell}$, is chosen to be the Polak-Ribière conjugate gradient directions, which is given by

[TABLE]

where, $\bm{g}^{\chi}_{\ell}$ is the preconditioned gradient of the modified multi-frequency cost functional for updating the contrast

[TABLE]

Here,

[TABLE]

[TABLE]

Specifically, $\bm{g}^{\chi}_{\ell}$ is given by

[TABLE]

where,

[TABLE]

The step size $\beta_{\ell}$ is determined by minimizing the updated cost function in the formulation of

[TABLE]

This is a problem of finding the minimum of a single-variable function, which can be solved efficiently by the Brent’s method [27, 28]. It is worth noting that the objects are assumed to be isotropic in this paper. Namely, we assume $\bm{\chi}_{x}=\bm{\chi}_{y}$ when TE-polarized data is processed.

2.2.4 Initialization

If no a priori information about the objects is available, the contrast source is initialized using the values obtained by back-propagation [11]

[TABLE]

and the starting value of the total field is

[TABLE]

The contrast is initialized by (see [25])

[TABLE]

Since free space is considered in the following examples, the contrast must have non-negative real part and non-positive imaginary part. Such properties are enforced in the following examples by simply setting the negative real part and the positive imaginary part of the contrast to zero following each update of the contrast for both MF-CC-CSI and MF-MR-CSI. In addition, the contrast sources and the contrast are initialized by Equation (20) and Equation (22), respectively, for fair comparison.

3 Numerical Simulation

In this section, both MF-CC-CSI and MF-MR-CSI are tested for comparison with a 2-D benchmark problem – the “Austria” profile, which was also used in [29, 30, 31, 24]. The objects consist of two disks and one ring. Let us first establish our coordinate system such that the $z$-axis is parallel to the axis of the objects. The disks of radius 0.2 m are centred at ($-0.3$, 0.6) m and ($0.3$, 0.6) m. The ring is centred at (0,$-0.2$) m, and it has an exterior radius of 0.6 m and an inner radius of 0.3 m. We assume that the three cylinders are made of same material. Two sets of dielectric parameters are considered, which are $\varepsilon_{r}=3$, $\sigma=5$ mS/m, and $\varepsilon_{r}=10$, $\sigma=10$ mS/m, respectively. Since we assume a free space environment, the contrast values are $\Delta\varepsilon_{r}=2$, $\Delta\sigma=5$ mS/m, and $\Delta\varepsilon_{r}=9$, $\Delta\sigma=10$ mS/m, respectively. The forward scattering problems are solved by a MATLAB-based “MaxwellFDFD” package [32] with a fine and non-uniform grid size of $\lambda/(45\sqrt{\varepsilon_{r}})$.

3.1 Measurement configuration

In order to approach the realistic situation, we selected a measurement configuration which is similar to the experiment conducted by the Remote Sensing and Microwave Experiments Team at the Institut Fresnel, France [33]. An emitter is fixed on the circular rail, while a receiver is rotating with the arm around a vertical cylindrical target. The targets rotated from 0*∘* to 350*∘* in steps of 30*∘, and the receiver rotated from 60∘* to 300*∘* in steps of 5*∘*. The distance from the emitter/receiver to the origin is $3$ m. Simulation is sequentially done without and with the objects at five frequencies of 0.1 GHz, 0.2 GHz, 0.3 GHz, 0.4 GHz, and 0.5 GHz. Scattered field data is obtained by subtracting the incident field data, $\bm{y}^{\text{inc}}_{p}$, from the total field data, $\bm{y}^{\text{tot}}_{p}$.

3.2 Inversion results

Although the incident fields in the inversion domain can be easily obtained in numerical simulations, it is not available in real applications. Therefore, we consider the modeling of the incident fields using the collected incident field data. We select the approach reported in [25]. Specifically, the transmitting antenna is approximated by line source parallel to the cylindrical objects. The incident field by the line source is calibrated by multiplying a complex ratio which is calculated using only the value of the field when the transmitting and the receiving antenna are in opposite directions. For each frequency and source position, one complex calibration factor is determined. To simulate the inevitable measurement error in real experiments, additive Gaussian random noise, $\bm{n}_{p}$, is added directly to the scattered field data. In addition, the noise is also added to the total field data and the incident field data by $\bm{y}^{\text{tot}}_{p}+\bm{n}_{p}/2$ and $\bm{y}^{\text{inc}}_{p}-\bm{n}_{p}/2$, respectively. In doing so, the disturbance of the noise is reflected not only in the scattered field data, but also in the modelling of the incident fields, which better matches the real situation. To appraise the inversion accuracy, let us define the reconstruction error as follows

[TABLE]

where, $\hat{\bm{\chi}}$ and $\bm{\chi}$ are the reconstructed contrast and the exact contrast, respectively. Since large values of frequency tolerates large errors of estimated conductivity (see Subsection 4.2 for more discussion), $\hat{\bm{\chi}}$ and $\bm{\chi}$ in Equation (23) correspond to the highest frequency to ensure that the imaginary part of the contrast error is not over amplified in the evaluation of the inversion accuracy. In the following two cases, the inversion domain is restricted to [$-1.20$, $1.20$] $\times$ [$-1.20$, $1.20$] m2 and discretized into uniform grids with a size of $0.30\times 0.30$ mm2.

3.2.1 Case 1: $\Delta\varepsilon_{r}=2$, $\Delta\sigma=5$ mS/m

Let us first consider the lower contrast case, i.e., $\Delta\varepsilon_{r}=2$, $\Delta\sigma=5$ mS/m. MF-CC-CSI and MF-MR-CSI are used to process the TM- and TE-polarized numerical data respectively with 2048 iterations. Different frequency bands of 0.1-0.5 GHz and 0.3-0.5 GHz and different SNRs of 30 dB and 10 dB are tested. Figure 1 shows the reconstruction error curves in term of the iteration number, which indicates obviously that MF-CC-CSI obtains less reconstruction errors compared to MF-MR-CSI when same frequency band and SNR are used. It is also easy to observe that noise indeed leads to a degradation of the inversion accuracy (see the error curves of MF-CC-CSI shown in the right subfigures of Figure 1 when processing the 10 dB SNR data). As the frequency goes up to 0.3-0.5 GHz, the reconstruction error of MF-MR-CSI remains close to 1, indicating that MF-MR-CSI completely fails to invert the data in this frequency band.

3.2.2 Case 2: $\Delta\varepsilon_{r}=9$, $\Delta\sigma=10$ mS/m

In this case, the contrast increases to $\Delta\varepsilon_{r}=9$, $\Delta\sigma=10$ mS/m, which is supposed to be more challenging. We have also considered the frequency band 0.1-0.5 GHz, and from Figure 2 we observe that both MF-CC-CSI and MF-MR-CSI obtained large reconstruction errors. However, the former obtained less and smoothly decreasing reconstruction errors than the latter. Now let us remove the high frequency data and consider the frequency band 0.1-0.2 GHz, one can see from the bottom two subfigures of Figure 2 that the inversion accuracy of MF-CC-CSI has been improved; The reconstruction error curve of MF-MR-CSI started to decrease when processing the 0.1-0.2 GHz 30 dB synthetic data. Both of the two cases demonstrate that MF-CC-CSI is a more reliable inversion method in comparison to MF-MR-CSI.

3.3 Analysis of global convergence capability

In this subsection, let us have more insights on the enhanced global convergence capability of MF-CC-CSI by analyzing the behavior of the modified cost function in the solution space. Specifically, let us first define the solution as follows

[TABLE]

The contrast sources are not considered here because, knowing the contrast and the total fields, it can be simply obtained. With two sampling parameters, $\beta_{1}$ and $\beta_{2}$, the solution space is defined as follows

[TABLE]

where, $\bm{x}^{\text{act}}$ is the actual solution, while $\bm{x}^{\text{opt}}_{\text{MF-CC-CSI}}$ and $\bm{x}^{\text{opt}}_{\text{MF-MR-CSI}}$ are the solutions of MF-CC-CSI and MR-MR-CSI, respectively. Now let us consider the 0.1-0.5 GHz TM example in Case 1 and plot the cost function value versus two solution-space sampling parameters, $\mathcal{C}_{\text{MF-CC-CSI}}\left(\beta_{1},\beta_{2}\right)=\mathcal{C}_{\text{MF-CC-CSI}}\left\{\bm{x}\left(\beta_{1},\beta_{2}\right)\right\}$, with $\beta_{1},\beta_{2}\in[-1.5,1.5]$. Figure 3 (a) and (b) show the behavior of the cost function, $\log_{10}\left\{\mathcal{C}_{\text{MF-CC-CSI}}\left(\beta_{1},\beta_{2}\right)\right\}$, with different SNRs, 30 dB and 10 dB, respectively. From Figure 3 we see that the MF-MR-CSI converges to a local minimum, while $\bm{x}^{\text{opt}}_{\text{MF-CC-CSI}}$ belongs to the attraction basin of the global optimum of the cost function, $\bm{x}^{\text{act}}$. In addition, Figure 3 indicates that the modified cost function still shows the multiminima property and a good initial guess is still critical for preventing the occurrence of false solutions. By comparison of Figure 3 (a) and Figure 3 (b) we can also observe that the behavior of the cost function is affected by the noise level in the measurement data.

Here we finally remark that a good termination condition is required for MF-CC-CSI. Otherwise, the inversion accuracy may deteriorate after a mount of iterations when the disturbance of noise to the measurement data is not negligible (see the top right figures in Figure 1 and Figure 2). Now let us consider the 0.1-0.5 GHz TE example in Case 1, Figure 4 (a) and (b) show the behavior of the reconstruction error and the cost function ($\log_{10}\left\{\mathcal{C}_{\text{MF-CC-CSI}}\right\}$) of MF-CC-CSI with SNR$=30$ dB and 10 dB, respectively. One can observe that a threshold can hardly be determined to terminate the inversion algorithm without considering the noise level. In addition, the value of the contrast and the measurement configuration also affect the cost function curve and the convergence rate of the iterative inversion algorithms.

One reasonable strategy in practice is to first estimate the noise level. Moreover, according to some a priori information, it is possible to estimate the range of the contrast value. Based upon these information, one can do simulation using a typical benchmark inverse problem for a specific application, then obtain approximately how many iterations are required to get a reasonable inverted result for a specific application. Take the following experimental data inversion as an example. Considering the fact that the measurement configuration, the noise level and the range of the contrast value are similar to the numerical simulation, let us run 2048 iterations in the following experimental data inversion.

4 Experimental Data

4.1 Configuration

In this section, we consider the experiment carried out by Institut Fresnel in the year of 2005 [26]. In this experiment, the receiver stays in the azimuthal plane (xoy) and rotates along two-thirds of a circle from 60∘ to 300∘ with the angular step being 1∘. The source antenna stays at the fixed location ($\theta=0^{\circ}$). The targets are rotated to obtain different illumination incidences. The distance from the transmitter and receiver to the centre of the target domain is 1.67 m.

To avoid the redundancy of our discussion, let us select the most complicated configuration that corresponds to the datasets FoamTwinDielTM and FoamTwinDielTE. The targets consist of one larger circular dielectric cylinder with a smaller one embedded inside and a smaller adjacent one outside (see Figure 5(c) in [26]). Two smaller circular dielectric cylinders have relative permittivity values of $\varepsilon_{r}=3\pm 0.3$ while the larger one has a relative permittivity value of $\varepsilon_{r}=1.45\pm 0.15$. In this configuration, the targets are rotated from 0∘ to 340∘ with angular step of 20∘. $241\times 18$ measurements are obtained at each frequency (in TE polarization, only the component orthogonal to both the invariance axis of the cylinder and the direction of illumination is measured). The measurement configuration is given by Figure 1 of [26]. To increase the inversion difficulty, let us assume the data in the low frequency band is not available anymore, and we only have the measurement data at 7 GHz, 8 GHz, 9 GHz, and 10 GHz.

4.2 Inversion Results

First, let us model the incident fields with the same approach as we did in the numerical simulation. To reduce the computational burden, we restrict the inversion domain to [$-82.5$, $82.5$] $\times$ [$-97.5$, $67.5$] mm2. The inversion domain is discretized with 1.5 $\times$ 1.5 mm2 grids. The multi-frequency datasets FoamTwinDielTM and FoamTwinDielTE at 7 GHz, 8 GHz, 9 GHz, and 10 GHz were processed by MF-CC-CSI and MF-MR-CSI, respectively. Both of them were terminated after 2048 iterations. In addition, we also did a simulation to generate the perfect data of the same targets in the same configuration. The only thing different is that the distance between antennas and the centre of the inversion domain was shortened from 1.67 m to 0.20 m. Otherwise, the scattering domain with the frequency up to 10 GHz is too huge for a standard desktop.

Figure 5 (FoamTwinDielTM) and Figure 6 (FoamTwinDielTE) show the inverted results by processing the multi-frequency dataset (Line 1: MF-CC-CSI; Line 2: MF-MR-CSI) and the noise-free synthetic data (Line 3: MF-CC-CSI; Line 4: MF-MR-CSI) at 7 GHz, 8 GHz, 9 GHz, and 10 GHz with 2048 iterations. One can see that MF-CC-CSI successfully reconstructed the three cylinders in good estimation accuracy with both the experimental data and the noise-free synthetic data, while MF-MR-CSI only obtained good inversion results with the noise-free TE-polarized synthetic data. The artifacts of the inverted contrast conductivity (see the right figure in Line 3 of Figure 5) is inevitable and so far can only be explained as a reconstruction error to the best of our knowledge. The higher the frequency is, the larger such artifacts would be. This is easy to understand by noting that $\bm{\chi}_{i}=\Delta\varepsilon_{r}-\text{i}\Delta\sigma/\omega_{i}$, i.e., a large value of the angular frequency, $\omega_{i}$, tolerates a large error of the contrast conductivity, $\Delta\sigma$. The synthetic data enables us to obtain the reconstruction error curves of both methods, which is shown in Figure 7. The error curves show us the inversion performance of both methods during the whole process of inversion, demonstrating again the advantage of MF-CC-CSI in more complicated scenarios. One can also see that the multiplicative regularization can improve the inversion accuracy indeed when the contrast can be reconstructed to a reasonable accuracy. The reason has been fully discussed in [24].

5 Conclusion

In this paper, a multi-frequency version of the cross-correlated contrast source inversion (CC-CSI) method is proposed. By processing the numerical data and experimental data in both transverse magnetic (TM) and transverse electric (TE) polarizations, we have demonstrated the advantage of the proposed multi-frequency CC-CSI (MF-CC-CSI) method over the multi-frequency multiplicative regularized contrast source inversion (MF-MR-CSI) method. In addition, the introduction of the cross-correlated error changed the classical design of the cost functional, which is of great significance for improving the inversion performance of other algorithms that are based on the classical cost functional design.

Acknowledgment

The experimental data used in the manuscript is available via this link: http://iopscience.iop.org/article/10.1088/0266-5611/21/6/S09/data. The numerical data is available via this link: https://surfdrive.surf.nl/files/index.php/s/pzy4U2i6caNauKd. The full code package of CC-CSI is now available at https://github.com/TUDsun/CC-CSI.

In the end, we would like to show our appreciation to the Fresnel Institute for providing open access to the experimental data.

Appendix A Derivation of the parameter, $\alpha_{p,i,\ell}$

First, let us rewrite the cost function $\left.\mathcal{C}_{\text{MF-CC-CSI},\ell-1/2}\right|_{\bm{j}_{p,i}=\bm{j}_{p,i,\ell-1}+\alpha_{p,i}\bm{\nu}_{p,i,\ell}}$ as follows

[TABLE]

Obviously, it can be further simplified in the form of

[TABLE]

Therefore, we have

[TABLE]

Note that

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

it is easy to obtain

[TABLE]

where, $a_{p,i,2}$, $b_{p,i,2}$, and $c_{p,i,2}$ are given by Equation (29), Equation (31), and Equation (33), respectively.