Efficient 3D Image Reconstruction for Near-Field Microwave Imaging Using Dynamic Metasurface Antenna

In this paper, rapid reconstruction of a three-dimensional (3D) image of a scene for near-field microwave imaging using a dynamic metasurface antenna (DMA) is addressed. The imaging system consists of a fixed transmitting DMA and a single moving antenna as the receiver. In the proposed approach, after pre-processing the raw data and converting it into an applicable form for Fourier computations, the converted data is processed by the presented algorithm (all with fast Fourier computations, without the need for the heavy Stolt interpolation step). Furthermore, after extracting a closed-form expression for image reconstruction, a minimization problem is defined using the compressive-sensing (CS) theory, which allows a sparse pattern to be sampled instead of collecting full data, thereby significantly improving the data acquisition rate (according to the CS rate). The performance of the proposed approach is examined by computer simulations and analytical discussions.


I. INTRODUCTION
In recent decades, the widespread applications of microwave imaging in areas such as security screening, medical imaging, nondestructive testing and evaluation, structural health monitoring, and through-wall imaging have led to significant advances in various layers of imaging radars, from the physical layer to the system and signal processing layers [1,2]. In imaging systems, gathering the scene data is usually done either by mechanical scanning or by electronic scanning (large phased arrays with independent antennas) [3,4]. Although the electronic scanning mechanism may be a suitable solution to improve the data acquisition rate, phased array antennas may be expensive and bulky. In addition, due to the need for complex control circuitry and a large number of radio frequency discrete components (such as amplifiers, attenuators, filters, phase shifters, etc.) they usually suffer from high power consumption.
As an alternative to conventional imaging techniques, some modern computational imaging schemes [5] such as coded aperture [6] and single-pixel [7] with relaxed hardware constraints have been used. Recently, dynamic metasurface antenna (DMA) has been proposed as an effective platform for computational microwave imaging [8,9]. DMA contains a waveguide with subwavelength metamaterial elements called unit-cells [10]. The DMA can produce different radiation patterns even at a single frequency [9]. This is because unit-cells can be actively modulated throughout the metasurface [11].
Although employing DMAs simplifies the hardware part, they do not produce uniform radiation patterns due to physical layer compression. As a result, the raw information provided by them cannot be directly processed by efficient computational Fourier transform (FT)-based scene reconstruction techniques. This can limit the advantage of fast Fourier calculations for an imaging system with a DMA transmitter (Tx) or receiver (Rx). To address this issue, in [12][13][14], an effective solution compatible with the range migration algorithm (RMA) [15] is presented that allows the expression of DMA-derived measurements in the spatialfrequency domain and the transformation of raw data on Fourier bases. The system layouts studied in [12][13][14] involve using a DMA panel aperture as Tx and a probe antenna as Rx, which essentially represents a panel-to-probe architecture. Therefore, to reconstruct a three-dimensional (3D) scene, the transceiver system needs to perform a full mechanical scanning in numerous uniform spatial-frequency steps. This is a major challenge for real-time applications in terms of the data acquisition rate as well as for implementation. Also, the image reconstruction process requires a Stolt interpolation step [16], which is the most complex computational part of RMA-based techniques.
To address the above two issues, in this paper, we propose an approach based on fast Fourier and sparse sampling without Stolt interpolation (which we call FFSSWSI). In the proposed approach, by defining a minimization problem based on compressive-sensing (CS), only a limited share of total spatial-frequency data will be required. This allows us to skip the full scanning without having to worry about violating the Nyquist criterion (either in the spatial or frequency domain), and thus significantly increase the data acquisition rate. Also, in the image reconstruction section for computational efficiency, we define a nonuniform inverse fast Fourier transform (NUIFFT)-based process that replaces the heavy Stolt interpolation step. The performance of the FFSSWSI approach will be analyzed by computer simulations.
The rest of this paper is organized as follows: Section II details the proposed approach, including the system and data model, and image reconstruction algorithm; Section III is devoted to presenting the results and discussion; finally, Section IV provides the concluding remarks.
Notation: Throughout the paper, symbols ( ) † . , . p and j stand for the pseudoinverse, p ℓ -norm and the imaginary unit, respectively. Fig. 1 shows an overview of the imaging system studied in this paper. This system uses a fixed DMA along the x-axis (horizontal) as the Tx. This DMA consists of a 1D microstrip waveguide with T N complimentary electric-LC resonators placed at a subwavelength spacing from each other [17]. Each element can be randomly turned on/off by an external stimulus. This can lead to the production of a set of masks (corresponding to diverse and distinct spatial radiation patterns). At a given frequency f multiple measurements can be provided by cycling through M masks (tuning states). As a Rx, a single antenna moves along the y-axis (vertical) to receive the backscattered signal. As a rule, to observe the Nyquist rate and also to benefit from fast Fourier calculations, it is necessary to perform this vertical scanning evenly and with limited intervals. However, this can greatly increase the data acquisition time. We assume here that both spatial and frequency sampling are performed non-uniformly (randomly) with a much smaller number of samples (with a specified CS rate) compared to full sampling. Fig. 1 shows the mechanical scanning points in blue and the omitted points in red. Since the DMA cannot be modeled as a single dipole, the measurement signal g (given below) has a more complex form than the total field that is conventionally calculated only by the Green's function G [12,13] ( ) ( ) ( ) ( ) ; , ; ,

