Fourier Compatible Near-Field Multiple-Input Multiple-Output Terahertz Imaging with Sparse Non-Uniform Apertures

Among image reconstruction methods, Fourier transform-based techniques provide computationally better performance. However, conventional Fourier-based reconstruction techniques require uniform data sampling at the radar aperture. In this paper, a multiple-input multiple-output (MIMO) scenario for near-field (NF) terahertz imaging systems is considered. A compressive-sensing-based method compatible with efficient fast Fourier-based techniques for image reconstruction is proposed. To reduce the error due to the multistatic array topology in the NF, a multistatic-to-monostatic conversion is used. Employing the proposed method significantly reduces the number of antennas and channels. This, in addition to saving hardware resources, can improve the overall performance of the system depending on the type of channel access scheme. The results based on both numerical and electromagnetic data, presented as reconstructed images of the scene, confirm the performance of the proposed method.


I. INTRODUCTION
In recent years, active terahertz (THz) imaging ranging from 0.1 to 10 THz has received increasing attention in the fields of security screening, aerial imaging, medical diagnostics and non-destructive testing [1,2]. THz imaging modalities offer non-ionizing radiation and can operate in all-weather conditions [3]. Although visible and infrared frequencies can provide very high resolution, they cannot penetrate through some materials such as clothing. Therefore, for applications such as security screening, THz frequencies are ideal and can be used to detect hidden objects under clothing. These waves provide a resolution commensurate with the size of the aperture [4].
Image reconstruction is a mathematical process that generates scene images from projection data acquired at the target position. The type of mechanism used in the image reconstruction process has a fundamental effect on image quality. In literature, there are many techniques for image reconstruction [5][6][7][8]. Among those, Fourier-based reconstruction techniques offer significant potential due to high computational efficiencies. However, such techniques suffer from several limitations, such as uniform sampling requirement, typically achieved at the Nyquist limit.
Compressive-sensing (CS) is a framework for reconstructing a signal from a reduced number of measurements [9]. According to CS theory, sparse or compressible data can be reconstructed from a small set of measurements fewer than usual [9]. Therefore, using this method makes a significant improvement in applications where the amount of data is large, or receiving and storing data is a time-consuming and costly process, or the time to receive all the data is limited and short. The idea of CS has already been used in imaging radar systems [10][11][12][13]. In [10], a three-dimensional (3D) synthetic aperture radar imaging system is introduced to reduce the number of measurement points in both frequency and space domains. However, the system introduced in [10] considers a monostatic scanning topology which, unlike the multistatic topology, is not challenging in the near-field (NF). In [11], the problem of 2 VOLUME XX, 2021 receivers position optimization by CS technique in a singleinput multiple-output scenario is considered. In [13], a superresolution technique is employed to reconstruct the image of a target which is acquired by using CS. This technique is based on scanning the scene by using randomly patterned masks with fixed pixel sizes. However, these methods are not compatible with efficient Fourier-based techniques. In [12], a CS approach using a planar scanning setup has been presented to improve the dynamic range compared with the conventional imaging technique based on time reversal of the measured fields. For 2D scanning in a multistatic scenario, a fixed 2D array, or a 1D array with 1D mechanical scanning can be used. However, to satisfy the Nyquist criterion, such a setup requires many antennas or spatial sampling points. For example, to mechanically scan a scene with a width of 0.3m at 220GHz and an inter-element spacing of 2  , 441 elements will be required. Moreover, if we have a fixed densely-populated 2D aperture, then we would need 441 441  elements. This may be a major challenge for THz imaging.
In this paper, a multiple-input multiple-output (MIMO) scenario for NF THz imaging systems is considered. Unlike many existing works focusing on target domain sparsity [14,15], our focus is not on the target domain sparsity but the physical layer sparsity. The CS technique is employed here to significantly reduce the number of physical antennas and spatial sampling points. In fact, many antennas (or channels) or sampling points that are evenly spaced may be eliminated. In addition to saving hardware resources, this can reduce data acquisition time in time-division-based methods [16,17], reduce bandwidth and sampling frequency in frequencydivision-based methods [18], and simplify implementation in coding-based methods [19]. To correct the error due to the multistatic array topology in the NF and consequently to improve the estimation accuracy, a multistatic-to-monostatic conversion is used. The application of this conversion technique to a sparse aperture at microwave frequencies has also been studied in [20]. Although the aperture in [20] is sparse, the transmitting (Tx) and receiving (Rx) antennas are uniformly sampled along the edges of each cluster, producing a uniformly sampled effective aperture pattern. In this work, we propose a Fourier-based processing technique that can work with non-uniformly distributed sparse aperture layouts. This gives us the freedom to choose an arbitrary set of Tx and Rx pairs (not necessarily uniformly sampled) and reduce the number of spatial sampling points (sparse sampling) to facilitate compressive imaging and relax the Nyquist sampling criterion.
The main contributions and novelties of this paper are summarized below: • Mathematical development of a Fouriercompatible method for NF MIMO THz imaging with sparse non-uniform apertures: Although Fourier transform (FT)-based methods are commonly known as efficient image reconstruction algorithms, they conventionally require uniform sampling. Sampling at the Nyquist rate ensures that raw data can sample a complete set of Fourier components, such that a DFT can be defined over the set of data. The two features in this paper's scenario that disrupt the uniform sampling structure are NF multistatic imaging and non-uniform apertures. Since the effective phase center principle is theoretically valid only under the far-field (FF) assumption, a phase compensation mechanism is required to adopt the Fourier-based image reconstruction technique for large MIMO apertures in the NF. We achieve such an adaptation using a multistatic-to-monostatic conversion. In this context, for the first time, the feasibility of using multistatic-to-monostatic conversion for THz NF imaging is demonstrated.
• In achieving the above, the ability to work with non-uniformly distributed apertures is of significant importance. Non-uniform apertures can be due to sparse spatial sampling or sparse arrangement of array elements or due to the inherent sparse configuration of the array. In this regard, we retrieve the data required for the image reconstruction process by solving a CS problem and forming a dictionary matrix in the Fourier domain. Also, for inherently non-uniform configurations, an interpolationbased solution is considered as a preprocessing step to improve the quality of the reconstructed image. As a result, the compatibility of the proposed method with non-uniform THz multistatic imaging in the NF is demonstrated for the first time.
• Validation of the performance of the proposed method is realized in various experiments with both numerical and electromagnetic data.
The rest of this paper is organized as follows: In Section II, the system model is presented; Section III briefly describes the main concepts and equations of the CS technique; Section IV presents the proposed CS-based method for NF imaging; In Section V, we present the simulation results; Section VI presents the concluding remarks.
Notation: Throughout the paper, superscript * represents the complex conjugate. Symbols . , . p and .ˆ stand for the absolute value, p -norm and estimation, respectively.

