*Article* **An Imaging Plane Calibration Method for MIMO Radar Imaging**

#### **Yuanyue Guo \*, Bo Yuan, Zhaohui Wang and Rui Xia**

Key Laboratory of Electromagnetic Space Information, Chinese Academy of Sciences,

University of Science and Technology of China, Hefei 230026, China; yuanb@mail.ustc.edu.cn (B.Y.);

wzh01@mail.ustc.edu.cn (Z.W.); xrke928@mail.ustc.edu.cn (R.X.)

**\*** Correspondence: yuanyueg@ustc.edu.cn

Received: 10 November 2019; Accepted: 27 November 2019; Published: 29 November 2019

**Abstract:** In two dimensional cross-range multiple-input multiple-output radar imaging for aerial targets, due to the non-cooperative movement of the targets, the estimated imaging plane parameters, namely the center and the posture angles of the imaging plane, may have deviations from true values, which defocus the final image. This problem is called imaging plane mismatch in this paper. Focusing on this problem, firstly the deviations of spatial spectrum fulfilling region caused by imaging plane mismatch is analyzed, as well as the errors of the corresponding spatial spectral values. Thereupon, the calibration operation is deduced when the imaging plane parameters are accurately obtained. Afterwards, an imaging plane calibration algorithm is proposed to utilize particle swarm optimization to search out the imaging plane parameters. Finally, it is demonstrated through simulations that the proposed algorithm can accurately estimate the imaging plane parameters and achieve good image focusing performance.

**Keywords:** two-dimensional radar imaging; multiple-input multiple-output (MIMO) radar; Particle Swarm Optimization (PSO); imaging plane calibration algorithm (IPCA)

#### **1. Introduction**

Aerial targets imaging is an important research direction in the field of radar imaging technology. Especially, it plays a crucial role in the military field, such as, aerial defense [1] and anti-missile defense [2] and so on. Multiple-input multiple-output (MIMO) radar is a new radar technique, which adopts multiple transmitters and receivers. By transmitting orthogonal space–time block codes [3] or frequency diversity signals [4], a MIMO radar with *M* transmitters and *N* receivers can eventually form a virtual array of the aperture length up to *M* times that of the receive array, which greatly saves the hardware cost. In addition, MIMO radar can obtain the images of the aerial targets with only one snapshot, and thus has enormous superiority in radar image acquisition time [5,6].

Since MIMO radar technique was proposed, it has been desired to build high performance imaging algorithms. The researches of MIMO radar imaging algorithms mainly focus on two aspects: The first is the wave-number domain imaging methods, and their imaging performances mainly depend on the spatial spectrum fulfilling region, which is generally required to be uniformly fulfilled, so that the fast Fourier transform (FFT) can be applied [7,8]. In this regard, Prof. Yarovoy et al. have conducted a large number of studies and verified the feasibility of the proposed algorithms in security check [9], wall penetrating [10] and ground penetrating imaging applications [7]. In addition to the traditional wave-number domain methods, the iterative optimization methods have also been applied in the MIMO radar imaging field. For example, Prof. Li's team of University of Florida proposed several iterative imaging algorithms, such as iterative adaptive approaches (IAA) algorithm [11] and the sparse learning via iterative minimization (SLIM) [12] algorithm, which have been well applied and verified in MIMO radar scenarios.

In MIMO radar imaging, the model mismatch caused by system errors or array spatial position errors will degrade the imaging quality. Therefore, the study of model error calibration algorithms is an important research direction of high-quality MIMO radar imaging. In this regard, a large number of studies have been used to calibrate the phase error [13,14], carrier frequency deviation [15], array position error [16], off-grid problem [17,18], and so on. The degradation of MIMO radar resolution under the condition of phase error from the perspective of point spread function (PSF) is analyzed in [13], and the sparse imaging via expectation maximization (SIEM) algorithm is proposed, which alternately estimates the phase errors and the target image, and obtains better imaging quality. Subsequently, the degradation of MIMO radar resolution under the condition of carrier frequency deviation from the perspective of PSF is analyzed in [15] as well, and an iterative algorithm employing iteration strategy similar to reference [13] is proposed, and good imaging results are achieved. MIMO imaging with array position errors is studied in reference [16], while in reference [17,18], the off-grid problem of MIMO radar imaging is studied. The algorithms proposed in [16–18] all employ sparse optimization by alternately estimating the target image and the errors during iterations, hence clear images are finally obtained.

Most of the above methods are proposed for two dimensional (2D) imaging in range and cross-range directions. However, for the target plane parallel to cross-range direction, these methods cannot be directly applied. In fact, it is necessary to set the origin of the coordinates at the center of the imaging plane [5,19] for the 2D cross-range imaging methods based on the spatial spectrum. In practical applications, as the target is non-cooperative, it is required to estimate the scene center and posture angles of the target plane, so as to make the final image focus on the image scene center and the target plane. However, there are always some deviations between the estimated imaging plane parameters and the real situations, which resulting in unfocused image and poor imaging quality. This problem is called the imaging plane mismatch problem in this paper.

To solve this problem, firstly, the deviations of spatial spectral fulfilling region caused by imaging plane mismatch are analyzed, and the location errors between estimated spatial spectral point and real spatial spectral point under the condition of imaging plane mismatch are deduced, as well as the errors of the corresponding spatial spectral values. Subsequently, in order to estimate imaging plane parameters and to be able to calibrate the locations and values of spatial spectral points, an imaging plane calibration algorithm (IPCA) is proposed. Aiming at minimizing the image entropy as well as promoting target sparsity, IPCA utilizes a particle swarm optimization (PSO) [20,21] algorithm to search out the parameters of imaging plane center deviation and pose angles deviations, and then calibrates the locations of spatial spectral points according to these parameters, so as to obtain images with better quality.

This paper is organized as follows. Section 2 introduces the spatial spectral imaging model of MIMO radar, and analyzes the problem of imaging plane mismatch, and the deviations between the estimated spatial spectral point and the real spatial spectral point position under the imaging plane mismatch are deduced. Section 3 provides the design and the detailed flow of IPCA. In Section 4, the validity of the proposed algorithm, robustness to noise, and tolerance to mismatching parameters are verified by simulations. Section 5 is the conclusion of this paper.

#### **2. Problem Formulation of Imaging Plane Mismatch**

In this section, the spatial spectral imaging model of MIMO radar is reviewed firstly, and then the model mismatch problem caused by the imaging plane mismatch is analyzed. Under the condition of imaging plane mismatch, the deviations between the positions of the obtained spatial spectral points and the positions of the real spectral points are analyzed, as well as the values of the corresponding spatial spectral points. Afterwards, the calibration operation is deduced to obtain the focused image after obtaining the imaging plane parameters. Notice that the radar system discussed here is the frequency diversity MIMO (f-MIMO) radar [4,22].

#### *2.1. Space Spectral Imaging Model of Multiple-Input Multiple-Output (MIMO) Radar*

In this subsection, the spatial spectral imaging model of MIMO radar is reviewed firstly.

Figure 1 illustrates the general 2D cross-range imaging scenario of MIMO radar. Let (*x*, *y*, *z*) be Cartesian coordinates with the origin *O* located at the center of the imaging plane and the 2D target is supposed to be on the imaging plane. The location of the *p*-th transmitting antenna and the *q*-th receiving antenna are denoted as *r<sup>p</sup>* = *rp*, *θp*, *ϕ<sup>p</sup>* and *r<sup>q</sup>* = *rq*, *θq*, *ϕ<sup>q</sup>* in spherical coordinate respectively.

**Figure 1.** Space spectral imaging model of multiple-input multiple-output (MIMO) radar.

Without loss of generality, regardless of the loss associated with the free-space propagation, the echo signal at the *q*-th receiver by the *p*-th transmitter is given by [5]:

$$s\_{p,q}(t) = \iint \sigma(\mathbf{x}\_T, y\_T) \exp\left(j2\pi \left[f\_p\left(t - \frac{R\_{p,T}}{c} - \frac{R\_{q,T}}{c}\right)\right]\right) d\mathbf{x}\_T d y\_T \tag{1}$$

where *σ*(*xT*, *yT*) denotes the reflectivity of the scatterer at (*xT*, *yT*) on the imaging plane, *fp* is the transmitting frequency of the *<sup>p</sup>*-th transmitter and *<sup>c</sup>* is the speed of light. *Rp*,*<sup>T</sup>* = *r<sup>p</sup>* <sup>−</sup> *<sup>r</sup><sup>T</sup>* , *Rq*,*<sup>T</sup>* = *r<sup>q</sup>* <sup>−</sup> *<sup>r</sup><sup>T</sup>* , *<sup>r</sup><sup>T</sup>* = (*xT*, *yT*, 0).

Then, down conversion is applied to the received signal, which can be achieved by multiplying it by the following reference signal:

$$s\_{ref}(t) = \exp\left(-j2\pi \left[f\_p\left(t - \frac{\left|r\_p\right|}{c} - \frac{\left|r\_q\right|}{c}\right)\right]\right) \tag{2}$$

In the case of the far field, approximate conditions can be used:

$$\begin{aligned} \left| R\_{p,T} = \left| r\_p - r\_T \right| = \left| r\_p \right| - r\_T \cdot \mathfrak{e}\_p \\ \left| R\_{q,T} = \left| r\_q - r\_T \right| = \left| r\_q \right| - r\_T \cdot \mathfrak{e}\_q \end{aligned} \right. \tag{3}$$

where ˆ*e<sup>p</sup>* <sup>=</sup> *<sup>r</sup>p*/|*rp*<sup>|</sup> and ˆ*e<sup>q</sup>* <sup>=</sup> *<sup>r</sup>q*/|*rq*|.

Thus it can get:

$$\begin{split} s\_{p,q}(t) \cdot s\_{ref}(t) &= \iint \sigma(\mathbf{x}\_{T}, \mathbf{y}\_{T}) \exp\left(j2\pi \left[f\_{p}\left(t - \frac{R\_{p,T}}{c} - \frac{R\_{q,T}}{c}\right)\right]\right) \\ &\times \exp\left(-j2\pi \left[f\_{p}\left(t - \frac{|r\_{p}|}{c} - \frac{|r\_{q}|}{c}\right)\right]\right) dx\_{T} dy\_{T} \\ &= \iint \sigma\left(\mathbf{x}\_{T}, \mathbf{y}\_{T}\right) \exp\left(j2\pi \mathbf{K}\_{p,q} \cdot \mathbf{r}\_{T}\right) dx\_{T} dy\_{T} \end{split} \tag{4}$$

where *Kp*,*<sup>q</sup>* = *kx <sup>p</sup>*,*q*, *k y p*,*q* . *k<sup>x</sup> <sup>p</sup>*,*<sup>q</sup>* and *k y <sup>p</sup>*,*<sup>q</sup>* are:

$$\begin{aligned} k\_{p,q}^x &= \frac{f\_p}{c} \left( \cos \theta\_p \sin \varphi\_p + \cos \theta\_q \sin \varphi\_q \right) \\ k\_{p,q}^y &= \frac{f\_p}{c} \left( \cos \theta\_p \cos \varphi\_p + \cos \theta\_q \cos \varphi\_q \right) \end{aligned} \tag{5}$$

Therefore, the value of the 2D spatial spectral point (*k<sup>x</sup> <sup>p</sup>*,*q*, *k y <sup>p</sup>*,*q*) of the imaging plane is obtained:

$$\mathcal{G}\left(k\_{p,q}^x, k\_{p,q}^y\right) = \iint \sigma\left(\mathbf{x}\_T, y\_T\right) \exp\left(j2\pi\left(\mathbf{x}\_T k\_{p,q}^x + y\_T k\_{p,q}^y\right)\right) d\mathbf{x}\_T d y\_T \tag{6}$$

Finally, the common algorithms can be applied on spatial spectral point value to obtain the target image, such as the inverse fast Fourier transform (IFFT) algorithm, the back-projection (BP) algorithm [23], and the non-uniform fast Fourier transform [8] and so on.

#### *2.2. Analysis of Model Mismatch Problem Caused by Imaging Plane Mismatch*

In the actual MIMO radar imaging applications, especially in aerial target imaging, the center and the posture angles of the imaging plane are uncertain due to the target's non-cooperative movement state. When the parameters of the imaging plane have deviations, the fulfilling region of the obtained spatial spectrum will deviate from the real region, which will cause images to be unfocused.

Figure 2 shows the imaging geometry of the scenario with imaging plane mismatch. Let *α* denote the estimated target plane, while let *β* denote the real target plane. Besides, set up the coordinate *O* − *xαyαz<sup>α</sup>* with the origin *O* located at the center of the estimated target plane and plane *xαOy<sup>α</sup>* coincides with the estimated target plane. The location vectors of *p*-th transmitting antenna and *q*-th receiving antenna are *r<sup>p</sup>* = *rp*, *θp*, *ϕ<sup>p</sup>* and *r<sup>q</sup>* = *rq*, *θq*, *ϕ<sup>q</sup>* respectively in Coordinate *O* − *xαyαzα*. Likewise, set up the coordinate *O* − *xβyβz<sup>β</sup>* with the origin *O* located at the center of the real target plane and plane *xβOy<sup>β</sup>* coincides with the real target plane. Meanwhile, the location vectors of *p*-th transmitting antenna and *q*-th receiving antenna are *r <sup>p</sup>* = *r <sup>p</sup>*, *θ <sup>p</sup>*, *ϕ p* and *r <sup>q</sup>* = *r <sup>q</sup>*, *θ <sup>q</sup>*, *ϕ q* respectively in coordinate *O* − *xβyβzβ*. In addition, the changes between plane *β* and plane *α* consist of the translation change and the posture angles' change. The direction vector *d* = *dβ*, *θβ*, *ϕβ* represents the translation of the origin of coordinate *O* − *xβyβz<sup>β</sup>* relative to the origin of coordinate *O* − *xαyαzα*. Define (*δ*, *μ*, *ξ*) as the posture angle of plane *β* in Coordinate *O* − *xαyαzα*, where *δ* denotes the angle between the axis *x<sup>β</sup>* of coordinate *O* − *xβyβz<sup>β</sup>* and the plane *xαOyα*, *μ* denotes the angle between the projection of the axis *x<sup>β</sup>* on the plane *xαOz<sup>α</sup>* and the axis *xα*, and *ξ* denotes the angle between the plane *yβO z<sup>β</sup>* and the plane *yαOzα*.

What is more, *θ <sup>p</sup>*, *θ <sup>q</sup>*, *ϕ <sup>p</sup>*, *ϕ <sup>q</sup>* and *θp*, *θq*, *ϕp*, *ϕ<sup>q</sup>* have relations:

$$\begin{aligned} \theta'\_p &= \theta\_p + \delta \sin \varphi\_p + \mu \cos \varphi\_p, & \quad \varphi'\_p &= \varphi\_p + \xi + \delta \cos \varphi\_p + \mu \sin \varphi\_p \\ \theta'\_q &= \theta\_\emptyset + \delta \sin \varphi\_\emptyset + \mu \cos \varphi\_\emptyset, & \quad \varphi'\_q &= \varphi\_\emptyset + \xi + \delta \cos \varphi\_\emptyset + \mu \sin \varphi\_\emptyset \end{aligned} \tag{7}$$

**Figure 2.** Imaging geometry of MIMO radar with imaging plane mismatch.

Hence, *r <sup>p</sup>* and *<sup>r</sup> <sup>q</sup>* can be computed by *<sup>r</sup>p*, *<sup>r</sup><sup>q</sup>* and *<sup>d</sup>*, namely

$$\begin{aligned} \left|r'\_p\right| &= \left|r\_p - d\right| = \sqrt{\left|r\_p\right|^2 + d\_\beta^2 - 2\left|r\_p\right| \cdot d\_\beta \left(\cos\theta\_p \cos\theta\_\beta \cos\left(\varphi\_p - \varphi\_\beta\right) + \sin\theta\_p \sin\theta\_\beta\right)} \\ \left|r'\_q\right| &= \left|r\_q - d\right| = \sqrt{\left|r\_q\right|^2 + d\_\beta^2 - 2\left|r\_q\right| \cdot d\_\beta \left(\cos\theta\_q \cos\theta\_\beta \cos\left(\varphi\_q - \varphi\_\beta\right) + \sin\theta\_q \sin\theta\_\beta\right)} \end{aligned} \tag{8}$$

Substitute Equation (8) into Equation (3), so that it can get:

$$\begin{aligned} \tau\_{p,T} &= \frac{R\_{p,T}}{c} = \frac{\left|\left.r\_p'\right| + x\_T \cos\theta\_p' \sin\theta\_p' + y\_T \cos\theta\_p' \cos\varphi'\right|}{c} \\ &= \frac{\sqrt{\left|\left.r\_p'\right|^2 + d\_{\beta}^2 - 2\left|\left.r\_p\right| \, d\_{\beta}\left(\cos\theta\_p \cos\theta\_{\beta} \cos\left(\varphi\_p - \varphi\_{\beta}\right) + \sin\theta\_p \sin\theta\_{\beta}\right)\right|}{c}}{+\frac{x\_T \cos\theta\_p' \sin\varphi'\_p + y\_T \cos\theta\_p' \cos\varphi'\_p}{c}} \\ &+ \frac{x\_T \cos\theta\_p' \sin\varphi'\_p + y\_T \cos\theta\_p' \cos\varphi'\_p}{c} \\ \tau\_{q,T} &= \frac{R\_{q,T}}{c} = \frac{\left|\left.r\_q'\right| + x\_T \cos\theta\_q' \sin\varphi'\_q + y\_T \cos\theta\_q' \cos\varphi'\_q\right|}{c} \\ &= \frac{\sqrt{\left|\left.r\_q\right|^2 + d\_{\beta}^2 - 2\left|\left.r\_q\right| \, d\_{\beta}\left(\cos\theta\_q \cos\theta\_{\beta} \cos\left(\varphi\_q - \varphi\_{\beta}\right) + \sin\theta\_q \sin\theta\_{\beta}\right)\right|}{c}}{+\frac{x\_T \cos\theta\_q' \sin\varphi'\_q + y\_T \cos\theta\_q' \cos\varphi'\_q}{c}} \\ \end{aligned} (9)$$

Since the real imaging plane parameters are not accurately known, the reference signal *sref*(*t*) is still the same as Equation (2). Substitute Equation (9) into Equation (4) and using far field approximation:

$$\begin{split} s\_{p,q}(t) \cdot s\_{\pi f}(t) &= \iint \sigma \left( \mathbf{x}\_{\Gamma \tau} y\_{\Gamma \tau} \right) \exp \left\{ j2\pi \frac{f\_p}{c} \left( -d\_{\beta} (\cos \theta\_p \cos \theta\_{\beta} \cos \left( \varphi\_p - \varphi\_{\beta} \right) + \sin \theta\_p \sin \theta\_{\beta} \right) \right. \\ &\left. - d\_{\beta} (\cos \theta\_q \cos \theta\_{\beta} \cos \left( \varphi\_q - \varphi\_{\beta} \right) + \sin \theta\_q \sin \theta\_{\beta}) \right. \\ &\left. + \mathbf{x}\_T \left[ \cos \left( \theta\_p + \delta \sin \varphi\_p + \mu \cos \varphi\_p \right) \sin \left( \varphi\_p + \xi + \delta \cos \varphi\_p + \mu \sin \varphi\_p \right) \right. \\ &\left. + \cos \left( \theta\_q + \delta \sin \varphi\_q + \mu \cos \varphi\_q \right) \sin \left( \varphi\_q + \xi + \delta \cos \varphi\_q + \mu \sin \varphi\_p \right) \right] \\ &\left. + \mathbf{y}\_T \left[ \cos \left( \theta\_p + \delta \sin \varphi\_p + \mu \cos \varphi\_p \right) \cos \left( \varphi\_p + \xi + \delta \cos \varphi\_p + \mu \sin \varphi\_p \right) \right] \right] \right\} \, d\mathbf{x}\_T t d\_T \\ &\left. + \cos \left( \theta\_q + \delta \sin \varphi\_q + \mu \cos \varphi\_q \right) \cos \left( \varphi\_q + \xi + \delta \cos \varphi\_q + \mu \sin \varphi\_q \right) \right) \right] \, d\mathbf{x}\_T t d\_T \end{split}$$

