A Linear Model for Microwave Imaging of Highly Conductive Scatterers

TL;DR

This paper introduces a linear, versatile model for microwave imaging of highly conductive objects, leveraging joint sparse regularization to improve resolution without prior information or model approximation.

Contribution

A novel linear model based on multiple measurement vectors for inverse scattering, utilizing joint sparse structure enforcement for enhanced imaging of conductive scatterers.

Findings

Higher resolving ability than linear sampling methods

Effective with synthetic and Fresnel data

No prior information or scattering model approximation needed

Abstract

In this paper, a linear model based on multiple measurement vectors model is proposed to formulate the inverse scattering problem of highly conductive objects at one single frequency. Considering the induced currents which are mostly distributed on the boundaries of the scatterers, joint sparse structure is enforced by a sum-of-norm regularization. Since no \textit{a priori} information is required and no approximation of the scattering model has been made, the proposed method is versatile. Imaging results with transverse magnetic and transverse electric polarized synthetic data and Fresnel data demonstrate its higher resolving ability than both linear sampling method and its improved version with higher, but acceptable, computational complexity.

Figures

Equations

$$\begin{array}[]{ccc}\begin{bmatrix}-\text{i}\omega\bm{\varepsilon}-\bm{\sigma}&\nabla\times\\ \nabla\times&\text{i}\omega\bm{\mu}\\ \end{bmatrix}&\begin{bmatrix}\bm{E}\\ \bm{H}\end{bmatrix}&=\begin{bmatrix}\bm{J}\\ -\bm{M}\end{bmatrix}\\ \end{array},$$

$$\nabla \times {\mathbf{\mu }}^{-1}\nabla \times \mathbf{E}-{\omega }^2\mathbf{ϵ}\mathbf{E}=-\text{i}\omega \mathbf{J}-\nabla \times {\mathbf{\mu }}^{-1}\mathbf{M},$$

$$\mathbf{ϵ}=\mathbf{\epsilon }-\text{i}\mathbf{\sigma }/\omega .$$

$${\mathbf{E}}_p^{\text{sct}}={\mathbf{E}}_p-{\mathbf{E}}_p^{\text{inc}},p=1,2,\dots ,P.$$

$$\nabla \times {\mu }_0^{-1}\nabla \times {\mathbf{E}}_p^{\text{sct}}-{\omega }^2{\mathbf{ϵ}}_{\text{b}}{\mathbf{E}}_p^{\text{sct}}={\omega }^2\mathbf{\chi }{\mathbf{E}}_p,p=1,2,\dots ,P,$$

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

$${\mathbf{e}}_p^{\text{sct}}=\left[\begin{array}{c}{\mathbf{e}}_{p,x_1}^{\text{sct}}\\ {\mathbf{e}}_{p,x_2}^{\text{sct}}\\ {\mathbf{e}}_{p,x_3}^{\text{sct}}\end{array}\right],{\mathbf{e}}_p=\left[\begin{array}{c}{\mathbf{e}}_{p,x_1}\\ {\mathbf{e}}_{p,x_2}\\ {\mathbf{e}}_{p,x_3}\end{array}\right],$$

$${\mathbf{y}}_p=\mathbf{\Phi }{\mathbf{j}}_p,p=1,2,3,\cdots ,P,$$

$$\mathbf{\Phi }={\mathcal{M}}_{\mathcal{S}}{\mathbf{A}}^{-1}{\omega }^2.$$

$${\mathbf{e}}_p={\mathbf{e}}_p^{\text{sct}}+{\mathbf{A}}^{-1}{\omega }^2{\mathbf{j}}_p,p=1,2,\dots ,P.$$

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

$$\mathbf{J}=\left[{\mathbf{j}}_1,{\mathbf{j}}_2,\cdots ,{\mathbf{j}}_P\right]$$

$$\mathbf{Y}=\left[{\mathbf{y}}_1,{\mathbf{y}}_2,\cdots ,{\mathbf{y}}_P\right]$$

$$({\text{SMV}}_{\text{TM}}{\text{BP}}_{\tilde{\sigma }}):\text{min}{\Vert {\mathbf{j}}_p\Vert }_1\text{s.t.}{\Vert \mathbf{\Phi }{\mathbf{j}}_p-{\mathbf{y}}_p\Vert }_2\le \tilde{\sigma }.$$

$$({\text{SMV}}_{\text{TM}}{\text{LS}}_{\tau }):\text{min}{\Vert \mathbf{\Phi }{\mathbf{j}}_p-{\mathbf{y}}_p\Vert }_2\text{s.t.}{\Vert {\mathbf{j}}_p\Vert }_1\le \tau .$$

$${\mathcal{P}}_{\tau }\left[{\mathbf{j}}_p\right]:=\left\{\underset{\mathbf{s}}{\mathrm{argmin}}{\Vert {\mathbf{j}}_p-\mathbf{s}\Vert }_2\text{s.t.}\Vert \mathbf{s}{\Vert }_1\le \tau \right\}.$$

$${\varphi }_{{\text{SMV}}_{\text{TM}}}(\tau )={\Vert \mathbf{\Phi }{\mathbf{j}}_{p,\tau }-{\mathbf{y}}_p\Vert }_2,$$

$${\tau }_{h+1}={\tau }_h+\frac{\tilde{\sigma }-{\varphi }_{{\text{SMV}}_{\text{TM}}}\left({\tau }_h\right)}{{\varphi }_{{\text{SMV}}_{\text{TM}}}^{\prime }\left({\tau }_h\right)},$$

$${\varphi }_{{\text{SMV}}_{\text{TM}}}^{\prime }({\tau }_{h+1})=-\frac{{\Vert {\mathbf{\Phi }}^H{\mathbf{r}}_{p,{\tau }_h}\Vert }_{\infty }}{\Vert {\mathbf{r}}_{p,\tau }{\Vert }_2}.$$

$$({\text{MMV}}_{\text{TM}}{\text{BP}}_{\tilde{\sigma }})\text{min}\kappa (\mathbf{J})\text{s.t.}\Vert \mathbf{\Phi }\mathbf{J}-\mathbf{Y}{\Vert }_F\le \tilde{\sigma },$$

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

$$({\text{MMV}}_{\text{TM}}{\text{LS}}_{\tau })\text{min}{\Vert \mathbf{\Phi }\mathbf{J}-\mathbf{Y}\Vert }_F\text{s.t.}{\Vert \mathbf{J}\Vert }_{1,2}\le \tau ,$$

$${\varphi }_{{\text{MMV}}_{\text{TM}}}(\tau )={\Vert \mathbf{\Phi }{\mathbf{J}}_{\tau }-\mathbf{Y}\Vert }_F,$$