II. SYSTEM MODEL
Consider the setups shown in Fig. 1. We assume that the radar measurements are obtained by a MIMO structure, either by a setup consisting of antenna elements filling a 2D aperture (Setup 1 in Fig. 1(a)) or by an array of 1D antenna elements that are mechanically scanned (Setup 2 in Fig.  1  According to the effective phase center principle, under the FF assumption, a multistatic array topology with TR NN + physical elements can be considered as a monostatic virtual array with TR NN elements (equal to the number of Tx-Rx channels) [23] (see Fig. 1(a)). Also, the Setup 2 configuration is equivalent to a denser linear virtual array of length  Fig. 1(b)). However, for NF imaging, we need a more accurate model to reconstruct the image. Details of this adaptation are given in Section IV.
The measured backscatter signal can be written as follows [24]: , , , , , , , is the wavenumber of the corresponding to the frequency f ,  is target reflectivity and c is the speed of light. Based on this, a raw 2D data of size captured over the xy-domain can be constructed for Setups 1 and 2, respectively, where s N is the number of frequency points. We refer to this raw data as S .

III. COMPRESSIVE-SENSING
According to CS theory, sparse or compressible signals can be reconstructed from a reduced number of measurements. Assume that signal x exhibits sparsity in certain orthonormal basis ψ . In matrix form, the signal x can be represented by using its sparse transform domain vector (vector of coefficients) θ as a K -sparse signal [25]: 11 .
In the CS scenario, a reduced set of measurements According to (2), the measurement procedure can be modeled by projections of x onto vectors   12 , ,..., where = A Θψ is the sensing matrix. Remark 1: In some computational imaging works [26,27], M refers to the total number of measurements and N refers to the number of pixels in the scene (unknown to be reconstructed, e.g., the reflectivity distribution). However, in this paper, both M and N refer to the number of antenna elements. In this work, the compression comes from the fact that we use a fewer number of spatial sampling points (and hence data acquisition channels) on the aperture in comparison to the regularly sampled (at Nyquist limit) case. In computational imaging works [26,27], it is possible to have a single channel (only one data acquisition channel) but the number of measurement modes, M , can still be large. It can even be made MN  by simply increasing the number of frequency points (each frequency point in this concept is a separate measurement).

