**Through-the-Wall Microwave Imaging: Forward and Inverse Scattering Modeling**

**Alessandro Fedeli 1, Matteo Pastorino 1, Cristina Ponti 2,3,\* Andrea Randazzo 1,\* and Giuseppe Schettini 2,3**


Received: 18 April 2020; Accepted: 16 May 2020; Published: 18 May 2020

**Abstract:** The imaging of dielectric targets hidden behind a wall is addressed in this paper. An analytical solver for a fast and accurate computation of the forward scattered field by the targets is proposed, which takes into account all the interactions of the electromagnetic field with the interfaces of the wall. Furthermore, an inversion procedure able to address the full underlying non-linear inverse scattering problem is introduced. This technique exploits a regularizing scheme in Lebesgue spaces in order to reconstruct an image of the hidden targets. Preliminary numerical results are provided in order to initially assess the capabilities of the developed solvers.

**Keywords:** electromagnetic scattering; buried objects; through-wall radar; microwave imaging; inverse scattering

#### **1. Introduction**

Microwave imaging of targets placed behind a wall is a topic of great interest in the remote detection of humans or objects in indoor environments, with applications in surveillance, rescue, and defense [1]. The instrument usually adopted in this class of surveys is the so-called through-wall (TW) radar, which includes a wide set of possible hardware architectures. For example, as regards the source of the interrogating field, a possible solution is to use wideband antennas radiating pulsed electromagnetic (EM) fields, or to adopt frequency-stepped sources, radiating an EM field at a limited number of frequencies.

Beyond the radar architecture, the processing of experimental data plays an important role. In most processing approaches, the goal is to localize hidden targets or produce a qualitative image of the indoor scenario [1–6], returning information about the target's shape or the presence of interfaces. Accurate forward scattering approaches may improve the reconstruction capabilities of radar surveys. First, the synthetic data obtained by the numerical modeling are helpful in providing a deeper physical insight on the fields scattered by typical TW targets, in several practical cases. Second, the theoretical solution to the forward scattering problem is a useful tool in imaging approaches, especially when aiming for a quantitative reconstruction of the target. In this respect, forward solvers find twofold applications. On the one hand, the numerical data they provide can be used as input data to validate inversion schemes. On the other hand, forward solvers can be employed as building blocks of non-linear inversion algorithms themselves. In this perspective, the forward scattering solver should develop a full-wave solution, where all the effects of the wall interfaces on the propagation of the scattered field are suitably modeled. As for the forward solvers in TW problems, the methods proposed in the literature are essentially numerical [2,3,7–10]. Due to its high flexibility, the finite-difference time-domain (FDTD)

method is mainly adopted, and, being directly developed in the time-domain, its application to TW radars with pulsed sources is straightforward [7–9,11]. When modeling stepped-frequency sources, frequency domain data should be computed through an inverse fast Fourier transform [2]. However, due to the large size of investigation domains in TW settings, the FDTD modeling for frequency domain analysis may be demanding in terms of execution times and memory requirements. As for the frequency-domain techniques, other methods such as the ones based on method-of-moments (MoM) approaches may be employed, which are still numerical [3]. Asymptotic techniques are also used for modeling very large regions [10]. However, when applied to imaging approaches, forward solvers are usually implemented through linearized formulations [2,3,12–14] or by using synthetic aperture schemes [4,5,15], thus leading to qualitative images of the target. Techniques to improve the spatial resolution [16,17] and human discrimination [18,19] have also been proposed, through frequency spectral analysis. The implementation of full-wave inverse scattering approaches for quantitative TW imaging is still an open issue, as most algorithms for solving the non-linear inverse scattering problem are usually developed for free-space applications [20–30] and must be properly adapted to include the presence of the wall.

In this paper, the issues relevant to the modeling of the forward scattering problem and to the quantitative imaging are both addressed. In particular, a single-frequency non-linear inversion procedure based on a regularization scheme in Lebesgue-spaces is proposed. This kind of inversion strategy was initially developed in free-space environments [31–35] and subsequently extended to the case of targets buried in a homogeneous soil [36], showing good reconstruction capabilities. It has been found that the geometrical properties of Lebesgue-space norms lead to regularized solutions endowed with less oversmoothing than classical Hilbert-space regularization schemes, improving both target localization and their shaping. In this work, the method is extended to the case of targets hidden behind a dielectric wall, by modifying the underlying scattering model in order to include the proper Green's function for layered media [7]. The validity of this inversion scheme is assessed in its forward scattering formulation, as well as in its application to the imaging procedure. For the validation of the proposed technique, an analytical solver has been employed, called the cylindrical wave approach (CWA) [37–41]. It provides an analytical/numerical solution to a layout with buried targets, given by circular cross-section cylinders, employing cylindrical waves as basis functions of the scattered field. Due to the presence of the interface, suitable cylindrical waves are introduced, defined through spectral integrals, to deal with reflection and transmission of the scattered field by the interfaces. A preliminary implementation of the CWA for TW scattering was proposed in [38], with an approach through multiple reflection fields to accurately describe the multipath inside the wall. In [39], the CWA has been developed through a non-iterative approach, where all the multiple interactions were included in suitable reflection and transmission coefficients. The same technique has been applied in [41], where a pulsed source field has been modeled. In this work, the results provided by the analytical solver are validated considering an integral formulation based on a Green's function approach, and they are used as reference data to test the reconstruction capabilities of the inversion procedure. It is worth remarking that the novelty of the paper is both on inverse and forward modeling. From the point of view of the inversion procedure, an efficient technique working in the framework of Lebesgue spaces (which was developed for free-space and half-space scenarios) is here extended to through-the-wall configurations. Specifically, the presence of the wall is considered by inserting the proper Green's function into the scattering model. In this way, the through-the-wall propagation phenomena are fully taken into account, even under near-field conditions (differently from synthetic aperture and beamforming schemes, where far-field conditions are usually assumed). Moreover, the adopted scattering model does not rely on approximations (e.g., the Born or Kirchhoff models often adopted in TW imaging). However, the presence of multiple interfaces and the availability of few limited-view measurements (i.e., only along one side of the wall) significantly increase the difficulty of the inverse problem. Consequently, the present paper is aimed at evaluating the regularization properties of the approach even in this more involved case. Moreover, an automatic criterion for

selecting the optimal Lebesgue-space norm parameter based on the entropy principle is proposed for the first time. The analytical solver used in the validation of the inversion procedure is implemented for a monochromatic line-current source and dielectric targets, applying the spectral-domain analysis developed in [39] to a more realistic source. In the proposed approach, the total field is decomposed into two different sets: scattered fields by the cylinders, and non-scattered fields, i.e., the field radiated by the line source and the fields excited by its reflection and transmission at the interfaces, in the absence of the cylinders. The non-iterative approach is applied to both sets and, through suitable reflection and transmission coefficients, all the multiple interactions of the fields inside the layer are collected in two contributions: an up-ward and a down-ward propagating wave. Therefore, a theoretical solution is developed through a very compact formulation, with the total field in each medium decomposed in a limited number of terms. This approach leads to a numerical implementation which is fast and efficient.

The paper is organized as follows: in Section 2, an overview of the proposed forward and inverse scattering approaches is reported. Forward and inverse scattering results are then presented in Section 3. Conclusions follow in Section 4.

#### **2. Theoretical Approach to the Through-Wall Imaging Problem** 7

The geometry of the TW imaging problem is shown in Figure 1. A two-dimensional layout is considered, with one lossless dielectric wall between two semi-infinite regions filled with air (i.e., characterized by the vacuum dielectric permittivity, ε0). The wall has relative permittivity ε*r*<sup>1</sup> and thickness *l*. The hidden investigation domain *Dinv*, highlighted by the dashed box in Figure 1, is located in the medium behind the wall, and contains one infinitely long cylinder with circular cross section having center in (*xc*, *yc*), radius *a*, and relative permittivity ε*rc*.

**Figure 1.** Configuration of the scattering problem, with one target placed behind a dielectric wall.

A set of *M* transmitting/receiving antennas placed along a line of length *Ls* parallel to the interface at a fixed distance *y* = *hs* in front of the wall is used. The transmitting antennas are modeled by monochromatic line-current sources with angular frequency ω, and it is assumed that a TM*z*-polarized incident electric field **E***inc* = *Einc*(*x*, *y*) **^ z** is excited. The expression of the field radiated by the transmitting antennas with center in (*ds*, *hs*) is given by [42]:

$$E\_{\rm inc}(\mathbf{x}, \mathbf{y}) = V\_0 H\_0^{(2)} \left( \sqrt{\left(\mathbf{x} - d\_s\right)^2 + \left(\mathbf{y} - h\_s\right)^2} \right) \tag{1}$$

where *<sup>H</sup>*(2) <sup>0</sup> (·) is the zero-th order second-kind Hankel function, and *V*<sup>0</sup> is the complex amplitude of the field. The *ej*ω*<sup>t</sup>* term is omitted throughout the paper.

Antennas are scanned in a multi-illumination multi-view configuration, i.e., each antenna is used in turn in transmission mode to radiate the incident electric field in (1), and the total field **E***tot* = *Etot*(*x*, *y*) **^ z** produced by the interaction of the EM wave with the investigation domain (including the wall and the target) is received by the remaining *M* − 1 antennas. It is worth remarking that the assumed scattering model is formally exact only when dealing with cylindrical targets (i.e., ideally infinite and invariant along the *z* direction) under a TM*<sup>z</sup>* illumination. In practical TW imaging applications, the inspected objects, as well as the wall, although not being infinite are usually elongated along the vertical direction (corresponding to the *z* axis in our settings). Consequently, the predicted fields are sufficiently accurate for solving the imaging problem at hand.

#### *2.1. Forward-Scattering Problem Formulation*

The theoretical method adopted to evaluate the scattered field in the layout of Figure 1 is the cylindrical wave approach [39,41]. The total field **E***tot* = *Etot*(*x*, *y*) **^ z** is given by the superposition of two sets of fields. The field radiated by the transmitting antenna in (1) and the fields relevant to its reflection and transmission from the interface (in the absence of the target) belong to the first set, representing known field contributions. The second set of fields is given by the scattered field by the target in the medium behind the wall, and by the scattered-reflected and transmitted fields through the wall interfaces. In the lowest medium, the scattered electric field is found from *Escatt*(*x*, *y*) = *Etot*(*x*, *y*) − *Et*2(*x*, *y*), where *Et*2(*x*, *y*) is the field related to the transmission of the incident field *Einc*(*x*, *y*) in the medium behind the wall. The scattered field in the medium behind the wall is given in turn by the superposition of three contributions, *Es*(*x*, *y*), *Esr*(*x*, *y*), *Esc*(*x*, *y*), i.e,

$$E\_{\text{scatt}}(\mathbf{x}, \mathbf{y}) = E\_{\text{s}}(\mathbf{x}, \mathbf{y}) + E\_{\text{s}\mathbf{r}}(\mathbf{x}, \mathbf{y}) + E\_{\text{sc}}(\mathbf{x}, \mathbf{y}) \tag{2}$$

where *Es* represents the field scattered by the target, *Esr* is the scattered-reflected field, describing the reflection of the field *Es* by the wall, and *Esc* is the contribution of scattered field that is transmitted inside the cylinder.

The scattered field *Es* in (2) is expressed through an expansion into a series of basis functions *CWm* [37]:

$$E\_s(\mathbf{x}, \mathbf{y}) = V\_0 \sum\_{m = -\infty}^{+\infty} c\_m \mathbf{C} \mathcal{W}\_m(\mathbf{x}, \mathbf{y}) \tag{3}$$

where *cm* are unknown expansion coefficients and the basis functions *CWm* are cylindrical waves, proportional to *m*-th order Hankel functions:

$$\text{CW}\_{\text{m}}(\mathbf{x}, \mathbf{y}) = H\_{\text{m}}^{(2)}(k\_0 r) e^{j m \theta} \tag{4}$$

where (*r*, θ) are polar coordinates centered on the cylinder.

The use of cylindrical waves as functions of expansion of the fields scattered by circular cross-section cylinders gives the analytical basis to the method. However, as the target is not in free space, but placed behind a dielectric wall, the interaction with the wall interfaces in terms of reflection and transmission must be suitably modeled. This is accomplished by expressing the cylindrical waves in (4) through an alternative definition, i.e., the plane-wave spectrum of a cylindrical wave:

$$\text{CW}\_{m}(x,y) = \frac{1}{2\pi} \int\_{-\infty}^{+\infty} F\_{m}(y,k\_{\parallel})e^{-jk\_{\parallel}x}dk\_{\parallel} \tag{5}$$

where *Fm <sup>y</sup>*, *<sup>k</sup>*|| is the plane-wave spectrum:

$$F\_m(y, k\_{\parallel}) = -\frac{2e^{-j|y|}\sqrt{1 - \left(k\_{\parallel}\right)^2}}{\sqrt{1 - \left(k\_{\parallel}\right)^2}} \begin{cases} \left.e^{jm\cos^{-1}n\_{\parallel}}\right|\_{\mathcal{I}} & y \ge 0\\ \left.e^{-jm\cos^{-1}n\_{\parallel}}\right|\_{\mathcal{I}} & y \le 0 \end{cases} \tag{6}$$

The expressions (5) and (6) are used to derive the scattered-reflected field *Esr*(*x*, *y*) in (2) and the scattered fields propagating inside the wall and in the first half-space [39]. In particular, in the half-space in front of the wall, where the field is probed by the receiving antennas, the scattered field is found as *Escatt*(*x*, *y*) = *Etot*(*x*, *y*) − *Einc*(*x*, *y*) − *Er*1(*x*, *y*), where *Er*1(*x*, *y*) is the contribution relevant to the reflection of incident field by the interface in *y* = 0. The scattered field *Escatt*(*x*, *y*) is defined through the following expansion [39]:

$$E\_{\rm scatt}(\mathbf{x}, \mathbf{y}) = V\_0 \sum\_{m = -\infty}^{+\infty} c\_m T \mathcal{W}\_m^0(\mathbf{x}, \mathbf{y}; \mathbf{y}\_\varepsilon) \tag{7}$$

where the basis functions *TW*<sup>0</sup> *<sup>m</sup>*(*x*, *y*; *yc*) are transmitted cylindrical waves, and they are expressed through spectral integrals:

$$T\mathcal{W}\_{m}^{0}(\mathbf{x},\mathbf{y};\boldsymbol{y}\_{\varepsilon}) = \frac{1}{2\pi} \int\_{-\infty}^{+\infty} T\_{10}(k\_{\parallel}) T\_{21}(k\_{\parallel}) F\_{m} \left[ n\_{2}(\boldsymbol{y}\_{\varepsilon} - \boldsymbol{I}), k\_{\parallel} \right] e^{-j\boldsymbol{y}\_{\varepsilon}\sqrt{k\_{0}^{2} - \left(n\_{2}k\_{\parallel}\right)^{2}}} e^{-j\boldsymbol{k}\_{\parallel} \left(\mathbf{x} - \mathbf{x}\_{\varepsilon}\right)} dk\_{\parallel} \tag{8}$$

In (8), *T*10 *n*|| and *T*21 *n*|| are the transmission coefficients from the wall to the upper medium and from the lowest medium to the wall, respectively. In the expression (8), all the multiple reflections excited by propagation of the scattered field *Es* in the wall are included through transmission and reflection coefficients related to the interaction of a plane wave with a dielectric slab [43]. A solution to the scattering problem is developed imposing the boundary conditions of continuity of the field components tangential to the cylinder's interface and deriving the unknown expansion coefficients *cm* in (3) and (8) [39].

#### *2.2. Inverse-Scattering Problem Formulation*

In the inversion procedure, the space-dependent dielectric properties of the investigation domain *Dinv* are described by the contrast function *c*(*x*, *y*) = ε(*x*, *y*)/ε<sup>0</sup> − 1, ε(*x*, *y*) being the dielectric permittivity in a generic point (*x*, *y*) ∈ *Dinv*, which represents the unknown to be retrieved. Such a quantity is related to the scattered field *Escatt* in the measurement points by means of the following integral relationship (data equation) [21]

