{'fraction_non_alphanumeric': 0.03825547513864085, 'fraction_numerical': 0.013159131497321177, 'mean_word_length': 4.539963551158553, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Compressive Sensing Imaging of 3-D Object by a Holographic Algorithm

Existing three-dimensional (3-D) compressive sensing-based millimeter-wave (MMW) imaging methods require a large-scale storage of the sensing matrix and immense computations owing to the high dimension matrix-vector model employed in the optimization. To overcome this shortcoming, we propose an efficient compressive sensing (CS) method based on a holographic algorithm for near-field 3-D MMW imaging. An interpolation-free holographic imaging algorithm is developed and used as a sensing operator, in lieu of the nominal sensing matrix typically used in the CS iterative optimization procedure. In so doing, the problem induced by the large-scale sensing matrix is avoided. With no interpolations required, both the computational speed and the image quality can be improved. Simulation and experimental results are provided to demonstrate the performance of the proposed method in comparison with those of the Omega-K based CS and the traditional Fourier-based imaging techniques.


I. INTRODUCTION
ILLIMETER-WAVE (MMW) has attractive characteristics compared with waves in the microwave band or lower frequency bands. These include higher carrier frequency and wider usable frequency band that enables higher target cross-range and down-range resolutions. Another important feature of MMW is the design of small and light systems and equipment. Accordingly, it is beneficial to adopt MMW for short-range broadband communications [1], [2], high-resolution sensing [3], [4], and radio astronomy [5]. MMW imaging techniques have been widely developed and applied to non-destructive testing [6], biomedical diagnosis [7], and personnel security inspection [8]- [11].
MMWs are capable of penetrating regular clothing to form an image of a person and concealed objects. Most importantly, MMWs are non-ionizing and, therefore, pose no known health Manuscript  hazard at moderate power levels. MMW imaging systems can be classified into two types: passive imaging systems and active imaging systems. Both types have their own offerings and challenges. This paper focuses on active MMW imaging techniques, which typically implement large scale antenna arrays to illuminate the whole human body, leading to high system cost. Compressive sensing (CS) and sparse reconstruction techniques, on the other hand, have been used to reduce the number of array elements, thereby the system cost, without degradation of image quality [12]- [14]. CS has been extensively studied in radar imaging [15]- [19], microwave imaging [20]- [23], array synthesis and diagnosis [24], [25], and direction-of-arrival estimation [26]- [28].
CS methods are typically based on the matrix-vector model, which employs a large-scale sensing matrix in 3-D MMW imaging for personnel inspection and security applications. Such model was incorporated in [29] for a 3-D compressive phased array imaging. CS was applied to single-frequency 2-D MMW holographic imaging in [30] and [31], where a Fourier-based imaging operator represented the sensing matrix. This replacement, in essence, avoided storing and processing of a large-scale sensing matrix which, in turn, simplified imaging and permitted its realization on an ordinary personal computer.
In [32], the 3-D K  algorithm, referred to as range migration algorithm, was used in combination with the CS principle for image reconstruction. Nevertheless, this algorithm includes the forward and inverse Stolt interpolations that entail require high computations and can lead to reduced image fidelity. An interpolation-free SAR imaging algorithm, named range stacking, was proposed in [33] and extended to 3-D imaging in [34]. The range stacking reconstruction method forms the target image at different range points by matched filtering the SAR signal in the spatial frequency domain. The result is integrated over frequencies to yield the marginal Fourier transform of the target function.
In [35] and [36], we considered a single-frequency based auto-focus holographic imaging algorithm. The auto-focus was obtained by calculating the amplitude integral values of the images reconstructed at different focusing range bins. These values draw to a minimum when the image is well focused. In this paper, we extend the above algorithm to the wideband signal, and change the integral variables of the imaging model to eliminate the Stolt interpolation. This changing is similar to the work in [33], [34]. However, unlike these references, we first apply the inverse Fourier transform over the azimuth-and elevation-frequencies with respect to the matched filtered data.

M
The results are then integrated over the fast-time frequencies, which can be considered as a coherent summation of the single-frequency holographic imaging results. These steps amount to a linear relationship between the scene and the measurements, defined by a sensing operator, and as such, the CS problem can be readily formulated. Due to the fact that interpolations are avoided in the iterations underlying the optimization algorithm, the computational speed and optimization solution can both be improved in comparison with the 3-D K  -based CS method [32].
In order to deal with compressed data, which corresponds to selecting few antennas, we utilize a uniform-random spatial undersampling scheme, in lieu of the totally random undersampling scheme. The benefits of the former over the latter sampling strategy can be made evident by analyzing the mutual coherence measure of the sensing operator. In imaging, this measure, in essence, represents the maximum sidelobe value of the point spread function (PSF) in Fourier-based imaging. Accordingly, it can be used to quantify the effect of different undersampling schemes on performance. Also, we compare the PSFs of the proposed sensing operator and that of the K  -based sensing operator, and show the superior performance of the former. The remainder of this paper is organized as follows. Section II presents the formulation of the 3-D interpolation-free holographic imaging algorithm, which is used to construct the operator considered as the sensing matrix of CS. In Section III, we provide the CS imaging method and present the uniform-random under-sampling scheme. The relationship between the mutual coherence and the point spread function is analyzed. Numerical simulations and experimental results are shown in Section IV. Finally, concluding remarks are presented in Section V.

A. Formulation of the interpolation-free holographic imaging algorithm
The configuration of the imaging system is shown in Fig. 1. For a frequency modulated continuous wave (FMCW) transceiver, the transmitted signal can be expressed as, where 0 f is the center frequency, t is the fast-time variable within one pulse repetition interval, and K is the frequency slope of the transmitted signal. The backscattered signal from a point target is given by, where f Kt  is considered as the fast-time frequency. We assume the residual video phase (RVP) has been removed. The time-delay is given by, Due to the target located in the near-field, which is illuminated by spherical waves, the square root of the above time-delay expression cannot be simplified, as in the case of far-field. For a volume target, extending in all three dimensions, propagation. The exponential term in (4) represents a spherical wave emanating from   ', ', x y Z . This term can be decomposed into a superposition of plane wave components [8], as xy . The distinction between the primed and unprimed coordinates is now dropped since the two coordinates coinciding. Typically, the data is sampled at uniform intervals k and, as such, require resampling at equally spaced positions which is referred to as Stolt interpolation [37]. To avoid this interpolation, and similar to the approach in [33], [34], we change the integral variable z dk in (6) to dk and modify the integrant as follows: where   ,, In [33], [34], the summation over fast-time frequencies was carried out prior to performing the inverse FFT over the spatial frequencies, with matched filter used to specify the particular range bin. On the opposite, in this paper, we first perform the inverse FFT with respect to the matched filtered outputs over the spatial frequencies x k and y k , then the results are integrated over all fast-time frequencies (or wavenumbers). In this way, imaging can be considered as an extension of 2-D single-frequency holographic imaging algorithm.
We set and substituting it into (7) where    , IFFT xy kk  represents the 2-D inverse FFT over   , xy kk . The above integral over k can be performed in parallel for all range bins, as also stated in [33] and [34]. The above steps eliminate Stolt interpolations. It is noted from (7) to (9) that the inner integral over k is identical to the single-frequency holographic imaging algorithm at a specified range bin [8]. Therefore, the above 3-D imaging is accomplished by the following procedure: We first perform a 2-D single-frequency based imaging in parallel at different specified range bins with respect to all individual fast-time frequencies. Then, these images are summed coherently over frequencies, and a well-focused 3-D image can be obtained. Therefore, the integral (9) can be considered as an extension of single-frequency holographic imaging algorithm, and is referred to as 3-D holographic imaging algorithm, for convenience. A byproduct of the above problem formulation is that we can freely choose the imaging zone of interest along the z direction. This is different from the traditional 3-D K  algorithm, which generates the image as follows [8]: k k k which is necessary for performing the 3-D inverse FFT. Due to the properties of FFT, imaging along the entire unambiguous range along the z direction must be implemented.

B. Forward and backward operators
Below, we construct a forward operator  which will be used as the sensing matrix for the sparse reconstruction presented in the next section. According to (9), the backward operator †  is as follows, where the term in   stands for the object on which the operator acts. We can then formulate the imaging problem (9) as the following matrix form: where G is a 3-D matrix representing the reconstructed image (13) Clearly, no interpolations are included in (11) or (13). Similarly, The complexity of the K  algorithm and the proposed algorithm is given in Table I  If we choose the resolution to be 5mm 5mm 3cm  (Usually, the range resolution is lower than the transverse resolution due to the limit on bandwidth), then the computation efficiency can be computed to be Although this efficiency appears modest, it is proven to be important when considering the overall optimization problem, as described in the next section.

III. CS APPROACH FOR NEAR-FIELD 3-D IMAGING
In this section, we perform CS-based 3-D MMW imaging based on the aforementioned algorithm incorporating the sensing operator.
The theory of compressive sensing states that sparse signals can be recovered using far fewer samples than that required by the Nyquist sampling. Considering spatial sampling, CS can be used to reduce the number of antenna array elements, thus reducing the overall system cost. In CS, the relationship between the measurements and the sparse image assumes a linear model, namely,   y where 1 l l    α denotes the 1 norm of α and    .
Eq. (14) can be solved efficiently using several linear or quadratic programming techniques [12]. In this paper, we assume that the image is sparse in its canonical basis, and as such, set the sparsifying basis to an identity matrix. For the 3-D MMW imaging, discussed in this paper, the elements of the sensing matrix  are given by, where 1,2, , nN  , 1, 2, , pP  , and 1, 2, , qQ  . The variables N , P , and Q , respectively, represent the total numbers of equivalent frequencies, spatial measurement samples, and pixels of the image scene. Accordingly, the size of sensing matrix  is NP Q  . For the 3-D human body imaging scene, Q could be of the order of 10 5 (assuming the cross-range resolution is 1cm and the range resolution is 3cm).
The number of measurements NP  could be in the order of 10 6 , under Nyquist sampling. Thus, the sensing matrix could be prohibitively large, causing challenges in both storage and processing using a personal computer. The computational complexity for one-time multiplication of this matrix with a 1 Q  vector should be   will be much higher than that of the operators. For a same image size and data size, as illustrated in Table I (15), to construct the model. The MMW images of human body are relatively smooth, and the concealed weapons usually are piecewise-constant objects for which the discrete gradient turns out sparse. This property invites the applications of total variation (TV) compressive sensing techniques [38]. Since there are also very small but lethal objects, such as razor and small lighter, which can be considered as point targets, we utilize the following unconstrained optimization model to reconstruct the image.
where 1  and 2  provide a tradeoff between fidelity to the measurements and noise sensitivity. In this model, the symbols Ĝ , G and S represent the reconstructed image, the target scene, and the scattered data, respectively, all being 3-D matrices. The symbol  denotes the sensing operator, as demonstrated in (13) The TV norm is obtained by, 3  12  1 2 3  1  2 3  1 2 3  1 2  3   1  2  3  1 2 3  1 2 3   12   1  2  3  1  2  3  1 2  3  1 2  3   12   3   1  2 3  1  2 3   3   TV   1  11 , , 1, , , Many optimization algorithms, such as the iterative shrinkage-thresholding based algorithms [39], Least Absolute Shrinkage and Selection Operator (LASSO) [40], can be used to solve the underdetermined equation (16). Because it is not the emphasis of this paper, we just adopt the algorithm in [41] 5 log to minimize (16), which is based on the conjugate gradient method. The gradient of The non-smooth functions in (18), i.e., the 1 norm and the TV norm, can be smoothed by using an approximate function * x x x   , where  is a small positive smoothing parameter. This enables the corresponding gradient to be calculated. The reader can refer to [41].
Another key factor in CS framework is the restricted isometry property (RIP) [12], which is widely used for analyzing the performance of sparse reconstruction algorithms. However, the RIP is often difficult. An alternative to RIP is the mutual coherence of the sensing matrix, which is a more practical approach for evaluating the CS recovery properties. The mutual coherence of a sensing matrix is defined as follows [16]: where i  represents the ith column of the matrix  . If  is small, we state that  is incoherent. However, in this paper, we cannot directly use (19) to evaluate the mutual coherence due to the fact that we represent  as a sensing operator as shown in (13). Instead, as has been shown in [41] and [16], the point spread function (PSF) could be used to measure the mutual coherence of a sampling scheme. It is defined as follows: where i e is the ith vector of the natural basis having 1 at the ith location and zeros elsewhere. The inner product of HH , i j j i      e e e e amounts to the selection of the   , ij -th element of H , which is exactly the inner product , ij  in (19) if  has a matrix form, where the superscript "H" represents the Hermitian operator. Accordingly, a simple measure to evaluate the incoherence is the maximum of the sidelobes of PSF: Although, with the proposed approach, the mutual coherence in (19) cannot be directly computed, the PSF can be readily obtained. It is desirable to have the   PSF , ij ij  to be as small as possible, and have random noise-like statistics for the random undersampling schemes. From (20), we can represent i e as a target scene with only one point target located at its ith element. Then, the imaging procedure is given by column of H  when  is a matrix, as illustrated in Fig. 2. Note that all elements, except for those on the diagonal of H , are sidelobes of the PSFs. The diagonal elements are the main-lobes of PSFs corresponding to the point targets located at all the possible positions in the entire target scene. In this respect, the maximum sidelobe value of ˆi e could be used to evaluate the mutual coherence.
Based on the above discussion, we analyze the effect of different undersampling schemes on the PSF. A uniform-random undersampling scheme is applied to choose the 2-D antenna positions. This approach has been used in [42] for stepped frequency waveform design, and yield a better performance than a totally random undersampling strategy. Specifically, in uniform-random undersampling, we first divide the 2-D antenna positions into a number of non-overlapping groups. Each group consists of the same number of contiguous antenna positions out of which few are selected. In so doing, we provide more uniform illumination of the target scene. In the simulation, we perform 200 independent runs with respect to different selections of the antenna positions according to their corresponding undersampling schemes. The mean of the PSFs is computed. We project the maximum values of the mean PSFs to one dimension, such as the azimuth dimension, as shown in Fig. 3, for both random and uniform-random undersampling schemes. The level of the horizontal line measures the maximal value of the sidelobes and aliasing artifacts. It is evident from Fig. 3 that the artifacts introduced by uniform-random undersampling are lower than those by totally random undersampling. Moreover, the variation of the artifacts of the uniform-random scheme is also lower than that of the totally random scheme. In Fig. 4, we compare the PSFs for different sensing operators, such as the K  based sensing operator and the proposed operator. It is evident that the level and variation of the artifacts and sidelobes for the proposed sensing operator are both lower than those for the K  based operator, which indicates better incoherence. The aliasing artifacts can be removed by using a nonlinear reconstruction technique improving sparsity as introduced in [16]. Next, we construct the sensing operators for CS optimization model for both algorithms, as described in Sec. III. We present the computational time for the Holography-CS and the K  -CS in Fig. 9. The size of the reconstructed image varies from 16      We perform 50 independent runs for each SNR and data ratio to obtain the error means and standard deviations. It is noted from Fig. 10 that the means and standard deviations of RMSEs of the operator-based CS methods are both much lower than those of the Fourier based method (except for the data ratio 0.2, the deviation of the Fourier method is smaller than that of the CS methods). The Holography-CS method has smaller errors than the K  -CS.
On the other hand, the errors of these two methods exhibit a small increase when the data ratio exceeds some extent. It could be caused by the combined effects of mutual coherence and data ratio.

B. Experimental results
For the real measured data, we provide a comparison between the Holography-CS and Fourier-based imaging method. The results of K  -CS method are not given because its computational time is extremely long. We construct a W-band (92GHz to 94 GHz) imaging system, as shown in Fig.   11. The antenna can be scanned at a 2-D planar aperture with 200 200  points. The imaging targets include ten small metal balls and a combination of knife and scissor, representing a point-like target and volume target, respectively, as illustrated in Fig. 12. The imaging results of Fourier-based method (specifically, the K  algorithm) are shown in Figs. 13 and 14, by using the full data set, 50%, and 30% of the data, respectively. Clearly, the random undersampling results in white noise like artifacts in the image domain. These artifacts cannot be removed by the Fourier-based methods. However, with CS, it is possible to remove the aliasing noisy artifacts without degrading the image quality. Figs. 15 and 16 illustrate the imaging results of the Holography-CS method. It is evident that the CS-based method can obtain better imaging results than the traditional Fourier-based method, even when using much reduced data.

V. CONCLUSIONS
This paper proposed a 3-D CS method to the near-field MMW imaging. Based on a 3-D holographic imaging algorithm, we constructed a sensing operator to avoid storing and processing of the large-scale sensing matrix. Most importantly, there are no interpolations for both the forward and backward operators in performing the optimization procedure iterations. We discussed the mutual coherence and PSF when dealing with sensing operator, and suggested a semi-random way of compressing the antenna elements. The proposed 3-D imaging technique has less computations, better performance and improved imaging quality compared to the K  based CS imaging algorithm. Simulations and experimental results also demonstrated that the proposed technique, not only improves over the Fourier-based imaging, but also outperforms the K  based CS method with lower RMSEs.