A Linear Method for Shape Reconstruction based on the Generalized Multiple Measurement Vectors Model

TL;DR

This paper introduces a linear shape reconstruction method using the generalized multiple measurement vectors model, leveraging joint sparsity and cross validation to improve accuracy and outperform traditional methods in focusing performance.

Contribution

It presents a novel linear reconstruction approach based on GMMV, incorporating joint sparsity and CV to enhance shape estimation without noise level prior.

Findings

Outperforms linear sampling method in focusing performance

Uses cross validation to avoid noise level estimation

Validates with TM experimental data

Abstract

In this paper, a novel linear method for shape reconstruction is proposed based on the generalized multiple measurement vectors (GMMV) model. Finite difference frequency domain (FDFD) is applied to discretized Maxwell's equations, and the contrast sources are solved iteratively by exploiting the joint sparsity as a regularized constraint. Cross validation (CV) technique is used to terminate the iterations, such that the required estimation of the noise level is circumvented. The validity is demonstrated with an excitation of transverse magnetic (TM) experimental data, and it is observed that, in the aspect of focusing performance, the GMMV-based linear method outperforms the extensively used linear sampling method (LSM).

Figures

Tables

Table 1. Table 1: Running times of the experimental examples

Data-set	Frequency number	GMMV /s	LSM /s
twodielTM_8f	1	2.5	0.0145
twodielTM_8f	4	12.7	0.0270
FoamDieIntTM	5	3.8	0.0693
rectTM_dece	4	2.7	0.0225
uTM_shaped	4	41.0	0.0498
FoamMetExtTM	7	15.6	0.0911

Equations

$$-{\nabla }^2{E}_{p,i}^{\text{sct}}-k_i^2{E}_{p,i}^{\text{sct}}={\omega }_i^2{\mu }_0{J}_{p,i},p=1,2,\dots ,P,i=1,2,\dots ,I,$$

$${\mathbf{A}}_i{\mathbf{e}}_{p,i}^{\text{sct}}={\omega }_i^2{\mathbf{j}}_{p,i},p=1,2,\dots ,P,i=1,2,\dots ,I,$$

$${\mathbf{y}}_{p,i}={\mathbf{\Phi }}_{p,i}{\mathbf{j}}_{p,i}^{\text{ic}},p=1,2,\dots ,P,i=1,2,\dots ,I,$$

$$\mathbf{Y}=\mathbf{\Phi }\cdot \mathbf{J}+\mathbf{U}$$

$$\mathbf{Y}=\left[\begin{array}{ccccccc}{\mathbf{y}}_{1,1}& {\mathbf{y}}_{2,1}& \dots & {\mathbf{y}}_{P,1}& {\mathbf{y}}_{1,2}& \dots & {\mathbf{y}}_{P,I}\end{array}\right],$$

$$\mathbf{J}=\left[\begin{array}{ccccccc}{\mathbf{j}}_{1,1}^{\text{ic}}& {\mathbf{j}}_{2,1}^{\text{ic}}& \dots & {\mathbf{j}}_{P,1}^{\text{ic}}& {\mathbf{j}}_{1,2}^{\text{ic}}& \dots & {\mathbf{j}}_{P,I}^{\text{ic}}\end{array}\right],$$

$$\mathbf{\Phi }\cdot \mathbf{J}=\left[\begin{array}{cccc}{\mathbf{\Phi }}_{1,1}{\mathbf{j}}_{1,1}^{\text{ic}}& {\mathbf{\Phi }}_{2,1}{\mathbf{j}}_{2,1}^{\text{ic}}& \dots & {\mathbf{\Phi }}_{P,I}{\mathbf{j}}_{P,I}^{\text{ic}}\end{array}\right],$$

$${\mathbf{\Phi }}^H\cdot \mathbf{Y}=\left[\begin{array}{cccc}{\mathbf{\Phi }}_{1,1}^H{\mathbf{y}}_{1,1}^{\text{ic}}& {\mathbf{\Phi }}_{2,1}^H{\mathbf{y}}_{2,1}^{\text{ic}}& \dots & {\mathbf{\Phi }}_{P,I}^H{\mathbf{y}}_{P,I}^{\text{ic}}\end{array}\right].$$

$$|\text{supp}(\mathbf{X})|<\frac{\text{spark}(\mathbf{A})-1+\text{rank}(\mathbf{X})}{2},$$

$$\text{minimize}\kappa (\mathbf{J})\text{subject to}\Vert \mathbf{\Phi }\cdot \mathbf{J}-\mathbf{Y}{\Vert }_F\le \tilde{\sigma },$$

$$\Vert \mathbf{J}{\Vert }_{\alpha ,\beta }:={\left(\sum _{n=1}^N{\Vert {\mathbf{J}}_{n,:}^T\Vert }_{\beta }^{\alpha }\right)}^{1/\alpha },$$

$$\text{minimize}{\Vert \mathbf{\Phi }\cdot \mathbf{J}-\mathbf{Y}\Vert }_F\text{subject to}{\Vert \mathbf{J}\Vert }_{1,2}\le \tau .$$

