Accurate Phase Calibration of Multistatic Imaging System for Medical and Industrial Applications

Abstract

Multistatic imaging systems are commonly used in radar systems and microwave imaging. In these systems, many antennas are arranged three-dimensionally and connected to RF switches. The length of each transmitter (Tx) and receiver (Rx) channel differs slightly, causing artifacts in high-resolution image reconstruction. This study presents a novel method for the phase calibration of multistatic systems. This method does not require system reconstruction and can automatically perform phase calibration in a short time. This method is expected to facilitate an accurate phase measurement in multistatic systems. The approach involves phase calibration by analyzing the reflection coefficients of antenna elements in the time domain. Imaging experiments were performed on a multistatic imaging system using this calibration method, and the position and shape of a metal rod with a diameter one-fourth of a wavelength were reconstructed by simple back-projection with an accuracy beyond the diffraction limit.

1. Introduction

Multistatic systems are an important technique for estimating the position and characteristics of a target object using multiple receivers. Multistatic systems are used in various fields such as radar and imaging. In the radar field, multistatic systems are used for multistatic radar and multistatic Synthetic Aperture Radar (SAR) [1]. In the imaging field, medical and industrial applications are expected. Research is being conducted in the medical field for the early diagnosis of breast cancer [2,3] and brain strokes [4,5]. When applied in the industrial field, it is used for the inspection of foreign objects in food [6] and security screening [7,8,9]. In these fields, a substantial amount of data from different perspectives of the same object are required for imaging with high sensitivity and accuracy. Multistatic systems are used in many studies [10,11,12] because they can efficiently switch between many Tx and Rx antenna combinations and reduce data acquisition time.

However, a challenge in the use of multistatic systems is the calibration of the system. The solid-state RF switches and coaxial cables used in these systems have different electrical lengths between ports due to mechanical tolerances. In addition, in actual use, phase drift occurs due to temperature changes and cable bends. For these reasons, the multistatic system needs calibration for each port. However, it is difficult to calibrate each port one by one without system reconfiguration because the reconfiguration takes a very long time and causes small phase shifts due to the unavoidable additional bending of RF cables. In addition, these systems are often measured in the frequency domain and need to be calibrated to take dispersion into account. Various methods of calibrating multistatic systems have therefore been devised. The calibration methods can be divided into those using the scattered field [13,14,15] and using Tx and Rx lines [16,17,18,19].

Majid Ostadrahimi et al. performed calibration using scattered fields by objects inside the antenna array [13]. Calibration coefficients were calculated using the electric field scattered by a phantom placed inside the antenna array and simulated data. In this study, the incident field radiated from the antenna was also used for the calibration. The calibration with the scattering field gave better image accuracy than the calibration with the incident field. It is noted that such a calibration method using simulation could be used for other systems.

Tegan Counts et al. calibrated 6ch ground-penetrating radar by connecting the Tx and Rx ports directly with a cable [17]. In this method, calibration was performed for delay and attenuation in the signal lines as well as direct coupling between the antennas. First, the antennas were disconnected from the cables, and each Tx-Rx channel was connected by a cable. A vector network analyzer (VNA) was then used to measure the response of the cables. And the attenuation and delay in the cables and the RF switches were calibrated. Next, the antennas were connected to the cables, and the free space bistatic responses were measured. This allowed for the calibration of the direct coupling between the antennas and the reflection from objects near the antenna array.

Manuel Kasper et al. performed a calibration of a multistatic system using an electronic calibration (ECal) module and a phantom with known dielectric properties [16]. The calibration process consisted of three steps. First, the reflection coefficient of one port of the antenna array was calibrated. The ECal module was connected between the VNA, and the antenna, and the standard reflection calibration procedure was performed. This procedure moved the calibration plane of the VNA from the VNA port to the connector of the ECal module. Next, a phantom was placed inside the antenna array, and then a calibration of the other ports was performed using the reflection coefficient of the antennas. Finally, the coupling between the antennas was calibrated using the transmission coefficients between each Tx and Rx antenna pair. It should be noted that this calibration method has limitations, such as requiring symmetrical antenna configurations and not using a RF switches. However, because of the low effort required for calibration, it is mentioned that it is possible to perform calibration over a short period.

