A Novel Near Field Image Reconstruction Method Based on Beamforming Technique for Real-Time Passive Millimeter Wave Imaging

A Ka-band 1024-channel passive millimeter wave (PMMW) imager (BHU-1024) based on synthetic aperture interferometric radiometer (SAIR) technique has been developed by Beihang University for security screening. BHU-1024 uses linear phased array to obtain resolution in the horizontal direction and uses aperture synthesis to obtain resolution in the vertical direction. This imager is designed for detecting concealed weapons on the human body and operated under the near-field condition of the antenna array. Thus, the conventional direct Fourier imaging theory, which is based on the far-field approximation, can’t be applied any longer. In this paper, a novel near-field image reconstruction method based on the beamforming technique is proposed. In this method, we derive the near-field imaging formula by beamforming theory and demonstrate the feasibility of obtaining spatial brightness temperature distribution by beam focusing. In the design of the beamformer, we further optimize the weight vector to suppress the sidelobe levels of the synthesized beam. The results of simulation and measurement demonstrate that the proposed method is an advantageous, effective imaging method for near-field passive millimeter wave imaging. At present, the proposed method has been applied to the actual passive millimeter wave imaging system BHU-1024.


I. INTRODUCTION
The synthetic aperture interferometric radiometer imaging (SAIR) is widely applied in radio astronomy and Earth observation [1]- [4]. The basic measurement of SAIR is the socalled "visibility function," obtained from the measured cross correlations between the signals received from each pair of antennas. According to the Van Cittert-Zernike theorem [5], the visibility function samples can be utilized to get the brightness temperature within the field of view (FOV) by the inverse Fourier transform.
In the recent past, it has been proved that the SAIR can be applied in passive millimeter wave (PMMW) imaging for security screening [6]- [9]. Beihang University has developed consecutive generations of PMMW demonstrators over the last few years [10]- [12], which are used to verify the ability to detect concealed objects and validate the advantages of high imaging rate and large FOV. At present, an improved PMMW imager using 1024 antenna-receiver channels with 1GHz bandwidth has been developed by Beihang University [13]- [15]. As shown in Figure 1, the 1024-channel PMMW imager (BHU-1024) uses linear phased array to obtain resolution in the horizontal direction and uses aperture synthesis to obtain resolution in the vertical direction. The linear phased array forms a fan beam to provide resolution in the horizontal direction by the beam scanning, and the synthetic aperture provides resolution in the vertical direction at each beam direction.
However, the imaging distance of BHU-1024 is 1-4m, typically in the near-field of the utilized antenna array. In the near-field condition, due to the cross-coupling between azimuth and distance, the conventional Fourier imaging method can't be applied any longer [16]. To overcome this problem, two methods have been reported: One is the modified FFT method, this method modifies the near-field visibility by introducing a correction phase term to obtain the equivalent far-field visibility but the images reconstructed by this method still have obvious noise pollution for the cases of extended targets [17]. Another imaging method is called the G-matrix method, which needs to measure the point source Geometry diagram of 1D near-field synthetic aperture imaging response of the imaging system to form the G-matrix, and then reconstruct image by solving the Moore-Penrose pseudo inverse of the G-matrix [18]- [19]. Nevertheless, for the PMMW imaging system applied to security screening, the measurement of the G matrix will take a long time (hundreds of hours), and the quality of the inverted image is very sensitive to the system errors [20]. Hence, it is difficult to apply the G-matrix method to the actual system. The beamforming technique of the antenna array has been widely used in many imaging applications, such as medical ultrasound imaging [21], seismic imaging [22], and sonar imaging [23]. In this article, inspired by the beamforming technique, we propose a novel near-field SAIR imaging method. We establish the relationship between near-field imaging and general array processing and then steer the weight vector to obtain the brightness temperature at each specified position using the beamforming technique. In the design of the beamformer, to suppress the sidelobe levels of the beam pattern, we optimize the weight vector.
The rest of this paper is organized as follows. In Section II, we deduce the near-field imaging principle of SAIR prove the unavailability to reconstruct millimeter wave image using the near-field visibility function directly. Section III describes the proposed imaging method. Numerical and experimental results are presented and analyzed in Section IV and Section V. Finally, Section VI summarizes the results and conclusions of this paper.

