A novel study on the bifocusing method for imaging unknown objects in two-dimensional inverse scattering problem

1.   Introduction

In this paper, we consider the localization of a set of small two-dimensional (2D) homogeneous inhomogeneities using various electromagnetic source transmitters and receivers. The problem is formulated as an inverse scattering problem for the 2D electromagnetic system. To solve this problem, various remarkable techniques have been investigated, most of which are classified as quantitative or qualitative methods.

In order to reconstruct the parameter (such as dielectric permittivity, electric conductivity, or magnetic permeability at a given frequency) distribution or identify an unknown object larger than the given wavelength, various quantitative methods (or iterative methods) have been developed, such as Newton's method ^[47,62], Gauss–Newton method ^[1,16], optimal control approach ^[3], Levenberg–Marquardt algorithm ^[13,27], Born iterative method ^[33,50], level-set technique ^[23,64], and subspace optimization method ^[18,69]. Although these algorithms have shown feasibility and applicability in inverse scattering problems, good initial guess must be generated before the iteration procedure because the success of iterative-based techniques significantly depend on the initial guess ^[48,61].

Alternatively, for obtaining a good initial guess or identifying the location or determining the outline shape of small objects, various qualitative methods (or non-iterative methods) have been investigated, such as the direct sampling method ^[34,37,39], factorization method ^[15,20,45], Kirchhoff and subspace migrations ^[5,54,59], MUltiple SIgnal Classification (MUSIC) algorithm ^[22,42,58], linear sampling method ^[9,19,46], topological derivative strategy ^[4,51,53], orthogonality sampling method ^[31,36,60], and variational algorithm based on the Fourier inversion ^[6,7,8].