$${\varphi }_{\text{GMMV}}(\tau )={\Vert \mathbf{\Phi }\cdot {\mathbf{J}}_{\tau }^{\text{opt}}-\mathbf{Y}\Vert }_F,$$

$${\tau }_{l+1}={\tau }_l+\frac{\tilde{\sigma }-{\varphi }_{\text{GMMV}}({\tau }_l)}{{\varphi }_{\text{GMMV}}^{\prime }({\tau }_l)},$$

$${\varphi }_{\text{GMMV}}^{\prime }({\tau }_l)=-\frac{{\Vert {\mathbf{\Phi }}^H\cdot (\mathbf{\Phi }\cdot {\mathbf{J}}_{{\tau }_l}^{\text{opt}}-\mathbf{Y})\Vert }_{\infty ,2}}{{\Vert \mathbf{\Phi }\cdot {\mathbf{J}}_{{\tau }_l}^{\text{opt}}-\mathbf{Y}\Vert }_F}.$$

$${\tau }_1=\frac{\tilde{\sigma }-\Vert \mathbf{Y}{\Vert }_F}{\Vert {\mathbf{\Phi }}^H\cdot \mathbf{Y}{\Vert }_{\infty ,2}}.$$

$$r_{\text{rec}}={\left(\sum _{i=1}^I\sum _{p=1}^P{\Vert {\mathbf{y}}_{p,i,r}-{\mathbf{\Phi }}_{p,i,r}{\mathbf{j}}_{p,i}\Vert }_2^2\right)}^{1/2}$$

$$r_{\text{CV}}={\left(\sum _{i=1}^I\sum _{p=1}^P{\Vert {\mathbf{y}}_{p,i,\text{CV}}-{\mathbf{\Phi }}_{p,i,\text{CV}}{\mathbf{j}}_{p,i}\Vert }_2^2\right)}^{1/2},$$

$$N_{\text{Iter}}>N_{\text{opt}}+\Delta N,$$

$${\mathbf{\gamma }}_{\text{GMMV},n}=\sum _{i=1}^I\sum _{p=1}^P{\left|{\mathbf{j}}_{p,i,n}^{\text{ic}}\right|}^2,n=1,2,\dots ,N,$$

$$\Delta \le \frac{\min \{{\lambda }_i\}}{15},i=1,2,\dots ,I,$$

$$\int E_i({\mathbf{x}}_r,{\mathbf{x}}_t)g_i({\mathbf{x}}_s,{\mathbf{x}}_t)d{\mathbf{x}}_t=\frac{{\omega }_i{\mu }_0}{4}{H}_0^{(1)}(-k_i\Vert {\mathbf{x}}_s-{\mathbf{x}}_r{\Vert }_2),$$

$${\mathbf{F}}_i{\mathbf{g}}_{i,{\mathbf{x}}_s}={\mathbf{f}}_{i,{\mathbf{x}}_s},i=1,2,\dots ,I,$$

$$\Vert {\mathbf{g}}_{i,{\mathbf{x}}_s}{\Vert }^2=\sum _{d=1}^D{\left(\frac{s_{i,d}}{s_{i,d}^2+a_i^2}\right)}^2{\left|{\mathbf{u}}_d^H{\mathbf{f}}_{i,{\mathbf{x}}_s}\right|}^2,$$

$${\mathbf{\gamma }}_{\text{LSM}}({\mathbf{x}}_s)=\frac{1}{\Vert {\mathbf{g}}_{{\mathbf{x}}_s}^{\text{MF}}{\Vert }^2},$$

$${\Vert {\mathbf{g}}_{{\mathbf{x}}_s}^{\text{MF}}\Vert }^2=\frac{1}{I}\sum _{i=1}^I\frac{{\Vert {\mathbf{g}}_{i,{\mathbf{x}}_s}\Vert }^2}{\underset{{\mathbf{x}}_s\in \mathcal{D}}{\max }\left(\Vert {\mathbf{g}}_{i,{\mathbf{x}}_s}{\Vert }^2\right)}.$$

$${\mathbf{\gamma }}_{\text{dB}}=10\times {\log }_{10}\left(\frac{\mathbf{\gamma }}{\max \{\mathbf{\gamma }\}}\right).$$

Full text

1 Introduction

Inverse electromagnetic (EM) scattering is a procedure of recovering the characteristics of unknown objects using the scattered fields probed at a number of positions. In many real applications, such as geophysical survey [1, 2, 3, 4], it is of great importance to retrieve the geometrical features of a system of unknown targets.