Based on these previous studies, this study proposed a method that enables automatic and quick calibration without system reconfiguration. The phase calibration without system reconfiguration reduces the time and effort required for calibration. In addition, errors are caused by the tightness of connectors and bends in cables. Therefore, the phase calibration can also reduce measurement errors. As a result, it is expected that the facilitation of accurate phase measurements in multistatic systems can be easily performed. This is expected to reduce the cost of image reconstruction calculations, improve image accuracy, and reduce artifacts in reconstructed images, especially in imaging systems. The phase calibration was performed because of the large phase contribution to the image reconstruction in the image reconstruction method used in this study. After calibrating the reference port using the ECal module, the reflection coefficients were analyzed in the time domain to calculate the electrical length difference in each port, and phase calibration was performed. Direct coupling of the Tx and Rx antennas and reflections at the antenna array were removed during reconstruction. Using the phase-calibrated system, image reconstruction experiments were performed using the simple back-projection method of waves.

2. Phase Calibration Method

In this study, phase calibration was performed using the difference in electrical length. The difference in electrical length was determined by the reflection coefficients from the measurement reference plane of the antenna element. The measurement reference plane is the position where the signal is reflected, such as at the conversion section of a transmission line. If there is a signal reflection, the position can be defined as a measurement reference plane for phase calibration. This calibration method was designed for a system like the one shown in Figure 1. The system consists of a VNA, RF switches, an ECal module, and antenna elements. ECal is an automated calibration module for VNA that simplifies traditional manual calibration procedures and improves accuracy and efficiency [20]. While traditional calibration requires multiple references (e.g., open, short, loaded, through, etc.) to be prepared and performed manually, ECal allows these references to be switched electronically within the module for a quicker calibration. The RF switches require one port to be connected to the ECal module, in addition to the ports used for measurements. This port is called the reference port. RF switches and coaxial cables have different electrical lengths for each channel. In Figure 1, the electrical length difference with the reference port is indicated by $l_{\mathrm{t}\mathrm{x}}\left(\omega , n_{\mathrm{t}\mathrm{x}}\right)$ and $l_{\mathrm{r}\mathrm{x}}\left(\omega , n_{\mathrm{r}\mathrm{x}}\right)$, while $\omega$ represents the angular frequency, and $n_{\mathrm{t}\mathrm{x}}$ and $n_{\mathrm{r}\mathrm{x}}$ denote the numbers of the Tx and Rx ports. If the electrical length of the port is longer than the electrical length of the reference port, $l_{\mathrm{t}\mathrm{x}}\left(\omega , n_{\mathrm{t}\mathrm{x}}\right)$ and $l_{\mathrm{r}\mathrm{x}}\left(\omega ,n_{\mathrm{r}\mathrm{x}}\right)$ are represented by a positive number. The green dotted line shows the calibration plane of the Ecal module; if only the ECal module is used for calibration, the VNA calibration plane is the green dotted line. The orange dotted line indicates the measurement reference plane of each channel. The calibration plane and the measurement reference plane are not the same for each channel, even if the reference port is configured with the ECal module. The phase calibration is performed by determining the electrical length deviation from the reflection coefficient in such a system. The electrical length is calculated by taking dispersion into account.

The calibration procedure is shown in Figure 2. The phase calibration is performed in three steps. First, the reference port connected to the ECal module is calibrated using ECal. The ECal module is connected only to the reference port for its calibration, and the reflection coefficients of other ports are measured by VNA in the frequency domain. In the second step, the reflection signal at the measurement reference plane is extracted from the measured reflection coefficients. The phase change is then obtained, considering the variance of each port. In the third step, phase calibration is performed using the phase change due to the difference in electrical length.

2.1. Reflection Coefficient Measurements

The reference ports of the RF switches are calibrated using the ECal module. After calibration with the ECal module, the calibration plane of the VNA does not match with the measurement reference plane due to different electrical lengths for different ports. Next, the reflection coefficients $S_{11}\left(\omega ,n_{\mathrm{t}\mathrm{x}}\right)$ and $S_{22}\left(\omega ,n_{\mathrm{r}\mathrm{x}}\right)$ of both the Tx and Rx ports are measured. This reflected signal includes the reflected signal from the phase origin, which indicates the difference in electrical length from the reference port. This signal is used for the phase calibration. However, the reflected signal also includes reflections at connector junctions (e.g., switch connections, or antenna connections). In the next step, the reflected signal at the phase origin is extracted from the reflected signal.

2.2. Time Domain Gating

