Real-time tracking of moving objects from scattering matrix in real-world microwave imaging

1.   Introduction

The real-time tracking of unknown objects using microwaves is an important and interesting inverse scattering problem that arises in fields such as physics, engineering, and military services, and is nowadays highly related to human life ^{[5,8,15,33,35,57]}. Most studies have focused on the development of detection algorithms for unknown objects located in a wide area; however, the detection of the movements of small objects or changes in permittivity/conductivity distributions in small or narrow regions has not been sufficiently investigated. Detecting small objects is not an easy problem to solve, but can be applied to various real-world problems, such as the diagnoses of cerebral hemorrhages ^[32], imaging of crack propagations in walls or bridges ^[16], through-wall imaging ^[6], breast cancer detection ^[13], and landmine detection ^[10]. In general, objects exhibit material properties (permittivities and conductivities) that are different from their background media; thus, most studies have focused on retrieving complete information about these properties. Unfortunately, owing to the intrinsic ill-posedness and nonlinearity of the inverse scattering problem, it cannot be successfully and satisfactorily resolved.

To solve this problem, various inversion techniques and corresponding computational environments have been investigated. The most popular and appropriate approaches in real-world applications are based on the Newton-type iteration scheme, which involves retrieving information about the shape, location, and material properties of unknown objects (minimizers); these unknown objects minimize the discrete norm (generally, $\ell^2-$norm) between the measured data in the presence of true and man-made objects. For examples, we refer to the Newton method for microwave tomography ^[56], Born approximation for compressive sensing imaging ^[34], level set technique for crack-like thin electromagnetic object ^[47], Gauss-Newton method for breast cancer monitoring ^[50], Levenberg-Marquardt method for retrieving permittivity distribution in microwave imaging ^[17], and Newton–Kantorovich technique for arm imaging ^[31].

Although iteration-based techniques have shown their feasibilities, some preceding conditions, such as good initial guesses that are close to the objects, a priori information about unknown objects, appropriate regularization terms significantly dependent on the problem, evaluations of complex Fréchet (or domain) derivatives must be fulfilled to guarantee successful iteration procedures. Furthermore, large computational costs are still incurred, and extensions to multiple objects still prove to be difficult. Hence, Newton-type iteration schemes are not appropriate in designing real-time detection of moving objects.

As alternatives, various non-iterative algorithms have been investigated, e.g., the bifocusing method ^{[24,26,27,53]}, direct sampling method ^{[22,23,39,52]}, MUltiple SIgnal Classification (MUSIC) algorithm ^[4,36,41,46], linear sampling method ^[1,11,19,29], topological derivatives ^[7,18,30,38], Kirchhoff and subspace migrations ^[2,3,37,45], and orthogonality sampling method ^{[20,43,44,48]}. It is worth mentioning that, although complete information about unknown objects cannot be retrieved via non-iterative algorithms, such algorithms are fast, effective, and stable in detecting multiple objects without any a priori information. So, instead of completely reconstructing the objects, it would be appropriate to design a real-time algorithm for detecting moving objects. Let us emphasize that, in order to apply MUSIC or subspace migration for identifying moving objects, a careful threshold of singular values of the scattering matrix is essential. Topological derivative based techniques or linear sampling methods are very stable and effective non-iterative techniques, but additional operations such as solving adjoint problems or nonlinear integral equations are required. Hence, these methods would be inappropriate as real-time detection techniques. Fortunately, the bifocusing method, direct and orthogonality sampling methods with multiple sources, and Kirchhoff migration require no additional operations; therefore, we believe that they can be applicable to design a real-time detection algorithm.

Let us emphasize that the imaging performance of direct and orthogonality sampling methods is significantly dependent to the location of source (see ^[44,52] for instance). In this paper, regarding this issue, we apply Kirchhoff migration for a real-time tracking of moving small objects from the scattering matrix, whose elements are measured scattering parameters at a fixed frequency. In contrast to the traditional studies in inverse scattering problems, the diagonal elements of a scattering matrix cannot be determined from a microwave machine, i.e., scattering parameters cannot be obtained when the transmitting and receiving antennas are at the same location. Considering such a limitation, an imaging function based on the Kirchhoff migration for tracking moving objects is designed, and its mathematical structure was rigorously analyzed by establishing a relation with an infinite series of Bessel functions of integer order and antenna configurations. Based on this relation, we can demonstrate that the designed tracking algorithm is fast and effective. Moreover, it is possible to guarantee unique determination of moving objects. To illustrate the feasibilities of the designed technique and to avoid committing inverse crimes, numerical simulations are performed with experimental data generated by a microwave machine ^[28].