The bifocusing method (BFM) is classified as a qualitative method using a multistatic measurement system. Recently, it has been successfully applied to various interesting problems, such as damage detection of inhomogeneous structures ^[26,44], radar imaging ^[35], and anomaly detection in microwave imaging ^[66]. {Let us emphasize that a reliable mathematical theory of BFM in the presence of small anomaly has been developed only for the transverse magnetic polarization (TM) case, refer to ^[43]. However, to the best of our knowledge, the application of BFM in inverse scattering problem for both transverse magnetic and electric (TE) polarizations is still heuristic; thus, developing an appropriate mathematical theory and performing numerical simulations with synthetic and experimental data are necessary to explain the applicability, fundamental limitations, and various intrinsic properties of the BFM in TM and TE polarizations.

In this paper, we design the imaging functions of the BFM for the dielectric permittivity (TM polarization) and magnetic permeability (TE polarization) contrast cases in full-aperture inverse scattering problem. Based on the asymptotic expansion formula in the presence of small inhomogeneities, we show that the designed imaging functions comprise the Bessel function of order zero and order one for the dielectric permittivity and magnetic permeability contrast cases, respectively, and the material properties (such as the total number, size, dielectric permittivity, and magnetic permeability) of the inhomogeneities. This analysis confirms the reason of the unique determination of inhomogeneities for the permittivity contrast case and the identification of two locations instead of true locations of inhomogeneities for the permeability contrast case. To demonstrate the theoretical results, various results of the numerical simulation with synthetic and experimental data are presented.

This study is organized as follows. In Section 2, we introduce the basic concept of the 2D direct scattering problem in the presence of a set of small, well-separated inhomogeneities and the asymptotic expansion formula for the scattered field. In Section 3, we design the imaging functions of the BFM for the dielectric permittivity and magnetic permeability contrast cases, explore the mathematical structures of the imaging functions, and discuss their applicability, limitations, and various properties. In Section 4, we present some numerical simulation results with synthetic and experimental data, and compared with those obtained by other qualitative methods such as direct sampling method, MUSIC, subspace migration, and topological derivative method. In Section 5, we present a short conclusion including an outline of future research.

Finally, let us emphasize that the imaging result via experimental data obtained through the BFM does not guarantee the exact shape of object. However, retrieved result can be adopted as a good initial guess and correspondingly, it will be possible to retrieve a better shape through the iterative-based algorithms. Let us mention some classical but very useful approaches of the retrieval of objects via Fresnel data; Bayesian approach ^[10], modified Born method ^[12], contrast source inversion method ^[14], conjugate gradient fast Fourier transform (CG-FFT) method ^[21], modified gradient method ^[24], image fusion approach ^[25], diffraction tomography, modified gradient, and real-coded genetic algorithms ^[52], level set method ^[63], linear spectral estimation technique ^[67], and distorted-wave Born approach ^[68].

2.   Forward problem and asymptotic expansion formula

We shall consider 2D direct scattering from a set of small electromagnetic inhomogeneities located in homogeneous space $\mathbb{R}^2$. We denote $\Omega$ be the region of interest (ROI) and assume that the inhomogeneity $\Sigma_s$, $s = 1, 2, \cdots, S$, is a small ball with radius $\alpha_s$ and location $\mathbf{x}_s$:

where $\mathbb{S}^1$ denotes the unit circle centered at the origin.

Let $\varepsilon_\mathrm{b}$ and $\mu_\mathrm{b}$ denote the values of the dielectric permittivity and magnetic permeability of $\mathbb{R}^2$, respectively, at the given angular frequency $\omega = 2\pi f$, where $f$ denotes the ordinary frequency measured in $\texttt{hertz}$. Similarly, we denote $\varepsilon_s$ and $\mu_s$ as the values of the dielectric permittivity and magnetic permeability of $\Sigma_s$. With this, we denote $\varepsilon(\mathbf{x})$ and $\mu(\mathbf{x})$ as the piecewise constants of dielectric permittivity and magnetic permeability, respectively, such that

With this, we denote $k$ as the background wavenumber that satisfies $k^2 = \omega^2 \varepsilon_\mathrm{b} \mu_\mathrm{b}$. Throughout this paper, we assume that all $\Sigma_s$ are well-separated from each other and, correspondingly, that $k$ satisfies the following:

We assume that transmitters are located at $\mathbf{t}_m$, $m = 1, 2, \cdots, M$, while receivers are located at $\mathbf{r}_n$, $n = 1, 2, \cdots, N$. That is, $\mathbf{t}_m$ and $\mathbf{r}_n$ are given by the following:

and

respectively. In this study, we consider the point-source illumination, letting

be the incident field at source $\mathbf{t}_m$. Here, $J_n$ and $Y_n$ denote the Bessel and Neumann functions of integer order $n$, respectively, and $H_0^{(1)}$ denotes the Hankel function or order zero of the first kind. Throughout this paper, we assume that the time dependence is of the form $e^{-i\omega t}$. Then, the governing equation for the time-harmonic total field $u(\mathbf{t}_m, \mathbf{x})$ emitted by a source at $\mathbf{t}_m$ is

with transmission conditions at the boundaries of $\Sigma_s$. With this, we denote $u_{{\rm{scat}}}(\mathbf{t}_m, \mathbf{r}_n)$ as the scattered field, which is obtained by subtracting the total and incident fields and satisfies the Sommerfeld radiation condition

uniformly in all directions $\mathbf{x}/| \mathbf{x}|$. Based on ^[7], $u_{{\rm{scat}}}(\mathbf{t}_m, \mathbf{r}_n)$ can be represented as an asymptotic expansion formula, which plays a key role in designing the imaging function.

Lemma 2.1 (Asymptotic formula). For sufficiently large $k$, $u_{{\rm{scat}}}(\mathbf{t}_m, \mathbf{r}_n)$ can be represented as follows:

where $\mathbb{M}(\mathbf{x}_s)$ denotes $2\times2$ diagonal matrix with diagonal element $2\pi \mu_\mathrm{b}/(\mu_s+ \mu_\mathrm{b})$.

3.   Designing imaging functions and qualitative analysis

In this section, we apply the asymptotic formula (2.5) to design an imaging function for identifying $\Sigma_s$ with no prior information and to establish a mathematical theory. Consequently, we examine the structure of the MSR matrix $\mathbb{K}$. The analysis is then led successively for two cases: dielectric permittivity ($\varepsilon(\mathbf{x})\ne \varepsilon_\mathrm{b}$ and $\mu(\mathbf{x}) = \mu_\mathrm{b}$) and magnetic permeability contrast cases ($\varepsilon(\mathbf{x}) = \varepsilon_\mathrm{b}$ and $\mu(\mathbf{x})\ne \mu_\mathrm{b}$). Throughout this section, we assume that there is no noise and the full MSR matrix $\mathbb{K}$ is available such that

3.1.   Dielectric permittivity contrast (transverse magnetic polarization) case

First, we assume that $\varepsilon(\mathbf{x})\ne \varepsilon_\mathrm{b}$ and $\mu(\mathbf{x}) = \mu_\mathrm{b}$. Based on the representation formula (2.5), $u_{{{\rm{scat}}}}^{(\varepsilon)}(\mathbf{t}_m, \mathbf{r}_n)$ can be written as follows:

Correspondingly, the MSR matrix becomes

Based on the above representation, we can design the imaging function forms every image point through the synthesis inversion of two Green's functions, $G(\mathbf{t}_m, \mathbf{x})$ and $G(\mathbf{r}_n, \mathbf{x})$, for $m = 1, 2, \cdots, M$ and $n = 1, 2, \cdots, N$. Correspondingly, the imaging function $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x})$ of BFM can be introduced as follows: for $\mathbf{x}\in\Omega$,

Here, $\Omega$ denotes the imaging area. The resulting plot of $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x})$ is expected to exhibit large peaks at $\mathbf{x}_s\in\Sigma_s$, $s = 1, 2, \cdots, S$ so that all inhomogeneities can be recognized through the map of $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x})$, refer to ^[35].

Although the location of the inhomogeneities can be retrieved from the map of $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x})$, further mathematical theories to explain the feasibilities, fundamental limitations, and some phenomena in the simulation results must be investigated. Here, we show that $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x})$ can be written by the Bessel function of order zero of the first kind. For a proper derivation, we introduce the following useful identity derived in ^[56].

Lemma 3.1. For sufficiently large $M$, $\boldsymbol{\theta}_m, \boldsymbol{\theta}\in\mathbb{S}^1$, and $\mathbf{x}\in\mathbb{R}^2$, the following relation holds

Then, we can obtain the following result.