The measured reflection signals are transformed into the time domain to extract only the reflection signals at the phase origin. This technique is called Time Domain Gating [21]. Initially, the discrete inverse Fourier transform is used to convert the frequency domain reflection signals $S_{11}\left(\omega ,n_{\mathrm{t}\mathrm{x}}\right)$ into the time domain signals $s_{11}\left(t,n_{\mathrm{t}\mathrm{x}}\right)$.

$$s_{11}\left(t\left[m_{\mathrm{t}}\right],n_{\mathrm{t}\mathrm{x}}\right)={\displaystyle \frac{1}{M}}\sum _{m_{\mathrm{f}}=0}^{M-1}{S}_{11}\left(\omega \left[m_{\mathrm{f}}\right],n_{\mathrm{t}\mathrm{x}}\right){\mathrm{e}}^{j{\displaystyle \frac{2\pi }{M}}{m}_{\mathrm{t}}{m}_{\mathrm{f}}}$$

(1)

where $M$ is the number of measurement data in the frequency domain; $m_{\mathrm{f}}$ is the discrete frequency index $m_{\mathrm{f}}=0, 1, 2, \dots ,M-1$; $m_{\mathrm{t}}$ is the discrete time index $m_t=0, 1, 2, \dots ,M-1;$ $t$ is the time axis in the time domain $t\left[m_{\mathrm{t}}\right]=\Delta t·m_{\mathrm{t}}$; and $\Delta t$ is the time resolution of the time axis.

$$\Delta t={\displaystyle \frac{2\pi }{{\omega }_{\mathrm{s}\mathrm{t}\mathrm{o}\mathrm{p}}-{\omega }_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{t}}}}$$

(2)

where ${\omega }_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{t}}$ is the lowest angular frequency ${\omega }_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{t}}=\omega \left[0\right]$ and ${\omega }_{\mathrm{s}\mathrm{t}\mathrm{o}\mathrm{p}}$ is the highest angular frequency ${\omega }_{\mathrm{s}\mathrm{t}\mathrm{o}\mathrm{p}}=\omega [M-1]$. To improve the accuracy of phase calibration, it is necessary to improve the time resolution $\Delta t$. The reflection coefficients of the Rx $s_{22}\left(t,n_{\mathrm{r}\mathrm{x}}\right)$ are also calculated the same way. The time resolution $\Delta t$ depends on the bandwidth in the frequency domain ${\omega }_{\mathrm{s}\mathrm{t}\mathrm{o}\mathrm{p}}-{\omega }_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{t}}$. Since it is difficult to ensure sufficient bandwidth in the experiments, the bandwidth was ensured by padding the data. In this case, a method of adding zeros to the measured frequency domain data to widen the bandwidth is used.

Next, the reflection signals at the measurement reference plane ${s^{\prime }}_{11}\left(t,n_{\mathrm{t}\mathrm{x}}\right)$ and $s_{22}^{\prime }\left(t,n_{\mathrm{r}\mathrm{x}}\right)$ are extracted from the time domain signals. The signals in the time domain were extracted by applying a gate function.

$${s^{\prime }}_{11}\begin{array}{c}\left(t\left[m_{\mathrm{t}}\right],n_{\mathrm{t}\mathrm{x}}\right)=s_{11}\left(t\left[m_{\mathrm{t}}\right],n_{\mathrm{t}\mathrm{x}}\right)·\mathrm{W}\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{o}\mathrm{w}\left(m_{\mathrm{t}}\right)\end{array}$$

(3)

The Kaiser–Bessel function was used for gate processing. The Kaiser–Bessel function is represented by the following equation:

$$\begin{array}{c}\mathrm{W}\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{o}\mathrm{w}\left(m_{\mathrm{t}}\right)=\left\{\begin{array}{cc}{\displaystyle \frac{{\mathrm{I}}_0\left(\beta \sqrt{1-\frac{4m_{\mathrm{t}}^2}{{\left(m_{\mathrm{t}\_\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}}-1\right)}^2}}\right)}{{\mathrm{I}}_0\left(\beta \right)}}& \mathrm{if} -{\displaystyle \frac{m_{\mathrm{t}\_\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}}-1}{2}}+{\mathrm{m}}_{\mathrm{t}\_\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}\le m_{\mathrm{t}}\le {\displaystyle \frac{m_{\mathrm{t}\_\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}}-1}{2}}+m_{\mathrm{t}\_\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}\\ 0& \mathrm{otherwise}\end{array}\right.\end{array}$$

(4)