The remainder of this paper is organized as follows. Section 2 briefly introduces the forward problem and scattering parameters caused by the presence of objects. In Section 3, we describe the Kirchhoff migration-based imaging algorithm without the diagonal elements of the scattering matrix, analyze the mathematical structure of the imaging function by establishing an infinite series of Bessel functions of integer order and antenna setting, and discuss some properties of the imaging function. Section 4 discusses the results of the simulation, which used real-world data to demonstrate the effectiveness of the algorithm. Section 5 concludes the paper.

2.   The forward problem and scattering parameters

In this section, the forward model and scattering parameters are briefly introduced in the presence of a set of objects with small diameter. Let $D_m(t)$, $m = 1, 2, \cdots, M$, be a (cross-sectional) object with location ${\bf r}_m(t)$ at time $t$, where $D(t)$ denotes the collection of $D_m(t)$. Throughout this study, all $D_m(t)$ are included in a homogeneous region of interest (ROI) $\Omega\subset\mathbb{R}^2$ and surrounded by several transmitting and receiving antennas ${\bf A}_n$ located at ${\bf a}_n$, $n = 1, 2, \cdots, N$, with $| {\bf a}_n| = R$. We denote $\mathcal{A}$ as the collection of antennas ${\bf A}_n$ and assume that all $D_m(t)$ and $\Omega$ are characterized by their dielectric permittivity and electric conductivity at a given angular frequency $\omega = 2\pi f$, i.e., the magnetic permeability of all objects are constant such that $\mu({\bf r}, t)\equiv \mu_\mathrm{b} = 4\pi\times10^{-7}{{\rm{H}}/{\rm{m}}}$, ${\bf r}\in\Omega$. We denote $\varepsilon({\bf r}, t)$ and $\sigma({\bf r}, t)$ as the piecewise constant permittivity and conductivity

respectively. With this, we denote $k_\mathrm{b}$ as the background wavenumber that satisfies

and assume that the following relationships hold for $m = 1, 2, \cdots, M$:

where $\text{diam}(D_m)$ denotes the diameter of $D_m$.

Let ${\bf E}_{{\rm{inc}}}({\bf r}, {\bf a}_q)$ be the time-harmonic incident electric field in a homogeneous medium resulting from a point current density at ${\bf A}_q$. Then, based on the Maxwell equation, ${\bf E}_{{\rm{inc}}}({\bf r}, {\bf a}_q)$ satisfies

Let ${\bf E}_{{\rm{tot}}}({\bf r}, {\bf a}_p)$ be the corresponding total electric field due to the existence of $D(t)$ measured at ${\bf A}_p$. Then, ${\bf E}_{{\rm{tot}}}({\bf r}, {\bf a}_p)$ satisfies

with the transmission condition on $\partial D_m$, $m = 1, 2, \cdots, M$. Here, the time-harmonic dependence $e^{-i\omega t}$ is assumed.

The scattering parameter ($S-$parameter) $S(p, q, t)$, $p, q = 1, 2, \cdots, N$, is defined as the ratio of the output voltage at ${\bf A}_p$ and the input voltage at ${\bf A}_q$ at time $t$. Let $S_{{\rm{tot}}}(p, q, t)$ and $S_{{\rm{inc}}}(p, q, t)$ be the scattering parameter data in the presence and absence of $D(t)$, respectively, and denote $S_{{\rm{scat}}}(p, q, t) = S_{{\rm{tot}}}(p, q, t)-S_{{\rm{inc}}}(p, q, t)$ as the measurement data^*. Based on the simulation setup introduced in recent works ^[28,42,45], only the $z$-components of ${\bf E}_{{\rm{inc}}}({\bf r}, {\bf a}_q)$ and ${\bf E}_{{\rm{tot}}}({\bf r}, {\bf a}_p)$ can be handled so that, by denoting $\mathrm{E}_{{\rm{inc}}}^{(z)}({\bf r}, {\bf a}_q)$ and $\mathrm{E}_{{\rm{tot}}}^{(z)}({\bf r}, {\bf a}_p)$ as the $z$-component of incident and total fields, respectively, $S_{{\rm{scat}}}(p, q, t)$ can be presented as the following integral equation formula:

^*This subtraction is essential and effective to remove unknown modeling error.

where $k_0$ denotes the lossless background wavenumber that satisfies $k_0^2 = \omega^2 \varepsilon_\mathrm{b} \mu_\mathrm{b}$, refer to ^[21].

3.   Kirchhoff migration for a real-time detection: introduction, analysis, and some properties

3.1.   Introduction to imaging function of the Kirchhoff migration

Herein, we introduce the traditional imaging function of the Kirchhoff migration and apply it for a real-time mornitoring of moving objects $D_m(t)$ generated from the scattering matrix $\mathbb{K}(t)$ such that

Unfortunately, it is very difficult to apply (2.2) directly to design an imaging function because exact expression of $\mathrm{E}_{{\rm{tot}}}^{(z)}({\bf r}, {\bf a}_p)$ is unknown. Since the condition (2.1) holds, every $D_m(t)$ can be regarded as an object with small diameter. Thus, based on ^[51], it is possible to apply the Born approximation $\mathrm{E}_{{\rm{tot}}}^{(z)}({\bf r}', {\bf a}_p)\approx \mathrm{E}_{{\rm{inc}}}^{(z)}({\bf r}', {\bf a}_p)$ so that (2.2) can be approximated as

and correspondingly, applying (3.2) to (3.1), $\mathbb{K}(t)$ can be decomposed as

Based on the structure of the above decomposition, let us generate the following unit vector: for each ${\bf r}\in\Omega$,

Now, let us define an inner product $\langle\cdot, \cdot\rangle_{\ell^2(\mathcal{A})}$ on the Lebesgue space $\ell^2(\mathcal{A})$ such that

