Fast Imaging of Thin, Curve-Like Electromagnetic Inhomogeneities without a Priori Information

Abstract

It is well-known that subspace migration is a stable and effective non-iterative imaging technique in inverse scattering problem. However, for a proper application, a priori information of the shape of target must be estimated. Without this consideration, one cannot retrieve good results via subspace migration. In this paper, we identify the mathematical structure of single- and multi-frequency subspace migration without any a priori of unknown targets and explore its certain properties. This is based on the fact that elements of so-called multi-static response (MSR) matrix can be represented as an asymptotic expansion formula. Furthermore, based on the examined structure, we improve subspace migration and consider the multi-frequency subspace migration. Various results of numerical simulation with noisy data support our investigation.

1. Introduction

There exists a considerable amount of interesting inverse scattering problems concerned with the retrieval of crack-like defects completely embedded in a medium from measurement data (for related works, see [1,2,3,4,5,6] and references therein). This problem is considered an interesting and important issue because it is closely related to human life. However, due to the intrinsic ill-posedness and nonlinearity, this problem has not yet been successfully resolved.

To solve this problem, various remarkable and effective inversion techniques have been developed. The famous and general approach for solving this problem is based on the Newton-type iteration scheme, i.e., estimating the solution (shape, material properties, etc.) that minimizes the discrete norm between the true measurement data and computed data in the presence of true and man-made targets [7,8,9,10,11,12,13,14,15]. For a successful application of iterative based schemes, one need a priori information of unknown target and must begin iteration procedure with a good initial guess that is close enough to the unknown target. If not, one will need large computational costs, encounter non-convergence issue, or obtain a local minimizer instead of true solution. For this reason, generation of a good initial guess without any a priori information of target must to be considered.

To obtain an outline shape of unknown targets as a good initial guess, various non-iterative techniques have been developed and successfully applied to various inverse problems. Among them, MUltiple SIgnal Classification (MUSIC) [16,17,18], topological derivative [19,20,21], linear and direct sampling methods [22,23,24], and factorization method [25,26,27] have been applied for identifying crack-like defects or thin electromagnetic inhomogeneities.

Subspace migration is a recently investigated non-iterative algorithm in inverse scattering problem for identifying or imaging of arbitrary shaped, unknown electromagnetic inhomogeneities. To our best knowledge, it was firstly designed and applied to the edge detection of volumetric reflectors by Borcea et al. [28]. After that, it has been applied to various inverse scattering problems, for example detection of a point-like reflector when the array response matrix is obtained in a noisy environment [29] and localization of a set of perfectly conducting cracks with small length [30]. In these studies, statistical approaches are used to demonstrate that subspace migration is a very fast, effective, and robust technique. In other studies, subspace migration has been applied to identify thin electromagnetic inhomogeneities [31] in full-view inverse scattering problem, completely embedded in a half-space [32], and arbitrary shaped extended perfectly conducting cracks [33] in the limited-view inverse scattering problem. Unlike the statistical approach, the feasibility was confirmedbecause the imaging function can be expressed by the Bessel function of integer order. Based on this fact, some undiscovered properties of subspace migration have been explored. In [34], the authors considered subspace migration for obtaining an outline shape of extended electromagnetic inclusions and used the result as an initial guess to the iterative algorithm. Recently, subspace migration is applied to the microwave imaging technique for identifying small and extended anomalies from synthetic [35] and real [36] scattering parameters. It is worth mentioning that subspace migration is closely related to the time-reversal migration imaging techniques (see [37] and references therein).

Although the subspace migration algorithm has been shown to be feasible for detecting various types of inhomogeneities, this has only been confirmed when a priori information of the target is known. However, for a proper application of subspace migration, a priori information of an unknown target must be known (as shown in Equation (8) below). For this reason, most of research is conducted on the assumption that information of the number, shape (small, extended, crack-like defects, etc., mostly unit outward normal vector), or material properties (penetrable, impenetrable, permittivity of permeability contrast, etc.) of unknown targets is known. Thus, estimation or identification of a priori information of the targets seems to be considered in the early stage of the application of subspace migration. Furthermore, some artifacts also appear in the imaging result; however, there is a specific method to eliminate them.

The purpose of this paper is to apply subspace migration for imaging of a thin crack-like electromagnetic inhomogeneity located in the two-dimensional homogeneous space without any a priori information of inhomogeneity. Based on the structure of left- and right-singular vectors of multi-static response (MSR) matrix and asymptotic expansion formula for the far-field pattern in the presence of thin inhomogeneity, we perform an analysis to explain that imaging function of the subspace migration can be expressed by the combination of the Bessel functions of order 0 and 1, and unit tangential and normal vectors on the supporting curve of the thin inhomogeneity. This also leads us to discover certain properties of imaging function and to give an idea of improvement of imaging performance by filtering small magnitudes or reducing the oscillation pattern of imaging related to the Bessel functions.

The authors of [29,30,38] confirmed that application of multi-frequency is advantageous not only to improve the imaging performance but also to crucially enable self-averaging in imaging of objects. In this paper, to upgrade the imaging quality and eliminate several unexpected artifacts, we consider the multi-frequency imaging of subspace migration and confirm that application of multiple frequency successfully guarantees the imaging performance. Contrary to existing results based on statistical approaches, we confirm the improvement of imaging performance of subspace migration based on the reduction of artifacts and oscillations in the imaging results.

This paper is organized as follows. In Section 2, we briefly survey the two-dimensional direct scattering problem in the presence of thin inhomogeneity and introduce the traditional subspace migration imaging technique. In Section 3, we analyze single- and multi-frequency subspace migration imaging functions without any a priori information of thin inhomogeneities by establishing a relationship with Bessel function of integer order of the first kind, discuss certain properties of subspace migration, and introduce an improved subspace migration. In Section 4, several results of numerical experiments with noisy data are presented to support our analysis. Finally, a short conclusion is mentioned in Section 5.

2. Direct Scattering Problem and Subspace Migration

2.1. Introduction to the Two-Dimensional Direct Scattering Problem

Let $\Gamma$ be a thin, curve-like homogeneous inhomogeneity located within a two-dimensional homogeneous space $\Omega \subset {\mathbb{R}}^2$. Throughout this paper, we assume that $\Gamma$ is localized in the neighborhood of a finitely long, smooth curve $\sigma$ such that

$$\Gamma =\left\{\mathbf{x}+\eta \mathbf{n}\left(\mathbf{x}\right):\mathbf{x}\in \sigma ,-h\le \eta \le h\right\},$$

(1)

where $\mathbf{n}\left(\mathbf{x}\right)$ is the unit normal to $\sigma$ at $\mathbf{x}$ and h is a strictly positive constant that specifies the thickness of the inhomogeneity (small with respect to the wavelength) (see Figure 1). Throughout this paper, we denote $\mathbf{t}\left(\mathbf{x}\right)$ as the unit tangent vector at $\mathbf{x}\in \sigma$.

Throughout the paper, we assume that $\Gamma$ is characterized by its dielectric permittivity and magnetic permeability at a given angular frequency $\omega$. Let $0<{\epsilon }_0<+\infty$ and $0<{\mu }_0<+\infty$ denote the permittivity and permeability of the embedding space $\Omega$, and $0<{\epsilon }_{\star }<+\infty$ and $0<{\mu }_{\star }<+\infty$ those of the inhomogeneity $\Gamma$. Then, we can define the following piecewise constant dielectric permittivity and magnetic permeability