II. THE NEAR FIELD IMAGING
The principle of a synthetic aperture radiometer is to use an antenna array composed of small-aperture antenna elements to achieve equivalent high-resolution imaging results of a large-aperture antenna through related measurements. Taking the case of one-dimensional (1D) interferometry as an example, the geometric relationship is shown in Figure 2.
The temperature distribution placed at a height h, and the near-field visibility function V j,k NF measured by two antennas A j and A k can be written as [24]: where T B (θ) is the brightness temperature of the target, k 0 is the wave number, Ω j,k are the antennas equivalent solid angles, ∆r = r k -r j is the path difference and F j,k are the normalized antenna patterns which observe the image from their own coordinate system. r̃j ,k (τ) is the fringe washing function, and it can be neglected in narrow-band imaging system [17]. The distance r j,k from the antennas to each pixel T B (θ) can be written as a function of the direction cosine (l = sin θ), pixel distance R and antenna position (x j,k , 0)as: The Taylor expansion approximation of (2) can be expressed as: and then the phase difference k 0 ∆r can be written as: where d = x k -x j is the antenna spacing.
In far-field condition (R ≫ d ), the image is at a large distance in comparison with the array dimensions, and the approximation θ j ≈ θ k ≈ θ, r j ≈ r k ≈ R, x k 2 -x j 2 R ≈ 0 ⁄ can be applied. In this case, the far-field visibility function can be expressed as [17]: where u = d/λ is called the baseline, λ is the wavelength.  It can be seen that the far-field visibility function is the inverse Fourier transforms of brightness temperature. For near-field application, the linear Taylor approximation of (5) will result in large approximation errors that indicate the direct Fourier theory is not suited any longer [17] [20]. To solve this problem, the modified Fourier method can be employed to correct the near-field visibility measurement of (5) by introducing the key item e -j(r k -r j ) and adding the farfield phase e -j2πul : However, for an extended source target, the near-field visibility function is related to different positions of the pixel, and the near-field visibility function just can be partially corrected by introducing the phase-corrected item e -j(r k -r j ) . Hence, the residual near-field model error cannot be avoided when the large FOV is required. To illustrate the impact of the phase-corrected item e -j(r k -r j ) on near-field imaging, we carry out some simulation experiments. The main parameters of these simulations are as follows. The imaging distance is 3m, the frequency is 34GHz, the antenna array is 32×32 and the element spacing is 1cm×1cm. The detailed simulation model will be demonstrated in Section IV. Figure 3(a) shows the original scene of brightness temperature distribution. Figure 3(b)-3(c) show the imaging results based on direct FFT method and modified FFT method. It can be seen that in near-field condition, the modified FFT method can significantly improve the image quality compared to the direct FFT method, but there is still obvious noise in the inverted image of extended targets [19].

A. SAIR ARRAY SIGNAL MODEL
As analyzed in the previous section, it is impossible to correct all the near-field visibilities to the far-field case. Considering the SAIR array can be viewed as a general antenna array [25], we can rewrite the near-field imaging formula by near-field beamforming theory. According to the conventional array processing [26], a signal received by an SAIR instrument with M sensor antennas can be expressed as: ) ) ) where r k, m is distance from the kth signal source to the mth antenna. The output of the beamformer y(t) is a weighted sum of the receiver outputs X: where W=[w 1 , w 2 , … w M ] H is the beamforming weight vector. And then the power (variance) of y(t) can be calculated as: where R=E[XX H ] is the covariance matrix of the received signal vector X. As shown in [27] and [28], the visibility samples are equivalent to the elements of the covariance matrix R. Therefore, we can directly use visibility samples to establish the covariance matrix R. For example, R 16,17 (the element at row 16 and column 17 of the matrix R) can be expressed as: 16 where visibility sample V 16,17 denotes the correlation result of the 16th and 17th receiver outputs.
Since the model of SAIR is equivalent to the general antenna array, we can focus the beam to a specified position by steering the weight vector W to obtain the estimation of the power at this position. And for the interferometric radiometer, the power of the signal is interpreted as the brightness temperature [29]. Therefore, through the above analysis, we can use the imaging method based on near-field beamforming technique to reconstruct the spatial brightness temperature distribution.

