# A Three-Dimensional Microwave Sparse Imaging Approach Using Higher-Order Basis Functions

Authors: Nebojsa Vojnovic; Lorenzo Crocco; Marija Nikolic Stevanovic
Journal: International Journal of Antennas and Propagation
 (2022)

---

## Abstract

This paper presents a novel method for three-dimensional microwave imaging based on sparse processing. To enforce the sparsity of the unknown function, we take advantage of the fact that arbitrary three-dimensional electromagnetic fields can be decomposed into two components with respect to the radial direction: one with transverse-magnetic polarization and the other with transverse-electric polarization. Each component can be further expressed as a sum of spherical harmonics, which provide the dictionary exploited by the sparse processing algorithm. Our measurement model relates the data and the parameters of the spherical harmonics’ sources, which are uniformly distributed on a grid sampling the imaging domain. By relying on the theory of degrees of freedom of electromagnetic fields, it can be shown that only a few harmonics are sufficient to accurately represent the measured scattered field from objects whose diameter is of the order of the wavelength, thus allowing reducing the dimension of the adopted dictionary. We analyze several imaging scenarios to assess the algorithm’s performance, including different object shapes, sensor orientations, and signal-to-noise ratios. Moreover, we compare the obtained results with other state-of-the-art linear imaging techniques. Notably, thanks to the adopted dictionary, the proposed algorithm can yield accurate images of both convex and concave objects.


---

## Body

## 1. Introduction

Microwave imaging (MWI) has many practical applications, such as ground-penetrating radar (GPR) [1, 2], through-the-wall-imaging (TWI) [3, 4], medical diagnostics [5, 6], and nondestructive evaluation (NDE) [7, 8]. Among many MWI applications, microwave medical imaging is one of the fastest progressing areas since it takes advantage of nonionizing radiation and affordable components. Examples are MWI systems for brain stroke monitoring [9], breast-cancer detection [10], and thermal ablation guidance [11]. Nevertheless, one of the drawbacks of MWI systems, which restricts their widespread adoption, is the limited spatial resolution [12] compared to other medical imaging modalities.Depending on the retrieved information, MWI techniques fall into two categories. The goal of the first group of algorithms, the so-called qualitative methods, is to find the shape and position of the inspected object. Examples of qualitative imaging methods are time reversal [13], MUSIC [14], linear sampling method (LSM) [15, 16], level set [17], truncated singular value decomposition (TSVD) [18, 19], sparse processing [20, 21], and many others. If additional information is needed, such as the permittivity and conductivity of the target, quantitative MWI methods are implemented [22, 23].In general, solving MWI tasks is difficult due to the inherent nonlinearity and ill-posedness. Different regularization techniques have been applied to tackle the ill-posedness of MWI problems [12]. Incorporating prior knowledge into the inverse scattering model, such as known physical/electrical object features, reduces the ambiguity of the reconstruction and moderates the occurrence of false solutions. If the unknown object is electrically small, one of the most effective regularization tools is sparse processing or compressive sensing (CS) [20, 21]. CS efficiently includes the prior knowledge about the target through the sparsity constraint or l1 regularization. The resulting images are sharp, with limited artifacts, and robust against noise. However, if the unknown object is electrically large, significant portions of the target shape remain unrecovered. Thus, the obtained images often appear as collections of disconnected points unless some additional knowledge about targets is available [24].This paper introduces a novel qualitative MWI algorithm, which achieves high spatial resolution by taking advantage of sophisticated electromagnetic modeling and the principles of sparse processing. We use the fact that with respect to the radial direction, the scattered field can be expressed as a superposition of a transverse-magnetic polarized wave (denoted as TMr) and a transverse-electric polarized wave (denoted as TEr) [25]. Both TMr and TEr fields are then written as a sum of spherical waves. Each spherical wave is determined by a pair of indices (n, m), related to radial and angular variation of the electromagnetic field, respectively. The sources of these waves are collocated in the target’s center.However, we further show that the scattered field can be represented as a sum of spherical waves of a single order (n, m), if the wave sources are sufficiently separated in space. This theoretical result makes the foundation of the proposed imaging algorithm. We compute different target images by using wave sources of a single order distributed on a uniform grid. Each order (harmonic) produces a unique image, and we obtain the final target image by superimposing these partial images.A further contribution of this work is the recognition of the so-called “dominant” harmonics, which have a major role in target imaging. Namely, we show that in certain cases, the dominant harmonics are sufficient for target reconstruction, even in the case of complex-shaped targets. Using the dominant harmonics remarkably simplifies the optimization process and reduces the computational time at the expenses of negligible loss in the reconstruction accuracy. Nevertheless, we also show that in the cases of electrically small targets, other harmonics have an indispensable role in improving the spatial resolution of the imaging results.Another novelty is the selection of leading orders, which is a nontrivial question in 3D space. Starting from the results of the theory of degrees of freedom [26, 27] and the condition m≤n, we determine the maximal order nmax, which can be used for reconstruction. We also show that images obtained with n>nmax produce pixels lying outside the target’s convex envelope so that they can be safely removed from the image.The present work completes and generalizes the previous work done by the authors, in which we dealt with the 2D case (i.e., targets infinitely extended along one direction) through dictionaries made with cylindrical basis functions for both TM [28] and TE [29] scenarios. Since cylindrical bases functions are determined by a single index, n, determining the relevant harmonics in the 2D space is a much simpler task than in 3D. Preliminary work on 3D reconstruction was presented in [30]. In that case, a simplified measurement model was considered, which assumed only vertically polarized antennas and, consequently, only vertically oriented low-order multipoles as the sources of the scattered field. Here, instead, we assume arbitrarily polarized measurement array and spherical harmonics of an arbitrary order as the sources of the scattered field, thus achieving an utterly general model.Another outcome of this work is the normalization of the system matrix, which allows for a straightforward selection of the regularization coefficient. Typically, the regularization coefficient is obtained by plotting the L-curve and selecting the value associated with the knee of the curve [31]. However, this means solving the minimization function several times, once for each value of the regularization parameter, which can be numerically consuming. The proposed normalization translates the regularization parameter in the range 0.1≤γ≤1, with the knee of the L-curve roughly being close to 0.5.Finally, it is worth noting that most CS/sparse microwave imaging algorithms in the literature deal with 2D problems. Examples of 2D algorithms that combine the first-order Born approximation and CS are [32, 33]. In contrast, we propose an algorithm that does not require Born approximation’s validity as it takes advantage of the linear relationship holding between the scattered field and the induced sources. Moreover, the algorithm is applicable for 3D imaging of arbitrary dielectric/metallic targets. In [34, 35], the authors use the sparsity constraint to regularize the distorted Born iterative method (DBIM). In [36], the standard 3D CS-based algorithm was utilized for the breast-cancer localization. However, the measurement model described in [36] uses jointly all the measurements, which is suitable only for small target imaging, such as early-stage cancer. In contrast, we divide the measurements into groups that we process separately and combine the obtained results at the end. In this way, it is possible to reconstruct complex-shaped targets and augment the reconstructed portion of the target, which is critical in sparse-based algorithms. In [37], the CS approach was applied to 2D extended targets with simple and concave cross sections.The proposed approach was assessed on simulated data involving various concave-shaped objects. Due to multiple scattering, imaging of such objects is significantly more challenging than that of simple, smooth objects. In most cases, standard qualitative algorithms retrieve a convex envelope, ignoring the concave parts of the target shape. Distinguishing between convex and concave targets is vital in many applications, such as distinguishing benign tumors from malignant ones [5].In addition, the performance of the proposed algorithm was compared to that of the gold standard methods in qualitative imaging, such as LSM and TSVD. We showed that the proposed algorithm outperforms both LSM and TSVD when electrically small targets are imaged. In addition, it requires fewer measurements and is less sensitive to noise. However, if the number of sensors is large, LSM may be the preferred choice due to short computational time and ease of implementation. Finally, we compared the imaging results of the presented algorithm and those obtained by another sparse processing algorithm [38]. Presented results show that the proposed algorithm restores better electrically small targets at the cost of more complex implementation. Both algorithms can retrieve complex shapes, yet the restored portions of the target were significantly larger when the spherical harmonics were utilized as opposed to the case when the standard sparse processing algorithm was applied. Consequently, the point-like nature of the sparse microwave images was considerably suppressed.The paper is organized as follows. The analytical background is given in Section2. The numerical model is described in Section 3. In Section 4, we describe the methods used for comparison. The results obtained for various scenarios are reported in Section 5. In Section 6, we give comparison of our method with other methods. Some final remarks are given in the Conclusion section.


## 2. Analytical Considerations

### 2.1. Scattered Field Decomposition