VOLUME XX, 2021
Under certain conditions (restricted isometry property) [25], the recovery problem can be described as the following minimization: y Aθ (4) However, the solution of (4) requires an exhaustive search and is NP-hard [25]. Hence, the 0 -minimization is replaced by convex 1 -minimization as follows: In the situation when the measurements are corrupted by the noise ( + = + y=Θx n Aθ n ), where n denotes the additive noise, the reconstruction problem can be defined as

IV. PROPOSED METHOD
Normally, to avoid aliasing in image reconstruction, the Nyquist criterion in spatial sampling must be met. This means that an inter-element spacing of x d and y d in the virtual array, along the x-and y-axes, respectively, must be implemented as follows [28]: where  is the wavelength. However, according to the description in the previous section and using the CS technique, it is no longer necessary to obtain samples based on the Nyquist rate and uniform sampling in the x-and y-axes. In fact, with sparse spatial sampling, the data required for the image reconstruction process can be retrieved based on (6). A summary of the main steps for implementing the proposed method is given in the block diagram of Fig. 2, which is described in more detail below. For each i , we formulate the uniformly sampled data ( ) Next, we construct an NN  dictionary matrix ψ (in practice, discrete cosine transform (DCT), DFT, or other dictionaries may be used [10]). Note that the issue of selecting the appropriate basis before retrieving information using CS has been discussed in the literature under the heading of best basis selection and dictionary learning [29,30]. Then we choose a random measurement matrix  . We refer to compressed data as S . To develop a CS approach for the MIMO NF imaging system, we need a more accurate system model for image reconstruction than that provided in Section II. This is because we are dealing with a short-range NF multistatic configuration, while the virtual arrays shown in Fig. 1 are theoretically valid for the FF. Therefore, a phase compensation mechanism is required to adopt the Fourierbased image reconstruction technique for large MIMO apertures in the NF. Suppose x y x y k correspond to the monostatic and multistatic reference signals, respectively [32]. Note that we still have a multistatic topology, however, this conversion provides a virtual (mathematical) approximation of the corresponding monostatic topology. We refer to converted data as S . By using (6), x (an estimation of uniformly sampled data) is retrieved. We refer to compressed data as Ŝ . Thus, although we have not performed spatial sampling evenly, the resultant data can be processed using Fourier-based image reconstruction techniques.
As an instance, 2D target reflectivity can be reconstructed as [33] ( ) A 3D target reflectivity can also be achieved only by modifying the image reconstruction step and using 3D IFFT [32].