For solving this problem, a wealth of reconstruction methods have been proposed over the recent decades. Due to their high efficiency, the linear focusing methods have been extensively used in real applications, among which are Kirchhoff migration [5], back-projection method [6], time-reversal (TR) technique [7, 8, 9, 10, 11, 12, 13, 14, 15], and so forth. However, as is well known the imaging resolutions of the linear focusing algorithms are bound by the diffraction limit [16]. As a variant of TR technique, time-reversal multiple signal classification (TR MUSIC) [17, 18, 19, 20] is capable to achieve a resolution that can be much finer than the diffraction limit by exploiting the orthogonality of the signal and noise sub-spaces. Linear Sampling Method (LSM) [21, 22] is a non-iterative inversion technique of finding an indicator function for each position in the region of interest (ROI) by first defining a far-field (or near-field [23]) mapping operator, and then solving a linear system of equations. LSM has been proven to be effective for impenetrable scatterers, and in some cases, also applicable for dielectric scatterers [24]. As a matter of fact, LSM can also be reinterpreted, apart from very peculiar cases, as a “synthetic focusing” problem [25], and more interestingly, an extension of the MUSIC algorithm [26]. There is another group of iterative surface-based inversion methods, which first parametrize the shape of the scatterer, then optimize the parameters by minimizing a cost functional iteratively [27]. The drawbacks of these methods are obvious. Firstly, they require a priori information about the position and the quantity of the scatterers. More research on this point can be found in [28, 29]. Secondly, it is intractable to deal with complicated non-convex objects. Quantitative inversion methods, such as contrast source inversion (CSI) [30, 31, 32, 33] and (Distorted) Born iterative methods (BIM and DBIM) [34, 35, 36, 37], can also be used for shape reconstruction. However, it is very time consuming due to the fact that the forward scattering problem needs to be solved in every iteration.

In this paper, a novel linear method using generalized multiple measurement vectors (GMMV) model [38, 39] is proposed for solving the problem of shape reconstruction. Specifically, as the objects are illuminated by EM waves from various incident angles at different frequencies, the contrast sources, i.e., the multiplication of the contrast and the total fields, are distributed in the same region with the objects. Therefore, the problem is consequently formulated as a GMMV model, and the contrast sources can be retrieved by solving multiple systems of linear equations simultaneously. In our method, the sum-of-norm of the contrast sources is used as a regularization constraint to address the ill-posedness. Finite difference frequency domain (FDFD) [40] is used to construct the scattering operator which enables simple incorporation of complicated background media, and the spectral projected gradient method, referred to as SPGL1 [41, 42], is selected to estimate the contrast sources by solving a sum-of-norm minimization problem. Sparse scatterer imaging has been studied in [43], in which the single measurement vector (SMV) model was used, but the joint sparsity was not considered. The application of joint sparsity in the field of medical imaging has been reported in [44], which is actually a hybridization of compressive sensing (CS) [45] and MUSIC based on a so-called generalized MUSIC criterion. In the aforementioned work, sparse targets (original or equivalently transformed) and their sparsest solutions are considered, and the problem of defining the best discretization grid and target number is critical for ensuring a level of sparsity that is recoverable. Equivalence principles have been considered in [46] for reconstructing the boundary of dielectric scatterers. In this paper, we use sum-of-norm as a regularization constraint and we demonstrate a regularized solution of the contrast sources is sufficient to recover the spatial profile of the non-sparse targets. In this paper, we only considered the transverse magnetic (TM) EM scattering problem, and we verified the validity of the proposed method with 2-D experimental data provided by the Institut Fresnel, France [47, 48] for three distinct cases – metallic objects, dielectric objects, and a hybrid one of both. Since the noise level is unknown in real applications, cross validation (CV) technique [49] is used to terminate the optimization process. Comparison of the inverted results indicates that the proposed method shows higher resolving ability than LSM.

The remainder of the paper is organized as follows: In Section 2, the problem statement is given. In Section 3, the proposed GMMV-based linear method111The GMMV-LIM package is available at https://github.com/TUDsun/GMMV-LIM. is introduced in detail. The validation of this method with experimental data is given in Section 4. Finally, Section 5 ends this paper with our conclusions.

2 Problem Statement

For the sake of simplicity, we consider the 2-D TM-polarized EM scattering problem. A bounded, simply connected, inhomogeneous background domain $\mathcal{D}$ contains unknown objects. The domain $\mathcal{S}$ contains the sources and receivers. The sources are denoted by the subscript $p$ in which $p\in\{1,2,3...,P\}$, and the receivers are denoted by the subscript $q$ in which $q\in\{1,2,3,...,Q\}$. We use a right-handed coordinate system, and the unit vector in the invariant direction points out of the paper.

Assume the background is known to a reasonable accuracy beforehand, and the permeability of the background and unknown objects is a constant, denoted by $\mu_{0}$. The contrast corresponding to the $i$-th frequency, $\chi_{i}$, is defined as $\chi_{i}=\epsilon_{i}-\epsilon^{\text{bg}}_{i}$, where $\epsilon_{i}=\varepsilon-\text{i}\sigma/\omega_{i}$ and $\epsilon^{\text{bg}}_{i}=\varepsilon^{\text{bg}}-\text{i}\sigma^{\text{bg}}/\omega_{i}$ are the complex permittivity of the inversion domain with and without the presence of the targets, respectively. Here, $\varepsilon$ and $\varepsilon^{\text{bg}}$ are the permittivity of the inversion domain with and without the presence of the targets, respectively; $\sigma$ and $\sigma^{\text{bg}}$ are the conductivity of the inversion domain with and without the presence of the targets, respectively; $\omega_{i}$ is the $i$-th angular frequency; i represents the imaginary unit. The time factor used in this paper is $\exp(\text{i}\omega_{i}t)$. For 2-D TM-polarized scattering problems, the electric field is a scalar and the scattering wave equation with respect to the scattered fields can be easily derived from Maxwell’s equations, which is given by