where $m_{\mathrm{t}\_\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}$ is a time index of the center of the peak, and $m_{\mathrm{t}\_\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}}$ is the index number of the peak width. ${\mathrm{I}}_0$ is the modified zeroth-order Bessel function. $\beta$ is the shape parameter of the function. When $\beta = 0$, the function is rectangular. As $\beta$ increases, the function’s half-width decreases. $\beta =6$, like the Han function, which is commonly used.

Finally, the discrete Fourier transform was used to convert the time domain signal into the frequency domain:

$${S^{\prime }}_{11}\left({\omega }^{\prime }\left[{m_{\mathrm{f}}}^{\prime }\right],n_{\mathrm{t}\mathrm{x}}\right)=\sum _{m_{\mathrm{t}}=0}^{M-1}{s^{\prime }}_{11}\left(t\left[m_{\mathrm{t}}\right],n_{\mathrm{t}\mathrm{x}}\right)e^{-j{\displaystyle \frac{2\pi }{M}}{m}_{\mathrm{t}}{m}_{\mathrm{f}}\prime }$$

(5)

where $m_{\mathrm{f}}^{\prime }$ is a new discrete frequency index $m_{\mathrm{f}}^{\prime }=0,1,2,\dots ,M-1$, and $\omega ’$ is the new frequency axis in the frequency domain ${\omega }^{\prime }\left[m_{\mathrm{f}}^{\prime }\right]=\Delta {\omega }^{\prime }·m_{\mathrm{f}}^{\prime }+{\omega }_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{t}}$. $\Delta {\omega }^{\prime }$ is the frequency resolution of the time axis $\Delta {\omega }^{\prime }={\displaystyle \frac{2\pi }{\Delta t·(M-1)}}$. Since zeros remain in the transformed frequency domain due to data padding, the data were cut to match the frequency domain used for the measurement. Also, since the frequency domain’s tick width is different from the $\Delta f$ of the measured reflection coefficient, then linear interpolation yields ${S^{\prime }}_{11}\left(\omega ,n_{\mathrm{t}\mathrm{x}}\right)$. ${S^{\prime }}_{22}\left(\omega ,n_{\mathrm{r}\mathrm{x}}\right)$ is also calculated.

2.3. Apply Phase Calibration

The round-trip phase is converted to a one-way phase for phase calibration. When converting the phase, it is necessary to consider the phase fringes. When the electrical length is longer than the wavelength, the measured phase has advanced by more than 360 degrees. When determining the one-way phase, it is necessary to have an absolute phase. To find the absolute phase, the phase wraps $p$, which is the number of times the phase has cycled $2\pi$ radians at the lower frequency limit, is required. There are two ways to obtain the phase wraps at lower frequency.

The first method is to adjust the lower frequency limit to $p=0$. If the wavelength is sufficiently long compared to the assumed electrical length, there is no need to consider phase fringes. The second method is to calculate $p$. If the assumed electrical length is known, it is possible to evaluate the absolute phase from the relative phase. The electrical length can be obtained from the peak of the reflection coefficient in the time domain, but attention must be paid to conditions such as the frequency to be used for the inverse Fourier transform, as it cannot consider dispersion. The absolute phase ${\varphi }_{\mathrm{t}\mathrm{x}}\left(\omega ,n_{\mathrm{t}\mathrm{x}}\right)$ is calculated using $p$:

$${\varphi }_{\mathrm{t}\mathrm{x}}\left(\omega ,n_{\mathrm{t}\mathrm{x}}\right)=\arg \left(S_{11}^{\prime }\left(\omega ,n_{\mathrm{t}\mathrm{x}}\right)\right)+2\pi p$$

(6)

The round-trip absolute phase ${{\varphi }^{\prime }}_{\mathrm{t}\mathrm{x}}$ is converted to the one-way phase for phase calibration. Dividing the unwrapped phase value by two yields the one-way phase. ${\varphi }^{{\prime }_{\mathrm{t}\mathrm{x}}}\left(\omega ,n_{\mathrm{t}\mathrm{x}}\right)$ represents the calibration data. ${\varphi }^{{\prime }_{\mathrm{r}\mathrm{x}}}\left(\omega , n_{\mathrm{t}\mathrm{x}}\right)$ is also calculated in the same way.