V. SIMULATION RESULTS
In this section, the performance results of the proposed method based on numerical and electromagnetic data simulated in MATLAB and FEKO, respectively, are presented. All computations are performed in MATLAB. The simulation parameters for different configurations are given in Tables I-III. All results are provided for signal-tonoise ratio (SNR) of 20dB to ensure a realistic channel response in the presence of additional loss factors. Our previous works in NF imaging [34][35][36] justify this SNR selection. According to the parameters of Configs. 1-5 given in Tables I-III, targets with distances less than approximately 5m, 2m, 5m, 54m and 54m from the imaging system, respectively, are located in the NF [37]. Θ is designed so that only one 1 is randomly placed in each row (without repeating the position in other rows) with a uniform distribution, and the rest of the entries are 0. For the dictionary matrix, we considered the DCT domain, which provided the best performance in our experiments. For data retrieval, wherever it is not mentioned, the smoothed L0 (SL0) algorithm [38] is considered, which has less complexity than other similar algorithms such as robust SL0 (RSL0) [39], orthogonal matching pursuit (OMP) [40] and spectral projected gradient for L1 minimization (SPGL1) [41]. This section provides qualitative comparisons as well as quantitative analyzes such as resolution and normalized mean squared error (NMSE) to examine performance. The computation formula for NMSE is as follows [42]:  in the virtual array. The simulated system has a theoretical cross-range resolution of 5mm (1.67 ) in both xand y-axis [28]. Fig. 4(a) shows a sparse MIMO array scenario in which the physical Tx and Rx antennas in Fig.  3(a) are randomly selected. Fig. 4(b) shows the corresponding virtual array. The rate of compressive sensing ( r M N = ) is approximately 12%. Note that this rate is defined by the size of the data, not just the number of physical antennas. Since in this scenario the Nyquist criterion is not followed and lacks a uniform pattern required for the conventional FFT-IFFT technique, the image reconstructed without using the proposed method is aliased and the object is not recognizable. This can be seen in Fig. 4(c). However, as Fig. 4(d) shows, by using the proposed method, we were able to reconstruct the image in such a scenario (with only 35% of the total number of antennas (12% of all channels)), so that the object is identifiable. In Fig. 5(b), we have increased the rate r to 25%. As Fig. 5(c) shows, the conventional FFT-IFFT technique still failed to properly reconstruct the image. As expected, using the proposed method with 0.25 r = , the image is reconstructed with a better resolution than the 12% rate. This advantage comes at the cost of increasing the complexity in the recovery step because the computational complexity of the SL0 algorithm is of order 2 () M . Note that the SPGL1 algorithm is able to estimate the noise power in the optimization problem-solving process. Therefore, it has the advantage that it provides higher accuracy in CS problems in the presence of noise. However, this algorithm has a very high computational complexity of order ( log ) NN . Also, OMP and RSL0 algorithms have computational complexities () KMN and 2 () M , respectively. Note that by employing the proposed method, the overall transmitting time (assuming the use of conventional time-division and code-division techniques) in Figs. 4(a) and 5(a) are reduced to 12% and 25% of the full array in the structure of Fig. 3(a), respectively [24]. Also, by employing the proposed method (assuming the use of the frequency-division technique), the bandwidth (and consequently the sampling rate) are reduced to 35% and 50%, respectively [24]. For further investigation, we increased the operating frequency to 220GHz. Other parameters are similar to Config. 1 (see Table I). The results of simulations based on Config. 2 are shown in Figs. 6-8. The most important limitation of Config. 2 compared to Config. 1 is a significant reduction in aperture size. This is due to the significant increase in frequency, while the inter-element spacing remained constant in terms of wavelength. This difference can be found by comparing Figs. 3(a) and 3(b) with Figs. 6(a) and 6(b), respectively. Such a system has a theoretical cross-range resolution of 3.67 in both x-and y-axis. As expected, even with full array data, the reconstructed image (Fig. 6(c)) does not have the good quality of the image obtained in Fig. 3(c) at 110GHz. The significant decrease in aperture size has also had a significant effect on the quality of images reconstructed with sparse data (see Fig. 7(e) and compare it with Fig. 5(d)). Table IV shows the NMSE values of some images reconstructed by the proposed method. In addition to the SL0 algorithm, we also used the OMP and RSL0 algorithms to retrieve the data, the results of which are shown in Figs. 7(d) and 7(f), respectively. Both qualitative results, i.e. reconstructed images, and quantitative results, i.e. NMSE values (see Table IV), indicate that the OMP algorithm provided poorer performance. Also, the RSL0 algorithm performed slightly better than the SL0. This is because RSL0 is more robust against noise than SL0. We also increased the compressive sensing rate from 25% to almost 50%. However, due to the aperture size limitation, the reconstructed images (Figs. 8(d), 8(e) and 8 (f)) still do not look desirable, despite the relative reduction of the sidelobes.  Note that despite the insufficient aperture size, the results obtained from the proposed method still provide a much better idea of the target than the images obtained from the conventional FFT-IFFT technique (Figs. 7(c) and 8(c)).
To improve the results at the latter frequency, it is necessary to increase the aperture size. For this purpose, we considered Config. 3 (see Table I). The corresponding results are shown in Figs. 9 and 10. Comparing Figs. 9(c) and 10 with Figs. 6(c) and 7, respectively, shows that by increasing the aperture size, we were able to significantly improve the quality of the reconstructed images. Such a system provides a theoretical cross-range resolution of 2.62 in both x-and yaxis. It is emphasized that the unsatisfactory quality of the reconstructed images in Figs. 7 and 8 is more due to the insufficient aperture size in the Config. 2 than to the performance of the proposed method. This can also be examined by comparing NMSEs. For example, the NMSE values in Figs. 7(e) and 10(e) are calculated as 0.0199 and 0.0258, respectively. This means that the similarity of Fig.  7(e) to its reference image (i.e. Fig. 6(c)) is even greater than the similarity of Fig. 10(e) to its reference image (i.e. Fig.  9(c)). Therefore, satisfactory image reconstruction with sparse data requires that the minimum conditions be met for ideal image reconstruction with full data; because according to information theory, no process can increase the information content of sparse data beyond the information content of the original data.
The rest of the results of this section are related to SPA, the structure of which was presented in Section II. Fig. 11 shows the results corresponding to the case in which Config. 4 is used (see Table II). In the case of Fig. 11(d), it is assumed that the target is located at a distance of 0.3m from the radar. The SPA, despite the advantages it provides (including large inter-element spacing [21,22,43]), creates a gap in the center of the virtual array, regardless of whether x N is even or odd (see Figs. 1(b), 11(b) and 11(c)). Although this gap does not pose a problem for non-FT-based image reconstruction techniques (such as generalized synthetic aperture focusing technique (GSAFT) [22]), it can affect the quality of results for Fourier-based techniques that require uniform spatial measurements. In [44], more details are provided and an interpolation-based solution to improve the results is presented. Fig. 11(e) shows the reconstructed image using the improved technique. As can be seen, Fig.  11(e) has less distortion than Fig. 11(d). We then increased the target distance to 1.1m (note that the target is still in the NF). It can be seen that with increasing range, the quality of the reconstructed image has improved. In fact, although increasing the range causes a tolerable degradation in resolution, by bringing the target closer to the FF region, the accuracy of the approximations used in multistatic to monostatic conversion and interpolation steps is improved [19,44,45]. Also, since the illumination footprint of the radar increases with the distance, the reconstructed image in Fig. 11(f) has larger dimensions than the one in Fig. 11(e).