[TABLE]

where, $\nabla^{2}$ is the Laplace operator, $k_{i}=\omega_{i}\sqrt{\epsilon_{b}\mu_{0}}$ is the $i$-th wavenumber, $J_{p,i}=\chi_{i}E_{p,i}^{\text{tot}}$ is the $p$-th contrast source at the $i$-th frequency, $E_{p,i}^{\text{sct}}$ and $E_{p,i}^{\text{tot}}$ are the scattered electric field and the total electric field at the $i$-th frequency, respectively. The inverse scattering problems discussed in this paper are to retrieve the geometrical features of the scatterers from a set of scattered field measurements.

3 The GMMV-based Linear Method

3.1 The GMMV Formulation

Following the vector form of the FDFD scheme in [40], we discretize the 2-D inversion space with $N$ grids and recast the scattering wave equation (1) into the following matrix formalism

[TABLE]

where $\bm{A}_{i}\in\mathbb{C}^{N\times N}$ is the FDFD stiffness matrix of the $i$-th frequency, which is highly sparse; $\bm{e}_{p,i}^{\text{sct}}\in\mathbb{C}^{N}$ and $\bm{j}_{p,i}\in\mathbb{C}^{N}$ are the scattered electric field and the contrast source in the form of a column vector, respectively. Obviously, the solution to Eq. (2) can be obtained by inverting the stiffness matrix $\bm{A}_{i}$, i.e., $\bm{e}_{p,i}^{\text{sct}}=\bm{A}_{i}^{-1}\omega_{i}^{2}\bm{j}_{p,i}$. For the inverse scattering problems discussed in this paper, the scattered fields are measured with a number of receivers at specified positions, yielding the data equations given by

[TABLE]

where, $\bm{\Phi}_{p,i}=\mathcal{M}^{\mathcal{S}}_{p}\bm{A}_{i}^{-1}\omega_{i}\in\mathbb{C}^{Q\times N}$ is the sensing matrix for the measurement $\bm{y}_{p,i}$, $\bm{j}^{\text{ic}}_{p,i}=\omega_{i}\bm{j}_{p,i}$ is the normalized contrast source proportional to the induced current $\text{i}\omega_{i}\mu_{0}\bm{j}_{p,i}$. Here, $\mathcal{M}^{\mathcal{S}}_{p}$ is a measurement matrix selecting the values of the $p$-th scattered field at the positions of the receivers.

In the rest of this subsection, a GMMV model [39] is constructed and solved by exploiting the joint sparsity of the normalized contrast sources. In doing so, the contrast sources can be well estimated by solving a sum-of-norm minimization problem, and consequently be used to indicate the shape of the scatterers. To do so, we reformulate the data equations, Eq. (3), as

[TABLE]

where

[TABLE]

[TABLE]

and $\bm{\Phi}\cdot\bm{J}$ is defined by

[TABLE]

and, correspondingly, $\bm{\Phi}^{H}\cdot\bm{Y}$ is defined as

[TABLE]

Here, $\bm{Y}\in\mathbb{C}^{Q\times PI}$ is the measurement data matrix, and the columns of $\bm{J}\in\mathbb{C}^{N\times PI}$ are the multiple vectors to be solved. $\bm{U}\in\mathbb{C}^{Q\times PI}$ represents the complex additive noises satisfying certain probability distribution. It is worth noting that for single frequency inverse scattering problem, if the positions of the receivers are fixed, i.e., $\bm{\Phi}_{1,1}=\bm{\Phi}_{2,1}=\cdots=\bm{\Phi}_{Q,1}$, Eq. (4) reduces to the standard multiple measurement vectors (MMV) model [38].

3.2 Guideline of the Measurement Configuration

Although the joint sparsity is used in this paper as a regularization constraint, an investigation on the uniqueness condition is still of much importance for two reasons: 1) It is of great interest to know how much we could benefit from a joint recovery; 2) It provides us a guideline of the measurement configuration to make the most of the joint processing.

According to the work of Chen and Huo [50] and Davies and Eldar [51], a necessary and sufficient condition for the measurements $\bm{Y}=\bm{A}\bm{X}$ to uniquely determine the row sparse matrix $\bm{X}$ is that

[TABLE]