Therefore Equation (6) actually gets:

$$\begin{aligned} G\left(\mathbf{x}\_{p,q}^{k}, \mathbf{y}\_{p,q}^{y}\right) &= \iint \sigma\left(\mathbf{x}\_{T}, \mathbf{y}\_{T}\right) \exp\left\{j2\pi\left(\mathbf{x}\_{T}\mathbf{x}\_{p,q}^{x} + \mathbf{y}\_{T}\mathbf{x}\_{p,q}^{y}\right) + j2\pi\frac{f\_{p}}{c}\right\} \\ &- d\_{\beta}(\cos\theta\_{p}\cos\theta\_{\beta}\cos\left(\varphi\_{p} - \varphi\_{\beta}\right) + \sin\theta\_{p}\sin\theta\_{\beta}) \\ &- d\_{\beta}(\cos\theta\_{q}\cos\theta\_{\beta}\cos\left(\varphi\_{q} - \varphi\_{\beta}\right) + \sin\theta\_{q}\sin\theta\_{\beta}) \\ &+ \mathbf{x}\_{T}\Big[\sin\left(\delta\sin\varphi\_{p} + \mu\cos\varphi\_{p}\right)\left(\sin\theta\_{p}\sin\varphi\_{p} + \sin\theta\_{q}\sin\varphi\_{q}\right) \\ &+ \sin\left(\beta + \delta\cos\varphi\_{q} + \mu\sin\varphi\_{q}\right)\left(\cos\theta\_{p}\cos\varphi\_{p} + \cos\theta\_{q}\cos\varphi\_{q}\right)\Big] \\ &+ \mathbf{y}\_{T}\Big[\sin\left(\delta\sin\varphi\_{p} + \mu\cos\varphi\_{p}\right)\left(\sin\theta\_{p}\cos\varphi\_{p} + \sin\theta\_{q}\cos\varphi\_{q}\right)\Big] \\ &+ \sin\left(\beta + \delta\cos\varphi\_{q} + \mu\sin\varphi\_{q}\right)\left(\cos\theta\_{p}\sin\varphi\_{p} + \cos\theta\_{q}\sin\varphi\_{q}\right)\Big]\Big\}\Big{]}\,\mathrm{d}\_{T}\mathrm{d}\_{l}\mathrm{y}\_{T}\end{aligned}$$

The location of space spectral points obtained by the above processing is not completely matched with the real spatial spectral fulfilling region. As a result, the entire image is not focused on the real target plane.

#### *2.3. The Calibration Operation*

In order to obtain more accurate imaging results, it is necessary to search out the parameters of the imaging plane, so as to eliminate the mismatch between the estimated spatial spectral fulfilling region and the real spatial spectral fulfilling region, and the according spatial spectral values. After the parameters of the imaging plane are obtained, the following calibration operation can be applied to achieve focused image.

When the parameters *d* = *dβ*, *θβ*, *ϕβ* and (*δ*, *μ*, *ξ*) have been obtained, the reference signal can be calibrated as:

$$s\_{ref}'\left(t, \mathbf{r}\_{p'}', \mathbf{r}\_q'\right) = \exp\left\{-j2\pi f\_p \left(t - \frac{\left|\mathbf{r}\_p'\right| + \left|\mathbf{r}\_q'\right|}{c}\right)\right\} \tag{12}$$

After the coherent processing in Equation (10), the final spatial spectral value is:

$$\mathbf{G}'\left(k\_{p,q}^{\prime x}, k\_{p,q}^{\prime y}\right) = \iint\_{T} \sigma\left(\mathbf{x}\_{T}, y\_{T}\right) \mathbf{e}^{j2\pi\left(\mathbf{x}\cdot\mathbf{k}\_{p,q}^{\prime x} + y\_{T}\mathbf{k}\_{p,q}^{\prime y}\right)} d\mathbf{x}\_{T} d\mathbf{y}\_{T} \tag{13}$$

where *T* denotes the imaging area and *k<sup>x</sup> <sup>p</sup>*,*q*, *k y <sup>p</sup>*,*<sup>q</sup>* are defined in Equation (14).

$$\begin{aligned} k\_{p,q}^{I\star} &= \frac{f\_p \left( \cos \theta\_p' \sin \varphi \theta\_p' + \cos \theta\_q' \sin \varphi \theta\_q' \right)}{c} \\ &= \frac{f\_p}{c} \left( \left( \cos \left( \theta\_p + \delta \sin \varphi\_p + \mu \cos \varphi\_p \right) \sin \left( \varphi\_p + \xi + \delta \cos \varphi\_p + \mu \sin \varphi\_p \right) \right) \\ &+ \cos \left( \theta\_q + \delta \sin \varphi\_q + \mu \cos \varphi\_q \right) \sin \left( \varphi\_q + \xi + \delta \cos \varphi\_q + \mu \sin \varphi\_q \right) \right) \\ &k\_{p,q}^{\vartheta} = \frac{f\_p \left( \cos \theta\_p' \cos \varphi\_p' + \cos \theta\_q' \cos \varphi\_q' \right)}{c} \\ &= \frac{f\_p}{c} \left( \cos \left( \theta\_p + \delta \sin \varphi\_p + \mu \cos \varphi\_p \right) \cos \left( \varphi\_p + \xi + \delta \cos \varphi\_p + \mu \sin \varphi\_p \right) \right) \\ &+ \cos \left( \theta\_q + \delta \sin \varphi\_q + \mu \cos \varphi\_q \right) \cos \left( \varphi\_q + \xi + \delta \cos \varphi\_q + \mu \sin \varphi\_q \right) \end{aligned} \tag{14}$$

Ultimately, the focused images can be obtained using *k<sup>x</sup> <sup>p</sup>*,*q*, *k y <sup>p</sup>*,*<sup>q</sup>* and the corresponding spatial spectral value *G k x <sup>p</sup>*,*q*, *k y p*,*q* .

#### **3. The Proposed Imaging Plane Calibration Algorithm (IPCA)**

When there are deviations between the estimated values of the imaging plane parameters (*dβ*, *θβ*, *ϕβ*, *δ*, *μ*, *ξ*) and the real values, it will lead to the problem of spatial spectrum mismatch. The mismatch of spatial spectrum results in the image not focusing at the real imaging plane and the whole image quality is deteriorated seriously.

In order to achieve better imaging performance, it is necessary to search out the real imaging plane parameters, namely the six parameters (*dβ*, *θβ*, *ϕβ*, *δ*, *μ*, *ξ*).

However, searching directions cannot be clear at beginning and always changed in the processing. Commonly, a good objective function will make the searching process in desirable directions.

Generally, image entropy can be used to measure the focusing performance of an image. In radar imaging, it is defined as follows:

$$E = -\sum\_{k=0}^{M-1} \sum\_{n=0}^{N-1} \frac{|\sigma\left(\mathbf{x}\_{k'} y\_n\right)|^2}{S} \ln \frac{|\sigma\left(\mathbf{x}\_{k'} y\_n\right)|^2}{S} \tag{15}$$

where, *S* = ∑*M*−<sup>1</sup> *<sup>k</sup>*=<sup>0</sup> <sup>∑</sup>*N*−<sup>1</sup> *<sup>n</sup>*=<sup>0</sup> <sup>|</sup>*<sup>σ</sup>* (*xk*, *yn*)<sup>|</sup> 2 , and the image is discretized into *M* × *N* grids. *σ*(*xk*, *yn*) denotes the scattering coefficient of the grid in *k*-th row at *n*-th column. At the same time, in the aerial imaging applications, the targets usually have sparse characteristics.

Therefore, both the focusing performance and the sparsity of the targets are considered, hence the objective function is to minimize the following function:

$$f(\mathbf{d}, \delta, \mu, \xi) = -\sum\_{k=0}^{M-1} \sum\_{n=0}^{N-1} \frac{|\sigma\left(\mathbf{x}\_k, y\_n\right)|^2}{S} \ln \frac{|\sigma\left(\mathbf{x}\_k, y\_n\right)|^2}{S} + \gamma \sum\_{k=0}^{M-1} \sum\_{n=0}^{N-1} |\sigma\left(\mathbf{x}\_k, y\_n\right)|\_1 \tag{16}$$

where *<sup>σ</sup>* (*xk*, *yn*) = *IFFT Gp*,*q kx*, *ky*, *d*, *δ*, *μ*, *ξ* and *γ* is a tunable parameter to adjust the weights of the image entropy and the sparsity.

In such a way, the estimation problem for (*d*, *δ*, *μ*, *ξ*) can be transformed into the following optimization problem:

$$\begin{aligned} \left(d, \delta, \mu, \xi\right) &= \underset{d, \delta, \mu, \xi}{\arg\min} \left(d, \delta, \mu, \xi\right) \\ \text{s.t.} \quad \sigma\left(\mathbf{x}\_{k}, y\_{n}\right) &= \text{IFFT}\left(\mathbf{G}\_{p, \emptyset}\left(\mathbf{k}\_{\mathbf{x}}, \mathbf{k}\_{\mathbf{y}}, \mathbf{d}, \delta, \mu, \xi\right)\right) \end{aligned} \tag{17}$$

However, the above optimization problem is not a convex problem. Commonly, among non-convex optimization algorithms, metaheuristic optimization [24] is becoming increasingly popular. Particle swarm optimization (PSO) algorithm is one of the most popular metaheuristic optimization algorithms, and has been successfully applied in radar imaging problems [25,26]. Hence, the imaging plane calibration algorithm (IPCA) is proposed to utilize PSO algorithm for solving the optimization problem in Equation (17). In IPCA, each particle *p<sup>i</sup>* = (*dβi*, *θβi*, *ϕβi*, *δi*, *μi*, *ξi*) denotes one estimation of (*dβ*, *θβ*, *ϕβ*, *δ*, *μ*, *ξ*), while the velocity *v<sup>i</sup>* of each particle denotes one search direction. In PSO, each particle remembers its personal best position. Meanwhile, global best position is also recorded. In each iteration, the velocity of each particle is updated combining the individual movement states and the group movement states. Particles acquire their new positions by constantly updating their speed, and eventually all particles will converge to the global optimal value [20,21].

The detailed algorithm of IPCA is described in Algorithm 1 and the flow chart is illustrated in Figure 3.

**Figure 3.** Flow chart of imaging plane calibration algorithm (IPCA).

#### **Algorithm 1:** IPCA


```
2 for k = 1 : Kmax do
3 for i = 1 : I do
4 compute the fitness of each particle: fitness(i) = f(pk
                                                            i );
5 if fitness(i) < fitness(Pbest(i)) then
 6 Pbest(i) = pk
                         i
7 end
8 end
9 Update Gbest with Gbest = argmin
                                    pi
                                        fitness(i), i = 1, 2, . . . , I;
10 Update the location and velocity of each particle using the following equations where
       wk, c1, c2 are parameters controlling the iteration step and randk
                                                                    1 and randk
                                                                               2 are random
       numbers in (0, 1).
                vk+1
                 i =wkvk
                           i + c1 · randk
                                       1 ·

                                          Pbest(i) − pk
                                                       i

                                                         + c2 · randk
                                                                    2 ·

                                                                       Gbest − pk
                                                                                 i

                pk+1
                 i =pk
                         i + vk+1
                             i
11 if Convergence condition reached then
12 break
13 end
14 end
  Output: (dβ, θβ, ϕβ, δ, μ, ξ) and the image σ
```
#### **4. Simulations**

In this section, several simulations are carried out to verify the effectiveness of IPCA.

The imaging distance is set to 1 km in the simulations and the target is a B727 airplane with its point scattering model shown in Figure 4a. The size of the imaging plane, approximately parallel to the radar antenna array, is 60 × 60 m. The radar consists of 31 × 31 transmitters and 4 receivers and all of the antennas are located on the same transceiver plane, as is illustrated in Figure 4b. Besides, the size of the radar array is 60 × 60 m and the radar works at C-band. The detailed simulation parameters are given in Table 1.

**Table 1.** Simulation parameters.


**Figure 4.** Target and antenna array (**a**) Scattering points distribution of the B727 airplane, (**b**) Locations of transmit antennas and receive antennas.

#### *4.1. Imaging Simulation*

To verify the effectiveness of the proposed method, the following simulation was carried out. System parameters in the simulations are shown in Table 1. The parameters of the center and pose angles of the imaging plane both have errors. The real parameters are given in Table 2.



According to the derivation in Section 2 and the deviation parameters of the imaging plane, the estimated and real fulfilling region of the spatial spectrum are illustrated in Figure 5a,b respectively. It can be seen that when there are deviations of the imaging plane parameters, the estimated fulfilling region of the spatial spectrum is quite different from the fulfilling region of the real spatial spectrum. The recovered image is illustrated in Figure 5c when deviations of the imaging plane parameters are not known, and the adopted algorithm is the IFFT algorithm. It can been seen that even when there are slight deviations of the imaging plane parameters, e.g., the deviations of the posture angles of the imaging plane are within 1◦, the image is badly defocused, and the target's contour is not clear.

**Figure 5.** Spatial spectrum and reconstructed image. (**a**) The estimated fulfilled region of the spatial spectrum, (**b**) The real fulfilled region of the spatial spectrum, (**c**) Image reconstructed using estimated spatial spectrum values.

In the following, the proposed IPCA is taken to search out the six parameters, namely *dβ*, *θβ*, *ϕβ*, *δ*, *μ*, *ξ* , and then the calibration operation is taken to obtain the final image. The number of

the particles, *I*, in IPCA is set to 100 and the maximum number of iterations, *Kmax*, is set to 100 too. Meanwhile some other parameters in PSO are set: *<sup>ω</sup><sup>k</sup>* <sup>=</sup> 0.8 <sup>−</sup> 0.5 <sup>×</sup> *<sup>k</sup> Kmax* , *c*<sup>1</sup> = *c*<sup>2</sup> = 1.

According to Figure 6, the image has clearer target outlines and better focusing performance with fewer noisy points around the target. In the final image, each scattering point of the target can be clearly identified. Meanwhile, the image entropy is smaller using IPCA than that using IFFT, as is shown in Table 3, which means better focusing performance.

**Figure 6.** Reconstructed images of different iterations, (**a**–**f**) images of the 10th, 20th, 30th, 50th, 80th, and 100th iteration respectively.


**Table 3.** Image entropy.

Figure 7 shows the value of the objective function and the image entropy during the IPCA iterations. It can been seen that when the number of iterations exceeds 80, the value of the objective function tends to be stable, that is, the search results gradually converge.

From what has been discussed above, the method proposed in this paper can effectively converge. Meanwhile, the entropy of the inversion image is lower and the image quality is higher.

Figure 8a shows the estimated *d* = (*dx*, *dy*, *dz*) during the iterations. It can be seen that, with the increase of the number of iterations, *dz* gradually approaches the real value, and finally converges to the real value. The values of *dx dy* are still deviated from the real values. This is because the deviations *dx*, *dy* of the image will only make the image translation in the imaging plane, which have no effect on the image focus effect and target sparsity. Therefore, *dx*, *dy* do not affect the imaging quality and IPCA cannot guarantee *dx*, *dy* convergence to the real value. Likewise, Figure 8b shows the estimated imaging plane posture angles (*δ*, *μ*, *ξ*) with the number of iterations. It can be seen that as the number of iterations increases, the value of *δ*, *μ* converges to the real value while *ξ* not. This is because the estimation error of the rotation angle *ξ*, whose rotating axis parallel to the line of sight, will make the image rotate in the imaging plane, which has no influence on the image entropy and target sparsity. Therefore *ξ* does not affect the imaging quality and IPCA cannot guarantee *ξ* convergence to the real value.

**Figure 7.** Value of objective function and image entropy during iterations. (**a**) Value of objective function during iterations (**b**) Image entropy during iterations.

**Figure 8.** Values of *d* = (*dx*, *dy*, *dz*) and (*δ*, *μ*, *ξ*) during iterations. (**a**) *d* = (*dx*, *dy*, *dz*), (**b**) (*δ*, *μ*, *ξ*).

Meanwhile, since the PSO algorithm is a metaheuristic optimization algorithm, the computation time of IPCA is hard to predict. Therefore, 10 Monte Carlo trials are taken to count the computation time. The computation times of 10 Monte Carlo trials vary from 954 s to 2017 s with average computation time is 1617 s.

To sum up, the algorithm proposed in this paper can obtain more accurate imaging plane parameters, resulting in better image focusing performance and better imaging quality.

#### *4.2. Simulations with Different Tunable Parameter γ*

In Equation (16), a tunable parameter *γ* is defined to adjust the weights of the image entropy and the sparsity. In order to analysis the influence of *γ*, the below simulations are taken.

The simulation parameters are set the same with that in Tables 1 and 2. *γ* is set to [0.1, 0.5, 1, 3, 5, 10, 50]. The imaging results are illustrated in Figure 9, meanwhile the image entropy are given in Table 4.

It can be seen from Figure 9 that the reconstructed images are all well focused and the image entropy are all lower than 4 when *γ* is chosen with different values. So it can conclude that *γ* has seldom influence on the final imaging results. Therefore, *γ* can be selected from a wide range in IPCA.

**Figure 9.** Reconstructed images with different *γ*, (**a**) *γ* = 0.1, (**b**) *γ* = 1, (**c**) *γ* = 3, (**d**) *γ* = 5, (**e**) *γ* = 10, (**f**) *γ* = 50.

**Table 4.** Image entropy with different *γ*.


#### *4.3. Simulations under Different Signal-to-Noise Ratios (SNRs)*

In order to verify the robustness to the noise of IPCA, the following simulation is carried out in different echo signal-to-noise ratio. The echo SNRs are set from 5 dB to 30 dB. The simulation parameters are the same with that in Tables 1 and 2. The simulation results are in Figure 10.

**Figure 10.** (**a**) Value of the objective function *f*(*d*, *δ*, *μ*, *ξ*) during iterations under different SNRs, (**b**) Image entropy during iterations under different SNRs.

The value of the objective function *f*(*d*, *δ*, *μ*, *ξ*) during the iterations is illustrated in Figure 10a. As it is shown, the value of the objective function *f*(*d*, *δ*, *μ*, *ξ*) decreases iteratively and finally converges. When the SNR is above 15 dB, the final values of *f*(*d*, *δ*, *μ*, *ξ*) are closed to each other.

The image entropy during the iterations is illustrated in Figure 10b. It can be seen that image entropy decreases iteratively. When the SNR is equal to or above 15 dB, the final image entropy is lower than 4 which means good image focusing performance.

#### *4.4. Simulations under Different Parameter Ranges*

In order to verify that the proposed method has a certain degree of tolerance for the deviations of the imaging parameters, the following simulations are conducted. According to the first subsection in this section, the center deviation parameter *dx*, *dy* of the imaging plane only cause the image to shift in the imaging plane, and has no influence on the imaging quality and focusing performance. In addition, the posture angle *ξ* makes the image rotate in the imaging plane, while the focusing performance is not affected. So in the following simulations, only *dz*, *δ* and *μ* are considered.

In the simulations, *dz*, *<sup>δ</sup>* and *<sup>μ</sup>* are uniform random selected in [−Δ*dz* <sup>2</sup> , <sup>Δ</sup>*dz* <sup>2</sup> ], [−Δ*<sup>δ</sup>* <sup>2</sup> , <sup>Δ</sup>*<sup>δ</sup>* <sup>2</sup> ] and [−Δ*<sup>μ</sup>* <sup>2</sup> , <sup>Δ</sup>*<sup>μ</sup>* <sup>2</sup> ]. Δ*dz* are set to [0.1, 0.2, 0.3, 0.5, 0.8, 1.0], while Δ*δ* are set to [1, 2, 3, 5, 8, 10]/180*π* whereas Δ*μ* are set to [1, 2, 3, 5, 8, 10]/180*π*. For each Δ*dz*, Δ*δ* and Δ*μ*, 10 Monte Carlo trials are taken. Meanwhile, when the value *dz*, *δ* and *μ* are changing, *dx*, *dy* and *ξ* are set to be zeros.