Theorem 3.1. Let $\mathbf{x}\in\Omega$, $\boldsymbol{\theta}_m = (\cos\theta_m, \sin\theta_m)$, $\boldsymbol{\vartheta}_n = (\cos\vartheta_n, \sin\vartheta_n)$, and $\mathbf{x}- \mathbf{x}_s = | \mathbf{x}- \mathbf{x}_s|(\cos\phi_s, \sin\phi_s)$, and assume that the values of $M$, $N$, $k$ are sufficiently large such that $| \mathbf{x}- \mathbf{x}_s|\gg0.25/k$. Subsequently, for $\mathbf{t}_m$ and $\mathbf{r}_n$ respectively defined in (2.2) and (2.3), $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x})$ can be represented as follows:

Proof. Since $| \mathbf{t}_m- \mathbf{x}|, | \mathbf{r}_n- \mathbf{x}|\gg0.25/k$ and $-2\pi < \text{arg}(k| \mathbf{t}_m- \mathbf{x}|), \text{arg}(k| \mathbf{r}_n- \mathbf{x}|) < \pi$ for all $m$ and $n$, the following asymptotic forms hold

and

Since

we have

Thus, $\Phi_{{\rm{TM}}}(\mathbf{x})$ can be written as follows:

Since $M$ and $N$ are sufficiently large, applying Lemma 3.1 yields

Hence, we can obtain

and correspondingly, (3.2) can be derived. □

Based on the identified structure (3.2), we can examine some properties of the imaging function $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x})$.

Remark 3.1 (Unique detectability and resolution). Since $J_0(0) = 1$, $\mathbf{x}_s$ can be identified uniquely because the local maxima of the $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x})$ will appear at $\mathbf{x}_s$. Moreover, several artifacts that depend on $k$ also appear in the map of $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x})$ due to the oscillating property of the $J_0$. That is, if the value of $k$ is sufficiently large, high-resolution results can be obtained, but several artifacts will be included in the map of $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x})$. Conversely, if the value of $k$ is small, a blurred image will appear, but less artifacts will be included. Thus, $\mathbf{x}_s$ cannot be identified using the map of $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x})$ if the value $k$ is extremely small (see Figure 11).

Remark 3.2 (Material properties). Based on (3.2), we can examine the value of $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x})$ as proportional to $\alpha_s^2(\varepsilon_s- \varepsilon_\mathrm{b})$. This means that the value of $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x})$ depends on the size and permittivity of the inhomogeneities. That is, if there are two inhomogeneities $\Sigma_1$ and $\Sigma_2$ satisfying $\alpha_1^2(\varepsilon_1- \varepsilon_\mathrm{b}) > \alpha_2^2(\varepsilon_2- \varepsilon_\mathrm{b})$, then $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x}_1) > \mathfrak{F}_{{\rm{TM}}}(\mathbf{x}_2)$ so that the identification of $\mathbf{x}_1$ is guaranteed. However, the identification of $\mathbf{x}_2$ is not guaranteed in some cases because the value of $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x}_2)$ cannot be distinguished from the artifacts. Hence, if the size or permittivity of an inhomogeneity is larger than that of another inhomogeneity, its location can be identified. If $\alpha_1^2(\varepsilon_1- \varepsilon_\mathrm{b})\gg\alpha_2^2(\varepsilon_2- \varepsilon_\mathrm{b})$, then the location $\mathbf{x}_1$ can be identified, while $\mathbf{x}_2$ cannot.

Based on the Remark 3.1, we can examine the following result of unique determination.

Corollary 3.1 (Unique determination in permittivity contrast case). With the same condition of the Theorem 3.1, the location of small inhomogeneities can be identified uniquely through the map of $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x})$.

3.2.   Magnetic permeability contrast (transverse electric polarization) case

Next, we assume that $\varepsilon(\mathbf{x}) = \varepsilon_\mathrm{b}$ and $\mu(\mathbf{x})\ne \mu_\mathrm{b}$. Based on the representation formula (2.5), $u_{{\rm{scat}}}^{(\mu)}(\mathbf{t}_m, \mathbf{r}_n)$ can be written as follows:

Correspondingly, the MSR matrix becomes

Based on the above representation, we can regard that the imaging function forms every image point through the synthesis inversion of $\nabla G(\mathbf{t}_m, \mathbf{x})\cdot\nabla G(\mathbf{r}_n, \mathbf{x})$ for $m = 1, 2, \cdots, M$ and $n = 1, 2, \cdots, N$. Correspondingly, the imaging function $f(\mathbf{x})$ can be introduced as follows: for $\mathbf{x}\in\Omega$,

Remark 3.3 Notably, since

sometimes $m, n\in\mathbb{N}$ and $\mathbf{x}\in\Omega$ exist such that

This means that the denominator of $\Phi(\mathbf{x})$ defined in (3.4) can be zero, we refer to Figure 1 for an illustration.

Based on the Remark 3.3, using the $\mathfrak{F}(\mathbf{x})$ of (3.4) to identify $\mathbf{x}_s$ is inappropriate. Thus, motivated by the permittivity contrast case, we consider the following imaging function of the BFM: for $\mathbf{x}\in\Omega$,

In order to explore some properties of $\mathfrak{F}_{{\rm{TE}}}(\mathbf{x})$, we explore its structure. To this end, we introduce the following useful identity derived in ^[56].

