On the identification of small anomaly in microwave imaging without homogeneous background information

1.   Introduction

Generally, one of the purposes of the microwave imaging is to localize or detect unknown objects from measured electromagnetic waves in the high-frequency regime (between ${300}\; {\rm{MHz}}$ and ${300}\; {\rm{GHz}}$). Due to this reason, identification of the location and shape of unknown anomaly whose values of dielectric permittivity and electric conductivity differ from the homogeneous background medium from scattering parameter data is an important and interesting research topic in microwave imaging. Many authors have proposed various remarkable algorithms for retrieving information on anomaly for example, Born iterative method for reconstructing permittivity distribution ^[13] and brain stroke detection ^[20], level-set method for shape reconstruction of unknown objects ^[16] and breast cancer detection ^[21], distorted iterated virtual experiments scheme for imaging unknown scatterers ^[28], conjugate gradient method for breast imaging ^[11], Levenberg-Marquardt technique for recovering parameter distribution ^[17]. We also refer to remarkable mathematical and experimental studies ^{[1,3,12,15,26,45,46,47,48]}. However, the success of iteration-based algorithms significantly depends on the priori information and initial guess, which must be close to the unknown anomaly, refer to ^[25,39].

For this reason, various non-iterative schemes in inverse scattering problems without a priori information about unknown anomaly have also been developed to retrieve the location and shape, such as a variational algorithm based on the inverse Fourier transform to retrieve small electromagnetic inhomogeneities ^[5,6], direct sampling method for localizing small electromagnetic inhomogeneities ^[22,23] and anomaly detection in real-world experiment ^[43], factorization method for crack detection ^[10], shape reconstruction of unknown obstacles ^[18] and numerical study for anomaly imaging ^[33], linear sampling method for imaging of unknown scatters in limited-aperture problem ^[8] and crack-like defects ^[24], MUltiple SIgnal Classification (MUSIC) algorithm for identifying small anomalies ^[38], fast imaging of small targets in limited-aperture measurement configuration ^[35] and real-world application of anomaly detection ^[34], orthogonality sampling method for imaging unknown targets ^[40], qualitative microwave imaging ^[2], anomaly detection in microwave imaging ^[37] and topological derivative strategy for imaging crack-like defects ^[31] and retrieving unknown scatterers in 3D ^[27].

The subspace migration (SM) algorithm is a well-known, non-iterative imaging technique in both inverse scattering problem and microwave imaging. It has been applied successfully to the various problems such as localization of small targets ^[4], identification of extended objects ^[9], fast imaging of curve-like cracks ^[29,30] and anomaly detection in microwave imaging ^[32,36]. Throughout several studies, it has been confirmed that SM is a fast, robust and effective technique for retrieving unknown anomaly from scattering parameter data. However, for successful application in microwave imaging, accurate values of background permittivity and conductivity must be known because the exact background wavenumber value must be applied. Generally, most researchers have used the statistical values of the background material instead of the true ones and obtained inaccurate locations and shapes of anomaly. This can be examined through the results of numerical simulations but no reliable mathematical theory explaining this phenomenon has yet been developed.

In this paper, we apply the SM to retrieve unknown anomaly from scattering parameter data without complete information about the background material; that is, inaccurate values of background permittivity and conductivity are applied. To explain the appearance of an inaccurate location and shape of an anomaly, we show that the imaging function of the SM can be written as the infinite series of Bessel function of integer order, antenna arrangement and applied inaccurate values of background permittivity and conductivity. This enables us to theoretically explain the appearance of an inaccurate location and shape of an anomaly through the SM. To confirm the theoretical results, simulation results with inaccurate values of background permittivity and conductivity are also presented.

The rest of this paper is organized as follows. In Section 2, the two-dimensional direct problem and imaging function of the SM are introduced. The structure of the imaging function with inaccurate values of background permittivity and conductivity is revealed in Section 3. In Section 4, a set of numerical simulation results with synthetic data generated by CST STUDIO SUITE is presented. Finally, a short conclusion is provided in Section 5.

2.   Scattering parameter and imaging function of the subspace migration