where, $\text{supp}(\bm{X})$ denotes the index set corresponding to non-zero rows of matrix $\bm{X}$, $|\text{supp}(\bm{X})|$ denotes the cardinality of $\text{supp}(\bm{X})$, the spark of a given matrix is defined as the smallest number of the columns that are linearly dependent. Thereafter, Heckel and Bölcskei have studied the GMMV problem and showed that having different measurement matrices can lead to performance improvement over the standard MMV case [39]. The above work about the uniqueness condition implies specifically in our method that in order to make the most of the joint processing, the column number of matrix $\bm{J}$ is supposed to be larger than the number of measurements, i.e., $P\times I>Q$. Moreover, with the same measurement configuration, the inversion performance can be further improved by exploiting the frequency diversity even for the case of $P>Q$. The latter is further demonstrated in Subsection 4.1.

3.3 Spectral Projected Gradient L1 method (SPGL1)

3.3.1 GMMV basis pursuit denoise (BPσ) problem

Suppose the noise level is known beforehand, the approach to finding the multiple vectors is based on solving a convex optimization problem (referred to as GMMV (BPσ) problem), which can be written as follows

[TABLE]

where, $\tilde{\sigma}$ represents the noise level; $\kappa(\bm{J})$ is the mixed ($\alpha,\beta$)-norm defined as

[TABLE]

where, $\bm{J}_{n,:}$ denotes the $n$-th row of $\bm{J}$; $\|\cdot\|_{\beta}$ is the conventional $\beta$-norm; $(\cdot)^{T}$ is the transpose operator; $\|\cdot\|_{F}$ is the Frobenius norm which is equivalent to the mixed (2,2)-norm $\|\cdot\|_{2,2}$. In this paper, we select the mixed norm $\|\cdot\|_{1,2}$ as a regularized constraint. Although $\|\cdot\|_{1,2}$ tends to enforce the row-sparsity of the matrix $\bm{J}$, sparsity is not a premise for this approach. The key point is the utilization of the joint structure for improving the focusing ability. As demonstrated in the following experiments, this approach is able to image objects which are not sparse by exploitation of the frequency diversity.

3.3.2 Multiple GMMV Lasso (LSτ) problems

Since it is not straightforward to solve the GMMV (BPσ) problem, we consider the GMMV (LSτ) problem formulated as [41]

[TABLE]

The GMMV (LSτ) problem is equivalent to the GMMV (BPσ) problem when $\tau=\tau_{\tilde{\sigma}}$. As the exact value of $\tau_{\tilde{\sigma}}$ is not available, a series of GMMV (LSτ) problems with different values of $\tau$ must be solved. Now let us first introduce the Pareto curve defined as follows

[TABLE]

where, $\bm{J}_{\tau}^{\text{opt}}$ is the optimal solution to the LSτ problem given by Eq. (12). When the optimal solution $\bm{J}_{\tau_{l}}^{\text{opt}}$ to the GMMV (LSτ) problem is found, $\tau_{l}$ is updated to $\tau_{l+1}$ by probing the Pareto curve. The searching procedure is terminated when $\phi_{\text{GMMV}}(\tau)=\tilde{\sigma}$. At the mean time, $\tau$ reaches $\tau_{\tilde{\sigma}}$.

3.3.3 Updating the parameter $\tau$

As the Pareto curve is proven to be a non-increasing convex function, Newton iteration is used for updating the parameter $\tau$.

Specifically, $\tau$ is updated by

[TABLE]

where,

[TABLE]

Here, $\|\cdot\|_{\infty,2}$ is the dual norm of $\|\cdot\|_{1,2}$. The searching procedure is illustrated in Fig. 1. Unless a good estimate of $\tau_{\tilde{\sigma}}$ can be obtained, we set $\tau_{\tilde{\sigma}}=0$, yielding $\phi(0)=\|\bm{Y}\|_{F}$ and $\phi^{\prime}(0)=\|\bm{\Phi}^{H}\cdot\bm{Y}\|_{\infty,2}$. With Eq. (14), it holds immediately that the next Newton iteration is

[TABLE]

We refer to [41, 42] for more details about SPGL1 and [52] for its application in inverse scattering problems.

3.4 CV-based Modified SPGL1

In real applications, the termination condition, $\phi_{\text{GMMV}}(\tau)=\tilde{\sigma}$, is not applicable, because the noise level, i.e., the parameter $\tilde{\sigma}$, is unknown in general. In order to deal with this problem, we modified the SPGL1 method based on the CV technique [49]. In doing so, the estimation of the noise level can be well circumvented.

Specifically, we separate the original scattering matrix to a reconstruction matrix $\bm{\Phi}_{p,i,r}\in\mathbb{C}^{Q_{r}\times N}$ and a CV matrix $\bm{\Phi}_{p,i,\text{CV}}\in\mathbb{C}^{Q_{\text{CV}}\times N}$ with $Q=Q_{r}+Q_{\text{CV}}$. The measurement vector $\bm{y}_{p,i}$ is also separated accordingly to a reconstruction measurement vector $\bm{y}_{p,i,r}\in\mathbb{C}^{Q_{r}}$ and a CV measurement vector $\bm{y}_{p,i,\text{CV}}\in\mathbb{C}^{Q_{\text{CV}}}$. The reconstruction residual and the CV residual are defined as