Lemma 3.2. For sufficiently large $M$, $\boldsymbol{\theta}_m, \boldsymbol{\theta}, \boldsymbol{\vartheta}_n, \in\mathbb{S}^1$, and $\mathbf{x}\in\mathbb{R}^2$, the following relation holds

Then, we can obtain the following result.

Theorem 3.2. Let $\mathbf{x}\in\Omega$, $\boldsymbol{\theta}_m = (\cos\theta_m, \sin\theta_m)$, $\boldsymbol{\vartheta}_n = (\cos\vartheta_n, \sin\vartheta_n)$, and $\mathbf{x}- \mathbf{x}_s = | \mathbf{x}- \mathbf{x}_s|(\cos\phi_s, \sin\phi_s)$, and assume that the values of $M$, $N$, $k$ are sufficiently large such that $| \mathbf{x}- \mathbf{x}_s|\gg0.25/k$. Then, for $\mathbf{t}_m$ and $\mathbf{r}_n$ respectively defined in (2.2) and (2.3), $\mathfrak{F}_{{\rm{TE}}}(\mathbf{x})$ can be represented as follows:

Proof. Since $| \mathbf{t}_m- \mathbf{x}|\gg0.25/k$, $| \mathbf{r}_n- \mathbf{x}|\gg0.25/k$, $-2\pi < \text{arg}(k| \mathbf{t}_m- \mathbf{x}|) < \pi$, and $-2\pi < \text{arg}(k| \mathbf{r}_n- \mathbf{x}|) < \pi$ for all $m$ and $n$, we have

Subsequently, using (3.3) and (3.7), $\Phi_{{\rm{TE}}}(\mathbf{x})$ can be written as follows:

Since $M$ is sufficiently large, applying Lemma 3.2 yields

Since $N$ is sufficiently large, applying Lemma 3.2 again, we can obtain

Therefore, by combining (3.9) and (3.10), $\Phi_{{\rm{TE}}}(\mathbf{x})$ of (3.8) becomes

and correspondingly, (3.6) can be derived. □

Based on the identified structure (3.6), we can examine some properties of the imaging function $\mathfrak{F}_{{\rm{TE}}}(\mathbf{x})$.

Remark 3.4 (Detectability and resolution). Since $J_1(0) = 0$, contrary to the permittivity contrast case, $\mathbf{x}_s$ cannot be identified through the map of $\mathfrak{F}_{{\rm{TE}}}(\mathbf{x})$. Instead, since $J_1(|x|)$ has its maximum value at $x\approx\pm1.8412$, the map of $\mathfrak{F}_{{\rm{TE}}}(\mathbf{x})$ has a large magnitude ring at each location $\mathbf{x}_s$ satisfies $k| \mathbf{x}- \mathbf{x}_s| = 1.8412$. Therefore, it can be confirmed that the center of the ring is the location of $\mathbf{x}_s$. Similar to the permittivity contrast case, several artifacts that depend on the value of $k$ are also included in the map of $\mathfrak{F}_{{\rm{TE}}}(\mathbf{x})$ due to the oscillation property of $J_1$.

If an extremely large value of $k$ is applied, since $| \mathbf{x}- \mathbf{x}_s| = 1.8412/k\longrightarrow0$, the ring with large magnitude is concentrated to the center $\mathbf{x}_s$ so that an almost accurate location $\mathbf{x}_s$ can be identified through the map of $\mathfrak{F}_{{\rm{TE}}}(\mathbf{x})$, refer to Figure 15 at $f = {10}{{\rm{GHz}}}$ and $f = {12}{{\rm{GHz}}}$.

Remark 3.5 (Material properties). Based on (3.6), we can examine the value of $\mathfrak{F}_{{\rm{TE}}}(\mathbf{x})$ as proportional to $\alpha_s^2/(\mu_s+ \mu_\mathrm{b})$. This means that the value of $\mathfrak{F}_{{\rm{TE}}}(\mathbf{x})$ depends on the size and permeability of the inhomogeneities. That is, if there are two inhomogeneities $\Sigma_1$ and $\Sigma_2$ satisfying $\alpha_1^2/(\mu_1+ \mu_\mathrm{b}) > \alpha_2^2/(\mu_2+ \mu_\mathrm{b})$, then $\mathfrak{F}_{{\rm{TE}}}(\mathbf{x}_1) > \mathfrak{F}_{{\rm{TE}}}(\mathbf{x}_2)$ so that the identification of $\mathbf{x}_1$ is guaranteed. However, the identification of $\mathbf{x}_2$ is not guaranteed in some cases because the value of $\mathfrak{F}_{{\rm{TE}}}(\mathbf{x}_2)$ cannot be distinguished from the artifacts. Hence, if the size of an inhomogeneity is larger than that of another inhomogeneity or the permeability of an inhomogeneity is smaller than that of another inhomogeneity, its location can be identified. If $\alpha_1^2/(\mu_1+ \mu_\mathrm{b})\gg\alpha_2^2/(\mu_2+ \mu_\mathrm{b})$, then the location $\mathbf{x}_1$ can be identified, while $\mathbf{x}_2$ cannot.

3.3.   Identification of accurate location in magnetic permeability contrast case