The errors between the estimated *dz*, *δ*, *μ* and the real values and the image entropy obtained by 10 Monte Carlo trials with different Δ*dz*, Δ*δ*, Δ*μ* are illustrated using boxplot in Figure 11. On each box, the central mark indicates the median, and the bottom and top edges of the box indicate the 25th and 75th percentiles, respectively. The whiskers extend to the most extreme data points not considered outliers, and the outliers are plotted individually using the red '+' symbol. According to Figure 11a–c, the errors between the final estimated *dz* and the real values are within ±0.01 m, and the errors between the final estimated *δ* and the real value are within ±0.35◦ when Δ*δ* is chosen little than 3◦. Moreover the errors between the final estimated *μ* and the real values are within ±0.35◦ when Δ*μ* is chosen little than 5◦. In addition, according to Figure 11d–f, the image entropy are all lower than 4 when Δ*dz*, Δ*δ*, Δ*μ* are chosen within 1 m, 3◦, 5◦ respectively.

**Figure 11.** Errors between the estimated *dz*, *δ*, *μ* and real values and Image entropy of 10 Monte Carlo trials, (**a**–**c**). Errors between the estimated *dz*, *δ*, *μ* and real values respectively when different Δ*dz*, Δ*δ* and Δ*μ* are chosen; (**d**–**f**) image entropy with different Δ*dz*, Δ*δ* and Δ*μ*.

In short, IPCA can obtain the final estimated *dz*, *δ* and *μ* with errors at most ±0.01 m, ±0.35◦ and ±0.35◦ when Δ*dz*, Δ*δ*, Δ*μ* are chosen within 1 m, 3◦, 5◦ respectively, and the image entropy is smaller than 4.

#### **5. Conclusions**

In this paper, the imaging plane mismatch problem of 2D cross-range MIMO radar imaging is analyzed, and the deviations between the estimated spatial spectral point location and the real spatial spectral point location are deduced as well as the corresponding spatial spectral values. To solve this problem, IPCA is proposed in this paper. Aiming at minimizing the image entropy and sparsity of the image, PSO is utilized in IPCA to obtain the image with better focusing performance. Simulation results verify the effectiveness of the proposed algorithm and the robustness to noise. At the same time, when the parameters of the imaging plane are different, the proposed algorithm can obtain the imaging plane parameters of the real values, and then carry out the imaging plane calibration operation to obtain the focused image. The method proposed in this paper can solve the imaging plane mismatch problem and obtain high quality MIMO images.

**Author Contributions:** All authors contributed extensively to the work presented in this paper. Y.G. proposed the original idea and designed the study; B.Y. and Z.W. performed the simulations and wrote the paper; Y.G. supervised the analysis and edited the manuscript; R.X. provided his valuable suggestions to improve this study.

**Funding:** This research was funded by the National Natural Science Foundation of China under contact No. 61771446 and No. 61431016.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Field Representation Microwave Thermography Utilizing Lossy Microwave Design Materials**

**Christoph Baer \*, Kerstin Orend, Birk Hattenhorst and Thomas Musch**

Institute of Electronic Circuits, Ruhr University Bochum, Universitaetsstr. 150, 44801 Bochum, Germany; Kerstin.Orend@rub.de (K.O.); Birk.Hattenhorst@rub.de (B.H.); Thomas.Musch@est.rub.de (T.M.) **\*** Correspondence: christoph.baer@rub.de; Tel.: +49-234-3227606

**Abstract:** In this contribution, we are investigating a technique for the representation of electromagnetic fields by recording their thermal footprints on an indicator material using a thermal camera. Fundamentals regarding the interaction of electromagnetic heating, thermodynamics, and fluid dynamics are derived which allow for a precise design of the field illustration method. The synthesis and description of high-loss dielectric materials is discussed and a technique for a simple estimation of the broadband material's imaginary permittivity part is introduced. Finally, exemplifying investigations, comparing simulations and measurements on the fundamental TE10-mode in an X-band waveguide are presented, which prove the above introduced sensing theory.

**Keywords:** microwave; thermography; field illustration; permittivity

#### **1. Introduction**

The measurement and visualization of electrical and magnetic fields has been an important task in the big scope of instrumentation, measurement, and education for many years as it enables a better understanding of physical phenomena. Here, several differentiations can be made for classifying methods and techniques. In terms of static and quasi-static field observation up to 100 MHz, optical sensors utilizing the Kerr effect, showing a two-dimensional field distribution, have been already reported in the 1970s [1]. Further instrumentation approaches utilize optical sensors [2,3]. However, when increasing the alternating frequency of the applied field up to wave propagation, sensor and probe design is often reduced to single point measurements that record a scalar value. Therefore, the illustration of complete field distributions is only possible by the time consuming scanning of the area of interest. Corresponding sensor and measurement techniques have been reported in [4–6] and include dielectric tips and dipole probes.

Anyway, the fast recording of microwave field distributions remains an interesting approach for several scientific and engineering fields including the investigation of field radiation in terms of antenna design, EMC testing and of course educational purposes. However, large scale and multi-dimensional visualization is challenging, because the desired resolution of the illustrated field image demands for a large amount of sensor array elements, which are used in near field sensor devices. Consequently, when applied to bigger investigation areas, these devices are costly, power consuming, and disturb the incident electromagnetic wave. Hence, far field measurements are often performed for antenna and electromagnetic compatibility (EMC) testing, which however give only information on radiated field components. Finally, nowadays simulation based illustrations are a key technology in antenna design and education. Yet, conformity between virtual and real world results must be verified. Another method for the illustration of electromagnetic fields and their interaction with matter is the microwave thermography (MWT) approach. Here, the thermal footprint of an incident electromagnetic wave on lossy dielectrics or metallic objects is recorded by means of a thermal camera and subsequently processed. While, materials for the illustration of two-dimensional laser beams of non-optical wavelengths,

**Citation:** Baer, C.; Orend, K.; Hattenhorst, B.; Musch, T. Field Representation Microwave Thermography Utilizing Lossy Microwave Design Materials. *Sensors* **2021**, *21*, 4830. https://doi.org/ 10.3390/s21144830

Academic Editors: Andrea Randazzo, Cristina Ponti and Alessandro Fedeli

Received: 18 June 2021 Accepted: 12 July 2021 Published: 15 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

e.g., infrared lasers, are very well known and have been already commercialized and are distributed by Macken Instruments or Cascadelasers, the illustration of microwave fields is challenging because field strenghts' are significantly lower. To date, MWT has been mostly used for non-destructive testing as described in [7]. Further recording techniques rely on thermo-elastic optical indicators [8,9], or solid state quantum sensing [10]. In [11–13], promising approaches for the microwave field illustration are presented that make use of films with magnetic losses in combination with thermo-fluoresecent indicators. However, due to the magnetic loss mechanism, the frequency range of the illustrated field is limited. Moreover, the presented materials demand for high feeding powers up to several watts, which result in comparably low heating and therefore low measurement effects.

Taking previous research results into consideration reveals that a frequency gap between several 100 MHz up to optical frequencies is present for field representation thermography, which we want to partly close in this paper. Next to a holistic modeling of the multi-physical background, we propose novel approaches on the description of lossy dielectric material mixing and further realizations for the accurate representation of electro-magnetic field in the microwave region between 1 GHz and 100 GHz. Obviously, for Field-Representation-MWT (FRMWT), an adequate indicator material must be found, which has to meet several requirements: its parameters regarding dielectric and thermal properties must be chosen wisely in order to increase the overall efficiency and resolution of the thermal footprint. This is desirable as it allows reduced excitation powers and the applicability of low-cost thermal cameras. Moreover, the chosen size of the indicator material plays an important role for its transparency in terms of field interference. Finally, special test benches for FRMWT must be constructed to reduce disturbing environmental influences, such as infrared reflections as well as forced and additional natural convection, which can tamper the field recording.

This manuscript is organized as follows: Section 2 explains necessary fundamentals regarding "wave-matter interaction", "dielectric heating", "thermal loss theory", and "dielectric design materials and mixing". In Section 3, we discuss and set up a novel transient FRMWT model, while Section 4 presents self-made design materials for an optimized FRMWT. Finally, Section 5 shows microwave thermography results in simulation and measurements. In this contribution, the exemplary field distribution under investigation is the fundamental TE10-mode of a rectangular X-band waveguide. Section 6 concludes the manuscript and gives an outlook on future modifications and investigations.

#### **2. Fundamentals**

Since the description of FRMWT demands for the combination of several physical phenomena from electromagnetics, thermodynamics, and flow dynamics, the following section gives a brief overview on the dependencies between the different disciplines.

#### *2.1. Wave Interaction and Dielectric Heating*

The description of the interaction between an electromagnetic wave with the indicator material is very fundamental for FRMWT. In this context we need to observe two properties: the dielectric power dissipation, and the total reflection coefficient. The dielectric power dissipation is essential for the heating process of the indicator material and therefore important for the later field illustration. The total reflection coefficient, however, is of great interest as it gives information on field interferences, which can lead to distorted field illustration. In the simplest case, a plane wave interacts with a dielectric material, which is oriented perpendicular to the Poynting-vector of the incident wave. At the interface between free space and material, reflection coefficients occur that depend on wave impedances of the free space (*Z*0) and inside of the material (*Z*1).

$$
\Gamma = \frac{Z\_1 - Z\_0}{Z\_1 + Z\_0} \tag{1}
$$

At this point, we need to distinguish between the wave's entrance and exit reflection of the material as their coefficients differ. Furthermore, it is obvious that a part of the incident wave's power *P*<sup>i</sup> remains inside of the material and is continuously attenuated by the material as illustrated in Figure 1. The power *P*di, which remains inside of the dielectric material can be calculated as follows:

$$P\_{\rm di} = P\_{\rm l} \cdot \left(1 - \Gamma\_1^2\right) \cdot a(f)^2 \cdot \frac{1}{1 - (\Gamma\_2 a(f))^2}.\tag{2}$$

In (2), the factors Γ<sup>1</sup> and Γ<sup>2</sup> stand for the reflection coefficient at the interfaces freespace-to-material and vice versa. The attenuation *α* depends on the material's thickness *d* and its dielectric properties, which are frequency dependent:

$$\mathfrak{a}(f) = \mathfrak{e}^{-\frac{2uf}{c\_0}\operatorname{Im}\left\{\sqrt{\mathfrak{k}\_\mathbf{r}(f)}\right\} \cdot d}.\tag{3}$$

The total reflection coefficient, which is important for the interference with the incident field, is the square root of the relation of the reflected power *P*refl to the incident power *P*<sup>i</sup> and can be described as the superposition of the reflections at the interfaces free spacematerial and material-free-space. Because of the material's thickness *d*, the contribution of the second reflection performs an additional phase shift, which must be considered for the superposition. By considering the refractive indices *<sup>n</sup>*<sup>0</sup> = 1 and *<sup>n</sup>*<sup>1</sup> = √*ε*r,1*μ*r,1 for free-space and a material, respectively. The total reflection Γtot is calculated by:

$$\Gamma\_{\text{tot}} = \sqrt{\frac{P\_{\text{refl}}}{P\_{\text{i}}}} = \frac{j\left(n\_0^2 - n\_1^2\right)\sin(kd)}{2n\_0n\_1\cos(kd) - j\left(n\_0^2 + n\_1^2\right)\sin(kd)}.\tag{4}$$

In (4), *j* is the imaginary unit, while *k* represents the corresponding wave number. As we consider the dielectric losses of the material, the heating process of the material is directly connected with the magnitude of the electrical field strength |*E*| of the incident wave. By considering the wave impedance *Z* and the length of the electric field vector *l*1, it can be described by:

$$|E| = \sqrt{P\_{\text{di}} \cdot Z} / l\_1 \tag{5}$$

Finally, the dissipated power depends on the imaginary part of the permittivity *ε* <sup>r</sup> [14], the frequency *f* and the material's volume *V*:

$$P\_{\rm loss} = 2\pi f \varepsilon\_0 \cdot \varepsilon\_\mathbf{r}^\prime \cdot |E|^2 \cdot V. \tag{6}$$

The actual heating, which can be regarded as the power transfer from electromagnetic to thermal power, is a transient process and depends on the material's mass *m* and its specific thermal capacity *c*w. It can be described by:

$$
\Delta T = \frac{P\_{\text{loss}}}{c\_W \cdot m} \cdot t. \tag{7}
$$

Equations (5)–(7) show that the heating of an indicator material is directly proportional to the microwave excitation power. However, for calculating the electrical microwave field distribution we need to consider the square root of the recorded thermal footprint.

**Figure 1.** Schematical presentation of the wave interaction with a thin dielectric material sheet.

#### *2.2. Thermal Losses*

As we imbalance the thermodynamic equilibrium by the dielectric heating, the consideration of thermal losses is mandatory at this point. Otherwise, the material's temperature would rise indefinitely, because the heating Equation (7) is a linear function with time. Consequently, different thermal loss mechanisms affect a cooling, which limits the heating process to a final temperature. Here, three main thermodynamic transfer processes must be taken into consideration: thermal conduction, radiation, and convection. In a first approximation, conduction effects can be neglected as we observe pure dielectric materials. Consequently, the main contribution to thermal losses are radiation (rad) and free convection (fc). Therefore, the total thermal power loss can be described by:

$$P\_{\text{th,tot}} = \Delta P\_{\text{rad}} + P\_{\text{fc}}.\tag{8}$$

Regarding radiation losses Δ*P*rad, we only need to consider the additional losses caused by the temperature difference to the surrounding. In good approximation, the phenomenon can be described using the Stefan–Boltzmann law for black bodied radiators, because for the kind of measurements carried out in this research, where temperature differences are measured, the choice of an emissivity is not fundamental. By using the Stefan–Boltzmann constant *<sup>σ</sup>*SB <sup>=</sup> 5.67 · <sup>10</sup>−<sup>8</sup> Wm−2K−4, we receive:

$$
\Delta P\_{\rm rad} = P\_{\rm rad,Tm} - P\_{\rm rad,T\infty} \tag{9}
$$

$$\sigma = \sigma\_{\rm SB} \cdot A\_{\rm rad} \cdot \left( T\_{\rm m}^4 - T\_{\infty}^4 \right). \tag{10}$$

In (9) and (10), *P*rad,Tm and *P*rad,T<sup>∞</sup> represent the thermal radiation powers caused by the temperatures *T*m and *T*∞ on the material and in the environment, respectively. The parameter *A*rad indicates the heated area of the indicator material. Free convection refers to a mechanism of fluid mechanics and is caused by a spatially non-uniform distribution of density. In this case, this non-uniformity is caused by the temperature difference between indicator material and environment [15]. In a first approximation, the transported heat in steady state can be regarded as:

$$P\_{\rm fc}(t \to \infty) = A\_{\rm rad} \cdot \mathfrak{a}\_{\rm th} \cdot \Delta T\_{\rm env}.\tag{11}$$

In (11), *A*rad is the heated area that affects free convection while Δ*T*env = *T*<sup>m</sup> − *T*<sup>∞</sup> describes the temperature difference between the actual indicator material's temperature and the environmental temperature. In case of FRMWT indicator foils, the heated area must be taken twice into consideration as free convection occurs in front and behind of the foil. The parameter *α*th is the heat transfer coefficient. Its determination leads to a problem of flow dynamics and therefore to the determination of the functional relationship between the non-dimensional group of Rayleigh (Ra), Prandtl (Pr), and Nusselt (Nu) numbers [16,17]. Unfortunately, this approach bases on fitted measurement models and the further investigation is not very precise. Anyway, for vertical material orientation, the literature approximates the following equations, which have been constructed from measurement series [18,19]:

$$\mathbf{a\_{th}} = \frac{\mathbf{Nu} \cdot \boldsymbol{\lambda}\_{\text{th}}}{l\_1},\tag{12}$$

$$\text{Nu} = \left( 0.825 + 0.387 \cdot \left( \text{Ra} \cdot \text{Pr}\_{\text{eff}} \right)^{1/6} \right)^2,\tag{13}$$

$$\text{Pr}\_{\text{eff}} = \left( 1 + \left( \frac{0.492}{Pr} \right)^{9/16} \right)^{-16/9} \text{ \AA} \tag{14}$$

$$\mathrm{Ra} = \frac{\mathrm{g} \cdot l\_2^3 \cdot \Delta T\_{\mathrm{env}} \cdot \mathrm{Pr}}{T\_{\infty} \cdot \nu^2}. \tag{15}$$

In (12)–(15), the parameters *λ*th and *ν* represent the surrounding gas' thermal conductivity and viscosity, respectively. Further, *l*<sup>1</sup> and *l*<sup>2</sup> represent geometrical lengths for thermal imbalance alongside the direction of convection as well as the characteristic length of the flow. Finally, *g* is the gravitational acceleration. Thermal losses are highly non-linear regarding time and geometrical setup. Therefore, the presented equations must be regarded as approximation but they facilitate setting up a model for the transient behavior, which is interesting in terms of FRMWT.

#### *2.3. Design Materials and Mixing*

As shown before, the successful FRMWT strictly depends on the right selection of the indicator material. Here, we prefer materials that on the one hand easily heat up and store the additional temperature, while the total reflection coefficient should be kept low in order to minimize field interferences on the other hand. Moreover, thermal conduction should be kept low as well, because it would negatively affect the field illustration's resolution. Consequently, several compromises must be met to serve these requirements. In order to meet as many requirements as possible, an indicator material was designed by mixing several raw ingredients and their corresponding physical parameters. In order to make accurate predictions, so-called mixing equations can be utilized in order to estimate final parameters. In terms of the thermal capacity *c*<sup>w</sup> and mass density *ρ*, the resulting quantity is directly dependent on the volume fractions *ζ* of the corresponding materials. In contrast to this, the description of the resulting, complex valued permittivity for dielectric mixing is not straightforward. However, the manufacturing of phantom materials for the substitution of, e.g., medical substances and surrogate explosives utilizing loss-free materials has been already published in [20,21], respectively. A general overview on dielectric mixing and corresponding models is given in [22]. In terms of FRMWT-materials, the most relevant material parameter are the losses in microwave range. Losses in dielectric materials can be caused by a small remaining conductivity, which leads to small conduction currents, or by replacement currents, which describe dipole orientation etc. Usually, for the general description of a material's complex valued permittivity those effects are tantamount. However, referring to the Wiedemann–Franz law [23], electrical and thermal conductivity are directly proportional so that an increased electrical conductivity will deteriorate the resolution of the FRMWT. Consequently, it is important to tune the material's losses by increasing replacement currents but keeping conductivity as low as possible. Since the electrical polarization strongly depends on the time-variation of the excitation, the permittivity depends on the frequency of those field variations. This effect is called dispersion and can be seen in both real and imaginary part of the permittivity. Part of a material's dispersion is caused by the different polarization effects like orientation polarization as well as atomic and electronic polarization, which cause a resonance behavior in the permittivity's real part associated with a peak in the imaginary part [22]. In any real material, several polarization mechanisms take place and superimpose, which contribute to different dispersion effects within the entire electromagnetic spectrum. Obviously, a holistic model is difficult to formulate. That is why in the literature several polarization response models can be found that describe the aforementioned dispersion and loss effects for limited frequency ranges. Because the proposed FRMWT approach in this manuscript handles the microwave frequency range, we choose the Debye model henceforth. The Debye model is an easy to understand model that considers single relaxations processes very well. As we only handle solid materials in this work, we do not expect various relaxation processes in the microwave region, but continuously decreasing permittivity's real part values. Therefore, the Debye model is a good choice for modeling our indicator materials. Because of the induced torque of the dipole moments, the polarization requires time to reach equilibrium. This relaxation time *τ* is essential for the description of the complex valued permittivity. Moreover, Debye model defines the edges of the observed spectrum. Therefore, *ε*r,<sup>∞</sup> and *ε*r,s describe the permittivities at optical and low frequencies, respectively. The description of a mixed material with Debye-type inclusions in a dispersion and loss-free background with permittivity *ε*r,e can be aligned to the basic Maxwell–Garnett mixing rule. The resulting material is also of Debye-type and can be formulated as follows:

$$
\varepsilon\_{\rm r,eff} = \varepsilon\_{\rm r,\infty,\rm eff} + \frac{\varepsilon\_{\rm r,s,eff} - \varepsilon\_{\rm r,\infty,eff}}{1 + j\omega\tau\_{\rm eff}}.\tag{16}
$$

Consequently, the parameters of the mixed material are:

$$
\varepsilon\_{\rm r,\infty,eff} = \varepsilon\_{\rm r,e} + \frac{3\zeta\varepsilon\_{\rm r,e}(\varepsilon\_{\rm r,\infty} - \varepsilon\_{\rm r,e})}{\varepsilon\_{\rm r,\infty} + 2\varepsilon\_{\rm r,e} - \zeta(\varepsilon\_{\rm r,\infty} - \varepsilon\_{\rm r,e})},\tag{17}
$$

$$
\varepsilon\_{\rm r,s,eff} = \varepsilon\_{\rm r,c} + \frac{3\zeta\varepsilon\_{\rm r,c}(\varepsilon\_{\rm r,s} - \varepsilon\_{\rm r,c})}{\varepsilon\_{\rm r,s} + 2\varepsilon\_{\rm r,c} - \zeta(\varepsilon\_{\rm r,s} - \varepsilon\_{\rm r,c})}\,,\tag{18}
$$

$$
\pi\_{\rm eff} = \pi \frac{(1 - \zeta)\varepsilon\_{\rm r,\infty} + (2 + \zeta)\varepsilon\_{\rm e}}{(1 - \zeta)\varepsilon\_{\rm r,s} + (2 + \zeta)\varepsilon\_{\rm e}}.\tag{19}
$$

Next to the Maxwell–Garnett averaging of the static and optical permittivity in (17) and (18), respectively, (19) reveals that the relaxation time of the including particles is also altered [22]. This is of great interest for the FRMWT as the material mixing allows for an optimization of the losses for a predefined operation frequency.

Since the measurement of lossy materials results in extremely low signal-to-noise ratios, especially the determination of the the loss-describing part, i.e., the imaginary part of the permittivity, can be crucial. Moreover, measurement probes are prone to delivering wrong loss values when the probe-material connection is not properly realized and generates radiation. However, obviously the connection between frequency behavior of the real and imaginary part can be formulated. Starting from the basic physical principle of causality, which means that the polarization response cannot lead the cause, the Kramers–Kronig relation must be fulfilled for the real and imaginary part of the permittivity [24]. As conse-

quence, this very fundamental connection allows for an estimation of the imaginary part when the broadband real part permittivity is known. The estimation itself can be done in various ways, such as solving the Kramers–Kronig integrals, applying Hilbert transformation [25], or simply using the complex valued completion (CVC) which is described further on. The CVC allows for the estimation of the permittivity's imaginary part, when the corresponding permittivity's real part is known. It makes use of the dependencies of complex valued time signals and their spectra, which are subdivided in their even and odd signal and spectral parts. Figure 2 illustrates the well known dependencies. Within the figure, the sign ❝ indicates the Fourier transformation from time (◦) to frequency domain (•). Moreover, the indices E and O represent even and odd signal and spectral parts, while the superscripts R and I indicate real and imaginary parts. As initial point, we consider a measured broadband, real part permittivity. As this quantity depends on frequency, we can regard it as a spectrum. Moreover, any commercial measuring probe will deliver measurement values for positive frequencies only. Therefore, the measured permittivity is the combination of its even and odd part:

$$
\varepsilon\_{\rm r}^{\rm K} = 1/2(\varepsilon\_{\rm r,E}^{\rm K} + \varepsilon\_{\rm r,O}^{\rm K}).\tag{20}
$$

By means of the following CVC-steps, illustrated in Figure 3, we can calculate the right-sided, complex permittivity. As initial step we form the evenly-shaped spectrum of the measured real part permittivity. Therefore, the spectrum *ε<sup>R</sup>* <sup>r</sup> is multiplied by the heaviside step function Θ and divided by 2. Additionally, the resulting data vector is mirrored to the negative frequency range in order to achieve an axis-symmetric spectrum, which is described as the evenly-shaped spectral part of the real part permittivity *ε<sup>R</sup>* r,E. According to Figure 2, the application of the inverse Fourier transformation delivers the evenly-shaped impulse response of the real part permittivity *E<sup>R</sup>* r,E(*t*). In [22], several interpretations of this pulse response signal are given, however in our case the time domain signal is a pure mathematical intermediate step and further interpretation is not relevant. As the next step, we transform the evenly-shaped impulse response to be oddly-shaped. Therefore, we once more apply the heaviside function and add the point symmetric signal for *t* < 0. By applying the Fourier transformation to the oddly-shaped time signal, we receive the oddly-shaped spectrum of the permittivity's imaginary part *εI* r,O(*ω*). By doubling this spectrum, applying another heaviside step function, and adding the right-sided spectrum of the real part permittivity, we finally receive the right sided complex valued permittivity *ε*<sup>+</sup> <sup>r</sup> (*ω*).

**Figure 2.** Dependencies of complex-valued time signal and the corresponding spectrum sub-divided into even (E) and odd (O) valued functions. The handle bar shows the intercorrespondencies.

It needs to be mentioned that the CVC is only applicable if the measured permittivity data is sufficiently broadband and covers all relevant polarization behaviors. In terms of FRMWT and the corresponding design material manufacturing, CVC plays an important role for the precise description of the high loss indicator materials as well as their corresponding mixing equations.

**Figure 3.** Processing steps of the complex valued completion (CVC). While rectangular boxes describe signal operation, circles indicate the current signal form. The applied function Θ indicates the heaviside function in frequency-domain (*ω*) and time-domain (*t*).

#### **3. Transient FRMWT Model**

In the following, we will observe the special case of FRMWT for mode characterization of a rectangular waveguide in X-band, as illustrated in Figure 4. The corresponding geometrical dimensions for the investigated WR90 rectangular waveguide and IEC 60154- 2:2016 flange are: *l*<sup>1</sup> = 10.16 mm, *l*<sup>2</sup> = 41.4 mm, and *w* = 22.86 mm. Moreover, the indicator material's thickness is considered to be *d* = 0.9 mm and it is positioned in close proximity to the waveguide's flange.

**Figure 4.** Sketch of the investigated waveguide flange.

However, it does not touch its metal boundary in order to prevent additional thermal conduction. After reaching steady state, i.e., when the microwave power which is transferred to heat equals the thermal losses, the heating-up process stops and the indicator material's end-temperature *T*end is reached. Consequently, by equating the dissipated microwave power with the thermal loss power, we can calculate *T*end and therefore the final temperature change Δ*T*end. Figure 5 combines the thermal heating model with the microwave scenario by showing the theoretical temperature change Δ*T*end as a function of the thermal loss power for the described setup and different surrounding gases. The determination of the temperature change is of great interest as it has direct influence on the sensitivity of the thermal sensor and the applied microwave power. Therefore, for a given thermal sensor sensitivity, the required microwave power can be reduced, which is beneficial in many scenarios. Apparently, a higher temperature change is reached when choosing surrounding gases like Helium or Xenon. However, the highest temperature change is achieved when the setup is kept in vacuum, because free convection vanishes and only thermal radiation contributes to the thermal loss. Anyway, this insight shows that the construction of an optimized FRMWT test stand that guarantees for a highly defined environment scenario is desirable. Yet, low cost setups also yield sufficient results in many applications.

Since the heating process of the indicator material can be interpreted as charging process, the transient description of microwave induced heating in FRMWT can be expressed as:

$$
\Delta T(t) = \Delta T\_{\rm end} \cdot \left(1 - \mathbf{e}^{-\frac{\beta\_{\rm th}}{\Delta T\_{\rm end}} \cdot \frac{P\_{\rm lens}}{r\_{\rm W} \cdot m} \cdot t} \right) = \Delta T\_{\rm end} \cdot \left(1 - \mathbf{e}^{-\frac{t}{\tau\_{\rm th}}} \right). \tag{21}
$$

In (21), the parameter *β*th is called thermal retardation, a form factor that compensates measurement setup geometries etc. Once the thermal retardation is evaluated the time constant *τ*th can be obtained. It combines all the previous theoretical considerations from electromagnetics in Section 2.1 and thermodynamics in Section 2.2. This time constant is of great interest, because it gives indication for the compromise that needs to be found for FRMWT-measurements: On the one hand, the induced temperature differences of the indicator should be as high as possible in order to meet the thermal resolution of the utilized recording device, such as a thermal camera. On the other hand, the footprint recording should be performed as fast as possible after initial excitation in order to prevent image blurring and further non-linear behavior such as conduction effects. Therefore, we propose to take the thermal snapshot for field illustration purposes within the linear part of the transient model described in (21), i.e., for

$$
\tau\_{\rm th} < t\_{\rm rec} < \Im \cdot \tau\_{\rm th}.\tag{22}
$$

So that the heating process already reached between 68% and 95% of the final temperature.

**Figure 5.** Theoretical temperature change for maximum thermal losses in steady state and different gaseous environments.

#### **4. Design Materials**

For the manufacturing of adequate indicator material foils, several material mixtures have been investigated. In order to meet the materials requirements, i.e., providing high dielectric losses by keeping the conductivity low and enabling the possibility of forming robust foils, we chose epoxy resin as matrix material and carbon black as additive. The utilized epoxy resin Breddepox E300 from Breddemann is a low viscosity, two component casting resin, providing a high UV resistance. As carbon black, we used the DUREX 0-Powder from manufacturer Orioncarbons. This amorphous carbon black provides low conductivity, offers very good processability, exhibits low compression set, and has good dynamic properties by keeping high elasticity. For tuning and fitting the mixing equation, we produced several test mixtures with different carbon black volume fractions *ζ*cb between 0.01 and 0.09, which than were characterized by means of a DAK1.2E dielectric probe kit from SPEAG. As mentioned earlier, the measurement of high loss materials can lead to inaccuracies, especially for the permittivity's imaginary part. Figure 6a–c show the broadband measured real and imaginary part as well as the loss tangent, respectively, for a volume fraction of *ζ*CB = 0.07. In (a) we added a fitting curve, which is the result of a least-square-fit of the Debye models permittivity from (16). While in (a) this curve perfectly balances the measured values, it absolutely mismatches the measured imaginary part in (b). Therefore, we added CVC-processed curves from both measurement and fitting data. Obviously, as the fitting curves and CVC-processed curve are nearly identical, the directly measured imaginary part is highly erroneous and will be ignored further on.

**Figure 6.** Exemplary results of the mixed material investigation at an carbon black volume fraction of *ζ*cb = 0.07 , indicating the real (**a**) and imaginary part (**b**) permittivity as well as the loss tangent (**c**) for (—) measured values, (—) Debye model fitted values, (—) CVC of measured values and (—) CVC of Debye model fitted values.

Consequently, we used the direct-fitted data for setting up the mixing equations. Figure 7 shows several 3D plots indicating the properties of processable indicator materials such as the complex permittivity (a) + (b), total reflection coefficient (c), and the expected, maximum heating (d) for the scenario described above. In the following investigation, we are using an indicator material foil with a volume fraction of *ζ*CB = 0.05, which allows for homogeneous mixing and a temperature increase of more than 2 ◦C at 10 GHz when applying 500 mW microwave power.

**Figure 7.** Resulting parameters of the fitted material mixing equation showing (**a**) real part permittivity, (**b**) imaginary permittivity, (**c**) total reflection coefficient, (**d**) maximum heating. Subfigures (**c**,**d**) apply for the above described indicator material and measurement scenario.

#### **5. Results**

In order to investigate the applicability of FRMWT and to validate the presented models, several simulations and measurements have been performed on a simple test scenario. The chosen scenario is the fundamental TE10-mode within a WR90 rectangular waveguide utilizing a IEC 60154-2:2016 flange as introduced in Figure 4. The indicator material was selected from the considerations made in Section 4 and provides a complex permittivity of *ε*<sup>r</sup> = 7.2 + j1.6; as operating frequency we chose 10 GHz and an excitation power of 500 mW in both, simulations and measurements. As the waveguide will operate in fundamental mode, we expect a *cos*-shaped E-field distribution along side the x-axis of the waveguide. However, since the square of the electrical field strength contributes the thermal heating as shown in (6) we expect a cos2-shaped heat distribution. It needs to be mentioned that because the waveguide is truncated it will excite evanescent modes as well, which will slightly interfere with the fundamental mode resulting in minimal distortion of the thermal foot print.

#### *5.1. Measurement Setup*

From theoretical approximations in Section 2.2 and 3, we can conclude that the observable temperature increase for FRMWT applications will be only several kelvin, when applying reasonable microwave power. Therefore, we need to build a measurement environment, which prevents thermal and environmental disturbances. In this work, the thermal footprint of the device under test will be recorded by means of a high resolution thermal camera. Usually, these kind of cameras are portable so that it can be integrated into a measurement chamber. Further key properties of this chamber are:

• Housing:

The measurement chamber needs to be closed and air tight for enabling different measurement environments such as gases. Moreover, it must prevent forced convection.

• Utilized Materials: The housing should be impervious for infrared radiation for preventing disturbing environmental reflections. The metal content should be kept low or needs a microwave absorber covering in order to avoid disturbing reflections.

• Peripheral Connections: Gas in- and outlets allow for FRMWT in different surrounding gases. Microwave connectors are necessary for feeding the structure under test.

The constructed measurement chamber was built of Plexiglas plates as these sufficiently attenuate infrared propagation within the spectrum of interest between 3 and 15 μm. The chamber measures a total size of (height, width, length) 50 cm × 55 cm × 45 cm. Gas inand outlets can be used for flooding the whole setup with gases, while an SMA-throughhole connector allows for the connection of microwave devices within a frequency range of 0.1–30 GHz. Furthermore, an optical breadboard as lower wall allows for an universal connection of mechanical components. Figure 8 shows a photography of the described measurement chamber. As thermal camera we used a VarioCAM HD research 900 from Infratec. This thermal camera records a spectral range from 7.5 to 14 μm by means of an uncooled microbolometer focal plane array, which results in an image resolution of 2048 × 1536 pixels and a thermal resolution of 20 mK. Because the indicator material was manufactured from epoxy resin as matrix material we used an emissivity of 0.96 and a spectral range of 8 μm to 14 μm for all recordings, which refers to the values recommended by the thermal camera manufacturer [26]. The influence of further parameters is neglected in this work, because we mainly focus on temperature differences and its evolution over time.

**Figure 8.** Photography of the realized setup containing: (**a**) Thermal camera, (**b**) through-hole microwave connection, (**c**) microwave device under test, (**d**) universal mechanical mounting board, (**e**) rear absorber wall.

#### *5.2. Simulation and Measurement Results*

In order to investigate the applicability of the indicator material the FRMWT-concept, multi-physics simulations as well as measurements with a high definition thermal camera were performed and compared. For all simulations we used the simulation software CST Microwave Studio 2020 , which allows for the coupling of three-dimensional full wave analysis with a thermal simulation of the corresponding loss distributions. In terms of the microwave investigation the time domain solver with a simulation accuracy of −30 dB was utilized. In order to ensure accurate simulation results, a hexahedral mesh type was used with a mesh of approximately 1.2 million meshcells. In a subsequent step, thermal simulations were carried out, basing on the microwave simulation results in order to track the heating behavior over time and to illustrate the temperature distribution caused by the microwave losses.

For this purpose, the convective heat transfer coefficient *α*th was set to 15.3 Wm−<sup>2</sup> K−1. Due to the gap between material and waveguide, convection was assumed on both sides. Thermal conduction was not considered further due to the experimental setup. In addition, the material parameters of the indicator foil were included. These contain the previously measured permittivity *ε*<sup>r</sup> = 7.2 + j1.6. Further material parameters were calculated from the raw materials' datasheets in order to agree with the measurement setup. Here, we set the heat capacity to be *c*<sup>w</sup> = 894 J (kg K)−1, the mass density to be *ρ* = 1100 kg m−<sup>3</sup> and the material thickness to *d* = 0.9 mm. The remaining geometrical dimensions apply as described above. All simulations were carried out with an microwave excitation power of 500 mW. The following thermography images were captured on a truncated WR90 rectangular waveguide surrounded by an IEC 60154-2:2016 flange. The test recording was performed at an operating frequency of 10 GHz and a power of 500 mW, for keeping the measurement results comparable to the simulation. For the verification of the proposed transient model from Equation (21), we recorded the setup's hotspot temperature over time in both simulation and measurement. By subtracting the initial temperature, the temperature increase over time is obtained, which is presented in Figure 9 and compared to the transient FRMWT model from (21) for different thermal retardation values *β*th. The steady state temperature difference was determined to be 2.5 K and 2.71 K in simulation and measurement, respectively. The discrepancy between simulation and measurement is 0.21 K, which is an acceptable, relative difference of 7.7% and can be explained by a different implementation of the thermal losses within the simulation tool. Because of the slightly different final temperature, the transient behavior also slightly differs between simulation and measurement as presented in Figure 9. However, the measurement fits very good to the transient FRMWT model from (21) when choosing the thermal retardation to *β*th = 0.35. After the determination of all relevant parameters, the thermal time constant is calculated to be *τ*th = 12 s. Consequently, the following thermal images are recorded at *t*rec = 35 s, which perfectly satisfies Equation (22). The recorded heat distributions, resulting from simulation and measurement, are shown in Figure 10a,b, respectively. In both images, the hotspot is located in the waveguide's center position as expected because of the E-field maximum of the fundamental mode. Towards the edges of the material the heating decreases as expected. However, in y-direction the expected rectangular shape is flattened due to a remaining thermal conductivity and fringe fields in the flange-material-interspace.

**Figure 9.** Transient temperature observation of the investigated waveguide's center, indicating the FRMWT model for different model parameters (solids), simulation results (– –) and measurement results (∗).

**Figure 10.** Thermography images of the investigated waveguide from (**a**) simulation and (**b**) thermal camera.

Figure 11 shows the temperature distributions alongside the waveguide's x-axis in comparison to a normalized cos2 function. Both, simulation and measurement results, clearly reveal the expected field distribution. However, the simulation result fit the cos2 reference a little better. An explanation for this could be the combination of evanescent fields at the truncated waveguide flange as well as the neglected thermal conduction in simulation.

**Figure 11.** Temperature distribution of the heating in the waveguide's center alongside the x-axis for simulation (– –) and measurements (–) compared to the theoretical *cos*2-shape (–).

#### **6. Discussion and Conclusions**

