Three-Dimensional Microwave Imaging for Concealed Weapon Detection Using Range Stacking Technique

1School of Information Engineering, Inner Mongolia University of Technology, Hohhot 010051, China 2Inner Mongolia Key Laboratory of Radar Technology and Application, Hohhot 010051, China 3School of Electronics and Information Engineering, Beihang University, Beijing 100191, China 4School of Communication and Information Engineering, University of Electronic Science and Technology of China, Chengdu 611731, China


Introduction
Recently, increasing terrors have been motivating the development of various concealed weapon detection techniques for better security monitoring at airports, courthouses, and other buildings [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15].The ability to detect dangerous hidden items is seen as a significant tool in defeating terrorism in today's environment.However, it is monumentally difficult to timely detect concealed weapons carried on an individual with conventional metal detectors.Since microwave possesses a unique property of passing transparently through materials such as heavy clothing, it is possible to identify the presence of hidden threatening objects under clothing with microwave imaging technique.In addition, microwave imaging systems eliminate the use of ionizing radiation such as those present in X-ray systems and, therefore, pose no known health hazard at a low power.In the past few years, an increased interest has been addressed to three-dimensional (3D) microwave imaging [16][17][18][19][20][21][22].The system has the capability of mapping the microwave reflectivity of an individual in 3D space and providing the basis for detection, identification, and further automatic recognition.Lately, the Smith Detection Inc. in England developed a flat-panel millimeter-wave technique for concealed weapon detection for the increasing people screening market [4].Two-dimensional (2D) planar scanning geometries are adopted in all the above systems and, consequently, a 3D reflectivity image of the human body can be formed.
As for the personnel surveillance applications, most of the previous image reconstruction algorithms based on 2D planar synthetic aperture mainly focus on the 3D range migration algorithm (RMA) [23], which requires a 3D STOLT interpolation before the 3D inverse discrete Fourier transform (DFT).The formulation of the 3D RMA is derived by using the stationary phase [24] or the Fourier decomposition of the 3D free space Green function.Here, we discuss yet another interpolation-free reconstruction algorithm which accounts precisely for not only the free space propagation losses but also the wavefront curvature.The algorithm, known as range stacking algorithm (RSA), is usually derived based on wave equations and has been widely and successfully used in seismic signal processing for oil exploration for decades [25].To date, it has only been used in its 2D version in the field of synthetic aperture radar (SAR) imaging.The discussed algorithm is an extension of the work of Soumekh, who primarily developed the 2D RSA for the conventional SAR imaging [26].Here we introduce an alternative viewpoint to reconstruct the image and derive the interpolationfree reconstruction algorithm from wavenumber domain instead of wave equations.This algorithm uses a different matched filtering reference for each range bin and then numerically integrates over the fast-time frequencies before inverse Fourier transform to get the final result.By explicitly performing these steps, the interpolation step required in existing algorithms is not required and, moreover, the final image does not suffer from the truncation and DFT wraparound errors in the wavenumber domain.The algorithm has been successfully extended to 3D microwave imaging for concealed weapon detection and the reflectivity images at several range bins within the illuminated range swath can be easily formed.
The remaining of the paper is organized as follows.An introduction to 3D microwave imaging geometry and an alternative viewpoint for image reconstruction under 2D planar aperture condition are described in Section 2. The extended 3D RSA for implementing the image reconstruction is presented, and the corresponding procedures are addressed in Section 3. Section 4 is dedicated to numerical simulations and experimental results with realistic measurement data.Finally, the conclusions are provided in Section 5.