[TABLE]

and

[TABLE]

respectively. In doing so, every iteration can be viewed as two separate parts: reconstructing the contrast sources by SPGL1 and evaluating the outcome by the CV technique. The CV residual curve turns to increasing when the reconstructed signal starts to overfit the noise. The reconstructed contrast sources are selected as the output on the criterion that its CV residual is the least one. To find the least CV residual, we initialize $\tilde{\sigma}$ as 0 and terminate the iteration when

[TABLE]

is satisfied, Here, $N_{\text{Iter}}$ is the current iteration number, and $N_{\text{opt}}$ is the iteration index corresponding to the least CV residual — the optimal solution. Namely, the CV residual is identified as the least one if the CV residual keeps increasing monotonously for $\Delta N$ iterations. In the following experimental examples, we set $\Delta N=30$.

Once the normalized contrast sources are obtained, one can achieve the shape of the scatterers defined as

[TABLE]

where $\bm{j}^{\text{ic}}_{p,i,n}$ and $\bm{\gamma}_{\text{GMMV},n}$ represent the $n$-th element of vector $\bm{j}^{\text{ic}}_{p,i}$ and $\bm{\gamma}_{\text{GMMV}}$, respectively.

4 Validation with Experimental Data

In order to validate the proposed GMMV-based linear method, we applied it to the experimental database provided by the Remote Sensing and Microwave Experiments Team at the Institut Fresnel, France, in the years of 2001 [47] and 2005 [48]. Three different cases of dielectric scatterers, metallic scatterers (convex and nonconvex), and a hybrid one of both, were considered, respectively. In order to guarantee the accuracy of the FDFD scheme, the inversion domain is discretized with a grid size $\Delta^{2}$ satisfying

[TABLE]

where, $\lambda_{i}$ is the wavelength of the $i$-th frequency.

We have also processed the same data by LSM for comparison. Since the background of the experiments is free space and only TM wave is considered, the LSM method consists in solving the integral equation of the indicator function $g_{i}(\bm{x}_{s},\bm{x}_{t})$ at the $i$-th frequency

[TABLE]

where, $E_{i}(\bm{x}_{r},\bm{x}_{t})$ is the scattered field probed at $\bm{x}_{r}$ corresponding to the transmitter at $\bm{x}_{t}$ and the $i$-th frequency. Here, $\bm{x}_{s}$ is the sampling point in the inversion domain, $H_{0}^{(1)}(\cdot)$ is the Hankel function of the first kind, $k_{i}$ is the wavenumber of the $i$-th frequency. Eq. (22) can be reformulated as a set of systems of linear equations

[TABLE]

where, $\bm{F}_{i}$ is the measurement data matrix, $\bm{g}_{i,\bm{x}_{s}}$ is the indicator function of the sampling point $\bm{x}_{s}$ in the form of a column vector, $\bm{f}_{i,\bm{x}_{s}}$ is the right side of Eq. (22) in the form of a column vector, the index $i$ represents the $i$-th frequency. Following the same approach of solving Eq. (23) in [25, 53], the indicator function $\bm{g}_{i,\bm{x}_{s}}$ is sought to be

[TABLE]

where, $s_{i,d}$ represents the singular value of matrix $\bm{F}_{i}$ corresponding to the singular vector $\bm{u}_{d}$; $(\cdot)^{H}$ is the conjugate transpose operator; $D=\min\{P,Q\}$; $a_{i}=0.01\times\underset{d}{\max}\{s_{i,d}\}$. The shape of the scatterers is defined by

[TABLE]

where, $\|\bm{g}^{\text{MF}}_{\bm{x}_{s}}\|^{2}$ is a multi-frequency modified indicator defined as the average of the normalized modified ones computed at each frequency [54]

[TABLE]

It is worth mentioning that both the normalized contrast sources and the indicator functions are proportional to the amplitude of the electric field. According to the definitions in Eq. (20) and Eq. (25), $\bm{\gamma}_{\text{GMMV}}$ and $\bm{\gamma}_{\text{LSM}}$ are proportional and inversely proportional to the power of electric fields, respectively. Therefore, the dB scaling shown in the following examples is defined as follows

[TABLE]

4.1 Dielectric Scatterers

4.1.1 Example 1

In the first example, we consider the twodielTM_8f data-set provided in the first opus of the Institut Fresnel’s database [47]. The targets consist of two identical circular cylinders, which are shown in Fig. 4a (a). All the cylinders have radius of 1.5 cm and relative permittivity $3\pm 0.3$. The emitter is placed at a fixed position on the circular rail, while a receiver is rotating around the center point of the vertical cylindrical target. The targets rotated from 0*∘* to 350*∘* in steps of 10*∘* with a radius of $720\pm 3$ mm, and the receiver rotated from 60*∘* to 300*∘* in steps of 5*∘* with a radius of $760\pm 3$ mm. Namely, we have $49$ $\times$ $36$ measurement data at each frequency when all the measurements are finished. The measurement configuration is shown in Fig. 2, from which we can see 9 red circles which represents the CV measurements and 40 green ones which represents the reconstruction measurements. The inversion domain is restricted to [$-75$, 75] $\times$ [$-75$, 75] mm2, and the size of the discretization grids is 2.5 $\times$ 2.5 mm2.