Although, it can be recognized the existence of small inhomogeneities through the appearance of rings in the map of $\mathfrak{F}_{{\rm{TE}}}(\mathbf{x})$, it is still difficult to identity their true locations. Here, we design a new imaging function of the BFM for accurate localization of small inhomogeneities. The main idea is to convert the factor $J_1(k| \mathbf{x}- \mathbf{x}_s|)^2$ of (3.6) to the form $J_0(k| \mathbf{x}- \mathbf{x}_s|)$. To this end, we consider the following result derived in ^[57].

Lemma 3.3. For sufficiently large $M$, $\boldsymbol{\theta}_m, \boldsymbol{\theta}\in\mathbb{S}^1$, and $\mathbf{x}, \mathbf{y}\in\mathbb{R}^2$, the following relation holds

On the basis of Lemma 3.3, we introduce the following imaging function

Then, the location of small inhomogeneities can be retrieved through the map of $\mathfrak{F}_{{\rm{ITE}}}(\mathbf{x})$. To support this assertion, we derive the following result.

Theorem 3.3. Let $\mathbf{x}\in\Omega$, $\boldsymbol{\theta}_m = (\cos\theta_m, \sin\theta_m)$, $\boldsymbol{\vartheta}_n = (\cos\vartheta_n, \sin\vartheta_n)$, and $\mathbf{x}- \mathbf{x}_s = | \mathbf{x}- \mathbf{x}_s|(\cos\phi_s, \sin\phi_s)$, and assume that the values of $M$, $N$, $k$ are sufficiently large such that $| \mathbf{x}- \mathbf{x}_s|\gg0.25/k$. Then, for $\mathbf{t}_m$ and $\mathbf{r}_n$ respectively defined in (2.2) and (2.3), $\mathfrak{F}_{{\rm{ITE}}}(\mathbf{x})$ can be represented as follows:

where

Proof. Based on the (3.3) and (3.7), $\Phi_{{\rm{ITE}}}(\mathbf{x})$ can be written as follows:

Since $M$ is sufficiently large, $\boldsymbol{\theta}_m\cdot(\mathbf{x}- \mathbf{x}_s) = | \mathbf{x}- \mathbf{x}_s|\cos\big(\theta_m-\phi_s)$, and $\boldsymbol{\theta}_m\cdot \boldsymbol{\vartheta}_n = \cos(\theta_m-\vartheta_n)$, applying Lemma 3.3 yields the following:

Hence,

Since $N$ is sufficiently large, applying Lemma 3.3 again, we can derive

Therefore, by combining (3.14) and (3.15), $\Phi_{{\rm{ITE}}}(\mathbf{x})$ of (3.13) becomes

and correspondingly, (3.12) can be derived. □

Based on the identified structure (3.6), we can examine some properties of the imaging function $\mathfrak{F}_{{\rm{TE}}}(\mathbf{x})$.

Remark 3.6 (Detectability). Since $\mathfrak{F}_{{\rm{ITE}}}(\mathbf{x})$ contains the factor $J_0(k| \mathbf{x}- \mathbf{x}_s|)^2$, similar to the permittivity contrast case, $\mathbf{x}_s$ can be identified through the map of $\mathfrak{F}_{{\rm{ITE}}}(\mathbf{x})$. Notice that the factor $J_2(k| \mathbf{x}- \mathbf{x}_s|)^2$ is also included in the $\mathfrak{F}_{{\rm{ITE}}}(\mathbf{x})$ but since (see ^[49] for instance)

it generates some artifacts with small magnitudes. Numerically, $J_2(x)^2$ has maximum value $0.2367$ at $|x| = 0.1940$. Therefore, by regarding local maxima of $J_0(k| \mathbf{x}- \mathbf{x}_s|)^2+J_2(k| \mathbf{x}- \mathbf{x}_s|)^2$, the factor $J_2(k| \mathbf{x}- \mathbf{x}_s|)^2$ does not disturb the identification significantly so almost accurate location $\mathbf{x}_s$ can be identified through the map of $\mathfrak{F}_{{\rm{ITE}}}(\mathbf{x})$. We refer to the Figure 2 for related illustration.

Correspondingly, we can obtain the following result of unique determination in the magnetic permeability contrast case.

Corollary 3.2 (Unique determination in permeability contrast case). With the same condition of the Theorem 3.3, the location of small inhomogeneities can be identified uniquely through the map of $\mathfrak{F}_{{\rm{ITE}}}(\mathbf{x})$.

4.   Numerical experiments

4.1.   Synthetic data experiments

In this section, a set of simulation results are shown to validate the investigated results of Theorems 3.1–3.3. To this end, we applied the background wavenumber of the form $k = 2\pi/\lambda$ with $\lambda = {0.6}\; {\rm{m}}$, ${0.4}\; {\rm{m}}$, and ${0.2}\; {\rm{m}}$, and selected $S = 3$ small circles with locations $\mathbf{x}_1 = ({0.7}\; {\rm{m}}, {0.5}\; {\rm{m}})$, $\mathbf{x}_2 = (-{0.7}\; {\rm{m}}, {0.0}\; {\rm{m}})$, and $\mathbf{x}_3 = ({0.4}\; {\rm{m}}, -{0.6}\; {\rm{m}})$. The background permittivity and permeability are set to $\varepsilon_\mathrm{b} = \mu_\mathrm{b} = 1$, and the size, permittivity, and permeability of the inhomogeneities $\Sigma_s$ are given in Tables 1 and 2. We selected $M = 64$ and $N = 32$, different transmitters and receivers with $| \mathbf{t}_m| = {3}\; {\rm{m}}$ and $| \mathbf{r}_n| = {2}\; {\rm{m}}$ (see Figure 3), and $\Omega$ as a square region $[-{1}\; {\rm{m}}, {1}\; {\rm{m}}]\times[-{1}\; {\rm{m}}, {1}\; {\rm{m}}]$. With this, the scattered field data $u_{{\rm{scat}}}(\mathbf{t}_m, \mathbf{r}_n)$ of $\mathbb{K}$ are generated by solving the Foldy-Lax formulation introduced in ^[32]. After generating the MSR matrix, ${20}{{\rm{dB}}}$ white Gaussian random noise is added to show the robustness.