Planar Scanning Geometry Signal Model and Image Reconstruction.
We begin by considering a 3D imaging geometry for the detection of concealed weapons or other hidden items carried on personnel, as shown in Figure 1.A wide-band stepped-frequency (or frequency-modulated) continuouswave (CW) signal is radiated with a linear array of antennas mounted on a fast mechanical scanner along the -axis or axis.When the scanner is moving along the corresponding axis, the sequentially switched transmitter receivers scan quickly over a large linear aperture to actively illuminate the target, so that a full 2D planar rectangular aperture of data can be gathered.The equivalent sample point of the quasimonostatic transmit and receive antenna pair is assumed to be at the position (0, , V), and the target is assumed to be at the position (, , ) with a spatial distribution of reflectivity function (, , ).Taking the free space propagation losses into consideration, the signal reflected from a given target at position (, , ) can be written as follows when measured at position (0, , V) [2,14,17]: where   ∈ [ min ,  max ] is the fast-time frequency wavenumber and related to the temporal angular frequency  and speed of light  via   = /;  min and  max are the wavenumbers at the minimum and maximum angular frequencies, respectively;  is the distance between the target   and the measurement point (0, , V); that is, where  ∈ [− 0 /2,  0 /2],   is the reference range, and  0 is the range swath along the -axis.
Let  denote the support band of the 3D illuminated area, and then the response measured at the planar aperture transceiver positions at frequency wavenumber   can be formulated as follows: where  =  and | − | ≤   /2 and | − V| ≤   /2 with   and   representing the aperture length along the -axis and -axis, respectively.Note that either of them subtends the corresponding half-power beamwidth in that direction.
The signal model (3) indicates that the measured response comprises echoes of different targets, any one of which has an amplitude inversely proportional to the square of the distance to the target.This amplitude, which decays with the range, is referred to as propagation loss.If left unhandled, it will degrade the image focusing quality.However, according to (3), the propagation losses are difficult to be directly removed in the frequency wavenumber domain.In this paper, we compensate by first transforming (  , , V) to the slant range spatial domain and then multiplying with an amplitude factor  2 ; that is, and, except for some additional computations, (4) brings about the benefit of being able to reduce the coupling between transmitter and receiver and undesired noise beyond the region of interest by a simple range gating in the spatial domain.The computational load and the data volume can be further reduced by subsampling (  , , V) without aliasing, and then the useful signal reflected from the human body could be picked up for subsequent processing at the subsampled set.
After the propagation loss compensation, the 3D microwave reflectivity image of the human body can be obtained by an integral resembling the inverse transform of (4) with respect to the variables   , , and V; namely, where  denotes the surface of the planar synthetic aperture and f(, , ) represents the reconstructed values of the ideal target function (, , ).The image reconstruction can be derived by using the stationary phase.For instance, both Sheen and Soumehk decomposed the spherical-wave exponential term in (5) into a superposition of plane-wave components with Fourier transform.To facilitate the understanding of the image reconstruction, we give an alternative derivation instead of that of (5) in the following subsection.

An Alternative Viewpoint for Image Reconstruction.
Firstly, let us consider the target function reconstruction in the 3D Fourier domain: where   ,   , and   are the frequency wavenumber variables along the -, -, and -axis, respectively;  denotes the support band of this 3D wavenumber domain.In order to reconstruct the 3D reflectivity image of the illuminated target area, the relationship between (  ,   ,   ) and the measured response (  , , V) can be found as follows: where (  ,   ,  V ) is the Fourier transform of the propagation loss compensated signal   (  , , V) with respect to (, V) and (  ,   ,  V ) is the matched filter for correcting the wavefront curvature of all target points with the same reference range   , which has the following form: When the 3D RMA is used, the wavefront curvature errors can be compensated by the known STOLT transformation; namely, where computationally intensive 3D interpolation is required to recover the rectilinear samples of (  ,   ,   ) from available (  ,   ,  V ) data samples.To avoid this, substituting (7)∼( 9) into ( 6) and making a change of variables from where |(  ,   ,  V )| is the Jacobian determinant to accommodate for the change of coordinates which can be formulated as The Jacobian determinant is a slowly fluctuating function and as a result its impact on the 3D image reconstruction is negligible.In the following discussions, we assume that the Jacobian function, or any other amplitude factor, is absorbed International Journal of Antennas and Propagation in the measured response spectrum (  ,   ,  V ).Suppose also that we are interested in forming the 3D image at range bins   s, where   ∈ [− 0 /2,  0 /2].Consider the above reconstruction for (, , ) at a fixed range  =   : It can be derived that the second exponential term in ( 12) is related to a reference signal    (  , , V) (also propagation loss compensated) at the reference range bin  =   via the following expression [25]: where the reference signal    (  , , V) is the response from a reflector located at (, , ) = (  , 0, 0) with    (  ,   ,  V ) representing its 2D Fourier transform with respect to (, V) and * denotes the complex conjugate operator.Then the target function reconstruction at ( =   , , ) can be rewritten as (14) indicates that the spatial distribution of target reflectivity at a specific range bin  =   could be obtained by multiplying the measured response with the reference signal at the same reference range bin   , followed by integrating the result over the available fast-time frequency wavenumbers   , as well as a 2D inverse Fourier transform with respect to (  ,  V ), rather than mapping the (  ,   ,  V ) domain into the (  ,   ,   ) domain through a computationally intensive interpolation as that of (9).It is worth noting that, before transforming to (  ,   ,  V ), the propagation loss compensated signal   (  , , V) is often upsampled along the /and /-axes through zero-padding to avoid any possible aliasing.Thereafter, a parallel implementation of the above procedures can simultaneously form the target function at the individual range bins   s and then stack the outcomes together to form a 3D image of the target function.
Next, we develop the detailed procedures of the 3D range stacking algorithm (3D RSA) for implementing the above reconstructions.