A. SYSTEM AND DATA MODEL
where dV x y z = ∂ ∂ ∂ , R y corresponds to the position of Rx, ρ represents the reflectivity of the target, and r is the position vector to a point in the scene. In the above equation, T U is the radiated field from the aperture, which can be expressed as the superposition of all metamaterial elements (see [12,13] for more details).

B. IMAGE RECONSTRUCTION
Assuming that the signal measured in the aperture plane can be expanded in terms of the fields corresponding to all the masks, we can write [12,13] where m Φ and s represents the field over the aperture for the m -th mask and total field (incident field), respectively. R N denotes the number of vertical scanning points. By creating a set of aperture modes with some degrees of orthogonality, using the delta function screening property and applying some mathematical simplifications, (2) can be written in matrix form as follows: The data is now in the form needed to introduce and apply a This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2022.3186363 fast Fourier-based technique. Therefore, according to the above and the geometry of the system, the received signal can be considered as x y z correspond to the positions of Tx and Rx. Also, T R and R R represent the distances between the positions of Tx and Rx to the target points, respectively. By taking the 2D FT from both sides of (4) on the aperture coordinates, employing the convolution theorem and method of stationary phase [18] and applying some mathematical approximations and simplifications, the signal representation in the wavenumber domain can be written as follows (see [12,19] for details): , .
x y y x x y x y In the above equation, the term acts as a filter in the Fourier domain. Also, from the above equation, it can be concluded that to retrieve the target reflectivity, a 3D inverse Fourier transform must be applied to the signal ( ) , , x y S k k k upon defining an equispaced distribution. Such mapping is known in the literature as Stolt interpolation [16]. Stolt interpolation is computationally the most expensive step of image reconstruction in RMA-based techniques. In the following, we intend to solve this problem by replacing the interpolation step with a less computationally expensive mechanism while improving the data acquisition rate according to the structure introduced in Section II-A. Suppose we denote the output signal obtained by applying x y S k k k ɶ . This filtering is performed by using the appropriate amplitude term mentioned above as well as backpropagating the signal to the center of the scene in the range 0 z (also called motion compensation). Motion compensation corrects the wavefront curvature of all scatterers at the same range as the scene center [20]. In this way, we can rewrite (5) as follows: and by discretizing it in spatial and spectral domains, we have Details of implementing a 1D NUIFFT operator can be found in detail in [21,22]. Therefore, according to the above, the relationship between the measurement data and the reconstructed image in the introduced system can be written in the following closed-form expression:  (10) where f N denotes the total number of frequency samples.
As mentioned in the previous sections, the introduced system does not provide the desired data acquisition rate due to numerous scanning points in the vertical direction and the need for many frequency samples. We address this issue by defining a minimization problem using CS theory and based on the relation extracted in (10) so that there is no need to This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. ), it can be represented by using its sparse transform domain vector θ as a compressible signal [23,24] 1 .
Let us denote the set of randomly selected reduced measurements by , where = A Γψ is the sensing matrix. The recovery problem can be described as the following convex 1 ℓ -minimization To solve (14), we define a cost function as follows [26]: where α is the regularization parameter that controls the relative weight of the two terms. Details of implementation are provided in Section III. By solving the above problem, s can be retrieved with much less data (according to CS theory), and then the scene image can be reconstructed by using (10). A summary of the main steps of implementing the proposed approach with the variables introduced above is given in the form of a block diagram in Fig. 2.