We consider a 3D scatterer, depicted, without loss of generality, as a yellow sphere in Figure1(a), located in a homogeneous nonmagnetic medium of permittivity ε. The scatterer is surrounded by an array of sensors, only two of which are shown in Figure 1(a). The field scattered from the target is a superposition of the transverse magnetic (TMr) and the transverse electromagnetic components (TEr) [25]. As detailed in the Appendix, we decompose both TMr and TEr components using spherical harmonics. Below, we present only the final expressions for the θ and ϕ components of the electric field:(1)EθTMr,θ,ϕ=∑n=1∞∑m=0ngnβrpn,mθDn,mcosmϕ+En,msinmϕ,(2)EφTMr,θ,ϕ=∑n=1∞∑m=1ngnβrqn,mθDn,mcosmϕ+En,msinmϕ,(3)EθTEr,θ,ϕ=∑n=1∞∑m=1nfnβrqn,mθZn,mcosmϕ+Yn,msinmϕ,(4)EφTEr,θ,ϕ=∑n=1∞∑m=0nfnβrpn,mθZn,mcosmϕ+Yn,msinmϕ,with(5)pn,mθ=mcosθsinθPnmcosθ+Pnm+1cosθ,qn,mθ=mPnmcosθsinθ,gnr=πβ2rβHn−1/22βr−nrHn+1/22βr,fnr=πβ2rHn+1/22βr.Figure 1
(a) Multipoles of various orders, depicted as small spheres collocated in the scatterer’s center. (b) Computation of the radiation pattern of a multipole of the order (m, n). (c) Spatial array of multipoles of the order (m, n).
(a)(b)(c)Herer is the receiver’s position vector; r,θ,φ are the spherical coordinates of r; Dn,m, En,m, Zn,m, and Yn,m are the expansion coefficients; Hn2βr is the cylindrical Hankel function of the second kind and order n; and Pnmcosθ is the associated Legendre function of the first kind, degree n, and order m.We assume that the antenna array is in the far field with respect to the scatterer, and, thus, the radial component of the electric field is negligible. Moreover, using the asymptotic expressions for the cylindrical Hankel functions, we have(6)gnr⟶βr⟶∞βjne−jβrr,fnr⟶βr⟶∞βjn+1e−jβrr.From the properties of the associated Legendre functions [25], it follows that n≥1, 0≤m≤n in (1) and (4) and n≥1, 1≤m≤n in (2) and (3). The number of the relevant terms in (1)–(4) is determined by the theory of degrees of freedom (DoF) of EM fields [26, 27]. This theory states that the scattered field of a bounded target, measured on a spherical domain enclosing it, can be approximated by βa2 independent pieces of information, where β=ωεμ is the phase coefficient and a is the radius of the minimum sphere enclosing the target. Thus, we can expect that the upper bound for relevant terms per each polarization is approximately βa2/2. Therefore, we can calculate the maximal required order n from the condition(7)∑n=1nmax∑m=0n1=nmaxnmax+32≤βa22.In Figure1(a), we visualize the sources of the spherical waves (multipoles) as small spheres collocated in the scatterer’s center (origin), where each color refers to one multipole order.


### 2.2. Alternative Scattered Field Decomposition

Let us consider wave sources having the same polarization and order (m, n), which are continuously distributed on a sphere of radius a. Figure 1(b) shows an example of such a source located at r′. The radial component of the potential generated by the source at the field point r is(8)Ar∝H^n2βRPnmcosΘcosmΦr0,R=r−r′,where R,Θ,Φ are the spherical coordinates of R, a,θ′,φ′ are the spherical coordinates of r′, and r0=R/R. In the far field, we have(9)H^n2βr⟶βr⟶∞jn+1e−jβr,Θ,Φ≈θ,φ,R≈r−ir⋅ir′,r0≈ir,where ir′=r′/r′ and ir=r/r. Thus, the potential is approximated as(10)Ar∝e−jβre+jβir⋅ir′Pnmcosθcosmφir,and the total potential is given by(11)Ar≈e−jβrPnmcosθcosmϕ∫0π∫02πgθ′,ϕ′exp+jβair⋅ir′dθ′dϕ′,where the function gθ′,φ′ determines the source magnitudes. Using the Taylor series, we obtain for the exponential in (11)(12)exp+jβair⋅ir′=1+∑k=1∞jβair⋅ir′kk!.Forβa≪1, only the first term in the series is important (i.e., exp+jβair⋅ir′≈1). Thus, the resulting potential is(13)Ar≈Ae−jβrPnmcosθcosmϕ,where A is the coefficient, which depends on gθ′,ϕ′. As expected, the sources of the order (m, n), located on an electrically small sphere, produce the electromagnetic wave of the same order (m, n). When a increases, more terms in (13) become significant. Thus, if we keep one more term, the resulting potential becomes(14)Ar≈e−jβrPnmcosθcosmϕ∫0π∫02πgθ′,ϕ′1+jβair⋅ir′dθ′dϕ′=e−jβrPnmcosθcosmϕ∫0π∫02πgθ′,ϕ′1+jβacosθcosθ′+sinθsinθ′cosφ−φ′dθ′dϕ′.After the integration and the application of recurrence formulas for associated Legendre polynomials, we can identify two groups of terms in the resulting potential. In the first group, we have the following terms:e−jβrPnmcosθcosmϕ, e−jβrPn+1m−1cosθcosm−1φ, e−jβrPn+1m−1cosθsinm−1φ, e−jβrPn+1m+1cosθcosm+1φ, and e−jβrPn+1m−1cosθsinm−1φ, which are the sources of the waves of the orders m−1,n−1 and m+1,n+1. In the second group, we have the following terms: e−jβrPnm+1cosθcosmφ, e−jβrPnm−1cosθcosmφ, e−jβrPn+1m−1cosθcosm+1φ, and e−jβrPn+1m+1cosθcosm−1φ, which do not satisfy the Helmholtz equation (A.3). If we further increase radius a, we have to keep additional terms in (16). Consequently, by adjusting the value of a, we can generate the spherical waves of arbitrary orders. Thus, we can model the field scattered from the target by the spatial array of multipoles of one order (e.g., (m, n)), as illustrated in Figure 1(c). The approximate expressions for the electric field in this case are(15)EθTMr,θ,ϕ≈∑l=1LgnβRlpn,mΘlan,mcosmΦl+bn,msinmΦl,(16)EϕTMr,θ,ϕ≈∑l=1LgnβRlqn,mΘlan,mcosmΦl+bn,msinmΦl,(17)EθTEr,θ,ϕ≈∑l=1LfnβRlqn,mΘlcn,mcosmΦl+dn,msinmΦl,(18)EϕTEr,θ,ϕ≈∑l=1LfnβRlpn,mΘlcn,mcosmΦl+dn,msinmΦl,where an,m, bn,m, cn,m, and dn,m are the multipole coefficients. Further, Rl=r−tl, where tl is the position vector of the lth multipole, Rl,Θl,Φl are the spherical coordinates of Rl, and L is the total number of the multipoles.The expansion coefficients must be selected in such a way as to boost the components belonging to the spectrum of the scattered field and to suppress the unwanted components. In this case, the multipoles will be distributed in the target’s interior. Thus, by reconstructing the locations of the multipoles, we can estimate the target’s shape. We note that the multipoles with orders (m, n) close to border values produce more unwanted components. Consequently, the sums (15)–(18) may diverge, yielding erroneous reconstructions. However, we will show later how to select optimal multipole orders for imaging.


## 2.1. Scattered Field Decomposition

We consider a 3D scatterer, depicted, without loss of generality, as a yellow sphere in Figure1(a), located in a homogeneous nonmagnetic medium of permittivity ε. The scatterer is surrounded by an array of sensors, only two of which are shown in Figure 1(a). The field scattered from the target is a superposition of the transverse magnetic (TMr) and the transverse electromagnetic components (TEr) [25]. As detailed in the Appendix, we decompose both TMr and TEr components using spherical harmonics. Below, we present only the final expressions for the θ and ϕ components of the electric field:(1)EθTMr,θ,ϕ=∑n=1∞∑m=0ngnβrpn,mθDn,mcosmϕ+En,msinmϕ,(2)EφTMr,θ,ϕ=∑n=1∞∑m=1ngnβrqn,mθDn,mcosmϕ+En,msinmϕ,(3)EθTEr,θ,ϕ=∑n=1∞∑m=1nfnβrqn,mθZn,mcosmϕ+Yn,msinmϕ,(4)EφTEr,θ,ϕ=∑n=1∞∑m=0nfnβrpn,mθZn,mcosmϕ+Yn,msinmϕ,with(5)pn,mθ=mcosθsinθPnmcosθ+Pnm+1cosθ,qn,mθ=mPnmcosθsinθ,gnr=πβ2rβHn−1/22βr−nrHn+1/22βr,fnr=πβ2rHn+1/22βr.Figure 1
(a) Multipoles of various orders, depicted as small spheres collocated in the scatterer’s center. (b) Computation of the radiation pattern of a multipole of the order (m, n). (c) Spatial array of multipoles of the order (m, n).
(a)(b)(c)Herer is the receiver’s position vector; r,θ,φ are the spherical coordinates of r; Dn,m, En,m, Zn,m, and Yn,m are the expansion coefficients; Hn2βr is the cylindrical Hankel function of the second kind and order n; and Pnmcosθ is the associated Legendre function of the first kind, degree n, and order m.We assume that the antenna array is in the far field with respect to the scatterer, and, thus, the radial component of the electric field is negligible. Moreover, using the asymptotic expressions for the cylindrical Hankel functions, we have(6)gnr⟶βr⟶∞βjne−jβrr,fnr⟶βr⟶∞βjn+1e−jβrr.From the properties of the associated Legendre functions [25], it follows that n≥1, 0≤m≤n in (1) and (4) and n≥1, 1≤m≤n in (2) and (3). The number of the relevant terms in (1)–(4) is determined by the theory of degrees of freedom (DoF) of EM fields [26, 27]. This theory states that the scattered field of a bounded target, measured on a spherical domain enclosing it, can be approximated by βa2 independent pieces of information, where β=ωεμ is the phase coefficient and a is the radius of the minimum sphere enclosing the target. Thus, we can expect that the upper bound for relevant terms per each polarization is approximately βa2/2. Therefore, we can calculate the maximal required order n from the condition(7)∑n=1nmax∑m=0n1=nmaxnmax+32≤βa22.In Figure1(a), we visualize the sources of the spherical waves (multipoles) as small spheres collocated in the scatterer’s center (origin), where each color refers to one multipole order.