3D Range Stacking Algorithm
where the constant amplitude factor in the Fourier integral is neglected.Comparing ( 15) with ( 14), we can get that the reconstruction for  , ( =   ,   ,   ) is where Practically, the above one-dimensional integral can be implemented by summing over the available discrete values of   and the digital implementation of the Fourier transform or its inverse; that is, the forward/inverse DFT can be realized through efficient fast Fourier transform (FFT) and inverse FFT (IFFT) algorithms.Thus, according to ( 14)∼( 17), the discrete form of the target function reconstruction at ( =   , , ) is expressed as In summary, the reconstruction procedure of the 3DRSA is illustrated in Figure 2, which contains the following steps.

Postprocessing
Step 1 Step 2 Step 3 Step 1. Perform an IFFT on the measured response (  , , V) with respect to the fast-time frequency wavenumbers   .Compensate the propagation losses, and reduce the coupling and noises beyond the region of interest by range gating, and then the compensated signal is upsampled along the and -axes through zero-padding to avoid any aliasing, followed by transforming to the 3D Fourier domain via a 3D FFT.
Step 2. The resultant 3D Fourier domain signal from Step 1 is first matched filtered with the reference signal spectrum at the range bin   ,  *   (  ,   ,  V ).The result is summed over the available fast-time frequency wavenumbers   to yield the marginal Fourier transform of the target function  , ( =   ,   ,   ) at the range   with   =   and   =  V .Then a 2D IFFT with respect to (  ,   ) is performed on  , ( =   ,   ,   ) to obtain the target reconstruction f( =   , , ) at the range bin   .For all available range bins   s, the above calculations are implemented.
Step 3. Stack all outcomes of Step 2 together to form a 3D image of the target function, f(, , ); then apply possible postprocessing procedures such as smoothing filtering, geometry regridding, and object recognition and classification [27,28] and finally display the image.

Sampling Criteria and Resolutions.
According to the Nyquist sampling criterion [29], the sampling interval along some specific direction is inversely proportional to the spatial frequency domain support band of the signal in that direction.
As for the wave propagation direction, the frequency domain support band of the transmitted signal is determined by Thus the following Nyquist sampling constraint should be met in order to avoid aliasing along the slant range direction: where Δ is the frequency sampling interval and  max denotes the maximum observed distance within the antenna beam pattern.
As for the sampling criterion along the / direction, the samples are realized via a real or synthetic antenna array along this direction, and the frequency wavenumber support band Ω   depends on the antenna beamwidth of the individual antenna Φ  through Thus the Nyquist sampling spacing Δ along the / axis for gathering the echoed signal must satisfy the following inequality: where  min is the smallest wavelength at the highest frequency  max and the antenna beamwidth Φ  is determined by the antenna size along the corresponding direction.The worst-case situation corresponds to the case of Φ  = 180 ∘ for which the Nyquist sampling spacing in the / dimension is International Journal of Antennas and Propagation As for the case along the / direction, a similar conclusion can be drawn that the Nyquist sampling criterion in the / dimension takes the form where Φ  denotes the antenna beamwidth along the / direction.
The frequency wavenumber support band and the derived sampling constraints depend on the dimensions of the volume occupied in the 3D wavenumber domain.Therefore, the spatial resolutions of the resulting 3D reflectivity image f(, , ) are determined by the signal bandwidth, the center frequency, and the effective synthetic aperture.The achievable −3 dB resolutions along the -, -, and -axes, are approximated, respectively, by where   is the wavelength at the center frequency   and the units of   ,   , and   are all meters.In practice, due to the sidelobe structure of the radiation patterns, the final resolutions are usually slightly poorer than those theoretical values given by ( 25)∼ (27).In spite of this, both simulation and experimental results provided in the next section show that the 3D RSA is well suited to reconstruct 3D microwave imaging in concealed weapon detection at a high-precision manner.