$${\varphi }_{{\text{MMV}}_{\text{TM}}}^{\prime }({\tau }_h)=-\frac{{\Vert {\mathbf{\Phi }}^H(\mathbf{\Phi }{\mathbf{J}}_{{\tau }_h}-\mathbf{Y})\Vert }_{\infty ,2}}{{\Vert \mathbf{\Phi }{\mathbf{J}}_{{\tau }_h}-\mathbf{Y}\Vert }_F},$$

$${\tau }_{h+1}={\tau }_h+\frac{\tilde{\sigma }-{\varphi }_{{\text{MMV}}_{\text{TM}}}\left({\tau }_h\right)}{{\varphi }_{{\text{MMV}}_{\text{TM}}}^{\prime }\left({\tau }_h\right)}.$$

$${\mathcal{P}}_{\tau ,{\text{MMV}}_{\text{TM}}}[\mathbf{J}]:=\left\{\underset{\mathbf{S}}{\mathrm{argmin}}{\Vert \mathbf{J}-\mathbf{S}\Vert }_F\text{s.t.}\Vert \mathbf{S}{\Vert }_{1,2}\le \tau \right\}.$$

$$({\text{MMV}}_{\text{TE}}{\text{BP}}_{\tilde{\sigma }})\text{min}{\kappa }_{\text{TE}}(\mathbf{J})\text{s.t.}\rho \left(\mathbf{\Phi }\mathbf{J}-\mathbf{Y}\right)\le \tilde{\sigma },$$

$${\kappa }_{\text{TE}}(\mathbf{J}):=\sum _{n=1}^N{\Vert {\left[{\mathbf{J}}_{2n-1,:}{\mathbf{J}}_{2n,:}\right]}^T\Vert }_2,$$

$$\rho (\cdot ):=\Vert \cdot {\Vert }_F,$$

$$({\text{MMV}}_{\text{TE}}{\text{LS}}_{\tau })\text{min}\rho \left(\mathbf{\Phi }\mathbf{J}-\mathbf{Y}\right)\text{s.t.}{\kappa }_{\text{TE}}(\mathbf{J})\le \tau .$$

Full text

1 Introduction

Inverse scattering is a procedure of recovering the characteristics of the objects from the scattered fields. It is of great importance because of its wide applications in different areas, such as seismic detection, medical imaging, sonar, remote sensing, and so forth. Most of the studies on the inverse scattering problems are focused on the frequencies of the resonant region, i.e., the wavelength is comparable to the dimension of the object. Challenges mainly lie in the nonlinearity and ill-posedness in the Hadamard sense [1]. There is a large variety of possible inverse scattering problems, for example, find the shape of the scatterer with the boundary condition already known, or find the space dependent coefficients of the object without any a priori information at all. The inverse scattering problem discussed in this paper is to determine the shape of the highly conductive scatterers with the scattered electromagnetic (EM) field for one or several incident fields at one single frequency of the resonant region.

Basically, there are two families of methods for solving this problem: the volume-based methods and the surface-based methods. The volume-based methods indicate the shape with space dependent coefficients. P. M. van den Berg proposed to retrieve the boundary of the highly conductive scatterer by doing the same with the iterative method of reconstructing the conductivity of an EM penetrable objects [2]. This idea was further extended to the mixed dielectric and highly conductive objects combined with contrast source inversion (CSI) method [3], see O Féron et al. [4], Chun Yu et al. [5]. Classical design of the cost functional consists of the data error and the state error without considering the cross-correlated mismatch of both errors. Recently, a so-called cross-correlated error function and a novel inversion method, referred to as cross-correlated CSI, have been proposed to enhance the inversion ability [6]. The idea of solving the nonlinear inverse scattering problem with linear algebra can also be found in [7]. An algorithmically efficient algorithm, the time-reversal multiple signal classification (TR-MUSIC) method [8, 9], is also of interest since it is capable to break through the diffraction limit. Linear sampling method (LSM) [10, 11, 12] is another typical volume-based method, which finds an indicator function for each voxel in the region of interest (ROI) by first defining a far-field mapping operator (or a near-field mapping operator [13]), then solving a linear system of equations. The norm of the indicator function approaches to zero when the position of the corresponding voxel approaches the highly conductive scatterer. Although LSM has been proved to be effective for highly conductive scatterer, and also applicable for dielectric scatterer in some cases [14], this method needs sufficient amount of measurements to guarantee the inversion performance [15]. Besides, it is very time-consuming to compute the dyadic Green function for an irregular inhomogeneous background [16], for instance, in the case of ground penetrating radar (GPR) [17]. The surface-based methods first parameterize the shape of the scatterer mathematically with several parameters, then optimize the parameters by minimizing a cost function iteratively [18]. The drawback of this method is obvious: 1) This type of methods require a priori information on the position and quantity of the scatterers, more research on this point can be found in [19, 20]; 2) It is intractable to deal with the complicated nonconvex objects. Apart from that, each iteration involves finding a solution to a forward scattering problem, which is extremely time-consuming for the large-scale inverse problems with an irregular background. As a matter of fact, this is a general drawback of the iterative inversion methods. In cases where the dimension of the solution space is not so huge, global optimization techniques [21, 22, 23] are good candidates to search for the global optimal solution. We also refer to [24] for a compressive sensing CSI method which solves the contrast source two-step formulation for detecting the non-radiating part of the equivalent currents.

Recently we have proposed a linear model to address the nonlinear highly conductive inverse scattering problem with transverse magnetic (TM) polarized incident fields [25]. The basic idea is to transfer the problem to a set of linear inverse source problems, formulate the set of problems with the multiple measurement vectors (MMV) model [26], and finally solve the problem with the sum-of-norm regularization constraint. We have also considered a cascade of the linear inverse source model and a linear optimization model for solving the 3-D inverse scattering problem with half-space configurations [27]. Although the feasibility of this idea has been demonstrated numerically, a theoretical framework has not been established yet. The extension in solving the inverse scattering problem of vectorial fields is not straightforward and therefore has to be proved theoretically. Moreover, the feasibility needs to be further demonstrated with experimental measurement data. In this paper, we first presented the theoretical framework in solving the highly conductive inverse scattering problem with TM-polarized incident fields. Based on the convex optimization theory and the spectral projected gradient method (SPGL1) [28, 29], we extended the theoretical framework to solve the transverse electric (TE) polarized (vectorial field) inverse scattering problems. Consequently, the extension in 3-D problems can be derived straightforwardly. Cross-validation (CV) technique [30, 31] is used to terminate the iteration such that the estimation of noise level is well circumvented. Both 2-D synthetic data generated by a MATLAB-based “MaxwellFDFD” package [32] and the experimental TM- and TE-polarized datasets of the Institut Fresnel, France [33] are processed for demonstrating the validity of the proposed model. What’s more, we have also presented an analysis on the computational complexity and the affect of transmitter/receiver number on the imaging quality of the proposed approach, which sheds a light on the design of the imaging system.