Let $D$ be a circle-like anomaly with radius $\alpha$ and center $\mathbf{r}_\star$ that is surrounded by a circular array of antennas $\mathbf{A}_n$ with location $\mathbf{a}_n$ and $| \mathbf{a}_n| = R$, $n = 1, 2, \cdots, N(>2)$. Throughout this paper, we use $\Omega$ to denote a homogeneous domain filled by matching liquid such that $D\subset\Omega$ and assume that all materials $D$ and $\Omega$ are nonmagnetic and classified by their value of dielectric permittivity and electric conductivity at a given angular frequency of operation $\omega$, i.e., the value of magnetic permeability is constant for every $\mathbf{r}\in\Omega$ say, $\mu(\mathbf{r})\equiv \mu_\mathrm{b} = {1.257 \times {\rm{1}}{{\rm{0}}^{{\rm{ - 6}}}}}\; {\rm{H/m}}$, refer to ^[41]. We use $\varepsilon_\mathrm{b}$ and $\sigma_\mathrm{b}$ to denote the background permittivity and conductivity, respectively. Analogously, let $\varepsilon_\star$ and $\sigma_\star$ be those of $D$. Then, we introduce the piecewise constant permittivity $\varepsilon(\mathbf{r})$ and conductivity $\sigma(\mathbf{r})$,

respectively. With this, let $k_{{\rm{b}}}$ be the background wavenumber that satisfies

and further assume that

With this, the time-dependent, homogeneous, linear Maxwell Equations take the form:

where $\mathcal{E}(\mathbf{r}, t)\in\mathbb{R}^3$ and $\mathcal{H}(\mathbf{r}, t)\in\mathbb{R}^3$ are the electric and magnetic fields, respectively. Here, we consider time-harmonic solutions to the (2.2) such that

Then, $\mathbf{E}(\mathbf{r})\in\mathbb{C}^3$ and $\mathbf{H}(\mathbf{r})\in\mathbb{C}^3$ satisfy

For a detailed description, we refer to ^[7].

Let $\mathbf{E}_{{\rm{inc}}}(k_{{\rm{b}}}, \mathbf{r}, \mathbf{a}_m)$ be the incident electric field due to the point current density $\mathbf{J}$ at $\mathbf{A}_m$. Then, based on (2.3), it satisfies

Analogously, let $\mathbf{E}_{{\rm{tot}}}(\mathbf{a}_n, \mathbf{r})$ be the total field measured at $\mathbf{A}_n$ in the presence of $D$ that satisfies

with transmission condition on $\partial D$.

Let $S(n, m)$ be the $S$-parameter (or scattering parameter), which is defined as the ratio of the output voltage (or reflected waves) at the $\mathbf{A}_n$ antenna and the input voltage (or incident waves) at the $\mathbf{A}_m$ (see ^[41] for instance). We use $S_{{\rm{inc}}}(n, m)$ and $S_{{\rm{tot}}}(n, m)$ to denote the incident and total field $S$-parameters, respectively, due to the absence and presence of $D$. Correspondingly, we let $S_{{\rm{scat}}}(n, m) = S_{{\rm{tot}}}(n, m)-S_{{\rm{inc}}}(n, m)$ be the scattered field $S$-parameter. Then, based on ^[19], $S_{{\rm{scat}}}(n, m)$ can be represented as follows:

where $k_0$ denotes the lossless background wave number that satisfies $k_0^2 = \omega^2 \mu_\mathrm{b} \varepsilon_\mathrm{b}$.

In this paper, we adopt the simulation configuration introduced in ^[32,43,44]. Notice that the height of microwave machine can be said to be long enough, only the $z$-component of the incident and total fields can be handled based on the mathematical treatment of the scattering of time-harmonic electromagnetic waves from thin infinitely long cylindrical obstacles. Correspondingly, by denoting $\mathrm{E}_{{\rm{inc}}}^{(z)}$ and $\mathrm{E}_{{\rm{tot}}}^{(z)}$ as the $z$-components of the incident and total fields, respectively, $S_{{\rm{scat}}}(n, m)$ can be written as follows:

Unfortunately, exact expression of the field $\mathrm{E}_{{\rm{tot}}}^{(z)}(k_{{\rm{b}}}, \mathbf{a}_n, \mathbf{r}')$ is unknown thus we cannot design the imaging function by using $S_{{\rm{scat}}}(n, m)$ of (2.4) directly. Since the condition (2.1) holds, we can apply the Born approximation such that (see ^[42] for instance)

where $H_0^{(1)}$ denotes the Hankel function of order zero of the first kind. Correspondingly, $S_{{\rm{scat}}}(n, m)$ of (2.4) can be written as

Now, we introduce the imaging function. To this end, let us generate the scattering matrix $\mathbb{K}$ such that

See ^[32] for an explanation of why the diagonal elements of $\mathbb{K}$ are set to zero. Since there is a single, small anomaly, the singular value decomposition (SVD) of $\mathbb{K}$ can be written as

where $*$ denotes the mark of Hermitian, $\tau_n$ are the singular values and $\mathbf{U}_n$ and $\mathbf{V}_n$ are the left and right singular vectors of $\mathbb{K}$, respectively. Then, on the basis of (2.5) and (2.7), we define a unit vector: For each $\mathbf{r}\in\Omega$,

With this, we introduce the following imaging function of the SM: For each $\mathbf{r}\in\Omega$,

where $\langle \mathbf{U}, \mathbf{V}\rangle = \mathbf{U}^* \mathbf{V}$ and $\overline{ \mathbf{V}}_1$ denotes the complex conjugate of $\mathbf{V}_1$. Then, based on ^[4]

and the orthonormal property of singular vectors, the value of $\mathfrak{F}(k_{{\rm{b}}}, \mathbf{r})$ will be close to $1$ when $\mathbf{r}\in D$ and less than $1$ at $\mathbf{r}\in\Omega\backslash D$, so the location and outline shape of $D$ can be identified through the map of $\mathfrak{F}(k_{{\rm{b}}}, \mathbf{r})$.

Let us emphasize that, to generate the test vector $\mathbf{W}(k_{{\rm{b}}}, \mathbf{r})$ of (2.8), the exact value of $k_{{\rm{b}}}$ must be known, i.e., a priori information of the $\varepsilon_\mathrm{b}$ and $\sigma_\mathrm{b}$ must be available. However, because these values are statistical, the exact value may not be unknown. For this reason, we assume that the exact values of $\varepsilon_\mathrm{b}$ and $\sigma_\mathrm{b}$ are unknown and apply an alternative value $k_{{\rm{a}}}$ instead of the true $k_{{\rm{b}}}$ and set a unit test vector $\mathbf{W}(k_{{\rm{a}}}, \mathbf{r})$ from (2.8). Then, by using the imaging function $\mathfrak{F}(k_{{\rm{a}}}, \mathbf{r})$ from (2.9), the exact location and shape of $D$ cannot be retrieved. Fortunately, we can recognize the existence of $D$ and the identified location is shifted in a specific direction.

3.   Analysis of the imaging function

In this section, we explore the structure of the imaging function $\mathfrak{F}(k_{{\rm{a}}}, \mathbf{r})$ to explain that retrieved location of $D$ is shifted in a specific direction and size is smaller or larger than the true anomaly. To explain this phenomenon, we explore the structure of the imaging function.

Theorem 3.1. Let $\boldsymbol{\theta}_n = \mathbf{a}_n/| \mathbf{a}_n| = (\cos\theta_n, \sin\theta_n)$ and $k_{{\rm{b}}} \mathbf{r}'-k_{{\rm{a}}} \mathbf{r} = |k_{{\rm{b}}} \mathbf{r}'-k_{{\rm{a}}} \mathbf{r}|(\cos\phi, \sin\phi)$. If $\mathbf{a}_n$ satisfies $| \mathbf{a}_n- \mathbf{r}|\gg\left\{{1/4|k_{{\rm{a}}}|, 1/4|k_{{\rm{b}}}|}\right\}$ for all $n = 1, 2, \cdots, N$, then $\mathfrak{F}(\mathbf{r}, k_{{\rm{a}}})$ can be represented as follows:

where $J_s$ denotes the Bessel function of order $s$ and

Proof. Since $\mathbb{K}\approx\tau_1 \mathbf{U}_1 \mathbf{V}_1^*$, we can examine that

Based on the assumption $| \mathbf{a}_n- \mathbf{r}|\gg\left\{{1/4|k_{{\rm{a}}}|, 1/4|k_{{\rm{b}}}|}\right\}$ for all $n = 1, 2, \cdots, N$, the following asymptotic forms of the Hankel function hold (see [14,Theorem 2.5], for instance)

Then, $\mathbf{W}(k_{{\rm{a}}}, \mathbf{r})$ and $\mathbb{K}$ can be represented as

and

respectively. Here the constant $C$ is given by

Note that, since the following JacobiAnger expansion formula holds uniformly

we have for $n = 1, 2, \cdots, N$

and correspondingly,

With this, we can evaluate

Hence,

Since $\langle \mathbf{W}(\mathbf{r}), \mathbf{U}_1\rangle\langle \mathbf{W}(\mathbf{r}), \overline{ \mathbf{V}}_1\rangle = 1$, $J_0(|k_{{\rm{b}}} \mathbf{r}'-k_{{\rm{a}}} \mathbf{r}|) = J_0(2|k_{{\rm{b}}} \mathbf{r}'-k_{{\rm{a}}} \mathbf{r}|) = 1$ and $\Psi(k_{{\rm{b}}}, k_{{\rm{a}}}, \mathbf{r}) = \Psi(2k_{{\rm{b}}}, 2k_{{\rm{a}}}, \mathbf{r}) = 0$ when $k_{{\rm{b}}} \mathbf{r}' = k_{{\rm{a}}} \mathbf{r}$, we have

Therefore,

With this, we can obtain the structure (3.1). □

From the derived structure (3.1), we can observe that since $J_0(|k_{{\rm{b}}} \mathbf{r}'-k_{{\rm{a}}} \mathbf{r}|) = 1$ and $\Psi(k_{{\rm{b}}}, k_{{\rm{a}}}, \mathbf{r}) = 0$ when $\mathbf{r} = (k_{{\rm{b}}}/k_{{\rm{a}}}) \mathbf{r}'$ for $\mathbf{r}'\in D$, an inaccurate location and shape of $D$ must be retrieved through the map of $\mathfrak{F}(k_{{\rm{a}}}, \mathbf{r})$. This is the theoretical reason why an inaccurate location and shape of the anomaly is retrieved when inaccurate values of $\varepsilon_\mathrm{b}$ and $\sigma_\mathrm{b}$ were applied. Further properties will be discussed in the simulation results.

4.   Simulation results

Here, we present simulation results to support the result in Theorem 3.1. To this end, a circular array of $N = 16$ antennas $\mathbf{A}_n$ is used to transmit and receive signals operated at $f = {1.24}\; {\rm{GHz}}$. The location of the antennas was set to

and the search domain $\Omega$ was selected as a square region $(-{0.1}{\; {\rm{m}}}, {0.1}{\; {\rm{m}}})\times(-{0.1}{\; {\rm{m}}}, {0.1}{\; {\rm{m}}})$ with $(\varepsilon_\mathrm{b}, \sigma_\mathrm{b}) = (20 \varepsilon_0, {0.2}{{\rm{S}}/{\rm{m}}})$. Here, $\varepsilon_0 = {8.854 \times {\rm{1}}{{\rm{0}}^{{\rm{ - 12}}}}}{{\rm{F}}/{\rm{m}}}$ is the vacuum permittivity. Correspondingly, the exact value of the background wavenumber is $k_{{\rm{b}}} = 116.5273+ 8.4020i$. For anomalies, we selected two small balls $D_1$ and $D_2$ with centers $\mathbf{r}_1 = ({0.01}{\; {\rm{m}}}, {0.03}{\; {\rm{m}}})$ and $\mathbf{r}_2 = (-{0.04}{\; {\rm{m}}}, -{0.02}{\; {\rm{m}}})$, same radii $\alpha = {0.01}{\; {\rm{m}}}$ and material properties $(\varepsilon_1, \sigma_1) = (55 \varepsilon_0, {1.2}{{\rm{S}}/{\rm{m}}})$ and $(\varepsilon_2, \sigma_2) = (45 \varepsilon_0, {1.0}{{\rm{S}}/{\rm{m}}})$. With these settings, the measurement data $S_{{\rm{scat}}}(n, m)$ of (2.6) and the incident field data of (2.8) were generated by CST STUDIO SUITE.

Example 4.1. (Only exact value of $\varepsilon_\mathrm{b}$ is unknown) First, we consider the case where only the exact value of $\varepsilon_\mathrm{b}$ is unknown. Instead, of the application of $\varepsilon_\mathrm{b}$, we applied alternative values $\varepsilon_\mathrm{a}$ and corresponding wavenumber

Note that we already assumed that $\omega \varepsilon_\mathrm{b}\gg \sigma_\mathrm{b}$. Thus, if the condition $\omega \varepsilon_\mathrm{a}\gg \sigma_\mathrm{b}$ is satisfied, the identified location becomes

Hence, identified anomalies will be concentrated at the origin and their retrieved sizes will be smaller than the true one when $\varepsilon_\mathrm{a} > \varepsilon_\mathrm{b}$. Otherwise, identified anomalies will be far from the origin and their retrieved sizes will be larger than the true one when $\varepsilon_\mathrm{a} < \varepsilon_\mathrm{b}$. See Figure 1 for a related illustration.

Figure 2 shows maps of $\mathfrak{F}(k_{{\rm{a}}}, \mathbf{r})$ with various selections of $\varepsilon_\mathrm{a}$ in the presence of $D_1$. As we discussed above, the location of the retrieved anomaly gets closer to the origin and it becomes smaller as the value of $\varepsilon_\mathrm{a}$ increases $($here, $\varepsilon_\mathrm{a} = 3 \varepsilon_\mathrm{b}, 10 \varepsilon_\mathrm{b})$. Otherwise, as the value of $\varepsilon_\mathrm{a}$ decreases $($here, $\varepsilon_\mathrm{a} = 0.5 \varepsilon_\mathrm{b})$, the identified location becomes far from the origin and the size becomes larger. If $\varepsilon_\mathrm{a} = 0.01 \varepsilon_\mathrm{b}$ i.e., the value of $\varepsilon_\mathrm{a}$ is very small, it is difficult to distinguish the $D_1$ and artifacts.

Notice that since

due to the oscillating property of the Bessel function, several artifacts will be included in the map of $\mathfrak{F}(k_{{\rm{a}}}, \mathbf{r})$ if $|k_{{\rm{a}}}|$ is large, i.e., the value of $\varepsilon_\mathrm{a}$ is large enough compared to $\varepsilon_\mathrm{b}$. In contrast, the map of $\mathfrak{F}(k_{{\rm{a}}}, \mathbf{r})$ will contain no artifacts but the imaging result will be blurred if the value of $\varepsilon_\mathrm{a}$ is small enough compared to $\varepsilon_\mathrm{b}$. This is the reason why several artifacts are included in the map of $\mathfrak{F}(k_{{\rm{a}}}, \mathbf{r})$ when $\varepsilon_\mathrm{a} = 10 \varepsilon_\mathrm{b}$ and why the obtained image is blurred when $\varepsilon_\mathrm{a} = 0.1 \varepsilon_\mathrm{b}$ and $\varepsilon_\mathrm{a} = 0.01 \varepsilon_\mathrm{b}$. We can observe the same phenomenon in the presence of multiple anomalies $D_1$ and $D_2$, as shown in the Figure 3.

Example 4.2. (Only exact value of $\sigma_\mathrm{b}$ is unknown) Next, we consider the case where only the exact value of $\sigma_\mathrm{b}$ is unknown and apply an alternative one $\sigma_\mathrm{a}$ such that

Same as the Example 4.1., if $\sigma_\mathrm{a}$ satisfies $\omega \varepsilon_\mathrm{a}\gg \sigma_\mathrm{a}$, the identified location becomes

Hence, it will be possible to retrieve almost accurate shapes and locations of anomalies when $\sigma_\mathrm{a}$ is sufficiently small.

Figure 4 shows maps of $\mathfrak{F}(k_{{\rm{a}}}, \mathbf{r})$ with various selections of $\sigma_\mathrm{a}$ in the presence of $D_1$. In contrast to the results in Example 4.1., almost the exact location and shape of $D_1$ were retrieved if $\sigma_\mathrm{a} < 3 \sigma_\mathrm{b}$, i.e., when $\sigma_\mathrm{a}$ was sufficiently small. The location of the retrieved anomaly gets closer to the origin and its size becomes smaller as the value of $\varepsilon_\mathrm{a}$ increases $($here, $\varepsilon_\mathrm{a} = 3 \varepsilon_\mathrm{b}, 10 \varepsilon_\mathrm{b})$. Unfortunately, it is very difficult to recognize $D_1$ due to the appearance of a huge artifact with a large magnitude if $\sigma_\mathrm{a}$ is not small $\sigma_\mathrm{a} = 10 \sigma_\mathrm{b} = {2}{{\rm{S}}/{\rm{m}}}$.

We can observe the same phenomenon in the presence of multiple anomalies $D_1$ and $D_2$, as shown in the Figure 5, and conclude that it will be possible to retrieve the accurate shape and location of anomalies by choosing a very small $($close to zero$)$ value of $\sigma_\mathrm{a}$ when the exact value of background permittivity is known.

Example 4.3. (Identification of circle and rectangular shaped anomalies) Here, we consider the imaging of anomalies with different shapes. To this end, we applied $f = {1.0}\; {\rm{GHz}}$, used $N = 36$ antennas $\mathbf{A}_n$, and selected $D_1$ as a ball of Examples 1 and 2 except $(\varepsilon_1, \sigma_1) = (45 \varepsilon_0, {1.0}{{\rm{S}}/{\rm{m}}})$, and $D_2$ as a square with vertices $(-0.05, -0.03)$, $(-0.03, -0.03)$, $(-0.03, -0.01)$ and $(-0.05, -0.01)$ with $(\varepsilon_2, \sigma_2) = (45 \varepsilon_0, {1.0}{{\rm{S}}/{\rm{m}}})$. With this configuration, the scattering parameter data were generated by using the FEKO (Feldberechnung für Körper mit beliebiger Oberfläche).

Figure 6 shows maps of $\mathfrak{F}(k_{{\rm{a}}}, \mathbf{r})$ when $\varepsilon_\mathrm{a}\ne \varepsilon_\mathrm{b}$ and $\sigma_\mathrm{a} = \sigma_\mathrm{b}$. Similar to the results in Example 4.1., identified anomalies are concentrated at the origin and their retrieved sizes are smaller than the true one when $\varepsilon_\mathrm{a} > \varepsilon_\mathrm{b}$. Moreover, identified anomalies located far from the origin and their retrieved sizes are larger than the true one when $\varepsilon_\mathrm{a} < \varepsilon_\mathrm{b}$. However, opposite to the results in Examples 1 and 2, it is hard to recognize the shape of $D_1$ and $D_2$ due to the appearance of several artifacts in the neighborhood of anomalies.

Figure 7 shows maps of $\mathfrak{F}(k_{{\rm{a}}}, \mathbf{r})$ when $\varepsilon_\mathrm{a} = \varepsilon_\mathrm{b}$ and $\sigma_\mathrm{a}\ne \sigma_\mathrm{b}$. Similar to the results in Example 4.2., we can examine that it is possible to recognize the outline shape of anomalies by choosing a very small value of $\sigma_\mathrm{a}$. However, exact shape of anomalies cannot be retrieved still.

5.   Conclusions

The structure of the imaging function of SM for retrieving small anomalies from scattering matrix is revealed when complete information of the background medium is not available. On the basis of its relationship with the infinite series of Bessel function of the first kind, we have theoretically confirmed why the accurate shape and location of anomalies cannot be retrieved.

The main subject of this paper is the imaging of small anomaly in two-dimensional microwave imaging. An extension to multiple, small anomalies will be carried out in forthcoming work. Moreover, the development of an effective algorithm for retrieving the exact value of background wavenumber will be an interesting research subject. Finally, we expect that the methodology presented in this paper could be applied to real-world microwave imaging with inhomogeneous background.

Use of AI tools declaration

The author declares he has not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The author would like to acknowledge anonymous reviewers for their comments that help to increase the quality of the paper. The author is also grateful to Sangwoo Kang, Kwang-Jae Lee and Seong-Ho Son for helping in generating scattering parameter data. This work was supported by the research program of the Kookmin University.

Conflict of interest

The author declares no conflicts of interest regarding the publication of this paper.