## 2.2. Alternative Scattered Field Decomposition

Let us consider wave sources having the same polarization and order (m, n), which are continuously distributed on a sphere of radius a. Figure 1(b) shows an example of such a source located at r′. The radial component of the potential generated by the source at the field point r is(8)Ar∝H^n2βRPnmcosΘcosmΦr0,R=r−r′,where R,Θ,Φ are the spherical coordinates of R, a,θ′,φ′ are the spherical coordinates of r′, and r0=R/R. In the far field, we have(9)H^n2βr⟶βr⟶∞jn+1e−jβr,Θ,Φ≈θ,φ,R≈r−ir⋅ir′,r0≈ir,where ir′=r′/r′ and ir=r/r. Thus, the potential is approximated as(10)Ar∝e−jβre+jβir⋅ir′Pnmcosθcosmφir,and the total potential is given by(11)Ar≈e−jβrPnmcosθcosmϕ∫0π∫02πgθ′,ϕ′exp+jβair⋅ir′dθ′dϕ′,where the function gθ′,φ′ determines the source magnitudes. Using the Taylor series, we obtain for the exponential in (11)(12)exp+jβair⋅ir′=1+∑k=1∞jβair⋅ir′kk!.Forβa≪1, only the first term in the series is important (i.e., exp+jβair⋅ir′≈1). Thus, the resulting potential is(13)Ar≈Ae−jβrPnmcosθcosmϕ,where A is the coefficient, which depends on gθ′,ϕ′. As expected, the sources of the order (m, n), located on an electrically small sphere, produce the electromagnetic wave of the same order (m, n). When a increases, more terms in (13) become significant. Thus, if we keep one more term, the resulting potential becomes(14)Ar≈e−jβrPnmcosθcosmϕ∫0π∫02πgθ′,ϕ′1+jβair⋅ir′dθ′dϕ′=e−jβrPnmcosθcosmϕ∫0π∫02πgθ′,ϕ′1+jβacosθcosθ′+sinθsinθ′cosφ−φ′dθ′dϕ′.After the integration and the application of recurrence formulas for associated Legendre polynomials, we can identify two groups of terms in the resulting potential. In the first group, we have the following terms:e−jβrPnmcosθcosmϕ, e−jβrPn+1m−1cosθcosm−1φ, e−jβrPn+1m−1cosθsinm−1φ, e−jβrPn+1m+1cosθcosm+1φ, and e−jβrPn+1m−1cosθsinm−1φ, which are the sources of the waves of the orders m−1,n−1 and m+1,n+1. In the second group, we have the following terms: e−jβrPnm+1cosθcosmφ, e−jβrPnm−1cosθcosmφ, e−jβrPn+1m−1cosθcosm+1φ, and e−jβrPn+1m+1cosθcosm−1φ, which do not satisfy the Helmholtz equation (A.3). If we further increase radius a, we have to keep additional terms in (16). Consequently, by adjusting the value of a, we can generate the spherical waves of arbitrary orders. Thus, we can model the field scattered from the target by the spatial array of multipoles of one order (e.g., (m, n)), as illustrated in Figure 1(c). The approximate expressions for the electric field in this case are(15)EθTMr,θ,ϕ≈∑l=1LgnβRlpn,mΘlan,mcosmΦl+bn,msinmΦl,(16)EϕTMr,θ,ϕ≈∑l=1LgnβRlqn,mΘlan,mcosmΦl+bn,msinmΦl,(17)EθTEr,θ,ϕ≈∑l=1LfnβRlqn,mΘlcn,mcosmΦl+dn,msinmΦl,(18)EϕTEr,θ,ϕ≈∑l=1LfnβRlpn,mΘlcn,mcosmΦl+dn,msinmΦl,where an,m, bn,m, cn,m, and dn,m are the multipole coefficients. Further, Rl=r−tl, where tl is the position vector of the lth multipole, Rl,Θl,Φl are the spherical coordinates of Rl, and L is the total number of the multipoles.The expansion coefficients must be selected in such a way as to boost the components belonging to the spectrum of the scattered field and to suppress the unwanted components. In this case, the multipoles will be distributed in the target’s interior. Thus, by reconstructing the locations of the multipoles, we can estimate the target’s shape. We note that the multipoles with orders (m, n) close to border values produce more unwanted components. Consequently, the sums (15)–(18) may diverge, yielding erroneous reconstructions. However, we will show later how to select optimal multipole orders for imaging.


## 3. Measurement Models

### 3.1. Inverse Model

We assume a spherical array of half-wavelength dipoles, whose positions are defined by the vectorsrk, k=1,…,M, where M is total number of dipoles. The array has two configurations: one with the dipoles parallel to the unit vector iθ, and the other with the dipoles parallel to the unit vector iϕ of the spherical coordinate system.Imaging is performed on uniform grids, as illustrated in Figure2. At each grid node, there is a multipole of the selected order and polarization. The locations of the grid nodes are defined by the position vectors, tl, l=1,…,L, where L is the overall number of the grid nodes. The measurement model follows from (15)–(18) and it relates the multipole coefficients and the measured scattered field. When the ith antenna is transmitting, we have(19)ei=Wn,mkn,mi.Figure 2
Simplified illustration of the grid used for the sampling of the imaging domain.Hereei is the measurement vector, Wn,m is the system matrix, and kn,mi is the unknown vector whose elements are the multipole coefficients. The system matrix consists of two dictionaries related to the cosine and sine expansion terms:(20)Wn,m=Cn,mSn,m.Similarly, vectorkn,mi consists of two parts:(21)kn,mi=cn,misn,mi.Herecn,mi and sn,mi are the unknown vectors corresponding to the dictionaries Cn,m and Sn,m, respectively. Depending on the array configuration and multipole polarization, Wn,m, kn,mi, and ei have different forms. Below, we list the possible combinations.Array withθ-oriented dipoles, TMr polarized multipoles:(22)Cn,m=gnβr11pn,mθ11cosmϕ11⋯gnβr1Lpn,mθ1Lcosmϕ1L⋮⋱⋮gnβrM1pn,mθM1cosmϕM1⋯gnβrMLpn,mθMLcosmϕML,Sn,m=gnβr11pn,mθ11sinmϕ11⋯gnβr1Lpn,mθ1Lsinmϕ1L⋮⋱⋮gnβrM1pn,mθM1sinmϕM1⋯gnβrMLpn,mθMLsinmϕML.(23)ei=Es,θr1,ri⋮Es,θrM,ri,cn,mi=a1⋮aL,sn,mi=b1⋮bL.Array withθ-oriented dipoles, TEr polarized multipoles:(24)Cn,m=fnβr11qn,mθ11cosmϕ11⋯fnβr1Lqn,mθ1Lcosmϕ1L⋮⋱⋮fnβrM1qn,mθM1cosmϕM1⋯fnβrMLqn,mθMLcosmϕML,Sn,m=fnβr11qn,mθ11sinmϕ11⋯fnβr1Lqn,mθ1Lsinmϕ1L⋮⋱⋮fnβrM1qn,mθM1sinmϕM1⋯fnβrMLqn,mθMLsinmϕML,(25)ei=Es,θr1,ri⋮Es,θrM,ri,cn,mi=c1⋮cL,sn,mi=d1⋮dL.Array withφ-oriented dipoles, TMr polarized multipoles:(26)Cn,m=gnβr11qn,mθ11cosmϕ11⋯gnβr1Lqn,mθ1Lcosmϕ1L⋮⋱⋮gnβrM1qn,mθM1cosmϕM1⋯gnβrMLqn,mθMLcosmϕML,Sn,m=gnβr11qn,mθ11sinmϕ11⋯gnβr1Lqn,mθ1Lsinmϕ1L⋮⋱⋮gnβrM1qn,mθM1sinmϕM1⋯gnβrMLqn,mθMLsinmϕML,(27)ei=Es,ϕr1,ri⋮Es,ϕrM,ri,cn,mi=a1⋮aL,sn,mi=b1⋮bL.Array withφ-oriented dipoles, TEr polarized multipoles:(28)Cn,m=fnβr11pn,mθ11cosmϕ11⋯fnβr1Lpn,mθ1Lcosmϕ1L⋮⋱⋮fnβrM1pn,mθM1cosmϕM1⋯fnβrMLpn,mθMLcosmϕML,Sn,m=fnβr11pn,mθ11sinmϕ11⋯fnβr1Lpn,mθ1Lsinmϕ1L⋮⋱⋮fnβrM1pn,mθM1sinmϕM1⋯fnβrMLpn,mθMLsinmϕML,(29)ei=Es,ϕr1,ri⋮Es,ϕrM,ri,cn,mi=c1⋮cL,sn,mi=d1⋮dL.In the above formulas,Rkl=rk−ti is the distance between the kth antenna and lth grid node; and rkl,θkl,ϕkl are the spherical coordinates of the vector Rkl. For better clarity, we omit the dependence of coefficients al, bl, cl, and dl, l=1,…,L, on the order of spherical harmonics (n, m) in (23), (25), (27), and (24).In the algorithm’s implementation, we compute the measurement vector in all four cases as(30)ei=Δs1,iΔs2,i⋮ΔsM,i,Δsj,i=sj,i−sj,i0,j=1,…,M.Heresj,i and sji0 are the transmission coefficients between the ith and jth antenna computed in the presence and the absence of the target, respectively. In addition, we define a vector whose elements are the square norms of the coefficients cn,mi and sn,mi at each grid node. For example, for the coefficients defined in (23) and (27), we have(31)pn,mi=a12+b12…aL2+bL2T.In the same way, we can definepn,mi for the coefficients defined in (25) and (29).


