Cross-correlated Contrast Source Inversion

TL;DR

This paper introduces a cross-correlated contrast source inversion method that enhances robustness and accuracy in electromagnetic imaging by interrelating state and data errors, demonstrated through 2-D benchmark problems.

Contribution

It proposes a novel cross-correlated cost functional for CSI, improving robustness and accuracy without significant computational overhead.

Findings

Better robustness and higher inversion accuracy than classical CSI and MR-CSI.

Effective in 2-D benchmark problems with transverse magnetic and electric waves.

Gradient calculation remains computationally efficient.

Abstract

In this paper, we improved the performance of the contrast source inversion (CSI) method by incorporating a so-called cross-correlated cost functional, which interrelates the state error and the data error in the measurement domain. The proposed method is referred to as the cross-correlated CSI. It enables better robustness and higher inversion accuracy than both the classical CSI and multiplicative regularized CSI (MR-CSI). In addition, we show how the gradient of the modified cost functional can be calculated without significantly increasing the computational burden. The advantages of the proposed algorithms are demonstrated using a 2-D benchmark problem excited by a transverse magnetic wave as well as a transverse electric wave, respectively, in comparison to classical CSI and MR-CSI.

Tables

Table 1. Table 1: The computational times

	CSI	MR-CSI	CC-CSI
TM	1361.2	1419.7	1429.6
TE	2931.9	3040.9	3146.2

Equations

$$\nabla \times {\mathbf{\mu }}^{-1}\nabla \times {\mathbf{E}}_p^{\text{sct}}-{\omega }^2{\mathbf{ϵ}}_b{\mathbf{E}}_p^{\text{sct}}={\omega }^2\mathbf{\chi }{\mathbf{E}}_p,$$

$$\mathbf{A}{\mathbf{e}}_p^{\text{sct}}={\omega }^2\mathbf{\chi }{\mathbf{e}}_p,p=1,2,\dots ,P,$$

$${\mathbf{f}}_p={\mathcal{M}}_{\mathcal{S}}{\mathbf{A}}^{-1}{\omega }^2\mathbf{\chi }{\mathbf{e}}_p,\mathbf{x}\in \mathcal{S},p=1,2,\dots ,P,$$

$${\mathbf{\rho }}_p={\mathbf{f}}_p-{\mathcal{M}}_{\mathcal{S}}{\mathbf{A}}^{-1}{\mathbf{j}}_p,p=1,2,\dots ,P,$$

$${\mathbf{\gamma }}_p=\mathbf{\chi }{\mathbf{e}}_p^{\text{inc}}+\mathbf{\chi }{\mathcal{M}}_{\mathcal{D}}{\mathbf{A}}^{-1}{\mathbf{j}}_p-{\mathbf{j}}_p,p=1,2,\dots ,P,$$

$$C_{\text{CSI}}^{\mathbf{j}}({\mathbf{\chi }}_{n-1},{\mathbf{j}}_p)={\eta }^{\mathcal{S}}\sum _{p=1}^P\Vert {\mathbf{\rho }}_p{\Vert }_{\mathcal{S}}^2+{\eta }^{\mathcal{D}}\sum _{p=1}^P\Vert {\mathbf{\gamma }}_p{\Vert }_{\mathcal{D}}^2.$$

$$C_{\text{CSI}}^{\mathbf{\chi }}(\mathbf{\chi },{\mathbf{j}}_{p,n})={{\eta }^{\mathcal{D}}\sum _{p=1}^P\Vert {\mathbf{\gamma }}_p{\Vert }_{\mathcal{D}}^2|}_{{\mathbf{j}}_p={\mathbf{j}}_{p,n}}.$$

$$C_{\text{MR-CSI}}^{\mathbf{\chi }}(\mathbf{\chi },{\mathbf{j}}_{p,n})={C_{\text{CSI}}^{\mathbf{j}}(\mathbf{\chi },{\mathbf{j}}_p)\times \frac{1}{V}{\int }_{\mathcal{D}}\frac{\Vert \nabla \mathbf{\chi }{\Vert }^2+{\delta }_{n-1}^2}{\Vert \nabla {\mathbf{\chi }}_{n-1}{\Vert }^2+{\delta }_{n-1}^2}d\mathbf{x}|}_{{\mathbf{j}}_p={\mathbf{j}}_{p,n}}.$$

$${\delta }_{n-1}^2=C_{\text{CSI}}^{\mathbf{\chi }}({\mathbf{\chi }}_{n-1},{\mathbf{j}}_{p,n}){\Delta }^2,$$

$$\kappa (\mathbf{\Phi }):=\frac{{\sigma }_{\max }(\mathbf{\Phi })}{{\sigma }_{\min }(\mathbf{\Phi })}.$$

$${\mathbf{\xi }}_p={\mathbf{f}}_p-\mathbf{\Phi }(\mathbf{\chi }{\mathbf{e}}_p^{\text{inc}}+\mathbf{\chi }{\mathcal{M}}_{\mathcal{D}}{\mathbf{A}}^{-1}{\mathbf{j}}_p).$$