Numerical Simulation.
In the computer simulations, a human body model consisting of 816 unit-amplitude point scatterers is numerically generated; as in Figure 3, the model has pixel spacings of 2.5 cm and 5.0 cm on the horizontal plane and in the height direction, respectively; therefore it is practically confined within a box of dimensions 25.0 cm (axis) × 50.0 cm (-axis) × 170.0 cm (-axis) (the definition of the coordinate is in accordance with that of Figure 1, with -, -, and -axes representing the range, azimuth, and height directions, resp.).Assuming that the human body is facing the planner scanning antenna aperture and considering that the microwave signals are not capable of penetrating the human body, therefore only 486 points from the model representing the surface of the front part of the body are picked out to generate the simulated radar echoes.The measurement parameters used in the simulations are listed in Table 1.Though a rectangular radiation pattern is assumed for each antenna element for simplicity, we will see in the latter real data experiments that the effects of the sidelobes from the practical antenna radiation pattern has a trivial influence on the performance of the proposed imaging algorithm.After performing the processing chain of Figure 2, a 3D image of the human body is reconstructed, as in Figure 4(a).The proposed 3D imaging algorithm shows good imaging performance; for example, the reflectivity values as well as the positions of the point targets are all in perfect agreement with those of the simulated model.To quantitatively evaluate the performance, a point target with the coordinate (, , ) = (0.5 m, 0, 0) is simulated independently, the reconstructed point target is shown in Figures 5 and 6, and the figures of merits like resolution, peak sidelobe ratio (PSLR), and integrated sidelobe ratio (ISLR) are evaluated and compared with their theoretic values, which are indicated in Table 2.All the above measured parameters indicate that the developed algorithm provides a proper and accurate implementation of the 3D reconstruction.
Actually, as for the people screening system, the resolution along the azimuth and height direction can be easily achieved, while being compared with that along the wave propagation (range) direction, so that the main lobe of the reconstructed impulse response along the range direction is much wider than the main lobes along other directions.Accordingly, if the 3D image is projected onto a 2D plane directly, the energy of the main lobe along the range direction would possibly be spread out across other directions due to the shift-varying property of the impulse response.To overcome this problem, the 3D image is collapsed using a maximum-value projection into a fully focused 2D image of the target for display, rather than being directly projected.