### 3.2. Sparse Processing

The goal of the algorithm is to minimize the number of the grid nodes with significant coefficients, that is, to reduce the number of nonzero elements of (31). Hence, the minimization function is(32)k^n,mi=minkn,miei−Wn,mkn,mi22+γpn,mi1,where γ is the regularization parameter that balances between the data fidelity (represented by the l2 norm of the error) and the solution sparsity (represented by the l1 norm of the solution vector). For this purpose, we utilize the convex programming package (CVX) [38].Typically, the regularization coefficientγ is obtained by plotting the L-curve and selecting the value associated with the knee of the curve [31]. That means solving (32) one time per each value of the regularization parameter, which is numerically consuming. However, we found that a simple normalization of the system matrix can significantly shorten the selection of the regularization coefficient. Namely, when we divide each column of (20) by its square norm, the range of useful values of the regularization coefficient falls in the range of 0.1<γ<1. Moreover, the knee of the L-curve corresponds approximately to γ=0.5.To suppress the noise, we process together data obtained from a few adjacent sensors in the array. The number of the combined transmitting antennas is denoted withP. For example, for P=2, the minimization function reads(33)k^n,mi=minkn,m,ieiei+1−Wn,mWn,mkn,mi22+γpn,mi1.The choice ofP is a trade-off between the image resolution and robustness to noise. Namely, large values for P are favorable in low SNR scenarios, but they produce smooth reconstructions (i.e., cause loosing nonconvex features of targets’ surfaces).


### 3.3. Image Construction

To obtain the complete image, we superimpose the partial results obtained for groups of incidences. In addition, we add results obtained for two polarizations. Hence, the image obtained using the multipoles of the order (m, n) equals(34)Im,nl=∑ipn,mil︸TMr+∑ipn,mil︸TEr,where l is the pixel index and I(l) is the corresponding image value. In the further analysis, we study images obtained using multipoles of different orders. We show that only a few orders are relevant, and that by summing the images produced by different orders the quality of the reconstruction is improved.


### 3.4. Digital Filtering

Estimation with higher-order multipoles may produce artifacts. To obtain clear images, we use a convex hull as a digital filter. To compute the convex hull, we adopted the approach described in [28, 29] for the 3D scenario. We first run the algorithm using the dominant harmonic (which we describe in detail later). The regularization parameter γ needs to have a high value associated with the vertical part of the L-curve. In our case, this is equivalent to γ>1(the upper limit of γ is insignificant.) The convex hull gets its final form after applying the binarization procedure in which all the pixels whose values are above a threshold are set to one, and the pixels whose values are below the threshold are set to zero.


### 3.5. Quality Measures

Typically, 3D reconstruction is performed in several 2D cuts, which are either parallel or orthogonal planes. In this work, we computed targets’ images in the two cuts coinciding withxOy and xOz planes of the Cartesian coordinate system. To quantify the precision of the method and its reconstruction abilities, we introduce two metrics. First, we define η1 parameter as(35)η1=NinNtot,where Nin denotes the number of pixels obtained inside the object’s contour and Ntot denotes the total number of the obtained pixels. In other words, η1 parameter appraises the method precision by measuring how many of the obtained pixels “fall” inside the object contour. The second metric, denoted by η2, assesses the ability of the method to reconstruct the object shape and is defined as(36)η2=NinNobj,where Nobj denotes the number of pixels inside the object’s actual contour and is fixed for a chosen object shape. In the calculations, we considered as relevant all pixels with the intensity higher than 5% of the maximal value. Both metrics are calculated for each cut separately. Ideally, they should be close to unity.


## 3.1. Inverse Model

We assume a spherical array of half-wavelength dipoles, whose positions are defined by the vectorsrk, k=1,…,M, where M is total number of dipoles. The array has two configurations: one with the dipoles parallel to the unit vector iθ, and the other with the dipoles parallel to the unit vector iϕ of the spherical coordinate system.Imaging is performed on uniform grids, as illustrated in Figure2. At each grid node, there is a multipole of the selected order and polarization. The locations of the grid nodes are defined by the position vectors, tl, l=1,…,L, where L is the overall number of the grid nodes. The measurement model follows from (15)–(18) and it relates the multipole coefficients and the measured scattered field. When the ith antenna is transmitting, we have(19)ei=Wn,mkn,mi.Figure 2
Simplified illustration of the grid used for the sampling of the imaging domain.Hereei is the measurement vector, Wn,m is the system matrix, and kn,mi is the unknown vector whose elements are the multipole coefficients. The system matrix consists of two dictionaries related to the cosine and sine expansion terms:(20)Wn,m=Cn,mSn,m.Similarly, vectorkn,mi consists of two parts:(21)kn,mi=cn,misn,mi.Herecn,mi and sn,mi are the unknown vectors corresponding to the dictionaries Cn,m and Sn,m, respectively. Depending on the array configuration and multipole polarization, Wn,m, kn,mi, and ei have different forms. Below, we list the possible combinations.Array withθ-oriented dipoles, TMr polarized multipoles:(22)Cn,m=gnβr11pn,mθ11cosmϕ11⋯gnβr1Lpn,mθ1Lcosmϕ1L⋮⋱⋮gnβrM1pn,mθM1cosmϕM1⋯gnβrMLpn,mθMLcosmϕML,Sn,m=gnβr11pn,mθ11sinmϕ11⋯gnβr1Lpn,mθ1Lsinmϕ1L⋮⋱⋮gnβrM1pn,mθM1sinmϕM1⋯gnβrMLpn,mθMLsinmϕML.(23)ei=Es,θr1,ri⋮Es,θrM,ri,cn,mi=a1⋮aL,sn,mi=b1⋮bL.Array withθ-oriented dipoles, TEr polarized multipoles:(24)Cn,m=fnβr11qn,mθ11cosmϕ11⋯fnβr1Lqn,mθ1Lcosmϕ1L⋮⋱⋮fnβrM1qn,mθM1cosmϕM1⋯fnβrMLqn,mθMLcosmϕML,Sn,m=fnβr11qn,mθ11sinmϕ11⋯fnβr1Lqn,mθ1Lsinmϕ1L⋮⋱⋮fnβrM1qn,mθM1sinmϕM1⋯fnβrMLqn,mθMLsinmϕML,(25)ei=Es,θr1,ri⋮Es,θrM,ri,cn,mi=c1⋮cL,sn,mi=d1⋮dL.Array withφ-oriented dipoles, TMr polarized multipoles:(26)Cn,m=gnβr11qn,mθ11cosmϕ11⋯gnβr1Lqn,mθ1Lcosmϕ1L⋮⋱⋮gnβrM1qn,mθM1cosmϕM1⋯gnβrMLqn,mθMLcosmϕML,Sn,m=gnβr11qn,mθ11sinmϕ11⋯gnβr1Lqn,mθ1Lsinmϕ1L⋮⋱⋮gnβrM1qn,mθM1sinmϕM1⋯gnβrMLqn,mθMLsinmϕML,(27)ei=Es,ϕr1,ri⋮Es,ϕrM,ri,cn,mi=a1⋮aL,sn,mi=b1⋮bL.Array withφ-oriented dipoles, TEr polarized multipoles:(28)Cn,m=fnβr11pn,mθ11cosmϕ11⋯fnβr1Lpn,mθ1Lcosmϕ1L⋮⋱⋮fnβrM1pn,mθM1cosmϕM1⋯fnβrMLpn,mθMLcosmϕML,Sn,m=fnβr11pn,mθ11sinmϕ11⋯fnβr1Lpn,mθ1Lsinmϕ1L⋮⋱⋮fnβrM1pn,mθM1sinmϕM1⋯fnβrMLpn,mθMLsinmϕML,(29)ei=Es,ϕr1,ri⋮Es,ϕrM,ri,cn,mi=c1⋮cL,sn,mi=d1⋮dL.In the above formulas,Rkl=rk−ti is the distance between the kth antenna and lth grid node; and rkl,θkl,ϕkl are the spherical coordinates of the vector Rkl. For better clarity, we omit the dependence of coefficients al, bl, cl, and dl, l=1,…,L, on the order of spherical harmonics (n, m) in (23), (25), (27), and (24).In the algorithm’s implementation, we compute the measurement vector in all four cases as(30)ei=Δs1,iΔs2,i⋮ΔsM,i,Δsj,i=sj,i−sj,i0,j=1,…,M.Heresj,i and sji0 are the transmission coefficients between the ith and jth antenna computed in the presence and the absence of the target, respectively. In addition, we define a vector whose elements are the square norms of the coefficients cn,mi and sn,mi at each grid node. For example, for the coefficients defined in (23) and (27), we have(31)pn,mi=a12+b12…aL2+bL2T.In the same way, we can definepn,mi for the coefficients defined in (25) and (29).