$$C_{\text{CC-CSI}}^{\mathbf{j}}({\mathbf{\chi }}_{n-1},{\mathbf{j}}_p)={\eta }^{\mathcal{S}}\sum _{p=1}^P\Vert {\mathbf{\rho }}_p{\Vert }_{\mathcal{S}}^2+{{\eta }^{\mathcal{D}}\sum _{p=1}^P\Vert {\mathbf{\gamma }}_p{\Vert }_{\mathcal{D}}^2+{\eta }^{\mathcal{S}}\sum _{p=1}^P\Vert {\mathbf{\xi }}_p{\Vert }_{\mathcal{S}}^2|}_{\mathbf{\chi }={\mathbf{\chi }}_{n-1}}.$$

$${\mathbf{g}}_{p,n}=-2{\eta }^{\mathcal{S}}{\mathbf{\Phi }}^H{\mathbf{\rho }}_p+2{\eta }^{\mathcal{D}}{\left(\mathbf{\chi }{\mathcal{M}}_{\mathcal{D}}{\mathbf{A}}^{-1}-\mathbf{I}\right)}^H{\mathbf{\gamma }}_p-{2{\eta }^{\mathcal{S}}{\left(\mathbf{\Phi }\mathbf{\chi }{\mathcal{M}}_{\mathcal{D}}{\mathbf{A}}^{-1}\right)}^H{\mathbf{\xi }}_p|}_{\mathbf{\chi }={\mathbf{\chi }}_{n-1},{\mathbf{j}}_p={\mathbf{j}}_{p,n-1}}.$$

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

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

$$C_{\text{CC-CSI}}^{\mathbf{\chi }}(\mathbf{\chi },{\mathbf{j}}_{p,n})={\eta }^{\mathcal{D}}\sum _{p=1}^P\Vert {\mathbf{\gamma }}_p{\Vert }_{\mathcal{D}}^2+{\eta }^{\mathcal{S}}\sum _{p=1}^P\Vert {\mathbf{\xi }}_p{\Vert }_{\mathcal{S}}^2$$

$${\mathbf{\chi }}_n={\mathbf{\chi }}_{n-1}+{\beta }_n{\mathbf{\nu }}_{\mathbf{\chi },n},$$

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

$${\mathbf{g}}_{\mathbf{\chi },n}=\frac{2{\eta }^{\mathcal{D}}{\sum }_{p=1}^P\overline{{\mathbf{e}}_{p,n}}{\mathbf{\gamma }}_p-2{\eta }^{\mathcal{S}}{\sum }_{p=1}^P\overline{{\mathbf{e}}_{p,n}}{\mathbf{\Phi }}^H{\mathbf{\xi }}_p}{{\sum }_{p=1}^P{\left|{\mathbf{e}}_{p,n}\right|}^2},$$

$$C_{\text{CC-CSI},n}^{\mathbf{\chi }}=\frac{\sum _{p=1}^P{\Vert ({\mathbf{\chi }}_{n-1}+{\beta }_n{\mathbf{\nu }}_{\mathbf{\chi },n}){\mathbf{e}}_{p,n}-{\mathbf{j}}_{p,n}\Vert }_{\mathcal{D}}^2}{\sum _{p=1}^P{\Vert ({\mathbf{\chi }}_{n-1}+{\beta }_n{\mathbf{\nu }}_{\mathbf{\chi },n}){\mathbf{e}}_p^{\text{inc}}\Vert }_{\mathcal{D}}^2}+{\eta }^{\mathcal{S}}\sum _{p=1}^P{\Vert {\mathbf{f}}_p-\mathbf{\Phi }({\mathbf{\chi }}_{n-1}+{\beta }_n{\mathbf{\nu }}_{\mathbf{\chi },n}){\mathbf{e}}_{p,n}\Vert }_{\mathcal{S}}^2.$$

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

$${\mathbf{\chi }}_0=\sum _{p=1}^P{\mathbf{j}}_{p,0}\overline{{\mathbf{e}}_{p,0}}/\sum _{p=1}^P{\mathbf{e}}_{p,0}\overline{{\mathbf{e}}_{p,0}},$$

$${\mathbf{A}}^T{\mathbf{\phi }}_m={{\mathcal{M}}_{\mathcal{S}}^m}^T,m=1,2,\dots ,M,$$

$$\Delta \le \frac{{\lambda }_0}{15\sqrt{{\epsilon }_r}},$$

$$\Delta \le \frac{{\lambda }_0}{15\sqrt{\max \{{\epsilon }_r\}}}.$$

$$\text{err}=\frac{\Vert \mathbf{\chi }-\hat{\mathbf{\chi }}{\Vert }^2}{\Vert \mathbf{\chi }{\Vert }^2}.$$

$${\mathbf{f}}_{p,noise}={\mathbf{f}}_p+\zeta \times \underset{m}{\max }\{|f_{p,m}|\}({\mathbf{n}}_1+\text{i}{\mathbf{n}}_2),$$

$$\frac{T_{\text{CC-CSI}}-T_{\text{CSI}}}{T_{\text{CSI}}}\times 100\%,$$