For the case of penetrable objects, the contrast sources are distributed everywhere in the interior of the object. Since the linear model is regularized with sum-of-norm constraint, the reconstruction algorithm tends to seek for a solution of the minimum sum-of-norm. According to the field equivalence principle [34], for penetrable object, there would more likely be several sparse solutions which not only generate the same scattered field pattern (satisfying the data equation) but also have the same non-zero structure (possessing smaller sum-of-norm than the real solutions). Since correct recovery cannot be guaranteed for penetrable objects, we restricted the discussion in this paper to highly conductive objects.

In a nutshell, major differences of the proposed approach in comparison to other existing methods are as follows:

In comparison to linear methods with linear model of weak scattering assumptions (such as Born/Rytov approximations), the proposed model is more applicable since no weak scattering approximation has been made. 2. 2.

In comparison to linear methods with linear model of no weak scattering approximations (such as linear sampling method), the joint sparse structure of the contrast sources is enforced in the proposed approach by use of sum-of-norm regularization constraint, resulting in higher resolving ability. 3. 3.

In comparison to linear iterative algebra of multiple levels with non-linear model (such as CSI), calculation of the total fields for the proposed approach is not needed, resulting in higher imaging efficiency. 4. 4.

In comparison to super resolution methods based on pseudo-spectrum analysis (such as Time-reversal MUSIC), the proposed approach does not need to estimate the scatterer number, nor does it need to care about how the imaging domain is discretized.

Since the proposed approach is also based on a linear model of no weak scattering approximation, LSM and an improved version are selected in this paper for comparison.

The remainder of the paper is organized as follows: In Subsection 2.1, the problem statement is given; In Subsection 2.2, the formulation of the linear model is presented; In Subsection 2.3, the SPGL1 method for solving the single measurement vector (SMV) model in TM case is introduced; In Subsection 2.4 and 2.5, we derived the sum-of-norm optimization method for solving the MMV model111The multi-frequency implementation of this method (referred to as GMMV-LIM) is available at https://github.com/TUDsun/GMMV-LIM. of TM case and TE case, respectively; In Subsection 2.6, a CV-based modified SPGL1 method is introduced; In Section 3, linear sampling method and its improved version are introduced; The inverted results with synthetic data and experimental data are given in Section 4 and Section 5, respectively; Finally, Section 6 ends the paper with our conclusions.