## 3.2. Sparse Processing

The goal of the algorithm is to minimize the number of the grid nodes with significant coefficients, that is, to reduce the number of nonzero elements of (31). Hence, the minimization function is(32)k^n,mi=minkn,miei−Wn,mkn,mi22+γpn,mi1,where γ is the regularization parameter that balances between the data fidelity (represented by the l2 norm of the error) and the solution sparsity (represented by the l1 norm of the solution vector). For this purpose, we utilize the convex programming package (CVX) [38].Typically, the regularization coefficientγ is obtained by plotting the L-curve and selecting the value associated with the knee of the curve [31]. That means solving (32) one time per each value of the regularization parameter, which is numerically consuming. However, we found that a simple normalization of the system matrix can significantly shorten the selection of the regularization coefficient. Namely, when we divide each column of (20) by its square norm, the range of useful values of the regularization coefficient falls in the range of 0.1<γ<1. Moreover, the knee of the L-curve corresponds approximately to γ=0.5.To suppress the noise, we process together data obtained from a few adjacent sensors in the array. The number of the combined transmitting antennas is denoted withP. For example, for P=2, the minimization function reads(33)k^n,mi=minkn,m,ieiei+1−Wn,mWn,mkn,mi22+γpn,mi1.The choice ofP is a trade-off between the image resolution and robustness to noise. Namely, large values for P are favorable in low SNR scenarios, but they produce smooth reconstructions (i.e., cause loosing nonconvex features of targets’ surfaces).


## 3.3. Image Construction

To obtain the complete image, we superimpose the partial results obtained for groups of incidences. In addition, we add results obtained for two polarizations. Hence, the image obtained using the multipoles of the order (m, n) equals(34)Im,nl=∑ipn,mil︸TMr+∑ipn,mil︸TEr,where l is the pixel index and I(l) is the corresponding image value. In the further analysis, we study images obtained using multipoles of different orders. We show that only a few orders are relevant, and that by summing the images produced by different orders the quality of the reconstruction is improved.


## 3.4. Digital Filtering

Estimation with higher-order multipoles may produce artifacts. To obtain clear images, we use a convex hull as a digital filter. To compute the convex hull, we adopted the approach described in [28, 29] for the 3D scenario. We first run the algorithm using the dominant harmonic (which we describe in detail later). The regularization parameter γ needs to have a high value associated with the vertical part of the L-curve. In our case, this is equivalent to γ>1(the upper limit of γ is insignificant.) The convex hull gets its final form after applying the binarization procedure in which all the pixels whose values are above a threshold are set to one, and the pixels whose values are below the threshold are set to zero.


## 3.5. Quality Measures

Typically, 3D reconstruction is performed in several 2D cuts, which are either parallel or orthogonal planes. In this work, we computed targets’ images in the two cuts coinciding withxOy and xOz planes of the Cartesian coordinate system. To quantify the precision of the method and its reconstruction abilities, we introduce two metrics. First, we define η1 parameter as(35)η1=NinNtot,where Nin denotes the number of pixels obtained inside the object’s contour and Ntot denotes the total number of the obtained pixels. In other words, η1 parameter appraises the method precision by measuring how many of the obtained pixels “fall” inside the object contour. The second metric, denoted by η2, assesses the ability of the method to reconstruct the object shape and is defined as(36)η2=NinNobj,where Nobj denotes the number of pixels inside the object’s actual contour and is fixed for a chosen object shape. In the calculations, we considered as relevant all pixels with the intensity higher than 5% of the maximal value. Both metrics are calculated for each cut separately. Ideally, they should be close to unity.


## 4. Comparison Methods

This section gives a brief description of the methods used to compare the reconstruction results. Since the proposed method is qualitative, we used LSM [15, 16] and TSVD [18, 19], which are the golden standard inversion methods adopted for qualitative real-time imaging. In addition, we compared the images obtained by the proposed algorithm, further denoted as SHSP (spherical harmonics sparse processing), and another sparse processing algorithm [39], denoted as SSP (standard sparse processing).


### 4.1. Truncated Singular Value Decomposition (TSVD)

The TSVD measurement model follows from the linearized scattering equation [12]. The target image is obtained from the linear system of equations(37)e=Lf,L=Einct1,r1⋅Einct1,r1…EinctL,r1⋅EinctL,r1Einct1,r1⋅Einct1,r2⋮EinctL,r1⋅EinctL,r2⋮⋮⋮Einct1,rM⋅Einct1,rM…EinctL,rM⋅EinctL,rM,e=Δs1,1Δs1,2⋮ΔsM,M,f=f1f2⋮fL.HereL is the system matrix, e is the known measurement vector, and f is the unknown vector (target image). The TSVD solution of the system is(38)f=∑i=1P1σiuiH⋅evi,where ui, vi are the singular vectors of matrix L, σi are the corresponding singular values, and P is the truncation index, obtained from the condition −20log10σP=TH, where TH is the adopted threshold.


### 4.2. Linear Sampling Method (LSM)

In this method, a system of equations is solved for each search point and each Cartesian component of the incident field:(39)Lfs=gstl,L=Δs1,1⋯Δs1,М⋮⋱⋮ΔsM,1⋯ΔsM,M,gstl=Einc,stl,r1⋮Einc,stl,rМ,fstl=f1f2⋮fL.HereL is the system matrix, gs is the known vector whose elements are s∈x,y,z components of the incident electric field, produced by the antenna array at tl, and fs is the unknown vector. The solution obtained by means of Tikhonov regularization is(40)fstl=∑i=1Mσiσi2+α2uiH⋅gtlvi,where α is the regularization coefficient, typically selected as α=0.01σ1. The resulting image at tl is computed as(41)Itl=1fxtl22+fytl22+fztl22.


### 4.3. Standard Sparse Processing (SSP)

The selected sparse microwave imaging algorithm is based on linearized scattering equation, as TSVD and, implicitly, LSM. In contrast to TSVD, partial systems of equations are formed for each transmitter. When theith antenna is transmitting, the system of equations reads(42)ei=Lifi,Li=Einct1,r1⋅Einct1,ri…EinctL,r1⋅EinctL,riEinct1,r2⋅Einct1,ri⋮EinctL,r2⋅EinctL,ri⋮⋮⋮Einct1,rM⋅Einct1,ri…EinctL,rM⋅EinctL,ri,ei=Δs1,iΔs2,i⋮ΔsM,i,fi=f1if2i⋮fLi,where Li is the partial system matrix, ei is the corresponding measurement vector, and fi is the unknown vector. Similar to (32), data obtained from P transmitting antennas are processed together. For example, if P = 2, the resulting system of equations becomes(43)eiei+1︸e≈LiLi+1︸Lfi.The minimization function is(44)f^i=minfie−Lfi22+γfi1.where γ is the regularization coefficient. The final image is obtained as the superposition of the partial results f^i.


## 4.1. Truncated Singular Value Decomposition (TSVD)