$$\begin{array}{rl}& C_{\text{CC-CSI}}^{\mathbf{j}}({\mathbf{\chi }}_{n-1},{\mathbf{j}}_{p,n-1}+{\alpha }_p{\mathbf{\nu }}_{p,n})={\eta }_{n-1}^{\mathcal{D}}\sum _{p=1}^P{\Vert {\mathbf{\gamma }}_{p,n-1}+{\alpha }_p({\mathbf{\chi }}_{n-1}{\mathcal{M}}_{\mathcal{D}}{\mathbf{A}}^{-1}-\mathbf{I}){\mathbf{\nu }}_{p,n}\Vert }^2+\\ & {\eta }^{\mathcal{S}}\sum _{p=1}^P{\Vert {\mathbf{\rho }}_{p,n-1}-{\alpha }_p\mathbf{\Phi }{\mathbf{\nu }}_{p,n}\Vert }^2+{\eta }^{\mathcal{S}}\sum _{p=1}^P{\Vert {\mathbf{\xi }}_{p,n-1}-{\alpha }_p\mathbf{\Phi }{\mathbf{\chi }}_{n-1}{\mathcal{M}}_{\mathcal{D}}{\mathbf{A}}^{-1}{\mathbf{\nu }}_{p,n}\Vert }^2\end{array}$$

$$\begin{array}{rl}{C}_{\text{CC-CSI}}^{\mathbf{j}}({\mathbf{\chi }}_{n-1},{\mathbf{j}}_{p,n-1}+{\alpha }_p{\mathbf{\nu }}_{p,n})=& a_{p,2}{\alpha }_{p,n}^2+a_{p,1}{\alpha }_{p,n}+a_{p,0}+b_{p,2}{\alpha }_p^2+b_{p,1}{\alpha }_p+b_{p,0}+\\ & c_{p,2}{\alpha }_{p,n}^2+c_{p,1}{\alpha }_p+c_{p,0}.\end{array}$$

Full text

1 Introduction

Inversion techniques have been applied extensively in many fields, e.g. radar imaging [1], seismic imaging [2], medical imaging [3, 4], and so forth. Developments in inversion techniques and research are focused on computational efficiency, the incorporation of a priori information to circumvent computational artifacts, and the calibration to the real antenna radiating pattern especially in near-field scenarios [5, 6, 7]. Methods to solve the inverse scattering problems include non-iterative methods, e.g. linear sampling method [8, 9], and iterative methods [10]. The Contrast Source Inversion (CSI) method is an iterative frequency domain inversion method to retrieve the value of the contrast of scattering objects, which was first proposed by van den Berg et al. [10], and was later applied to subsurface object detection in combination with integral equations based on the Electric Field Integral Equation (EFIE) formulation, see Kooij et al. [11]. A priori information was introduced in the form of mathematical regularization constraints like the positivity constraints of the material properties and the Total Variation (TV) constraint in [12] to further enhance the performance. In [12], a multiplicative regularized CSI (MR-CSI) method is proposed, in which the estimation of the tuning parameter is avoided. In addition, Crocco et al. [13, 14], applied the so-called contrast source-extended Born-model to 2-D subsurface scattering problems. Later the CSI technique was introduced in combination with the finite-difference frequency domain (FDFD) scheme by Abubakar et al. [15, 16]. The scheme based on the FDFD technique turned out to have computational advantages compared to EFIE scheme, especially if a non-homogeneous background, like the half-space configuration in ground penetrating radar (GPR), is required in the inversion. For a more accurate representation of complex geometry, finite-element method (FEM) was introduced and combined with CSI by Zakaria et al., and the 2-D inversion results with the transverse magnetic (TM) wave and the transverse electric (TE) wave can be found in [17] and [18], respectively. FEM was applied as well to the contrast source-extended Born method in [19].

One obvious drawback of the iterative methods is that a good initial guess must be provided beforehand to ensure the iterative inverting process converges to the global optimal solution. The reason is that there is more than one variable needed to be estimated during the inverting process. More specifically, the contrast and the contrast source are both unknown, and a less accurate initial guess is more likely to give a false gradient and thus leads the iterative inverting process to a local optimal solution. To overcome this drawback, the hybrid inversion schemes have been considered, which first recovered the shape of the scatterers faithfully by sampling-type technique, and then estimated their dielectric properties or improved the shape iteratively. The same idea can be found in the recently published papers [20, 21, 22, 23].

For iterative inversion methods, a cost function is normally needed which consists of the data error and the state error [24]. On one hand, it ensures that the algorithm fits the measurement data. On the other hand, it tends to optimize the estimation of the contrast and the contrast source to satisfy the Maxwell equations. In this paper, we show that a minor state error can become large when mapped into the measurement domain because of the ill-posedness of the inverse scattering problem. In another word, there might be cases in which both the state error and data error are minimized sufficiently, while the state error is still large when mapped into the measurement domain, which indicates that the estimated contrast is not the global optimal solution. Inspired by this fact, we introduce a new error equation that interrelates the state error and the data error in the measurement domain and modifies the cost functional accordingly. In doing so, the state error and the data error are cross-correlated, and the inverting process is stabilized by minimizing the state error not only in the field domain, but also in the measurement domain. We refer to the proposed algorithm as cross-correlated contrast source inversion (CC-CSI) method. In addition, we also show how the gradient of the new cost functional can be calculated without significantly increasing the computational complexity. The performance of the proposed method is investigated based on a 2-D benchmark problem excited by a TM-polarized wave and a TE-polarized wave, respectively. As can be observed from the results, the CC-CSI method shows better robustness and higher inversion accuracy than classical CSI and MR-CSI. Since the Maxwell’s equations are formulated in a three-dimensional FDFD formulation, it is already applicable to the reconstruction of future 3-D scattering objects.