Then, based on the orthogonal property of the $\ell^2(\mathcal{A})$, we can examine that the value of $\langle {\bf F}({\bf r}), {\bf G}({\bf r}', t)\rangle_{\ell^2(\mathcal{A})}$ will reach its maximum value when ${\bf r} = {\bf r}'\in D(t)$. Correspondingly, the following imaging function of the Kirchhoff migration can be introduced: for each ${\bf r}\in\Omega$,

Based on (3.3) and (3.4), the value of $\mathfrak{F}_{{\rm{KIR}}}({\bf r}, t)$ is expected to reach its maximum value when ${\bf r} = {\bf r}'\in D(t)$. For a detailed description of the imaging function, we refer to ^[3].

Contrary to the traditional simulation setup, each $N$ antenna is used for signal transmission, whereas the remaining $N-1$ antennas are used for signal reception. Therefore, the value of $S_{{\rm{scat}}}(p, q, t)$ for $p = q$, i.e., the diagonal elements of $\mathbb{K}(t)$, cannot be determined so that generated scattering matrix is of the following form

For related discussions, we recommend the following references ^[40,41,54]. In this study, we set $S_{{\rm{scat}}}(p, q, t) = 0$ instead of the unknown measurement data, and the following scattering matrix is considered for the design of the imaging function:

We refer to ^[40,44] for an explanation of why the diagonal elements of scattering matrix were set to zero. Although the diagonal elements are missing, the following imaging function can be introduced in a similar way to (3.5): for each ${\bf r}\in\Omega$,

The theoretical reason of the applicability of real-time object detection is discussed next.

3.2.   Structure of the imaging function

In order to explain the availability of real-time detection of a moving object, the mathematical structure of the designed imaging function is carefully explored by establishing a relationship with an infinite series of Bessel functions.

Theorem 3.1 (Structure of imaging function). Let $\boldsymbol{\theta}_n = {\bf a}_n/R = (\cos\theta_n, \sin\theta_n)^{\mathtt{T}}$ and ${\bf r}- {\bf r}' = | {\bf r}- {\bf r}'|(\cos\phi, \sin\phi)^{\mathtt{T}}$. If $| k_\mathrm{b}({\bf a}_n- {\bf r})|\geq0.25$ for all ${\bf r}\in\Omega$, then $\mathfrak{F}({\bf r}, t)$ can be represented as follows:

where $J_s$ is a Bessel function of the order $s$ of the first kind,

Proof. Since $| k_\mathrm{b}({\bf a}_n- {\bf r})|\geq0.25$ for all $n = 1, 2, \cdots, N$, the following asymptotic form of the Hankel function holds (see [12, Theorem 2.5], for instance)

Then, since

and

we can derive

where $\mathcal{N}_p = \left\{{1, 2, \cdots, N}\right\}\backslash\left\{{p}\right\}$ and $C = (e^{2iR k_\mathrm{b}}\omega \varepsilon_\mathrm{b})/(32R k_\mathrm{b}\pi)$.

Since $\boldsymbol{\theta}_n\cdot({\bf r}- {\bf r}') = | {\bf r}- {\bf r}'|\cos(\theta_n-\phi)$ and the following Jacobi-Anger expansion formula holds uniformly

we can evaluate

and correspondingly,

Now, applying (3.10) again, we can evaluate

and

Then,

Finally, by taking the absolute value, we can obtain the result (3.8). This completes the proof. □

Based on Theorem 3.1, we can briefly derive the following important result of the availability and unique determination of moving objects.

Corollary 3.1 (Availability and unique determination of moving objects). For given angular frequency $\omega$, the moving objects $D_m(t)$ at each time $t$ can be detected uniquely through the map of $\mathfrak{F}({\bf r}, t)$ with the same condition of Theorem 3.1.

Proof. Since $J_0(0) = 1$ and $J_s(0) = 0$ for $s\ne0$, the map of $\mathfrak{F}({\bf r}, t)$ will contain peaks of large magnitudes when ${\bf r} = {\bf r}'\in D(t)$. Hence, although some unexpected artifacts with small magnitudes will be included in the map of $\mathfrak{F}({\bf r}, t)$, the locations of the moving objects can be imaged uniquely through the map of $\mathfrak{F}({\bf r}, t)$. □

3.3.   Various properties of imaging function

On the basis of the result (3.8), we can explore some properties of the imaging function as follows.

Property 3.1 (Dependence of the material properties). If ${\bf r} = {\bf r}'\in D_m(t)$, then, since $J_0(k_\mathrm{b}| {\bf r}- {\bf r}'|) = 1$ and $\mathcal{E}(k_\mathrm{b}| {\bf r}- {\bf r}'|, t) = \mathcal{E}(2 k_\mathrm{b}| {\bf r}- {\bf r}'|, t) = 0$, we have

Hence, we can conclude that the value of $\mathfrak{F}({\bf r}, t)$ significantly depends on the size, permittivity, and conductivity of the object. This means that if the size, permittivity, and conductivity of an object $D_m(t)$ is considerably larger than those of the others, peaks of large magnitude will appear on the map of $\mathfrak{F}({\bf r}, t)$ at the location of this object because $\mathfrak{F}({\bf r}, t) > \mathfrak{F}({\bf r}', t)$ for ${\bf r}\in D_m(t)$ and ${\bf r}'\in D(t)\backslash\overline{D_m(t)}$.

Property 3.2 (Ideal conditions for a proper detection). Based on the structure of the factor $\mathcal{E}(k_\mathrm{b}| {\bf r}- {\bf r}'|, t)$, it disturbs the detection of objects. Thus, eliminating the factors $\mathcal{E}(k_\mathrm{b}| {\bf r}- {\bf r}'|, t)$ and $\mathcal{E}(2 k_\mathrm{b}| {\bf r}- {\bf r}'|, t)$ will guarantee good results. Notice that if one can increase $N$ as much as possible ($N\longrightarrow+\infty$) or apply extremely high frequency ($\omega\longrightarrow+\infty$), then their effects can be reduced, but this is inappropriate for real-world applications.

Property 3.3 (Practical condition for a proper detection). To reduce the adverse effect of $\mathcal{E}(k_\mathrm{b}| {\bf r}- {\bf r}'|, t)$, a condition must be found to satisfy the following equation:

Since we have no a priori information of objects, $i^sJ_{s}(k_\mathrm{b}| {\bf r}- {\bf r}'|)$ cannot be eliminated. This means that we must find a condition $\theta_n$ such that

Based on the periodic property of the cosine and sine functions, the effect of $\mathcal{E}(k_\mathrm{b}| {\bf r}- {\bf r}'|, t)$ can be reduced when an even number of antennas are uniformly distributed on a circular array and the total number of antennas $N$ is greater than $8$. This means that the array configuration of the antennas affects the imaging quality and this is the theoretical reason for the even number of antennas and their symmetric location with respect to the origin in general. We refer to ^[25,38,42] for a similar phenomenon in various imaging techniques.

Property 3.4 (On the imaging of objects close to an antenna). On the imaging of objects, the distance between each object or between objects and antennas significantly influences the imaging performance because the integral equation formula (2.2) holds for well-separated objects and Theorem 3.1 holds when $| k_\mathrm{b}({\bf r}- {\bf a}_n)|\geq0.25$. If an object $D_m(t)$ is very close to an antenna ${\bf A}_n$ at time $t$ such that

then, based on the asymptotic form of the Hankel function

some elements of $\mathbb{G}(t)$ will be very large valued so that unexpected imaging results (appearance of several artifacts with large magnitudes, invisible objects, etc.) will be obtained. Here, $\gamma = 0.57721\ldots$ denotes the Euler–Mascheroni constant.

Property 3.5 (Similarity with bifocusing method). The comparison of the derived structure (3.8) and the result in ^[27] indicates that the imaging performance of Kirchhoff migration and the bifocusing method will be very similar because their structures are almost the same.

4.   Simulation results with experimental data

In this section, various numerical simulation results with experimental data are shown for demonstrating the feasibility of the Kirchhoff migration and supporting the theoretical result. For the simulation, $N = 16$ dipole antennas equally distributed on a circle with a diameter of ${0.18}\ {{\rm{m}}}$ are placed into a cylindrical tank with a height of ${0.3}\ {{\rm{m}}}$ and a diameter of ${0.2}\ {{\rm{m}}}$ to satisfy the condition in Property 3.3. The tank was filled with water such that $(\varepsilon_\mathrm{b}, \sigma_\mathrm{b}) = (78 \varepsilon_0, 0.2 \mathrm{~S} / \mathrm{m})$ at $f = {925}\ {{\rm{MHz}}}$, where $\varepsilon_0=8.854 \times 10^{-12} \mathrm{~F} / \mathrm{m}$ denotes the vacuum permittivity.

Under the current simulation configuration, since $k_\mathrm{b}\approx171.27- 4.26i$ (i.e., $| k_\mathrm{b}|\approx171.3237$), we set the ROI as a circle centered at the origin with radius of ${0.085}\ {{\rm{m}}}$ to satisfy the condition $| k_\mathrm{b}({\bf a}_n- {\bf r})|\geq0.25$ for all $n$; refer to the condition of Theorem 3.1 and discussion of Property 3.4. Throughout this paper, the elements $S_{{\rm{scat}}}(p, q, t)$ of $\mathbb{G}(t)$ are generated using a microwave machine manufactured by Electronics and Telecommunications Research Institute (ETRI); refer to ^[28]. For describing objects, the cross-section of four long objects $D_m$, $m = 1, 2, 3, 4$, are chosen. Table 1 presents the material properties of each object and Figure 1 exhibits the manufactured microwave machine, antenna arrangements, and selected objects. All the numerical experiments were implemented on a 3.6GHz 8 Core i9-10910 CPU with 64GB RAM, running on macOS 13.6.5 and MATLAB R2023b. With these computational resources, each of the experiments took approximately 15 seconds.

Example 4.1 (Tracking of a single moving object). Figure 2 shows the maps of $\mathfrak{F}({\bf r}, t)$ for a single, moving object $D_3(t)$. Although some artifacts are included on the map, it is possible to track the moving object. Note that at some time $t$ (e.g., $t = {4.5}\ {{\rm{s}}}$, ${5.0}{\ {\rm{s}}}$, ${10.5}{\ {\rm{s}}}$, etc.), the artifacts were eliminated. We cannot explain the exact reason of this phenomenon, but we believe that at those moments, $\mathcal{E}(k_\mathrm{b}| {\bf r}- {\bf r}'|, t)\approx0$ of (3.8).

Example 4.2 [Tracking of moving objects: slightly different radii and same material properties] Figure 3 shows the maps of $\mathfrak{F}({\bf r}, t)$ for two moving objects $D_2(t)$ and $D_3(t)$ with slightly different radii and same material properties. Let us emphasize that object $D_2(t)$ has a diameter of ${6.55}{\ {\rm{mm}}}$ in the thick part and ${6.40}{\ {\rm{mm}}}$ in the narrow part, and we varied the size by moving it up and down; refer to Figure 1. Due to this reason, as we discussed in Property 3.1, the magnitudes in several maps of $\mathfrak{F}({\bf r}, t)$ during $t = {8.5}{\ {\rm{s}}}$ and $t = {14.5}{\ {\rm{s}}}$ are different such that $\mathfrak{F}({\bf r}_2, t) < \mathfrak{F}({\bf r}_3, t)$ because the size of cross-section of $D_2(t)$ is slightly smaller than that of $D_3(t)$. Similar to the result in Example 4.1, the track of the moving objects can be recognized even though some artifacts degrade the imaging quality. Fortunately, in contrast to Example 4.1, the artifacts do not disturb the recognition of the moving objects, but at certain instances $t = {1.5}{\ {\rm{s}}}$ and $t = {2.5}{\ {\rm{s}}}$, significantly large amounts of artifacts appear. Unfortunately, the truth behind this phenomenon cannot be estimated at present, but the distance between the object and antennas may possibly be the reason; refer to Property 3.4.

Example 4.3 (Tracking of moving objects: different radii and material properties). Now we consider the tracking of moving objects $D_1(t)$ and $D_2(t)$ with different radii and material properties. On the basis of the result in Figure 4, we can observe that a peak of large magnitude appears when ${\bf r}\in D_1(t)$ at $t = {5.0}{\ {\rm{s}}}$, ${5.5}{\ {\rm{s}}}$, ${12.0}{\ {\rm{s}}}$, etc., because the size of $D_1(t)$ is larger than that of $D_2(t)$, i.e., $\mathfrak{F}({\bf r}_1, t) > \mathfrak{F}({\bf r}_2, t)$ for ${\bf r}_1\in D_1(t)$ and ${\bf r}_2\in D_2(t)$; refer to Property 3.1. Similar to the results in Example 4.2, some artifacts are included on the map of $\mathfrak{F}({\bf r}, t)$; however, they disappeared at certain instances, e.g., $t ={1.0}{\ {\rm{s}}}$, ${4.0}{\ {\rm{s}}}$, and ${5.5}{\ {\rm{s}}}$.

Example 4.4 (Tracking of moving objects: same radii but different material properties). For the final example, we present the simulation result for moving objects $D_2(t)$ and $D_4(t)$ with the same size but different material properties such that $\mathfrak{F}({\bf r}_2, t)\gg\mathfrak{F}({\bf r}_4, t)$ for ${\bf r}_2\in D_2(t)$ and ${\bf r}_4\in D_4(t)$. Figure 5 shows the maps of $\mathfrak{F}({\bf r}, t)$, and we can observe that, in contrast to Example 4.2, it is impossible to recognize the movement of $D_4(t)$ because the value of $\mathfrak{F}({\bf r}, t)$ at ${\bf r}\in D_4(t)$ is significantly smaller than that of $\mathfrak{F}({\bf r}, t)$ at ${\bf r}\in D_2(t)$. It is interesting to observe that when an object moves quickly, few ghost replicas are exhibited on the map, e.g., two peaks of large magnitude appeared at $t = {8.0}{\ {\rm{s}}}$ when $D_2(t)$ moved quickly.

5.   Conclusions and perspectives

Owing to the existence of small, moving objects and the non-iterative Kirchhoff migration technique in inverse scattering problems, we designed a real-time algorithm for imaging the moving objects on the basis of the representation formula of the scattering parameters. To examine the feasibility and explore some properties of the designed algorithm, we proved that the imaging function could be represented as the total and array configuration of the antennas and an infinite series of Bessel functions of integer order.

Numerical simulations were performed using experimental data generated by the ETRI to demonstrate that the Kirchhoff migration is very effective for real-time detection of moving objects in microwave imaging. However, the algorithm's application is currently restricted to the detection of small objects; therefore, further applications to the detection of the movements of arbitrary shaped extended objects or the evolution of crack-like defects must be considered. Forthcoming studies will focus on designing appropriate imaging algorithms, performing related mathematical analyses, and conducting related simulations. Following several studies ^[9,49,58,60], imaging techniques based on terahertz wave propagation are promising and effective. We believe that extension to the terahertz imaging and establishment of reliable mathematical theories will be an interesting research subject. Finally, inspired by several remarkable studies ^{[14,55,59,61]}, development of machine learning based inversion techniques with experimental datasets will be a remarkable research topic.

Use of AI tools declaration

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

Acknowledgments

This research was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. NRF-2020R1A2C1A01005221), the Institute of Information & Communications Technology Planning & Evaluation (IITP) grant funded by the Korea government (MSIT) (No. 2021-0-00731), and the Soonchunhyang University research fund.

Conflict of interest

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