2; (a) sparse MIMO array (50% of all antennas), (b) virtual array, (c) image reconstructed by conventional FFT-IFFT technique, (d) image reconstructed by the proposed method (with the OMP algorithm), (e) image reconstructed by the proposed method (with the SL0 algorithm), (f) image reconstructed by the proposed method (with the RSL0 algorithm).
Author Name: Preparation of Papers for IEEE Access (November 2021) VOLUME XX, 2021 9 The SPA itself has a physically sparse form and due to its properties [21,22,43], compared to the URA, with a smaller number of physical elements, a larger aperture size can be achieved. So here for Setup 2, we consider sparsity in spatial sampling. However, at the end of this section, we will also consider another scenario. Figs. 12 and 13 show the results obtained from the conventional FFT-IFFT technique, the non-uniform FFT (NUFFT) technique [46] and the proposed method for spatial sampling rates of 60% and 70%, respectively. As can be seen, the conventional FFT technique still failed to reconstruct the image properly (see Figs. 12 12(d) and 13(d)), the target can be identified (with more clarity in Fig. 13(d), as expected). Sparse reconstruction analyzes [47,48] show that approximately at least 15% of the signal samples are required for reconstruction by sparse recovery algorithms. Here we perform a numerical study to determine an approximate limit for the CS rate that provides a reconstructed image of satisfactory quality using the proposed method. For this purpose, we consider Configs. 1, 3 and 4. As we saw in the discussions above, the images reconstructed by the full data in these configurations were of good quality. Fig. 14 shows the average values of NMSEs obtained by the proposed method (using OMP and SL0 algorithms) versus r in 100 independent experiments for various configurations. First, let   Fig. 15(a) compared to Fig. 15(b) is consistent with the findings of Fig. 14 To validate the performance of the proposed method, in addition to the numerical data simulated in MATLAB, we used the electromagnetic data simulated in FEKO. The simulation parameters are given in Table III. Fig. 16 shows images reconstructed by using full data. Fig. 17 shows images reconstructed by the conventional FFT-IFFT method, the matched filtering method [49], the NUFFT technique, the GSAFT method [22], and the proposed method in a sparse scenario. It is observed that we can identify the target image only by using the proposed method because images reconstructed by other methods suffer from severe distortion.    This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Note that in the proposed method, in the last step, we use (9) to reconstruct the image from Ŝ (see Fig. 2). Here, instead of using (9) in the proposed mechanism, we reconstruct the images from Ŝ by using matched filtering technique and GSAFT. The related results are shown in Fig.  18. As can be seen, employing the matched filtering technique (which is based on FT) in the proposed mechanism provides a more acceptable image of the scene than employing GSAFT, which is not based on FT. Figs. 18(a) and 18(b) may be compared to Figs. 17(i) and 17(j), respectively. It is observed that the images reconstructed by (9) have a better quality than those reconstructed by employing the matched filtering technique. The reason for this is that the matched filtering technique, despite its more straightforward implementation, has limitations in terms of resolution [50].
For a case like Setup 2, where a 1D array is combined with a mechanical scanning, it is ideal for the sparse scenario to skip some scanning steps (lines) randomly (as in experiments we did the above for Setup 2). This reduces data acquisition time and completes mechanical scanning faster. However, here we consider another scenario and that instead of skipping some scanning steps, we randomly turn off some SPA elements in each step (see Fig. 19). As can be seen, compared to Figs. 17(i) and 17(j), a slight quality improvement is observed in the reconstructed images. Note that this comparison is made with similar compressive sensing rates. The reason for this relative improvement is that in the latter scenario, the data is sparser, which is more desirable for CS.

VI. CONCLUSION
In this paper, a CS-based method compatible with Fourierbased techniques for NF mm-wave imaging was presented in a practical MIMO scenario. To reduce the error due to the multistatic array topology in the NF, we used a multistatic-tomonostatic conversion. By using both numerical and electromagnetic data, we showed in various experiments that we were able to successfully reconstruct the scene images, despite the non-observance of the Nyquist condition (which is a prerequisite for most Fourier-based methods). Employing the proposed method can reduce the complexity of data acquisition in imaging applications. This advantage is particularly important at THz frequencies where conventional Nyquist sampling of the aperture can require an excessive amount of sampling points. Moreover, the capability of the developed Fourier-based reconstruction technique to work with non-uniformly sampled data introduces an additional degree of freedom in designing sparse apertures for imaging at THz frequencies.