In this contribution, we present a method called Field Representation Microwave Thermography (FRMWT) for the fast and accurate recording of electromagnetic fields in the microwave range. Theoretical background from microwave engineering, thermodynamics, and fluid mechanics has been derived in order set up a holistic theory, which allows the prediction of the desired field illustration. Regarding the necessary indicator material, the so called complex valued completion (CVC) has been introduced, which allows for the measurement validation of high loss materials. Investigations on mixed high-loss materials showed the necessity of the CVC for dielectrics, as straight forward measurements delivered erroneous values for the material permittivity's imaginary part. Further, within simulations and practical experiments, we demonstrated the applicability of FRMWT using the example of the fundamental TE10-mode of rectangular waveguides. Within these investigations, we proved both, the theoretical transient FRMWT model as well as the ability of the field illustration by recording the predicted *cos*2-shape of the TE10 mode's thermal foot print. Compared to previous work, which already showed the applicability of field representation in low frequency ranges up to 100 MHz and in the optical region, this work accesses the gap-region of micro- and mmWaves from 1 GHz up to 100 GHz. Although the proposed theoretical FRMWT model makes use of several approximations, it provides an excellent and holistic foundation for future measurement setups. The introduced lossy material modeling and complex valued completion combines previous work on dielectric material mixing and gives easy access to future material modeling. Of course, its accuracy is limited to the preciseness of a priori material-knowledge and the accuracy of the necessary broad material characterization. The presented measurement setup for the FRMWT already showed a good functionality. However, as we still make use of an extremely expensive thermal camera, the current setup is only useful for laboratory measurements. In order to establish the FRMWT in practical fields, the thermal camera must be substituted. In future work, further investigations on more complex propagation modes at higher frequencies will be performed and a post-processing concept will be evaluated, which delivers key facts on the investigated mode such as the mode purity.

**Author Contributions:** Conceptualization, C.B.; methodology, C.B. and T.M.; software, K.O.; validation, B.H. and C.B.; formal analysis, C.B., K.O. and B.H.; investigation, C.B.; resources, T.M.; data curation, C.B.; writing—original draft preparation, C.B.; writing—review and editing, C.B.; visualization, C.B.; supervision, C.B. and T.M.; project administration, C.B. and T.M.; funding acquisition, C.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by MERCUR Research Center within the TORONTO project (AN-2019-0026).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

FRMWT Field Representation Microwave Thermography CVC Complex Valued Completion

#### **References**


## *Article* **Model-Based Systems Engineering Applied to Trade-Off Analysis of Wireless Power Transfer Technologies for Implanted Biomedical Microdevices**

**Juan A. Martínez Rojas 1,\*, José L. Fernández 2, Rocío Sánchez Montero 1, Pablo Luis López Espí <sup>1</sup> and Efren Diez-Jimenez <sup>1</sup>**


**Abstract:** Decision-making is an important part of human life and particularly in any engineering process related to a complex product. New sensors and actuators based on MEMS technologies are increasingly complex and quickly evolving into products. New biomedical implanted devices may benefit from system engineering approaches, previously reserved to very large projects, and it is expected that this need will increase in the future. Here, we propose the application of Model Based Systems Engineering (MBSE) to systematize and optimize the trade-off analysis process. The criteria, their utility functions and the weighting factors are applied in a systematic way for the selection of the best alternative. Combining trade-off with MBSE allow us to identify the more suitable technology to be implemented to transfer energy to an implanted biomedical micro device.

**Keywords:** trade-off analysis; medical MEMS; wireless power transfer

#### **1. Introduction**

New sensors and actuators based on MEMS technologies are increasingly complex and quickly applied to new products. MEMS-based devices demand system engineering approaches, which were previously limited to very large projects, and it is expected that this need will increase in the future, particularly in biomedical applications.

During product research and development, one important step is the powering of the implantable microdevices [1]. For this choice, there are many different technologies that should be analyzed in order to properly select the best option for each product. This trade-off is normally carried out based on previous experience and/or using a classical weighted approach. However, when the trade-off is started from current literature reviews, the amount of data is so large that it becomes necessary to develop more complex and systematic tools for reaching the optimum selection for microdevice applications.

In this work, we propose the application of Model Based Systems Engineering (MBSE) to systematize and optimize the trade-off analysis process. The criteria, their utility functions and the weighting factors are applied in a systematic way for the selection of the best alternative. Combining trade-off studies with MBSE allows us to identify the more suitable technology to be implemented to transfer energy to an implanted biomedical MEMS device. At present, a fully automated algorithm able to perform the complete decision process in the design of a new device is unfeasible. MBSE and trade-off analysis involve a wise combination of science, engineering and art. However, this cognitive endeavor can be made more rigorous by applying certain rules and mathematical techniques. In this work, we show how such techniques are applied to a very difficult challenge, which tries to push the limits of the MEMS engineering to design an implantable medical device with

**Citation:** Rojas, J.A.M.; Fernández, J.L.; Sánchez Montero, R.; Espí, P.L.L.; Diez-Jimenez, E. Model-Based Systems Engineering Applied to Trade-Off Analysis of Wireless Power Transfer Technologies for Implanted Biomedical Microdevices. *Sensors* **2021**, *21*, 3201. https://doi.org/ 10.3390/s21093201

Academic Editors: Andrea Randazzo, Cristina Ponti and Alessandro Fedeli

Received: 25 March 2021 Accepted: 1 May 2021 Published: 5 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

size in the order of 1 mm3. Logically, the selected criteria will be severely constrained by this demand.

In Section 2, the systems engineering methodology used in this study is briefly explained. It is a very general, comprehensive and powerful approach, but in this case, we show its application to the trade-off analysis of the best present Wireless Power Transfer (WPT) alternatives for the intended biomedical implantable microdevice.

In Section 3, the state-of-the-art WPT for medical applications is reviewed. In Section 4, we detail the different steps which permit a complete and accurate assessment of different engineering alternatives. This trade-off analysis is described in the context of a larger systems engineering design effort, which is explained in Section 5.

In Section 6, we comment on the results of the trade-off analysis and its implications. The information obtained in this study permits a logical selection of the best WPT system for our application and provides valuable clues to other kind of important studies, such as sensitivity analysis.

#### **2. Trade-Off and Heuristics in the MBSE Methodology ISE&PPOOA**

At present, there are several MBSE methodologies, more or less complete or general in scope [2]. The ISE&PPOOA methodology is a logically consistent approach to MBSE which combines the best features of the traditional and modern tools for optimal design; for example, N2 charts with SysML diagrams. ISE&PPOOA proposes two views of the architecture of the modeled system [3]. These views are the functional architecture and the physical architecture. Below, we briefly describe each of them.

The functional architecture uses diverse SysML diagrams [4] and tables. We understand a function as a transformation to be performed by the system that consumes mass energy or data and generates new ones or transforms them.

Combining diagrams and tables, the functional architecture represents the functional hierarchy using a SysML block definition diagram. This diagram is complemented with activity diagrams for the main system functional flows to represent the system behavior. The N2 chart is a table used for an interface description where the main functional interfaces are identified. A textual description of the system functions is provided, as well.

The physical architecture, or architecture of the solution, is developed in two main steps, the results of which are the so-called modular architecture and refined architecture. The modular architecture represents a first version of the solution architecture representing its main logical blocks. These logical block or modules are blocks allocating functions based on the principles of maximum cohesion and minimum coupling between them.

The transition from the modular architecture to the refined architecture is where we apply a combination of trade-off analysis and design heuristics. Although the scope of systems engineering application of trade-off studies is wider, trade-off analysis is prescribed here for choosing and ranking alternative solutions to be applicable to the system component level. Instead of trade-off analysis we recommend the use of design heuristics [3] to implement those non-functional requirements that apply to the architecture design so that they are implemented as design patterns or at the level of system connectors.

Solution architecture is represented by the system decomposition into subsystems and parts using a SysML block definition diagram. This diagram is complemented with SysML internal block diagrams representing the system physical blocks with either logical or physical connectors for each identified subsystem, and activity and state diagrams for behavioral description as needed. A tabular description of the system parts may be provided as well. Functional allocation may be represented either in tabular form or at the system blocks, allocation by definition, or as partitions in the activity diagrams, allocation by usage, represented using SysML notation.

We get benefit of this systematic approach for system engineering to develop an optimal trade-off selection method. In this case, we apply this methodology to WPT technologies for the specific case of an implantable MEMS device. The benefits of integrating the design process for a suitable medical WPT in the range of a 1 mm<sup>3</sup> into a complete MBSE

design framework are manifold, due to its complexity, which demands the exploration of multiple alternatives that have not been already tested for such a small size.

#### **3. State-of-the-Art Wireless Power Transfer Technology for Biomedical Applications**

Many patients can take advantage of implanted biomedical devices. New microelectromechanical systems (MEMS) technologies allow the development of increasingly complex and miniaturized devices. This reduces many risks associated with the treatment, from surgery to tissue compatibility, due to reduced invasiveness. Another very important factor for implanted devices is energy consumption. Permanent wires inside the body are an unacceptable solution in most cases, while batteries can be dangerous and must be replaced using surgery.

The potential of MEMS for low energy consumption opens the possibility of remote powering without external contacts. The technologies that permit remote powering are collectively known as Wireless Power Transfer (WPT). These technologies are relatively recent and are under development. Although there are many publications exploring detailed aspects of some technologies, there are relatively few reviews comparing different alternative approaches in a systematic and weighted way.

A very recent and complete review was written by Khan et al. in 2020 [5], who compared several WPT technologies using the following parameters: implant type, implant WPT system size, distance from power source to device, type of radiation and frequency, input power, efficiency, test model, SAR (specific absorption rate), safety considerations and technological maturity. The different WPT approaches studied in this article were: NRCC (Non-Radiative Capacitive Coupling), NRIC (Non-Radiative Inductive Coupling), NRMRC (Non-Radiative Magnetic Resonance Coupling), NRRMF (Non-Radiative and Radiative Mid-Field), RFF (Radiative Far-Field), APT (Acoustic Power Transfer) and OPT (Optical Power Transfer). The general conclusion of this review, after a qualitative study of the performance of these technologies, using a Low–Medium–High scale, was that NRIC and NRMRC were better than other WPT techniques due to their moderate size, range and higher PTE performance. Additionally, more complete studies on tissue safety exist for NRIC and NRMRC WPT. APT was comparable to NRIC and NRMRC WPT in terms of performance; the other technologies (NRCC, NRRMF, RFF and OPT) were still too immature from a technological point of view.

Another recent review was provided by Zhou et al. [6], who compared several WPT systems, investigating their key performances such as power transfer capability, power level, efficiency, safety requirements and some others. They organized their review by medical applications, instead of energy types or technologies, but their study was restricted to electromagnetic near-field WPT devices.

Moore et al., in 2019, reviewed the state-of-the-art electromagnetic WPT, especially magnetic resonance, in medicine [7]. They found 17 relevant journal papers and/or conference papers and separated them into defined categories: Implants, Pumps, Ultrasound Imaging and Gastrointestinal (GI) Endoscopy. They found no strong correlation between the system parameters and biomedical applications.

Mahmood et al., in 2019, proposed a new ultrasound sensor-based WPT for lowpower medical devices [8]. A 40 kHz ultrasound transducer was used to supply power to a wearable heart rate sensor for medical application. The system consisted of a power unit and a heart rate measurement unit. The power unit included an ultrasonic transmitter and receiver, rectifier, boost converter and super-capacitors. At 4 F, the system achieved 69.4% transfer efficiency and 0.318 mW power at 4 cm. They also compared their work with previous ultrasound WPT systems. They remarked that power and efficiency decreased as the air gap increased to more than 4 cm.

Kakkar in 2018 designed an ultra-low power system architecture for implantable medical devices based on an embedded processor platform chip [9]. This work explored the partitioning of a chip, as well as the trade-offs associated with design choices, especially intelligent power management. He described several design requirements implying

ultra-small volume, improved resolution for both sensing and stimulation, integrated electronics design, autonomous operation without batteries and intelligent power management. Concrete values for power operation were provided in figures, being in the order of 250 microwatts.

Shadid and Noghanian, in 2018, wrote a literature survey on WPT for biomedical devices based on inductive coupling, concentrating on the applications using near-field power transfer methods [10]. They compared different systems with the following parameters: frequency, output power, transmitter dimensions, receiver dimensions, gap and efficiency. They plotted efficiency versus delivered power and efficiency versus frequency for the reviewed systems.

Taalla et al., in 2018, presented a comparison of inductive and ultrasonic WPT techniques used to power implantable devices [11]. The inductive and ultrasonic techniques were analyzed studying their sizes, operating distance, power transfer efficiency, output power and overall system efficiency standpoints. They concluded that the inductive coupling approach can deliver more power with higher efficiency compared to the ultrasonic technique, but the ultrasonic technique can transmit power to longer distances.

Agarwal et al., in 2017, presented a comparison of various power transfer methods based on their power budgets and WPT range [12]. Power requirements of specific implants such as cochlear, retinal, cortical and peripheral were also considered. Patient's safety concerns with respect to electrical, biological, physical, electromagnetic interference and cyber security were also explored. Their conclusion was that EM, NRIC and NCC were better for high-power devices, but for low-power in the order of fewer milliwatts, ultrasonic, mid-field, or far-field technologies are promising.

Dinis et al., in 2017, presented a review of the state-of-the-art implantable electronic devices with wireless power capabilities, ranging from inductive coupling to ultrasounds [13]. They compared the different power transmission mechanisms and showed that the power that current technologies can safely transmit to an implant is reaching its limit. In order to overcome these difficulties, they proposed a new approach, capable of multiplying the available power inside a brain phantom for the same specific absorption rate (SAR) value. They compared previous devices using WPT link distance, antenna/transducer size, received power at the implant, link efficiency and the calculated power density (obtained by dividing the received power by the antenna/transducer size), and ultrasound seemed to be the best WPT solution. Moreover, biological energy harvesters for implantable devices using biologically renewable energy sources, such as muscle movement, vibrations or glucose, were briefly discussed.

The review of Kim et al., published in 2017, focused on Near-Field Wireless Power and Communication for biomedical applications [14]. They proposed that near-field magnetic wireless systems had advantages in water-rich environments, such as biological tissues, due to lower power absorption. However, various issues in near-field magnetic systems remained, such as transmission range, misalignment and limited channel capacity for communications. They suggested that mid-field coupling based wireless powering was convenient for smaller-sized implants using the sub-GHz range. Finally, they indicated that the Q-factor of the coils and their cross coupling are the primary factors that need to be taken into account for system performance optimization.

Lu and Ma (2016) made a review of the best architectures for efficient WPT, allowing further device miniaturization and higher power loss reduction [15]. The main contribution of this article was its best design guidelines for WPT systems based on near-field inductive coupling wireless power transfer. They remarked that operating at a higher WPT frequency can lead to significant size reduction of passive components and better Q values with smaller inductance. However, higher frequencies produced higher tissue absorption inside the human body. They explained that the selection of a correct WPT system architecture for portable or implantable biomedical applications was a trade-off between device volume, efficiency, regulation accuracy, speed and functionality.

Altawy and Youssef, in 2016, studied the trade-off between security, safety and availability in implantable medical devices [16]. They discussed the challenges and constraints

associated with securing such systems and focused on the tradeoff between security measures required for blocking unauthorized access to the device and the safety of the patient in emergency situations where such measures must be dropped to allow access. They analyzed the up-to-date proposed solutions and discussed their strengths and limitations.

Safety and thermal aspects of implanted medical devices were described in [17] by Campi et al. in 2016. They studied a WPT system based on magnetic resonant coupling applied to a pacemaker for recharging its battery using a low operational frequency (20 kHz). Other safety considerations for an implantable rectenna for far-field WPT can be found in [18]. Other very informative specific monographs dedicated to WPT for biomedical devices can be found in [19–22].

All these reviews and monographs articles give exhaustive and detailed information about current trends in WPT. However, even with this information, so many alternatives make it difficult to properly select the most adequate solution for a specific application. Here, we propose a rigorous engineering approach for a trade-off study based on systems engineering and particularly using the models of the system developed by the ISE&PPOOA MBSE methodology. This new comprehensive approach permits conclusion with the most convenient choice and it also allows the ranking of the rest of the solutions for an eventual selection.

#### **4. Steps of Trade-Off Studies in the ISE&PPOOA MBSE Context**

Based on diverse processes for trade-off analysis found in the literature, we use here a trade-off analysis process that can be integrated with the ISE&PPOOA architecting design using its outputs and producing inputs to the ISE&PPOOA MBSE process presented in Chapter 4 of the ISE&PPOOA book [3]. Traditional approaches such as those found in the NASA report [23] in 1994 do not use SysML system models. However, recent approaches such as IBM use SysML notation and diagrams [24]. The steps of the proposed trade-off subprocess of ISE&PPOOA are presented in Figure 1.

**Figure 1.** Trade-off and heuristics subprocess for obtaining the refined physical architecture.

#### *4.1. Identify the System Modules for the Trade-Off Study*

From the modular architecture obtained in step 4.2 of the ISE&PPOOA process described elsewhere [3], modules (blocks) are selected, which are the logical building elements clustering cohesive functionality. These may be implemented using alternative technical solutions that are identified in the next step.

#### *4.2. Identify Credible Alternative Technical Candidates for Implementing the System Logical Building Blocks or Modules under Consideration*

The list of technical alternatives selected during brainstorming sessions may be reduced if the system requirements are considered to have been met [25]. Some alternatives may be discarded based either on cost, technology readiness or other criteria used to eliminate alternatives. The remaining alternatives that need to be assessed should be described in detail.

#### *4.3. Define Trade-Off Criteria*

Objectives related to stakeholder needs and system requirements are transformed into a set of performance, cost and other criteria to be used for the trade-off study.

#### *4.4. Assign Relative Weightings to the Criteria*

There are diverse methods for deriving numerical values to the weights to be assigned to the criteria. If a pair-wise comparison is used, we recommend the Analytical Hierarchical Process (AHP) [23] to establish the relative weights to the criteria at the same level, so that all weights sum to 1.0. Another method we recommend and use in the present WPT case is the "swing weight matrix" [26]. Swing weights are assigned to the criteria based not only on importance but on the variation of their scales as well. The reason is that it does not make sense to consider one criterion more important than another without considering the degree of variation among the consequences for the alternatives under trade-off analysis.

The swing matrix is a matrix where the top row defines the value measure importance and the left column represents the range of value measure variation. As recommended by Parnell, weights should descend in magnitude as we move in the diagonal from top left cell to the bottom right cell of the swing weight matrix. Multiple criteria can be placed in the same cell [27].

#### *4.5. Generate the Utility Function for Each Criterion*

For the trade study, it is necessary to represent, as part of the system model, the selected assessment criterion and the utility or value functions associated to each criterion. Utility functions can be discrete or continuous. Utility functions follow three basic shapes: linear, curve and S shape curve.

For example, when an increasing utility function is created for a particular criterion, the systems engineer ascertains whether the project stakeholders consider it as the minimum value of the measure to be accepted, mapping it to the 0 value on the score scale (*y*-axis). The measure beyond which an alternative provides no additional value is mapped to the highest score scale (*y*-axis). It is important to pick the appropriate inflection points for drawing the curve, which may be either convex or concave.

#### *4.6. Assess Each Alternative*

Every alternative should be estimated for a given criterion in terms of its score, based on the applied utility function. Then, the relative weights assigned to each criterion are used to compute the objective function that combines the weights and scores. The sum combining function is frequently used. The assessment is based on the created utility functions, using criteria values from each alternative, obtained from test data, vendor provided data, simulations, prototypes, engineering practice or literature.

#### *4.7. Show the Trade-Off Results*

Generally, a summary table of criteria versus alternatives is presented to summarize the results from the preceding steps. Based on these results, a decision is made. A sensitivity analysis is recommended to determine the robustness of the alternatives selected based on their highest rank in the trade-off analysis [27].

Concurrently, design heuristics are selected and used to refine the architecture. The use of heuristics is recommended to implement non-functional requirements that cannot be allocated to some building blocks but are to be implemented as design patterns or layouts of connectors between the building elements.

#### **5. Application of Trade-Off Analysis within the ISE&PPOOA MBSE Approach to Select the Best WPT Alternative for an Implanted Biomedical Device**

In this work, we are interested in the design of a challenging WPT system able to power a newly created micromotor for medical intravascular surgery. This demands the previous study of the best technological alternatives for its implementation. Many severe constraints are associated with this problem, including lack of data in the intended size range of 1 mm3. Due to this, a comprehensive design framework such as the ISE&PPOOA methodology is needed. We shall study the most promising WPT alternatives to power a 1 mm<sup>3</sup> MEMS device inside a human artery, using the information reviewed in the Section 3.