The remainder of the paper is organized as follows: Section 2 gives the problem statement and the introduction of classical CSI and MR-CSI. The cross-correlated CSI method is introduced in Section 3. Simulation results based on a 2-D benchmark problem are given in Section 4, in which the performance investigation of CC-CSI in comparison to classical CSI and MR-CSI is fully discussed111The CC-CSI package is available at https://github.com/TUDsun/CC-CSI, in which MR-CSI and CSI are also contained for comparison.. Finally, we give our conclusions in Section 5.

2 Problem statement and Classical CSI Based on FDFD

2.1 Problem statement

We consider a scattering configuration as depicted in Figure 1,

which consists of a bounded, simply connected, inhomogeneous background domain $\mathcal{D}$, which is also referred to as field domain in this paper. The domain $\mathcal{D}$ contains an object, whose location and complex permittivity profile are unknown. The measurement domain $\mathcal{S}$ contains the sources and receivers. The sources are denoted by the subscript $p$ in which $p\in\{1,2,3...,P\}$, the receivers are denoted by the subscript $q$ in which $q\in\{1,2,3,...,Q\}$. Sources and receivers that have equal subscripts are located at 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 $\text{e}^{\text{i}\omega t}$ is considered in this paper, where $\text{i}^{2}=-1$.

In our notation for the vectorial quantities, we use a bold notation which represents a vector with three components. The general mathematical representations presented are consistent with any 3-D configuration, in which the 2-D TE and TM excitations are a special case, resulting in vectors containing zero elements. According to the linear relation among the incident electric field $\bm{E}^{\text{inc}}_{p}$, the scattered electric field $\bm{E}^{\text{sct}}_{p}$ and the total electric fields $\bm{E}_{p}$, $\bm{E}^{\text{sct}}_{p}=\bm{E}_{p}-\bm{E}^{\text{inc}}_{p}$, the scattering equation with respect to the scattered electric field $\bm{E}^{\text{sct}}_{p}$ can be easily obtained which is [25]

[TABLE]

with $p=1,2,\dots,P$. Here, $\bm{\mu}$ represents the permeability of the background and is assumed to be equal to the permeability of the free space in this paper; $\omega$ is the angular frequency; the permittivity $\bm{\varepsilon}_{b}$ and the conductivity $\bm{\sigma}_{b}$ of the background are incorporated into the complex permittivity $\bm{\epsilon}_{b}$ satisfying $\bm{\epsilon}_{b}=\bm{\varepsilon}_{b}-\text{i}\bm{\sigma}_{b}/\omega$.

The scattering equation can be further formulated based on FDFD scheme by [25]

[TABLE]

where $\bm{A}$ is the FDFD stiffness matrix; $\bm{e}_{p}^{\text{sct}}$ and $\bm{e}_{p}$ are the scattered electric field and the total electric field in the form of a column vector, respectively; and $\bm{\chi}$ is the contrast consisting of the difference of the permittivity $\bm{\varepsilon}_{\text{c}}$ and the difference of the conductivity $\bm{\sigma}_{\text{c}}$ with the relation of $\bm{\chi}=\bm{\varepsilon}_{\text{c}}-\text{i}\bm{\sigma}_{\text{c}}/\omega$. Then the solution of (2) is obtained by inverting the stiffness matrix $\bm{A}$, which yields $\bm{e}_{p}^{\text{sct}}=\bm{A}^{-1}\omega^{2}\bm{\chi}\bm{e}_{p}$. This leads to the data equation

[TABLE]

where $\mathcal{M}_{\mathcal{S}}$ is an operator that interpolates field values defined at the finite-difference grid points to the appropriate receiver positions.

In the remainder of this paper, $\omega^{2}$ is incorporated into $\bm{A}$ for the sake of conciseness. The inverse scattering problem is to reconstruct the contrast $\bm{\chi}$ as a function of space from the incomplete measured field data $\{\bm{f}_{p}\},\ p=1,2,\cdots,P$, which is full of challenges because of the nonlinearity and the ill-posedness.

2.2 Classical CSI and MR-CSI

Classical CSI is a method of iteratively minimizing a cost functional consisting of the data error and the state error for reconstructing the contrast source. The contrast is updated during the iterations. Specifically, the multiplication of the contrast $\bm{\chi}$ and the total field $\bm{e}_{p}$ is referred to as the contrast source, which is represented by $\bm{j}_{p}$. Then, the data error $\bm{\rho}_{p}$ and the state error $\bm{\gamma}_{p}$ are defined by

[TABLE]

and

[TABLE]

respectively. Here, $\mathcal{M}_{\mathcal{D}}$ is an operator that selects fields only inside the field domain $\mathcal{D}$. The cost functional is given by

[TABLE]

with $\bm{\chi}=\bm{\chi}_{n-1}$. Here, $1/\eta^{\mathcal{S}}=\sum_{p}\|\bm{f}_{p}\|^{2}_{\mathcal{S}}$, $1/\eta^{\mathcal{D}}=\sum_{p}\|\bm{\chi}\bm{e}_{p}^{\text{inc}}\|^{2}_{\mathcal{D}}$, and $\|\cdot\|_{\mathcal{S}}$ and $\|\cdot\|_{\mathcal{D}}$ represent the norms on the measurement space $L^{2}(\mathcal{S})$ and the field space $L^{2}(\mathcal{D})$, respectively. The contrast source $\bm{j}_{p}$ is iteratively optimized by minimizing the cost functional $C^{\bm{j}}(\bm{\chi},\bm{j}_{p})$, followed by the update of the contrast $\bm{\chi}$ which is done by minimizing the cost functional