Let us first process the single frequency data at 4 GHz by the GMMV-based linear method and the LSM method. The data matrix $\bm{F}_{i}$ for the LSM is a 72 $\times$ 36 matrix in which the data entries that are not available are replaced with zeros. The reconstruction residual curve and the CV residual curve are shown in Fig. 3 (a), from which we see the CV residual decreases before the 52-nd iteration and starts to increase thereafter. The solutions at the turning point correspond to the optimal ones. In addition, the reconstruction residual corresponding to the turning point gives an estimation of the noise level $\tilde{\sigma}\approx 0.05\|\bm{Y}\|_{F}$. Fig. 4a (b) and Fig. 4a (c) show the images achieved by the two methods at 4 GHz in a dynamic range of [$-25$, [math]] dB. As we can see the GMMV image is more clear than the LSM image. However, there is obvious shape distortion in the former. Note that $Q=49$, $P=36$ and $I=1$, we have $P\times I<Q$. Recalling the guideline of the measurement configuration discussed in Subsection 3.2, the reconstruction performance can be further improved via exploiting the frequency diversity. It is worth mentioning that an obvious position mismatch of the true objects and the reconstructed result can be observed. The reason is very likely to be the minor displacement and tilt occurred in the placement of the objects while doing this measurement, because the same phenomenon can be observed as well in the inverted results reported in [55].

Now let us process the data at 2 GHz, 4 GHz, 6 GHz, and 8 GHz, simultaneously. The residual curves are shown in Fig. 3 (b) and the reconstructed images are shown in Fig. 4a (d) and Fig. 4a (e). By comparison of Fig. 4a (b) and Fig. 4a (d), one can see that the reconstruction performance of the proposed GMMV-based linear method is improved by exploiting the frequency diversity. One can also observe that the GMMV-based linear method achieves lower sidelobes than LSM in the case of dielectric scatterers.

4.1.2 Example 2

In the second example, we consider the FoamDieIntTM data-set provided in the second opus of the Institut Fresnel’s database. The targets consist of a circular dielectric cylinder with a diameter of 30 mm embedded in another circular dielectric cylinder with a diameter of 80 mm. The smaller cylinder has a relative permittivity value of $\varepsilon_{r}=3\pm 0.3$, while the larger cylinder has a relative permittivity value of $\varepsilon_{r}=1.45\pm 0.15$. Fig. 7a (a) shows the true objects, and we refer to [48] for more description of the targets. The experiment is carried out in 2005, in which the receiver stays in the azimuthal plane ($xoy$) and is rotated 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}$), and the object is rotated to obtain different illumination incidences from 0∘ to 315∘ with angular step of 45∘. Namely, we have $241\times 8$ measurements at each frequency. The distance from the transmitter/receiver to the centre of the targets has increased to 1.67 m. The measurement configuration is shown in Fig. 5, in which the blue one represents the emitter, the $4\times 9$ red ones represent the CV measurements, and the green ones are the reconstruction measurements.

The inversion domain is restricted to [$-60$, $60$] $\times$ [$-60$, $60$] mm2, and the discretization grid size is 2.5 $\times$ 2.5 mm2. Let us process the multi-frequency data at 2 GHz, 4 GHz, 6 GHz, 8 GHz, and 10 GHz simultaneously by the GMMV-based linear method and the LSM method, respectively. The data matrix $\bm{F}_{i}$ for LSM is a $360\times 8$ matrix in which the data entries that are not available are replaced with zeros. The reconstruction residual curve and the CV residual curve are shown in Fig. 6, from which we see the CV residual decreases during the first 62 iterations and starts to increase thereafter. Fig. 7a (b) and Fig. 7a (c) show the reconstructed images by the GMMV-based linear method and LSM, respectively. One can observe that the profile of the objects is reconstructed by the proposed method with high resolution, while in the LSM image the objects cannot be distinguished at all.

4.2 Metallic Scatterers

In this subsection, we applied the proposed method to the rectTM_dece and uTM_shaped data-sets provided in the first opus of the Institut Fresnel’s database [47], which correspond to a convex scatterer – a rectangular metallic cylinder, and a non-convex scatterer – a “U-shaped” metallic cylinder, respectively. The dimensions of the rectangular cross section are 24.5 $\times$ 12.7 mm2, while those of the “U-shaped” cylinder are about 80 $\times$ 50 mm2. The measurement configuration is same with that in Subsection 4.1. More details about the description of the targets can be found in [47].