Below, we describe the main outcomes produced by performing the trade-off steps described in the previous section.

#### *5.1. Identify the System Modules for the Trade-Off Study*

The modular architecture previously obtained using ISE&PPOOA can be seen as an internal block diagram represented using SysML notation in Figure 2. This architecture was designed for a complex autonomous implanted system with sensors, actuators and communication capabilities. In order to simplify the trade-off study, we limit ourselves to the system parts related to WPT, functionality that can be reduced to the blocks "Power Source" and "Internal Electrical Power Generator".

**Figure 2.** Logical blocks of a biomedical implanted system with WPT.

*5.2. Identify Credible Alternative Technical Candidates for Implementing the System Logical Building Blocks or Modules under Consideration*

After a comprehensive literature review of the state-of-the-art WPT, summarized in the Introduction, possible candidates to power our MEMS-based device in the one cubic millimeter range are:


#### *5.3. Define Trade-Off Criteria*

Based on the system requirements [25] related to energy transfer and energy harvesting by the device inside the patient body, we select the following trade-off criteria:


costs. The higher the number of parts, the smaller their size, the more complex their shapes and the costlier their materials, the higher will be the mechanical complexity;

• Technical maturity (low–medium–high), an abstract measure of the degree of consolidation and performance of a technological solution.

Other safety criteria besides SAR and resilience impact the architecture at the connectors level, as well, so we recommend the use of heuristics instead of trade-off studies to implement these safety and resilience requirements.

Both the criteria selection process and the application of heuristics cannot be fully automated without human intervention, so a certain degree of subjectiveness, but also creativity, is unavoidable. Additionally, many unknown variables are present in the microdevice design, due to the exploratory and frontier research nature of this project. Thus, the trade-off criteria have been selected trying to constrain the space of alternatives (trade-space) as much as possible in order to produce a reasonable solution for a physical architecture, avoiding the introduction of spurious information that cannot be confirmed at present in the intended range of sizes. Another important consideration in the selection process was that the criteria values obtained from the literature revision could be reasonably extrapolated to our design objectives to build the utility curves. The final result of this process is the refined physical architecture of the best WPT alternative (APT) after the trade-off study, which can be seen in Figure 3.

It is not easy or even possible to know a priori which set of trade-off variables will be optimal for every problem, because it is the result of many constraints, especially the stakeholders' needs [25], which are usually presented in an ambiguous or less than desirable rigorous way. A very useful piece of advice to confirm that the criteria selection process is correct is based on the exploration of the resulting utility curves. This is a powerful approach to refine the criteria selection process itself in an iterative way, which is inherent to the MBSE design approach. From a mathematical point of view, the curves should be smooth (continuous and differentiable) and restricted to a few reasonable types, such as polynomial, logarithmic or exponential functions. The extremal and inflection points should indicate critical values of the corresponding trade-off criterion. Any lack of regularity, bijectivity (with suitable domain restrictions if necessary, like in exponential or even degree polynomials, for example, in order to be invertible) or oscillations indicate some serious problem with the utility curve and the associated trade-off variable, which must be corrected either by changing the criterion completely or revising the prescribed values. In our case, it can be observed that all utility curves in Figures 4–10 fulfill these mathematical conditions. These mathematical-based "metacriteria" for trade-off analysis verification are not explicitly described in the systems engineering literature, as far as we know, and can be considered a valuable contribution.

#### *5.4. Assign Relative Weightings to the Criteria*

We apply the swing weight matrix for our trade-study because it considers variation in the measured range as well as importance. Thus, the swing weight matrix is more complete than other approaches that only consider importance. The swing matrix obtained for the trade-off criteria selected in the previous step is shown in Table 1. Weights are assigned based on the technical literature review, previous experience with MEMS projects and stakeholders' needs. One common criticism of trade-off analysis based on weighting criteria outside the systems engineering community is that it involves some degree of subjectiveness. However, this human assessment is unavoidable, because there are many competing stakeholder interests, generally formulated in an incomplete and ambiguous manner, in addition to very complex physical and technical constraints. At present, no algorithm or artificial intelligence approach is able to complete this step automatically without human assistance. Similar observations can be applied to other widely used decision-making approaches such as Analytical Hierarchy Process.


**Table 1.** Swing matrix of the studied criteria. Weights are assigned based on technical literature values and stakeholders needs.

#### *5.5. Generate the Utility Function for Each Criterion*

Utility functions for the selected criteria can be seen in Figures 4–10. These utility functions are built and represented using the recommendations described in ISE&PPOOA book [3] admitting the mathematical conditions described in more detail in Section 5.3.

#### *5.6. Assess Each Alternative and Show the Trade-Off Results*

The results of this trade-off analysis can be seen in Table 2. Weights have been normalized with respect to the total weight sum of the swing matrix values (weight sum = 4.65), converting the swing weights into measure weights (or global weights), so that their sum is now 1. The final weighted sum for a given technological alternative is calculated as the sum of the products of the measure weights of each criterion by the values obtained from the utility curves for the same criterion. In other words, the values of the weights column are multiplied by the corresponding utility function values of a given alternative column and then summed to produce the final score. Comparing the obtained scores, we can select the preferred alternative. In this case, the Acoustic Power Transfer solution is the best, followed by the Radiating Far Field technology.

**Table 2.** Summary of results of the trade-off analysis. Weights are normalized with respect to the total weight sum of the swing matrix values (weight sum = 465).


#### **6. Results and Discussion**

The results indicate that the most appropriate technology is APT (Acoustic Power Transfer) due to its transfer effectiveness, size, effective operation distance and SAR. After we identified the technology, we were able to decompose in more detail the system blocks that have to implement it, and we used a BDD diagram made with standard SysML notation (Figure 3).

**Figure 3.** Refined physical architecture of the best WPT alternative (APT) after the trade-off study.

From the detailed decomposition of these blocks, it is possible to build, by using design heuristics and design patterns, the IBD diagrams that would define the refined architecture of the solution which, for brevity reasons, are not shown here.

Utility curves were obtained compiling relevant data from the revised scientific publications, as shown in Section 3, and verifying its mathematical correctness following the conditions explained in Section 5.3.

Specifically, data for the utility curves were collected from:

Table 2 (page 46) of reference [5];


Table 1 (page 2101) and Table 2 (page 2102) of reference [11];

Table 1 (page 9) of reference [13].

The values correspond to the best results of the state-of-the-art technologies combined. We used a quantitative safety criterion (SAR), but safety heuristics could be applied to refine the architecture as well.

A sensitivity analysis is recommended to determine the robustness of the alternatives. A complete sensitivity analysis is not performed here, but a very useful qualitative analysis can be performed from the information derived from the utility curves, which are far richer in content than their mere use as trade-off tools could suggest.

Input power (Figure 4) is a decreasing linear function in a semilogarithmic scale, which implies that its dynamic range is large. Thus, a small variation in other variables can produce a large variation in the required input power, making this variable very sensitive to small changes in the design for all the studied alternatives. This variable is also one of the most important, so that it can be deduced that the overall performance of the WPT system will be very sensitive. Thus, the input power values can vary for the same alternative with different values of the other variables. Moreover, input power is very relevant from a safety point of view, because it is limited to levels which cannot produce any damage, pain or discomfort on the skin or inside human tissues. Other safety considerations related to input power would be possible electromagnetic interferences and damage to other devices or to human operators.

**Figure 4.** Utility function for Input Power.

Power Transmission Effectiveness (PTE) (Figure 5) is a key variable with the largest range of variation of all studied criteria, even larger than the input power; it is the most sensitive value and the latest technological advances are critical to fix it. Its mathematical behavior is an increasing linear function in a semilog scale. For example, ultrasound-based energy harvesters only a few years ago could have a PTE as low as 0.001%, while the most recent ones can achieve 40%. Obviously, such a large increment has changed this alternative from a feasible WPT technology for medical implants to the most promising one in our trade-off analysis for MEMS. Thus, the two most sensitive variables, input power and PTE, have opposite contributions to the performance of the WPT system, because we want the largest PTE with the lowest input power.

**Figure 5.** Utility function for Power Transmission Effectiveness.

Implant size (Figure 6) is one of the most important constraints in our example due to an extreme constraint: it must be fit in a volume of the order of 1 mm3. However, PTE decreases exponentially and power input increases exponentially when the device size is decreased, even by small amounts, due to the small range of variation of this variable with a decreasing linear utility function. A sensitivity analysis would reveal that even small building tolerances in the manufacturing process with current technologies could radically change the performance of the final implant, which is unacceptable for a medical device.

**Figure 6.** Utility function for Implant WPTRx Size.

The effective operation distance (Figure 7) has a nonlinear response curve. This variable determines the depth of the implant. It can be observed that present technologies do not allow very deep implants in practice if millimeter size devices are desired. The range of variation is small, so that small variations in depth or distance from the external power source or both can produce very large changes in input power and PTE. However, in this case, the sensitivity is even worse than in the case of the device size, because the utility function increase is not linear and the steepest variation can be seen in the smallest distance values, in the range of a few mm. This implies that in order to have a stable power supply, very strict positioning mechanisms should be used, which could be impractical or even unfeasible if the implant can or must be moved.

**Figure 7.** Utility function for Effective Operation Distance.

SAR (Figure 8) is one of the most important safety related metrics, although the input power has to be limited as well. We can see that its range is severely constrained and most of the state-of-the-art devices are dangerously near the limit for the best values of the other variables. Thus, a sensitivity analysis would indicate that even a modest decrease in the SAR values in order to comply with legal restrictions would enormously impact the values of input power and PTE, surely forcing the implant size to be considerably larger, frustrating the goal of 1 mm3 total volume if radically new design ideas are not found.

**Figure 8.** Utility function for Specific Absorption Rate.

The last two variables, mechanical complexity (Figure 9) and technical maturity (Figure 10), are very different from the previous ones. They are abstract metrics that try to describe the effort and cost of building the WPT system without using economic or monetary values, presently unknown, because we are dealing with a present research project pushing the limits of near-future MEMS technology. In fact, their values are more categorical than numerical, although an easy conversion can be done in order to draw the utility curves, which are simply linear in these cases. In this study, these two variables are not as critical, but within a constrained project budget, much more detailed utility curves should be obtained in this regard. Surely, with accurate costs, their influence would severely constrain the feasible technological options and the conclusions of the trade-off analysis could be very different.

**Figure 9.** Utility function for Mechanical complexity.

**Figure 10.** Utility function for Technical maturity.

#### **7. Conclusions**

The combination of MBSE and trade-off analysis is a very useful engineering best practice for the selection of the most suitable technology to be applied to solve a particular problem. Here, MBSE allows us to identify the system level where the application of tradeoff should be accomplished. In this way, we identified the system building blocks and which functions should be implemented to solve the problem of a wireless energy transfer and energy harvesting system for a new type of micromotor inside the patient's body.

The results obtained by this trade-off study are consistent with the results obtained by other researchers and published in the literature, but expand and detail them within the context of a real research systems engineering effort with many unknowns because it pushes the limits of present MEMS technologies to sizes of the order of 1 mm3. The methodological approach we propose here helps to inform better design decisions in the early phases of the project, saving costly redesign efforts in later phases. Moreover, we have detailed the mathematical conditions that the utility functions must have in order to be useful for trade-off studies. Finally, we have shown how the information from the trade-off analysis can be used to make a valuable qualitative sensitivity study of the system.

**Author Contributions:** Conceptualization, J.A.M.R. and J.L.F.; methodology, J.L.F.; formal analysis, J.A.M.R. and J.L.F.; investigation, J.A.M.R. and P.L.L.E.; resources, J.A.M.R. and R.S.M.; data curation, R.S.M. and P.L.L.E.; writing—original draft preparation, E.D.-J. and J.A.M.R.; visualization, P.L.L.E.; supervision, E.D.-J. and J.L.F.; project administration, E.D.-J.; funding acquisition, E.D.-J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work has been conducted for the UWIPOM2 project, which received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 857654.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors want to recognize the work of Alba Martínez Pérez during preparation of the figures and text edition.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Detecting Axial Ratio of Microwave Field with High Resolution Using NV Centers in Diamond**

#### **Cui-Hong Li 1,2,\*, Deng-Feng Li 3, Yu Zheng 3, Fang-Wen Sun 3, A. M. Du 1,2 and Ya-Song Ge 1,2**


Received: 7 April 2019; Accepted: 14 May 2019; Published: 21 May 2019

**Abstract:** Polarization property characterization of the microwave (MW) field with high speed and resolution is vitally beneficial as the circularly-polarized MW field plays an important role in the development of quantum technologies and satellite communication technologies. In this work, we propose a scheme to detect the axial ratio of the MW field with optical diffraction limit resolution with a nitrogen vacancy (NV) center in diamond. Firstly, the idea of polarization selective detection of the MW magnetic field is carried out using a single NV center implanted in a type-IIa CVD diamond with a confocal microscope system achieving a sensitivity of 1.7 <sup>μ</sup>T/√Hz. Then, high speed wide-field characterization of the MW magnetic field at the submillimeter scale is realized by combining wide-field microscopy and ensemble NV centers inherent in a general CVD diamond. The precision axial ratio can be detected by measuring the magnitudes of two counter-rotating circularly-polarized MW magnetic fields. The wide-field detection of the axial ratio and strength parameters of microwave fields enables high speed testing of small-scale microwave devices.

**Keywords:** MW magnetic field; axial ratio; polarization; NV center

#### **1. Introduction**

Building and measuring circularly-polarized microwave (MW) fields are of great importance in many fields like magnetic resonance [1,2] and mobile satellite communications [3–5]. Conventionally, precision measurement of MW field strength is realized by converting the field strength to easilymeasurable electrical signals like using the calorimetric method [6] and the technique of peak demodulation [7]. The antenna axial ratio measurement is conducted with the help of a linearlypolarized antenna [4]. However, the calorimetric method is susceptible to environment temperature, and the technique of peak demodulation does not perform well at weak field strength. The assistant antenna perturbs the properties of the antenna under test (AUT) and usually has a large volume. Quantum-based MW field detection technologies, such as super quantum interference device [8] and cold atoms [9], which offer high sensitivity, require extreme measurement conditions. Vapor cell devices [10,11] that are capable of high sensitivity MW field sensing at near room temperature have a spatial resolution confined to sub-100 μm.

Recently, quantum sensing based on a point defect, the nitrogen vacancy (NV) center, in diamond has attracted wide attention for its long coherence time at room temperature [12]. It has already been deeply researched and applied to magnetic field, electric field, and temperature sensing [13–15] in many subjects. As its electron spin state can be coherently manipulated by the MW field, it can also

be applied to MW field sensing. In previous works, the NV-based MW magnetic field sensing schemes mainly focused on detecting MW magnetic field strength [16,17].

In this paper, we propose a scheme to detect the axial ratio of an MW device with high resolution using NV centers in diamond. The axial ratio detection is based on MW magnetic field strength detection of the two counter-rotating circularly-polarized components of an MW field. An NV center is a negatively charged point defect in diamond consisting of a substitutional nitrogen atom and an adjacent carbon vacancy. Its ground electron-spin triplet state can be selectively and coherently manipulated by a specified circularly-polarized MW field. The MW magnetic field strength of specific polarization can be derived from Rabi frequency or magnetic resonance spectra contrast [16,18]. Firstly, the MW axial ratio detection is demonstrated with a single NV center in diamond with a conventional confocal microscope. Then, wide-field detecting of the MW magnetic field, which greatly promotes sensing speed, is conducted with a wide-field microscopy system. The MW axial ratio detection technique facilitates the test and design of MW devices aimed at quantum spin manipulation [19–21] and satellite communication.

#### **2. Principle**

The axial ratio is an important parameter to describe a polarized MW field. Broadly, an arbitrarilypolarized MW field can be treated as an elliptically-polarized wave. The arbitrarily- polarized MW field spread along the z direction can be described as a combination of two circularly-polarized microwave fields, **B**<sup>+</sup> and **B**−, spread along the z direction with the same frequency as the source MW field that rotate in different directions in the xoy plane (Figure 1).

$$\mathbf{B} = \mathbf{B}\_{+} + \mathbf{B}\_{-}, \tag{1}$$

where *<sup>B</sup>*<sup>+</sup> (*B*−) indicates the magnetic field strength of the *<sup>σ</sup>*<sup>+</sup> (*σ*−) circularly-polarized field. The axial ratio can be described as the ratio of the magnitudes of the major and minor axes defined by the magnetic field vector. It can be figured out that the length of the major (minor) axis of the ellipse equals 2|*B*<sup>+</sup> + *B*−| (2|*B*<sup>+</sup> − *B*−|).

$$AR = |\frac{B\_+ + B\_-}{B\_+ - B\_-}|\,. \tag{2}$$

Obviously, *AR* = 1 for a circularly-polarized MW field, and *AR* = ∞ for a linearly-polarized MW field.

**Figure 1.** At a fixed point in space (or for fixed z), the magnetic vector **B** of a polarized microwave (MW) field traces out an ellipse in the xoy plane.

Axial ratio detection of MW fields with NV centers in diamond is possible as the strength of two counter-rotating MW magnetic fields can be selectively acquired by a spin-1 system. NV electron spin in its ground state can be optically pumped to the excited state and emit red fluorescence while decaying to the ground spin state.The spin state selective intersystem crossing (ISC) process [22] in the NV center enables the state read out and state initialization of NV electron spin. The NV electron in the *ms* = 0 state emits more photons then the NV center in the *ms* = ±1 state. Successive illumination initializes the NV center to the *ms* = 0 state. As shown in Figure 2a, the electron spin ground triplet state of NV centers has an energy splitting of *D* = 2.87 GHz between the *ms* = ±1 state and *ms* = 0 state under zero magnetic field. The spin transition *ms* = 0 ↔ *ms* = −1(+1), which can be simply treated as a two-level system (TLS) and can be selectively driven by *σ*−(*σ*+) circularly-polarized microwave field perpendicular to the NV axis [20]. According to the semi-classical theory of light–matter interaction, the transition frequency, i.e., Rabi frequency Ω*R*, and the optically-detected signal contrast *CR* that corresponds to the electron spin state population of the TLS under a driving MW magnetic field can be written as:

$$
\Omega\_R = \sqrt{\triangle^2 + \Omega\_{0^\prime}^2} \quad \mathcal{C}\_R = \mathcal{C}\_{R0} \frac{\Omega\_0^2}{\Omega\_0^2 + \triangle^2} \tag{3}
$$

where = *ω*<sup>0</sup> − *ω<sup>m</sup>* is the frequency difference between the frequency corresponding to the energy separation *h*¯ *ω*<sup>0</sup> of the TLS and the frequency of the driving microwave field *ωm*. Ω<sup>0</sup> = *γeBmw* is the Rabi frequency at resonance, i.e., *ω*<sup>0</sup> = *ωm*; *γ<sup>e</sup>* = 2.8 MHz/G is the gyromagnetic ratio of NV electron spin. *Bmw* is the amplitude of the resonant component of the driving microwave magnetic field. *CR*<sup>0</sup> is the maximum signal contrast between the *ms* = 0 state and *ms* = ±1 state. Thus, the MW magnetic field strength of particular polarization can be deduced from Rabi frequency.

Furthermore, due to the Zeeman effect, the resonance frequency of transition *ms* = 0 ↔ *ms* = ±1 moves to *ω*<sup>±</sup> = *Dgs* ± *γeBz* under a bias magnetic field *Bz* along the NV axis (Figure 2a). This enables wide-band detecting of the microwave magnetic field. Simply by turning over the permanent magnet to apply bias magnetic fields in reversed direction, the strength of both *σ*<sup>−</sup> and *σ*<sup>+</sup> circularly-polarized components of the MW magnetic field at a certain frequency perpendicular to the NV axis that rotate in different directions can be extracted by detecting Rabi frequencies. Therefore, the axial ratio of the MW magnetic field at the location of the NV center perpendicular to the NV axis can be deduced.