B. OPTIMIZATION OF THE WEIGHT VECTOR W
According to the analysis in section III.A, the near-field beamforming method can be applied to the near-field synthetic aperture radiometer imaging, and the weight vector W is the key for the near-field beamforming method. Generally, there are two beamforming methods to obtain the weight vector W: One beamforming method is to adaptively control the weight vector according to the received signal data, and the minimum variance distortionless response (MVDR) method is the most popular approach to adaptive beamforming [30]- [31], which minimizes the array output power while maintaining a distortionless main-lobe response toward the desired. However, there are still some drawbacks in the application of this method to near-field SAIR imaging. First, the signal-to-noise ratio (SNR) of the imaging system determines the performance of the adaptive beamforming method, which indicates this method lacks robustness and is sensitive to the system errors. Second, the adaptive beamforming method includes large-scale matrix inversion operations, which is not conducive to real-time imaging. Therefore, the adaptive beamforming method is not suitable for the near-field SAIR imaging system. Another beamforming method is independent of the received signal data, and the most widely used method is the conventional beamforming (CBF) method [26]. The CBF method is also called spatial matching beamforming, which is to steer the weight vector W to maximize the output power at the desired spatial position. Assuming that the desired beam focus point is (h,θ p ), we can design the optimization problem of the CBF method: By means of the Lagrange multiplier method [32], we can get the solution of (12) easily: (13) where W CBF is the CBF weight vector. The CBF is a fast and effective beamforming method, but this method tends to have a high sidelobe level of the beam pattern. The high sidelobe level degrades the contrast between the human body and dangerous objects [33], which is detrimental to the detection of dangerous objects in the security screening. We need to exert more constraints on (12) to suppress the sidelobe levels of the synthesized beam to obtain a lowsidelobe level beam pattern for SAIR imaging. Therefore, we design the weight vector with the following constraints: We keep the beam focus point as (h,θ p ) firstly. Secondly, we use multiple quadratic inequality constraints outside the mainlobe area to keep the sidelobe level below the preset value. At last, we minimize the norm of the weight vector to increase the robustness of the beamformer. The above analysis leads to the following optimization problem:  (14) where Θ SL is the sidelobe region, (h,θ s ) is the field point in Θ SL and the δ s is the prescribed sidelobe level. The optimization problem (14) is a typical convex optimization problem, and it can be solved by the MATLAB CVX toolbox [32] [34]. Once we obtain the optimized weight vector W, the spatial brightness temperature distribution T can be reconstructed by the following formula: where W opt (h,θ i ) H is the optimized weight vector corresponding to the focus point (h,θ i ), R is the covariance matrix constructed from the visibility samples and P is the number of pixels in the brightness temperature distribution scene. In this section, the performance of the proposed method is verified by numerical simulation. In the simulation, we establish the simulation model based on BHU-1024. As is shown in Figure 4, the BHU-1024 consists of 32×32=1024 receiving elements and the element spacing is about 1.2 λ×1.2 λ at 34GHz. Each row of the receiver array consists of four 8-channel receiver modules, and the entire array has a total of 128 receiver modules. To perform the simulation, the specific simulation parameters have been listed in Table  I.

A. POINT SPREAD FUNCTION
The point spread function (PSF) is a good indicator to characterize the system spatial resolution and side-lobe suppression. According to [2], the PSF is equivalent to the imaging result of the point source. In our simulation, two methods are used to calculate the PSF of the antenna array for comparison. The first one is the modified FFT method, and the simulation result is shown in Figure 5(a). The second one is the proposed method, which uses (14) to calculate the optimized weight vector W of the antenna array. To suppress the sidelobe level while maintaining the width of the main lobe, we use the following constraints to calculate the weights: We keep the sidelobe region of the PSF calculated by the proposed method the same as the sidelobe region of PSF calculated by the modified FFT method, so we set the sidelobe region Θ SL to (-90°, -1.4°)∪(1.4°, 90°) . Meanwhile, we set the sidelobe level δ s to -20dB.
The PSF of the proposed method is shown in Figure 5(b). As mentioned in section I, we use the linear phased arrays to obtain resolution in the horizontal direction (x-direction) and use aperture synthesis to obtain resolution in the vertical direction (y-direction). That is, we carry out 1-D synthetic aperture imaging at each beam direction. Therefore, we need the vertical section of the imaging result for further comparison. Figure 6 shows the vertical section of the PSF. The sidelobe level of the modified FFT method and the proposed method are -13.3dB and -20.0dB, and the half power beam width (HPBW) are 1.14° and 1.20°, respectively. It is found that the proposed method effectively reduces the sidelobe level of the PSF of the antenna array. The sidelobe level of the proposed method is much lower than that of the modified FFT method while the HPBW of the two methods are similar.