$$\epsilon \left(\mathbf{x}\right)=\left\{\begin{array}{ccc}{\epsilon }_0& \mathrm{for}& \mathbf{x}\in \Omega \backslash \overline{\Gamma }\hfill \\ {\epsilon }_{\star }& \mathrm{for}& \mathbf{x}\in \Gamma \hfill \end{array}\right.\mathrm{and}\mu \left(\mathbf{x}\right)=\left\{\begin{array}{ccc}{\mu }_0& \mathrm{for}& \mathbf{x}\in \Omega \backslash \overline{\Gamma }\hfill \\ {\mu }_{\star }& \mathrm{for}& \mathbf{x}\in \Gamma ,\hfill \end{array}\right.$$

respectively. Note that, if there is no inhomogeneity, i.e., in the homogeneous space, $\mu \left(\mathbf{x}\right)$ and $\epsilon \left(\mathbf{x}\right)$ are equal to ${\mu }_0$ and ${\epsilon }_0$, respectively. In this paper, we set ${\epsilon }_{\star }>{\epsilon }_0$ and ${\mu }_{\star }>{\mu }_0$ for convenience but exact values of ${\epsilon }_{\star }$ and ${\mu }_{\star }$ are assumed unknown.

At strictly positive angular frequency $\omega$, we let k be the background wavenumber that satisfies $k^2={\omega }^2{\epsilon }_0{\mu }_0$ and $u_{tot}(\mathbf{x};k)$ be the time-harmonic total field which satisfies the Helmholtz equation

$$\nabla ·\left(\frac{1}{\mu \left(\mathbf{x}\right)}\nabla u_{tot}(\mathbf{x};k)\right)+{\omega }^2\epsilon \left(\mathbf{x}\right)u_{tot}(\mathbf{x};k)=0\mathrm{in}\Omega$$

(2)

with transmission condition on the boundary $\partial \Gamma$. Similarly, the incident field $u_{inc}(\mathbf{x};k)$ satisfies the homogeneous Helmholtz equation

$$▵u_{inc}(\mathbf{x};k)+k^2{u}_{inc}(\mathbf{x};k)=0\mathrm{in}\Omega .$$

Throughout this paper, we consider the illumination of plane waves

$$u_{inc}(\mathbf{x};k)=exp(ik\mathbf{\theta }·\mathbf{x})\mathrm{for}\mathbf{x}\in \Omega ,$$

where $\mathbf{\theta }$ is a two-dimensional vector, which characterizes the propagation direction, on the two-dimensional unit circle centered at the origin ${\mathbb{S}}^1$. As is usual, the total field $u_{tot}(\mathbf{x};k)$ divides itself into the incident field $u_{inc}(\mathbf{x};k)$ and the corresponding scattered field $u_{scat}(\mathbf{x};k)$ that satisfies $u_{tot}(\mathbf{x};k)=u_{inc}(\mathbf{x};k)+u_{scat}(\mathbf{x};k)$ and the following Sommerfeld radiation condition

$$\underset{\left|\mathbf{x}\right|\to \infty }{lim}\sqrt{\left|\mathbf{x}\right|}\left(\frac{\partial u_{scat}(\mathbf{x};k)}{\partial \left|\mathbf{x}\right|}-ik{u}_{scat}(\mathbf{x};k)\right)=0\mathrm{uniformly}\mathrm{in}\mathrm{all}\mathrm{directions}\mathbf{\vartheta }=\frac{\mathbf{x}}{\left|\mathbf{x}\right|}.$$

The far-field pattern $u_{\infty }(\mathbf{\vartheta },\mathbf{\theta };k)$ of the scattered field $u_{scat}(\mathbf{x};k)$ is defined on the ${\mathbb{S}}^1$ that satisfies

$$u_{\infty }(\mathbf{\vartheta },\mathbf{\theta };k)=\frac{exp\left(ik\right|\mathbf{y}\left|\right)}{\sqrt{\left|\mathbf{y}\right|}}{u}_{scat}(\mathbf{x};k)+o\left(\frac{1}{\sqrt{\left|\mathbf{y}\right|}}\right)\mathrm{as}\left|\mathbf{y}\right|⟶\infty \mathrm{uniformly}\mathrm{on}\mathbf{\vartheta }=\frac{\mathbf{y}}{\left|\mathbf{y}\right|}.$$

Note that $u_{\infty }(\mathbf{\vartheta },\mathbf{\theta };k)$ can be represented as the single-layer potential with unknown density function $\phi (\mathbf{x},\mathbf{\theta })$:

$$u_{\infty }(\mathbf{\vartheta },\mathbf{\theta };k)=\frac{k^2(1+i)}{4\sqrt{k\pi }}{\int }_{\Gamma }exp(-ik\mathbf{\vartheta }·\mathbf{x})\phi (\mathbf{x},\mathbf{\theta })d\mathbf{x}.$$

(3)

The purpose of this kind of inverse scattering problem is to identify outline shape of $\Gamma$ from a set of measured far-field pattern data $\{u_{\infty }({\mathbf{\vartheta }}_j,{\mathbf{\theta }}_l;k):j,l=1,2,\cdots ,N\}$ with various observation ${\mathbf{\vartheta }}_j$ and incident ${\mathbf{\theta }}_l$ directions. Since the structure of density function $\phi (\mathbf{x},\mathbf{\theta })$ of Equation (3) is unknown, it is very hard to design an algorithm for identifying $\Gamma$ so that we need an approximation form of $\phi (\mathbf{x},\mathbf{\theta })$. Following Beretta and Francini [39], $u_{\infty }(\mathbf{\vartheta },\mathbf{\theta };k)$ can be written as an asymptotic expansion formula.

Lemma 1

(Asymptotic expansion formula). For $\mathbf{\vartheta },\mathbf{\theta }\in {\mathbb{S}}^1$ and $\mathbf{x}\in \Omega \backslash \overline{\Gamma }$, the far-field pattern $u_{\infty }(\mathbf{x},\mathbf{\theta };k)$ can be represented as

$$u_{\infty }(\mathbf{\vartheta },\mathbf{\theta };k)=h\frac{k^2(1+i)}{4\sqrt{k\pi }}{\int }_{\sigma }\left(\frac{{\epsilon }_{\star }-{\epsilon }_0}{\sqrt{{\epsilon }_0{\mu }_0}}-2\mathbf{\vartheta }·\mathbb{M}\left(\mathbf{y}\right)·\mathbf{\theta }\right)exp\left(ik(\mathbf{\theta }-\mathbf{\vartheta })·\mathbf{y}\right)d\mathbf{y}+o\left(h\right).$$

(4)

Here,$\mathbb{M}\left(\mathbf{y}\right)$is a$2\times 2$symmetric matrix defined as follows: let$\mathbf{t}\left(\mathbf{y}\right)$and$\mathbf{n}\left(\mathbf{y}\right)$denote unit tangent and normal vectors to σ at$\mathbf{y}$, respectively. Then,

• $\mathbb{M}\left(\mathbf{y}\right)$has eigenvectors$\mathbf{t}\left(\mathbf{y}\right)$and$\mathbf{n}\left(\mathbf{y}\right)$.
• The eigenvalues corresponding to$\mathbf{t}\left(\mathbf{y}\right)$and$\mathbf{n}\left(\mathbf{y}\right)$are$2\left(\frac{1}{{\mu }_{\star }}-\frac{1}{{\mu }_0}\right)$and$2\left(\frac{1}{{\mu }_0}-\frac{{\mu }_{\star }}{{\mu }_0^2}\right)$, respectively.

2.2. Introduction to Subspace Migration for Imaging of Thin Inhomogeneity

Now, we introduce the traditional imaging function of subspace migration for identifying the shape of $\Gamma$. A detailed description can be found in [30,33]. Let $\mathbb{K}\left(k\right)\in {\mathbb{C}}^{N\times N}$ be the MSR matrix whose elements are $u_{\infty }({\mathbf{\vartheta }}_j,{\mathbf{\theta }}_l;k)$ at observation direction ${\mathbf{\vartheta }}_j$ for the incident direction ${\mathbf{\theta }}_l$, $j,l=1,2,\cdots ,N$, such that

$$\mathbb{K}\left(k\right)=\left(\begin{array}{cccc}{u}_{\infty }({\mathbf{\vartheta }}_1,{\mathbf{\theta }}_1)& u_{\infty }({\mathbf{\vartheta }}_1,{\mathbf{\theta }}_2)& \cdots & u_{\infty }({\mathbf{\vartheta }}_1,{\mathbf{\theta }}_N)\\ u_{\infty }({\mathbf{\vartheta }}_2,{\mathbf{\theta }}_1)& u_{\infty }({\mathbf{\vartheta }}_2,{\mathbf{\theta }}_2)& \cdots & u_{\infty }({\mathbf{\vartheta }}_2,{\mathbf{\theta }}_N)\\ ⋮& ⋮& \ddots & ⋮\\ u_{\infty }({\mathbf{\vartheta }}_N,{\mathbf{\theta }}_1)& u_{\infty }({\mathbf{\vartheta }}_N,{\mathbf{\theta }}_2)& \cdots & u_{\infty }({\mathbf{\vartheta }}_N,{\mathbf{\theta }}_N)\end{array}\right).$$

In this paper, we assume that ${\mathbf{\vartheta }}_j=-{\mathbf{\theta }}_j$, i.e., we have the same incident and observation directions configuration. It is worth emphasizing that, for a given wavenumber $\lambda$ such that $k=2\pi /\lambda$, based on the resolution limit, any detail less than one-half of the wavelength cannot be retrieved. Hence, if we divide thin inhomogeneity $\Gamma$ into M different segments of size of order $\lambda /2$, only one point, say ${\mathbf{x}}_m$, $m=1,2,\cdots ,M$, at each segment will affect the imaging (see [40,41]). If total number of incident and observation directions N is sufficiently large and $N>3M$, the elements of MSR matrix can be represented as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill u_{\infty }({\mathbf{\vartheta }}_j,{\mathbf{\theta }}_l;k)=h\frac{k^2(1+i)}{4\sqrt{k\pi }}& {\int }_{\sigma }\left(\frac{{\epsilon }_{\star }-{\epsilon }_0}{\sqrt{{\epsilon }_0{\mu }_0}}-2\left(\mathbf{\vartheta }·\mathbb{M}\left(\mathbf{y}\right)·\mathbf{\theta }\right)\right)exp\left(ik(\mathbf{\theta }-\mathbf{\vartheta })·\mathbf{y}\right)d\mathbf{y}\hfill \\ \hfill \approx h\frac{k^2(1+i)}{4\sqrt{k\pi }}& \frac{\left|\sigma \right|}{M}{\displaystyle \sum _{m=1}^M}\left[\frac{{\epsilon }_{\star }-{\epsilon }_0}{\sqrt{{\epsilon }_0{\mu }_0}}+\left(\frac{1}{{\mu }_{\star }}-\frac{1}{{\mu }_0}\right)\left({\mathbf{\theta }}_j·\mathbf{t}\left({\mathbf{y}}_m\right)\right)\left({\mathbf{\theta }}_l·\mathbf{t}\left({\mathbf{y}}_m\right)\right)\right.\hfill \\ \hfill & +\left.\left(\frac{1}{{\mu }_0}-\frac{{\mu }_{\star }}{{\mu }_0^2}\right)\left({\mathbf{\theta }}_j·\mathbf{n}\left({\mathbf{y}}_m\right)\right)\left({\mathbf{\theta }}_l·\mathbf{n}\left({\mathbf{y}}_m\right)\right)\right]exp\left(ik({\mathbf{\theta }}_j+{\mathbf{\theta }}_l)·{\mathbf{y}}_m\right),\hfill \end{array}\end{array}$$

(5)

where $\left|\sigma \right|$ denotes the length of $\sigma$.

Based on the representation in Equation (5), $\mathbb{K}\left(k\right)$ can be decomposed as follows:

$$\mathbb{K}\left(k\right)=\mathbb{H}\left(k\right)\mathbb{B}\left(k\right)\overline{\mathbb{H}}\left(k\right)=\sum _{m=1}^{3M}{\tau }_m\left(k\right)\mathbb{H}\left(k\right)\overline{\mathbb{H}}\left(k\right),$$

(6)

where $\mathbb{B}\left(k\right)\in {\mathbb{R}}^{3M\times 3M}$ is a block diagonal matrix with components

$$h\frac{k^2(1+i)\left|\sigma \right|}{4M\sqrt{k\pi }}\left(\begin{array}{cc}{\displaystyle \frac{{\epsilon }_{\star }-{\epsilon }_0}{\sqrt{{\epsilon }_0{\mu }_0}}}& {\mathbb{O}}_{1\times 2}\\ {\mathbb{O}}_{2\times 1}& \mathbb{M}\left({\mathbf{x}}_m\right)\end{array}\right)$$

and $\mathbb{H}\left(k\right)\in {\mathbb{C}}^{N\times 3M}$ is written as

$$\mathbb{H}\left(k\right)=\left({\mathbb{H}}_1\left(k\right),{\mathbb{H}}_2\left(k\right),\cdots ,{\mathbb{H}}_M\left(k\right)\right).$$

Here, ${\mathbb{O}}_{p\times q}$ denotes the $p\times q$ zero matrix and vectors and

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathbb{H}}_m\left(k\right)& =\left({\mathbf{H}}_m^{\left(1\right)}\left(k\right),{\mathbf{H}}_m^{\left(2\right)}\left(k\right),{\mathbf{H}}_m^{\left(3\right)}\left(k\right)\right)\hfill \\ \hfill & =\left(\begin{array}{ccc}exp(ik{\mathbf{\theta }}_1·{\mathbf{x}}_m),& \left({\mathbf{\theta }}_1·\mathbf{t}\left({\mathbf{x}}_m\right)\right)exp(ik{\mathbf{\theta }}_1·{\mathbf{x}}_m),& \left({\mathbf{\theta }}_1·\mathbf{n}\left({\mathbf{x}}_m\right)\right)exp(ik{\mathbf{\theta }}_1·{\mathbf{x}}_m)\\ exp(ik{\mathbf{\theta }}_2·{\mathbf{x}}_m),& \left({\mathbf{\theta }}_2·\mathbf{t}\left({\mathbf{x}}_m\right)\right)exp(ik{\mathbf{\theta }}_2·{\mathbf{x}}_m),& \left({\mathbf{\theta }}_2·\mathbf{n}\left({\mathbf{x}}_m\right)\right)exp(ik{\mathbf{\theta }}_2·{\mathbf{x}}_m)\\ ⋮& ⋮& ⋮\\ exp(ik{\mathbf{\theta }}_N·{\mathbf{x}}_m),& \left({\mathbf{\theta }}_N·\mathbf{t}\left({\mathbf{x}}_m\right)\right)exp(ik{\mathbf{\theta }}_N·{\mathbf{x}}_m),& \left({\mathbf{\theta }}_N·\mathbf{n}\left({\mathbf{x}}_m\right)\right)exp(ik{\mathbf{\theta }}_N·{\mathbf{x}}_m)\end{array}\right).\hfill \end{array}\end{array}$$

(7)

Now, let us perform the Singular Value Decomposition (SVD) of $\mathbb{K}\left(k\right)$

$$\mathbb{K}\left(k\right)=\mathbb{U}\left(k\right)\mathbb{S}\left(k\right)\overline{\mathbb{V}}{\left(k\right)}^T=\sum _{m=1}^N{\rho }_m\left(k\right){\mathbf{U}}_m\left(k\right){\overline{\mathbf{V}}}_m{\left(k\right)}^T\approx \sum _{m=1}^{3M}{\rho }_m\left(k\right){\mathbf{U}}_m\left(k\right){\overline{\mathbf{V}}}_m{\left(k\right)}^T,$$

where ${\rho }_m\left(k\right)$, $m=1,2,\cdots ,3M$, are nonzero singular values such that

$${\rho }_1\left(k\right)\ge {\rho }_2\left(k\right)\ge \cdots \ge {\rho }_{3M}\left(k\right)>0\mathrm{and}{\rho }_m\left(k\right)\approx 0\mathrm{for}m\ge 3M+1,$$

and ${\mathbf{U}}_m\left(k\right)$ and ${\mathbf{V}}_m\left(k\right)$ are left- and right-singular vectors of $\mathbb{K}\left(k\right)$, respectively. Based on the structure of Equation (7), define a test vector $\mathbf{T}(\mathbf{z};k)\in {\mathbb{C}}^{N\times 1}$ as

$$\mathbf{T}(\mathbf{z};k)={\left({\mathbf{c}}_1·(1,{\mathbf{\theta }}_1)exp(ik{\mathbf{\theta }}_1·\mathbf{z}),{\mathbf{c}}_2·(1,{\mathbf{\theta }}_2)exp(ik{\mathbf{\theta }}_2·\mathbf{z}),\cdots ,{\mathbf{c}}_N·(1,{\mathbf{\theta }}_N)exp(ik{\mathbf{\theta }}_N·\mathbf{z})\right)}^T$$

and corresponding unit vector

$$\mathbf{F}(\mathbf{z};k):=\frac{\mathbf{T}(\mathbf{z};k)}{\left|\mathbf{T}\right(\mathbf{z};k\left)\right|},$$

(8)

where the selection of ${\mathbf{c}}_n\in {\mathbb{R}}^3\backslash \left\{\mathbf{0}\right\}$, $n=1,2,\cdots ,N$, is depending on the $\mathbf{t}\left({\mathbf{x}}_m\right)$ and $\mathbf{n}\left({\mathbf{x}}_m\right)$, i.e., shape of $\sigma$ (see [41] for a detailed discussion). Then, for a proper choice of ${\mathbf{c}}_n$, we can observe that

$${\mathbf{U}}_m\left(k\right)=exp\left(i{\gamma }_m^{\left(1\right)}\right)\mathbf{F}({\mathbf{y}}_m;k)\mathrm{and}{\mathbf{V}}_m\left(k\right)\approx exp(-i{\gamma }_m^{\left(2\right)})\mathbf{F}({\mathbf{y}}_m;k),{\gamma }_m^{\left(1\right)}+{\gamma }_m^{\left(2\right)}=arg\left({\tau }_m\left(k\right)\right).$$

Since singular vectors are orthonormal to each others, the first $3M$ columns of $\mathbb{U}\left(k\right)$ and $\mathbb{V}\left(k\right)$ are orthonormal. Furthermore, since ${\rho }_m\left(k\right)\in \mathbb{R}$, it follows that

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \langle \mathbf{F}(\mathbf{z};k),{\mathbf{U}}_m\left(k\right)\rangle \langle \mathbf{F}(\mathbf{z};k),{\overline{\mathbf{V}}}_m\left(k\right)\rangle \approx 1\mathrm{if}\mathbf{z}={\mathbf{y}}_m\hfill \\ \hfill & \langle \mathbf{F}(\mathbf{z};k),{\mathbf{U}}_m\left(k\right)\rangle \langle \mathbf{F}(\mathbf{z};k),{\overline{\mathbf{V}}}_m\left(k\right)\rangle <1\mathrm{if}\mathbf{z}\ne {\mathbf{y}}_m,\hfill \end{array}\end{array}$$

(9)

where $\langle \mathbf{a},\mathbf{b}\rangle =\overline{\mathbf{a}}·\mathbf{b}$ for $\mathbf{a},\mathbf{b}\in \mathbb{C}$.

Hence, we can introduce subspace migration for imaging of thin inhomogeneity at a given wavenumber k as

$$f_{SF}(\mathbf{z};k):=\left|\sum _{m=1}^M\langle \mathbf{F}(\mathbf{z};k),{\mathbf{U}}_m\left(k\right)\rangle \langle \mathbf{F}(\mathbf{z};k),{\overline{\mathbf{V}}}_m\left(k\right)\rangle \right|.$$

(10)

Based on the properties in Equation (9), map of $f_{SF}(\mathbf{z};k)$ should exhibit peaks of magnitude 1 at $\mathbf{z}={\mathbf{x}}_m\in \sigma$, and of small magnitude at $\mathbf{z}\in \Omega \backslash \overline{\Gamma }$. This is the reason thin inhomogeneity can be imaged via subspace migration if a priori information of target is known.

3. Analysis of Subspace Migration without a Priori Information of Inhomogeneity

Based on the above, defining a vector $\mathbf{T}(\mathbf{z};k)$ in Equation (8) is very important to obtain a good result. This means that a proper selection of test vectors ${\mathbf{c}}_n$ seems to be considered beforehand. Notice that, based on Equations (5) and (7), ${\mathbf{c}}_n$ must be a linear combination of tangential $\mathbf{t}\left({\mathbf{x}}_m\right)$ and normal $\mathbf{n}\left({\mathbf{x}}_m\right)$ vectors at ${\mathbf{x}}_m\in \sigma$. Unfortunately, because we have no a priori information of shape of $\Gamma$, it is very hard to define an optimal vector $\mathbf{T}(\mathbf{z};k)$. Thus, ${\mathbf{c}}_n$ is chosen as a fixed vector [41,42] or selected from a set of various directions [33,43].

We now identify the structure of Equation (10) without consideration of shape of $\Gamma$. Since we have no information of $\mathbf{t}\left(\mathbf{y}\right)$ and $\mathbf{n}\left(\mathbf{y}\right)$ for $\mathbf{y}\in \sigma$, we cannot select ${\mathbf{c}}_n$ of Equation (8) so that by neglecting ${\mathbf{c}}_n$ such that ${\mathbf{c}}_n·(1,{\mathbf{\theta }}_n)=1$ for all n, we consider the following test vector

$$\mathbf{W}(\mathbf{z};k)=\frac{1}{\sqrt{N}}{\left(exp(ik{\mathbf{\theta }}_1·\mathbf{z}),exp(ik{\mathbf{\theta }}_2·\mathbf{z}),\cdots ,exp(ik{\mathbf{\theta }}_N·\mathbf{z})\right)}^T,$$

(11)

and consider the corresponding single-frequency imaging function of subspace migration

$$f_{SUB}(\mathbf{z};k):=\left|\sum _{m=1}^M\langle \mathbf{W}(\mathbf{z};k),{\mathbf{U}}_m\left(k\right)\rangle \langle \mathbf{W}(\mathbf{z};k),{\overline{\mathbf{V}}}_m\left(k\right)\rangle \right|.$$

(12)

For starting analysis, we introduce two useful identities derived in [33].

Lemma 2.

Let$\mathbf{x}={\left|\mathbf{x}\right|(cos\varphi ,sin\varphi )}^T,\mathbf{\xi }=\left|\mathbf{\xi }\right|{(cos\xi ,sin\xi )}^T\in {\mathbb{R}}^2$, and${\mathbf{\theta }}_n={(cos{\theta }_n,sin{\theta }_n)}^T\in {\mathbb{S}}^1$,$n=1,2,\cdots ,N$. Then, the following relations hold uniformly:

$$\frac{1}{N}\sum _{n=1}^{N}exp(ik{\mathbf{\theta }}_n·\mathbf{x})=J_0\left(k\right|\mathbf{x}\left|\right)+\frac{2}{\pi }\sum _{s=1}^{\infty }\frac{i^s}{s}{J}_s\left(k\right|\mathbf{x}\left|\right)cos\frac{s({\theta }_N+{\theta }_1-2\varphi )}{2}sin\frac{s({\vartheta }_N-{\vartheta }_1)}{2}$$

and

$$\begin{array}{c}\hfill \frac{1}{N}\sum _{n=1}^N({\mathbf{\theta }}_n·\mathbf{\xi })exp(ik{\mathbf{\theta }}_n·\mathbf{x})=i\left(\frac{\mathbf{x}}{\left|\mathbf{x}\right|}·\mathbf{\xi }\right)J_1\left(k\right|\mathbf{x}\left|\right)+\frac{1}{\pi }{J}_0\left(k\right|\mathbf{x}\left|\right)sin\frac{{\vartheta }_N-{\vartheta }_1}{2}cos\frac{{\theta }_N+{\theta }_1-2\xi }{2}\\ \hfill +\frac{i}{2\pi }{J}_1\left(k\right|\mathbf{x}\left|\right)sin({\theta }_N-{\theta }_1)cos({\theta }_N+{\theta }_1-\xi -\varphi )\\ \hfill +\frac{1}{\pi }\sum _{s=2}^{\infty }{i}^s{J}_s\left(k\right|\mathbf{x}\left|\right)[\frac{1}{1-s}sin\frac{(1-s)({\theta }_N-{\theta }_1)}{2}cos\frac{(1-n)({\theta }_N+{\theta }_1)+2s\varphi -2\xi }{2}\\ \hfill +\frac{1}{1+s}sin\frac{(1+s)({\theta }_N-{\theta }_1)}{2}cos\frac{(1+s)({\theta }_N+{\theta }_1)-2s\varphi -2\xi }{2}],\end{array}$$

where$J_s(·)$denotes the Bessel function of integer order s of the first kind.

3.1. Analysis of Single-Frequency Imaging Function

Based on Equation (11) and Lemma 2, we obtain the following main result.

Theorem 1.

Let${\mathbf{x}}_m-\mathbf{z}=|{\mathbf{x}}_m-\mathbf{z}|{(cos{\varphi }_m,sin{\varphi }_m)}^T$,$\mathbf{t}\left({\mathbf{x}}_m\right)={(cos{\xi }_m^{\left(1\right)},sin{\xi }_m^{\left(1\right)})}^T$,$\mathbf{n}\left({\mathbf{x}}_m\right)={(cos{\xi }_m^{\left(2\right)},sin{\xi }_m^{\left(2\right)})}^T$, and${\mathbf{\theta }}_n={(cos{\theta }_n,sin{\theta }_n)}^T\in {\mathbb{S}}^1$. Then, for sufficiently large$N>3M$,$f_{SUB}(\mathbf{z};k)$can be represented as follows:

$$\begin{array}{c}{f}_{SUB}(\mathbf{z};k)=\left|\sum _{m=1}^M{\left(J_0\left(k\right|{\mathbf{x}}_m-\mathbf{z}\left|\right)+{\Psi }_1(\mathbf{z};k)\right)}^2\right.\hfill \\ \hfill \left.+\sum _{m=1}^M{\left(\sqrt{2}i\left(\frac{{\mathbf{x}}_m-\mathbf{z}}{|{\mathbf{x}}_m-\mathbf{z}|}·\left(\mathbf{t}\left({\mathbf{x}}_m\right)+\mathbf{n}\left({\mathbf{x}}_m\right)\right)\right)J_1\left(k\right|{\mathbf{x}}_m-\mathbf{z}\left|\right)+\frac{1}{\pi }\sum _{r=1}^2{\Psi }_2^{\left(r\right)}(\mathbf{z};k)\right)}^2\right|,\end{array}$$

(13)

where

$${\Psi }_1(\mathbf{z};k)=\frac{2}{\pi }\sum _{s=1}^{\infty }\frac{i^s}{s}{J}_s\left(k\right|{\mathbf{x}}_m-\mathbf{z}\left|\right)cos\frac{s({\theta }_N+{\theta }_1-2{\varphi }_m)}{2}sin\frac{s({\theta }_N-{\theta }_1)}{2}$$

and

$$\begin{array}{cc}\hfill {\Psi }_2^{\left(r\right)}(\mathbf{z};k)=& J_0\left(k\right|{\mathbf{x}}_m-\mathbf{z}\left|\right)sin\frac{{\theta }_N-{\theta }_1}{2}cos\frac{{\theta }_N+{\theta }_1-2{\xi }_m^{\left(r\right)}}{2}\hfill \\ \hfill & +\frac{i}{2}{J}_1\left(k\right|{\mathbf{x}}_m-\mathbf{z}\left|\right)sin({\theta }_N-{\theta }_1)cos\left({\theta }_N+{\theta }_1-{\xi }_m^{\left(r\right)}-\varphi \right)\hfill \\ \hfill & +\sum _{s=2}^{\infty }{i}^s{J}_s\left(k\right|{\mathbf{x}}_m-\mathbf{z}\left|\right)[\frac{1}{1-s}sin\frac{(1-s)({\theta }_N-{\theta }_1)}{2}cos\frac{(1-n)({\theta }_N+{\theta }_1)+2s{\varphi }_m-2{\xi }_m^{\left(r\right)}}{2}\hfill \\ \hfill & +\frac{1}{1+s}sin\frac{(1+s)({\theta }_N-{\theta }_1)}{2}cos\frac{(1+s)({\theta }_N+{\theta }_1)-2s{\varphi }_m-2{\xi }_m^{\left(r\right)}}{2}].\hfill \end{array}$$

Furthermore, if$k⟶+\infty$, then

$$f_{SUB}(\mathbf{z};k)=\delta ({\mathbf{x}}_m-\mathbf{z}),$$

where δ is the Dirac delta function.

Proof. 

Since N is large, ${\mathbf{H}}_m^{\left(1\right)}$, ${\mathbf{H}}_m^{\left(2\right)}$ and ${\mathbf{H}}_m^{\left(3\right)}$ of Equation (7) are orthogonal to each other for all $m=1,2,\cdots ,M$ (see Appendix A). Then, by applying Equation (11) and Lemma 2, we can examine

$$f_{SUB}(\mathbf{z};k)=\left|\sum _{m=1}^{3M}〈\mathbf{W}(\mathbf{z};k),{\mathbf{U}}_m〉〈\mathbf{W}(\mathbf{z};k),{\overline{\mathbf{V}}}_m〉\right|\approx \left|\sum _{m=1}^M\sum _{r=1}^3{〈\mathbf{W}(\mathbf{z};k),{\widehat{\mathbf{H}}}_m^{\left(r\right)}\left(k\right)〉}^2\right|,$$

where

$$\begin{array}{cc}\hfill {\widehat{\mathbf{H}}}_m^{\left(1\right)}\left(k\right)& =\frac{{\mathbf{H}}_m^{\left(1\right)}\left(k\right)}{\left|\right|{\mathbf{H}}_m^{\left(1\right)}\left(k\right)\left|\right|}=\frac{1}{\sqrt{N}}{\mathbf{H}}_m^{\left(1\right)}\left(k\right),\hfill \\ \hfill {\widehat{\mathbf{H}}}_m^{\left(2\right)}\left(k\right)& =\frac{{\mathbf{H}}_m^{\left(2\right)}\left(k\right)}{\left|\right|{\mathbf{H}}_m^{\left(2\right)}\left(k\right)\left|\right|}=\frac{\sqrt{2}}{\sqrt{N}}{\mathbf{H}}_m^{\left(2\right)}\left(k\right),\hfill \\ \hfill {\widehat{\mathbf{H}}}_m^{\left(3\right)}\left(k\right)& =\frac{{\mathbf{H}}_m^{\left(3\right)}\left(k\right)}{\left|\right|{\mathbf{H}}_m^{\left(3\right)}\left(k\right)\left|\right|}=\frac{\sqrt{2}}{\sqrt{N}}{\mathbf{H}}_m^{\left(3\right)}\left(k\right).\hfill \end{array}$$

Then, we can evaluate

$$\begin{array}{c}〈\mathbf{W}(\mathbf{z};k),{\widehat{\mathbf{H}}}_m^{\left(1\right)}\left(k\right)〉=\frac{1}{N}\sum _{n=1}^{N}exp\left(ik{\mathbf{\theta }}_n·({\mathbf{x}}_m-\mathbf{z})\right)\hfill \\ \hfill =J_0\left(k\right|{\mathbf{x}}_m-\mathbf{z}\left|\right)+\frac{2}{\pi }\sum _{s=1}^{\infty }\frac{i^s}{s}{J}_s\left(k\right|{\mathbf{x}}_m-\mathbf{z}\left|\right)cos\frac{s({\theta }_N+{\theta }_1-2{\varphi }_m)}{2}sin\frac{s({\theta }_N-{\theta }_1)}{2},\end{array}$$

(14)

$$\begin{array}{c}\begin{array}{c}〈\mathbf{W}(\mathbf{z};k),{\widehat{\mathbf{H}}}_m^{\left(2\right)}\left(k\right)〉=\frac{\sqrt{2}}{N}\sum _{n=1}^N\left({\mathbf{\theta }}_n·\mathbf{t}\left({\mathbf{x}}_m\right)\right)exp\left(ik{\mathbf{\theta }}_n·({\mathbf{x}}_m-\mathbf{z})\right)\hfill \\ =\sqrt{2}i\left(\frac{{\mathbf{x}}_m-\mathbf{z}}{|{\mathbf{x}}_m-\mathbf{z}|}·\mathbf{t}\left({\mathbf{x}}_m\right)\right)J_1\left(k\right|{\mathbf{x}}_m-\mathbf{z}\left|\right)+\frac{1}{\pi }{J}_0\left(k\right|{\mathbf{x}}_m-\mathbf{z}\left|\right)sin\frac{{\theta }_N-{\theta }_1}{2}cos\frac{{\theta }_N+{\theta }_1-2{\xi }_m^{\left(1\right)}}{2}\hfill \\ +\frac{i}{2\pi }{J}_1\left(k\right|{\mathbf{x}}_m-\mathbf{z}\left|\right)sin({\theta }_N-{\theta }_1)cos\left({\theta }_N+{\theta }_1-{\xi }_m^{\left(1\right)}-\varphi \right)\hfill \\ +\frac{1}{\pi }\sum _{s=2}^{\infty }{i}^s{J}_s\left(k\right|{\mathbf{x}}_m-\mathbf{z}\left|\right)[\frac{1}{1-s}sin\frac{(1-s)({\theta }_N-{\theta }_1)}{2}cos\frac{(1-n)({\theta }_N+{\theta }_1)+2s{\varphi }_m-2{\xi }_m^{\left(1\right)}}{2}\hfill \\ +\frac{1}{1+s}sin\frac{(1+s)({\theta }_N-{\theta }_1)}{2}cos\frac{(1+s)({\theta }_N+{\theta }_1)-2s{\varphi }_m-2{\xi }_m^{\left(1\right)}}{2}],\hfill \end{array}\end{array}$$

(15)

and

$$\begin{array}{c}\begin{array}{c}〈\mathbf{W}(\mathbf{z};k),{\widehat{\mathbf{H}}}_m^{\left(3\right)}\left(k\right)〉=\frac{\sqrt{2}}{N}\sum _{n=1}^N\left({\mathbf{\theta }}_n·\mathbf{n}\left({\mathbf{x}}_m\right)\right)exp\left(ik{\mathbf{\theta }}_n·({\mathbf{x}}_m-\mathbf{z})\right)\hfill \\ =\sqrt{2}i\left(\frac{{\mathbf{x}}_m-\mathbf{z}}{|{\mathbf{x}}_m-\mathbf{z}|}·\mathbf{n}\left({\mathbf{x}}_m\right)\right)J_1\left(k\right|{\mathbf{x}}_m-\mathbf{z}\left|\right)+\frac{1}{\pi }{J}_0\left(k\right|{\mathbf{x}}_m-\mathbf{z}\left|\right)sin\frac{{\theta }_N-{\theta }_1}{2}cos\frac{{\theta }_N+{\theta }_1-2{\xi }_m^{\left(2\right)}}{2}\hfill \\ +\frac{i}{2\pi }{J}_1\left(k\right|{\mathbf{x}}_m-\mathbf{z}\left|\right)sin({\theta }_N-{\theta }_1)cos\left({\theta }_N+{\theta }_1-{\xi }_m^{\left(2\right)}-\varphi \right)\hfill \\ +\frac{1}{\pi }\sum _{s=2}^{\infty }{i}^s{J}_s\left(k\right|{\mathbf{x}}_m-\mathbf{z}\left|\right)[\frac{1}{1-s}sin\frac{(1-s)({\theta }_N-{\theta }_1)}{2}cos\frac{(1-n)({\theta }_N+{\theta }_1)+2s{\varphi }_m-2{\xi }_m^{\left(2\right)}}{2}\hfill \\ +\frac{1}{1+s}sin\frac{(1+s)({\theta }_N-{\theta }_1)}{2}cos\frac{(1+s)({\theta }_N+{\theta }_1)-2s{\varphi }_m-2{\xi }_m^{\left(2\right)}}{2}].\hfill \end{array}\end{array}$$

(16)

Hence, Equation (13) can be derived by Equations (14)–(16).

Now, assume that $k⟶+\infty$. Then, it is clear that $f_{SUB}(\mathbf{z};k)=1$ when $\mathbf{z}={\mathbf{x}}_m$. If $\mathbf{z}\ne {\mathbf{x}}_m$, then the following asymptotic form of Bessel function holds for $k|{\mathbf{x}}_m-\mathbf{z}|\gg |s^2-0.25|$,

$$J_s\left(k\right|{\mathbf{x}}_m-\mathbf{z}\left|\right)\approx \sqrt{\frac{2}{k\pi |{\mathbf{x}}_m-\mathbf{z}|}}cos\left\{k|{\mathbf{x}}_m-\mathbf{z}|-\frac{s\pi }{2}-\frac{\pi }{4}+O\left(\frac{1}{k|{\mathbf{x}}_m-\mathbf{z}|}\right)\right\}⟶0,$$

where s denotes a positive integer. Based on this asymptotic form, we can easily observe that $f_{SUB}(\mathbf{z};k)=0$ when $\mathbf{z}\ne {\mathbf{x}}_m$. Hence,

$$f_{SUB}(\mathbf{z};k)=\delta ({\mathbf{x}}_m-\mathbf{z}).$$

This completes the proof. □

Based on the identified structure in Equation (13), we can examine certain properties of $f_{SUB}(\mathbf{z};k)$:

Property 1.

The terms${\Psi }_1(\mathbf{z};k)$and${\Psi }_2^{\left(r\right)}(\mathbf{z};k)$cause the generation of unexpected artifacts so that they must be eliminated to obtain a good result. The simplest method is to apply even, large number N and select the set of incident directions$\left\{{\mathbf{\theta }}_n:n=1,2,\cdots ,N\right\}$to satisfy${\theta }_N-{\theta }_1=2\pi$. This is the theoretical reason an even, large number of incident and observation directions is selected and uniformly distributed on${\mathbb{S}}^1$in most studies. If this condition is satisfied,$f_{SUB}(\mathbf{z};k)$can be written by

$$f_{SUB}(\mathbf{z};k)=\left|\sum _{m=1}^M(J_0\left(k\right|{\mathbf{x}}_m-{\mathbf{z}\left|\right)}^2-2{\left(\frac{{\mathbf{x}}_m-\mathbf{z}}{|{\mathbf{x}}_m-\mathbf{z}|}·\left(\mathbf{t}\left({\mathbf{x}}_m\right)+\mathbf{n}\left({\mathbf{x}}_m\right)\right)\right)}^2{J}_1\left(k\right|{\mathbf{x}}_m-\mathbf{z}{\left|\right)}^2\right|,$$

(17)

Property 2.

Based on recent work [41], the dominant eigenvectors of matrix$\mathbb{M}\left(\mathbf{x}\right)$are$\mathbf{t}\left(\mathbf{x}\right)$and$\mathbf{n}\left(\mathbf{x}\right)$for${\mu }_{\star }<{\mu }_0$and${\mu }_{\star }>{\mu }_0$, respectively. Generally (and based on our assumption), since${\mu }_{\star }>{\mu }_0$,$f_{SUB}(\mathbf{z};k)$of Equation (17) can be written as

$$f_{SUB}(\mathbf{z};k)\approx \left|\sum _{m=1}^M\left\{J_0\left(k\right|{\mathbf{x}}_m-{\mathbf{z}\left|\right)}^2-2{\left(\frac{{\mathbf{x}}_m-\mathbf{z}}{|{\mathbf{x}}_m-\mathbf{z}|}·\mathbf{n}\left({\mathbf{x}}_m\right)\right)}^2{J}_1\left(k\right|{\mathbf{x}}_m-\mathbf{z}{\left|\right)}^2\right\}\right|.$$

Property 3.

Since$J_0\left(0\right)=1$and$J_n\left(0\right)=0$for$n=1,2,\cdots ,$the terms

$$J_0\left(k\right|{\mathbf{x}}_m-{\mathbf{z}\left|\right)}^2\mathrm{and}2{\left(\frac{{\mathbf{x}}_m-\mathbf{z}}{|{\mathbf{x}}_m-\mathbf{z}|}·\mathbf{n}\left({\mathbf{x}}_m\right)\right)}^2{J}_1\left(k\right|{\mathbf{x}}_m-\mathbf{z}{\left|\right)}^2$$

contribute to and disturb the detection of Γ, respectively.

Property 4.

Based on the properties of$J_0\left(x\right)$and$J_1\left(x\right)$, plots of$f_{SUB}(\mathbf{z};k)$will show peaks of magnitude 1 at$\mathbf{z}={\mathbf{x}}_m\in \Gamma$and a small one at$\mathbf{z}\notin \Gamma$. Note that, since$J_1{\left(x\right)}^2$has its maximum value at two points, say$x_1$and$x_2$, symmetric with respect to$x=0$, two curves with large (but less than 1) magnitude and many artifacts with small magnitude will included in the map of$f_{SUB}(\mathbf{z};k)$. This tells us that one can recognize the shape of thin inhomogeneities via the map of$f_{SUB}(\mathbf{z};k)$without a priori information of the shape of the thin inhomogeneities.

Property 5.

If one can increase the value of k, it will be possible to recognize the shape of thin inhomogeneities very clearly. This means that application of high frequency will guarantee a good result but due to the oscillating property of Bessel function, several artifacts will be included in the imaging result also. Note that$f_{SUB}(\mathbf{z};k)=1$at$\mathbf{z}\in \Gamma$for small value of k but due to the spreading effect of oscillation of$J_0$, a result with very low resolution will be obtained. Figure 2 shows the oscillation pattern of the Bessel function related to the applied frequency.

It is worth mentioning that, based on Property 3, $f_{SUB}(\mathbf{z};k)$ has its maximum value at $\mathbf{z}={\mathbf{x}}_m\in \Gamma$. Hence, we can immediately examine following result of unique determination.

Corollary 1.

Let the applied frequency ω be sufficiently high. If the total number N of incident and observation directions is sufficiently large, then the shape of supporting curve σ of thin inclusion Γ can be obtained uniquely via the map of$f_{SUB}(\mathbf{z};k)$.

3.2. Improvement of Imaging Performance Part 1: Filtering

Now, we consider the method of improvement. Based on the structure in Equation (13), we can examine that good results can be obtained via the map of $f_{SUB}(\mathbf{z};k)$ when $k⟶+\infty$. However, this is a theoretically ideal situation. Hence, to obtain good results, an alternative method must be considered.

Based on Property 4, one can obtain good results by eliminating two curves with large magnitude and correspondingly many artifacts with small magnitude. Hence, to design a filtering strategy, let us consider the maximum value of the following term:

$$2{\left(\frac{{\mathbf{x}}_m-\mathbf{z}}{|{\mathbf{x}}_m-\mathbf{z}|}·\mathbf{n}\left({\mathbf{x}}_m\right)\right)}^2{J}_1\left(k\right|{\mathbf{x}}_m-\mathbf{z}{\left|\right)}^2.$$

Since $J_1{\left(x\right)}^2$ has its maximum at $x\approx 0.147$, based on the numerical computation,

$$\underset{\mathbf{z}\in \Omega }{max}{J}_1\left(k\right|{\mathbf{x}}_m-\mathbf{z}{\left|\right)}^2\approx {(0.58186522\dots )}^2\approx 0.338567\mathrm{and}2{\left(\frac{{\mathbf{x}}_m-\mathbf{z}}{|{\mathbf{x}}_m-\mathbf{z}|}·\mathbf{n}\left({\mathbf{x}}_m\right)\right)}^2\le 2,$$

we can say that

$$2{\left(\frac{{\mathbf{x}}_m-\mathbf{z}}{|{\mathbf{x}}_m-\mathbf{z}|}·\mathbf{n}\left({\mathbf{x}}_m\right)\right)}^2{J}_1\left(k\right|{\mathbf{x}}_m-\mathbf{z}{\left|\right)}^2\approx 0.677134268\dots <0.678.$$

This means that the magnitude of two curves in the neighborhood of $\Gamma$ will be less than $0.678$ (see Figure 3). Hence, let us introduce a filtering function ${\mathfrak{F}}_{SF}$

$${\mathfrak{F}}_{SF}\left(x\right)=\left\{\begin{array}{ccc}x& \mathrm{if}& 0.678\le x\le 1\\ 0& \mathrm{if}& 0\le x<0.678.\end{array}\right.$$

Then, better imaging results of thin inhomogeneity can be obtained via the map of ${\mathfrak{F}}_{SF}\left(f_{SUB}(\mathbf{z};k)\right)$. Note that this method can be applied in the imaging of single inhomogeneity (see Figure 9).

Remark 1.

Theoretically, the maximum value of$f_{SUB}(\mathbf{z};k)$is equal to 1 but, generally, the maximum value is smaller than 1 in the results of numerical simulations (see Figure 9). There are many reasons for this, e.g., the algorithm is based on the asymptotic expansion formula, the influence of random noise, and computational errors. Thus, instead of ${\mathfrak{F}}_{SF}\left(f_{SUB}(\mathbf{z};k)\right)$, we introduce the following normalized imaging function

$$f_{NSUB}(\mathbf{z};k):=\frac{f_{SUB}(\mathbf{z};k)}{{\displaystyle \underset{\mathbf{z}\in \Omega }{max}\left|f_{SUB}(\mathbf{z};k)\right|}}$$

and consider the filtered map ${\mathfrak{F}}_{SF}\left(f_{NSUB}(\mathbf{z};k)\right)$ for identifying shape of inclusions.

3.3. Improvement of Imaging Performance Part 2: Application of Multi-Frequency

The authors of [30,33] confirmed that applying multi-frequency offers better results than applying single frequency. However this fact holds with an optimal choice of ${\mathbf{c}}_n$ in Equation (8). Now, we consider the multi-frequency subspace migration and examine whether it improves single-frequency one. Let $\{k_f=2\pi /{\lambda }_f:f=1,2,\cdots ,F\}$ be the set of $F-$different wavenumbers such that

$$k_1<k_2<\cdots <k_F,\mathrm{i}.\mathrm{e}.,{\lambda }_1>{\lambda }_2>\cdots >{\lambda }_F,$$

and let $\mathbb{K}\left(k_f\right)$ be the collected MSR matrix at $k_f$. Then, by performing SVD of $\mathbb{K}\left(k_f\right)$ as

$$\mathbb{K}\left(k_f\right)=\mathbb{U}\left(k_f\right)\mathbb{S}\left(k_f\right)\overline{\mathbb{V}}{\left(k_f\right)}^T=\sum _{m=1}^N{\rho }_m\left(k_f\right){\mathbf{U}}_m\left(k_f\right){\overline{\mathbf{V}}}_m{\left(k_f\right)}^T\approx \sum _{m=1}^{3M_f}{\rho }_m\left(k_f\right){\mathbf{U}}_m\left(k_f\right){\overline{\mathbf{V}}}_m{\left(k_f\right)}^T,$$

we can introduce multi-frequency subspace migration,

$$f_{MF}(\mathbf{z};F):=\frac{1}{F}\left|\sum _{f=1}^F{f}_{SF}(\mathbf{z};k)\right|=\frac{1}{F}\left|\sum _{f=1}^F\sum _{m=1}^{M_f}\langle \mathbf{F}(\mathbf{z};k_f),{\mathbf{U}}_m\left(k_f\right)\rangle \langle \mathbf{F}(\mathbf{z};k_f),{\overline{\mathbf{V}}}_m\left(k_f\right)\rangle \right|,$$

(18)

and the corresponding multi-frequency subspace migration without a priori information,

$$f_{MSUB}(\mathbf{z};F):=\frac{1}{F}\left|\sum _{f=1}^F{f}_{SUB}(\mathbf{z};k_f)\right|=\frac{1}{F}\left|\sum _{f=1}^F\sum _{m=1}^{M_f}\langle \mathbf{W}(\mathbf{z};k_f),{\mathbf{U}}_m\left(k_f\right)\rangle \langle \mathbf{W}(\mathbf{z};k_f),{\overline{\mathbf{V}}}_m\left(k_f\right)\rangle \right|,$$

(19)

where $\mathbf{F}(\mathbf{z};k_f)$ and $\mathbf{W}(\mathbf{z};k_f)$ are defined in Equations (8) and (11), respectively.

Based on the results in several works [29,30,37,38], Equation (18) is an improved imaging function compared to Equation (10). Opposite to the statistical approach, we identify the reason of improvement by establishing a relationship with Bessel functions of integer order as follows.

Theorem 2.

For sufficiently large$N(>3M)$and k,$f_{MSUB}(\mathbf{z};k)$can be represented as follows:

$$\begin{array}{c}{f}_{MSUB}(\mathbf{z};F)\approx \frac{1}{k_F-k_1}\left|\sum _{m=1}^M\{\Phi ({\mathbf{x}}_m-\mathbf{z};k_F)-\Phi ({\mathbf{x}}_m-\mathbf{z};k_1)\right.\hfill \\ \hfill +\left.\left(1-2{\left(\frac{{\mathbf{x}}_m-\mathbf{z}}{|{\mathbf{x}}_m-\mathbf{z}|}·\mathbf{n}\left({\mathbf{x}}_m\right)\right)}^2\right){\int }_{k_1}^{k_F}{J}_1\left(k\right|{\mathbf{x}}_m-\mathbf{z}{\left|\right)}^{2}dk\}\right|,\end{array}$$

(20)

where

$$\Phi (x;k):=k\left(J_0{\left(k\right|x\left|\right)}^2+J_1{\left(k\right|x\left|\right)}^2\right).$$

Proof. 

For the sake of simplicity, we assume that $M_f=M$ for $f=1,2,\cdots ,F$, i.e., the difference ${\lambda }_1-{\lambda }_F$ is small enough (in the numerical experiments, we set ${\lambda }_1=0.6$ and ${\lambda }_F=0.3$, i.e., ${\lambda }_1-{\lambda }_F=0.3$ is small enough; see Section 3). Then, based on the structure in Equation (13), we can say that

$$\begin{array}{c}{f}_{MSUB}(\mathbf{z};F)\approx \frac{1}{F}\left|\sum _{f=1}^F\sum _{m=1}^M\left\{J_0(k_f|{\mathbf{x}}_m-{\mathbf{z}\left|\right)}^2-{\left(\frac{{\mathbf{x}}_m-\mathbf{z}}{|{\mathbf{x}}_m-\mathbf{z}|}·\left(\mathbf{t}\left({\mathbf{x}}_m\right)+\mathbf{n}\left({\mathbf{x}}_m\right)\right)\right)}^2{J}_1(k_f|{\mathbf{x}}_m-\mathbf{z}{\left|\right)}^2\right\}\right|\hfill \\ \hfill \approx \frac{1}{k_F-k_1}\left|\sum _{m=1}^M\left\{{\int }_{k_1}^{k_F}{J}_0\left(k\right|{\mathbf{x}}_m-{\mathbf{z}\left|\right)}^{2}dk-{\left(\frac{{\mathbf{x}}_m-\mathbf{z}}{|{\mathbf{x}}_m-\mathbf{z}|}·\left(\mathbf{t}\left({\mathbf{x}}_m\right)+\mathbf{n}\left({\mathbf{x}}_m\right)\right)\right)}^2{\int }_{k_1}^{k_F}{J}_1\left(k\right|{\mathbf{x}}_m-\mathbf{z}{\left|\right)}^{2}dk\right\}\right|.\end{array}$$

Since the following relation holds for $x\in \mathbb{R}$,

$$\int J_0{\left(x\right)}^{2}dx=x\left(J_0{\left(x\right)}^2+J_1{\left(x\right)}^2\right)+\int J_1{\left(x\right)}^{2}dx,$$

we can evaluate

$$\begin{array}{c}{\int }_{k_1}^{k_F}{J}_0\left(k\right|{\mathbf{x}}_m-\mathbf{z}{\left|\right)}^{2}dk=k_F\left(J_0(k_F|{\mathbf{x}}_m-{\mathbf{z}\left|\right)}^2+J_1(k_F|{\mathbf{x}}_m-\mathbf{z}{\left|\right)}^2\right)\hfill \\ \hfill -k_1\left(J_0(k_1|{\mathbf{x}}_m-{\mathbf{z}\left|\right)}^2+J_1(k_1|{\mathbf{x}}_m-\mathbf{z}{\left|\right)}^2\right)+{\int }_{k_1}^{k_F}{J}_1\left(k\right|{\mathbf{x}}_m-\mathbf{z}{\left|\right)}^{2}dk.\end{array}$$

with which we can obtain Equation (20). □

Now, let us compare the results in Theorems 1 and 2. Based on the structures in Equations (13) and (20), imaging functionals are composed with contributing and disturbing terms for imaging. First, the contributing terms of Equations (13) and (20) are

$$J_0\left(k\right|{\mathbf{x}}_m-\mathbf{z}\left|\right)\mathrm{and}\frac{1}{k_F-k_1}\left(\Phi ({\mathbf{x}}_m-\mathbf{z};k_F)-\Phi ({\mathbf{x}}_m-\mathbf{z};k_1)\right),$$

respectively. Since $\Phi ({\mathbf{x}}_m-\mathbf{z};k_F)-\Phi ({\mathbf{x}}_m-\mathbf{z};k_1)$ oscillates less than $J_0\left(k\right|{\mathbf{x}}_m-\mathbf{z}\left|\right)$ (see [33]), the identified shape of the supporting curve via the map $f_{MSUB}(\mathbf{z};F)$ will be better than the one via the map $f_{SUB}(\mathbf{z};k)$. Next, the disturbing terms of Equations (13) and (20) are

$$2{\left(\frac{{\mathbf{x}}_m-\mathbf{z}}{|{\mathbf{x}}_m-\mathbf{z}|}·\mathbf{n}\left({\mathbf{x}}_m\right)\right)}^2{J}_1\left(k\right|{\mathbf{x}}_m-{\mathbf{z}\left|\right)}^2\mathrm{and}\left(1-2{\left(\frac{{\mathbf{x}}_m-\mathbf{z}}{|{\mathbf{x}}_m-\mathbf{z}|}·\mathbf{n}\left({\mathbf{x}}_m\right)\right)}^2\right){\int }_{k_1}^{k_F}{J}_1\left(k\right|{\mathbf{x}}_m-\mathbf{z}{\left|\right)}^{2}dk,$$

respectively. Similar to the comparison of contributing terms, we can observe that, since the term

$${\int }_{k_1}^{k_F}{J}_1\left(k\right|{\mathbf{x}}_m-\mathbf{z}{\left|\right)}^{2}dk$$

oscillates less than $J_1\left(k\right|{\mathbf{x}}_m-\mathbf{z}{\left|\right)}^2$ and the factor

$$1-2{\left(\frac{{\mathbf{x}}_m-\mathbf{z}}{|{\mathbf{x}}_m-\mathbf{z}|}·\mathbf{n}\left({\mathbf{x}}_m\right)\right)}^2$$

will reduce magnitude of disturbing term, the disturbing term of Equation (20) will affect imaging performance less than the one of Equation (13). Thus, we can conclude Corollary 2.

Corollary 2.

Maps of multi-frequency subspace migration$f_{MSUB}(\mathbf{z};F)$yields better imaging results owing to less oscillation than single-frequency one$f_{SUB}(\mathbf{z};k)$. This means that unexpected artifacts in the map of$f_{MSUB}(\mathbf{z};F)$are mitigated when F is sufficiently large.

As in the single-frequency case, we can immediately conclude that the following result of uniqueness holds.

Corollary 3.

Suppose that the values of$k_f$are sufficiently high. If the total number N of incident and observation directions and total number F of applied frequencies are sufficiently large, then the shape of supporting curve σ of thin inclusion Γ can be obtained uniquely via the map of$f_{MSUB}(\mathbf{z};F)$.

Remark 2

(Filtering). Similar to the case of single-frequency, let us consider the method of filtering. Since

$$\underset{\mathbf{z}\in \Omega }{max}\left|1-2{\left(\frac{{\mathbf{x}}_m-\mathbf{z}}{|{\mathbf{x}}_m-\mathbf{z}|}·\mathbf{n}\left({\mathbf{x}}_m\right)\right)}^2\right|=1,$$

and

$$\frac{1}{F}\sum _{f=1}^F{J}_1(k_f|\mathbf{z}-{\mathbf{x}}_m{\left|\right)}^2\approx \frac{1}{k_F-k_1}{\int }_{k_1}^{k_F}{J}_1\left(k\right|\mathbf{z}-{\mathbf{x}}_m{\left|\right)}^{2}dk\le 0.338567\approx 0.340$$

for multi-frequency imaging, we can define a filtering function such that

$$\mathfrak{F}\left(x\right)=\left\{\begin{array}{ccc}x& \mathit{if}& 0.340\le x\le 1\\ 0& \mathit{if}& 0\le x<0.340.\end{array}\right.$$

Then,$\mathfrak{F}\left(f_{MSUB}(\mathbf{z};F)\right)$will be an improved version of$f_{MSUB}(\mathbf{z};F)$.

Remark 3.

Similar to Remark 1, in this paper, we consider the following normalized multi-frequency imaging function

$$f_{NMSUB}(\mathbf{z};F):=\frac{f_{MSUB}(\mathbf{z};F)}{{\displaystyle \underset{\mathbf{z}\in \Omega }{max}\left|f_{MSUB}(\mathbf{z};F)\right|}}$$

and consider the filtered map$\mathfrak{F}\left(f_{NMSUB}(\mathbf{z};F)\right)$instead of$f_{MSUB}(\mathbf{z};F)$.

4. Results of Numerical Simulations

In this section, we exhibit some results of numerical simulations to support Theorems 1 and 2. To describe thin inclusions ${\Gamma }_j$, $j=1,2$, two supporting smooth curves are selected as follows:

$$\begin{array}{cc}\hfill {\sigma }_1& =\{{(s-0.2,-0.5s^2+0.4)}^T:-0.5\le s\le 0.5\}\hfill \\ \hfill {\sigma }_2& =\{{(s+0.2,s^3+s^2-0.6)}^T:-0.5\le s\le 0.5\}.\hfill \end{array}$$

The thickness h of thin inclusions ${\Gamma }_j$ is equally set to $0.015$ and the imaging domain $\Omega$ is selected as $\Omega =[-1,1]\times [-1,1]$. We denote ${\epsilon }_j$ and ${\mu }_j$ as the permittivity and permeability of ${\Gamma }_j$, respectively, and set parameters ${\mu }_j$, ${\mu }_0$, ${\epsilon }_j$, and ${\epsilon }_0$ as 5, 1, 5, and 1, respectively. Since ${\mu }_0$ and ${\epsilon }_0$ are set to unity, the applied frequencies reads as $k_f=2\pi /{\lambda }_f$ at wavelength ${\lambda }_f$ for $f=1,2,\cdots ,F(=10)$, which is varied in the numerical examples between ${\lambda }_1=0.6$ and ${\lambda }_{10}=0.3$. For single frequency imaging, $k=2\pi /0.4$ is applied. Total number of incident and observation directions is set to $N=32$ and

$${\mathbf{\vartheta }}_n=-{\mathbf{\theta }}_n=\left(cos\frac{2(n-1)\pi }{N},sin\frac{2(n-1)\pi }{N}\right),n=1,2,\cdots ,N.$$

Figure 4 shows an illustration of ${\Gamma }_j$ and ${\mathbf{\vartheta }}_n$. The far-field pattern data $u_{\infty }({\mathbf{\vartheta }}_j,{\mathbf{\theta }}_l;k)$ is generated by solving a second-kind Fredholm integral equation along the supporting curve, linked to the thinness of the inclusion [44]. To show robustness, 15 dB Gaussian random noise is added to the unperturbed data.

Example 1.

First, let us examine the effect of the selection${\mathbf{c}}_n$of Equation (8). Figure 5 shows maps of$f_{SF}(\mathbf{z};k)$for${\mathbf{c}}_n={(1,1,0)}^T$,${\mathbf{c}}_n={(1,cos\left(45{}^{\circ }\right),sin\left(45{}^{\circ }\right))}^T$, and${\mathbf{c}}_n={(1,cos\left(30{}^{\circ }\right),sin\left(30{}^{\circ }\right))}^T$when the thin inhomogeneity is${\Gamma }_1$. Throughout the result, we can observe that one cannot identify the shape of${\Gamma }_1$at this moment. Note that, based on Property 1, the dominant eigenvectors are$\mathbf{n}\left({\mathbf{x}}_m\right)$, thus${\mathbf{c}}_n$must be of the form${(1,\mathbf{n}{({\mathbf{x}}_m)}^T)}^T$. Hence, from now on, we apply${\mathbf{c}}_n={(1,0,1)}^T$for$f_{SF}(\mathbf{z};k)$.

Example 2.

Now, let us consider the imaging results of$f_{SF}(\mathbf{z};k)$,$f_{SUB}(\mathbf{z};k)$, and$\mathfrak{F}\left(f_{NSUB}(\mathbf{z};k)\right)$. On the basis of the results in Figure 6, we can observe that the shape of${\Gamma }_1$can be recognized via the maps of$f_{SF}(\mathbf{z};k)$and$f_{SUB}(\mathbf{z};k)$but the imaging seems rather coarse for the traditional method$f_{SF}(\mathbf{z};k)$, better for the proposed one$f_{SUB}(\mathbf{z};k)$. Furthermore,$\mathfrak{F}\left(f_{NSUB}(\mathbf{z};k)\right)$exhibits very accurate shape of${\Gamma }_1$so that suggested filtering method seems very effective.

Example 3.

Here, we examine the effect of total number of incident and observation direction N for imaging of$f_{SUB}(\mathbf{z};k)$. From the results in Figure 7a, we can see that the shape of${\Gamma }_1$cannot be recognized when N is even but small. If one wants to recognize the shape, as discussed for Property 1, N must be increased, but, if N is an odd, the map of$f_{SUB}(\mathbf{z};k)$contains some unexpected artifacts (see Figure 7b). To remove the artifacts, a large N must be applied (see Figure 7c), and this result supports the theoretical result and Property 1.

Example 4.

Maps of$f_{SF}(\mathbf{z};k)$,$f_{SUB}(\mathbf{z};k)$, and$\mathfrak{F}\left(f_{NSUB}(\mathbf{z};k)\right)$are shown in Figure 8 when the thin inhomogeneity is${\Gamma }_2$. Similar to the results in Figure 6, the shape of${\Gamma }_1$can be recognized but the obtained shape of${\Gamma }_1$via$f_{SUB}(\mathbf{z};k)$looks like an anchor due to the appearance of unexpected artifacts. Although some artifacts are still visible, the result via$f_{SUB}(\mathbf{z};k)$seems better than the one via$f_{SF}(\mathbf{z};k)$. Similar to the previous result, the filtering method still seems effective.

Example 5.

It is well-known that one of advantage of subspace migration is its straightforward application to the imaging of multiple inhomogeneities. Figure 9 shows the maps of$f_{SF}(\mathbf{z};k)$,$f_{SUB}(\mathbf{z};k)$, and$\mathfrak{F}\left(f_{NSUB}(\mathbf{z};k)\right)$for imaging multiple thin inhomogeneities${\Gamma }_1\cup {\Gamma }_2$with the same permittivity${\epsilon }_1={\epsilon }_2=5$and permeability${\mu }_1={\mu }_2=5$. Similar to the imaging of single inhomogeneity, we can observe that$f_{SUB}(\mathbf{z};k)$is an effective method. However, due to the appearance of artifacts, it is hard to identify true shape of inhomogeneities. Thus, in contrast to the imaging of single inhomogeneity, filtering method is not effective for imaging of multiple inhomogeneities.

Example 6.

Figure 10 shows the maps of$f_{SF}(\mathbf{z};k)$,$f_{SUB}(\mathbf{z};k)$, and$\mathfrak{F}\left(f_{NSUB}(\mathbf{z};k)\right)$under the same configuration as the previous result in Figure 9, except for different material properties,${\epsilon }_1={\mu }_1=10$and${\epsilon }_2={\mu }_2=5$. Then, based on the results in [31,41], the values of$f_{SF}(\mathbf{z};k)$and$f_{SUB}(\mathbf{z};k)$for$\mathbf{z}\in {\Gamma }_2$are smaller than$f_{SF}(\mathbf{z};k)$and$f_{SUB}(\mathbf{z};k)$for$\mathbf{z}\in {\Gamma }_1$, respectively. Furthermore, due to the unexpected artifacts, the shape of${\Gamma }_2$cannot be identified, while${\Gamma }_1$can be identified. Correspondingly, only the shape of${\Gamma }_1$can be identified in the map of$\mathfrak{F}\left(f_{NSUB}(\mathbf{z};k)\right)$.

Example 7.

From now on, we consider the multi-frequency imaging. Maps of$f_{MF}(\mathbf{z};F)$,$f_{MSUB}(\mathbf{z};F)$, and$\mathfrak{F}\left(f_{NMSUB}(\mathbf{z};F)\right)$are exhibited in Figure 11 when the thin inhomogeneity is${\Gamma }_1$. Although the map of$f_{MSUB}(\mathbf{z};F)$contains more artifacts than$f_{MF}(\mathbf{z};F)$, the identified shape via the map of$f_{MSUB}(\mathbf{z};F)$seems close to the true shape of${\Gamma }_1$. It is interesting to observe that, opposite to Remark 2, the unexpected peak of small (but cnon-negligible) magnitude remains in the neighborhood of the tip of${\Gamma }_1$but the result seems very nice. A similar phenomenon can be examined through the results in Figure 12 when the thin inhomogeneity is${\Gamma }_2$.

Example 8.

Figure 13 shows the maps of$f_{MF}(\mathbf{z};F)$,$f_{MSUB}(\mathbf{z};F)$, and$\mathfrak{F}\left(f_{NMSUB}(\mathbf{z};F)\right)$for imaging multiple thin inhomogeneities${\Gamma }_1\cup {\Gamma }_2$with same permittivity${\epsilon }_1={\epsilon }_2=5$and permeability${\mu }_1={\mu }_2=5$. Similar to the single-frequency imaging, it seems that$f_{MSUB}(\mathbf{z};F)$performs better imaging accomplishment than$f_{MF}(\mathbf{z};F)$. However, due to the appearance of artifacts, it is hard to identify true shape of inhomogeneities. Furthermore, in contrast to the single-frequency imaging, the filtering method seems very effective for imaging, although some peaks of small magnitudes remain.

Example 9.

Figure 14 shows the maps of$f_{MF}(\mathbf{z};F)$,$f_{MSUB}(\mathbf{z};F)$, and$\mathfrak{F}\left(f_{NMSUB}(\mathbf{z};F)\right)$under the same configuration as the previous result in Figure 13, except for different material properties,${\epsilon }_1={\mu }_1=10$and${\epsilon }_2={\mu }_2=5$. Similar to the previous example,$f_{MSUB}(\mathbf{z};F)$can be regarded as an improved version of$f_{MF}(\mathbf{z};F)$.

Example 10.

For the final example, let us consider the application of filtering for imaging of multiple inhomogeneities when their permittivities and permeabilities are different to each other. A simple method is to divide the search domain into two (or more) disjoint areas and applying the filtering function to each area. Figure 15 shows corresponding results. By comparing the results in Figure 10 and Figure 12, the identified shapes are more accurate than the traditional ones.

5. Conclusions

In this paper, we consider the subspace migration imaging functional without any a priori information of thin inhomogeneities. We derive a relationship between the imaging functional and the Bessel functions of integer order of the first kind. By comparing with traditional research, the results obtained are not considered to be better than the ones of previous studies conducted with a priori information of inhomogeneities. Nevertheless, the derived results indicate that, although some unexpected artifacts appear, it is possible to identify the outline shape of thin electromagnetic inhomogeneities without any a priori information. Furthermore, based on the established theoretical results, the designed filtering technique will be effective to remove unnecessary artifacts in the imaging results. For a further improvement, a multi-frequency based imaging algorithm is also suggested and successfully applied.

Here, we focus on the imaging of thin, curve-like electromagnetic inhomogeneities. In the same line of thought, the analysis of subspace migration for the imaging of perfectly conducting cracks in Transverse Magnetic (TM) and Transverse Electric (TE) cases would be an interesting research topic. In this paper, we consider an imaging of a two-dimensional thin electromagnetic inclusions. In three-dimensional microwave imaging, the measurement data (scattering parameter) in the presence of small object $\Sigma$ can be approximated by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill S_{\mathrm{scat}}(m,n)& =\frac{i{k}^2}{4\omega {\mu }_0}{\int }_{\Omega }\left(\frac{\epsilon \left(\mathbf{r}\right)-{\epsilon }_0}{{\epsilon }_0}+i\frac{\sigma \left(\mathbf{r}\right)-{\sigma }_0}{\omega {\epsilon }_0}\right){\mathbf{E}}_{inc}({\mathbf{r}}_n,\mathbf{r})·{\mathbf{E}}_{tot}(\mathbf{r},{\mathbf{r}}_m)d\mathbf{r}\hfill \\ \hfill & \approx \frac{i{k}^2\mathrm{area}(\Sigma )}{4\omega {\mu }_0}\left(\frac{{\epsilon }_{\star }-{\epsilon }_0}{{\epsilon }_0}+i\frac{{\sigma }_{\star }-{\sigma }_0}{\omega {\epsilon }_0}\right){\mathbf{E}}_{\mathrm{inc}}({\mathbf{r}}_n,{\mathbf{r}}_{\star })·{\mathbf{E}}_{\mathrm{inc}}({\mathbf{r}}_m,{\mathbf{r}}_{\star }),\hfill \end{array}\end{array}$$

(21)

where ${\mathbf{E}}_{tot}\in {\mathbb{C}}^3$ and ${\mathbf{E}}_{inc}\in {\mathbb{C}}^3$ denote the total and incident fields, respectively; ${\mathbf{r}}_n$ is the nth signal transmitter or receiver; and ${\sigma }_{\star }$ and ${\sigma }_0$ are conductivities of $\Sigma$ and $\Omega$, respectively (see [35,36,45]. As we can see, Equation (21) has some similarities to Equation (4). Hence, we expect that the analysis could be extended to a three-dimensional inverse scattering problem and microwave imaging. Following Park [32], subspace migration is successfully applied to the half-space problem. The application to the detection of anti-personnel mines buried in the ground or cracks in concrete walls or bridges would be also an interesting research subject.