[TABLE]

MR-CSI is the CSI method regularized with a multiplicative weighted total variation (TV) constraint, which was first introduced by van den Berg et al. [12]. In comparison to CSI, the contrast $\bm{\chi}$ in MR-CSI is now updated by minimizing the cost functional

[TABLE]

Here, $V=\int_{\mathcal{D}}d\bm{x}$, and $\delta^{2}_{n-1}$ are introduced for restoring the differentiability of the TV factor [12]. The value of $\delta^{2}_{n-1}$ is chosen to be large in the beginning of the optimization and small towards the end, which is given by

[TABLE]

where $\Delta$ denotes the mesh size of the discretized domain $\mathcal{D}$. It is worth noting that the contrast $\bm{\chi}$ is assumed to be isotropic in this paper. Namely, for the TE case, only one component of $\bm{\chi}$ is used in the TV penalty function. More details of classical CSI and MR-CSI can be found in [15, 25, 12, 26].

3 Cross-Correlated CSI

3.1 Motivation

As aforementioned, a good initial guess is very critical for ensuring that the iterative methods can successfully converge to the global optimal solution. This can be explained firstly by the fact that there are two unknown variables — the contrast and the contrast source. The less accurate the initial guess is, the more inaccurate the gradient with respect to the contrast source will be. Secondly, although classical CSI is able to minimize the data error $\bm{\rho}(\bm{j}_{p})$ by constraining the state error $\bm{\gamma}(\bm{\chi},\bm{j}_{p})$ at the same time, a global optimal solution can still not be guaranteed because of the ill-posedness of the inverse scattering problem. For simplicity, let us first define the measurement matrix as $\bm{\Phi}:=\mathcal{M}_{\mathcal{S}}\bm{A}^{-1}$. The condition number of matrix $\bm{\Phi}$ is further defined as

[TABLE]

where $\sigma_{\max}(\bm{\Phi})$ and $\sigma_{\min}(\bm{\Phi})$ are maximal and minimal singular values of $\bm{\Phi}$, respectively. As discussed in Subsection 4.2, the measurement matrix has a large condition number, which means a minor state error in the field space $L^{2}(\mathcal{D})$ may cause a large error in the measurement space $L^{2}(\mathcal{S})$. This potential mismatch cannot be reflected by the cost functional of the classical CSI method. Inspired by this fact, we came up with the idea of introducing the so-called cross-correlated cost functional.

3.2 Cross-correlated CSI

In this subsection, a new cost functional is proposed, which interrelates the mismatch of the state equation and the data error in the measurement space. This proposed algorithm is referred to as the cross-correlated contrast source inversion method.

Specifically, the cross-correlated error $\bm{\xi}_{p}$ is defined as

[TABLE]

Note that if the state error $\bm{\gamma}_{p}$ is zero, then theoretically we have $\bm{\rho}_{p}=\bm{\xi}_{p}$. In classical CSI, sufficiently minimizing the cost function of Eq. (6) does not necessarily mean that the cross-correlated error is sufficiently minimized. Therefore, the cost functional of the contrast source in the proposed CC-CSI method is modified and defined as

[TABLE]

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

[TABLE]

Here, $\bm{I}$ represents the identity matrix, and $(\cdot)^{H}$ is the conjugate transpose operator. Now suppose $\bm{j}_{p,n-1}$ and $\bm{\chi}_{n-1}$ are known, then we update $\bm{j}_{p}$ by

[TABLE]

where $\alpha_{p,n}$ is constant and the update directions $\bm{\nu}_{p,n}$ are functions of the position. The update directions are chosen to be the Polak-Ribière conjugate gradient directions, which are given by

[TABLE]

where $\left\langle\cdot\right\rangle_{\mathcal{D}}$ represents the inner product defined in the field space $L^{2}(\mathcal{D})$. The step size $\alpha_{p,n}$ can be explicitly found by minimizing the cost functional $C_{\text{CC-CSI}}^{\bm{j}}(\bm{\chi}_{n-1},\bm{j}_{p,n-1}+\alpha_{p,n}\bm{\nu}_{p,n})$ (see Appendix A for the derivation).

Once the contrast source $\bm{j}_{p,n}$ is determined, we update the contrast $\bm{\chi}$ by minimizing the cost functional of the contrast which is defined by

[TABLE]

with $\bm{j}_{p}=\bm{j}_{p,n}$. Specifically, $\bm{\chi}$ is updated via

[TABLE]

where $\beta_{n}$ is constant and the update directions $\bm{\nu}_{p,n}$ are chosen to be the Polak-Ribière conjugate gradient directions, which are given by

[TABLE]

where $\bm{g}_{\bm{\chi},n}$ is the preconditioned gradient of the contrast cost functional $C_{\text{CC-CSI}}^{\bm{\chi}}(\bm{\chi})$ defined as

[TABLE]

with $\bm{\chi}=\bm{\chi}_{n-1}$, $\bm{j}_{p}=\bm{j}_{p,n}$, where $\overline{(\cdot)}$ represents the conjugate operator. The step size $\beta_{n}$ is determined by minimizing the cost function

[TABLE]

This is a problem of finding the minimum of a single-variable function, and can be solved efficiently by the Brent’s method [27, 28].