$$E\_{\rm scatt}(\mathbf{x}, y) = \mathcal{G}\_{\rm uv}^{\rm ext}(cE\_{\rm lat})(\mathbf{x}, y) = -k\_0^2 \int\_{D\_{\rm inv}} c(\mathbf{x'}, y') \mathcal{E}\_{\rm lat}(\mathbf{x'}, y') \mathcal{g}\_{\rm uv}(\mathbf{x}, y, \mathbf{x'}, y') d\mathbf{x'} dy' \tag{9}$$

where *k*<sup>0</sup> = ω(ε0μ0) 0.5 is the vacuum wavenumber and <sup>G</sup>*ext <sup>w</sup>* is a linear integral operator whose kernel is the two-dimensional Green's function of the considered three-layer background, *gw*, which is given by [7]

$$g\_{w}(\mathbf{x}, y, \mathbf{x}', y') = \frac{j}{4\pi} \int\_{-\infty}^{+\infty} \frac{e^{j\zeta(\mathbf{x} - \mathbf{x}')}}{\gamma\_{1}} \begin{cases} e^{-j\gamma\_{0}y - y')} + Re^{-j\gamma\_{0}y \cdot (y + y')}, & y \ge 0\\ Te^{j\gamma\_{0}(y + l - y')}, & y \le -l\_{w} \end{cases} d\zeta \tag{10}$$

where γ<sup>0</sup> = *k*2 <sup>0</sup> <sup>−</sup> <sup>ζ</sup><sup>2</sup> and

$$R = \rho\_w \frac{1 - e^{-2j\gamma\_1 l}}{1 - \rho\_w^2 e^{-2j\gamma\_1 l}}, \quad T = \frac{(1 - \rho\_w^2)e^{-j\gamma\_1 l}}{1 - \rho\_w^2 e^{-2j\gamma\_1 l}}\tag{11}$$

with γ<sup>1</sup> = *k*2 <sup>1</sup> <sup>−</sup> <sup>ζ</sup>, *<sup>k</sup>*<sup>1</sup> <sup>=</sup> *<sup>k</sup>*0ε0.5 *<sup>r</sup>*<sup>1</sup> being the wavenumber in the wall, and ρ*<sup>w</sup>* = (γ<sup>0</sup> − γ1)/(γ<sup>0</sup> + γ1). For the sake of simplicity, a single view case is considered in this Section. The total electric field *Etot*(*x* , *y* ) inside the integral in (9) depends itself on the contrast function *c* and can be expressed

by means of a second integral equation similar to (9) (the so-called state equation), i.e., *Etot*(*x*, *y*) = *Einc*(*x*, *y*) + G*int <sup>w</sup>* (*cEtot*)(*x*, *y*), where G*int <sup>w</sup>* is again a linear integral operator whose kernel is the Green's function for the through-wall configuration [21]. By combining the data and state equations, the inverse scattering problem can be finally formulated as [21]

$$E\_{\rm scatt}(\mathbf{x}, \mathbf{y}) = \mathcal{T}(\mathbf{c})(\mathbf{x}, \mathbf{y}) = \mathcal{G}\_{\rm uv}^{\rm ext} \mathbf{c} \{ \mathbf{I} - \mathcal{G}\_{\rm uv}^{\rm int} \mathbf{c} \} \mathbf{E}\_{\rm inc}(\mathbf{x}, \mathbf{y}) \tag{12}$$

The non-linear problem at hand is solved in a regularized sense by using an inversion procedure developed in the framework of Lebesgue spaces *Lp*, i.e., function spaces endowed with the norm *u p <sup>L</sup><sup>p</sup>* <sup>=</sup> *u*(*x*, *<sup>y</sup>*) *p dxdy* (u being a generic function belonging to *Lp*). It is worth noting that the norm exponent *p* represents an additional parameter that can be tuned in order to enhance the reconstruction performance. In particular, the developed procedure is based on an iterative outer–inner Newton scheme, which can be summarized by the following steps [31,33]:


$$\mathcal{L}\_{n,l+1} = \mathcal{J}\_q \Big( \mathcal{J}\_p(\boldsymbol{\xi}\_{n,l}) - \beta \mathcal{T}\_n^{'\*} \mathcal{J}\_p(\mathcal{T}\_n^{'} \boldsymbol{\xi}\_{n,l} - \mathcal{E}\_{\text{scatt}}(\mathbf{x}, y) + \mathcal{T}(\boldsymbol{c}\_n)) \Big) \tag{13}$$

where ξ*n*,0 = 0, β = T *n* −2 <sup>2</sup> is a relaxation coefficient, *q* = *p*/(*p* − 1) is the Hölder conjugate of *p*, and the duality map J*<sup>p</sup>* is defined as J*p*(*e*) = *e* 2−*p <sup>p</sup>* |*e*| *<sup>p</sup>*−1sign(*e*), with sign(*e*)<sup>=</sup> *<sup>e</sup>*/|*e*<sup>|</sup> (if *<sup>e</sup>* - 0, otherwise it is equal to zero).


#### **3. Numerical Results**

#### *3.1. Validation of the Forward Methods*

A comparison between the analytical TW solver (Section 2.1) and the forward scattering model embedded inside the inversion procedure (Section 2.2) is reported here, for a cross-validation of the two forward approaches. A multistatic and multiview configuration has been simulated, with *M* = 15 transmitting and receiving antennas aligned in front of the wall along a line of length *Ls* = 1.5 m, with spacing *d* = *Ls*/(*M* − 1), and parallel to the wall at distance *hs* = 30 cm. The *s*-th transmitting antenna (*s* = 1, ... , *M*) is placed along the horizontal axis in the following position:

$$\mathbf{x}\_s^{TX} = -\frac{L\_s}{2} + (s-1)d\tag{14}$$

whereas the scattered field is probed at the remaining *M* − 1 positions along *Ls*. The working frequency has been fixed equal to 1 GHz. As a first case, a single dielectric cylinder with center in (−20 cm, 60 cm), radius *a* = 10 cm, and relative permittivity ε*rc* = 2, placed behind a wall of relative permittivity ε*r*<sup>1</sup> = 4 and thickness *l* = 20 cm, has been considered. The actual distribution of the relative dielectric

permittivity in the investigation domain is shown in Figure 2a. In the MoM solver, the target has been discretized into *N* = 900 square subdomains of side about equal to 6.7 mm. In the CWA, the order *m* of the cylindrical waves in Equation (7) has been truncated to *Mt* = 9, being the total number of terms in the cylindrical expansions equal to 2*Mt* + 1. The truncation order has been determined applying the rule *Mt* = 3 <sup>√</sup>*<sup>r</sup>*1*a*(2π)/λ, that allows a compromise between accuracy and computational heaviness. Figure 3 shows the amplitude and phase of the fields computed by the two approaches for some of the considered views, at the *M* − 1 measurement receiving points. Plots are evaluated for different values of the index *s*, which denotes the antenna used in the transmission mode, according to Equation (14). As can be seen, a very good agreement between the analytical solver used in the forward approach and the solver employed in the inversion procedure is obtained.

**Figure 2.** Actual configuration. (**a**) Single dielectric cylinder and (**b**) two separate dielectric cylinders.

**Figure 3.** (**a**) Amplitude and (**b**) phase of the scattered fields in some of the considered view computed by the analytical forward solver based on the cylindrical wave approach (CWA) and by the integral equation formulation used in the inverse scattering procedure. Single dielectric cylinder.

As a second test case, two dielectric cylinders have been considered inside the investigation domain. The first one is the same considered above, whereas the second one is a dielectric cylinder with center in (20 cm, 60 cm), radius *a* = 10 cm, and relative permittivity ε*rc* = 2. The corresponding distribution of the relative dielectric permittivity in the investigation domain is shown in Figure 2b. Figure 4 reports the amplitude and phase of the field computed by using the CWA and the MoM approaches. In this more complex case, too, there is a good agreement between the two solving schemes, confirming the correctness and suitability of the analytical and numerical solvers adopted in the data generation and inversion steps.

**Figure 4.** (**a**) Amplitude and (**b**) phase of the scattered fields in some of the considered view computed by the analytical forward solver based on the CWA and by the integral equation formulation used in the inverse scattering procedure. Two separate dielectric cylinders.

#### *3.2. Inversion Scheme*

Some preliminary examples of reconstructions provided by the previously described inversion procedure are reported in this Section. In particular, the two configurations adopted for the comparison in the previous Section are considered. In order to simulate a more realistic scenario, a Gaussian noise with zero mean value and variance corresponding to a signal-to-noise ratio of 20 dB has been added to the computed scattered electric field data. The following values of the algorithm's parameters have been used: *p* ∈ [1.1, 2.5]; maximum number of outer iterations, 10; maximum number of inner iterations, 50; iterations stopped when the relative variation of the residual falls below the threshold 0.005. Such values have been empirically selected, based on the previous experience on Lebesgue-space inversion in the free-space scenario. The optimal value of the norm parameter *p* has been found by performing a sweep in the assumed range of values and by selecting the one providing the maximum entropy. Such a choice has been made since, for this particular application, it is expected that localized targets are usually present inside the inspected scenario, and consequently the sharpness of the image, which is related to its entropy, may represent a discriminating feature [44,45].

Figure 5 shows the reconstructed distribution of the relative dielectric permittivity retrieved by the developed inverse scattering procedure. In particular, the reconstruction obtained with the optimal value of the norm parameter, i.e., *p* = *popt* = 1.3, is reported in Figure 5a. This result evidences a correct localization of the target. Indeed, the estimated center of the cylinder is (−19.9 cm,−61.4 cm), which corresponds to an average percentage error of 1.3%. Moreover, the reconstructed value of the dielectric permittivity is close to the actual one. Specifically, the maximum value of the estimated permittivity is 1.81, which compares very well with the actual value of 2. Nevertheless, the target shape, which is a circular one, is elongated along the range direction and the cross-range size is underestimated. Such a behavior can be ascribed to the use of data collected at a single frequency, as well as to their aspect-limitedness. However, it is worth noting that even with such a small number of available data and considering just a single working frequency, the approach is able to effectively provide a quite accurate indication about the target. For comparison purposes, the reconstruction obtained by using a standard inversion procedure in Hilbert spaces (corresponding to *p* = 2) is provided in Figure 5b. In this case, the target is still visible (the average percentage error on the center position is 2.9%), but the dielectric permittivity is significantly underestimated (the maximum value is 1.35). Moreover, stronger artifacts are present in the background.

**Figure 5.** Reconstructed distribution of the relative dielectric permittivity inside the through-wall (TW) investigation domain. Single dielectric cylinder. (**a**) Optimal value of the norm parameter (*popt* = 1.3) and (**b**) standard Hilbert-space approach (*p* = 2 ).

As a second test case, the scattering data from the two dielectric cylinders considered in the previous Section have been used. In this case, too, the scattered field has been corrupted with a Gaussian noise with zero mean value and variance corresponding to *SNR* = 20 dB. The parameters of the inversion procedure are the same as in the previous case. The reconstructed distribution of the relative dielectric permittivity is shown in Figure 6. In particular, Figure 6a shows the results obtained by considering the optimal value of the norm parameter, which is equal to *popt* = 1.3 in this case. Even in this situation, the two targets are correctly localized, although the dielectric permittivity is slightly underestimated. Indeed, the estimated centers of the cylinders are (−19.4 cm,−57.4 cm) and (−19.0 cm,−57.0 cm), which correspond to the mean percentage errors of 3.7% and 5%, respectively, whereas the maximum values of the dielectric permittivity are both equal to 1.7. The corresponding reconstruction obtained by using the standard Hilbert-space reconstruction technique is shown in Figure 6b. Similarly to the preceding configuration, the two targets are visible (the mean percentage errors on the position are 11.2% and 7.3%), although their properties are strongly underestimated and with amplitude comparable to the background artifacts (the maximum value of the dielectric permittivity is equal to 1.26). The criterion for the selection of the optimal reconstruction has been assessed in this case by comparing the behavior of the scaled entropy with the one of the reconstruction errors (defined as NMSE = *c* − *cact* 2/*cact* 2, *cact* being the actual distribution of the contrast function). As can be seen from Figure 7, the scaled entropy (defined as in [45]) has a maximum corresponding to the value of the norm parameter for which the lowest reconstruction error is obtained.

**Figure 6.** Reconstructed distribution of the relative dielectric permittivity inside the TW investigation domain. Two dielectric cylinders. (**a**) Optimal value of the norm parameter (*popt* = 1.3) and (**b**) standard Hilbert-space approach (*p* = 2 ).

**Figure 7.** Behavior of the scaled entropy and of the reconstruction error versus the norm parameter.

#### **4. Conclusions**

In this work, the forward scattering problem modeling and the quantitative imaging in through-wall scenarios have been addressed. For the solution of the forward EM problem, an analytical solver based on the cylindrical wave approach has been presented. The inverse problem is solved by using a non-linear regularization technique developed in the framework of Lebesgue-spaces. The suitability of the forward and inverse solvers for the problem at hand has been evaluated with preliminary numerical simulations. Future works will be mainly devoted to the extension to multi-frequency processing, in order to increase the reconstruction accuracy, and to three-dimensional configurations. Moreover, the developed forward and inverse schemes will be validated by considering experimental measurements obtained with a real hardware setup.

**Author Contributions:** Conceptualization, A.F., M.P., C.P., A.R., and G.S.; formal analysis A.F., M.P., C.P., A.R., and G.S.; methodology A.F., M.P., C.P., A.R., and G.S.; validation, A.F., M.P., C.P., A.R., and G.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partly supported by the Italian Ministry for Education, University, and Research under the project PRIN2015 U-VIEW, grant number 20152HWRSL.

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

#### **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/).

## *Article* **Experimental Validation of a Microwave Imaging Method for Shallow Buried Target Detection by Under-Sampled Data and a Non-Cooperative Source**

**Adriana Brancaccio 1,2,\*, Giovanni Leone 1,2, Rocco Pierri 1,2 and Raffaele Solimene 1,2,3**


**Abstract:** In microwave imaging, it is often of interest to inspect electrically large spatial regions. In these cases, data must be collected over a great deal of measurement points which entails long measurement time and/or costly, and often unfeasible, measurement configurations. In order to counteract such drawbacks, we have recently introduced a microwave imaging algorithm that looks for the scattering targets in terms of equivalent surface currents supported over a given reference plane. While this method is suited to detect shallowly buried targets, it allows one to independently process all frequency data, and hence the source and the receivers do not need to be synchronized. Moreover, spatial data can be reduced to a large extent, without any aliasing artifacts, by properly combining single-frequency reconstructions. In this paper, we validate such an approach by experimental measurements. In particular, the experimental test site consists of a sand box in open air where metallic plate targets are shallowly buried a (few *cm*) under the air/soil interface. The investigated region is illuminated by a fixed transmitting horn antenna, whereas the scattered field is collected over a planar measurement aperture at a fixed height from the air-sand interface. The transmitter and the receiver share only the working frequency information. Experimental results confirm the feasibility of the method.

**Keywords:** radar imaging; target detection; experimental measurements; microwave imaging

#### **1. Introduction**

Microwave imaging, and in general radar imaging, is a mature research field that finds application in a number of different contexts where pursuing non-destructive investigation is convenient or mandatory [1–11].

Ground Penetrating Radar (GPR) is a radar system that is properly conceived to address non-destructive imaging. Generally, GPRs work in contact with the interface between the air and the medium under investigation. However, there is great interest in achieving target detection through non-contact measurement layouts, for example, with GPRs mounted on a flying platform [12,13]. Indeed, stand-off distance configurations allow for the investigation of regions that are not easily (or safely) accessible, as it happens, for instance, when one has to deal with mine or unexploded device detection [14]. Moreover, a flying GPR can allow for inspecting large areas quickly [15,16].

In this framework, however, the system cost and the achievable performance must be traded-off [17]. This requires finding a compromise between the time needed to collect data, the number of sensors to be simultaneously deployed, and the way transmitter and receivers "communicate". In this regard, a single-view/multistatic configuration seems convenient, since only one sensor transmits and the others act as mere passive receivers,

**Citation:** Brancaccio, A.; Leone, G.; Pierri, R.; Solimene, R. Experimental Validation of a Microwave Imaging Method for Shallow Buried Target Detection by Under-Sampled Data and a Non-Cooperative Source. *Sensors* **2021**, *21*, 5148. https:// doi.org/10.3390/s21155148

Academic Editor: Giacomo Oliveri

Received: 21 June 2021 Accepted: 27 July 2021 Published: 29 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/).

hence with reduced weight and cost. However, synchronization between TX and RXs is a critical issue when they are mounted on different platforms, since, differently from multi-monostatic arrangement, TX and RXs are no more co-located and hence do not share the same electronic system. In addition, the number of measurement points is directly linked to the number of flying platforms and hence must be reduced as much as possible.

Actually, a number of different processing schemes have been proposed in the literature to address subsurface imaging [18]. Recently, we have introduced a reconstruction scheme that allows one to mitigate the previously mentioned drawbacks [19,20]. This method relies on a certain scattering formulation where, according to the equivalence principle [21], the scattered field is modeled as being radiated by equivalent surface currents that are supported over the air/soil interface, or at some other reference plane whose depth is chosen. A related method that uses equivalence principles for target shape reconstruction is reported in [22], where, however, the multiple experiments arise from using different illuminations (i.e., a multistatic/multiview configuration is employed).

In our approach, basically, the main idea is that if the reference plane is close to the scattering target, then the spatial support of the surface current gives information concerning the transverse location of the target. If it is a priori known that the targets are in close proximity to the air/soil interface, e.g., for detection of a mine [23] or unexploded improvised device [24] (where the targets are often very shallowly buried just to hide from sight), then the reference plane can be set just at the air/soil interface. In this case, the surface currents radiate in free-space, and the related simple Green function can be considered as a propagator. The reconstruction is cast as a 2D inverse problems, since only the targets detection and their transverse locations are looked for, and hence single-frequency data can be employed. Accordingly, RXs do not need information about the TX, except the working frequency. Hence, the source can be considered as being non-cooperative (it does not share information with the receivers). However, it is not opportunistic as in most passive radar, since it is deliberately deployed in the scene. Multi-frequency data can be processed separately (i.e., incoherently) and then combined to counteract aliasing artifacts that can arise if the spatial measurement points are reduced (not properly sampled) [20]. Finally, depth can be explored by performing the reconstructions at different reference planes.

In previous contributions (see [19,20]), we have shown the feasibility of the method by employing synthetic numerical data and some measured data collected under lab conditions for a free-space scattering scenario. However, the method still need to be validated in a realistic scattering scenario.

In this contribution, we aim at pursuing such a task. To this end, the method and the related achievable performances are checked for a more realistic scattering scenario where the background medium is actually not homogeneous. Indeed, the test site mimics a realistic on-field situation, as it consists of an open-air sand box. As to the RXs and the Txs, they are not mounted on flying platforms. However, they are at a stand-off distance from the air/soil interface, and the TX signal is unknown to the RXs but the working frequency is known.

In sum, the advancements that we are conveying in this contribution concern:


The rest of the paper is organized as follows. In Section 2 the scattering model is generalized to deal with the new scattering configuration, whereas in Section 3, the related reconstruction algorithm is briefly described. In Section 4, the experimental set-up as well as a few experimental results are presented and discussed. Conclusions end the paper.

#### **2. Scattering Model**

In this section, we introduce the adopted scattering model and the necessary notation that is required in the following. A detailed derivation of the model can be found

in [19]. Here, we adapt that derivation to the configuration used in the experimental set up (described later).

The background medium is two-layered with the upper half-space being air, while the lower one models the soil. The two half-spaces are separated by a planar interface (i.e., the air/soil interface) located at *z* = 0. The scattering targets are located in the lower half-space (i.e., for *z* < 0) and buried in close proximity to the separation interface. Moreover, the "transverse" investigation domain (i.e., the spatial region where the targets can belong to) is denoted as *D* = [−*xM*, *xM*] × [−*yM*, *yM*]. The 2D investigation domain is considered at a fixed depth *zT*. Generally, we will consider *zT* = 0 (which is just at the separation interface). Reconstructions at different *zT* < 0 can be considered as well in order to explore the depth.

The scattering scene is probed by a single source located in the upper half-space at some stand-off distance from the air/soil interface, *ht*, whereas the field scattered by the buried targets is collected over a set of sensors still located on the air side and all at the same height *hO*. Accordingly, the spatial measurement positions lay over a plane; say **r***<sup>n</sup>* = (*xn*, *yn*, *hO*), *n* = 1, 2, ... , *NO*, with their positions, *NO* being the number of spatial measurements. Figure 1 shows a pictorial view of the scattering configuration along with the adopted reference frame.

**Figure 1.** Reference system: the air/soil interface is at *z* = 0, the investigation plane at *zT*.

After the scene is illuminated by the incident field, the scattered field arises. Since all the targets are located in the half-space *z* < *zT*, by invoking the equivalence theorem [21], the scattered field can be considered as being radiated by equivalent surface currents supported over the plane at *z* = *zT*. In particular, by filling the half-space *z* < *zT* with a perfect electric conductor, only the magnetic surface current survives. Such a current can be expressed as

$$\mathbf{J}\_m(\mathbf{r};k) = \mathbf{J}\_{eq}(x,y;k)\delta(z-z\_T). \tag{1}$$

Note that the magnetic equivalent current depends on the scattering scene and on the incident field. As such, it depends on the working frequency, which is indeed highlighted in (1) in terms of the wavenumber *k*. In particular, if *zT* = 0 (i.e., just at the air/soil interface), the current in (1) radiates in free space, and hence the scattered field (in the upper half-space *z* > 0) can be written as

$$\mathbf{E}\_{\rm S}(\mathbf{r};k) = \int \nabla\_{\rm r} \mathbf{g}(\mathbf{r}, \mathbf{r}';k) \times \mathbf{J}\_{\rm ll}(\mathbf{r}';k) d\mathbf{r}' = \int\_{D} \underline{\mathbf{G}}(\mathbf{r}, \mathbf{r}';k) \cdot \mathbf{J}\_{\rm eq}(\mathbf{r}';k) d\mathbf{r}'\_{l\prime} \tag{2}$$

where **r** = **r** *<sup>t</sup>* + *zTz*ˆ is the source point, with **r** *<sup>t</sup>* = (*x* , *y* ), **r** is the field point and

$$g(\mathbf{r}, \mathbf{r}'; k) = -\frac{e^{-jk|\mathbf{r} - \mathbf{r}'|}}{4\pi|\mathbf{r} - \mathbf{r}'|},\tag{3}$$

is the free-space 3D scalar Green function. *G*(**r**,**r** ; *k*) is the magnetic to electric dyadic Green function, whose expression is given as

$$
\underline{\mathbf{G}}(\mathbf{r},\mathbf{r}';k) = \begin{bmatrix}
0 & -\frac{\partial \underline{\chi}}{\partial \mathbf{z}}(\mathbf{r},\mathbf{r}';k) & \frac{\partial \underline{\chi}}{\partial \mathbf{w}}(\mathbf{r},\mathbf{r}';k) \\
\frac{\partial \underline{\chi}}{\partial \mathbf{z}}(\mathbf{r},\mathbf{r}';k) & 0 & -\frac{\partial \underline{\chi}}{\partial \mathbf{x}}(\mathbf{r},\mathbf{r}';k) \\
\end{bmatrix}.\tag{4}
$$

It must be remarked that, to be rigorous, surface integration in (2) should run over the entire plane *z* = 0. However, since the targets are very close to the air/soil interface (or in general to the reference plane *zT*), it can reasonably be assumed that the current support is very similar to the target's cross section. Hence, *D* is chosen according to the size of the spatial region to be investigated. Again, in (2), we considered the free-space Green function. This is correct for *zT* = 0. When this is not the case, because the reconstruction at a different depths is required, we will still use the same Green function. Indeed, using the free-space Green function avoids dealing with the computation of the Green function pertaining to a layered background medium. What is more, the layered medium Green function requires the knowledge of the electromagnetic features (dielectric permittivity and conductivity) of the soil, which are in general not available. Herein, such background medium electromagnetic parameters are assumed not to be known. By contrast, using the free-space Green function leads to the targets appearing more deeply located, because soil is electromagnetically denser than air. However, this is not a serious drawback if the targets of interest are shallowly buried.

The magnetic surface current **J***eq* has no component along *z*ˆ. Accordingly, (2) particularizes as

$$\mathbf{E}\_{S}(\mathbf{r};k) = \int\_{D} \begin{bmatrix} 0 & -\frac{\partial \chi}{\partial \mathbf{z}}(\mathbf{r}, \mathbf{r}'; k) \\ \frac{\partial \chi}{\partial \mathbf{z}}(\mathbf{r}, \mathbf{r}'; k) & 0 \\ -\frac{\partial \chi}{\partial \mathbf{y}}(\mathbf{r}, \mathbf{r}'; k) & \frac{\partial \chi}{\partial \mathbf{x}}(\mathbf{r}, \mathbf{r}'; k) \end{bmatrix} \cdot \mathbf{J}\_{\epsilon \eta}(\mathbf{r}'\_{t}; k) d\mathbf{r}'\_{t}. \tag{5}$$

It is seen that *ESx* is solely linked to *Jeqy* and *ESy* to *Jeqx*. Therefore, if one collects separately such field components, then the inverse problem in (5) splits in two identical scalar inverse problems from which one can reconstruct the two source components independently. Then, these reconstructions can be combined as in [20]. However, in general, this is not the case, even in view of the receiving antenna response. More precisely, what one can actually measure is the antenna output voltage. Therefore, in place of (5), the following equation should be considered

$$V(\mathbf{r};k) = \mathcal{T}(\mathbf{E}\_{\mathbb{S}})(\mathbf{r};k),\tag{6}$$

where *V* is the voltage data, and T is a linear operator schematizing the antenna response. Eventually, this is the scattering model from which to start in order to perform the current reconstructions. A few details concerning the reconstruction algorithm along with some further simplifications to achieve such a task are provided in the next section.

#### **3. Reconstruction Algorithm**

According to (6), the magnetic current components cannot in general be separately reconstructed. Moreover, the antenna response should be known and put into the model. It would be useful to avoid both these issues. Indeed, looking for simultaneously both the source components means to deal with a doubled number of unknowns. Furthermore, antenna response must be measured/known in advance.

As to the first question, if the receiving antenna is linearly polarized, for example, in the *x* − *z* plane, then its plane-wave spectrum vector belongs to the same plane. Accordingly, the main contribution to the voltage is due to *Jeqy*, which is the one that contributes to *ESx*. Similar considerations apply if the other source component is considered. Therefore, we make the problem scalar by approximating (6) as

$$V(\mathbf{r}\_{n\prime};k) = \int\_{D} H(\mathbf{r}\_{n\prime}, \mathbf{r}\prime; k) f\_{\rm{eq}y}(\mathbf{r}\prime\_{t\prime}; k) d\mathbf{r}\prime\_{t\prime} \tag{7}$$

where **r***<sup>n</sup>* are the measurement positions (antenna phase center) and *H*(**r***n*,**r** ; *k*) is the kernel of the integral operator in (7), which is actually the approximation of the composition between the antenna response operator and the propagator, once the contribution due to *Jeqx* has been neglected. However, the antenna response is still there.

It is noted that both the data and the unknown depend on the frequency *k*. Hence, employing all the available multiple-frequency data to perform the reconstruction is not formally allowed. On the one hand, single-frequency reconstructions do not allow one to estimate targets' depths. This is a minor drawback for shallowly buried targets. Moreover, depth at which reconstruction is achieved can be changed. On the other hand, at a single frequency, the antenna response basically introduces a complex weight, which is the same for each measurement position. This means that it does not affect source localization, which is instead related to the phase term that depends on **r***<sup>n</sup>* − **r** . Hence, since we are mainly interested in the detection and the localization of the targets (i.e., we do not aim at quantitative reconstructions, even in view of the other approximations), antenna response is neglected in (7) while achieving single-frequency reconstructions. Note that this would not be the case for multi-frequency reconstructions, since the complex weight in general changes with frequency.

In order to perform the reconstruction, the presented model (at a given frequency *kl*) is discretized by representing the unknown current component *Jeqy* as a truncated Fourier series, that is,

$$J\_{eqy}(x, y; k\_l) = \sum\_{n=-N}^{N} \sum\_{m=-M}^{M} I\_{mn} \exp\left[-j\pi(\frac{mx}{\varkappa\_M} + \frac{ny}{y\_M})\right].\tag{8}$$

where *Inm* are the Fourier coefficients and the exponentials represent the two-dimensional Fourier spatial harmonics over extent of the domain of investigation . Accordingly, the unknowns of the problem now become the expansion coefficients *Inm*. The choice of *N* and *M* is linked to the so-called number of degrees of freedom of the problem, reflects the ill-posedness of the inverse problem at hand, and depends on the operating frequency as well as the investigation and the observation domain extensions [19,20]. In general, however, this discretization scheme allows one reduce to a large extent the number of unknowns, as compared to a pixel based representation. The corresponding discrete model is then obtained as

$$\mathbf{V}(k\_l) = \mathbf{H} \cdot \mathbf{I},\tag{9}$$

where **<sup>V</sup>**(*kl*) <sup>∈</sup> *<sup>C</sup>NO* is the data vector at the *<sup>l</sup>*-th frequency, **<sup>H</sup>** <sup>∈</sup> *<sup>C</sup>NO*×(2*N*+1)(2*M*+1) is the matrix version of the scattering operator, and **<sup>I</sup>** <sup>∈</sup> *<sup>C</sup>*(2*N*+1)(2*M*+1) is the (vectorized) expansion coefficient matrix.

Equation (9) is inverted for **I** by a standard Truncated Singular Value Decomposition (TSVD) [25] of the relevant matrix operator. This allows to counteract the ill-posedness of the problem and to obtain a stable reconstruction.

Once the coefficients *Inm* have been recovered, the corresponding equivalent current *Jeqy*(*x*, *y*; *kl*) is computed by means of (8). Then, the support of such an equivalent surface current is provided by the image in the *x*, *y* investigation domain,

$$I(\mathbf{x}, y; k\_l) = |I\_{\text{eqy}}(\mathbf{x}, y; k\_l)|. \tag{10}$$

In order to limit the system complexity, the number of spatial data must be reduced. This in general can lead to a reconstruction that is corrupted by aliasing artifacts that are difficult to distinguish from the actual current. To cope with this drawback, a simple strategy based on the combination of single-frequency reconstructions has been introduced in [20]. In more detail, suppose that *Nk* is the number of adopted frequencies; then the final reconstruction is obtained as

$$I(\mathbf{x}, \mathbf{y}) = \Pi\_{l=1}^{N\_k} I(\mathbf{x}, \mathbf{y}; k\_l). \tag{11}$$

The very basic idea behind (11) is that aliasing artifacts change positions with the working frequency, whereas the actual source reconstruction does not. Therefore, (11) tends to mitigate all those peaks in the reconstruction that do not overlap (or overlap only partially) while the frequency changes. A criterion for the choice of the frequencies is provided in [20].

In sum, the algorithm presents the following steps:


#### **4. Experimental Results**

In this section, we check the proposed algorithm against some experimental measurements collected in a semi-controlled scattering scenario. In particular, we first describe the test site and then show some reconstructions aiming to highlight the role of the number of spatial measurements and of the employed frequencies.

#### *4.1. Test Site*

The test site consisted of a tank full of sand of about 3.5 m (length) 2.5 m (width) and 1.5 m (depth) in size. The tank was placed in the open air so that the sand appeared wet, apart from the very surface layer, which was dried by sun. The electromagnetic features of the sand were unknown and wer not estimated for detection purposes.

The transmitting antenna was a horn positioned at a *ht* = 1.5 m height from the sand floor and located at one of the end sides of the tank. It was tilted to point to the spatial region under investigation. The receiving antenna was still a horn and was located on the other side of the tank. In particular, it was mounted on a wooden slide that allowed it to synthesize a planar measurement aperture at a fixed height from the air/sand interface (in the following examples *hO* = 80 cm or *hO* = 130 cm). Furthermore, the receiving antenna was tilted toward the investigated spatial region and was linearly polarized. Figure 2 shows a schematic view of the measurement configuration along with some pictures of the test site. As a target, a metallic rectangular plate 17.5 cm × 48 cm in size is considered.

A vector network analyzer was connected to the antennas by means of coaxial cables. Standard calibration at the end of each channel was performed at the beginning of each measurement stage in order to avoid mismatch between VNA and cables. Data have been acquired in the frequency band [2–9] GHz (201 equispaced frequencies) for each different position of the receiving antenna.

We performed measurements under two different conditions: *flat* and *rough* air/sand interface. Flatness was obtained by manually using a shovel for smoothing the sand floor. Of course, the obtained sand surface, though smooth, was far from "ideally" flat. Roughness interface was instead obtained by turning over the sand. For such scenarios, we took measurements with and without targets (background measurements) for comparison purposes.

It is worth remarking that the measurement scenario actually contained many features that are not accounted for in the scattering model used to develop the detection algorithm. For example, though antennas were tilted towards the scattering scene, they still presented a direct link, which implies direct coupling between the receiving and the transmitting antennas. Furthermore, because of the finite dimensions of the tank, the two-half-space medium assumption clearly does not correspond to the actual background medium. This entails that the received signal actually consisted of different contributions besides the one expected from the targets. In addition, note that, even though the medium was a perfect two-layered medium with a flat separation interface, the air/sand interface reflection always superimposes the target signals. Nonetheless, in the following reconstructions, we did not mitigate such unwanted contributions by data pre-processing.

**Figure 2.** Schematic view of the measurement configuration (**top**) and photos of the test site (**bottom**).

#### *4.2. Detection Results for Flat Air/Soil Interface*

We start by considering the case of flat air/sand interface in the sense clarified above. For this case, we collected data over a grid of 5 × 7 positions. To this end, the receiving antenna scanned the measurement aperture with a spatial step (along both the *x* and *y* directions) of 20 cm at the height and *hO* = 80 cm.

The investigation domain is a rectangle of 160 cm along the *x* and 100 cm along the *y* direction, and it is located at *zT* = 0. We also introduce the two parameters *offsetx* and *offsety*, which indicate the displacement, along the *x* and *y* directions, respectively, of the center of the investigation domain with respect to the central point of the first measurement line (see Figure 3a for the reference system). Basically, after acquiring the data, changing the investigation domain center location (by varying the offset parameters) entails looking for the targets in different spatial regions. Accordingly, the same target will appear at different relative positions. This can be considered as a way to check reconstructions' consistency and stability. A detection can be considered successful if the target localization point moves coherently with the change of the investigation domain center position.

It must be remarked that the number of employed spatial data is already below the one required if the field has to be properly spatially sampled (see [20]). This, of course, limits the highest frequency that can be used in the reconstructions. Indeed, we tested the inversion algorithm against different bands inside the available one of [2–9] GHz. As expected, when dealing with higher frequencies, even by our approach, detection is impaired because the reconstruction results were crowded by a number of artifacts. Hence, we ended up processing data collected only within the band [2–4] GHz. In particular, in such a band, *Nk* = 11 frequencies were considered in order to be as close as possible to frequency selection criterion provided in [20].

Finally, the target was located as shown in Figure 3b and roughly buried 2 cm below the air/sand interface.

**Figure 3.** (**a**) Investigation domain and measurement points; the reference system is centred on the investigation domain. (**b**) Different positions of the target with respect to the first (numbered from the left) measurement line.

In particular, Figure 4 reports the reconstructions for a target located as in position 1 sketched in Figure 3b, for different investigation domain center offsets. As can be seen, a good detection was achieved with no significant artifacts corrupting the reconstructions. What is more, the reconstructed spot moved coherently with the investigation domain displacement.

Some comments are in order here.

Firstly, we remark that the inversion algorithm does not aim at providing the shape of the targets, but rather, it is intended for target detection. This is mainly due to the strategy we adopted to combine the different single frequency reconstructions. Indeed, the proposed multiplicative combination highlighted in (11) allows us to mitigate aliasing artifacts but at the same time tends to enhance the strongest part of the reconstructions so that targets actually appear as hot spots.

Secondly, as we mentioned above, we did not process data to counteract the direct link (from the transmitting antenna to the receiving one) or the air/sand reflection. However, reconstructions do not suffer from such spurious signals. This can be justified by observing

that we have performed the reconstruction over a given spatial region: the investigation domain. Hence, the direct link should appear located outside such a region. As to the air/sand interface reflection, it actually enters in the investigation domain. However, it is not localized as the target contribution and tends to be spread over the whole investigation domain. The multiplicative combination strategy hence also enhances the target reconstruction against such a contribution. In Figures 5 and 6, reconstructions corresponding to the target (still approximately buried at 2 cm) located as in position 2 (see Figure 3b) are reported. In particular, while Figure 5 has been obtained using the same source polarization as in the previous case, in Figure 6, the transmitting antenna has been rotated 90◦ so that the scene is illuminated by an orthogonal polarization. As can be seen, in this case as well, the target is clearly detected, regardless of the transmitting antenna features (in this case polarization). Actually, according to the formulation and related approximations presented in Section 2 and under Equation (7), what matters is that the equivalent source contributes to the polarization and that the receiving antenna is sensitive. Hence, even by changing the transmitting antenna polarization, it is expected that the method works as long as the previous statement holds true. This is basically what happened in Figure 6. In other words, unless the equivalent current has rigorously no component to which the receiving antenna is sensitive, the method is expected to work.

**Figure 4.** Normalized reconstructions for target located at position 1 shown in Figure 3b for different investigation domain offsets (reported in the figures). The red rectangles show the actual target positions.

**Figure 5.** Normalized reconstructions for target located at position 2 shown in Figure 3b for different investigation domain offsets (reported in the figures). The red rectangles show the actual target positions.

**Figure 6.** The same as in Figure 5 but incident field polarized orthogonal with respect to the previous case.

#### *4.3. Detection Results for Rough Air/Sand Surface*

We now turn to consider the case in which the air/sand interface was not smoothed. In this case, we consider data collected over a grid of 8 × 7 positions with the same spatial step as above but at a height *hO* = 130 cm. Note that the spatial data are slightly greater than the previous case but still under-sampled [20].

First, we show the reconstruction of a shallowly buried target (the target depth and type are the same as above) by employing all the available frequencies. These results are reported in the left column of Figure 7. As can be appreciated, the target is clearly detected and the related hot spot indicator changes position accordingly to the investigation domain center offset. On the same figure (right column), instead, we report the reconstructions obtained by processing background data, i.e., in absence of the target. Differently from the target case, now the reconstructions do not exhibit a clear hot spot. Moreover, the reconstruction changes as the investigation domain center offset varies. This confirms the previous discussion that processing air/soil interface reflection does not return a focused hot spot. Furthermore, the reconstruction corresponding to this contribution is in general different when the investigated spatial region changes. This is in particularly true here because roughness entails that the air/soil interface has different spatial details. Eventually, these results suggest a possible strategy to recognize actual targets against surface clutter. Indeed, comparing images obtained using different investigation domain offsets, the actual

targets are those ones for which the reconstructions "move" coherently with the change in the investigation domain center.

(**b**) *offsetx* = 2.0 m, *offsety* = 0.3 m

**Figure 7.** Normalized reconstructions of a shallowly buried targets for two different investigation domain offsets along with the actual target location denoted as red rectangles (figures un the left column). Normalized reconstructions background medium data, i.e., in absence of target, (figures on the right column).

In Figure 8, as well as in Figure 7, by keeping fixed the investigation domain with *offsetx* = 2.0 m, *offsety* = 0.3 m is considered. However, different numbers of frequencies are employed. As can be seen, when the number of frequencies is increased, the aliasing artifacts actually tend to disappear and the hot spots narrows. This is, of course, expected and perfectly consistent with the theoretical arguments discussed above.

The simple proposed strategy for reducing spatial data hence works very well. To further check this procedure, we consider an even more challenging case by reducing the spatial data employed to achieve the reconstructions. In more detail, we collected data over a 4 × 4 measurement grid, with twice the spatial step, that is, 40 cm. The height of the measurement aperture was still *hO* = 130 cm, and the frequencies were taken within the same band as above. Note that in this case the number of spatial data is even lower than the ones used in the first example reported at the beginning of this section.

**Figure 8.** The same case as in Figure 7 with *offsetx* = 2.0 m; *offsety* = 0.3 m. Comparison of reconstructions obtained by employing different number of frequencies *Nk*.

It can be seen in Figure 9 that even under this more challenging situation where the spatial data have been reduced further, the proposed method works very well in detecting and localizing the target and in counteracting aliasing artifacts, though the number of frequencies is actually low as well.

As a final example, we consider the case the target is buried more deeply. In particular, in this case, the plate target is buried approximately 10 cm below the air/sand interface. Figure 10 shows the reconstructions obtained by considering the image plane at different depths. In all the examples, *Nk* = 7 frequencies (which were shown to be enough for

mitigating aliasing artifacts) have been employed, and the same spatial grid measurement as in Figure 9 is considered. These results show that the proposed method is still effective in detecting the target. Moreover, as discussed above, changing the depth at which reconstruction is achieved allows one to get information about the target depth as the one for which the reconstruction is more focused. Of course, since the soil permittivity is not known (and hence not accounted for in the reconstruction algorithm), the estimated target depth does not coincide with the actual one. In particular, since the relative dielectric permittivity of a wet sand likely stands between 4 and 9, the "apparent" target depth can range from 20 cm to 30 cm, which is perfectly consistent with the result shown in Figure 10.

A brief comment is in order about the actual location of the reconstructed spots that, in most of the figures above, appear focused at the edges of the true shape. As a matter of fact, one can expect that the detected hot spot should roughly appear in correspondence with the target center. Indeed, this circumstance has been observed in [19,20], where synthetic data were employed. However, in those cases, we have considered a reflectionmode configuration (i.e., the TX (plane wave incidence) and the RXs were on the same side), the investigation domain was centered with respect to the measurement aperture, and the target was precisely parallel to the reconstruction plane. For the present case, because of the experimental set up, of course, the target position and its orientation are not precisely known; uncertainty is a few cm. In addition, the configuration is quite different since a transmission-mode set up is under consideration and the measurement aperture is actually side-looking the investigation plane. All these aspects could have contributed to the obtained results. However, we believe that such reconstructions are most likely due to the illumination, which is not uniform across the target and depends on the relative position between the TX and the target itself. In order to obtain a clear understanding of this effect, one should go through the study of the equivalent current behavior. This can be numerically achieved and is beyond the aim of the paper.

**Figure 10.** Normalized reconstructions of a target approximately buried 10 cm below the air/sand interface. Data have been collected over 4 × 4 measurement grid with a spatial step of 40 cm and *offsetx* = 1.7 m, *offsety* = 0.0 m. The actual target location is denoted as red rectangles. Comparison of reconstructions obtained using *Nk* = 7 at different depths.

#### **5. Conclusions**

In this paper, we experimentally validated a method for detecting and localizing shallowly buried targets from under-sampled multi-frequency data. The method relies on two main ingredients: a suitable scattering model and a simple procedure to process multi-frequency under-sampled data. The scattering model is based on the equivalence theorem, which allows one to cast the detection as the reconstruction of equivalent sources supported over a reference plane. This way, detection becomes a 2D problem. Furthermore, the incident field is embodied into the unknown equivalent sources. Therefore, coherence between the TXs and Rxs is not necessary when single-frequency data are employed. Reconstructions obtained at different frequencies are then suitably combined. This allows one to mitigate aliasing artifacts and hence to achieve the reconstructions with a very reduced number of spatial measurements.

Experimental results confirm the feasibility of the method even under a rather complex scattering scenario, such as the one addressed herein. Indeed, it is shown that the targets can be detected by using very few spatial measurement points and a number of properly selected frequencies. In addition, the results confirm the robustness of the method against the reflection occurring at the air/soil interface, even when this exhibits some degree of roughness.

The obtained results can be meant as a proof of concept. However, they encourage the application of the proposed measurement configuration and of the related inversion algorithm to more challenging scenarios, where many targets may be present inside a non-homogeneous terrain and the sensors may be deployed at a larger stand-off distance.

**Author Contributions:** Conceptualization, A.B., R.S., and G.L.; methodology, R.S.; software, A.B. and R.S.; validation, A.B.; formal analysis, A.B. and R.S.; investigation, A.B. and G.L.; resources, A.B.; data curation, A.B.; writing—original draft preparation, A.B. and R.S.; writing—review and editing, A.B.; visualization, A.B.; supervision, G.L. and R.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partially funded by VALERE: VAnviteLli pEr la RicErca research program by Universita degli Studi della Campania Luigi Vanvitelli.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** No new data were created or analyzed in this study. Data sharing is not applicable to this article.

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

#### **References**


## *Article* **Low Cost, High Performance, 16-Channel Microwave Measurement System for Tomographic Applications**

**Paul Meaney 1,\*, Alexander Hartov 1, Timothy Raynolds 1, Cynthia Davis 2, Sebastian Richter 3, Florian Schoenberger 3, Shireen Geimer <sup>1</sup> and Keith Paulsen <sup>1</sup>**


Received: 17 August 2020; Accepted: 17 September 2020; Published: 22 September 2020

**Abstract:** We have developed a multichannel software defined radio-based transceiver measurement system for use in general microwave tomographic applications. The unit is compact enough to fit conveniently underneath the current illumination tank of the Dartmouth microwave breast imaging system. The system includes 16 channels that can both transmit and receive and it operates from 500 MHz to 2.5 GHz while measuring signals down to −140 dBm. As is the case with multichannel systems, cross-channel leakage is an important specification and must be lower than the noise floors for each receiver. This design exploits the isolation inherent when the individual receivers for each channel are physically separate; however, these challenging specifications require more involved signal isolation techniques at both the system design level and the individual, shielded component level. We describe the isolation design techniques for the critical system elements and demonstrate specification compliance at both the component and system level.

**Keywords:** microwave imaging; breast; multipath; dynamic range; software defined radio; leakage

#### **1. Introduction**

Interest in microwave imaging for medical applications has grown significantly since the early 1980s, but more prominently over the past decade. Much of the recent expansion has resulted from advances in technology that facilitate implementation of novel ultrawideband approaches [1–3] and the advent of powerful computing that can accommodate time-consuming inverse problems [4–6]. Microwave imaging offers important advantages over other emerging technologies. For instance, microwave methods achieve far greater penetration depths relative to optical imaging techniques, which allow broader applications to be pursued [7–9]. Compared to lower frequency electrical impedance approaches, microwave imaging systems provide superior spatial resolution [10–13]. The list of medical microwave imaging applications is substantial and growing, and currently includes breast cancer imaging [14,15], stroke diagnosis [16,17], bone health assessment [18], and thermal therapy monitoring [19–21], among others. In all of these cases, microwave methods either exploit significant tissue dielectric property contrast [16,17,22–25], or inherent property dependencies [26,27]. While some debate over breast tumor/normal tissue property contrast remains, most studies have demonstrated that substantial contrast exists between tumor properties and those of both the associated adipose and fibroglandular tissue [22–29]. Further, several clinical publications indicate effective tumor detection, diagnosis, and treatment monitoring based on endogenous dielectric property contrast

between tumor and normal tissue [14,15,30]. Although the potential of medical microwave imaging is significant, tangible progress in terms of physical implementations has been slow.

Hardware challenges arise from the need to devise a near-field configuration that illuminates the body from multiple directions and receives scattered waves in either monostatic or multistatic modes. In multistatic cases, the difficulties are compounded by antenna mutual coupling, multipath signal distortion within the illumination array, and possible cross-channel leakage within the measurement electronics, themselves [31]. Various strategies have been developed to contend with these signal corruption problems, for example, by multiplexing high-end vector network analyzers (e.g., [32,33]), utilizing lossy coupling media (e.g., [6]), and/or translating small numbers or single source–detector pairs (e.g., [1]).

An alternative and increasingly popular hardware option [34,35] is centered around a measurement system design based on inexpensive, commercially available software defined radios (SDRs). We have described a 4-channel prototype using the B210 USRP boards developed by Ettus Research (Santa Clara, California) [36]. These units incorporate the major functional elements of a sophisticated transmitter and receiver into a single chip (Analog Devices AD9361, Norwood, MA), and operate over an impressively large bandwidth of 70 MHz to 6 GHz. They are finding use in a range of radar applications [34]. Measurement requirements for our log transformed tomographic image reconstruction algorithm along with the associated calibration process have been stable and robust over several generations of hardware systems. Our microwave tomography implementation that does not require a priori information for convergence to the desired solution, and does so with less measurement data relative to competing implementations [37]. These innovations are integral to our success in translating this technology into the clinic [14,15]. Our software advances are described elsewhere, and a comprehensive summary of the most recent innovations appears in Meaney et al. [37]. The system design includes a synchronization strategy, a technique to enhance dynamic range, and a calibration process that compensates for uneven amplitude and phase steps associated with internal variable gain amplifiers [36]. In this paper, we focus on the significant challenges associated with channel-to-channel leakage, packaging SDRs into a compact and light form factor that fits underneath the antenna mounting fixture for our existing microwave breast imaging system, and instrumentation control through existing software platforms. In the case of channel-to-channel leakage, we investigate the problem in terms of (a) leakage through existing transmission lines, and (b) ambient leakage from signals escaping their immediate structures and propagating along the overall structure as surface waves.

In this context, the novel aspects of the system include: (a) enhanced dynamic range, (b) multichannel operation with high degree of channel-to-channel isolation, and (c) data acquisition acceleration via reordering of the data recording sequence. Integration of these attributes enables the realization of an SDR-based acquisition system with data frames rates short enough for use in clinical breast examinations. We demonstrate a compact and efficient packaging of a complete multichannel system, and present results that confirm excellent dynamic range and cross-channel isolation.

#### **2. Methods**

#### *2.1. Overall System Configuration*

Figure 1 shows a schematic diagram for a 16 channel system. It includes: (a) a single B210 USRP board operating as a dedicated transmitter, (b) a single-pole, 16 throw switch (SP16T) for multiplexing the transmission channel, (c) a 6-bit, digital attenuator (0−63 dB) for controlling the reference signal amplitude, (d) a series of 16 single-pole, double-throw switches (SPDT) for toggling individual channels between transmit and receive modes, (e) a series of 16 single-pole, single-throw switches (SPST) for extra isolation between the SP16T and SPDT's (not shown), (f) a series of low noise amplifiers (LNAs) for increasing received signal strength, (g) a series of 16 monopole antennas, (h) eight B210 USRP boards for dedicated reception, (i) a SPST in combination with an 8-way power divider for synchronizing reference signals to each receiver board, and (j) an OctoClock-G (Ettus

Research, Santa Clara, CA) for providing coherent clock signals and an accurate reference for B210 signal synthesis. Control, bias lines, USB communication lines, power supply, and controlling computer are not shown. Since our imaging system approach only utilizes transmission data, incorporation of reflection measurements is unnecessary, and would increase system complexity substantially. Substantially greater channel-to-channel isolation was achieved by separating the transmit and receive functions physically and utilizing a dedicated B210 (Ettus Research, Santa Clara, CA) for transmission. The signal from Channel 1 of the transmit board drove individual antennas as each radiated into the illumination tank (illumination tank shown later in Figure 5). Part of the signal was separated via a power divider and passed through a digital attenuator after which it acted as a variable reference signal for synchronization—a process described in more detail in Meaney et al. [36].

**Figure 1.** Schematic diagram of the complete system illustrating: (**a**) transmitting B210 board, (**b**) transmitting SP16T switch, (**c**) 16 switch/amplifier modules, (**d**) eight dedicated receive B210s, (**e**) switch module for the reference signal, and (**f**) 1-by-8 power dividers for the reference signal, respectively.

Sets of SP16T, SPDT, and SPST switches direct the transmission signal to a single antenna (one at a time) while allowing the reception of response signals by the remaining 15 channels. Transmitting signals must not leak into receiver ports where they will be amplified and detected erroneously by the receive B210 modules. Here, maintaining adequate isolation is challenging because transmitted signals are on the order of 0 dBm while received signals on the order of −140 dBm need to be detected. This level of isolation is uncommon in most measurement systems and demands significant attention. In particular, signal attenuation from transmitters to receivers directly opposing the signal source often exceed −120 dBm and can reach −140 dBm or more at higher frequencies. Lossiness in the coupling liquid confers substantial benefits, specifically, in suppression of surface waves that can corrupt desired signals [31]. Achieving a large dynamic range with concomitant channel-to-channel isolation is challenging but essential for high fidelity data acquisition. As a result, measurement of signals at these low power levels is critical to system success, especially in the breast imaging context. The SP16T switch was purchased in a connectorized configuration from Universal Microwave Components Corporation (UMCC SR-J010-16S, Alexandria, VA) while the SPDT and SPST switches

(Peregrine Semiconductor PE4240 and PE4246, San Diego, CA) were integrated into a custom housing (one per channel) with the LNA.

LNAs (Mini-Circuits TSS-53LNB+, Brooklyn, NY) were added to the design to improve the receiver noise figure and dynamic range. The dynamic range of the B210 s by themselves is limited by the internal A/D discretization error at the low end instead of the actual noise level. By adding extra gain in front of these boards, the low end of the dynamic range is determined by the noise floor, which can be controlled by altering the sampling bandwidth—i.e., increasing the sampling time decreases the noise floor. This design approach was discussed in depth in Meaney et al. [36].

The B210 SDRs are driven primarily by the AD9361 (Analog Devices, Norwood, MA) as an agile transceiver, which incorporates two channels of transmission and four ports of reception, respectively (only two of the latter can be used simultaneously). As discussed in detail in Meaney et al. [36], multiport synchronization is possible and leads to a coherent, multichannel system. Variable gain was included in order to measure signals at very low received power depending on sampling bandwidth. In this configuration, Rx ports were utilized for signal reception and one of the Tx/Rx ports was used for reference signal handling. Since the on-board LO oscillators were synchronized, sampling the reference signal on only one of the two remaining ports was required. Furthermore, reference signal sampling could not be performed simultaneously with the signal measurements from the antenna on the same channel, since the Tx/Rx port shared the same pin on the AD9361 chip with the corresponding Rx port. In fact, as discussed in Section 2.5.3., the reference signal was sampled simultaneously with the antenna detection signal from the complementary channel on the same board (only their phase differences were coherent), and then a second scan samples both antenna signals (multiple subtractions produced coherent versions of both antenna—see Section 2.5.3.). In addition, more B210 shielding is critical relative to cross channel leakage, and is discussed in detail below.

The coherent reference signal was fed into a module with three SPSTs placed in series after a digital attenuator that was inserted for extra isolation, which was necessary because the reference signal is not turned OFF during antenna transmission and becomes a potential multipath propagation conduit. The signal was then fed into an 8-way power divider (Pulsar Microwave PS8-12-454/5S, Clifton, NJ) from which the Channel 1 Tx/Rx ports of the receive B210s were driven. Feeding eight receive B210s allowed reception to be performed simultaneously on all 15 channels. The system synchronization utilizing this external reference signal configuration is described in detail in Meaney et al. [36].

The OctoClock-G CDA-2990 (Ettus Research, Santa Clara, CA, USA) provided eight pairs of outputs to the B210 SDR. The first signal is a stable, 10 MHz reference that was used to synchronize the B210s, and the second is a 1 Hz PPS clock signal that was utilized for B210 coordination. In our system configuration, eight signal pairs service nine B210s: one pair is fed into two low frequency power splitters that subsequently drive two B210s.

The list below summarizes the overall system costs. Costs for instruments with comparable performance (such as the Keysight M9800A and Rohde and Schwarz ZNBT) are in excess of \$120,000. Utilizing the Ettus B210 SDRs as the foundational building block offers a substantial price advantage.


#### *2.2. Isolation and Shielding*

Two components of the hardware require significant shielding to attenuate signals propagating into the environment and back into the system—the B210 USRP circuit boards and the switch/amplifier modules.

#### 2.2.1. B210 USRP Housing Design

B210 USRP circuit boards can be purchased with no shield or one provided by the manufacturer as shown in Figure 2a,b. Our system was comprised of nine B210 boards—one serving as a dedicated transmitter and the remaining eight acting as dedicated receivers (two channels per board to produce a 16 channel system).

(**a**) (**b**)

**Figure 2.** Photographs of the B210 USRP circuit board (**a**) without and (**b**) with a commercial cover.

In order to develop a suitable shield, we explored several options before converging on the final design shown in Figure 3. An initial design, based on housing walls configured to accommodate existing B210 connectors, was unsuccessful (signals easily escaped through gaps between connector flanges and the housing wall). Instead, coaxial SMA connectors were removed and panel mount flange connectors were utilized to achieve a significantly better seal. The power connector (lower left in Figure 2) was also replaced with a panel mount coaxial cable connector to increase isolation further. The USB3.0 connector was retained on the board and was accommodated through a cutout in the housing wall around its casement.

**Figure 3.** Photograph of the B210 USRP circuit board mounted inside its custom shielded housing and associated cover, which exhibits a central ridge that isolates RF fields from the digital portion of the circuitry.

The shielded housing cover was designed with a raised surface that extended 0.5 mm into the main housing chamber, and fit snugly within the bottom of the housing—adding isolation from ambient fields escaping through seals between housing parts and aligning components to minimize damage during assembly. The board shown in Figure 3 was logically divided into two functional areas: (a) the RF componentry (right) and (b) the digital interface, power supplies, and FPGA controlling unit (left). A gold band was printed on the top surface with plated-through vias to the ground plane below. We incorporated a ridge into the housing cover that closely contacted the gold band to prevent leakage of RF signals from the right portion of the board.

#### 2.2.2. Switch/Amplifier Housing Design

The switch/amplifier module (see Figure 4) incorporated a single-pole/single-throw (SPST) switch (lower left), a single-pole/double-throw (SPDT) switch (lower right), and a low noise amplifier (LNA; top right). The SPDT selects whether a particular antenna operates in transmit or receive mode. In the transmit mode, the signal from the transmit B210 passes through both switches directly to the antenna port for propagation into the illumination tank. In the receive mode, the signal passes from the antenna port, through the SPDT and LNA to the receive port. This SPDT selects the receiver mode, and the LNA added overall gain to the composite receive channel (including the receiver portion of the B210 boards) and reduced the effective channel noise figure [36]. The SPST also added extra isolation so that signals did not leak from the transmitter to the receiver through the circuitry, but instead propagated along the desired path into the tank and through intervening object. The design goal for minimizing leakage from the transmit to receive ports was 80 dB attenuation. The design compartmentalized components such that openings between each acted as quasi-cutoff waveguides, which attenuated unwanted signal propagation between partitions—leaving the coplanar waveguide transmission line as the only viable signal path. Bias and control lines were connected via resin sealed, bolt-in filters (API Technologies Corporation—Spectrum Control—part # 51-729-312, Schwabach, Germany) threaded into the housing floor. Previous analysis of connectorized, single component devices demonstrated that bias and control lines created significant leakage pathways. The extra filtering impacted ON/OFF switching speed; however, elimination of signal corruption was deemed more important even though the design degraded data acquisition speed slightly.

**Figure 4.** Photograph of the switch/amplifier module illustrating compartmentalization of the single-pole/single-throw (SPST; left), single-pole/double-throw (SPDT; lower right), and low noise amplifier (LNA; upper right), and associated cover with raised surfaces.

#### 2.2.3. Antenna Mutual Coupling

While not part of the microwave transmit/receive electronics, antenna mutual coupling warrants some discussion. The system feeds an array of 16 monopole antennas submerged in a glycerin: water bath. Figure 5 shows (a) a photograph of a test illumination chamber and (b) a schematic diagram of the imaging field-of-view with a test object present. In our measurement system, a complement of 13 antennas directly opposite the transmitter receives signals from each transmitting antenna. We have studied the effects of the presence of adjacent antennas on the signals received. In Paulsen and Meaney [38] and Meaney et al. [39], we identified the signal perturbations caused by adjacent antennas and developed a strategy for compensating for the effects by modeling each adjacent antenna as an electromagnetic sink. The compensation is effective and produces substantially improved images. Interestingly, at frequencies above 1000 MHz, these perturbations were not observed, primarily because the lossiness of the glycerin-based coupling liquid increased substantially with frequency, and effectively suppressed the interantenna interactions.

**Figure 5.** (**a**) Photograph of a test illumination chamber and (**b**) schematic diagram of the imaging field-of-view with 16 monopole antennas and the presence of a yellow test object.

#### *2.3. Packaging*

Key features that directly addressed expansion of our 4-channel prototype [36] were: (1) physical separation of receive modules, (2) separation of the transmit module from receive modules, (3) SPST addition to the switch/amplifier modules to minimize signal leakage from the transmitter when operating as a receiver, and (4) component compartmentalization in the switch/amplifier modules to reduce internal multipath propagation of unwanted surface waves.

Figure 6a,b shows two views of the microwave subsystem as (a) fully assembled and in a (b) separated configuration. The unit was mounted below the antenna array in the operational mode as shown in Figure 6c. Components were organized to minimize wiring and cabling lengths and complexity. For example, groups of eight switch/amp modules had their associated digital I/O cards and power supply screw terminals located directly above to keep bias and control lines as short as possible. Connections to external computers and power supplies were limited to two USB cables, a single 110 V 60 Hz power cable, and two DC power supply wires.

**Figure 6.** Photographs of the microwave electronic subsystem showing: (**a**) a complete system fully assembled external to the imaging system, (**b**) electronics with the grouping of eight shielded B210s (bottom) and switch/amplifier modules (top; USB hubs removed to expose componentry), respectively, and (**c**) a complete system integrated below the imaging tank and supported by the antenna array mounting plate.

#### *2.4. Data Acquisition Sequencing and Time Considerations*

Several techniques were combined to accomplish efficient data acquisition: (1) transmission runs continuously, (2) time stabilization occurs at start-up, and (3) transmission antenna switching occurs through the SP16T whose rise time is well under 1 μsec.

The data acquisition sequence was arranged so that a single receive module was selected first and then transmission was incremented over the remaining antennas. Delays because receiver gain levels need to be adjusted did occur (depending on which transmitting antenna was broadcasting), but only involved one or two acquisitions before settling. Here, a single data set was comprised of 10,000 measurements obtained with a sampling rate of 10 MHz for an acquisition time of 8.0 msec. Since stabilization time for each module was identical, the process was initiated for a subsequent module several seconds before the data acquisition was completed from the current module. Data acquisition time for a single receive module from the 14 complementary transmitting antennas and for 5 frequencies required 15 s—total time to acquire data from all 8 receive modules was 2.2 min for a given antenna array position. For each transmitting antenna, a closely coupled channel was always available within the system design, which allowed associated switch/amplifier modules to be connected to the same receive B210 circuit board. Since complete isolation of the transmitting signal from neighboring receive channels was difficult, we omitted the adjacent data points, leaving 14 receive signals for each transmitter (7 boards with two receivers per board). Our previous experience has shown that image quality is not reduced significantly, if at all, when 13 receiver antennas are used instead of the full complement of 15 channels (i.e., the two receivers closest to the transmitter are omitted).

#### *2.5. Software and Performance Considerations*

Software challenges associated with maintaining data acquisition performance included: (1) retaining coherence between transmit channels on a single board—especially when changing operating frequency, (2) utilizing the receive function on Tx/Rx ports, (3) acquiring signals from different ports on the same board coherently, and (4) delay times associated with each receive B210 board when switching between boards.

#### 2.5.1. Transmit Channel Coherence

To maintain transmit channel coherence, only signals from one of the transmitting ports were used by incorporating a power divider to create two coherent waveforms, and one of which was attenuated digitally and served as the reference signal, which guaranteed coherence (Figure 1). In general, maintaining coherence overcomes limitations that occur when the two transmit signals are not phase locked and different boards, themselves, are not phase locked. Without coherence, the measurement phase becomes arbitrary. We exploited the B210 design feature that the receive channels within a single board are coherent (i.e., differences in their phase are constant), and performed combinations of signal acquisitions involving a transmitter reference signal to maintain coherence between the transmitter board and all receiver boards. The details of this design strategy are discussed in Meaney et al. [36].

#### 2.5.2. Tx/Rx Port Receive Function

Each Ettus B210 board was comprised of two channels—an Rx port for receiving signals and a Tx/Rx port for both transmitting and receiving signals. Signal leakage occurred within the B210 board, itself, from an unconnected Tx/Rx port to its associated Rx port such that a reference signal could be connected to the Tx/Rx port, and sampled on the Rx port. Respective signal power levels were optimized to mitigate cross-channel contamination and associated interference through a calibration process.

#### 2.5.3. Coherent Signal Acquisition Across Same-Board Ports

Receive signals from the two channels (Rx and Tx/Rx) on each board were coherent only when sampled simultaneously. Accordingly, reference signal amplitude was tailored to minimize signal corruption from channel crosstalk.

#### 2.5.4. Set-up Time Minimization

Relative to earlier prototypes [40], data acquisition was reconfigured to loop over receive modules, and then loop over associated transmitting antennas for each module. Transmitter channel selection was performed by an external SP16T switch (not by the B210 boards), which offered rise times on the order of 10 nsec, which were negligible compared to receive board start up times. Only a few receive boards can be operational at any given time without overloading the control computer. Board set up times were minimized by starting each subsequent board before its previous counterpart had completed its data acquisition cycle. Under these conditions, the majority of next-board set up times was performed while data acquisition continued, which allowed sets of data consisting of 5 frequencies and 7 vertical antenna positions to be completed in approximately 15 min per breast.

#### 2.5.5. System Calibration

Our imaging algorithm utilizes measurement data formed as differences between field values acquired when an object is present in the illumination chamber relative to when the tank is empty [40]. The subtraction was performed in terms of the logarithm of field values, and therefore, processed both log magnitudes and phases. It also cancelled phase and amplitude variations associated with the individual antennas, cables, and measurement channels, and preserved measurement coherence. Coherence for this configuration was described and discussed in detail in Meaney et al. [36]. Processing log magnitudes and phases required phase unwrapping, and unwrapping strategies are described in detail in Meaney et al. [37].

#### **3. Results**

#### *3.1. Isolation of Individual B210s*

To test isolation characteristics of representative configurations, all B210 coaxial connectors were terminated with shielded matched loads and Channel 1 of each board was programmed to transmit 0 dBm signals at frequencies from 0.9 to 1.7 GHz. Configurations under the test included (a) B210 with no shielding and (b) B210 mounted in the custom shielded housing (see Figure 3). Measurements were recorded at representative locations (Figure 7) to assess shielding benefits, and are reported in Tables 1 and 2 for transmission from TX/Rx Channel 1 and 2 ports, respectively. An Electro-Metrics (Johnstown, NY, USA) hand-held probe (EM-6992) connected to Keysight Technologies (formerly HP) 8563E spectrum analyzer (26.5 GHz) sensed emissions outside of the housing at selected locations.

**Figure 7.** Photograph of the enclosed B210 with labels indicating probe measurement sites.

**Table 1.** Measured leakage signal values at 900, 1300, and 1700 MHz for circuit board locations with: (**a**) no shield and (**b**) shielded housing, respectively, for Channel 1.


**Table 2.** Measured leakage signal values at 900, 1300, and 1700 MHz for circuit board locations with: (**a**) no shield and (**b**) shielded housing, respectively, for Channel 2.


Predictably, leakage values were reduced for shielded B210s relative to unshielded boards, and advantages were more pronounced at the higher frequencies. For channel transmissions with no shield, leakage signals were high relatively at measurement points 1–4 compared to 5–8 (see Figure 7), likely due to the RF circuitry being immediately adjacent to these connectors compared to the more distant measurement sites (5–8). For points 1–4, shielded values were uniformly superior. For positions 5–8, leakage was also superior uniformly with shielding for the highest frequency, but more uneven for lower frequencies. Interestingly, instances occurred where unshielded values were nearly identical or slightly better than their shielded counterparts. In cases where shielded and unshielded values were similar, individual measurements were promising (above −94 dBm once).

Predictably, leakage values were reduced for shielded B210s relative to unshielded boards, and advantages were more pronounced at the higher frequencies. For channel transmissions with no

shield, leakage signals were high relatively at measurement points 1–4 compared to 5–8 (see Figure 7), likely due to the RF circuitry being immediately adjacent to these connectors compared to the more distant measurement sites (5–8). For points 1–4, shielded values were uniformly superior. For positions 5–8, leakage was also superior uniformly with shielding for the highest frequency, but more uneven for lower frequencies. Interestingly, instances occurred where unshielded values were nearly identical or slightly better than their shielded counterparts. In cases where shielded and unshielded values were similar, individual measurements were promising (going above −94 dBm only once).

#### *3.2. Performance of the Switch*/*Amplifier Module*

Figure 8 shows plots of insertion loss (or gain) for cases where (a) signals were transmitted to an antenna, (b) signals were received from an antenna and amplified prior to being measured by the B210 receiver, and (c) leakage signals from a transmitter port to receiver port when the switch/amplifier module operated in the receive mode. The latter leakage signal is also shown when the housing cover was removed (i.e., no shielding). For (a), insertion loss was mild and less than 2.4 dB up to 2 GHz, which is nominally the highest frequency used in the imaging system. The result implies that the SPST and SPDT switches were operating efficiently and power loss was minimal in the system transmission mode. Similarly, gain for signals received by antennas (and directed to receiver boards) was relatively flat and ranged from 20.0 to 21.3 dB with modest monotonically decreasing gain up to 2 GHz. Isolation was greater than 80 dB over the frequency range of interest when the housing enclosure was available, and remained above 70 dB up to 3 GHz. In the unshielded case, values were substantially worse and uneven across the frequency range likely due to excitation of surface wave modes that propagated over the circuit board and housing surfaces and recombined in unpredictable patterns with the desired signals at the output. For the enclosed case, small channels between compartments acted as cutoff waveguides and provided isolation for all modes other than the desired Quasi-TEM transmission line mode that led to well-behaved transfer characteristics.

**Figure 8.** Plots of the switch/amplifier insertion loss (or gain) for the transmission mode (between the Tx and Ant ports), and receive mode (between the Ant and Rx ports). Leakage between transmit and receive ports while operational in the receive mode is shown for the completely shielded housing and the compartmentalized housing without cover, respectively.

We also tested housing isolation for ambient signals. Figure 9 shows a SolidWorks rendering of the switch/amplifier module with four associated measurement sites: (a) connected to the antenna (Ant), (b) connected to the transmission network (Tx), (c) connected to the receive B210 boards (Rx), and (d) not directly associated with a connector. In one test, a 0 dBm signal was supplied to the Tx port while the other ports were terminated with shielded 50 ohm loads. Switch control lines were set so that the signal was transmitted to the antenna port. Measurement data are presented in Table 3. Interestingly, no significant difference occurred in shielded and unshielded values for the three frequencies reported in the table. Overall, measurements were lower than −90 dBm. Signals measured near Rx ports were noticeably lower. While signals normally emanating from this port were amplified by the low noise

amplifier by about 20 dB, the SPDT switch restricted signal propagation to the port, which more than compensated for the additional gain.

**Figure 9.** SolidWorks 3D rendering of the switch amplifier housing with four measurement sites.

**Table 3.** Average transmission network (Tx) mode leakage signal measurements at designated locations for 900, 1300, and 1700 MHz, respectively.


In a second test, a −10 dBm signal was fed into the antenna port while the remaining ports were terminated with shielded matched loads. The SPDT control line was set to direct a signal coming from the antenna port to pass through the low noise amplifier and Rx port. The measurement data are presented in Table 4. Substantial improvement in leakage occurred with the shielded housing relative to the unshielded case for all measurements. Interestingly, signals measured at the Rx site were about 10 dB higher than those from the other sites, which most likely reflected the fact that these signals have been amplified previously by the LNA (which does not occur for the other ports).

**Table 4.** Average receive B210 board (Rx) mode leakage signal measurements at designated locations for 900, 1300, and 1700 MHz, respectively.


#### *3.3. System Isolation Specifications*

Antenna ports in the switch/amplifier modules were terminated with 50 ohm loads and a 0 dBm signal was the output from the transmit B210 and response signals were measured at each of 16 channels with receiver gain levels set to their maximum (76 dB). Measurements were repeated for frequencies from 700 to 1900 MHz in 200 MHz increments. Figure 10 reports the isolation level measured at each port when the transmit channel was selected to be Channel 1, and remaining ports were activated in the receive mode. Under these conditions, isolation levels were largely below −140 dB and closer to −150 dB. A few measurements for the 700 MHz case rose above the −130 dB level. A few of the 1100 MHz values were also above the −140 dB level.

**Figure 10.** Isolation levels measured at receivers for signal transmission from Channel 1 at 7 frequencies when remaining channels were activated in the receive mode and antenna ports were terminated with a 50 Ω matched load. Except for the 700 MHz case, all values are less than or equal to −135 dB.

#### *3.4. Images Reconstructed from Measurement Data*

For the imaging experiment, we utilized a 2.7 cm × 2.7 cm square cylindrical phantom filled with 60% glycerin and offset 3.2 cm from the center. Its properties were ε<sup>r</sup> = 51.53 and σ = 1.18 S/m at 1100 MHz. The background medium was an 80% glycerin bath having ε<sup>r</sup> = 29.87 and σ = 1.25 S/m at 1100 MHz. The height of the liquid was adjusted to be the same as when the object is present and when it was absent. Figure 11a shows measured magnitudes for a single transmitting antenna when the object signals were calibrated, i.e., they had the measurement values for the homogeneous bath subtracted (phases were also shown for completeness). Here, leakage was assumed to arise from any one of the situations described above. Measured data were converted from the log magnitude/phase format to the associated complex numbers. We generated random numbers between −1 and 1 for the real and imaginary leakage signal components and multiplied them by amplitude factors, which were equivalent to those associated with the prescribed leakage simulation levels utilized in the results (in these cases, phases and amplitudes of multipath contributions are unknown. Accordingly, we utilized this random number approach to control the amplitude within a worst-case scenario for the phase. Ultimately, the experiment assesses trends in the presence of progressively distorting signals). These numbers were then added to the complex-valued versions of the original data. Once these terms were summed, they were transformed back to their log magnitude and phase forms, and used directly in the reconstruction algorithm. Here, deviations between the non-leakage case and the background were more pronounced (because these values were normalized to those of the background, the background values were effectively zero). Plots confirm the observations in (a), i.e., that the field values varied most for the most remote antennas. They also show signal deviations increased as leakage levels increased for both magnitude and phase.

Figure 12a reports reconstructed relative permittivity and conductivity images at 1100 MHz with no leakage signal. Images were reconstructed with a Gauss–Newton iterative algorithm under log transformation [37]. Regularization consisted of a standard Levenberg–Marquardt scheme, and the initial image estimate corresponded to the properties of the homogeneous bath. The algorithm completed 50 iterations with an iteration step size of 0.2, which required 5.5 min of computation time to produce each image pair. The permittivity inclusion was recovered at the correct location and with the right size. Artifact levels were minimal in the permittivity images but higher in the conductivity case. In addition, the squared shape of the inclusion cross section is also evident. Conductivity images have the inclusion recovered with the correct size and location, approximately, and the recovered inclusion property values were lower than the background as expected.

**Figure 11.** Plots of measurement data for antenna 1 transmission including cases with varying levels of added leakage signal: (**a**) raw magnitude, (**b**) calibrated magnitude, and (**c**) calibrated phases. The dashed black line represents the measured signal for the homogeneous bath whereas the solid black line indicates results when the object is present, respectively. Colored lines symbolize signals when the object is present but with progressively increasing leakage added.

**Figure 12.** Reconstructed relative permittivity (**top**) and conductivity (**bottom**) at 1100 MHz for (**a**) no signal leakage, and (**b**–**e**) signal leakage of −130 dB, −120 dB, −110 dB, and −100 dB, respectively, of the square object depicted in Figure 5 for the 14.2 cm diameter field of view.

Figure 12b–e shows images from the same experiment but with data with increased levels of leakage signals. For the −130 through to −110 dB leakage, the permittivity object was recovered properly but with increasing degrees of artifacts. For the −100 dB case, the image reconstruction diverged. The conductivity object was visible in the −130 dB case, and not discernable at all once signal leakage reached −110 dB. Artifacts increased with increasing leakage signal.

Figure 13a,b shows the plots of permittivity and conductivity values along transects through the center of the imaging zone (the actual property distributions are shown for reference). These results confirm previous observations that the algorithm recovers correct trends in the actual permittivity distribution, although peak values in the object are attenuated because of smoothing inherent in the regularization process. The conductivity profile is flatter as expected because contrast between the object and background is low, and the algorithm recovers values slightly lower than the background in the proper location, even though the object properties are only 5.6% lower than the background. Artifact levels are higher in the conductivity images, which is a consistent feature of our approach as well as other systems [6]. Corresponding plots are also shown for the different leakage levels. These results confirm that when leakage is suppressed sufficiently, quality images can be recovered from the data. As leakage levels increase, image quality progressively decreases.

**Figure 13.** Horizontal transects through the 1100 MHz (**a**) permittivity and (**b**) conductivity images shown in Figure 12.

#### **4. Discussion**

We examined leakage sources from critical elements in the system—the B210 transceivers and the switch/amplifier modules—because they were the most likely to cause signal contamination. In the case of B210 transceivers, we examined their behavior in the transmission mode while we tested the switch/amplifier modules in both transmit and receive modes. Leakage tests compared instances in which boards were operated with and without the shielding. We also measured leakage effects when boards were isolated with a commercial enclosure supplied by Ettus Research, and these results were not noticeably different from the non-shielded case. Leakage results obtained with custom shields were superior and the leakage signals that did appear were most likely derived from seams between the top and bottom housing parts. While overall system design specifications called for cross-channel isolation greater than 140 dB, the Electro-Metrics probe we used measured levels near −100 dB, which were considered acceptable because these data only reflect a one-way loss through the shielded housing.

Results for the switch/amplifier module exhibited similar trends. Here, we examined both unwanted signals propagating through existing transmission lines, and those caused by surface waves that leaked between shielded housing structures. For the former, leakage signals from the Tx to Rx port were generally below −80 dB. Given this loss would be combined with about 65 dB of attenuation from the SP16T switch associated with the dedicated transmit B210 channel, the maximum unwanted signal reaching the receive B210 was expected to be on the order of −145 dBm, which is acceptable because the lowest level signals getting to the B210 were expected to be on the order of −120 dBm. Corresponding surface wave-related leakage results were similar. Interestingly, leakage did not appear to be significantly better in the shielded case in the transmit mode. Nonetheless, overall leakage levels were low and suggest these signals were not problematic. When the probe was placed above uncovered circuits at random positions, we observed transient values often 30–40 dB higher than those measured near connectors. In the receive mode, the shielded housing offered a clear advantage and isolation data reported here were informative relative to overall system performance.

Field measurement power levels varied when no breast was present relative to when breasts of different sizes and densities were placed in the illumination tank. For example, the magnitude plot in Figure 11 for the case when no object or breast was present represented a lower bound on signal amplitude for the particular liquid and operating frequency given the lossy coupling medium. When the breast was present, signal strengths were higher, and could be as much as 30 dB higher for a fatty breast and 5–10 dB higher for denser breasts. The example emphasizes the fact that central antenna signals were impacted most by leakage levels, similarly to the phantom results, and that the data acquisition needs to handle a broad range of signal levels. Our imaging algorithm was able to recover images of quality similar to those reported in the past for standard phantoms with relatively high property contrast relative to the background. These results were consistent with image reconstruction performance obtained with data acquired with previous hardware implementations of our imaging system, and confirmed the multichannel software defined radio approach described here was operating properly.

While classical imaging defines resolution in terms of the Rayleigh criteria, we were typically able to recover objects considerably smaller than the λ/2 standard. Our experience indicates shows we could detect objects as small as λ/12, although we could not characterize their size or electrical properties, accurately. Experiments with varying levels of noise indicate that resolution was inversely related to the noise level.

#### **5. Conclusions**

We developed a high performance, multichannel measurement system that could be used for microwave imaging of biological targets. Novel aspects of the system include: (a) an unusually high dynamic range that is critical in medical applications involving lossy coupling media; (b) multichannel SDR-based design with multipath signals reduced to levels having minimal impact on reconstructed image quality; (c) reorganization of data acquisition sequencing to reduce data recording time significantly; and (d) complete integration of these features that makes clinical microwave exams possible. The system operated from 500 MHz to 2.5 GHz, allowed transmission and reception from 16 channels, had a dynamic range of roughly 140 dB, and achieved excellent cross-channel signal isolation. The system design is compact and can be located below the associated antenna array. It also included design features required to ensure coherence between transmit and receive signals. By exploiting emerging technology designed into commercially available software defined radios as a fundamental system building blocks, overall system costs and complexity were reduced dramatically while still leading to a high-performance measurement system.

**Author Contributions:** The following summarizes the contributions for each author: Conceptualization, P.M., C.D. and K.P.; methodology, P.M., T.R. and S.G.; software, A.H., T.R., S.R., and F.S.; validation, T.R., C.D., S.R., F.S. and S.G.; formal analysis, P.M. and S.G.; writing—original draft preparation, P.M. and K.P.; writing—review and editing, P.M. and K.P.; supervision, P.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by NIH grant # R01-CA240760.

**Acknowledgments:** We thank Dr. Selaka Bulumulla for his technical support in designing the microwave electronics while he worked at GE Global Research in Niskayuna, NY, USA.

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

#### **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/).

## *Article* **A Prototype Microwave System for 3D Brain Stroke Imaging**

**Jorge A. Tobon Vasquez 1, Rosa Scapaticci 2, Giovanna Turvani 1, Gennaro Bellizzi 3, David O. Rodriguez-Duarte 1, Nadine Joachimowicz 4, Bernard Duchêne 5, Enrico Tedeschi 6, Mario R. Casu 1, Lorenzo Crocco <sup>2</sup> and Francesca Vipiana 1,\***


Received: 6 April 2020; Accepted: 30 April 2020; Published: 3 May 2020

**Abstract:** This work focuses on brain stroke imaging via microwave technology. In particular, the open issue of monitoring patients after stroke onset is addressed here in order to provide clinicians with a tool to control the effectiveness of administered therapies during the follow-up period. In this paper, a novel prototype is presented and characterized. The device is based on a low-complexity architecture which makes use of a minimum number of properly positioned and designed antennas placed on a helmet. It exploits a differential imaging approach and provides 3D images of the stroke. Preliminary experiments involving a 3D phantom filled with brain tissue-mimicking liquid confirm the potential of the technology in imaging a spherical target mimicking a stroke of a radius equal to 1.25 cm.

**Keywords:** microwave imaging; brain stroke; monitoring; antenna array

#### **1. Introduction**

Brain stroke is one of the main causes of permanent injury and death worldwide, with an incidence of over 5 million annual deaths [1]. Since prompt intervention (such as the administration of specific drugs that can affect the acute stage of the stroke) can significantly improve the prognosis, a rapid diagnosis of the disease and continuous monitoring after its onset represent important clinical goals.

Currently, the most effective tool for brain stroke diagnosis is the imaging-based diagnostics performed in an emergency department after the recognition of stroke-like symptoms. In this respect, magnetic resonance imaging (MRI) or X-ray based computerized tomography (CT) are the clinically adopted techniques. However, although they are continuously evolving, these technologies still have several limitations. In particular, despite its high resolution and accuracy, MRI is not widely available in emergency settings and is therefore actually used only as a secondary diagnostic tool [2,3]. On the other hand, non-contrast CT may be limited by the fact that the early signs of ischemia may not be easily recognizable by non-experienced personnel. Moreover, due to the use of ionizing radiation, CT is

not suitable for repeated examinations, which are especially useful for post-acute monitoring purposes. Furthermore, both MRI and CT equipment is bulky and so not currently suited for ambulance use or as bedside devices.

The above circumstances have led to increased interest in the development of different diagnostic imaging techniques [3]. Among others, microwave imaging (MWI) [4] has emerged as a complementary technique which is able to address the different needs arising in stroke diagnosis and management, namely the early—possibly prehospital—diagnosis of the kind of stroke (ischemia or hemorrhage), bedside brain imaging and continuous brain monitoring for stroke in the post-acute stage. MWI takes advantage of the different electric properties (electric permittivity and conductivity) that human tissues exhibit at microwave frequencies depending on their kind (e.g., blood versus gray or white matter) and pathological status. These differences permit a functional map of the inspected anatomical region to be obtained. The benefits of MWI mainly stem from the non-ionizing nature of microwave radiation and the reduced intensity required to obtain reliable imaging (at an intensity comparable to that currently used for mobile phones), which make it completely safe and suitable for repeated applications. Moreover, MWI technology is cost-effective and benefits from a reduced size, as it makes use of miniaturized, low-cost, off-the-shelf components that are available in the microwave frequency range for signal generation and acquisition [5] and low-cost accelerators to speed up processing [6].

Recently, several MWI devices and prototypes have been proposed [7–17]. Among them, the two most prominent examples (which are already being tested on humans) are the Strokefinder, developed by Medfield Diagnostics [7,8], and the EMTensor BrainScanner [10]. The Strokefinder is a device which aims to discriminate between ischemic and hemorrhagic strokes in the early stage of patient rescue, based on an automated classification which is carried out by comparing the measured data to a database (obtained by data collected from already examined patients). This device is characterized by its very simple and compact hardware, consisting of a small number of printed antennas mounted on a support that can be adapted to the patient's head. Some initial clinical trials have been reported for the Strokefinder [7], but it should be remarked that it does not provide images; thus, its intended role is to complement standard imaging tools. The EMTensor BrainScanner aims to perform brain stroke tomography. The system is characterized by a high complexity, as it consists of a large number of radiating elements (177 truncated waveguides, loaded with ceramic material of appropriate permittivity [10]), which considerably increase its cost and size, thus reducing to some extent the advantages of its use. In addition, the image reconstruction task involves the processing of a considerable amount of measured data and has to face the pitfalls of dealing with a non-linear and ill-posed inverse problem. This entails long elaboration times and possibly results in false solutions; i.e., producing images which fulfill the underlying optimization but are different from the ground truth.

In this paper, we describe the realization, characterization and initial experimental validation of a prototype device representing a different approach to dealing with a still open issue in stroke management; that is, continuous monitoring during the hospitalization of the patient in order to evaluate the effectiveness of the administered therapies [18]. This specific application aims to image only a "small" variation occurring in the brain and not its overall structures and features. As a consequence, it is possible to keep the device complexity low, and therefore also its size and cost, as well as to rely on the Born approximation to model the scattering phenomenon, thus enabling reliable real-time imaging. Accordingly, the proposed device is based on the low-complexity architecture designed with the rigorous procedure as described in [18]. Moreover, it adopts a differential imaging approach, where data gathered at two different acquisition times are processed [19] with simple and fast imaging algorithms based on the distorted Born approximation [18,19].

The proposed system provides 3D images of the head by relying on data measured through an array of 24 printed monopole antennas organized as an anatomically conformal shape mimicking a wearable and adaptable helmet. Each antenna is enclosed in a box of graphite-silicon material, acting as the coupling medium, and connected to a two-port vector network analyzer (VNA) through a 24 × 24 switching matrix, which allows the whole differential scattering matrix required for imaging

to be acquired. The use of a semi-solid matching medium is a distinct feature of the system which allows for an increased simplicity of operation and repeatability, as compared to other arrangements that make use of a coupling liquid [10]. Finally, as detailed below, the device presented here is equipped with a "digital twin" based on a proprietary electromagnetic (EM) solver that allows us to properly characterize and foresee its behavior, as well as to provide the building blocks needed for the imaging.

In the following sections, the different components of the device are described and discussed, and a first experimental assessment on an anthropomorphic head phantom is presented. This phantom consists of a plastic shell with the shape and size of a human head, which is filled up with a homogeneous material whose dielectric properties are equal to an average value of the properties of the different tissues present in the brain. The reported experimental results provide an initial demonstration of the capabilities of the developed device.

It is worth noting that the presented system is not the only example of a low-complexity device for brain imaging, as other devices, using a low number of antennas (8–16) arranged in a circular array, have been proposed [12–16]. However, these devices only provide 2D maps of the transverse cross-section of the head in the array plane , whereas the device herein presented provides a full 3D image of the head.

#### **2. Material and Methods**

#### *2.1. Three-Dimensional Microwave Imaging System Design*

In this section, we present the choices made in the design of the proposed imaging system (i.e., the operating frequency, coupling medium, antenna number and arrangement).

The operating frequency as well as the properties of the selected coupling medium were set according to previous findings obtained from theoretical formulations and experimentally validated [20] and [4] (Chapter 2). In particular, a working frequency of around 1 GHz and a coupling medium with a relative permittivity of around 20 were determined to be optimal and therefore chosen for the realization of the hardware. As all the simulations and measurements were made at 1 GHz, this operating frequency will be considered as implicit in the following.

To design the layout of the array of antennas (the number and position of the radiating/measuring elements), a recently proposed rigorous procedure was adopted [18,21]. This procedure is based on the analysis of the spectral properties of the discretized scattering operator [22], assuming the antennas to be located on a surface conformal to the head (a helmet) and taking into account the dynamic range and signal-to-noise ratio (SNR) of the adopted measurement device as well as the actual size of the antennas. The result of this study allowed us to identify a 24-element array as the suitable candidate to perform the desired imaging task while keeping the system complexity as low as possible. The expected performances were confirmed by a preliminary numerical analysis [18].

As far as the choice of the imaging algorithm was concerned, the key aspect was that the targeted application was to monitor the time evolution of the stroke. Hence, a differential approach was a suitable choice [18,23]. In particular, the input data of the imaging problem, denoted as Δ*S* in the following, were represented by the difference between the scattering matrices measured at two different times, while the output was a 3D image showing the possible variation of the electric contrast of the brain tissues—i.e., Δ*χ*—occurring between these two different times. The differential electric contrast function Δ*χ* is defined as Δ/*b*, where *<sup>b</sup>* is the complex permittivity of the non-homogeneous background at the reference instant (e.g., a map of the brain at the first diagnosis) and Δ is the complex permittivity variation between the two time instants.

Since the contrast variation Δ*χ* was localized in a small portion of the imaging domain, it was possible to take advantage of the distorted Born approximation [22], so that a linear relationship held between Δ*S* and Δ*χ*:

$$
\Delta \mathcal{S} \left( r\_p, r\_q \right) = \mathcal{S} \left( \Delta \chi \right), \tag{1}
$$

where <sup>S</sup> is a linear and compact integral operator, whose kernel is <sup>−</sup>*jωb*/4 *<sup>E</sup>*<sup>b</sup> *rm*,*r<sup>p</sup>* · *E*b *rm*,*r<sup>q</sup>* , with *<sup>r</sup><sup>m</sup>* <sup>∈</sup> *<sup>D</sup>*, *<sup>r</sup><sup>p</sup>* and *<sup>r</sup><sup>q</sup>* denote the positions of the transmitting and receiving antennas, and *<sup>r</sup><sup>m</sup>* shows the positions of the points in which the imaging domain *D* is discretized. *E*<sup>b</sup> is the background field in the unperturbed scenario; that is, the field radiated inside the imaging domain by each element of the array. The symbol "·" denotes the dot product between vectors, *ω* = 2*π f* is the angular frequency, and *j* the imaginary unit.

As a reliable and well-established method to invert (1), we exploit the truncated singular value decomposition (TSVD) scheme [22], which allows us to obtain the unknown differential contrast function through the explicit inversion formula:

$$
\Delta \chi = \sum\_{n=1}^{L\_l} \frac{1}{\sigma\_n} \left< \Delta S\_\prime \left[ u\_n \right] \right> \left[ v\_n \right] \,, \tag{2}
$$

where *σn*, [*un*] and [*vn*] are the singular values and the right and left singular vectors of the discretized scattering operator S, respectively. *Lt* is the truncation index of the SVD, which acts as a regularization parameter and was chosen to obtain a good compromise between the stability and accuracy of the reconstruction [22].

#### *2.2. Three-Dimensional Microwave Imaging System Realization*

The realized 3D microwave imaging system prototype is shown in Figure 1. It consists of several parts, described in the following sub-sections, all of which are controlled by a laptop.

**Figure 1.** Overview of the realized 3D microwave imaging system prototype. From left to right: vector network analyzer (VNA), switching matrix and head phantom wearing a helmet interconnected to the switching matrix by means of coaxial cables.

#### 2.2.1. Vector Network Analyzer and Switching Matrix

First, all the signals are generated and received by a standard VNA (Keysight N5227A, 10 MHz–67 GHz) where the input power is set to 0 dBm and the intermediate frequency (IF) filter to 10 Hz. The two ports of the VNA are connected, via flexible coaxial cables, to the two input ports of the 2 × 24 switching matrix. The switching matrix has been realized with two single-pole-four-throw (SP4T), eight single-pole-six-throw (SP6T), and 24 single-pole-double-throw (SPDT) electro-mechanical coaxial switches. The internal connections between the switches are made with semi-rigid coaxial

cables to maximize the isolation and minimize the insertion losses. Then, the 24 output ports of the switching matrix are connected, via flexible coaxial cables, to the 24 antennas placed on the helmet as supports hosting the 3D anthropomorphic head phantom. As detailed in [23], the switching matrix has been realized so that there are 24 paths from VNA port 1 to the corresponding 24 antennas, as well as 24 paths back to VNA port 2; all the paths were designed to have the same electrical length. In this way, each antenna can work as a transmitter (TX) or as a receiver (RX), and while one antenna is transmitting, the signals collected by the other 23 antennas are received in sequence by the VNA.

#### 2.2.2. Brick Antenna Array

The antenna array is the core of the imaging system. The antenna numbers, positions and orientations were determined according to the design procedure described in Section 2.1. In the prototype shown in Figure 2, the 24 antennas are placed on a 3D printed plastic (acrylonitrile butadiene styrene, ABS) support with the shape of a helmet conformal to the head phantom. This prototype configuration allows us to easily change or remove the antennas, if needed.

**Figure 2.** Brick antenna array conformal to the head phantom.

Each antenna, denoted as "brick" antenna in the following, was a monopole antenna printed on a standard FR4 substrate, with a thickness equal to 1.55 mm, and embedded in a dielectric brick. The dielectric brick was made of a mixture of urethane rubber and graphite powder and was designed in order to reach a relative dielectric permittivity of *<sup>r</sup>* ∼= 20 and to minimize the losses. The actual EM properties obtained with this mixture were *<sup>r</sup>* = 18.3 and *σ* = 0.19 S/m. The graphite powder was needed to increase the relative dielectric constant of the urethane rubber; however, it also increased the conductivity, although still within acceptable levels compared to other matching media used in medical microwave imaging (e.g., [23,24]). Moreover, the adopted matching medium is usually liquid [10], which is inherently inconvenient for a helmet-like device. Here, instead, the implemented matching medium was solid rubber, which can be easily placed conformally to the head, as shown in Figure 2. The overall dimensions of each brick were 5 <sup>×</sup> <sup>5</sup> <sup>×</sup> 7 cm<sup>3</sup> to accommodate the need to place 24 brick antennas around the head.

Figure 3 reports the 24 × 24 scattering matrix measured by the microwave imaging system in the presence of the head phantom; the self-terms were set to zero in order to highlight the range variation of the measured coupling coefficients, which are the input data of the used TSVD imaging algorithm.

**Figure 3.** The 24 × 24 scattering matrix measured in the presence of the head phantom; the self-terms are forced to zero.

It can be noticed that the measured signals are above −100 dB; considering that the VNA noise floor is −110 dBm (at 1 GHz and with an intermediate frequency (IF) filter equal to 10 Hz), with an input power of 0 dBm the measured data are then above the VNA noise floor [25]. Moreover, as expected, the 24 × 24 matrix is approximately symmetric, which confirms the reciprocity of the realized system.

#### 2.2.3. Anthropomorphic Head Phantom

The 3D anthropomorphic head phantom used for the validation and testing of the microwave imaging system was made of polyester casting resin. It was realized by additive manufacturing from a stereo-lithography (STL) file derived from MRI scans. The STL file was obtained with computer-aided design (CAD) software by modifying an original file from the Athinoula A. Martinos Center for Biomedical Imaging at Massachusetts General Hospital [26].

The phantom consisted of a cavity, shown in Figure 4a, whose wall thickness and height were equal to 3 mm and 26 cm, respectively. Its maximum cross section was approximately an ellipse whose minor and major axes were equal to 20 cm and 26 cm, respectively. The cavity was filled with a liquid mixture, made of Triton X-100, water and salt, which mimicked the average value of the dielectric characteristics of different brain tissues (white matter, gray matter) [27]. The measured dielectric characteristics of the mixture are reported in Figure 4b.

**Figure 4.** Anthropomorphic head phantom; (**a**) structure, and (**b**) dielectric characteristics of the liquid mixture, mimicking brain tissues as averages, which fills the phantom.

#### 2.2.4. Digital Twin of the Device

The developed device was equipped with a digital twin, namely an accurate full-wave numerical model simulating its behavior, thanks to the use of proper CAD and EM simulation softwares. This tool had a two-fold purpose. First, it allowed us to foresee the expected outcomes of the planned experiments and to analyze them before running the actual experiments. Second, it provided the EM fields *E*<sup>b</sup> needed to build the imaging kernel, as described in Section 2.1.

Figure 5 shows the CAD model of the device which includes both the antenna array and the 3D anthropomorphic head phantom.

**Figure 5.** Computer-assisted design (CAD) model of the 3D microwave imaging system together with the anthropomorphic head phantom; two printed monopole antennas within the coupling medium bricks are highlighted in yellow.

Each brick, placed on the head, represented one TX/RX antenna together with the dielectric coupling medium; the antenna ports were placed at the end of the coaxial cables leading out from the bricks. The dielectric characteristics of the bricks were the same as the nominal ones used in the realized system (see Section 2.2.2), i.e., *<sup>r</sup>* = 18.5 and *σ* = 0.2 S/m. The head phantom was made of a dielectric medium representing the average brain tissues with *<sup>r</sup>* = 42.5 and *σ* = 0.75 S/m, according to the properties of the medium adopted in the experimental validation (Section 2.2.3).

To perform the EM simulation, which provided both the S-parameters at the antenna ports and the EM fields in the whole scenario, the CAD model was introduced into an in-house full-wave software, based on the finite element method (FEM) [28]. The implemented software used the standard "curl–curl" formulation for the electric field and Galerkin testing. The whole volume, including the CAD model, was discretized with edge-basis functions, defined over a mesh of tetrahedral cells. The antenna metal parts were modeled as perfect electric conductors (PEC), and all the dielectrics were modeled via sub-volumes with given *<sup>r</sup>* and *σ* values. The tetrahedral mesh was terminated at the volume boundaries with appropriate absorbing boundary conditions (ABC). Each antenna port was modeled as a coaxial cable section where, if excited, a tangential electric field distribution was enforced.

In the numerical modeling of the MWI system, the 24 antennas were excited one at a time, setting all the others as receivers in order to generate the 24 × 24 scattering matrix. Moreover, the field radiated by each antenna was evaluated at different locations within the head to generate the discretized scattering operator S in (1).

#### 2.2.5. Experiment Set-Up

The microwave imaging system was validated using the 3D anthropomorphic head phantom in which a 1.25 cm-radius plastic sphere was inserted, as shown in Figure 6.

**Figure 6.** (**a**) Experimental set-up; (**b**) 1.25 cm-radius plastic sphere.

The sphere was fixed via a monofilament polymer line (fishing line) to an external support located above the head and was immersed in the liquid mixture at a location and height that could be varied. The sphere density was higher than the liquid mixture filling the head phantom, and the sphere was therefore not floating on the liquid and could be easily positioned in different locations within the head phantom. Three positions were considered, as shown in Figure 7. However, as their exact location could be affected by inaccuracies, it is obvious that the actual positions may have slightly differed from the expected ones.

**Figure 7.** Three sphere locations (color dots) considered within the 3D head phantom: (**a**) 3D view; (**b**) 2D sagittal view.

The aim of these experiments was to identify the 3D shape and location of the sphere that was supposed to represent the region of the brain affected by the stroke. In this respect, it is important to remark that, while the plastic sphere adopted here for the sake of simplicity was not intended to mimic a hemorrhage or a clot, it showed a dielectric contrast with brain tissues which was comparable to a hemorrhage but with the opposite sign. Thus, the experiment was expected to provide a meaningful, though initial, validation of the imaging capabilities of the prototype device.

#### **3. Results**

In the following, we describe the validation of the realized 3D microwave imaging system (Section 2.2), whose experimental set-up is detailed above in Section 2.2.5. First, by means of the digital twin, the outcomes of the experiments are foreseen and assessed. Then, the results of the actual experiments are reported. All the imaging results were obtained by using the TSVD algorithm described in Section 2.1. Figure 8 shows the singular values calculated for the relevant scattering operator, which were computed with the digital twin. The truncation index *Lt* in (2) was set to −20 dB for all the reported cases. To quantitatively assess the quality of the reconstruction images, the root mean square error (RMSE) was evaluated as

$$\text{RMSE} = \sqrt{\frac{\sum\_{n=1}^{N\_t} (\Delta \hat{\chi} - \Delta \chi)^2}{N\_s}},\tag{3}$$

where *Ns* is the number of samples of the discretized domain, <sup>Δ</sup>*χ* the retrieved differential contrast and Δ*χ* is the actual contrast.

**Figure 8.** Singular value behavior of the scattering operator.

#### *3.1. Numerical Assessment*

To confirm the validity of the above statement, the digital twin of the system described in Section 2.2.4 was exploited. To this end, the experiment for the target corresponding to the blue sphere in Figure 7 was simulated by including a 1.25 cm-radius plastic sphere with *<sup>r</sup>* = 2.1 in the CAD head phantom. The simulation was repeated for the case of noiseless measurements (to provide an ideal benchmark) and for two SNR levels, namely 65 dB and 55 dB. Noise was modeled as an additive Gaussian white noise, which was added to the simulated data given by the S matrices computed with and without the target, respectively. Figure 9 shows the outcomes of the simulations. In particular, the first row shows the amplitude of the scattering matrices, where, in agreement with the actual experimental situation, the self-contributions have been omitted. It appears that the considered noise level severely affects the data matrix. The middle and bottom rows of Figure 9 show the normalized amplitude of the reconstructed differential contrast Δ*χ* arising from the TSVD algorithm. It is worthwhile to note that the target is properly imaged, even when the SNR is comparable to or higher than the maximum values of the differential S matrix and the corresponding matrices appear very noisy (see Figure 9b,c), as the unavoidable occurrence of some artifacts does not restrict the interpretation of the result. The RMSE values obtained for the reconstructions at the three considered noise levels are 0.04, 0.06 and 0.11, respectively.

**Figure 9.** Digital twin: case of a target plastic sphere; the exact sphere location and shape are indicated by red circles; (**a**–**c**) differential scattering matrices, and (**d**–**f**) cross-sections of the reconstructed images at the sphere center and (**g**–**i**) at various levels, for different values of the signal-to-noise ratio (SNR) (left column: no noise, middle column: SNR = 65 dB , right column: SNR = 55 dB).

Once the kinds of results that the device was expected to provide in the experiments had been characterized, the second goal was to show whether—and to what extent—the experiment was relevant for the considered clinical scenario. To this end, the same simulation as before was repeated considering a spherical target with the properties of blood (*<sup>r</sup>* = 63.4 and *σ* = 1.6 S/m) instead of plastic, thus mimicking a hemorrhagic stroke. The results of the simulations are shown in Figure 10. Figure 10d–i show the imaging results, which are comparable with the previous ones but for an expected increased presence of (not meaningful) artifacts. In agreement with this, the RMSE values are essentially the same as for the case of the plastic target, at 0.04, 0.06 and 0.13 for the different considered levels of noise, respectively. According to the above results, we can conclude that the developed prototype will be able to pass the planned experimental validation and that the considered experiments provide a relevant, though initial, test-bed for the designed microwave imaging system.

**Figure 10.** Digital twin: case of a target blood sphere; the exact sphere location and shape are indicated by red circles. (**a**–**c**) Differential scattering matrices; (**d**–**f**) cross-sections of the reconstructed images at the sphere center and (**g**–**i**) at various levels, for different values of the signal-to-noise ratio (left column: no noise, middle column: SNR = 75 dB , right column: SNR = 65 dB).

#### *3.2. Experimental Validation*

To perform the experimental validation of the system, for each sphere location, two measurement sets were taken at different times. The first data-set was measured in the absence of the sphere, and the second one was taken when the sphere was positioned inside the phantom. Each pair of measured 24 × 24 scattering matrices were then differentiated and given as an input to the TSVD algorithm. It can be observed that VNA port calibration was not needed as possible systematic measurement errors were cancelled out via the use of differential data or did not affect the obtained qualitative imaging of the differential contrast function. The time needed to perform the total measurement with the current prototype was less than 4 min for each sphere location, and the elapsed time between the two measurement sets was around 10 min. The kernel simulated with the digital twin (Section 2.2.4) was exploited to generate the images, and the data processing time needed by the TSVD algorithm was negligible (less than 1 s).

The 24 × 24 differential scattering matrices obtained for the three cases are shown in Figure 11a–c. The second and third rows of Figure 11 display the images obtained for the different sphere locations

indicated in Figure 7. The second row corresponds to the horizontal cross-section passing through the sphere center; the expected sphere location and size are highlighted with a red circle. The third row displays horizontal cross-sections at the different levels indicated in Figure 7b. In all cases, the targets are very well retrieved, and the images agree with the simulation results, which confirms our expectations. The RMSE values for these reconstructions are 0.16, 0.19 and 0.18 for the three considered positions of the target, respectively. Moreover, the last row of Figure 11 shows the 3D rendering of the imaged stroke obtained by plotting the values of the normalized differential contrast amplitude, which are above −3 dB. Finally, Figure 12 reports horizontal cross-sections at various levels of the 3D reconstructed image obtained by differentiating two different sets of measurements of the same scenario; all the images are normalized with respect to the maximum value of Figure 11b, which represents a case when the target is present.

**Figure 11.** Measurement results: (**a**–**c**) differential scattering matrices; (**d**–**f**) horizontal cross-sections at sphere center and (**g**–**i**) at the various levels as highlighted in Figure 7 (the expected sphere location and shape are indicated by red circles); from left to right, images corresponding to the blue, red and yellow spheres of Figure 7, respectively; (**j**–**l**) the 3D rendering of the imaged stroke.

**Figure 12.** Horizontal cross-sections at various levels of the obtained 3D image, differentiating two sets of different measurements of the same scenario; images normalized with respect to the maximum value of Figure 11b.

#### **4. Discussion**

The main goal of the developed device was to image (qualitatively) possible anomalies (clots or hemorrhages) in the head to support clinicians in the evaluation of the effectiveness of the administrated therapies. The results shown in the previous section, although preliminary, confirm the potential of the technology in providing reliable results, as it is capable of imaging a target as small as 1.25 cm in radius.

A second important achievement of the analysis carried out here is represented by the validation of the proposed system through its digital twin, which provides simulated data which are quantitatively consistent with the measured data. As such, the adopted modeling tool provides a reliable representation of the device, making it possible both to build the imaging kernel and to synthetically reproduce laboratory experiments. As a matter of fact, the results from the simulations and measurements appear to be essentially the same, except for a slight deterioration of the images in the case of measured data. In particular, the worse RMSE results in the experimental tests, with respect to the corresponding numerical cases, are due to the possible inaccuracies in the expected positions of the target as well as to model inaccuracies in the digital twin.

A very important feature of the presented device is its robustness against false positives, assessed through a specific experiment, the result of which is shown in Figure 12. We can observe that the reconstructed values, which represent only the overall noise between different sets of measurements, are significantly lower than 0 dB (the maximum value in the 3D reconstructed image is equal to −4.75 dB) and are therefore clearly different from the cases when the target is present (see Figure 11).

To the best of our knowledge, this paper presents the first system based on a low-complexity antenna arrangement conformal to the head which is able to provide full 3D images; other imaging systems available in the literature only provide 2D images and often exploit a large number of antennas.

With this work representing a preliminary validation of the developed hardware, it is worthy of note that, for practical reasons, the target used in the experimental test does not exhibit the same dielectric properties as a stroke. On the other hand, the digital twin can help us to predict what will happen in an experiment dealing with a target mimicking a hemorrhage, for example. As a matter of fact, by comparing results from simulations and measurements, it can be observed that the differential scattering matrices exhibit a similar pattern but are lower in amplitude (by about 10 dB) in the case of simulations due to the lower maximum amplitude of the differential contrast between blood and the average brain with respect to plastic and the average brain (0.49 and −0.95, respectively). This implies that, in an experiment dealing with a target mimicking a hemorrhage, slightly weaker useful signals should be collected. However, this is not a significant limitation, as the amplitudes of the differential scattering matrices are well above the VNA noise floor, which represents the ultimate limitation for accurate measurements [25].

Finally, while the repeatability of the experiment has not been tested in this paper with respect to the possible misalignment of the phantom in the two gathered data sets, from our previous studies, we expect that such uncertainties will produce "structured" artifacts in the final image, which are easily attributable to positioning errors [19].

#### **5. Conclusions and Future Work**

In this paper, a prototype of a novel low-complexity device, dedicated to brain stroke monitoring in the post-acute stage, has been presented. The reported experiments aimed to provide an initial validation of the device and confirm that there is an agreement between the design of the system [18] and its performance. This is assessed by the comparison of the achieved results with those obtained from its digital twin. The availability of such a model also allows us to show the kinds of outcomes that are expected to be obtained by the system when operated with more realistic targets.

The next steps of the system validation and development will involve an assessment with the more realistic anthropomorphic head phantom described in [27], which includes an additional cavity modeling a stroke [29]. As this phantom is derived from an STL file, a numerical validation using the digital model of the system will also be carried out in this case. In addition, improvements of the prototype will be performed by considering image reconstruction procedures using multi-frequency data calibration techniques as well as sparsity-promoting regularization schemes; furthermore, there will be a refinement of the antenna array support to make it wearable.

**Author Contributions:** All authors contributed substantially to the paper. Conceptualization, J.A.T.V., R.S., M.R.C., L.C. and F.V.; Software, J.A.T.V., R.S., G.T. and D.O.R.-D.; Validation, J.A.T.V., R.S., G.B., D.O.R.-D., N.J., B.D. and E.T.; Writing—Original Draft Preparation, J.A.T.V., R.S., G.B., L.C., and F.V.; Writing—Review and Editing, M.R.C., L.C. and F.V.; Supervision, M.R.C., L.C. and F.V.; Project Administration, F.V.; and Funding Acquisition, L.C. and F.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by the Italian Ministry of University and Research under the PRIN project "MiBraScan—Microwave Brain Scanner for Cerebrovascular Diseases Monitoring", and by the European Union's Horizon 2020 research and innovation program under the EMERALD project, Marie Skłodowska-Curie grant agreement No. 764479.

**Acknowledgments:** The authors would like to acknowledge the valuable help of Eng. Gianluca Dassano for the realization of the switching matrix.

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

#### **Abbreviations**

The following abbreviations are used in this manuscript:

MRI Magnetic resonance imaging CT Computerized tomography MWI Microwave imaging VNA Vector network analyzer EM Electromagnetic SNR Signal-to-noise ratio TSVD Truncated singular value decomposition SP4T Single-pole-four-throw SPDT Single-pole-double-throw IF Intermediate frequency TX Transmitter RX Receiver ABS Acrylonitrile butadiene styrene STL Stereo-lithography FEM Finite element method PEC Perfect electric conductor ABC Absorbing boundary conditions CAD Proper computer-aided design RMSE Root mean square error

#### **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/).

*Article*