Real-Field Experiment Result and Discussion
. A wideband imaging system has been developed that utilizes a scanned linear array of antennas to quickly gather the image data for subsequent reconstruction and display.The system consists of a vertically scanned linear antenna array which consists of 281 virtual elements which are formed by the multiple-input multiple-output (MIMO) technology [30,31].The antenna array is driven by a wide-band millimeterwave transceiver which transmits and receives step-frequency continuous-wave (SF-CW) signals and is demodulated to baseband through analogue dechirping.The array antenna is mounted on a large fast horizontal mechanical scanner, which allows a full-image data set to be gathered in less than 1 s.During a scan, the array is sequentially switched at a very high speed to gather each column of data as the scanner is moving from left or right, thereby gathering a full 2D aperture of data.
The system operates in Ka band and each antenna element operates in the HH polarization state and has a highest sidelobe level of −25 dB.A typical set of measurement parameters are listed in Table 3, where the current target scene consists of a human body model carrying a concealed knife near the left waist and a small metallic sphere around the neck.The wide-band data are sampled vertically across the array and horizontally over the synthetic aperture and are digitized by an A/D converter for subsequent storage in the computer.The proposed 3D range stacking algorithm is programmed with CUDA C language to run on an NVIDIA TITAN X GPU platform parallel, which allows it to be executed at very high speed with an average period of 0.7 seconds to generate a frame of 3D image of the human body, compared to the average period of 5 seconds or more for a typical STOLTinterpolated counterpart.So the proposed algorithm is more  Figure 7(a) shows the optical image of the human body model to be imaged.A 2 cm diameter metallic sphere fixed by the adhesive tape is exposed near the neck in the front of the body; however, a knife composed of a stain-less steel blade and a plastic handle is also concealed under the model's jacket near its left waist.Figure 7(b) is the front view of the 3D microwave image result.Due to the high sensitivity of the transceiver design, as well as the highprecision focus of the imaging algorithm, a significantly high dynamic range is apparent in the microwave image.It can be seen that microwaves readily pass through the clothing materials and are reflected by the body and the concealed knife.In this experiment, the 3D range stacking imaging techniques incorporating large apertures and wide angle illumination achieve the high-resolution reconstruction of the targets, and the theoretical resolution for the image is approximately 5 mm, where the differences in the shapes and reflectivity of the targets, such as the metallic sphere, the knife handle, blade, and the model's body, are readily apparent in the reconstructed imagery.Given the remarkable quality of the reconstructed image, it can also be concluded that the effects of the sidelobes of the antenna radiation pattern on the final image are negligible.
It is worth noting that, given a sufficient small point-like target (as compared to the system resolution), the resolution, PSLR, and ISLR can be measured from the reconstructed image, and this kind of point target usually works as a calibrator to calibrate any errors (such as the multichannel amplitude-phase errors of the array) of the imaging system.Other metrics of the image may include the image contrast relating the level of the signal to noise-clutter ratio, which is important to enhance the succeeding concealed weapon detection performance (such as the detection probability and the precision).
This experiment demonstrates microwave imaging can effectively view objects through barriers that are opaque in the optical band.It also demonstrates that the proposed 3D imaging algorithm can achieve the high resolution and image fidelity.

Conclusions
This paper presented a 3D imaging algorithm especially tailored for accurately reconstructing the 3D image of the human body for concealed weapon detection under planner scanning geometry.The algorithm, which is referred to as 3D range stacking algorithm, allows for the propagation losses compensation, wavefront reconstruction, and aliasing mitigating and without interpolations.Moreover, it suffers no approximations and truncation errors.The essence of the algorithm is the multiplication of the reference signal spectrum at each specific range bin, followed by the coherent summation over the whole frequency band, instead of the 3D frequency interpolation.In addition, the sampling criteria and the resolutions for planner scanning geometry are also derived.The imaging results using simulated data and experimental data demonstrated the high resolution and good image fidelity of the proposed algorithm.

Figure 2 :
Figure 2: Flowchart of 3D range stacking algorithm for concealed weapon detection.

Figure 3 :
Figure 3: Numerically generated human body model with an ensemble of 816 point scatterers, where 486 of the scatterers from the front of the body (red points) are utilized to generate the simulated radar echoes, whereas the left ones (blue crosses) represent the body part that cannot be penetrated by microwave: (a) 3D view; (b) side view; (c) front view.

Figure 7 :
Figure 7: Photograph of the human body model carrying a concealed knife and a metallic sphere (a); orthogonal projection of the 3D microwave image result onto the azimuth-elevation (-) plane for display (b).
3.1.Reconstruction Procedure.Let  , (,   ,   ) denote the marginal Fourier transform of the reconstructed target function f(, , ) with respect to the variables (, ), and then the target function at range  =   can be recovered by a 2D inverse Fourier transform of  , ( =   ,   ,   ) with respect to (  ,   ), which is

Table 1 :
Measurement parameters used in the simulation experiment.

Table 2 :
The measured values of PSLR, ISLR and resolutions.

Table 3 :
Measurement parameters of the real data experiment.