$${\varphi }^{{\prime }_{\mathrm{t}\mathrm{x}}}\left(\omega , n_{\mathrm{t}\mathrm{x}}\right) ={\displaystyle \frac{{\varphi }_{\mathrm{t}\mathrm{x}}\left(\omega , n_{\mathrm{t}\mathrm{x}}\right)}{2}}$$

(7)

Finally, calibration data are applied to the measured data using De-Embedding [22,23]. The phase shift is eliminated using the difference in electrical lengths determined from the transmission coefficients measured between the calibration planes. The circuit between the calibration planes in Figure 1 is connected in series, with three circuits as shown in Figure 3.

$S\left(\omega ,{ n}_{\mathrm{t}\mathrm{x}}, n_{\mathrm{r}\mathrm{x}}\right)$ represents the data between the measurement reference planes of the Tx and Rx antennas. $S_{\mathrm{t}\mathrm{x}}\left(\omega ,n_{\mathrm{t}\mathrm{x}}\right)$ is the S-parameter between the calibration plane and the measurement reference plane of the Tx antenna side, and $S_{\mathrm{r}\mathrm{x}}\left(\omega , n_{\mathrm{t}\mathrm{x}}\right)$ is the S-parameter between the calibration plane and the measurement reference plane of the Rx antenna side. Assuming that these two circuits are lossless, they are represented by the following S-parameters:

$${\mathrm{S}}_{\mathrm{t}\mathrm{x}}\left(\omega ,n_{\mathrm{t}\mathrm{x}}\right)=\left[\begin{array}{cc}0& e^{j{\varphi }_{\mathrm{t}\mathrm{x}}^{\prime }\left(\omega ,{\mathrm{n}}_{\mathrm{t}\mathrm{x}}\right)}\\ e^{j{\varphi }_{\mathrm{t}\mathrm{x}}^{\prime }\left(\omega ,{\mathrm{n}}_{\mathrm{t}\mathrm{x}}\right)}& 0\end{array}\right], {\mathrm{S}}_{\mathrm{r}\mathrm{x}}\left(\omega ,{\mathrm{n}}_{\mathrm{t}}\right)=\left[\begin{array}{cc}0& e^{j{\varphi }_{\mathrm{r}\mathrm{x}}^{\prime }\left(\omega ,n_{\mathrm{r}\mathrm{x}}\right)}\\ e^{j{\varphi }_{\mathrm{r}\mathrm{x}}^{\prime }\left(\omega ,{\mathrm{n}}_{rx}\right)}& 0\end{array}\right]$$

(8)

These circuits are expressed using T-parameters, as follows:

$${\mathrm{T}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}\left(\omega ,n_{\mathrm{t}\mathrm{x}}, n_{\mathrm{r}\mathrm{x}}\right) = {\mathrm{T}}_{\mathrm{t}\mathrm{x}}\left(\omega , n_{\mathrm{t}\mathrm{x}}\right)\mathrm{T}\left(\omega , n_{\mathrm{t}\mathrm{x}}, n_{\mathrm{r}\mathrm{x}}\right) {\mathrm{T}}_{\mathrm{r}\mathrm{x}}\left(\omega , {\mathrm{n}}_{\mathrm{t}\mathrm{x}}\right)$$

(9)

S-parameters and T-parameters are transformed by the following equations:

$$\left[S\right]={\displaystyle \frac{1}{T_{22}}}\left[\begin{array}{cc}{T}_{12}& T_{11}{T}_{22} - T_{12}{T}_{21}\\ 1& -T_{21}\end{array}\right]$$

(10)

$$\left[T\right]={\displaystyle \frac{1}{S_{21}}}\left[\begin{array}{cc}{S}_{12}{S}_{21}-S_{11}{S}_{22}& S_{11}\\ -S_{22}& 1\end{array}\right]$$

(11)

Multiply the transmission coefficient ${\mathrm{T}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}\left(\omega , n_{\mathrm{t}\mathrm{x}}, n_{\mathrm{r}\mathrm{x}}\right)$ between the calibration planes by the inverse of ${\mathrm{T}}_{\mathrm{t}\mathrm{x}}\left(\omega ,n_{\mathrm{t}\mathrm{x}}\right)$ and ${\mathrm{T}}_{\mathrm{r}\mathrm{x}}\left(\omega ,n_{\mathrm{r}\mathrm{x}}\right)$ to obtain the transmission coefficient $\mathrm{T}\left(\omega ,n_{\mathrm{t}\mathrm{x}},n_{\mathrm{r}\mathrm{x}}\right)$ between the measurement reference planes.