The TSVD measurement model follows from the linearized scattering equation [12]. The target image is obtained from the linear system of equations(37)e=Lf,L=Einct1,r1⋅Einct1,r1…EinctL,r1⋅EinctL,r1Einct1,r1⋅Einct1,r2⋮EinctL,r1⋅EinctL,r2⋮⋮⋮Einct1,rM⋅Einct1,rM…EinctL,rM⋅EinctL,rM,e=Δs1,1Δs1,2⋮ΔsM,M,f=f1f2⋮fL.HereL is the system matrix, e is the known measurement vector, and f is the unknown vector (target image). The TSVD solution of the system is(38)f=∑i=1P1σiuiH⋅evi,where ui, vi are the singular vectors of matrix L, σi are the corresponding singular values, and P is the truncation index, obtained from the condition −20log10σP=TH, where TH is the adopted threshold.


## 4.2. Linear Sampling Method (LSM)

In this method, a system of equations is solved for each search point and each Cartesian component of the incident field:(39)Lfs=gstl,L=Δs1,1⋯Δs1,М⋮⋱⋮ΔsM,1⋯ΔsM,M,gstl=Einc,stl,r1⋮Einc,stl,rМ,fstl=f1f2⋮fL.HereL is the system matrix, gs is the known vector whose elements are s∈x,y,z components of the incident electric field, produced by the antenna array at tl, and fs is the unknown vector. The solution obtained by means of Tikhonov regularization is(40)fstl=∑i=1Mσiσi2+α2uiH⋅gtlvi,where α is the regularization coefficient, typically selected as α=0.01σ1. The resulting image at tl is computed as(41)Itl=1fxtl22+fytl22+fztl22.


## 4.3. Standard Sparse Processing (SSP)

The selected sparse microwave imaging algorithm is based on linearized scattering equation, as TSVD and, implicitly, LSM. In contrast to TSVD, partial systems of equations are formed for each transmitter. When theith antenna is transmitting, the system of equations reads(42)ei=Lifi,Li=Einct1,r1⋅Einct1,ri…EinctL,r1⋅EinctL,riEinct1,r2⋅Einct1,ri⋮EinctL,r2⋅EinctL,ri⋮⋮⋮Einct1,rM⋅Einct1,ri…EinctL,rM⋅EinctL,ri,ei=Δs1,iΔs2,i⋮ΔsM,i,fi=f1if2i⋮fLi,where Li is the partial system matrix, ei is the corresponding measurement vector, and fi is the unknown vector. Similar to (32), data obtained from P transmitting antennas are processed together. For example, if P = 2, the resulting system of equations becomes(43)eiei+1︸e≈LiLi+1︸Lfi.The minimization function is(44)f^i=minfie−Lfi22+γfi1.where γ is the regularization coefficient. The final image is obtained as the superposition of the partial results f^i.


## 5. Numerical Results

The antenna array consisted of 72 half-wavelength dipoles, divided into six circles uniformly distributed along theθ-direction and operating at 2 GHz. We considered arrays with φ- and θ-oriented dipoles, as illustrated in Figure 3. The array response was computed numerically using the 3D EM solver WIPL-D Pro [40]. The data were corrupted by additive white Gaussian noise (AWGN). The signal-to-noise ratio (SNR) was determined as the ratio of the power of the differential signals to that of noise. We computed the target images in xOy and xOz planes of the Cartesian coordinate system. In each plane, we used a 50×50 search grid.Figure 3
Antenna array configurations. (a)φ-Oriented dipoles and (b) θ-oriented dipoles.


### 5.1. Cross-Shaped PEC Target

As the first example, we considered a cross-shaped PEC object depicted in Figure4. The dimensions of the objects were chosen in such a way that βa≈4, where a is the radius of the minimal sphere encompassing the target. According to (7), for adopted target size, nmax=3. The images were computed for each harmonic in the range (m, n), n≤nmax, m≤n. We noted that for φ-oriented dipoles, the TM11 harmonic was almost sufficient for the target reconstruction, particularly in the xOy plane (Figures 5(a) and 5(b)). For θ-oriented dipoles, the major contribution came from the TM10 harmonic (Figures 5(c)–5(d)). As we will show later, the same harmonics appeared “dominant” in other examples that we studied.Figure 4
3D model of the cross-shaped target.Figure 5
Cross-shaped PEC target images obtained using dominant harmonics. (a)xOy plane and φ-oriented dipoles, (b) xOz plane and φ-oriented dipoles, (c) xOy plane and θ-oriented dipoles, and (d) xOz plane and θ-oriented dipoles. The adopted signal-to-noise ratio was SNR = 10 dB.
(a)(b)(c)(d)In Table1, we show the contributions of all harmonics, described by the reconstruction ability, η2 (the values in each column do not sum up to 1, because in many cases the results obtained from different harmonics overlap). The contributions of harmonics are unique, but some general conclusions can be drawn. For φ-oriented dipoles, the most valuable information came from the harmonics with the index m = 1 (i.e., harmonics: TM11, TE11, TM21, TE11, TE31, and TM31). The harmonics with the index m = 0 (TE10, TE20, and TE30) can be excluded from the analysis with φ-oriented dipoles. These harmonics neglect the φ-dependence, which is critical for this array configuration. Harmonics with borderline values of the indices, m=2,3, tend to produce many artifacts; hence, they can be also eliminated. For θ-oriented dipoles, similar conclusions hold, except that harmonics with m=0 can be used in the image reconstruction, due to predominantly vertical dipole orientation. To simplify the image processing, in the further analysis, we used harmonics shown in red in Table 1.Table 1
Contributions of each harmonic for the cross-shaped PEC target.
φ-dipoles xOy plane<>φ-dipoles xOz plane<>θ-dipoles xOy plane<>θ-dipoles xOz plane<>TE1000.12——TE110.310.130.250TE200.130——TE210.190.250.250.25TE22000.060.13TE300.100——TE310.210.030.250.06TM10——0.630.69TM110.880.720.250.31TM20——0.060.31TM210.190.160.130.06TM220.050.1500.06TM30——0.240.13TM310.050.150.250Figure6 shows the complete target reconstructions, obtained with both arrays. The target shape was accurately estimated, regardless of the plane and dipole orientation. The parameters η1 and η2 were listed in Table 2. In particular, the reconstruction ability (η2) was close to one, whereas the precision η1 was slightly lower, mainly due to the pixels that encircle the target contour from the outside, which actually make the target “thicker.” For the purpose of comparison, we have also included the accuracy parameters calculated for the noiseless case.Figure 6
Cross-shaped PEC target images obtained using multiple harmonics. (a)xOy plane and φ-oriented dipoles. (b) xOz plane and φ-oriented dipoles. (c) xOy plane and θ-oriented dipoles. (d) xOz plane and θ-oriented dipoles. The adopted signal-to-noise ratio was SNR = 10 dB.
(a)(b)(c)(d)Table 2
Parametersη1 and η2 calculated for the cross-shaped PEC target.
PlaneSensor orientationSNRη1η2xOyφ∞0.731.0010 dB0.731.00θ∞0.640.8810 dB0.600.88xOzφ∞0.801.0010 dB0.611.00θ∞0.680.9410 dB0.620.88We note that the dominant harmonics have simple form:(45)EϕTMr,θ,ϕ≈Eϕ,1,1TMr,θ,ϕ=πβ2rβH1/22βr−1rH3/22βrDcosϕ+Esinϕ,(46)EθTMr,θ,ϕ≈Eθ,1,0TMr,θ,ϕ=Dcosθπβ2rβH1/22βr−1rH3/22βr,thus allowing the significant simplification of the reconstruction process at the cost of some accuracy loss. The obtained formula (45) resembles 2D case. However, the more search planes deviate from the xOy plane, the use of higher-order spherical harmonics becomes more important, and those harmonics do not have their counterpart in the 2D case. In addition, higher-order harmonics have their role in obtaining the fuller image (i.e., reducing point-like nature of sparse images).


### 5.2. Cross-Shaped Dielectric Target

In our next example, we investigated the algorithm’s ability to reconstruct the cross-shaped object made of dielectric with the complex permittivityεr=10−j1. Figure 7 shows the reconstruction results obtained using the same harmonics as in previous example, for SNR = 10 dB. The shape of the target was clearly visible in the xOz plane. In the xOy plane, the quality of the reconstruction was somewhat lower. Table 3 displays the values of η1 and η2 parameters, calculated for both noiseless and noisy scenarios. The accuracy parameters were slightly lower than those calculated for the cross-shaped PEC target. Again, for φ-oriented dipoles, the dominant harmonic was TM11. For θ-oriented dipoles, the major contribution came from TM10 harmonic.Figure 7
Cross-shaped dielectric target images obtained using multiple harmonics. (a)xOy plane and φ-oriented dipoles. (b) xOz plane and φ-oriented dipoles. (c) xOy plane and θ-oriented dipoles. (d) xOz plane and θ-oriented dipoles. The adopted signal-to-noise ratio was SNR = 10 dB.
(a)(b)(c)(d)Table 3
Parametersη1 and η2 calculated for the cross-shaped dielectric target.
PlaneSensor orientationSNRη1η2xOyφ∞0.650.6310 dB0.610.63θ∞0.670.7510 dB0.650.59xOzφ∞0.740.8810 dB0.710.84θ∞0.750.8110 dB0.700.78