The CC-CSI method is given in Algorithm 1, where $\Re\{\cdot\}$ represents the real part operator and, correspondingly, the imaginary part operator is represented by $\Im\{\cdot\}$. Since $\mathcal{M}_{\mathcal{D}}$ always exists together with the stiffness matrix $\bm{A}$, it is neglected for better readability in the remainder of this paper. It is worth noting that $\bm{\chi}$ is assumed to be isotropic in this paper. Therefore, we average the two components of the contrast $\bm{\chi}$ for TE case, and the three components of the contrast $\bm{\chi}$ for 3-D case, after each update of $\bm{\chi}$.

3.3 Initialization

If no a priori information about the objects is available, the contrast sources are initialized by (see [29, 10])

[TABLE]

which are obtained by back-propagation, multiplied by a weight to ensure that the data error is minimized. The contrast is initialized by (see [12])

[TABLE]

with $\bm{e}_{p,0}=\bm{e}_{p}^{\text{inc}}+\mathcal{M}_{\mathcal{D}}\bm{A}^{-1}\bm{j}_{p,0}$.

3.4 Computational complexity

In this subsection, we show that the CC-CSI method can be implemented without significantly increasing the computational complexity compared to the classical CSI method. Note that since the selecting matrix $\mathcal{M}_{\mathcal{S}}\in\mathbb{C}^{M\times N}$ has only $M<<N$ rows, the matrix $\bm{\Phi}$ can be calculated iteratively by solving $M$ linear systems of equations,

[TABLE]

where, $\mathcal{M}^{m}_{\mathcal{S}}$ is the $m^{th}$ row of the selecting matrix $\mathcal{M}_{\mathcal{S}}$, and $(\cdot)^{T}$ represents the transpose operator. The matrix $\bm{\Phi}$ is assembled by $\bm{\Phi}=[\bm{\varphi}_{1},\bm{\varphi}_{2},\dots,\bm{\varphi}_{M}]^{T}$. Since $\bm{\Phi}$ has only $M<<N$ rows, it is computationally much more efficient than the LU decomposition of the stiffness matrix $\bm{A}$ (if we use LU decomposition). This feature makes it suitable to be computed and stored beforehand, which is of great importance for real applications, especially for 3-D inverse scattering problems. Although CC-CSI requires more matrix-vector multiplications, the extra computational cost is not significant by noting that the matrix $\bm{\Phi}$ has only $M<<N$ rows. This is further demonstrated in Subsection 4.5.

4 Simulation results

4.1 Configuration

In this section, the proposed algorithm is tested with a 2-D benchmark problem – the “Austria” profile, which was also used in [30, 31, 26, 32]. Based on the benchmark problem, the performance of CC-CSI is analyzed in comparison to classical CSI and MR-CSI.

Specifically, the objects to be inverted 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. Belkebir and Tijhuis [30] and Litman et al. [31] have used 64 sources and 65 receivers on a circle of radius 3 m centred at (0, 0), while the inverting domain was discretized into 30 $\times$ 30 cells. Van den Berg et al. [26, 32] have taken 48 source/receiver stations, while the inverting domain was discretized into 64 $\times$ 64 cells. In our simulation, 36 source/receiver stations are used and uniformly distributed on the same circle, which means we have under-sampled this problem further. In our simulation, $P=Q=36$, viz., we have $36\times 36$ measurement data for TM case and $36\times 72$ measurement data for TE case. Same objects but of different relative permittivity have been considered. The conductivity is fixed at 10 mS/m, so we have the same attenuation, while the relative permittivity attains values – 2.0, 2.5, 3.0, and 3.5, respectively. The operating frequency is 300 MHz, therefore, the corresponding values of the contrast are $\chi=1.0-0.6\text{i}$, $\chi=1.5-0.6\text{i}$, $\chi=2.0-0.6\text{i}$ and $\chi=2.5-0.6\text{i}$, respectively. The original “Austria” profile is given in Fig. 2, in which the green dots represent the 36 source/receiver stations.

The forward EM scattering problem is solved by a MATLAB-based 3-D FDFD package “MaxwellFDFD” [33]. The $x$- and $y$-normal boundaries are covered by perfect matching layers (PML) to simulate the anechoic chamber environment (see the gray layers of Fig. 2 at the boundaries of the test domain), while the two $z$-normal boundaries are subject to periodic boundary conditions (PBC) to simulate the 2-D configuration. Line sources parallel to the $z$-axis are used to generate TM-polarized and TE-polarized incident wave. Non-uniform meshes are used to generate the scattered data, which means the testing domain is discretized with different mesh sizes determined by the distribution of the permittivity, viz., coarse meshes for low permittivity and fine meshes for high permittivity. The accuracy of the FDFD scheme is ensured by the following criterion [34]

[TABLE]

where, $\lambda_{0}$ is the wavelength in free space. Non-uniform meshes greatly reduce the computational burden for solving the forward scattering problem. In contrast, uniform meshes are used to invert the scattered data, since we do not know the distribution of the permittivity beforehand. To guarantee the inverting accuracy, the following condition is satisfied

[TABLE]

The scattered field is obtained by subtracting the incident field from the total field.

4.2 Condition number of the measurement matrix