$$\mathrm{T}\left(\omega ,n_{\mathrm{t}\mathrm{x}},n_{\mathrm{r}\mathrm{x}}\right)={\mathrm{T}}_{\mathrm{t}\mathrm{x}}^{-1}\left(\omega ,n_{\mathrm{t}\mathrm{x}}\right) {\mathrm{T}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}\left(\omega ,n_{\mathrm{t}\mathrm{x}}, n_{\mathrm{r}\mathrm{x}}\right) {\mathrm{T}}_{\mathrm{r}\mathrm{x}}^{-1}\left(\omega ,n_{\mathrm{r}\mathrm{x}}\right)$$

(12)

The following equation is the S-parameter between phase centers calculated from the S-parameter between calibration planes:

$$\mathrm{S}\left(\omega ,n_{\mathrm{t}\mathrm{x}},n_{\mathrm{r}\mathrm{x}}\right)=\left[\begin{array}{cc}{S}_{11}\left(\omega , n_{\mathrm{t}\mathrm{x}}\right)·e^{-2j{\varphi }_{\mathrm{t}\mathrm{x}}^{\prime }\left(\omega ,n_{\mathrm{t}\mathrm{x}}\right)}& S_{12}\left(\omega ,n_{\mathrm{t}\mathrm{x}},n_{\mathrm{r}\mathrm{x}}\right)·e^{-j\left({\varphi }_{\mathrm{t}\mathrm{x}}^{\prime }\left(\omega ,n_{\mathrm{t}\mathrm{x}}\right)+{\varphi }_{\mathrm{r}\mathrm{x}}^{\prime }\left(\omega ,n_{rx}\right)\right)}\\ S_{21}\left(\omega ,n_{\mathrm{t}\mathrm{x}},n_{\mathrm{r}\mathrm{x}}\right)·e^{-j\left({\varphi }_{\mathrm{t}\mathrm{x}}^{\prime }\left(\omega ,n_{\mathrm{t}\mathrm{x}}\right)+{\varphi }_{\mathrm{r}\mathrm{x}}^{\prime }\left(\omega ,n_{\mathrm{r}\mathrm{x}}\right)\right)}& S_{22}\left(\omega , n_{\mathrm{r}\mathrm{x}}\right)·e^{-2j{\varphi }_{\mathrm{r}\mathrm{x}}^{\prime }\left(\omega ,n_{\mathrm{r}\mathrm{x}}\right)}\end{array}\right]$$

(13)

3. Experiment

The microwave imaging system used in this study is illustrated in Figure 4. This device consists of a VNA, RF switches, an antenna array, and coaxial cables. The two ports of the VNA (Keysight P9371A, 300 kHz to 6.5 GHz) are connected to Tx and Rx switches. Each switch port is connected to antenna elements fixed on the antenna array and an ECal module (Keysight N7551A). The RF switches for Tx and Rx have 32 ports each. The Tx switch consists of an SP4T (Single-Pole Four-Throw, Mini-Circuit USB-SP4T-63) and four SP8T (Mini-Circuit USB-1SP8T-63H) switches, and the Rx switch consists of an SP4T and two SP16T (Mini-Circuit USB-1SP16T-83H) switches. The RF switch is connected to antennas and the ECal module via 3 m coaxial cables (CRYSTEK CCCSMA18-MM-086F-120). The antenna array is cylindrical in shape, with antenna elements fixed along the circumference. The antenna array has an outer diameter of 150 mm and an inner diameter of 280 mm, with a height of 85 mm. This antenna array consists of three rows of antennas. Ten Tx antennas and ten Rx antennas are fixed per row. The Tx and Rx antennas are arranged alternately every 18 degrees.

The Tx and Rx antennas used are balanced Vivaldi antennas (BVAs) [24,25]. This type of antenna has a small error in phase measurement and is thus suitable for this study. A schematic of a BVA is shown in Figure 5. The BVA is a three-layer structure consisting of two ground layers and a feed line layer. A connector and the feed line are sandwiched and implemented between a two-layer Printed Circuit Board (PCB). The feed line of this antenna is the stlip line (SL) that is sandwiched between both ground layers. The microwave propagated through the SL is converted into a transmission line, propagating double-sided slotline, i.e., the bilateral slotline (BSL). By exponentially widening the width of the BSL, the microwave is irradiated. The BVA, by changing the SL into the BSL, offers the following advantages [25,26]:

• Using SL to BSL reduces the amount of radiation from the feed line.
• BSL has less dispersion compared to the slotline.

Figure 5. Balanced Vivaldi antenna (BVA): (a) exploded view of the BVA construction; (b) geometry and parameters of the BVA.

Figure 5. Balanced Vivaldi antenna (BVA): (a) exploded view of the BVA construction; (b) geometry and parameters of the BVA.

In this study’s phase calibration, we confirm the reflection coefficients of the BVA. Figure 6a shows the phase of the reflection coefficients of the BVA. The graph has frequency on the x-axis and phase on the y-axis. In the frequency sweep from 1.0 to 6.5 GHz, the phase decreases at a constant rate. Then, it is considered that the position of reflections in the BVA is constant regardless of frequency.

Next, we examine the positions of reflections within the BVA. Since it was assumed that signals reflect at the transmission line transitions of the SL and BSL shown in Figure 3, we checked the reflection coefficients of five different BVAs with varying BSL widths, $d$. Figure 6b illustrates the result of transforming the BVA’s $S_{11}$ into the time domain using discrete inverse Fourier transform. The x-axis represents time, and the y-axis represents reflection intensity, with peaks appearing at the positions where signals reflect. The time at which the first peak appears is $t_1$ and the time at which the second peak appears is $t_2$.

The reflection coefficients of the five antennas showed two peaks. The peak at $t_1$ matched the delay time when a signal entered from the BSL connector and was reflected at the transmission line transition. Comparing the peak heights, the peak at $t_1$ varied with the width of the BSL; the narrower the BSL, the higher the peak height. On the other hand, the peak at $t_2$ remained the same height regardless of the BSL width. From these results, it was inferred that the peak at $t_1$ represented signals reflecting at the transmission line transition of the BSL. The antenna is suitable for the phase calibration method used in this study because it allows us to calculate the distance from the antenna’s reflection coefficients to the phase origin.

The phase calibration was performed using the measured reflection coefficients. Python was used to calculate the calibration data. The discrete inverse Fourier transform, and the discrete Fourier transform were calculated using the ifft and the fft functions from the Numpy library. Kaiser–Bessel functions were computed using the Kaiser function in the signal module of the Scipy library.

4. Calibration Result

Figure 7 shows the transmission coefficients measured at the antenna array. The transmission coefficients were measured at ten pairs of opposite Tx and Rx antennas in the antenna array. To reduce multipath reflections within the antenna array, the cylindrical cavity in the center of the antenna array is filled with the same resin as the antenna array. Figure 7a shows the coefficients among the calibration planes before phase calibration, while Figure 7b shows the coefficients among measurement reference planes after phase calibration. These graphs have frequency on the horizontal axis and phase on the vertical axis. In Figure 7a, the slopes of the phase waveforms for frequency sweeps are different for each channel. This phase is measured between calibration planes, and the reference planes are different from each channel. Therefore, the slopes of the phase waveforms are different from each channel. In Figure 7b, where phase correction is applied, the slopes of the phase waveforms match. The slope of the phase waveform decreased at a constant rate with frequency sweep. Around 3.0 GHz and 6.0 GHz, the phase increased and then decreased again. There were also fringe jumps in phase around 3.0 GHz in some channels. These were thought to be due to multiple reflections between transmitting and receiving antennas in the central cavity of the antenna array.

The phase difference between channels of the calibrated phase was less than 8° at 1.5–6.5 GHz. The phase difference between channels was obtained by unwrapping the calibrated phase and calculating the standard deviation from the difference with one channel. As the phase waveform contains ripples, the difference was calculated using the 100 MHz wide moving average phase, and the phase difference around 3 GHz was large, but this was due to fringe jumps. The standard deviation was calculated by subtracting the transmission coefficients of 10 pairs of transmitting and receiving antennas at each frequency. The standard deviation at 1.5~6.5 GHz was about 8 degrees. The phase difference between channels was approximately 1/45 of the wavelength, confirming a highly accurate phase calibration.

5. Imaging Result

The differences in reconstructed images with and without phase calibration were compared using the simple back-projection method of waves. The simple back-projection method of waves reconstructs the object from the amplitude and phase of the scattered waves measured at the Rx antenna. Microwaves emitted from a Tx antenna are scattered by objects, and the scattered waves spread with the objects as the wave source. The scattered electric field is measured with the Rx antenna so that the space containing the wave source is surrounded by a closed surface. Calculating the time reversal of this measured scattered electric field converges to the wave source. In this way, the object that is the wave source is back-projected. This method is based on the Cauchy problem [27], and back-projection by time reversal using pulses has been used in many studies [28,29,30]. The reconstructed image of the object is obtained by subtracting the reconstructed images calculated for the presence and absence of the object, respectively.