### 5.3. Two Targets

We studied the reconstruction of two cross-shaped PEC objects, placed in thexOy plane. The dimensions of each target were the same as in previous examples. Due to the increased complexity of the overall scenario, we used a denser antenna array, consisting of 104 dipoles. Figures 8(a) and 8(b) show the obtained results for two target positions using only the dominant harmonic (TM11). As can be observed, both objects were clearly distinguishable. In Figure 8(c), we show the reconstruction results for the spheres, whose radii differ for 50%. Again, the dominant harmonic was also TM11. Interestingly, in this case, the higher-order harmonics detected only the larger sphere.Figure 8
Reconstruction of two objects inxOy plane and φ-oriented dipoles. The adopted signal-to-noise ratio was SNR = 10 dB.
(a)(b)(c)


### 5.4. Algorithm Summary

The performed extensive numerical analysis allowed us to simplify the proposed algorithm whose final implementation consists of the following steps:(i)
Compute the target image using the dominant harmonic TM11 for φ-oriented dipoles or TM10 for θ-oriented dipoles. For the regularization, use γ slightly less than the value corresponding to the knee of the L-curve (i.e., γ = 0.4). Set the size of the transmission array to P = 1 for larger arrays or P = 2 for smaller arrays (results are very similar, but smaller P reduces the computational time).(ii)
Compute the target images using higher-order harmonics (shown in red in Table1). Use γ slightly higher than the value corresponding to the knee of the L-curve (i.e., γ = 0.6).(iii)
Superimpose the partial results.(iv)
Compute the convex hull using the dominant harmonic andγ = 2. Binarize the obtained image to obtain a digital filter. The procedure for selecting threshold was described in [28].(v)
Apply the binarized convex hull as a digital filter to remove artifacts, if necessary.The rationale behind the described regularization parameter selection is the following: images obtained using the dominant harmonic are stable in a wide range ofγ. Hence, by selecting γ that is slightly smaller than the knee-curve value, we increase the portion of the target that is being restored. On the other hand, when the suboptimal harmonics are used, selecting a value of γ slightly higher than the knee-curve value allows suppressing the artifacts.


## 5.1. Cross-Shaped PEC Target

As the first example, we considered a cross-shaped PEC object depicted in Figure4. The dimensions of the objects were chosen in such a way that βa≈4, where a is the radius of the minimal sphere encompassing the target. According to (7), for adopted target size, nmax=3. The images were computed for each harmonic in the range (m, n), n≤nmax, m≤n. We noted that for φ-oriented dipoles, the TM11 harmonic was almost sufficient for the target reconstruction, particularly in the xOy plane (Figures 5(a) and 5(b)). For θ-oriented dipoles, the major contribution came from the TM10 harmonic (Figures 5(c)–5(d)). As we will show later, the same harmonics appeared “dominant” in other examples that we studied.Figure 4
3D model of the cross-shaped target.Figure 5
Cross-shaped PEC target images obtained using dominant harmonics. (a)xOy plane and φ-oriented dipoles, (b) xOz plane and φ-oriented dipoles, (c) xOy plane and θ-oriented dipoles, and (d) xOz plane and θ-oriented dipoles. The adopted signal-to-noise ratio was SNR = 10 dB.
(a)(b)(c)(d)In Table1, we show the contributions of all harmonics, described by the reconstruction ability, η2 (the values in each column do not sum up to 1, because in many cases the results obtained from different harmonics overlap). The contributions of harmonics are unique, but some general conclusions can be drawn. For φ-oriented dipoles, the most valuable information came from the harmonics with the index m = 1 (i.e., harmonics: TM11, TE11, TM21, TE11, TE31, and TM31). The harmonics with the index m = 0 (TE10, TE20, and TE30) can be excluded from the analysis with φ-oriented dipoles. These harmonics neglect the φ-dependence, which is critical for this array configuration. Harmonics with borderline values of the indices, m=2,3, tend to produce many artifacts; hence, they can be also eliminated. For θ-oriented dipoles, similar conclusions hold, except that harmonics with m=0 can be used in the image reconstruction, due to predominantly vertical dipole orientation. To simplify the image processing, in the further analysis, we used harmonics shown in red in Table 1.Table 1
Contributions of each harmonic for the cross-shaped PEC target.
φ-dipoles xOy plane<>φ-dipoles xOz plane<>θ-dipoles xOy plane<>θ-dipoles xOz plane<>TE1000.12——TE110.310.130.250TE200.130——TE210.190.250.250.25TE22000.060.13TE300.100——TE310.210.030.250.06TM10——0.630.69TM110.880.720.250.31TM20——0.060.31TM210.190.160.130.06TM220.050.1500.06TM30——0.240.13TM310.050.150.250Figure6 shows the complete target reconstructions, obtained with both arrays. The target shape was accurately estimated, regardless of the plane and dipole orientation. The parameters η1 and η2 were listed in Table 2. In particular, the reconstruction ability (η2) was close to one, whereas the precision η1 was slightly lower, mainly due to the pixels that encircle the target contour from the outside, which actually make the target “thicker.” For the purpose of comparison, we have also included the accuracy parameters calculated for the noiseless case.Figure 6
Cross-shaped PEC target images obtained using multiple harmonics. (a)xOy plane and φ-oriented dipoles. (b) xOz plane and φ-oriented dipoles. (c) xOy plane and θ-oriented dipoles. (d) xOz plane and θ-oriented dipoles. The adopted signal-to-noise ratio was SNR = 10 dB.
(a)(b)(c)(d)Table 2
Parametersη1 and η2 calculated for the cross-shaped PEC target.
PlaneSensor orientationSNRη1η2xOyφ∞0.731.0010 dB0.731.00θ∞0.640.8810 dB0.600.88xOzφ∞0.801.0010 dB0.611.00θ∞0.680.9410 dB0.620.88We note that the dominant harmonics have simple form:(45)EϕTMr,θ,ϕ≈Eϕ,1,1TMr,θ,ϕ=πβ2rβH1/22βr−1rH3/22βrDcosϕ+Esinϕ,(46)EθTMr,θ,ϕ≈Eθ,1,0TMr,θ,ϕ=Dcosθπβ2rβH1/22βr−1rH3/22βr,thus allowing the significant simplification of the reconstruction process at the cost of some accuracy loss. The obtained formula (45) resembles 2D case. However, the more search planes deviate from the xOy plane, the use of higher-order spherical harmonics becomes more important, and those harmonics do not have their counterpart in the 2D case. In addition, higher-order harmonics have their role in obtaining the fuller image (i.e., reducing point-like nature of sparse images).


## 5.2. Cross-Shaped Dielectric Target

In our next example, we investigated the algorithm’s ability to reconstruct the cross-shaped object made of dielectric with the complex permittivityεr=10−j1. Figure 7 shows the reconstruction results obtained using the same harmonics as in previous example, for SNR = 10 dB. The shape of the target was clearly visible in the xOz plane. In the xOy plane, the quality of the reconstruction was somewhat lower. Table 3 displays the values of η1 and η2 parameters, calculated for both noiseless and noisy scenarios. The accuracy parameters were slightly lower than those calculated for the cross-shaped PEC target. Again, for φ-oriented dipoles, the dominant harmonic was TM11. For θ-oriented dipoles, the major contribution came from TM10 harmonic.Figure 7
Cross-shaped dielectric target images obtained using multiple harmonics. (a)xOy plane and φ-oriented dipoles. (b) xOz plane and φ-oriented dipoles. (c) xOy plane and θ-oriented dipoles. (d) xOz plane and θ-oriented dipoles. The adopted signal-to-noise ratio was SNR = 10 dB.
(a)(b)(c)(d)Table 3
Parametersη1 and η2 calculated for the cross-shaped dielectric target.
PlaneSensor orientationSNRη1η2xOyφ∞0.650.6310 dB0.610.63θ∞0.670.7510 dB0.650.59xOzφ∞0.740.8810 dB0.710.84θ∞0.750.8110 dB0.700.78


## 5.3. Two Targets

We studied the reconstruction of two cross-shaped PEC objects, placed in thexOy plane. The dimensions of each target were the same as in previous examples. Due to the increased complexity of the overall scenario, we used a denser antenna array, consisting of 104 dipoles. Figures 8(a) and 8(b) show the obtained results for two target positions using only the dominant harmonic (TM11). As can be observed, both objects were clearly distinguishable. In Figure 8(c), we show the reconstruction results for the spheres, whose radii differ for 50%. Again, the dominant harmonic was also TM11. Interestingly, in this case, the higher-order harmonics detected only the larger sphere.Figure 8
Reconstruction of two objects inxOy plane and φ-oriented dipoles. The adopted signal-to-noise ratio was SNR = 10 dB.
(a)(b)(c)