In the following simulations, the inversion domain is restricted to the region [-1.5,1.5] $\times$ [-1.5,1.5] m2. The dimension of the mesh grid is 30 $\times$ 30 mm2. Thus we have specifically in this simulation, $\bm{\Phi}\in\mathbb{C}^{36\times 10000}$ for TM polarization, and $\bm{\Phi}\in\mathbb{C}^{72\times 20000}$ for TE polarization. The condition numbers of $\bm{\Phi}$ are $\kappa(\bm{\Phi})_{TM}=5.19\times 10^{2}$ for TM polarization and $\kappa(\bm{\Phi})_{TE}=8.25\times 10^{7}$ for TE polarization. As one can see the condition number of the matrix $\bm{\Phi}$ is large for both TM and TE polarization, indicating that an error in the contrast sources $\bm{j}_{p}$ may cause an increased error in the measurement data $\bm{f}_{p}$. In addition, compared to TM polarization, TE polarization is more ill-conditioned because $\kappa(\bm{\Phi})_{TE}$ is much larger than $\kappa(\bm{\Phi})_{TM}$, due to the different formulation of the scattering equations. This implies that the introduction of the cross-correlated cost functional has higher influence on a TE case than a TM case, which is demonstrated by the following simulation results. It is worth noting that in the formulation of TE scattering problems, the operators involved have the same form, but one spatial dimension lower, compared to full 3-D scattering problems. Hence, the performance gain with CC-CSI in future 3-D inversion problems can be compared to the performance gain in the TE case.

4.3 Noise-free data

For fair comparison, in the following simulations, the contrast sources are initialized by Eq. (21) and the contrast is initialized by Eq. (22) for all the three algorithms. Since the background of this benchmark problem is free space, viz., $\Re\{\bm{\chi}\}\geq 0$ and $\Im\{\bm{\chi}\}\leq 0$, we exploit this a priori information by simply enforcing the negative real part and the positive imaginary part of the contrast to zero after each update of the contrast [26]. Let us first investigate the inversion performance to the noise-free data. Both the TM-polarized data and the TE-polarized data are processed by classical CSI, MR-CSI, and CC-CSI, respectively. The relative permittivity and conductivity of the reconstructed contrast after 2048 iterations are shown in Fig. 3 for the TM case and Fig. 5 for the TE case. From Fig. 3 and Fig. 5 we see that MR-CSI generates blocky images because of the introduction of the total variation constraint, while classical CSI and CC-CSI have obvious variation in the reconstructed images. As we can see CC-CSI show better robustness by noting from Fig. 3 that the images of the contrast $\chi=2.0-0.6\text{i}$ obtained by classical CSI and MR-CSI show more distortion than those of CC-CSI, and that classical CSI and MR-CSI fail to reconstruct the contrast of $\chi=2.5-0.6\text{i}$ using the TM-polarized data. We can also see From Fig. 5 that MR-CSI fail to reconstruct the contrast $\chi=2.0-0.6\text{i}$ using the noise-free TE-polarized data. In addition, we see that the interior hollow tube is better reconstructed by CC-CSI in the TE case, and the two smaller tubes are better distinguished by CC-CSI compared to classical CSI and MR-CSI.

To quantitatively investigate the reconstruction accuracy, the reconstruction error of the three inversion methods is defined in the following as

[TABLE]

Fig. 4 and Fig. 6 give the comparison of the reconstruction error curves in terms of the iteration number of the three methods in the TM case and the TE case, respectively. As we can see quantitatively that, in the TM case, the three methods reach the same reconstruction errors in reconstructing the contrasts $\chi=1.0-0.6\text{i}$, $1.5-0.6\text{i}$, and $2.0-0.6\text{i}$. However, the reconstruction errors of classical CSI and MR-CSI do not decrease in reconstructing the contrast $\chi=2.5-0.6\text{i}$. In the contrast, the decreasing tendency of the reconstruction error curve of CC-CSI is not obviously affected by increasing the value of the contrast. In the TE case, we can see from Fig. 6 that CC-CSI can achieve lower reconstruction errors than classical CSI and MR-CSI, indicating the higher inversion accuracy of CC-CSI compared to classical CSI and MR-CSI. The reconstruction error of MR-CSI in reconstructing the contrast $\chi=2.0-0.6\text{i}$ does not decrease, and the reconstruction error curve of classical CSI in reconstructing the contrast $\chi=2.5-0.6\text{i}$ shows non-monotonicity. This demonstrates the poor robustness of both classical CSI and MR-CSI, and the better robustness of CC-CSI.

It is worth noting that CC-CSI shows higher inversion accuracy than classical CSI and MR-CSI in the TE case, but similar inversion accuracy with classical CSI and MR-CSI in TM case. Recalling the previous subsection, we know that the matrix $\bm{\Phi}$ of TE polarization has a much larger condition number than that of TM polarization. Therefore, compared to the TM case, same level of cross-correlated error presented in the measurement domain in TE case corresponds to a smaller reconstruction error in the field domain.

In addition, as we can see from the simulation results, MR-CSI shows not only poor robustness, but also unstable performance with respect to inversion accuracy compared to classical CSI and CC-CSI. As is well known, total variation regularization was originally proposed for noise removal in the digital image processing [35]. Obviously, the feasibility condition for applying the total variation regularization is that this noisy image is suitable for processing. However, this is apparently not the case in CSI, because the image of the contrast in CSI is optimized iteratively. In the design of MR-CSI, total variation constraint is very likely to be applied to a seriously distorted image of the contrast in the beginning, and thus may mislead and degrade the optimization process. Therefore, benefits can be possibly obtained from MR-CSI only if the contrast can be reliably reconstructed with CSI. Namely, the benefits from MR-CSI are not guaranteed. This perfectly explains the instable performance of MR-CSI shown in Fig. 4 and Fig. 6.