III. SIMULATION RESULTS AND DISCUSSION
All computations are performed on a MATLAB R2020b platform running on a 64-bit Windows 10 operating system with 12GB of random-access memory and a Core-i7 central processing unit at 2.7GHz. To implement NUIFFT, a fast Gaussian gridding-based method is used, which includes the main steps of initialization, convolution, IFFT and deconvolution [27]. Γ is designed so that only one 1 is randomly (with a uniform distribution) placed in each row, and the rest of the entries are 0. For the dictionary matrix, we considered the discrete Fourier transform domain, which provided the best performance compared to the discrete cosine transform and wavelet domains in our experiments. To minimize the cost function in (15), we used the two-step iterative shrinkage/thresholding (TwIST) algorithm [28], which converges much faster than similar algorithms such as IST [29].
The simulation parameters are given in Table I, where Q is the quality factor of the metamaterial elements [30] and λ is the wavelength corresponding to the highest frequency in free space. In each mask, half of the elements are randomly turned on. In all experiments, a 3D B-shaped distributed target (see Fig. 1) is considered. Due to the size of the aperture ( 2 0.71 0.71m × ) and the operating frequency, targets with a range of less than approximately 66m are located in the near-field (NF) region [31]. First, we compare the performance of the image reconstruction algorithm presented in this paper with the adapted RMA in [12,13] without sparse data (with full data). This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2022.3186363 Note that since the geometry of the imaging system in this work is somewhat different from the work [12,13], we adapted the algorithm [12,13] based on the geometry of Fig.  1 to make a fair comparison. Figs. 3(a) and 3(b) show the outputs of the algorithm [12,13] as isosurfaces of the normalized image when using the spline interpolation [32] and an interpolation method with a simpler implementation based on the minimum distance criterion (given below), respectively: x q y q z q x q y q z q x q y q l z q l x q l y q q S k k k where 1,..., z q N ′′ = and Ŝ is initialized to zero. The reconstructed images are displayed in 3D and 2D (xy-plane) views. Note that the images are distance color-coded and the colorbar represents the range in meters.
It can be seen that Fig. 3(b) has fewer sidelobes than Fig.  3(a), however, it requires more computational time. The average computational times for the interpolation and 3D IFFT steps for Figs. 3(a) and 3(b) are calculated as 1.03 and 0.86 seconds, respectively. Fig. 3(c) shows the corresponding outputs obtained by the FFSSWSI approach. The target image is well reconstructed and completely recognizable.
Although the reconstructed image is less focused along the zaxis than in Figs. 3(a) and 3(b), the average computational time of alternative operations (1D NUIFFT + 2D IFFT instead of Stolt interpolation + 3D IFFT) is 0.63 seconds, which shows the computational superiority of the proposed approach.
To validate this observation analytically, we extracted the major computational complexity of the above operations.  Fig. 4 shows a comparison of the computational complexity between the two schemes versus different numbers of frequency and spatial samples. As can be seen, the scheme employed in the proposed approach shows much better computational performance than the conventional interpolation-based scheme, especially when the number of samples increases. In addition, in Table II, we calculated the normalized mean squared error (NMSE) values [4,24] and This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. summarized them with the computational times mentioned above. To calculate the NMSE, Fig. 3(a) is taken as the reference image. The NMSE value calculated for Fig. 3(c) is very close to that for Fig. 3(b), which can be verified by visual inspection of the figures. Furthermore, Table II provides the image contrast and entropy values [33][34][35] averaged from 2D reconstructed images (on the xy-plane) focused in different ranges. These results agree with the above findings in terms of the performance of various techniques.  In addition to the experiments performed on the full data, we also examined the performance of the FFSSWSI approach on the spatial-spectral sparse data and compared the results with the adapted RMA method. , respectively. These values are equivalent to CS rates of 25% and 42.25%, respectively. As expected, due to non-compliance with the Nyquist criterion in both spatial and frequency sampling, the resulting images are aliased and the target is not recognizable. Also, as can be seen, the deviation from the target position is much more severe in the y-and z-axes, because the sparse data are related to mechanical scanning samples in the y-direction and frequency samples (which contain the range information), but in the horizontal direction, the full data (DMA related data) is used. The corresponding images, when the proposed approach is used, are shown in Figs. 5(c) and 5(d). Although these images contain some undesirable sidelobes, unlike Figs. 5(a) and 5(b), the target image is well recognizable. Also, by comparing Figs. 5(a) and 5(b), it is clear that increasing the CS rate can be effective in improving the quality of the reconstructed image; however, this comes at the cost of increasing the number of measurements (reducing the data acquisition rate). Table III lists the NMSE, contrast and entropy values, which are consistent with the findings in Fig. 5 and the description above. As the CS rate increases, the error value decreases, the contrast improves, and the entropy value, which is interpreted here as the degree of irregularity (heterogeneity), decreases [33][34][35].

IV. CONCLUSION
An effective approach to 3D image reconstruction in a NF microwave imaging system including a DMA was addressed.
The most important achievement of this work is to formulate a solution to match the sparse data (instead of the full measurements) and to eliminate the Stolt interpolation step by replacing it with a more efficient computational process (fast Fourier type). The results of simulations and computational analyses showed that using the proposed approach, a 3D image of the scene can be reconstructed faster than conventional RMA-based methods. Adapting the solution with sparse data can lead to an improved data acquisition rate (proportional to CS rate) in the introduced system.   FIG. 5 (FIG. 3(A)