## 5.4. Algorithm Summary

The performed extensive numerical analysis allowed us to simplify the proposed algorithm whose final implementation consists of the following steps:(i)
Compute the target image using the dominant harmonic TM11 for φ-oriented dipoles or TM10 for θ-oriented dipoles. For the regularization, use γ slightly less than the value corresponding to the knee of the L-curve (i.e., γ = 0.4). Set the size of the transmission array to P = 1 for larger arrays or P = 2 for smaller arrays (results are very similar, but smaller P reduces the computational time).(ii)
Compute the target images using higher-order harmonics (shown in red in Table1). Use γ slightly higher than the value corresponding to the knee of the L-curve (i.e., γ = 0.6).(iii)
Superimpose the partial results.(iv)
Compute the convex hull using the dominant harmonic andγ = 2. Binarize the obtained image to obtain a digital filter. The procedure for selecting threshold was described in [28].(v)
Apply the binarized convex hull as a digital filter to remove artifacts, if necessary.The rationale behind the described regularization parameter selection is the following: images obtained using the dominant harmonic are stable in a wide range ofγ. Hence, by selecting γ that is slightly smaller than the knee-curve value, we increase the portion of the target that is being restored. On the other hand, when the suboptimal harmonics are used, selecting a value of γ slightly higher than the knee-curve value allows suppressing the artifacts.


## 6. Comparison with Other Methods

When juxtaposing two methods, there are many critical criteria to consider: reconstruction accuracy, resolution, robustness against the noise, minimal array size, computational time, and ease of implementation. As a first criterion, we studied the resolution of the proposed method (SHSP), contrasted to those of TSVD, LSM, and SSP. We utilized the cross-shaped object from Figure3, whose dimensions were scaled so that the radius of the minimum encompassing ball changed from a=0.25λ to a=λ. The images were computed using ϕ-oriented dipoles in the xOz plane. Figure 9 shows the reconstruction results for a=0.3λ obtained by TSVD (TH = – 20 dB), LSM α=0.1σ1, SSP (γ = 0.4), and SHSP with the settings described in Section 4.4. The SHSP was the only method to recover any parts of the vertical arm (for better clarity, the SSP and SHSP results were shown in logarithmic scale). Since nmax≈1, the relevant harmonics were TM11 and TE11. The harmonic TM11 yielded similar results as SSP. The harmonic TE11 was responsible for the pixels at the endpoints of the vertical arm. No digital filtering was applied. The higher-order harmonics produced images with pixels outside the object’s convex envelope.Figure 9
Reconstruction of the cross-shaped objecta = 0.3λ in xOz plane and φ-oriented dipoles. (a) TSVD, (b) LSM, (c) SSP, and (d) SHSP. SNR = 10 dB.
(a)(b)(c)(d)Figure10 shows the reconstruction results for a=0.46λ. Again, the vertical part was still invisible for all other algorithms apart from SHSP. Since nmax≈2, utilized harmonics were TM11, TM21, TE11, and TE21. The harmonics TM21 and TE11 helped in restoring the vertical arm. No digital filtering was applied. Again, the higher-order harmonics produced images with pixels outside the object’s convex envelope.Figure 10
Reconstruction of the cross-shaped objecta = 0.46λ in xOy plane and φ-oriented dipoles. (a) TSVD, (b) LSM, (c) SSP, and (d) SHSP. SNR = 10 dB.
(a)(b)(c)(d)Figure11 shows the reconstruction results for a=0.62λnmax≈3. The cross shape was invisible in the TSVD image (Figure 11(a)); it starts to show in the LSM image (Figure 11(b)), where it is evident in SSP and SHSP images (Figures 11(c) and 11(d)). Using the rules adopted in Section 4.4, we used TM11, TM21, TE11, and TE21 harmonics. However, only the dominant harmonic (TM11) was sufficient. The higher-order harmonics, such as TE31 and TM31, required digital filtering. However, they were not used since they did not contain any new information.Figure 11
Reconstruction of the cross-shaped objecta = 0.6λ in xOy plane and φ-oriented dipoles. (a) TSVD, (b) LSM, (c) SSP, and (d) SHSP. SNR = 10 dB.
(a)(b)(c)(d)Finally, Figure12 shows the reconstruction results for a=0.77λnmax≈4. Apart from TSVD, all algorithms restored the target’s shape. In SHSP, we used the same harmonics as in the previous case (TM11, TM21, TE11, and TE21). The same conclusions hold as for a=0.62λnmax≈3.Figure 12
Reconstruction of the cross-shaped objecta = 0.75λ in xOy plane and φ-oriented dipoles. (a) TSVD, (b) LSM, (c) SSP, and (d) SHSP. SNR = 10 dB.
(a)(b)(c)(d)Further, we compared the influence of the array size and SNR on the reconstruction accuracy. Thus, we decreased the number of antennas toM = 32. Figures 13 and 14 show the images obtained by LSM and SHSP, respectively, for SNR = 10 dB and SNR = 5 dB. The LSM results were severely affected by noise. In contrast, when SHSP was used, the shape of the cross was visible even for SNR = 5 dB.Figure 13
LSM images of the cross-shaped PEC target inxOy plane. (a) θ-Oriented dipoles, SNR = 10 dB; (b) θ-oriented dipoles, SNR = 5 dB; (c) φ-oriented dipoles, SNR = 10 dB; and (d) φ-oriented dipoles, SNR = 5 dB.
(a)(b)(c)(d)Finally, we performed one more comparison of the two sparse processing algorithms. We considered the star-shaped object shown in Figure15(a), illuminated by the array of ϕ-oriented dipoles (M = 62). Figure 15(b) shows the image obtained by SSP. For better clarity, the image was binarized with a 5% threshold. Although the concave nature of the target was visible, less than 30% of the target’s cross section was identified. Generally, restoring only a portion of the targets’ domain is a significant drawback of the sparsity-based algorithms. In contrast, Figure 15(c) shows the image obtained by SHSP, which restored about 60% of the target’s cross section. The accuracy metrics for both methods were given in Table 4. SHSP yielded better accuracy parameters than SSP, even when only the dominant harmonic was used; see Figure 15.Figure 14
Sparse processing images (SHSP) of the cross-shaped PEC target inxOy plane. (a) θ-Oriented dipoles, SNR = 10 dB; (b) θ-oriented dipoles, SNR = 5 dB; (c) φ-oriented dipoles, SNR = 10 dB; and (d) φ-oriented dipoles, SNR = 5 dB.
(a)(b)(c)(d)Table 4
Parametersη1 and η2 calculated for the star-shaped target from Figure 11(a).
Methodη1η2SHSP (TM11)0.920.42SHSP (TM11 et al.)0.910.57SSP [23]0.820.27Figure 15
(a) Star-shaped object. Reconstruction of the star-shaped object inxOy plane, obtained using f-oriented dipoles, SNR = 5 dB: (b) SSP and (c) SHSP.
(a)(b)(c)In general, SHSP yielded more accurate reconstructions than LSM and TSVD when the target was electrically small or when the amount of data was low. Thus, using SHSP reduces the array size and measurement system complexity without sacrificing the resolution. In addition, SHSP was more robust against the noise than LSM and TSVD. However, LSM and TSVD are computationally less expensive and straightforward to implement. Also, the computational time of the sparse processing is longer than that of LSM and TSVD. Obtaining the image with a dominant harmonic for the array with 72 sensors and 2500 grid points takes about 3 min with the SHSP, less than 2 min with TSVD, and about 0.7 min with LSM. When two sparse processing algorithms are compared, the main advantage of the SHSP is the higher resolution. In addition, the SHSP restores larger portions of a target at the cost of more complex implementation.


## 7. Conclusions

This paper describes a 3D MWI qualitative algorithm specifically designed to analyze complex-shaped objects. The algorithm utilizes spherical wave functions to model the scattered field and incorporate them into the sparse processing framework.The algorithm was tested on a vast number of examples. Due to the sophisticated electromagnetic modeling and the utilization of the L1 regularization, the algorithm was able to capture the abrupt changes in the target shape.The extensive study showed that the presented method is suitable for imaging metallic/dielectric targets even in low SNR scenarios. The performance of the algorithm was compared to that of LSM and TSVD. The analysis showed that the sparse processing outperforms the other two algorithms in imaging electrically small targets or low SNR scenarios with few measurements. However, the computational cost of the sparse processing was higher than that of TSVD and LSM.In addition, the performance of the developed algorithm was compared to that of another sparse processing algorithm, which uses the linearized scattering equation as the kernel. The obtained results showed that the utilization of the spherical harmonics significantly increased the reconstructed portion of the target, which is one of the critical points when the sparsity constraint is utilized.The conducted analysis has also shown the existence of dominant harmonics, which have a major influence on the target reconstruction. Hence, utilizing only the dominant harmonics reduces the computational time with slight performance loss [41–49].



---

*Source: 1839842-2022-03-23.xml*