**Figure 2.** (**a**) The ground electron spin triplet state of the NV center with zero field splitting of *D* = 2.87 GHz between the *ms* = 0 state and *ms* = ±1 state. A magnetic field of *Bz* along the direction of the NV axis leads to an energy level splitting of 2*γeBz*, where *γ<sup>e</sup>* = 2.8 MHz/G is the gyro-magnetic ratio of NV electron spin. The spin transition *ms* = 0 ↔ *ms* = −1(+1) can be selectively addressed by the *σ*−(*σ*+) polarized MW field. (**b**) Scanning picture of a diamond sample containing a single NV center. The marked bright spot is the single NV center used for MW field detection in this experiment. The other bright signals below are from ensemble NV centers. (**c**) Optically-detected magnetic resonance spectra of nitrogen vacancy (NV) center under a bias magnetic field of about 500 G. The data are fitted with the Gauss function. The first (second) dip at *f*<sup>1</sup> = 1.44 GHz (*f*<sup>2</sup> = 4.30 GHz) corresponds to the transition *ms* = 0 ↔ *ms* = −1 (+1) state.

#### **3. Experiments and Results**

#### *3.1. MW Magnetic Field Detection Using a Single NV Center in Diamond*

The detection of the MW magnetic field is firstly demonstrated with a single NV center in diamond. The experiment setup was based on a home-built confocal microscope [23]. A 532-nm green laser was used to initialize and detect the NV center. A microscope objective (NA = 0.9, Olympus, Kyoto, Japan) was used to focus the laser and collect fluorescence photons. The diamond sample was a single crystal synthetic (100)-oriented electronic-grade type-IIa diamond with dimension of <sup>2</sup> <sup>×</sup> <sup>2</sup> <sup>×</sup> 0.5 mm3 from the Element Six company. The diamond sample typically had less than a 0.03 ppb NV concentration before implantation. The investigated individual NV center was generated by <sup>14</sup>*N* ion implantation. The implantation dosage was 10<sup>11</sup> cm−2, and the estimated average depth of NV was about 20 nm. The MW field (R&S SMB 100A signal generator, Munich, Germany) was delivered to the sample via a coplanar waveguide (CPW) fabricated on the top surface of the diamond substrate. The diamond sample was mounted on a NanoCube Piezo System (PI 611.3S, Karlsruhe, Germany) enabling nano-scale scanning. A permanent magnet mounted on a three-axis translation stage applied a bias magnetic field of between zero and roughly 600 G, aligning with the NV axis (z axis). Figure 2b is a fluorescence picture obtained by scanning the diamond sample. The marked bright spot in the picture is a single NV center chosen for MW magnetic field sensing. Second order photon correlation measurement of the bright spot with strong antibunching observed at zero delay *g*2(0) ∼ 0.11 < 0.5 is proof that a single NV center was investigated [24] (data not shown).

To conduct polarization selective MW magnetic field sensing, a bias magnetic field was applied to split the resonant spectra lines corresponding to *ms* = 0 ↔ *ms* = −1 and *ms* = 0 ↔ *ms* = +1 transitions. In our experiment, a bias magnetic field of about 500 G was applied. As shown in Figure 2c, the optically-detected magnetic resonance (ODMR) spectra revealed an energy level splitting of 2.86 GHz between the *ms* = −1 state and *ms* = +1 state. Although the MW field from the signal generator was linearly polarized, the transition *ms* = 0 ↔ *ms* = −1 was only sensitive to the *σ*<sup>−</sup> MW field when the MW magnetic field strength was not too strong according to the rotating wave approximation in quantum optics [19]. We implemented the measurement of the strength of the *σ*− circularly-polarized component of the MW magnetic field with the NV center by detecting Rabi frequency. The MW field frequency was adjusted in resonance to the *ms* = 0 ↔ *ms* = −1 transition. The Rabi frequency detection sequence is shown in Figure 3a. A laser pulse was firstly applied to initialize the NV electron to the ground *ms* = 0 state. Then, an MW manipulation pulse was inserted. Finally, the population of the electron spin state was detected by another laser pulse. To obtain the manipulation frequency of the MW field, the MW pulse varied in time duration. Moreover, a photon count read pulse-2was applied after the NV center was polarized to the *ms* = 0 state again to provide a reference for drifts of the NV center. Figure 3b shows the Rabi oscillation detected with the MW field source power of about 7.8 mW applied. The oscillation signal was fitted by 1 <sup>−</sup> *<sup>C</sup>*0*e*−*t*/*t*<sup>0</sup> *sin*(*wt* <sup>+</sup> *<sup>φ</sup>*). The manipulation frequency of 0.63 MHz indicated the *σ*− polarized MW magnetic field strength of *B*<sup>−</sup> = 0.22 G at the NV center spot perpendicular to the NV axis. According to the fitting error of *δB*<sup>−</sup> = 0.0024 G and detecting time of *tmeas* = 1000 s of the Rabi nutation curve, the sensitivity for MW field detection was *η* = *δB* <sup>√</sup>*tmeas* <sup>=</sup> 1.7 <sup>μ</sup>T/√Hz. The Rabi oscillation frequency can also be derived by Fourier transformation (Figure 3c).

**Figure 3.** (**a**) Rabi frequency measurement sequence. A first green pulse of 3 μs was applied to initialize the NV electron spin to the *ms* = 0 state; then, the MW field was applied to interact with the NV electron spin; finally, the electron spin state was determined by simultaneously applying the laser pulse and photon counts read pulse (0.3 μs). To avoid the drift effect of the NV center on determining the NV state, the spin state photon counts were normalized to another photon counts read pulse that was applied in succession to the first one after the NV center was polarized again. (**b**) Rabi oscillation signal driven by a resonant MW field of 1.44 GHz. Each data point on the diagram was obtained by repeating the measurement sequences one million times. The oscillation frequency of 0.63 MHz indicates an MW magnetic field strength of *Bσ*<sup>−</sup> = 0.22 G at the NV spot. The fitting error yields a sensitivity of 1.7 <sup>μ</sup>T/√Hz. (**c**) Fourier transformation of (**b**).

As *BMW* <sup>∝</sup> <sup>√</sup>*P*, *<sup>ω</sup>MW* <sup>=</sup> *<sup>γ</sup>eBMW*, the Rabi frequency increased linearly with the square root of the resonant MW field power. The linear dependence was verified by detecting Rabi oscillation signals with different MW field powers applied. Figure 4 shows the linear dependence of Rabi frequency on microwave magnetic field strength. This is consistent with the measurement principle.

Turning over the permanent magnet to apply a bias magnetic field in the reverse direction, the *σ*<sup>+</sup> component of the MW magnetic field perpendicular to the NV axis can be detected similarly. Therefore, the axial ratio of the MW magnetic field at the NV point perpendicular to the NV axis can be derived.

**Figure 4.** Rabi frequency versus MW source power applied. The data are fitted by a linear function.

#### *3.2. Wide-Field MW Magnetic Field Imaging with Ensemble NV Centers*

A single NV center is capable of detecting the strength of the particular circularly-polarized component of the MW magnetic field with high spatial resolution and high sensitivity. Combining the wide-field microscopy system and ensemble NV centers, the detection speed can be crucially improved, and the spatial resolution is only confined to the optical diffraction limit, which is generally at the micron or sub-micron degree. For many practical applications, wide-field microscopy that offers MW field strength imaging at a short time scale and requires simpler experimental technologies could be highly beneficial.

For wide-field detection of the MW magnetic field, the experiment was based on a home-built wide-field back focal imaging microscope system. As shown in Figure 5a, the green laser was focused at the back focal area of the object (*NA* = 0.7) to illuminate the sample with collimated light. The red fluorescence emitted by ensemble NV centers at the focal plane was collected by the same objective and was separated from the source light by a dichroic mirror. Then, the light was filtered and focused on an EMCCD (Andor Ultra iXon 897, Oxford, UK) for imaging. The spatial resolution of the microscope system was about 700 nm, as the wavelength of fluorescence emitted from the NV center was roughly between 600 nm and 750 nm. The sample used in this experiment was a 2.6 <sup>×</sup> 2.6 <sup>×</sup> 0.3 mm<sup>3</sup> general single crystal (100)-oriented CVD diamond from Element Six company with inherent ensemble NV centers. The inherent NV concentration in the diamond was estimated to be about 3 <sup>×</sup> 1014 cm−<sup>3</sup> (2 ppb) by comparing the fluorescence photon counts of the NV centers in the diamond to that of a single NV center in type-IIa diamond under the confocal microscope system. The fluorescence emitted from the diamond was ascribed to the NV− center according to its zero photon line at 637 nm detected by the spectrometer FHR 640 (data not shown). A metal micro-strip line with a width of 8 μm was fabricated on the top surface of the diamond to deliver the microwave field.

**Figure 5.** (**a**) Schematic diagram of wide-field fluorescence microscopy. The green light was focused at the back focal are of the object to excite NV centers with collimated light. The red fluorescence emitted by the NV centers was collected by the same object. The dichroic mirror was used to separate the laser and fluorescence. At last, the fluorescence was further filtered and focused on the EMCCD. (**b**) Wide-field fluorescence imaging of the diamond sample without the MW field fed in. The dark strip indicates the position of the metal MW strip line.

Figure 5b is a fluorescence picture taken by EMCCD with an exposure time of 2 *s*. The metal strip line on the top surface blocked the fluorescence from the NV center, as shown in Figure 5b. Rough detection of the MW magnetic field was carried out by comparing the signals with and without the resonant MW magnetic field applied [18]. At a zero magnetic bias field, the NV ground electron spin energy splitting was *Dgs* ∼ 2.87 GHz. Taking pictures of the same region as in Figure 5b with an MW field of *f* = 2.87 GHz applied, the NV fluorescence brightness decreased in the vicinity of the metal strip line (data not shown). The signal contrast decreased with the distance between NV centers and the MW strip line.

Taking pictures of the same region while scanning the MW frequency applied to the strip line, the ODMRsignal of the region can be extracted. To verify the wide-field MW magnetic field sensing

method, ODMR signals of Spot A were extracted with the MW magnetic fields of different source powers P transmitted to the MW strip line (Figure 6a). Obviously, the magnetic resonance signal contrast increased with the power of the MW field applied, as the red square in Figure 6b indicates [18]. Thus, simply by tuning the NV electron spin transition frequency in resonance with the MW field, the qualitative MW magnetic field strength of specific polarization and spread direction [19] can be determined from the NV fluorescence signal contrast. Besides, the blue dot in Figure 6b shows that the full width at half maximum (FWHM) of the NV magnetic resonance signal decreased with reduced MW source power till a threshold value. The minimum FWHM of NV ODMR spectra reflects the NV electron spin decoherence property [18]. Moreover, the NV axial directions in the diamond were the same as the four tetrahedral diamond axes [25]. The MW magnetic field spread along each NV axis can be detected. Thus, the axial ratio of an MW field spread along an arbitrary direction can be extracted. Further applying the Rabi measurement sequence in the wide-field microscopy system, precise detection of MW magnetic field strength and the axial ratio is attainable [26,27].

**Figure 6.** (**a**) ODMRspectra of Spot A in Figure 5b extracted with four different MW sources' power fed in. The line width and dip contrast increase with the MW source power. (**b**) ODMR signal contrast (red square) and FWHM (blue dot) versus MW source power extracted from ODMR spectra of Spot A with different MW source powers fed in.

#### **4. Conclusions**

In summary, we proposed a new method that can detect the axial ratio of microwave fields with high resolution for wide-field and a short time scale. The measurement method relied on detection of two circularly-polarized MW magnetic fields that rotate in different directions in the plane perpendicular to the NV axis. The axial ratio detecting of the MW magnetic field was firstly demonstrated using a single NV center in diamond, achieving a sensitivity of 1.7 <sup>μ</sup>T/√Hz. Then, wide-field detection of the MW magnetic field with an imaging field size of 85 <sup>×</sup> <sup>85</sup> <sup>μ</sup>m<sup>2</sup> was displayed. The single NV center in the type-IIa CVD diamond investigated was generated by ion implantation, and the ensemble NV centers in the general CVD diamond were inherent. In a future study, we will apply this technique to test the polarization-controllable MW devices that we designed for selective manipulation of NV sublevels.The method we proposed can effectively benefit the design of MW fields with special polarization properties at the sub-millimeter to millimeter scale and thus benefit the development of techniques relying on quantum spin manipulation [28–30] and a small-sized antenna.

**Author Contributions:** Conceptualization, C.-H.L.; methodology, C.-H.L.; software, Y.Z.; validation, D.-F.L., A.M.D.; formal analysis, C.-H.L.; investigation, C.-H.L.; resources, F.-W.S.; data curation, D.-F.L.; writing—original draft preparation, C.-H.L.; writing—review and editing, C.-H.L., A.M.D.; visualization, C.-H.L.; supervision, F.-W.S., A.M.D.; project administration, Y.-S.G.; funding acquisition, Y.-S.G., A.M.D.

**Funding:** This research was funded by "Development of High Sensitivity Diamond NV Magnetic Sensor" in Chinese Academy of Sciences Research Equipment Development Project, the Strategic Priority Program in Space

Science of CAS grant number XDA1502030105, the Strategic Priority Research Program of the Chinese Academy of Sciences grant number XDA14040204, the National Key Research and Development Program of China grant numbers 2016YFB0501300 and 2016YFB0501304 and the APC was funded by "Development of High Sensitivity Diamond NV Magnetic Sensor" in Chinese Academy of Sciences Research Equipment Development Project.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Letter* **Terahertz Frequency-Scaled Di**ff**erential Imaging for Sub-6 GHz Vehicular Antenna Signature Analysis**

**Jose Antonio Solano-Perez 1,\*, María-Teresa Martínez-Inglés 2, Jose-Maria Molina-Garcia-Pardo 1, Jordi Romeu 3, Lluis Jofre-Roca 3, Christian Ballesteros-Sánchez 3, José-Víctor Rodríguez <sup>1</sup> and Antonio Mateo-Aroca <sup>4</sup>**


Received: 25 August 2020; Accepted: 28 September 2020; Published: 2 October 2020

**Abstract:** The next generation of connected and autonomous vehicles will be equipped with high numbers of antennas operating in a wide frequency range for communications and environment sensing. The study of 3D spatial angular responses and the radiation patterns modified by vehicular structure will allow for better integration of the associated communication and sensing antennas. The use of near-field monostatic focusing, applied with frequency-dimension scale translation and differential imaging, offers a novel imaging application. The objective of this paper is to theoretically and experimentally study the method of obtaining currents produced by an antenna radiating on top of a vehicular platform using differential imaging. The experimental part of the study focuses on measuring a scaled target using an imaging system operating in a terahertz band—from 220 to 330 GHz—that matches a 5G frequency band according to frequency-dimension scale translation. The results show that the induced currents are properly estimated using this methodology, and that the influence of the bandwidth is assessed.

**Keywords:** frequency-dimension scale; terahertz; measurements; differential imaging

#### **1. Introduction**

The future generation of connected vehicle, along with the autonomous vehicle it evolves into, will require significantly increasing the number of antennas on its surface operating at different frequency bands from sub-6 GHz to millimeter-wave (mmWave), as effective communication and environment sensing are ensured in this way in respect of other vehicles and different base-stations [1].

The operation of these antennas and a focus on a 3D spatial angular response (3D radiation pattern) may be significantly influenced (perturbed) by the nearby vehicular structure. The task of studying these effects may require ascertaining the distribution of the currents induced by the antenna on the vehicular surface by means of numerical or experimental methods. The process of obtaining these currents for realistic vehicle geometries, especially when using experimental techniques, may require complex and bulky setups.

The new imaging systems operating in terahertz frequencies provide a very interesting tool that is applicable for solving the stated problem in relation to the estimation of the induced currents on the vehicular structure by the antenna. This imaging system provides the reconstructed image and shape obtained via the scattered field of metallic or non-metallic targets [2].

The concept of using terahertz frequencies and imaging techniques has been reviewed in order to confirm its applicability to the present problem. In [3], the ultrawide-band (UWB) imaging radar system used in this paper is tested at mmWave and terahertz bands, validating the spatial resolution on the order of millimeters and its imaging capabilities. The work presented in [4] is an imaging system that inspects a polymer radome in the frequency range of 70–170 GHz, reporting good performance in resolution and proposing to increase the frequency up to 300 GHz for improved resolution. Continuing with frequencies around 300 GHz, the three-dimensional (3D) imaging system presented in [5] is a synthetic aperture radar (SAR) operating at 340 GHz, utilizing the Fourier transform in two dimensions for the reconstruction process. In [6], the operating frequency is increased to 522 GHz. Finally, to conclude the background review, a novel reconstruction process in time domain operating in the microwave and mmWave frequency bands is presented in [7].

Extensive discussion has focused on the structural scattering and the scattering of an antenna's radiation mode [8]. It is intended to reconstruct the currents associated with scattering of the radiation mode of an antenna.

The main objective of this paper is to investigate the distribution of currents induced by the antenna on the vehicular surface by means of differential image reconstruction processes, and using a previously developed multi-frequency system [3] operating in the terahertz frequency band and extending this concept to explore an area. This experimental method is a non-invasive way of obtaining the induced currents without placing the antenna in the supporting structure. A frequency-dimension scale translation is applied to reduce complexity from an experimental point of view. As the defining aspect, a differential imaging concept is exploited to process the images captured during the experiment. In this specific case of measurements, it could be applied to the study and optimization of the co-location of several antenna systems in vehicles. To the authors' best knowledge, the uniqueness of this paper concerns dealing with a differential imaging method to obtain the current distribution without influencing the measure.

As discussed previously, the terahertz imaging radar system makes it possible to obtain the spatial distribution of the reflectivity of an object for the analytical and experimental assessment of induced currents in the vehicle structure. This system uses a double focal system with an independent transmitter antenna and receiver antenna, being connected to both ports of a vector network analyzer (VNA) to obtain the scattering parameters over a defined frequency range [3]. Using frequency-dimension scale translation, and as a proof-of-concept, the behavior of a car with a length of 4.5 m operating at the new 6 GHz (3.5 GHz) 5G frequency band was studied through a 1:43 scale metallic car prototype of 10.5 cm length, measured at the R-band (220 GHz–330 GHz). The 1:43 scale metallic car corresponds with a frequency range from 2.5 to 3.8 GHz. The results are differential measurements calculated from the difference between the scattering image reconstructed by the imaging system with the antenna (ON state or short-circuited monopole antenna) and without the antenna (OFF state or open-circuited monopole antenna). Then, an analysis of the differential image provides key information about the current car surface distribution.

This paper is structured as follows: Section 2 jointly develops the theoretical formulation—related to the total image and differential image reconstruction process with the measurement description. Section 3 includes the total image reconstruction for the different measurement cases. Section 4 shows the results related to differential imaging. Finally, Section 5 presents the conclusions.

#### **2. The Imaging Process**

In this section, the two-step process of imaging is formulated, first in terms of total image, and then in a novel differential form, to apply to the reconstruction of the currents created by a radiating antenna on top of a vehicle.

#### *2.1. Formulation of the Total Image Reconstruction Process*

The ultrawide-band (UWB) near-field imaging radar system used for the measurement is a flat SAR system composed of a transmitter antenna (*Tx*) and a receiver antenna (*Rx*) with fixed spatial separation, named as *X* offset distance, being a monostatic configuration but with a bi-static angle. The *Tx*–*Rx* system is fixed and explores the environment by moving the object in an *XY* plane using a positioning table. The system points toward the target located at a distance *rt* in the *Z*-axis, as shown in Figure 1.

**Figure 1.** Measurement arrangements to explore the object in an *XY* plane using a positioning table.

Figure 1 shows the measurement arrangements, depicting the *Tx*–*Rx* and the object.

The objective is to get the spatial distribution image of the electric contrast associated with the object under test. The electric contrast is written as *c* → *r* = ε*t* → *r* −ε*<sup>b</sup>* <sup>ε</sup>*<sup>b</sup>* , with ε*<sup>t</sup>* → *r* being the object's complex dielectric permittivity relative to its environment (background) constant value ε*<sup>b</sup>* → *r* .

The imaging system generates an illuminating field with polarization in the *y*-axis (transversal) that is incident on the object under investigation, represented by a metallic car in Figure 1. → *E* → *r*, *f*; → *rti* denotes the electric field generated in point <sup>→</sup> *rti* by the *Tx* antenna placed at <sup>→</sup> *rt* and operating at frequency *f*.

$$
\overrightarrow{J}\_{eq}(\overrightarrow{r},f;\overrightarrow{r}\_{t\_i}) = j\omega\varepsilon\_b \ c(\overrightarrow{r})\overrightarrow{E}(\overrightarrow{r},f;\overrightarrow{r}\_{t\_i})\tag{1}
$$

The equivalent current <sup>→</sup> *<sup>J</sup> eq* is the scattering field source generated by the object. Then, <sup>→</sup> *J eq* is understood as a "trace" of the original object (image).

Figure 1 shows the *Tx* and *Rx* antennas fixed in position, and the object is moved using *XY* plane geometry located at a distance of <sup>→</sup> *rt*, on a set of positions *Nx*, *Ny* across *x*-axis and *y*-axis, to gather a 2D array that includes monostatic measurements.

For targets with highly conductive behavior (σ ωε0ε *<sup>r</sup>*), like the metallic car used in this paper, the equivalent currents tend to produce the conduction current <sup>→</sup> *J eq* = σ → *E* [9].

The equivalent current <sup>→</sup> *<sup>J</sup> eq*<sup>→</sup> *rt*, *f*; → *rti* generates a scattered field (<sup>→</sup> *Es*) in the point <sup>→</sup> *rti* of the object, which is measured by the *Rx* antenna at each point of measurement. <sup>→</sup> *Es* has the following equation:

$$\overrightarrow{E}\_s(\overrightarrow{r}\_{R\_j \prime} f; \overrightarrow{r}\_{t\_i}) = -j\omega\mu\_0 \int\_{V\_0} \overrightarrow{J}\_{eq}(\overrightarrow{r}\_{t \prime} f; \overrightarrow{r}\_{t\_i}) G\left(\left|\overrightarrow{r}\_{R\_i} - \overrightarrow{r}\_{t}\right|, f\right) dV\tag{2}$$

where <sup>ω</sup> is the angular frequency, <sup>μ</sup><sup>0</sup> is the magnetic permeability in vacuum, <sup>→</sup> *J eq* is the equivalent current, and *G* → *<sup>r</sup> Ri* <sup>−</sup> <sup>→</sup> *rt* , *f* is the Green's function that matches the geometry of the problem. For a 3D measurement geometry, Green's function is written as *G*(*r*) = *<sup>e</sup>* −*jkbr <sup>r</sup>* , where *kb* <sup>=</sup> <sup>ω</sup> <sup>√</sup>μ0ε0ε*b*.

For the case where the object is uniformly illuminated using high-directivity antennas, <sup>→</sup> *Es* is written as:

$$\overrightarrow{E}\_s(\overrightarrow{r}\_{R\_j}, f; \overrightarrow{r}\_{t\_i}) = -k\_b^2(f) A\_d \int\_{V\_0} c(\overrightarrow{r}', f) G\left( \left| \overrightarrow{r}\_t - \overrightarrow{r}\_{t\_i} \right| \epsilon \right) G\left( \left| \overrightarrow{r}\_{R\_i} - \overrightarrow{r}\_{t\_i} \right| \epsilon \right) dV \tag{3}$$

where the parameters of the transmitter and receiver antennas are included in the complex constant named *Ad*.

The total image reconstruction process introduced by the previous works [3,9–13] is applied, characterized as a bi-focusing technique using multi-frequency. The contrast factor *c* → *r* is averaged through the complete terahertz band and to the entire reconstruction space, which can be expressed as

$$\mathcal{L}\left(\overrightarrow{r}\right) = A\_o \sum\_{f\_{\text{min}}}^{f\_{\text{max}}} \sum\_{i=1}^{N\_{\text{T}-R}} \frac{\overrightarrow{E}\_s(\overrightarrow{r}\_{R\_i}, f; \overrightarrow{r}\_{I\_i})}{k\_b^2(f)} e^{j\mathbf{k}\_b(f)\left(|\overrightarrow{r}\_{I} - \overrightarrow{r}\_{I\_i}| + |\overrightarrow{r}\_{R\_i} - \overrightarrow{r}\_{I}|\right)} \tag{4}$$

where the parameters of the transmitter and receiver antennas are included in the complex constant called *Ao*. In order to reduce the computation times, fast Fourier transform (FFT) is used to solve Equation (4).

#### *2.2. Formulation of Di*ff*erential Image Reconstruction Process*

The two-step process, described in the previous paragraph, is composed of a first step—the direct problem of obtaining the scattering fields produced by the external illumination responsible for the creation of the currents into the object—and a second step—an inverse problem focusing back on the previously obtained scattered fields to obtain an image of the currents created by the nearby illuminating geometry.

Here, the goal is to obtain the currents created not by the external illuminating geometry but by the vehicle antenna close to the object, but still using the same external illuminating geometry. These vehicular antenna currents, as described below, should be considered as a differential contribution for two different states: an ON state, or short-circuited monopole antenna, installing monopole in contact with vehicle surface, and an OFF state, or open-circuited monopole antenna, removing the monopole antenna.

To extract the differential impact of in-vehicle antenna radiation, the scattering fields will be collected and the corresponding image of the currents on the surface of the vehicle will be obtained for two successive states of the vehicle antenna: short-circuited and open-circuited [14].

In the first case (the vehicle antenna in the ON, short-circuited state), the total scattered fields → *E total ant*−*on* are produced by the combination of the scattered fields due to the currents created on the surface of the object by the external illuminating geometry <sup>→</sup> *E ext*−*illun ant*−*on* , the scattered fields due to the structural currents flowing into the vehicle antenna structure without interacting with its port <sup>→</sup> *E ant*−*str ant*−*on* , and the scattered fields due to the antenna's radiating mode, <sup>→</sup> *E ant*−*rad ant*−*on* [15]:

$$
\stackrel{\rightarrow}{E}\_{ant-ou}^{total} = \stackrel{\rightarrow}{E}\_{ant-ou}^{ext-illun} + \stackrel{\rightarrow}{E}\_{ant-ou}^{ant-str} + \stackrel{\rightarrow}{E}\_{ant-ou}^{ant-rad} \tag{5}
$$

In the second case (vehicle antenna in the OFF, open-circuited state), the total scattered fields → *E total ant*−*off* are produced by the combination of the scattered fields due to the currents created on the surface of the object by the external illuminating geometry <sup>→</sup> *E ext*−*illun ant*−*off* and the scattered fields due to the structural currents flowing into the vehicle antenna structure without interacting with its port <sup>→</sup> *E ant*−*str ant*−*off* .

$$
\stackrel{\rightarrow}{E}\_{ant-off}^{total} = \stackrel{\rightarrow}{E}\_{ant-off}^{ext-illun} + \stackrel{\rightarrow}{E}\_{ant-off}^{ant-str} \tag{6}
$$

For the compact vehicle antenna, resonant antennas are utilized in most vehicles. Then, it may be considered: (i) that interactions of a high order, such as reflections among external illuminating geometry and vehicle antenna or vehicular platform, may be neglected since they are much smaller than the rest; and (ii) that the scattered fields produced by the unchanged parts, like the vehicular platform and the structure of the antenna, are approximately the same for both states of the vehicle antenna: <sup>→</sup> *E ext*−*illun ant*−*on* <sup>→</sup> *E ext*−*illun ant*−*off* and <sup>→</sup> *E ant*−*str ant*−*on* <sup>→</sup> *E ant*−*str ant*−*off* .

Based on the previous assumption, the reconstructed currents above the vehicle surface reconstructed from the differential scattered fields <sup>→</sup> *E di f ant* can be established as:

$$
\stackrel{\rightarrow}{E}\_{\text{ant}}^{\rightarrow df} = \stackrel{\rightarrow}{E}\_{\text{ant-on}}^{\leftarrow total} - \stackrel{\rightarrow}{E}\_{\text{ant-off}}^{\rightarrow ant-rad} \cong \stackrel{\rightarrow}{E}\_{\text{ant-on}}^{\rightarrow ant-rad} \tag{7}
$$

being related to the currents produced by the radiating mode of the vehicle antenna, and Equation (4) becomes:

$$c\_{ant}^{diff}(\overrightarrow{r'}) = f(\overrightarrow{E\_{ant}}) \cong c\_{ant-on}^{total}(\overrightarrow{r'}) - c\_{ant-off}^{total}(\overrightarrow{r'}) \tag{8}$$

#### *2.3. Measurement Description and Configuration*

The measurement device was a vector network analyzer (VNA) (Rohde Schwarz model ZVA 67). One VNA port was connected to the frequency converter and the *Tx* antenna in which the converter was mounted. In the same way, the other VNA port was connected to the frequency converter and the *Rx* antenna was mounted in this converter as well. The frequency converters allow us to extend the VNA operating frequency band up to 220–330 GHz [16]. The setup configuration developed in [3] measured the S21 parameter (scattering parameter) in the terahertz band to reconstruct the image using a UWB near-field multi-frequency bi-focusing algorithm.

The geometry for the measurement depicted in Figure 1 is used to estimate the distribution of the induced currents on a metal car of 10 cm length—using a terahertz band for imaging—through frequency-dimension scale translation to apply the novel concept of differential imaging.

The *XY* planar measurement geometry is used for illuminating the metallic vehicle using the ultrawide-band (UWB) imaging radar system [3]. This illumination provides optimization of the induced currents on the vehicular surface, although the vertical monopole antenna is installed perpendicular but connected to the vehicle surface. One of the reasons for using a monopole antenna oriented in the direction of the *Tx*/*Rx* antennas is because, even when the excitation of the monopole mode may be reduced, its mutual coupling with the *Tx*/*Rx* antennas is also reduced and so the multiple reflections between them may certainly be neglected and, at the same time, the surface currents on top of the vehicle due to the antenna should be visible—especially when the *Tx*/*Rx* antennas move out of the perpendicular direction.