Example 4.1 (Permittivity contrast case) Figure 4 shows the maps of $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x})$ for Cases 1–4. Based on the results, we determined that every location $\mathbf{x}_s$ can be identified for Case 1. For Case 2, since $\varepsilon_3 > \varepsilon_2 > \varepsilon_1$, the location $\mathbf{x}_3$ can be recognized clearly, but due to the appearance of several artifacts, identifying $\mathbf{x}_1$ is very difficult. Similarly, since $\alpha_1 > \alpha_2 > \alpha_3$ for Case 3, the location $\mathbf{x}_1$ can be identified through the map of $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x})$, but $\mathbf{x}_3$ cannot be recognized. Notice that the value of $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x}_2)$ is greater than the values in the neighborhood of $\mathbf{x}_2$, but it is difficult to distinguish from the artifacts. For Case 4, since

the locations $\mathbf{x}_1$ and $\mathbf{x}_2$ can be recognized. Similar to Case 3, although the value of $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x}_3)$ is greater than the values in the neighborhood of $\mathbf{x}_3$, it is difficult to distinguish from the artifacts.

Example 4.2 (Comparison with different algorithms). Figure 5 shows the imaging results obtained by MUSIC algorithm ^[55], subspace migration technique ^[56], and direct sampling method ^[39] to compare the imaging performances. Based on the simulation results, it is possible to recognize the existence of small inhomogeneities by using MUSIC (map of $\mathfrak{F}_{{\rm{MUSIC}}}(\mathbf{r}))$ and subspace migration (map of $\mathfrak{F}_{{\rm{SUB}}}(\mathbf{r}))$ very clearly. However, it is impossible to distinguish the characteristics (relative proportions of the size and permittivity) of inhomogeneities. It is interesting to observe that the imaging results by using the direct sampling method with multiple sources (map of $\mathfrak{F}_{{\rm{DSM}}}(\mathbf{r}))$ are almost similar to the results via the maps of $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x})$ but map of the $\mathfrak{F}_{{\rm{DSM}}}(\mathbf{r}))$ contains more artifacts than the one of the $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x})$.

Example 4.3 (Permeability contrast case). Figure 6 shows the maps of $\mathfrak{F}_{{\rm{TE}}}(\mathbf{x})$ for Cases 1–4. Similar to the permittivity contrast case, we determined that every location $\mathbf{x}_s$ can be identified for Case 1. For Case 2, since $\mu_3 > \mu_2 > \mu_1$, the rings in the neighborhood of all $\mathbf{x}_s$ can be recognized, but the magnitudes of each ring are different. In contrast, since $\alpha_1 > \alpha_2 > \alpha_3$ for Case 3, the rings with large magnitudes in the neighborhood of $\mathbf{x}_1$ and $\mathbf{x}_2$ can be identified through the map of $\mathfrak{F}_{{\rm{TE}}}(\mathbf{x})$. However, it is difficult to find a ring in the neighborhood of $\mathbf{x}_3$ so that the existence of $\Sigma_3$ cannot be recognized. For Case 4, since

the ring in the neighborhood of $\mathbf{x}_1$ can be recognized, but it is difficult to distinguish the ring in the neighborhood of $\mathbf{x}_2$ and the artifacts in the map. Furthermore, the ring in the neighborhood of $\mathbf{x}_3$ cannot be recognized because the value of $\alpha_3^2/(\mu_3+ \mu_\mathrm{b})$ is close to zero.

Example 4.4 (Comparison with different algorithms). Figure 7 shows the imaging results obtained by MUSIC algorithm, subspace migration technique, and direct sampling method. Similar to the results in Example 4.2, it is possible to recognize the existence of small inhomogeneities through the maps of $\mathfrak{F}_{{\rm{MUSIC}}}(\mathbf{r})$ and $\mathfrak{F}_{{\rm{SUB}}}(\mathbf{r})$ by regarding the rings with large magnitudes but it is still impossible to distinguish size and permeability of inhomogeneities. Map of $\mathfrak{F}_{{\rm{TE}}}(\mathbf{x})$ contains less artifacts than the map of $\mathfrak{F}_{{\rm{DSM}}}(\mathbf{r})$ however, opposite to the results in Example 4.2, imaging results via the maps of $\mathfrak{F}_{{\rm{DSM}}}(\mathbf{r})$ are quite different.