For the rectangular cylinder, the inversion domain is restricted to [$-25$, 25] $\times$ [15, 65] mm2 and the multiple frequency data at 10 GHz, 12 GHz, 14 GHz, and 16 GHz are processed simultaneously. While for the larger “U-shaped” cylinder, the inversion domain is restricted to [$-70$, 70] $\times$ [$-70$, 70] mm2 and the multiple frequency data at 4 GHz, 8 GHz, 12 GHz, and 16 GHz are processed simultaneously. The size of the discretization grids is 1.3 $\times$ 1.3 mm2. Fig. 8 (a,b) give the residual curves and the reconstructed images are shown in Fig. 9a and Fig. 10a, respectively, from which we can see that the focusing performance of LSM is poor in the rectangular cylinder case, and it is even worse in retrieving the shape of the non-convex “U-shaped” cylinder, while the rectangular shape and “U” shape are well reconstructed by the proposed GMMV-based linear method, indicating that the latter shows higher resolving ability than the former in both the convex metallic target case and the non-convex metallic target case.

4.3 Hybrid Scatterers

In this subsection, we applied the proposed method to hybrid scatterers consisting of a foam circular cylinder (diameter = 80 mm, $\varepsilon_{r}=1.45\pm 0.15$) and a copper tube (diameter = 28.5 mm, thickness = 2 mm), which was considered in the FoamMetExtTM data-set provided in the second opus of the Institut Fresnel’s database. We refer to [48] for more description of the targets. The measurement configuration is the same with the one shown in Fig. 5. In contrast to the FoamDieIntTM data-set, this data-set is obtained using 18 transmitters, while other settings are kept the same. Specifically, the source antenna stays at the fixed location ($\theta=0^{\circ}$), and the object is rotated to obtain different illumination incidences from 0∘ to 340∘ in steps of 20∘.

Let us first restrict the inversion domain to [$-90$, 60] $\times$ [$-75$, 75] mm2 and discretize this domain with 2.5 $\times$ 2.5 mm2 grids. The multi-frequency data at 7 frequencies, 2 GHz, 3 GHz, $\dots$, and 8 GHz, are jointly processed. The data matrix $\bm{F}_{i}$ for LSM is a $360\times 18$ matrix in which the data entries that are not available are replaced with zeros. Fig. 11 gives the normalized residual curves of the GMMV-based linear method, and the reconstructed images by both methods are shown in Fig. 12a. As we can see both the metallic cube and the circular foam cylinder can be well reconstructed by the GMMV-based linear method with high resolution, but for a slight part lost in between. In addition, one can also see from the GMMV image that the metallic cube obviously has larger intensity than the foam cylinder, showing a potential ability of distinguishing dielectric objects and metallic objects. In contrast, LSM shows a poor focusing ability in the hybrid scatterer case, indicating once again that the proposed GMMV-based linear method is able to achieve higher resolution image than LSM in this case.

4.4 Computation Time

In this subsection, we discuss the computational complexity of the GMMV-based linear method. Since the sensing matrices can be computed (or analytically given for the experiments in homogeneous backgrounds) and stored beforehand, the GMMV-based linear method only involves a number of matrix-vector multiplications. The codes for reconstructing the contrast sources are written in MATLAB language. We ran the codes on a desktop with one Intel(R) Core(TM) i5-3470 CPU @ 3.20 GHz, and we did not use parallel computing. The running times of the GMMV-based linear method and LSM are listed in Table 1, from which we see that, on one hand, all the reconstructions by the GMMV-based linear method require less than 1 minute (or even a couple of seconds for some examples); on the other hand, LSM shows overwhelmingly high efficiency in comparison to the GMMV-based linear method, because singular value decomposition (SVD) in LSM is done only once, then all of the indicator functions can be obtained simultaneously by several matrix-matrix multiplications. However, in view of the higher resolving ability of the proposed method, the extra computational cost is worth to pay. It is also worth mentioning that the proposed method is faster than the iterative shape reconstruction methods which solve the forward scattering problem in each iteration. In addition, parallel computing can be straightforwardly applied to the proposed method for acceleration.

5 Conclusion

In this paper, a novel linear method for shape reconstruction based on the generalized multiple measurement vectors (GMMV) model is proposed. The sum-of-norm of the contrast sources at multiple frequencies was used for the first time as a regularization constraint in solving the electromagnetic inverse scattering problems. We applied this method to process 2-D transverse magnetic (TM) experimental data, and the results demonstrated that a regularized solution of the contrast sources by the sum-of-norm constraint is sufficient to recover the spatial profile of the non-sparse targets. Comparison results indicated that the GMMV-based linear method outperforms LSM in all the three cases of dielectric scatterers, convex and non-convex metallic scatterers, and hybrid scatterers in the shape reconstruction quality and the level of the sidelobes in the images. In view of the resolving ability and computational efficiency, the proposed method looks very promising in the application to three-dimensional imaging problems. Besides, the outcome of the GMMV-based linear method — the contrast sources, can be directly used for quantitative imaging when the incident fields are known with a reasonable accuracy.