B. TWO-DIMENSIONAL IMAGING
In this section, two-dimensional (2D) imaging simulation experiments of four imaging methods are performed . Figure  7(a) shows the original brightness temperature distribution, and the number of horizontal and vertical pixels are 80×160. The brightness temperature is about 300 K for the human body and 160 K for the concealed object. The observation distance is set as 3 m and the observation range is 1m×2m. Under this condition, the images are reconstructed by the direct FFT method, modified FFT method, G-matrix method, and the proposed method. The corresponding 2D simulation results are shown in Figures 7(b)-(e). The vertical sections of these reconstructed images at the plane (x = 0 m) are illustrated in Figure 8.
To compare the results of different imaging methods objectively, we calculate the root-mean square error (RMSE) of the imaging results as follows [20]: where X 0 is the original image, X R is the reconstructed image. We give the RMSE results for these reconstructed VOLUME XX, 2017 6 The RMSE results show that the FFT method is not applicable in the near-field case. The modified FFT method has a better performance than the direct FFT method, but there still exist significant errors. The G-matrix method and proposed method have the best performance, but the latter does not require additional measurements, so it is more suitable for engineering applications of actual imaging systems.
For the purpose of comparing the imaging results of these  (17) where D denotes the array size. According to the simulation model based on the antenna array of BHU-1024, the general far-field is 27.8 m, and the absolute far-field is 278 m. Hence, we perform the simulation at the imaging distance of 1 m to 1000 m, and the RMSE results are shown in Figure 9. It can be seen that the direct FFT method is only applicable to farfield conditions, and the modified FFT method has unacceptable error under extreme near-field condition (1m). G-matrix method and proposed method have good performance in near-field and far-field conditions. Finally, the imaging results of these methods tend to be consistent with the increase of distance.
In the above simulation, the accurate imaging distance of the three near-field imaging methods (modified FFT method, G-matrix method and proposed method) is our known condition. However, in the actual imaging system, the measured imaging distance often has errors. To illustrate the  robustness of the proposed method, the distances errors are simulated. As shown in Figure 10, we calculate the RMSE results of modified FFT method, G-matrix method, and proposed method at the imaging distance of 2.5m to 3.5m, and the corresponding focusing distance error is (-0.5m, 0.5m). From Figure 10, we can see that the RMSE of the modified FFT method is relatively large, but the modified FFT method presents a low sensitivity to imaging distance errors. The G-matrix method is extremely sensitive to the imaging distance error, the RMSE result even exceeds the modified FFT method with the increase of imaging distance error. The RMSE result of the proposed method is below to the above two methods, which indicates the effectiveness and robustness of the proposed method.

V. MEASUREMENT EXPERIMENT
To further validate the proposed imaging method, the measurement experiments have been carried out by using the prototype of BHU-1024 developed by Beihang University. The Ka-band 1024-channel prototype with 1 GHz bandwidth is shown in Figure 11(a). The imaging experiments including observing a noise point source and the human with a concealed metal knife are carried out to test the performance of the proposed method. The main parameters of this prototype are listed in Table II.

A. POINT SPREAD FUNCTION
As analyzed in section IV, the PSF is a key tool to characterize the system spatial resolution and side-lobe suppression, and the PSF is equivalent to the imaging result of the point noise source. Figure 11(b) shows the experimental scene of the point noise source, in which we use BHU-1024 to observe the point noise source at a distance of 2m. The point noise source consists of matched load, low noise amplifier, and rectangular horn antenna. The measured PSF results are shown in Figure 12(a)-(b), and the vertical section of the PSF results are shown in Figure 13. The HPBW of the modified FFT method and the proposed method are 1.16° and 1.21°, respectively, which agree well with the system design. The sidelobe level of the modified FFT method and the proposed method are -12.0dB and -18.8dB, which indicates that the proposed method effectively reduces the sidelobe level of the PSF of the antenna array. However, the sidelobe level of measured PSF is worse than expected. This could be a comprehensive result of antenna pattern imperfection, and calibration residual errors.

B. TWO-DIMENSIONAL IMAGING
The real outdoor imaging experiments have been carried out to verify the presented imaging method. Figure 14 shows an optical image of the outdoor experimental scene. In the experiment, a metal plate with a 45incline to the ground is fixed behind the person. The metal plate can reflect the spontaneous radiation from the sky and form a uniform cold background.
As shown in Figure 15(a), the imaging target is a person carrying a concealed metal knife about 3m in front of the BHU-1024. Figure 15(b) shows the imaging result based on the direct FFT method, it can be seen that the metal knife is almost unrecognizable. To further compare the performance of the imaging method, the reconstructed image using the modified FFT method is shown in Figure 15(c), and the imaging result based on the proposed method is shown in Figure 15(d). From the imaging results in Figure 15(d), we can see that the outline of the human body is clear, and the concealed metal knife is in sharp contrast with the human body. Clearly, we can see that the image reconstructed by the proposed method is better than the one reconstructed by the modified FFT method. The image reconstructed by the proposed method has less noise and the contrast between the concealed metal knife and the human body is more obvious.

VI. CONCLUSION
In this paper, we first analyze the synthetic aperture imaging theory under near-field condition, and illustrate the shortcomings of the conventional near-field imaging method. And then we propose a novel near-field imaging method based on the beamforming theory. The numerical simulation and imaging experiment prove its validation. Finally, we have applied this proposed method to the actual passive millimeter wave imaging system BHU-1024.
In future work, we will further optimize the weight vector W to improve the imaging quality and expand the application of the proposed imaging method. In particular, due to the ease of implementation and high imaging quality, the proposed imaging method is well suited for a variety of applications.