The object is placed on a radiofrequency-absorbing material located 0.5 m from the *Tx*/*Rx* antennas. The separation between the *Tx* antenna with a frequency converter and the *Rx* antenna with a frequency converter is 0.1 m, installed in a fixed position in order to avoid any perturbation in the measurement associated with coaxial cables and converter movements. The object under measurement is placed on a positioning table that is capable of moving in two axes (*x*-axis and *y*-axis) with steps of 1 mm [16]. Figure 2 displays the frequency converter with the horn antenna mounted and the positioning table with the metallic car.

**Figure 2.** Measurement setup.

The S21 parameter measurements were carried out by sweeping the frequency band from 220 GHz to 330 GHz. The two-frequency converter equipment was manufactured by Rohde & Schwarz (Munich, Germany)®, with the model being the ZVA-Z325 Frequency Converter, J-Band WR-03 [17]. Table 1 includes the VNA configuration parameters.



The *Tx* and *Rx* antennas were horns produced by Flann Microwave LTD (Cornwall, United Kingdom) [18]. The selected horn model was the 32240-25, with a bandwidth (BW) from 217 GHz to 330 GHz. The *E* and *H* plane gain was 23.70 dBi at 217 GHz and 26.99 dBi at 330 GHz. The E and H plane beamwidth was 11.9◦ at 217 GHz and 7.8◦ at 330 GHz.

The selected object is a metallic car with a 1:43 scale and 10 cm length. This is a good model for a real car, in terms of both shape and materials. The scale of the car corresponds with the frequency range between 2.5 and 3.8 GHz. A short 15 mm metal wire was placed on top of the car to act as an antenna, as shown in Figure 3. The measurement was performed with (ON state) and without the antenna (OFF state) mounted in the metallic car with the aim of estimating the induced current on the car after using the imaging system to illuminate.

**Figure 3.** A metallic car with an antenna for differential imaging.

#### **3. Results**

In the first step, the total images were obtained based on the two-step process by measuring the scattered fields first and then obtaining the focused field by application of Equation (4) for different frequency ranges and focusing the car with (ON state, short-circuited monopole antenna) and without the antenna (OFF state, open-circuited monopole antenna), as depicted in Figures 4–7.

**Figure 4.** Focused field amplitude of the metallic car without an antenna. Frequency range: 295 GHz– 305 GHz (bandwidth (BW) = 10 GHz).

The bandwidth used to reconstruct the image is an important imaging parameter. The first image is shown in Figure 4, which depicts the total image reconstruction without an antenna on top of the metallic car with a frequency of 295 GHz–305 GHz. The top of the metallic car has a curvature that generates a specific focused field distribution.

Then, Figure 5 presents the total image reconstructed without an antenna using a bandwidth of 110 GHz, sweeping from 220 GHz to 330 GHz. This bandwidth provides better resolution performances, meaning the top of the metallic car can be identified in the image as a very clear reconstruction. Additionally, Figure 5b highlights the focused field amplitude value and the distribution due to the reflection in the top of the car, including some minor contributions around the top of the car.

**Figure 5.** The focused field amplitude of a metallic car without an antenna. Frequency range: 220 GHz– 330 GHz (BW = 110 GHz): (**a**) *XY* plane representation; (**b**) 3D representation of the focused field.

Figure 6 shows the total image reconstructed with the antenna installed on the top and using a bandwidth of 110 GHz, from 220 GHz to 330 GHz. The maximum of the focused fields is identified around the center of the top of the metallic car, where the antenna is installed. The antenna provides reflectivity, but it is masked with the reflectivity of the metallic car top in the graph.

**Figure 6.** The focused field amplitude of a metallic car with an antenna. Frequency range: 220 GHz– 330 GHz (BW = 110 GHz).

Finally, Figure 7 depicts the total image with the antenna installed on the top, using a bandwidth of 10 GHz around 300 GHz. The shape of the focused field is like that in Figure 4 due to the antenna contribution being added to the reflectivity of the metallic car top.

**Figure 7.** *Cont*.

**Figure 7.** The focused field amplitude of a metallic car with an antenna. Frequency range: 295 GHz– 305 GHz (BW = 10 GHz). (**a**) *XY* plane representation; (**b**) 3D representation of the focused field.

#### **4. Discussion**

The measured fields are focused on different image planes along the *z*-axis (*z* = 0.4787 m) where the capability of the system for focusing on different longitudinal distances (along the *z*-axis) and different frequency ranges and bandwidths, over different transversal planes (*XY* planes), may be observed. The experimental results show that using a wide bandwidth (110 GHz) is not required in order to obtain a good performance. The performance was of sufficient quality using 10 GHz bandwidth around 300 GHz. The resolutions would be:

$$
\Delta z \sim 1/\_{\text{AB}}\tag{9}
$$

$$
\Delta \mathbf{x}\_\prime \Delta \underline{y} \sim^{\lambda\_0} \sqrt{(2\sin \theta\_0)} \tag{10}
$$

$$
\tan \Theta\_0 \sim \chi\_0 / \mathbf{z}\_0 \text{ or } \mathbf{y} \llcorner \mathbf{z}\_0 \tag{11}
$$

where Δ*B* is the frequency bandwidth of the measurement, λ<sup>0</sup> is the central frequency wavelength, θ<sup>0</sup> is the target angle visualization seen from the measurement distance, and *Z*<sup>0</sup> and *X*0, are the scanning distances along the *x*-axis and *y*-axis, respectively.

Differential imaging is obtained by means of calculating the difference of the experimental data (total image obtained from focused field amplitude) measured using the car, including with a monopole antenna short-circuited (ON state) and without the antenna (OFF state or open-circuited antenna). As stated in Section 2, the short-circuited antenna represents the antenna placed on the top of the car, which is in contact with the car's metal surface.

The objective is to obtain the currents created by the antenna over the surface of the vehicle without getting dazzled ("flash-out") by the higher intensity of the currents on the antenna itself. A differential image may be produced due to the difference in module between the focused field amplitude images with the antenna (see Figures 6 and 7) and without the antenna (see Figures 4 and 5), resulting in the current distribution of the antenna that is shown in Figures 8 and 9.

**Figure 8.** Differential focused field amplitude between the metallic car with and without an antenna. Frequency range: 220 GHz to 330 GHz (BW = 110 GHz).

**Figure 9.** Differential focused field amplitude between the metallic car with and without an antenna. Frequency range: 295 GHz to 305 GHz (BW = 10 GHz).

In Figure 8, it can be identified that the strong central spot corresponding to the antenna has disappeared and, instead, the nearby currents to the antenna are visualized. The maximum can be clearly identified.

Figure 9 depicts the current distribution calculated as a differential focused field amplitude using a frequency range from 295 GHz to 305 GHz (BW = 10 GHz) to assess influence of bandwidth in differential imaging.

The differential focused field amplitude is higher with a bandwidth of 10 GHz around 300 GHz compared with a BW of 110 GHz from 220 GHz to 330 GHz. This effect is due to the gain variation of the horns across the frequency range. The imaging system performs an averaging through all the bandwidths, and the horn antenna gain variation generates this change in the amplitude, with higher amplitude when the BW is lower. However, the image resolution is better when the BW is greater [3]. If the application requires high resolution, then the BW parameter should be maximized.

#### **5. Conclusions**

In this work, an application of the approach in the near-field exploration of metallic objects explained in [3] is used to perform the focusing of a metallic car with and without a monopole antenna in order to estimate the induced currents on the vehicle surface using a differential methodology. This method makes it possible to estimate the distribution of the currents induced by the antenna on the vehicle's surface using a very simple method, which is also combined with the frequency-dimension scale translation.

The technique presented in [3] has been extended to a different measurement scenario in a 3D scan in the terahertz band and by using the differential current concept as both a novel form and a major contribution.

The findings verify that, while distance discrimination increased with the bandwidth, the transverse resolution (*XY* plane) did not change significantly.

The results show promising performance from this technique using differential currents, as depicted in the 300 GHz band. Gaining knowledge of these vehicle-antenna-created currents may offer hints regarding the impact produced by the shapes (borders and corners), discontinuities (slots and windows), or materials (composition or roughness) introduced into the vehicular platform.

**Author Contributions:** Conceptualization, J.A.S.-P., J.-M.M.-G.-P., J.R., and L.J.-R.; Methodology, J.A.S.-P. and J.-M.M.-G.-P.; Software, J.A.S.-P. and J.-M.M.-G.-P.; Validation, J.A.S.-P., J.-M.M.-G.-P., J.R., and L.J.-R.; Resources, J.-M.M.-G.-P. and M.-T.M.-I.; Data curation, M.-T.M.-I.; Writing—Original Draft Preparation, J.A.S.-P.; Writing—Review and Editing, J.A.S.-P., J.-M.M.-G.-P., J.R., L.J.-R., C.B.-S., and J.-V.R.; Supervision, J.-M.M.-G.-P., J.R., and L.J.-R., M.-T.M.-I., J.-V.R., and A.M.-A.; Project Administration, J.-M.M.-G.-P. and J.-V.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Ministry of Education and Science, Spain (PID2019-107885GB-C33/ AEI/10.13039/501100011033, TEC2016-78028-C3-1-P, TEC2016-78028-C3-2-P and TEC2016-78028-C3-3-P, MDM2016- O6OO), Catalan Research Group 2017 SGR 219, and the European FEDER funds.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Sensors* Editorial Office E-mail: sensors@mdpi.com www.mdpi.com/journal/sensors

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com

ISBN 978-3-0365-4204-1