Example 4.5 (Accurate localization in permeability contrast case). Figure 8 shows the maps of $\mathfrak{F}_{{\rm{ITE}}}(\mathbf{x})$ for Cases 1–4. As we already discussed in the Remark 3.6, locations of all inhomogeneities are successfully identified for Cases $1$ and $2$. Moreover, similar to the Example 4.3, it is very difficult to recognize the location $\mathbf{x}_3$ for Cases $3$ and $4$. Notice that opposite to the results in Example 4.1, rings of large magnitudes are included in the neighborhood of all $\mathbf{x}_s$ due to the effect of the factor $J_2(k| \mathbf{x}- \mathbf{x}_s|)^2$ but they do not disturb the recognition of inhomogeneities.

Example 4.6 (Application of multiple frequencies). On the basis of several studies ^{[2,5,28,29,30,38,56,61,65]}, we consider the multi-frequency imaging to improve the quality of the reconstruction. For this purpose, we introduce the following multi-frequency imaging functions

and

where $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x}, k_p)$, $\mathfrak{F}_{{\rm{TE}}}(\mathbf{x}, k_p)$, and $\mathfrak{F}_{{\rm{ITE}}}(\mathbf{x}, k_p)$ are the imaging functions at $k_p$ defined in (3.1), (3.5), and (3.11), respectively. Here, applied wavenumbers $k_p$ are equidistributed in the interval $[k_1, k_P]$ with $k_1 = 2\pi/0.6$ and $k_P = k_{10} = 2\pi/0.2$.

Figure 9 exhibits multi-frequency imaging results. On the basis of the results, we can examine that several artifacts in the single-frequency imaging were successfully eliminated. Nevertheless, it is very difficult to recognize the inhomogeneity $\Sigma_3$ for the Case $3$ via the maps of $\mathfrak{F}_{{\rm{MTM}}}(\mathbf{x})$, $\mathfrak{F}_{{\rm{MTE}}}(\mathbf{x})$, and $\mathfrak{F}_{{\rm{MITE}}}(\mathbf{x})$, and for the Case $4$ via the maps of $\mathfrak{F}_{{\rm{MTE}}}(\mathbf{x})$ and $\mathfrak{F}_{{\rm{MITE}}}(\mathbf{x})$. Hence, development of imaging technique for a further improvement is required.

4.2.   Real data experiments in limited-aperture configuration

Subsequently, we show the simulation results using the Fresnel data introduced in ^[11] at various frequencies. The transmitters and receivers are placed on the circles with radii ${0.72}\; {\rm{m}}$ and ${0.76}\; {\rm{m}}$, respectively, and the ROI $\Omega$ was chosen as the square $[-{0.15}\; {\rm{m}}, {0.15}\; {\rm{m}}]\times[-{0.15}\; {\rm{m}}, {0.15}\; {\rm{m}}]$. Notably, due to the limitation of the experimental setting, the range of receivers is restricted from ${60}°$ to ${300}°$, with a step size of ${5}°$ based on each direction of the transmitters. The transmitters are evenly distributed with step sizes of ${10}°$ from ${0}°$ to ${350}°$ (Figure 10). Thus, the simulation results are obtained with limited-aperture data, even though the theoretical results are derived under the full-aperture measurement hypothesis.

Example 4.7 (Imaging of dielectric objects: transverse magnetic case). Figure 11 shows maps of $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x})$ at various frequencies in the presence of single and multiple dielectric circular objects with the same radii in transverse magnetic (TM) polarization. As we already discussed in Remark 3.1, a blurred imaging result was obtained when $f = {1}{{\rm{GHz}}}$. Moreover, the presence of the two objects at $f = {1}{{\rm{GHz}}}$ cannot be recognized because a peak of large magnitude appeared at their center. Fortunately, every location of objects can be identified by increasing the value of frequency e.g., $f = {4}{{\rm{GHz}}}$ and ${8}{{\rm{GHz}}}$.

Example 4.8 (Imaging of metallic objects: transverse magnetic case). Here, we consider the imaging of a small rectangular-shaped and U-shaped metallic objects in TM polarization. Figure 12 includes maps of $\mathfrak{F}_{{\rm{TM}}}(\mathbf{x})$ with various frequencies. Based on the results, the center of the small rectangular-shaped object can be recognized for any frequency. Interestingly, unlike the other results, a considerable number of artifacts with small magnitudes are included in the map when $f = {8}{{\rm{GHz}}}$. Notice that it is very hard to retrieve the shape of U-shaped object at any frequency. Hence, further improvement is required still.

Example 4.9 (Comparison with different methods: transverse magnetic case). Figures 13 and 14 exhibit the imaging results obtained by MUSIC algorithm and direct sampling method, respectively. Same as the results in Example 4.7 and Figure 11, it is impossible to recognize the existence of multiple objects when applied frequency is low ($f = {1}{{\rm{GHz}}}$). Moreover, the shape of U-shaped object cannot be recognized through the MUSIC and direct sampling method for any frequency. It is interesting to examine the opposite to the synthetic data experiments (see Example 4.2 and Figure 5), it is very difficult to recognize the outline shape of multiple objects through the maps of MUSIC when a high frequency $f = {8}{{\rm{GHz}}}$ was applied.

Additionally, let us consider the imaging results of U-shaped object obtained by the topological derivative (TD) method introduced in [17,Figure 10]. Similar to the results in Figure 12, it is possible to recognize the existence of object at low frequency $f = {2}{{\rm{GHz}}}$. Opposite to the results via the BFM, although some artifacts are also included, one can recognize some features of the shape via the TD at high frequencies $f = 10$ and ${12}{{\rm{GHz}}}$.