2 MMV Linear Inversion Model

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}$. The domain $\mathcal{D}$ contains a highly conductive cylinder $\Gamma$, whose surface is represented by $\partial\Gamma$. The dielectric properties of the background are known beforehand. The domain $\mathcal{S}$ contains the sources and receivers. The sources are denoted by the subscript $p$ (where $p\in\{1,2,3...,P\}$), and the receivers are denoted by the subscript $q$ (where $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.

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. In this paper, we consider the time factor $\exp(\text{i}\omega t)$. Here, i represents the imaginary unit, then the vectorial Maxwell’s equations in the frequency domain can be written as

[TABLE]

where $\bm{E}$ and $\bm{H}$ are the electric and magnetic fields, respectively; $\bm{J}$ and $\bm{M}$ are the electric and magnetic current source densities, respectively; $\bm{\sigma}$, $\bm{\varepsilon}$ and $\bm{\mu}$ are the electric conductivity, electric permittivity and magnetic permeability, respectively. For most of the real problems, $\bm{\mu}$ can be reasonably assumed to be the permeability of free space, $\mu_{0}$, while $\bm{\sigma}$ and $\bm{\varepsilon}$ are functions of both the position vector $\bm{x}=[x_{1},x_{2},x_{3}]^{T}$ and the angular frequency $\omega$.

Since the solution of the $E$-field is two orders of magnitude larger than the $H$-field, it is for better numerical accuracy (see [32]) to eliminate either the $E$-field or the $H$-field from equation (1). In this paper, we assume the electric field is measured, so we rewrite our equations in terms of the electric field $\bm{E}$ according to

[TABLE]

where $\bm{\epsilon}$ is the complex permittivity given by

[TABLE]

The total electric field, $\bm{E}_{p}$, and the incident electric field, $\bm{E}_{p}^{\text{inc}}$, are excited by the $p$-th external source, the scattered electric field $\bm{E}_{p}^{{\text{sct}}}$ is then found by

[TABLE]

Here, the subscript $p$ corresponds to the $p$-th source. According to the above relation and the electric field equation (2), it is easy to obtain the basic equation of the inverse scattering problem, which is

[TABLE]

where, the contrast $\bm{\chi}$ is the difference of the complex permittivity of the test domain, $\bm{\epsilon}$, and that of the background, $\bm{\epsilon}_{\text{b}}$, i.e., $\bm{\chi}=\bm{\epsilon}-\bm{\epsilon}_{\text{b}}$. The problem we are going to resolve is to find the shape of the highly conductive scatterer $\partial\Gamma$ from the measurement of the scattered electric fields $\bm{E}_{p}^{{\text{sct}}}$. Since the total electric field $\bm{E}_{p}$ is a function of the contrast $\bm{\chi}$, this is obviously a nonlinear problem.

2.2 Formulation

First, let us formulate the inverse problem following the vector form of the FDFD scheme in [32], and rewrite (5) as follows

[TABLE]

where, $\bm{A}$ is the stiffness matrix, and

[TABLE]

and $\bm{\chi}\bm{e}_{p}$ is the component-wise multiplication of the two vectors, $\bm{\chi}$ and $\bm{e}_{p}$. The scattered fields $\bm{e}^{\text{sct}}_{p}=\bm{A}^{-1}\omega^{2}\bm{\chi}\bm{e}_{p}$ are probed and the measurements are used to estimate the unknown $\bm{\chi}$. Now let us use a measurement operator, $\mathcal{M}_{\mathcal{S}}$, to select the field values at the receiver positions, then we can formulate the data equations as $\bm{y}_{p}=\mathcal{M}_{\mathcal{S}}\bm{A}^{-1}\omega^{2}\bm{\chi}\bm{e}_{p}$, where $\bm{y}_{p}$ is the measurement vector of the $p$-th scattered fields. Let us further define the contrast source as $\bm{j}_{p}:=\bm{\chi}\bm{e}_{p}$, then we have

[TABLE]

where, $\bm{\Phi}$ is the sensing matrix defined by

[TABLE]

In TM case, $\bm{\Phi}\in\mathbb{C}^{Q\times N}$, while in TE case, $\bm{\Phi}\in\mathbb{C}^{2Q\times 2N}$. Here, $N$ is the grid number of the discretized inversion domain.

Since the contrast source $\bm{j}_{p}$ shows sparsity due to the fact that the induced current only exists on the surface of the highly conductive objects, the ill-posedness of the inverse scattering problem can be overcome by exploiting the sparsity of the contrast sources. Further, although the contrast sources $\bm{j}_{p}$ excited by the illumination of the incident fields $\bm{e}^{\text{inc}}_{p}$ are of different values, the non-zero elements are located at the same positions — the boundary of the highly conductive scatterers. This inspired us to improve the inversion performance by enhancing the joint sparse structure, so the linear data model (8) is further formulated as an multiple measurement vectors model. In the following part of this paper, a linear model is constructed and a sum-of-norm optimization problem is derived for TM and TE cases, respectively. In doing so, the nonlinear inverse scattering problem can be simplified and addressed by a linear optimization scheme without considering the state equations (i.e., the calculation of the total fields is circumvented)

[TABLE]

Specifically, Eq. (8) is rewritten as follows

[TABLE]

where, $\bm{\Phi}$ is the joint sensing matrix,

[TABLE]

is the contrast source matrix, and

[TABLE]

is the measurement data matrix, and $\bm{U}$ represents the additive complex measurement noise matrix.

2.3 Solving the SMV model: TM case

First, consider the single source configuration illuminated by TM-polarized wave. The inverse scattering problem is formulated as a basis pursuit denoise (BP${}_{\widetilde{\sigma}}$) problem [35]:

[TABLE]

where, $\widetilde{\sigma}$ represents the noise level, and the contrast source is regularized with the $\ell_{1}$-norm constraint. Solving this problem means searching for a solution $\bm{j}_{p}$ which is of the smallest $\ell_{1}$-norm and meanwhile satisfies the inequality condition. Although the BP${}_{\widetilde{\sigma}}$ problem is straightforward for understanding the inverse problem, it is not easy to solve directly even if we exactly know the value of $\widetilde{\sigma}$. An equivalent problem that is much simpler to solve is the Lasso (LSτ) problem [36], which is formulated as

[TABLE]

The LSτ problem can be solved using a spectral projected gradient (SPG) method that is proposed based on convex optimization theory [37, 38, 39]. Details of the SPG method for solving the (LSτ) problem (13) is given by [28, Algorithm 1], in which $\mathcal{P}_{\tau}[\cdot]$ is a projection operator defined as

[TABLE]

$\mathcal{P}_{\tau}[\cdot]$ gives the projection of a vector $\bm{j}_{p}$ onto the one-norm ball with radius $\tau$.

In practice, $\tau$ is usually not available. For solving this problem, a Pareto curve is defined in SPGL1 algorithm [28] by

[TABLE]

where, $\bm{j}_{p,\tau}$ is the optimal solution to the LSτ problem. It is easy to find that the (BP${}_{\widetilde{\sigma}}$) problem is equivalent to the (LSτ) problem when $\phi_{\text{SMV}_{\text{TM}}}(\tau)=\widetilde{\sigma}$ is satisfied.

The Pareto curve is proved to be differentiable under some conditions, and the root of the nonlinear equation $\phi_{\text{SMV}_{\text{TM}}}(\tau)=\widetilde{\sigma}$ can be reached by Newton iterations [28]

[TABLE]

where,

[TABLE]

where, $\bm{\Phi}^{H}$ is the conjugate transpose of matrix $\bm{\Phi}$, $\bm{r}_{p,\tau_{h}}=\bm{\Phi}\bm{j}_{p,\tau_{h}}-\bm{y}_{p}$ is the residual vector. The update of $\tau$ by probing the Pareto curve is illustrated in Figure 2. This procedure requires computing successively more accurate solutions of LSτ. The Newton root-finding framework for solving the ($\text{SMV}_{\text{TM}}\ \text{BP}_{\widetilde{\sigma}}$) problem is given in [29, Algorithm 1].

2.4 Solving the MMV model: TM case

Let us now first consider the 2-D multi-source configuration with TM-polarized illumination. As all the contrast sources are focused on the boundary of the scatterers, the contrast source matrix $\bm{J}$ shows row-sparsity. Therefore, the inverse scattering problem with multi-source configurations can be formulated as an ($\text{MMV}_{\text{TM}}\ \text{BP}_{\widetilde{\sigma}}$) problem regularized by a sum-of-norm constraint

[TABLE]

where, $\kappa(\bm{J})$ is the mixed ($\alpha,\beta$)-norm defined as

[TABLE]

with $\bm{J}_{i,:}$ denoting the $i$-th row of $\bm{J}$, and $\|\cdot\|_{\beta}$ the conventional $\beta$-norm. $\|\cdot\|_{F}$ is the Frobenius norm which is equivalent to the mixed (2,2)-norm $\|\cdot\|_{2,2}$. In this problem, we consider $\alpha=1$ and $\beta=2$, which is a sum-of-norm constraint. Accordingly, the $\text{MMV}_{\text{TM}}\ \text{LS}_{\tau}$ problem is reformulated as

[TABLE]

and the Pareto curve for the MMV model is defined as

[TABLE]

where, $\bm{J}_{\tau}$ is the optimal solution to the LSτ problem (20).

According to [29, Theorem 2.2] and [40, Chapter 5], $\phi_{\text{MMV}_{\text{TM}}}(\tau)$ is continuously differentiable and

[TABLE]

where, $\|\cdot\|_{\infty,2}$ is the dual norm of $\|\cdot\|_{1,2}$ [29, Corollary 6.2]. Similarly, the root of the nonlinear equation $\phi_{\text{MMV}_{\text{TM}}}(\tau)=\widetilde{\sigma}$ can also be reached by Newton iterations

[TABLE]

The projection operator $\mathcal{P}_{\tau}[\cdot]$ is replaced by an orthogonal projection onto $\|\cdot\|_{1,2}$ balls, $\mathcal{P}_{\tau,\text{MMV}_{\text{TM}}}[\cdot]$, which is defined as follows

[TABLE]

We refer to [29, Theorem 6.3] for the implementation of the projection operator. The (MMV${}_{\text{TM}}$ BP${}_{\widetilde{\sigma}}$) problem is solved by Algorithm 1 and Algorithm 2 with the Pareto curve, $\phi_{\text{MMV}_{\text{TM}}}(\tau)$, its derivative with respective to $\tau$, $\phi^{\prime}_{\text{MMV}_{\text{TM}}}(\tau)$, and the projection operator, $\mathcal{P}_{\tau,\text{MMV}_{\text{TM}}}[\cdot]$, defined by (21), (22), and (24), respectively.

2.5 Solving the MMV model: TE case

For the TE polarization case, the electric field is not a scalar anymore. Therefore, care must be given to the formulation of the (MMV${}_{\text{TE}}$ BP${}_{\widetilde{\sigma}}$) problem. Considering the two components of electric field, $E_{x}$ and $E_{y}$, the inverse scattering problem for the TE case can be formulated as

[TABLE]

where,

[TABLE]

and

[TABLE]

are gauge functions [41]. The $\text{MMV}_{\text{TE}}\ \text{LS}_{\tau}$ problem is formulated accordingly as

[TABLE]

2.5.1 Derivation of the dual

Let us rewrite (28) in terms of $\bm{J}$ and an explicit residual term $\bm{R}$

[TABLE]

The dual to this equivalent problem is given by [42, Chapter 5]

[TABLE]

where $\bm{Z}\in\mathbb{C}^{(2M)\times P}$ and $\lambda\in\mathbb{C}$ are dual variables, and $\mathcal{G}$ is the Lagrange dual function, given by

[TABLE]

where Tr represents the trace of a matrix. By separability of the infimum over $\bm{J}$ and $\bm{R}$ we can rewrite $\mathcal{G}$ in terms of two separate suprema,

[TABLE]

It is easy to see that the first supremum is the conjugate function of $\rho$ and the second supremum is the conjugate function of $\kappa_{\text{TE}}$ [42, Chapter 3.3], by noting that

[TABLE]

and

[TABLE]

respectively. Here, $\text{vec}\{\cdot\}$ is the vectorization of a matrix, $\widetilde{\bm{Z}}=\bm{\Phi}^{H}\bm{Z}\in\mathbb{C}^{(2N)\times P}$, and $\kappa_{\text{TE}}\left(\text{vec}\{\bm{J}\}\right)$ is defined equivalently as $\kappa_{\text{TE}}(\bm{J})$ in (26). Therefore, we have

[TABLE]

and

[TABLE]

where, the polar of $\rho$ and $\kappa_{\text{TE}}$ are defined by

[TABLE]

and

[TABLE]

respectively. If the gauge function is a norm, the polar reduces to the dual norm [42, Section 3.3.1], i.e., $\rho^{o}(\bm{Z})=\|\bm{Z}\|_{F}$ and

[TABLE]

(for more details see [29, Corollary 6.2]). Substitution of (35) and (36) into (32) yields

[TABLE]

In the case $\rho(\cdot)=\|\cdot\|_{F}$, the dual variable $\bm{Z}$ can be easily derived from

[TABLE]

which is $\bm{Z}=\frac{\bm{R}}{\|\bm{R}\|_{F}}$. To derive the optimal $\lambda$, we can observe from (40) that as long as $\tau>0$, $\lambda$ must be at its lower bound $\kappa_{\text{TE}}^{o}(\widetilde{\bm{Z}})$, otherwise one can increase the objective $\text{Tr}\left\{\bm{Y}^{H}\bm{Z}\right\}-\tau\lambda$. Therefore, we obtain

[TABLE]

According to [40, Theorem 5.2], we know that, on the open interval $\tau\in(0,\tau_{0})$, where

[TABLE]

the Pareto curve $\phi_{\text{MMV}_{\text{TE}}}(\tau)=\rho(\bm{R})$ is strictly decreasing, and continuously differentiable with

[TABLE]

The projection operator $\mathcal{P}_{\tau,\text{MMV}_{\text{TM}}}[\cdot]$ is replaced by an orthogonal projection onto $\kappa_{\text{TE}}(\cdot)$ balls, $\mathcal{P}_{\tau,\text{MMV}_{\text{TE}}}[\cdot]$, which is defined as follows

[TABLE]

With a simple matrix transformation of $\overline{\bm{J}}_{n,:}=\left[\bm{J}_{2n-1,:}\ \bm{J}_{2n,:}\right]$ and $\overline{\bm{X}}_{n,:}=\left[\bm{X}_{2n-1,:}\ \bm{X}_{2n,:}\right]$, we can rewrite (44) as follows

[TABLE]

In doing so, the projection operator in TE case satisfies [29, Theorem 6.3]. The (MMV${}_{\text{TE}}$ BP${}_{\widetilde{\sigma}}$) problem is solved by Algorithm 1 and Algorithm 2 with the Pareto curve, its derivative with respect to $\tau$, and the projection operator replaced by $\phi_{\text{MMV}_{\text{TE}}}(\tau)$, $\phi^{\prime}_{\text{MMV}_{\text{TE}}}(\tau)$, and $\mathcal{P}_{\tau,\text{MMV}_{\text{TE}}}[\cdot]$, respectively.

2.6 CV-based Modified SPGL1

In real applications, the noise level, i.e., the parameter $\widetilde{\sigma}$, is generally unknown, which means the termination condition, $\phi_{\text{MMV}}(\tau)=\widetilde{\sigma}$, does not work anymore. In order to deal with this problem, we modified the SPGL1 method based on the CV technique [30, 31], in which $\widetilde{\sigma}$ is set 0 and the iteration is terminated using CV technique. In doing so, the problem of estimating the noise level, i.e., the parameter $\widetilde{\sigma}$, can be well circumvented.

Cross validation is a statistical technique that separates a data-set into a training (estimation) set and a testing (cross validation) set. The training set is used to construct the model and the testing set is used to adjust the model order so that the noise is not over-fitted. The basic idea behind this technique is to sacrifice a small number of measurements in exchange of prior knowledge. Specifically, when CV is utilized in SPGL1 method, we separate the original scattering matrix to a reconstruction matrix $\bm{\Phi}_{p,r}\in\mathbb{C}^{Q_{r}\times N}$ and a CV matrix $\bm{\Phi}_{p,CV}\in\mathbb{C}^{Q_{CV}\times N}$ with $Q=Q_{r}+Q_{CV}$. The measurement vector $\bm{y}_{p}$ is also separated accordingly, to a reconstruction measurement vector $\bm{y}_{p,r}\in\mathbb{C}^{Q_{r}}$ and a CV measurement vector $\bm{y}_{p,CV}\in\mathbb{C}^{Q_{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 trend of CV residual in iteration behaves abruptly different (turns from decreasing to increasing) comparing to that of reconstruction residual, as soon as 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 smallest one.

In order to find the smallest CV residual, a maximum number, $N_{\max}$, is needed and set a large value to guarantee the smallest CV residual occurs in the range of the $N_{\max}$ iterations. In this case, a large number of iterations are performed in vain, which decreases the efficiency of the algorithm. Therefore, we consider an alternative termination condition given by

[TABLE]

where, $N_{\text{Iter}}$ is the current iteration number, $N_{\text{opt}}$ is the iteration index corresponding to the smallest CV residual — the optimal solution. The idea behind this criterion is that the CV residual is identified as the smallest one if the CV residual keeps increasing monotonously for $\Delta N$ times of iteration. In the following experimental examples, this termination condition works well with $\Delta N=30$.

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

[TABLE]

or

[TABLE]

with $n=1,2,\dots,N$, where $\bm{j}_{p,n}$, $\bm{\gamma}_{\text{MMV}_{\text{TM}}}[n]$, and $\bm{\gamma}_{\text{MMV}_{\text{TE}}}[n]$ represent the $n$-th element of vector $\bm{j}_{p}$, $\bm{\gamma}_{\text{MMV}_{\text{TM}}}$, and $\bm{\gamma}_{\text{MMV}_{\text{TE}}}$, respectively. In the end of this section, we remark that as the regularized solution corresponds to the least sum-of-norm, the non-measurable equivalent contrast sources [43] tend to be ignored.

3 Linear sampling method and its improved version

In this section, the proposed method is tested with both synthetic data and experimental data. In the meanwhile, we have also processed the same data using linear sampling method for comparison. Since the background of the experiments is free space, the LSM method consists in solving the integral equation of the indicator function

[TABLE]

and

[TABLE]

for the TM and TE cases, respectively, where, $E_{1}(\bm{x}_{r},\bm{x}_{t})$, $E_{2}(\bm{x}_{r},\bm{x}_{t})$, and $E_{3}(\bm{x}_{r},\bm{x}_{t})$ represent $x_{1}$-, $x_{2}$-, and $x_{3}$-components of the scattered field probed at $\bm{x}_{r}$ corresponding to the transmitter at $\bm{x}_{t}$, respectively; $\bm{x}_{s}$ is the sampling point in the inversion domain; $E_{i,j}(\bm{x}_{s},\bm{x}_{r})$ is $i$-th component of the electric field at $\bm{x}_{r}$ generated by an ideal electric dipole located at $\bm{x}_{s}$ with the polarization vector parallel to $x_{j}$-axis, which are given by

[TABLE]

where, $R=\|\bm{x}_{s}-\bm{x}_{r}\|_{2}$. Eq. (51) and (52) can be reformulated as a set of linear systems of equations

[TABLE]

where, $\bm{Y}$ is the measurement data matrix, $\bm{g}_{\bm{x}_{s}}$ is the indicator function of the sampling point $\bm{x}_{s}$ in the form of a column vector, $\bm{f}_{\bm{x}_{s}}$ is the right side of Eq. (51) in the form of a column vector. Following the same approach of solving Eq. (54) in [44, 45], the shape of the scatterers is defined by

[TABLE]

where, $\|\bm{g}_{\bm{x}_{s}}\|^{2}$ is given by

[TABLE]

where, $s_{d}$ represents the singular value of matrix $\bm{Y}$ corresponding to the singular vector $\bm{u}_{d}$, $D=\min\{P,Q\}$, and $a=0.01\times\underset{d}{\max}\{s_{d}\}$.

We have also considered, in the TM cases, the improved linear sampling method proposed in [12] in the comparison of the proposed method and LSM. The indicator function of improved LSM is defined as

[TABLE]

where, $a$ is the radius of a smallest ball that covers the targets, the power $\frac{1}{2I}$ is normalization factor, and $\bm{g}^{x}_{i,\bm{x}_{s}}$ and $\bm{g}^{y}_{i,\bm{x}_{s}}$ are obtained by replacing $E_{3,3}(\bm{x}_{s},\bm{x}_{r})$ in Eq. (51) with $\varphi^{x}_{i}(\bm{x}_{s},\bm{x}_{r})$ and $\varphi^{y}_{i}(\bm{x}_{s},\bm{x}_{r})$, respectively,

[TABLE]

where, $\phi_{r}$ and $\phi_{s}$ are the angular components of $\bm{x}_{s}$ and $\bm{x}_{r}$ in polar coordinate system, respectively. We refer to [12] for more details of this indicator function.

It is worth mentioning that both the contrast source $\bm{j}_{p}$ and the indicator function $\bm{g}_{\bm{x}}$ are proportional to the amplitude of the electric field. According to the definition in Eq. (49), Eq. (50), and Eq. (55), $\bm{\gamma}_{\text{MMV}}$ and $\bm{\gamma}_{\text{LSM}}$ are proportional and inversely proportional to the power of the electric field, respectively. Therefore, the dB scaling shown in the following examples is defined as

[TABLE]

4 Synthetic Data Imaging

In this subsection, the proposed method is tested with synthetic data. The transmitting antenna is simulated for simplicity with an ideal electric dipole (TM-polarization case) and an ideal magnetic dipole (TE-polarization case). Coordinate system is established such that the dielectric parameters are variable along the $x$- and $y$-axes, but invariable along $z$-axis. The transmitting antenna rotates on a circular orbit of 3 m radius centering at the origin (0,0). The receiving positions are taken on the same orbit without any position close than 30∘ from the transmitting antenna. The measurement configuration of Simulation 1 and 2 is shown in Figure 3, in which the selection of CV measurements and reconstruction measurements is illustrated. Empirically, an arc length $\geq\lambda/3$ is a good selection. The number of the CV receivers on each arc depends on how dense the receiver positions are, and the total CV receiver number is around 20% of the total measurement number [31]. The operating frequency is 500 MHz. Two configurations of different objects are considered. One is combined with two circular metallic cylinders and the other one is a metallic cylinder with a “crescent-shaped” cross section. The radius of the circular cross section is $0.2$ m ($=\lambda/3$), and the centers of the two circles are ($-0.45$, $0.6$) and ($0.45$, $0.6$), respectively. The crescent is the subtraction of two circles of radius $0.6$ m ($=\lambda$) centering at ([math], [math]) and ($0.4$, [math]). See Figure 6(a) and Figure 10(a) for their true geometry. The forward EM scattering problem is solved by a MATLAB-based 3-D FDFD package “MaxwellFDFD” [32]. The technique of non-uniform staggered grids is used to reduce the computational burden, while for inverting the measurement data, we consider uniform discretization such that an inverse crime is circumvented. In the forward solver, we consider a fine grid size of $\lambda/(45\sqrt{\epsilon_{r}})$. The data for inversion is obtained by subtracting the incident field from the total field. Periodic boundary conditions are imposed on the design of the FDFD stiffness matrix in order to simulate the 2-D configuration. Perfect matching layer (PML) is used to simulate an anechoic chamber environment.

4.1 Determine the measurement configuration

To determine the measurement configuration, we need to investigate the relationship between the transceiver numbers and the imaging quality. Let us first consider 120 transmitters and 151 receivers, i.e., the transmitter rotates on the circular orbit with a step of 3∘, and the receiver rotates on the measurement arc of 300∘ with a step of 2∘. The CV receivers are selected in the same way shown in Figure 3, but there are 4 continuous CV receiver positions in each CV arc (equivalent to 8 ∘). Now let us disturb the measurement data (the scattered fields) with Gaussian additive random noise of 30 dB signal-to-noise ratio (SNR), and then process the data by the proposed method. If we use the reconstructed image as the reference image, denoted by $\bm{\gamma}_{\text{ref}}$, then a correlation coefficient can be defined as

[TABLE]

where, $\bm{\gamma}$ denotes the MMV image with different measurement configurations and noise levels, $\overline{\gamma_{\text{ref}}}$ and $\overline{\gamma}$ are the mean values of $\bm{\gamma}_{\text{ref}}$ and $\bm{\gamma}$, respectively. The correlation coefficient reflects the similarity degree of two images. The minor negative correlation coefficients are forced to zeros, as negative correlation does not make any sense for two amplitude images.

Now we first fix the receiver number to 151, and calculate the correlation coefficients of Simulation 1 and Simulation 2 with different transmitter numbers. Fig. 4 (a) and (b) show the correlation coefficient curves in terms of transmitter number by processing the TM-polarized data and TE-polarized data, respectively. Two SNRs, 10 dB and 30 dB, are considered. From Fig. 4 we observe that obvious decrease of correlation coefficient occurs at 18 transmitters, indicating that the image quality gets worse when transmitter number is less than 18. The correlation coefficient curves of 10 dB and 30 dB maintain the same trend, and the correlation coefficients of 10 dB maintains above 0.95 when more than 18 transmitters are used, indicating the proposed method is robust against the Gaussian additive random noise. Then we fix the transmitter number to 18 and 120, respectively, and image the targets in Simulation 1 and 2 with different receiver numbers. Since CV technique needs enough amount of measurements, the noise level is assumed exactly known when the receiver number is less than or equal to 31. Fig. 5(a) and (b) show the correlation coefficient curves in terms of receiver number by processing the TM-polarized data and TE-polarized data, respectively. Two SNRs, 10 dB and 30 dB, are considered. From Fig. 5 we observe that the smallest receiver number for ensure a stable imaging quality is 16. The correlation coefficient of 18 transmitters and 10 dB SNR maintains $r_{\text{corr}}\geq 0.90$ when the receiver number $\geq 16$.

Since the reference image in the definition of the correlation coefficient is not the real shape of the targets, it is actually an asymptotic measure of the imaging quality. Discussion of the imaging results are given in the next subsection in comparison to the LSM images to further investigate the imaging performance. In the following subsection of imaging results and the next section of experimental data imaging, we select 18 transmitters (equivalent to an interval of 20∘) for two reasons: 1) the proposed method works well in the numerical simulations with 18 transmitters; 2) more experiments and targets are required to demonstrate the good imaging performance when 18 transmitters are used. Note that the noise level is not available in real applications, 61 receivers (equivalent to an interval of 5∘) are selected in the numerical experiments for the use of CV technique.

4.2 Imaging results

4.2.1 Simulation 1

To reduce the computational cost, we restrict the inversion domain to $[-1.0,1.0]\times[-0.4,1.6]$ m2. The inversion domain is discretized with equal grid size of $0.01$ m ($=\lambda/60$). Let us first consider the TM-polarized data disturbed with Gaussian additive random noise of 30 dB SNR. The residual curves are shown in Figure 8(a). The trend of the residual curves is like staircases, and each step corresponds to one update of the parameter $\tau$. The CV residual starts to increase after 80 iterations, and $\Delta N=30$ more iterations are performed before termination. The solution of the minimum CV residual is selected as the optimal solution. The scatterer shape reconstructed by MMV, LSM, and improved LSM with $I=7$ is shown in Figure 6(b), (c) and (d), respectively. By comparison of Figure 6(c) and (d), it is observed that the artifacts between the two circular cylinders are suppressed by improved LSM. However, the average amplitude of the sidelobes in the region of no targets increases from $-15$ dB to $-10$ dB. From Figure 6(b) we observe that the proposed method shows higher resolution and lower sidelobes in comparison to 6(c) and (d), indicating that the resolving ability of the proposed method is better than LSM. To study the imaging performance with different SNRs, let us decrease the SNR to 10 dB, and the images results are shown in Figure 7. By comparing Figure 6 and Figure 7, we can observe obvious degradation in the LSM images, while there is no obvious degradation in the MMV images. Figure 8(b) shows the residual curves, from which we can see the reconstruction residual of the optimal solution, $r_{\text{rec}}\approx 0.105$, is higher than that of Figure 8(a), $r_{\text{rec}}\approx 0.025$.

Now let us process the TE-polarized data of different SNRs, 30 dB and 10 dB. Figure 9 shows the scatterer shape reconstructed by MMV and LSM, respectively. The imaging results demonstrate again that, in the perspective of resolving ability, the proposed method outperforms the linear sampling method. In addition, the proposed method maintains good imaging performance for different SNRs.

4.2.2 Simulation 2

In the second simulation, we restrict the inversion domain to $[-1.0,1.0]\times[-1.0,1.0]$ m2, in which the target is fully covered. The inversion domain is discretized with equal grid size of $0.01$ m ($=\lambda/60$). First, we process the TM-polarized data of 30 SNR by MMV and LSM, respectively. Figure 10(b), (c), and (d) show the reconstructed shape by MMV, LSM, and improved LSM, respectively. We can see from the results that the boundary at the left side is well reconstructed by the three methods, while the arch at the right side shows more artifacts in both Figure 10(b) and (c), because the arch at the right side is concave and multi-path scattering is more severe than the left side which is convex. Comparison of Figure 10(c) and (d) shows minor suppression of artifacts in the interior of the cylinder by improved LSM. The imaging results of 10 dB SNR data are shown in Figure 11. Apart from some minor artifacts, no obvious degradation occurs in MMV image, while we can observe severe degradation of image resolution in LSM images.

The MMV image and LSM image obtained by processing the TE-polarized data of 30 dB SNR and 10 dB SNR are shown in Figure 12. From the results, we can observe that the proposed method is able to reconstruction the scatterer’s shape with some artifacts occur at the concave side, while LSM fail to give the basic profile of the target. Considering the length of this paper, the residual curves in this simulation are not given.

4.3 Analysis of computational complexity

The sensing matrix $\bm{\Phi}$ can be computed (or analytically given for the experiments in homogeneous background) and stored beforehand. It is easy to see from Algorithm 1 and Algorithm 2 that, the computational complexity of the GMMV-based linear method primarily depends on the number of matrix-vector multiplications, $\bm{\Phi}\bm{J}$ and $\bm{\Phi}^{H}\bm{R}$. Empirically, each iteration involves maximally 2 times of $\bm{\Phi}\bm{J}$ and 1 time of $\bm{\Phi}^{H}\bm{R}$. In order to study the computational complexity of the proposed algorithms, we use one complex data multiplication as the measurement unit. The computational complexity for computing $\bm{\Phi}\bm{J}$ and $\bm{\Phi}^{H}\bm{R}$ is $QN$ for the TM case and $4QN$ for the TE case. Let us use $L$ to denote the iteration number, then the computational complexity of the proposed method in the TM and TE case is

[TABLE]

respectively. When the mesh gets finer or the inversion domain gets larger, the iteration number, $L$, almost keeps unchanged, and the running time therefore increases linearly with the grid number.

In our experiments, the imaging algorithms are implemented with MATLAB language. We ran the program on a desktop with one Intel(R) Core(TM) i5-3470 CPU @ 3.20 GHz, and we did not use parallel computing. Table 1 lists the running times of the proposed method, LSM, and improved LSM in the two simulations. As one can see that the computation time of the proposed method is hundreds of times longer than that of LSM and tens of times longer than improved LSM. The most time-consuming part of the proposed method is the matrix-vector multiplication in each iteration, while LSM only calls singular value decomposition to the measurement data matrix for once. However, the running times of the proposed method are still acceptable in view of the gain of resolving ability.

5 Experimental Data Imaging

In this subsection, we applied our method to the experimental data provided by the Remote Sensing and Microwave Experiments Team at the Institut Fresnel, France, using an HP8530 network analyzer [33]. The experimental setup consists of a large anechoic chamber, $14.50$ m long, $6.50$ m wide and $6.50$ m high, with a set of three positioners to adjust antennas or target positions. A 2-D bistatic measurement system is considered, with an emitter placed at a fixed position 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 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. In TE polarization, only the component orthogonal to both the invariance axis of the cylinder and the direction of illumination is measured. Figure 13 gives the measurement configuration, in which the selection of CV measurements and reconstruction measurements is illustrated. The time dependence in this experiment is $\exp(\text{i}\omega t)$. Therefore, after subtracting the incident field from the total field, the measurement data can be directly used for inversion. The targets we consider here are a rectangular metallic cylinder and a “U-shaped” metallic cylinder, which have been shown by Figure 14(a) and (b). Three data-sets are processed: rectTM_cent at 16 GHz, uTM_shaped at 8 GHz, and rectTE_8f at 16 GHz. The selected frequencies correspond to wavelengths which are comparable to the dimension of the targets, i.e., the resonance frequency range.

As we have discussed in Subsection 4.1, 18 transmitters of $20^{\circ}$ arc interval are enough for the proposed method to give good resolution. Therefore, in the following experiments, only 18 transmitter positions (half of the measurement data) are used for imaging. First, let us process the TM-polarized data-set: rectTM_cent at 16 GHz. The inversion domain for imaging the rectangular metallic cylinder is restricted to $[-50,50]\times[-50,50]$ mm2, and the inversion domain is discretized with equal grid size of $0.5$ mm ($=\lambda/37.5$). Figure 15a(a), (b), and (c) show the scatterer shape reconstructed by MMV, LSM, and improved LSM, respectively. The residual curves are shown in Figure 18(a). From the imaging results we observe that Figure 15a(a) shows higher resolution and less artifacts than Figure 15a(b) and (c). And we also observe that there is no big difference between Figure 15a(b) and (c).

Now let us consider the “U-shaped” metallic cylinder. The inversion domain is restricted to $[-100,100]\times[-100,100]$ mm2, and the inversion domain is discretized with equal grid size of $1$ mm ($=\lambda/37.5$). Figure 16a(a), (b), and (c) give the scatterer shape reconstructed by MMV, LSM, and improved LSM, respectively. The residual curves are shown in Figure 18(b). Severe artifacts can be observed in LSM image and improved LSM image. Furthermore, the suppression to the artifacts is not obvious in the improved LSM image. In the contrary, the “U-shaped” cylinder is reconstructed by the proposed method with the boundary well distinguished. Some artifacts can be observed vertically aligned in the interior and below the opening which are caused by the complicated scattering in the opening area.

Finally, let us process the TE-polarized data-set: rectTE_8f at 16 GHz. The scatterer shape reconstructed by MMV and LSM is shown in Figure 17(a) and (b), respectively. The residual curves are shown in Figure 18(c). It can be observed that the boundary of the rectangular metallic cylinder is not distinguishable in Figure 17(b), while in Figure 17(a) the rectangular boundary can be clearly distinguished. The data processing is performed on the same computing platform, and the running times of all the methods are listed in Table 2. In the end, we summarize this section as follows: 1) the proposed method is able to obtain higher resolution than LSM and improved LSM; 2) the proposed method is more computationally expensive than LSM and improved LSM; 3) the suppression to the artifacts by improved LSM is not obvious when less transmitters are used.

6 Conclusion

In this paper, we addressed the nonlinear inverse scattering problem of highly conductive objects with a linear model. Multiple measurement vectors model is exploited and sum-of-norm constraint is used as a regularization approach. A CV-based modified SPGL1 method is proposed to invert the data without estimating the noise level. Numerical results and experimental results of both TM polarization and TE polarization demonstrate that the proposed method shows higher resolving ability than both LSM and improved LSM. In some cases where the latter two methods fail, the proposed method can still successfully recover the profile of the targets. Numerical experiments also demonstrate that the MMV method maintains good imaging performance in the disturbance of 10 dB SNR Gaussian random noise. The running time of the proposed method is hundreds of times longer than LSM and tens of times longer than improved LSM. However, it is still promising in view of the gain of resolving performance and the linear relationship between the computational complexity and the size of the inversion domain. In addition, it might fail in the presence of not conductive scatterers, which is an obvious limitation of the proposed method.