4.4 Noise-disturbed data

In real applications, the measurement data are very likely to be disturbed by noises. Apart from that, there is always error in modeling the incident fields. In this subsection, we investigate the inversion performance of the three methods to noise-disturbed data while the incident fields are assumed to be exactly known. Random white noise is added to the measurement data following the same procedure used in [32],

[TABLE]

with $p=1,2,\cdots,P$, $m=1,2,\cdots,M$. Here, $\bm{n}_{1}$ and $\bm{n}_{2}$ are two random numbers varying from $-1$ up to 1, $\zeta$ = 10% is the amount of noise, and $\underset{m}{\max}\{|f_{p,m}|\}$ gives the largest value among the amplitudes of the $M$ measurement data, which means the noise is scaled by the largest amplitude of the measurement data. $M=36$ in the TM case and $M=72$ in the TE case.

Fig. 7 and Fig. 9 show the inverted results after 2048 iterations by the three methods using the TM-polarized data and the TE-polarized data, respectively. In comparison to Fig. 3 and Fig. 5, we can see obvious distortion in the reconstructed images because of the disturbance by the additive random noise. What’s worse, we see from Fig. 9 that the hollow tube is not distinguishable any more in the images obtained by MR-CSI using the noise-disturbed TE-polarized data. As mentioned previously in Subsection 4.3, since the hollow tube cannot be well recognized by CSI (see the images obtained by CSI in Fig. 9), we lose the basic feasibility condition for applying the total variation constraint and therefore the hollow tube is reconstructed by MR-CSI to a solid one. In the contrast, CC-CSI is still capable to distinguish the hollow tube even though the measurement data has been disturbed by 10% additive random noise.

The reconstruction error curves of classical CSI, MR-CSI and CC-CSI in the noise-disturbed cases are given in Fig. 8 for the TM case and Fig. 10 for the TE case. As we can see, the error curves of classical CSI and MR-CSI show instability as the contrast goes higher, while the error curves of CC-CSI are always monotonously decreasing, indicating again the better robustness of CC-CSI. We can also see from Fig. 10 that CC-CSI has lower reconstruction error than classical CSI and MR-CSI, which is consistent with the simulation results in the noise-free cases.

What’s more, from Fig. 8(a) we see that the reconstruction errors of classical CSI and CC-CSI turn out to increase as the iteration goes on. Same phenomenon also occurs in Fig. 10(a,b) for CC-CSI. By comparison to Fig. 4(a) and Fig. 6(a,b) in the noise-free case, this phenomenon can been easily explained by the introduction of the additive random noise. Therefore, in real applications, a good termination condition is critical for preventing the methods from over-fitting the noise and for saving computation time.

4.5 Computational performance

As mentioned previously, CC-CSI can be implemented without significantly increasing the computational complexity. To demonstrate this point, we ran the MATLAB codes on a desktop with one Intel(R) Core(TM) i5-3470 CPU @ 3.20 GHz, and we did not use parallel computing. The computation times of classical CSI, MR-CSI, and CC-CSI, running for 2048 iterations are given in Table 1.

As we can see CSI is the most efficient, MR-CSI is in the middle, and CC-CSI runs slightly longer. If we define the increment percentage of the running times of CC-CSI as

[TABLE]

we have that the running time of CC-CSI is slightly longer than that of classical CSI by around 5.0% in the TM case and 7.3% in the TE case.

5 Conclusion

In this paper, a cross-correlated contrast source inversion (CC-CSI) method is proposed by modifying the cost functional of the CSI method to interrelate the state error and the data error. The proposed algorithm is tested with a 2-D benchmark problem which has also been tested by Belkebir and Tijhuis [30], Litman et al. [31], and Berg et al. [26, 32]. The simulation results with both TM-polarized wave and TE-polarized wave show that CC-CSI outperforms classical CSI and MR-CSI with respect to robustness and inversion accuracy, especially in the TE case. Which shows to be promising for the robustness and inversion accuracy in full 3-D inversion problems. We have also shown that CC-CSI can be implemented without significantly increasing the computational burden. As the Maxwell equations are formulated within a 3-D finite difference frequency domain (FDFD) scheme, it is straightforward to extend the proposed inversion scheme to future 3-D inverse scattering problems. Numerical results of 3-D scattering objects, including the application of the proposed method to experimental data will be published in future work.

Appendix A Derivation of the step size $\alpha_{p,n}$

First, let us rewrite the cost function $C_{\text{CC-CSI}}^{\bm{j}}(\bm{\chi}_{n-1},\bm{j}_{p,n-1}+\alpha_{p}\bm{\nu}_{p,n})$ as follows

[TABLE]

Obviously, it can be further simplified in the form of

[TABLE]

Therefore, we have

[TABLE]

Note that

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

it is easy to obtain that

[TABLE]

where, $a_{p,2}$, $b_{p,2}$, and $c_{p,2}$ are given by Eq. (32a), Eq. (33a), and Eq. (34a), respectively.