Figure 8 shows the reconstructed image of a metal rod placed in the cylindrical cavity in the center of the antenna array by the simple back-projection method of waves. The diameter of the metal rod is 12 mm. This diameter is less than 1/4 of the wavelength at the 6.5 GHz frequency used for reconstruction and is beyond the diffraction limit. The frequency used for the reconstructed images was 6.5 GHz and the tomograms were taken at a height of 37.5 mm from the bottom of the antenna array. This height of 37.5 mm corresponds to the midpoint of the antenna array. Figure 8a,b are calculated from data without phase calibration, and Figure 8c,d are calculated from data with phase calibration. Figure 8a,c show the amplitude distributions, while Figure 8b,d show the phase distributions. The horizontal axis of the reconstructed image shows the x-axis, and the vertical axis shows the y-axis. The white solid circles in the reconstructed image show the shape of the antenna array. The outer circle is the outer edge of the antenna array, and the inner circle is the inner edge of the array. The color bar in the reconstructed image of amplitude is proportional to the amplitude of the scattered wave, with higher values indicating greater scatter. The color bar in the reconstructed image of the phase shows the phase. The color bars of the amplitude distribution are set about the reconstructed image in the presence of the object.

The phase calibration enabled the identification of the object from the reconstructed image using the simple inverse projection method of the wave motion. Initially, the differences in amplitude distribution with and without the phase calibration are compared. In the absence of the object, the object could not be identified from the reconstructed image. In contrast, the object could be seen in the reconstructed image with phase calibration. The areas with large values appeared at the position where the object was placed. The diameter and shape of the object in the reconstructed image were almost identical to the actual object. Artifacts appeared around the object in a circular area centered on the object. These artifacts were due to multiple reflections in the metal and on the inner surface of the antenna array. In the region where $y>0$, artifacts appeared on the body of the antenna array. This artifact appears at the location of the Rx antenna and is considered to have been generated during the subtraction of the reconstructed image. Next, the phase distribution of the reconstructed image is compared. The phase distribution without phase calibration is asymmetric for the y-axis. In contrast, when phase calibration is performed, the distribution is symmetrical about the y-axis. Symmetrical phase distributions were obtained by performing phase calibration on all channels.

Figure 9 shows a 3D representation of the reconstructed images of Figure 8 at different heights. The graph shows that the amplitude distribution and the frequency used for reconstruction is 6.5 GHz. The x- and the y-axes show the same axes as in Figure 8. The z-axis shows the height at which the reconstructed image was calculated, from the bottom of the antenna array to a height of 75 mm, divided at 12.5 mm intervals. The color bar values indicate the intensity of scattering, and a transparency is set to the color of the color map for visualization. The black circles in the reconstructed image at each height indicate the outer diameter of the antenna array, the inner diameter, and the position of the metal column. The location and shape of the metal columns were confirmed and are presented in Figure 9. The metal columns are visible in the reconstructed images at four different heights from 25 mm to 62.5 mm. In the reconstructed images at 0 mm, 12.5 mm, and 75 mm, the metal columns are obscured. Several previous studies have shown that measurements at the Tx and Rx antennas located in different planes (cross-planes) are required to improve the image accuracy of microwave imaging systems [31,32,33,34]. It has been reported that measurements that do not include cross-planes reduce the vertical resolution of the reconstructed image. The fact that some of the reconstructed images in this study are blurred may indicate the same tendency as in the aforementioned previous studies.

6. Conclusions

The phase calibration in a multistatic system has been developed. Since no system reconfiguration is required, the system can be calibrated automatically, and the phase calibration can be performed in a short time. The phase calibration was performed by analyzing the reflection coefficient of the antenna in the time domain. The difference in electrical length of each port was determined. In the phase-calibrated system, the phase difference between the ports was found to be approximately 1/45 of the wavelength. An image reconstruction experiment was performed on the calibrated system using a simple inverse projection of the waves. The image reconstructed at 6.5 GHz showed the location and shape of the metal rod, which is one-fourth the diameter of the wavelength and beyond the diffraction limit.