Example 4.10 (Imaging of metallic object: transverse electric case). Now, we consider the imaging of a rectangular-shaped metallic object in transverse electric (TE) polarization. Figure 15 shows maps of $\mathfrak{F}_{{\rm{TE}}}(\mathbf{x})$ at various frequencies. Similar to the results in Figure 6, we can observe that two peaks of large magnitudes are included in the map when $f = 2, 4, 6, {8}{{\rm{GHz}}}$ instead of the true location of the object. As discussed in Remark 3.3, an almost exact location of the object is identified through the map when the value of the applied frequency is sufficiently large enough such as $f = {10}{{\rm{GHz}}}$ and ${12}{{\rm{GHz}}}$.

Example 4.11 (Comparison with different methods: transverse electric case). Figures 16 and 17 show the imaging results for a rectangular-shaped metallic object in transverse electric case through the MUSIC algorithm and direct sampling method, respectively. Similar to the results in Example 4.10, two peaks of large magnitudes were included in the map of $\mathfrak{F}_{{\rm{MUSIC}}}(\mathbf{x})$ and $\mathfrak{F}_{{\rm{DSM}}}(\mathbf{x})$ when $f = 2$ and ${4}{{\rm{GHz}}}$. Moreover, an almost exact shape of the object is retrieved via the direct sampling method when a high frequency $f = {8}{{\rm{GHz}}}$ was applied. Opposite to the TM case in Example 4.9 and Figure 13, it is very difficult to recognize the shape of the object via MUSIC algorithm when high frequencies $f = 8$, $10$, and ${12}{{\rm{GHz}}}$ were applied.

Example 4.12 (Accurate imaging of metallic object: transverse electric case). Here, we apply $\mathfrak{F}_{{\rm{ITE}}}(\mathbf{x})$ for identifying a rectangular-shaped metallic object in transverse electric (TE) polarization. Similar to the results in Example 4.5, the location of the object was successfully identified for various frequencies, refer to Figure 18. However, at the $f = {2}{{\rm{GHz}}}$, it seems to be difficult to recognize the object due to the blurring effect in the neighborhood of the object. Similarly with the results in Example 4.10, some artifacts are also included in the map but there is no difficulty to recognize the location of object.

Example 4.13 (Multi-frequency imaging: both transverse magnetic and electric cases). For the final example, we consider the multi-frequency imaging to improve the quality of the reconstruction. For this purpose, we adopted multi-frequency imaging functions $\mathfrak{F}_{{\rm{MTM}}}(\mathbf{x})$, $\mathfrak{F}_{{\rm{MTE}}}(\mathbf{x})$, and $\mathfrak{F}_{{\rm{MITE}}}(\mathbf{x})$ defined in Example 4.6. We applied frequency ranges from $1$ to ${8}{{\rm{GHz}}}$ in steps of ${1}{{\rm{GHz}}}$ for imaging single and multiple circular objects in TM polarization, from $2$ to ${16}{{\rm{GHz}}}$ in steps of ${2}{{\rm{GHz}}}$ for imaging rectangular-shaped object in both TM and TE polarizations and U-shaped metallic objects in TM polarization.

Figure 19 shows multi-frequency imaging results. Throughout the results, we can observe that most of objects were successfully retrieved in both TM and TE polarizations. Moreover, by regarding the large magnitudes in the map, an outline shape of U-shaped metallic object can be recognized. Hence, same as the previous studies, application of multiple frequencies on the BFM guarantees better imaging results.

We compare the imaging results via multi-frequency topological derivative (MFTD) method in TM polarization case, refer to Figures 11–14 in ^[17]. By comparing the results, we can observe that obtained imaging results are very close to the ones via MFTD hence, multi-frequency BFM will guarantee successful reconstruction of shape, size and number of objects.

5.   Conclusions

We considered the BFM for identifying small 2D inhomogeneities when their permittivities or permeabilities differ from the homogeneous background. Through careful analysis, we show that the imaging function can be expressed by the characteristics (size, permittivity, and permeability) of the inhomogeneities and the Bessel function of orders zero and one for the permittivity and permeability contrast cases, respectively. Based on this, we observed and verified the applicability and several properties of the BFM. Numerical simulations with synthetic and experimental data were provided to support the theoretical result.

To achieve the best imaging of inhomogeneities in the permeability contrast case, an improved version of the BFM is needed. Improving the BFM to identify the exact locations of inhomogeneities for permeability contrast case will be the focus of future studies. Based on the experimental data experiments, we determined that the BFM is highly effective in the limited-aperture inverse scattering problem. The theoretical study of BFM in limited-aperture inverse scattering problem will be a significant research topic. Moreover, imaging results from experimental data indicate that the BFM is still effective in situations not reduced to the case of penetrable objects. Hence, extension to the imaging of metallic objects will be also an interesting research. Furthermore, application to the BFM with monostatic ^[40] and bistatic ^[41] measurement configuration will be the forthcoming work.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors would like to acknowledge the anonymous reviewers for their comments that help to increase the quality of the paper. This work was supported by the research program of the Kookmin University.

Conflict of interest

The authors declare no conflicts of interest regarding the publication of this paper.