# Modified wavenumber domain algorithm for three-dimensional millimeter-wave imaging.

1. INTRODUCTION

The increasing threat of terrorism has made personnel surveillance increasingly important. Conventional security systems have been effective for detecting metal targets. However, some modern threats, such as plastic or ceramic handguns and knives, as well as extremely dangerous liquid explosives, etc., cannot be detected with conventional security technology, such as metal detectors. X-ray systems can present an effective solution to the problem. However, potential health risks due to exposure to ionizing radiation make it less acceptable to public safety.

Thus, airports, the traveling public, and other secure locations have a need for a new personnel surveillance system that is capable of detecting concealed body-worn threats that are undetectable with conventional security technology. One of the choices is to use millimeter waves which are capable of penetrating common clothing barriers to form an image of a person as well as any concealed items with different reflectivity or emissivity [1-3]. Millimeter waves are nonionizing and, therefore, pose no known health hazard at moderate power levels.

The combination of frequency-modulated continuous-wave (FMCW) and MMW imaging techniques leads to compact, light-weight, cost-effective, low-power operation, and high-resolution imaging systems, which are especially suitable for security and detection application. FMCW signal has been widely used in imaging radar [4-7]. Typically, FMCW systems use long pulses. Therefore, conventional SAR imaging algorithms that use the stop-and-go approximation need to be modified for FMCW image processing. Several conventional algorithms, such as range-Doppler algorithm, frequency-scaling algorithm, and chirp-scaling algorithm, etc., have been modified to focus the FMCW SAR data [8-10]. Especially, an accurate received two-dimensional (2-D) signal model was proposed in [6], in which the effect of the variation of the instantaneous slant range on the transmitted and received signal was accurately represented, and a range-azimuth coupling term was formulated for the first time in the FMCW SAR community.

The 3-D MMW imaging techniques for concealed weapon detection were discussed in [1,11-14], etc., whereas the effects of the continuous motion of the array on the transmitted and echoed signal were not addressed. In this paper, we begin with an analysis of the 3D signal model which accurately includes the effects of the continuous motion of the array on the transmitted and received signal during the pulse time. The main effects of the motion on the image result turn out to be a range walk and out of focus. These effects can be corrected by a matched filter multiplication in the 3-D frequency domain.

After the corrections, the Stolt mapping can be implemented to compensate for the range curvature of all scatterers by an appropriate warping of the wavenumber domain backscattered data. For the 3D imaging, the Stolt mapping is implemented by a 1-D interpolation process in the range-frequency domain, which is repeated for every azimuth-frequency and elevation-frequency point. The complex image can then be reconstructed by the 3-D inverse FFT (IFFT). However, the repeated interpolation process reduces the computational efficiency of the algorithm. In [15,16], the Stolt interpolation and IFFT were substituted by a nonuniform fast Fourier transform (NUFFT). The NUFFT approaches have been developed to overcome the limitation of equally spaced sampling needed by FFT [17-23]. NUFFTs have wide applications in magnetic resonance image reconstruction [24], ground penetrating radar with the range migration algorithm [16], synthetic aperture imaging radiometers [25], and near-field imaging [26]. In this paper, the 3-D NUFFT is used to replace the Stolt interpolation and the followed 3-D IFFT to improve the computational efficiency.

This paper is organized as follows. In Section 2, an accurate backscattered signal model is derived with the corresponding spectrum analysis. The 3-D imaging procedure is given in Section 3. Section 4 shows the simulation results to verify the effectiveness of the motion compensation. The comparison of the computational time between the interpolation-FFT and NUFFT are also shown. Section 5 summarizes the conclusions.

2. FMCW THREE-DIMENSIONAL SIGNAL MODEL AND SPECTRUM ANALYSIS

This section derives an analytical model of the FMCW backscattered signal in the 3-D spatial-frequency domain. We consider the 3-D imaging geometry, as shown in Figure 1. The 1-D antenna array aligns with the x direction, and the array is scanned along the minus y direction.

For an FMCW transceiver, the transmitted signal can be expressed as [27]:

[S.sub.T] (t) = exp [j2[pi] ([f.sub.0]t + 1/2 K[t.sup.2])], (1)

where [f.sub.0] is the center frequency, t is the time variable varying within one cycle of signal transmitting, and K is the frequency sweep rate of the transmitted signal.

The signal is transmitted at an arbitrary time [rho]. And the corresponding instantaneous range between the antenna element and the target is R ([rho]). Let the time [[rho].sub.d] be the round-trip delay time. Then the backscattered signal is received at time t + Td with the corresponding instantaneous range R([rho] + [[rho].sub.d]). Thus, the round-trip delay time can be given by

[[rho].sub.d] R ([rho]) + R([rho] + [[rho].sub.d])/c (2)

where

R([rho]) = [square root of [(x'-x).sup.2] + [(y'-y).sup.2] + [(z'-z).sup.2]] = [square root of [(x'-x).sub.2] + ([y.sub.0]-[[v.sub.[rho]]-y).sup.2] + [(z'- z).sup.2]]0, R([rho] + [[rho].sub.d]) = [square root of [(x'-x).sup.2] + [(y'-y).sup.2] + [(z'-z).sup.2]] = [square root of [(x'-x).sup.2] + [[[y'.sub.0]-v ([rho] + [[rho].sub.d])- y].sup.2] + [(z'-z).sup.2],

and c is the speed of light. (x, y, z) is the location of a point target, and (x', y', z') is the position of the antenna element, both within the same coordinates, as shown in Figure 1. And y' = [y'.sub.0]-v[rho], where [y'.sub.0] is the initial position of the array, v is the scanning velocity, and t = [nT.sub.x] + [mT.sub.y] + t is the continuous time; [T.sub.x] is the period of signal transmitting along the array elements, assuming equal to the signal duration time of one cycle, thus t [member of] [0,[T.sub.x]), and [T.sub.y] is the period of sampling time along y direction; n = 0, 1, ..., N-1, and m = 0, 1,..., M-1, where N is the number of the array elements (corresponding to the number of samples along x direction), and M is the number of samples along y direction.

[FIGURE 1 OMITTED]

Due to the fact that the antenna array is near the target for the personnel surveillance system, (2) can be accurately approximated as

[[rho].sub.D] [approximately equal to] 2R([rho])/c

Omitting the time scaling influences on the envelope, the backscattered signal from a point target is given by

[S.sub.R] (x',y'; t) = [sigma] (x,y,z) [S.sub.T] (t-[[rho].sub.d]) (4)

where [sigma] (x, y, z) is the backscattering coefficient of the point target.

The dechirp-on-receive technology is generally used in the FMCW systems [27]. The received signal and a delayed version of the transmitted signal are mixed in order to reduce the required sampling rate. The intermediate frequency signal after mixing can be written as

[S.sub.IF] (x',y'; t) = [sigma](x, y, z)exp[-j2[pi][f.sub.0]([[rho].sub.d]- [[rho].sub.c])]exp[-j2[pi]K([[rho].sub.d]- [[rho].sub.c])([rho]-[[rho].sub.c])] exp [-j[pi]K ([[rho].sub.d]-[([[rho].sub.c]).sup.2] (5)

where [[rho].sub.c] is the delayed time of the transmitted signal. The last exponential term of (5) is known as the residual video phase (RVP), which is caused by the difference between any two displaced linear FM waveforms [27]. The compensation for RVP is conducted in the frequency domain, which is not the emphasis of this paper. Assuming the RVP has been removed in the following, i.e.,

[S.sub.IF] (x', y'; t) = [sigma] (x, y, z) exp [-j 2[pi][f.sub.0] ([[rho].sub.d]-[[rho].sub.c])] exp [-j2[pi]K ([[rho].sub.d]-[[rho].sub.c]) ([rho]-[t.sub.c])] (6)

Using the substitution/= K(t-[[rho].sub.c]), (6) can be rewritten as

[S.sub.IF] (x', y';f) = [sigma] (x, y, z) exp [-j2[pi] (f + [f.sub.0]) ([[rho].sub.d]-[[rho].sub.c])] (7)

Substituting t = [mT.sub.y] + [nT.sub.x] + t = [[rho].sub.m] + [[rho].sub.n] + t into (7) yields

[S.sub.IF] (x', y', [??]; [[rho].sub.m] [[rho].sub.n]; f) = [sigma]([??])exp[- j4[pi](f + [f.sub.0]) (T([rho])/c- [R.sub.c]/c = [sigma] ([??]) exp[-j4[pi](f + [f.sub.0])(R([[rho].sub.m] + [[rho].sub.n] + t/c-[R.sub.c]/c] (8)

where [??] = (x, y, z) and [R.sub.c] = [ct.sub.c]/2.

Then the 2-D Fourier transform of (8) with respect to the spatial variables x' and [y'.sub.m] ([y'.sub.m] = [v[rho].sub.m] = vmTy) is given by

[S.sub.IF] ([k'.sub.x],[k'.sub.y], f;[[rho].sub.n])

= 1/v [integral] [S.sub.IF](x', y', [??]; [[rho].sub.m], [[rho].sub.n];f) exp[- [jk.sub.x']x') exp(- [jk.sub.y'][y'.sub.m]) dx'[dy'.sub.m]

= 1/v[sigma]([??])[integral] exp[-j[PHI]([k.sub.x'], [k.sub.y'], f; x', [y.sub.m'], [[rho].sub.n])]dx'[d'.sub.ym] (9)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and [k.sub.r] = 2[pi](f + [f.sub.0])/c.

The integral in (9) can be solved by means of the principle of the stationary phase method [27, 28]. At the point of stationary phase, the first partial derivatives of the phase [PHI] ([k.sub.x]',[k.sub.y]', f;x',[y'.sub.m], [[rho].sub.n]) are equal to zeros, such as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (10)

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (11)

Solving (10) and (11) for the stationary phase points, i.e., [y'.sub.m0] and [x.sub.0]', respectively, yields

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (12)

and

[x.sub.0]' = x-[k.sub.x]'(z-z')/[square root of [4k.sup.2.sub.r]- [k.sup.2.sub.x']-[k.sup.2.sub. y']. (13)

Then the 3-D spectrum of the backscattered signal can be expressed as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (14) where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (15)

Based on the relation between frequency and time of FMCW signal, we have t = f/K + [2R.sub.c]/c. Consequently, the phase of (14) can also be expressed as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (15)

where [R.sub.0] =-z'.

Note that the first three terms of (15), i.e., [k.sub.x']x, [k.sub.y']y and [square root of [4k.sup.2.sub.r]-[k.sup.2.sub.x']-[k.sup.2.sub.y']z are linearly dependent on the target position x, y and z, respectively. Thus the image can be reconstructed by using inverse Fourier transform. The fourth term [k.sub.y'][y.sub.0]' is introduced by the initial position of the array. The term [k.sub.y']v[[rho].sub.n] represents the phase variation corresponding to the elevation-range walk caused by the motion of the nth antenna element. [k.sub.y']vf/K is a space-invariant term also caused by the motion of the array within the signal duration time. The terms [k.sub.y']v([2R.sub.c]/c) and [2k.sub.r][R.sub.c] refer to the constant azimuth and range shifts, respectively, and are introduced by the dechirp-on-receive approach. The last term [square root of [4k.sup.2.sub.r]-[k.sup.2.sub.x']- [k.sup.2.sub.y'] [R.sub.0] represents the constant range shift.

3. THREE-DIMENSIONAL IMAGING PROCEDURE

Based on the spectrum analysis of the 3-D FMCW backscattered signal in Section 2, we present a wavenumber domain algorithm for the image reconstruction. The wavenumber domain algorithm (or called range migration algorithm) [27-29] is an accurate algorithm to process the near-field data.

The processing steps of the algorithm are shown in Figure 2. Removing RVP is not the emphasis of this paper. So it is not addressed. As can be seen from Figure 2, the phase corrections are the key parts for the image reconstruction. To remove the unwanted phase in (15), a 3-D matched filter is given by

[G.sub.F] ([k.sub.x'],[k.sub.y'],[k.sub.r];[[rho].sub.n], [R.sub.ref]) = [G.sub.F1]([k.sub.y'],[k.sub.r]; [[rho].sub.n]) x [G.sub.F2] {[k.sub.x'], [k.sub.y'],[k.sub.r]; [R.sub.ref]), (16)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and

[G.sub.F2] ([k.sub.x'],[k.sub.y'], [k.sub.r]; [R.sub.ref]) = exp[j([square root of [4k.sup.2.sub.r]-[4k.sup.2.sub.y']- [k.sub.r]; [R.sub.ref]]],

where [R.sub.ref] is the reference range for focus processing, which is generally defined as the range from the scene center to the aperture center.

Note that the term [G.sub.F1] must be first multiplied exactly after the elevation Fourier transform and before the azimuth Fourier transform, as shown in Figure 2, due to the fact that the term [k.sub.y']V[[rho].sub.n] exists in the elevation-frequency domain and the azimuth-time domain. Then the term [G.sub.F2] is multiplied to correct the range curvature of all scatterers having minimum range [R.sub.ref]. Residual range curvature is still present on those scatterers with a minimum range that is different with [R.sub.ref] [27]. The residual range curvature can be removed by Stolt mapping, which will be shown in the next.

After matched filtering, the remaining signal becomes

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (17)

where

[[PHI].sub.F] ([k.sub.x'],[k.sub.y'], [k.sub.r]; [R.sub.ref]) = [k.sub.x']x- [k.sub.y']y + [square root of [4k.sup.2.sub.r]-[k.sup.2.sub.x']- [k.sup.2.sub.y']z] + [square root of [4k.sup.2.sub.r]-[k.sup.2.sub.x']- [k.sup.2.sub.x']]([R.sub.0]-[R.sub.ref]). (18)

The minus in "-[k.sub.y']y " is introduced by the scanning direction of the array. The distinction between the primed and unprimed coordinate systems can now be dropped since the coordinate systems coincide.

[FIGURE 2 OMITTED]

As can be seen from the phase [[PHI].sub.F], (17) suffices for the image reconstruction if the data are defined continuously in [k.sub.x], [k.sub.y], and [k.sub.r]. However, for practical imaging configurations, the data are discretely sampled at uniform intervals of position and frequency. After 2-D Fourier transform with respect to x and y, and matched filtering, a sampled version of [S.sub.F] ([k.sub.x], [k.sub.y], [k.sub.r]; [R.sub.ref]) is obtained, which are uniformly spaced in [k.sub.x], [k.sub.y] and [k.sub.r]. However, the samples are nonuniformly spaced in [square root of [4k.sup.2.sub.r]-[k.sup.2.sub.x]- [k.sup.2.sub.y]]. In order to reconstruct the image via 3-D IFFT in (17), the data [S.sub.F] will need to be resampled to uniformly spaced positions in [4k.sup.2.sub.r]-[k.sup.2.sub.x]-[k.sup.2.sub.y]. This is accomplished by Stolt interpolation, which converts [S.sub.F] into a 3-D linear phase grating for targets at all ranges [27]. The appropriate change of variables is [square root of [4k.sup.2.sub.r]-[k.sup.2.sub.x]-[k.sup.2.sub.y]] [right arrow] [k.sub.z] (19)

where [k.sub.z]s are uniformly spaced.

After Stolt interpolation, (18) becomes

[PHI]F ([k.sub.x], [k.sub.y], [k.sub.z]; [R.sub.ref]) = [k.sub.x]x-[k.sub.y]y + [k.sub.z]z + [k.sub.z] (R-[R.sub.ref]) (20)

The aforementioned description is focused on a point target. For a volume target, we have

[S.sub.F] ([k.sub.x],[k.sub.y],[k.sub.z]; [R.sub.ref]) = [integral] (1/v) f (x, y, z) exp [-j[k.sub.z] ([R.sub.0]- [R.sub.ref])]

exp [-j ([k.sub.x]x-[k.sub.y]y + [k.sub.z]z)] dxdydz (21)

where f (x, y, z) is the backscattering coefficient of a volume target.

Clearly, a 3-D IFFT can compress the signal in range, azimuthrange and elevation-range, i.e.,

f (x,y,z) = v [integral] [S.sub.F] ([k.sub.x], [k.sub.y], [k.sub.z] ; [R.sub.ref]) exp [j[k.sub.z] ([R.sub.0]- [R.sub.ref])] exp [j ([k.sub.x]x-[k.sub.y]y + [k.sub.z]z)] [dk.sub.x] [dk.sub.y] [dk.sub.z] (22)

The constant factor relative to 2n is ignored in (22).

However, the change of variables and the interpolation procedures will introduce approximation errors [16] and are time-consuming. In this paper, we substitute the Stolt interpolation and IFFT with the 3D inverse NUFFT to reconstruct the image. The NUFFT idea was first presented by Dutt and Rokhlin [17]; later a new approach was proposed by Liu and Nguyen, using the regular Fourier matrices [18, 19]. A minmax method for fast and accurate approximation of the NUFFT was developed by Fessler and Sutton [21], which is optimal in the min-max sense of minimizing the worst-case approximation error over all signals of unit norm.

Here, the NUFFT algorithm developed by Fessler and Sutton is implemented. For the sake of simplicity, we consider only the 1-D case:

[x.sub.n] = [M-1.summation over (m-0)] [y.sub.m] exp (j[[omega].sub.m]n, n = 0,1. (23)

where [y.sub.m] represents the nonuniform frequency domain sample, [[omega].sub.m] is the nonuniform digital frequency in radian, and [x.sub.n] is the uniform output signal in time domain.

The computation of (23) can be summarized as follows:

1) The "gridding" step: [X.sub.k] = [m-1.summation over (m=0)] [v.sub.mk][y.sub.m], where [V.sub.mk]s denote interpolation coefficients. This step requires O (JM) operations, where J is the number of the nearest neighbors to wm.

2) K-point inverse discrete Fourier transform (DFT) of [X.sub.k] scaled

by K: [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], 0 [less than or equal to] n [less than or equal to] N-1. The output

signal [[??].sub.n]s are the first N signal values. This step requires O (K [log.sub.2] N) operations by using the reduced FFT.

3) Scaling each [[??].sub.n] by [s.sup.*.sub.n] to get [x.sub.n], where the sns are the scaling factors, and "*" denotes conjugate. The overall operation count per NUFFT is O(JM) + O(K[log.sup.2] N). Typically, we choose K = 2N and J [less than or equal to] 10.

Therefore, the overall computational requirements are similar to an FFT but with a larger constant. For more details, the readers are referred to [21]. The computational complexity of the 3-D NUFFT is known as O ([N.sup.3] [log.sub.2] N).

The flow diagram of the image reconstruction process is demonstrated in Figure 2 in which the 3-D inverse NUFFT can also be replaced by 1-D inverse NUFFT and 2-D IFFT due to the fact that [k.sub.x]s and [k.sub.y]s are uniformly spaced.

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

4. SIMULATION RESULTS

In this section, simulations are performed to validate the performance of the proposed method. The simulation parameters are listed in Table 1.

For a first validation of the derivation, the frequency domain data in (14) is multiplied by the conventional and the proposed matched filters, respectively, for a target in the scene center.

Under the stop-and-go approximation, the conventional matched filter [27] is given by

[G.sub.F] ([k.sub.x],[k.sub.y],[k.sub.r]; [R.sub.ref]) = exp [j([square root of [4k.sup.2.sub.r],-[k.sup.2.sub.x]- [k.sup.2.sub.y], [R.sub.ref] + [k.sub.y] [y.sub.0]'-[2k.sub.r][R.sub.c], (24)

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

Figure 3 shows the phase of the signal after matched filtering with a reference function that makes the stop-and-go approximation, i.e., the phase of the remaining signal after (14) is multiplied by (24). The phase of the remaining signal in the elevation-frequency and range-frequency domain, i.e., in the [k.sub.y]-[k.sub.r] domain, is shown in Figure 3(a), and the phase in the elevation-frequency and azimuth-frequency domain, i.e., in the [k.sub.y]-[k.sub.x] domain, is illustrated in Figure 3(b). The phase in the azimuth-frequency and range-frequency domain is not demonstrated due to the fact that it is not affected by the motion of the array, as shown in Figure 1.

Figure 4 illustrates the phase results produced by filtering the signal (14) with the proposed function in (16), i.e., the phase of (17).

[FIGURE 8 OMITTED]

If the signal (14) is completely compensated for by the matched filter, the phase of the remaining signal, i.e., the phase of (17), should be constant and does not vary with [k.sub.x], [k.sub.y] and [k.sub.r] in the case that the simulated point target is in the scene center, i.e., (x, y, z) = (0, 0, 0) and [R.sub.0] = [R.sub.ref]. Clearly, only in Figure 4 does the phase in the support band is practically constant, indicating that the matched filter does match the received signal.

Figures 5(a) to (c) show the slices of the reconstructed 3-D image of a point target located in the scene center. These slices correspond to the x-y (with z = 0), y-z (with x = 0), and x-z (with y = 0) section planes, respectively, with the matched filtering function that makes the stop-and-go approximation. Clearly, the motion of the array has a significant effect on the azimuth-elevation (the x-y plane) image, as shown in Figure 5(a). It can be seen that the nulls of the image are raised which make the slice image, i.e., the 2-D impulse response (or called point spread function) departure from the ideal 2-D sinc function. In addition, a range walk exists in the elevation direction which is caused by the motion, as shown in Figures 5(a) and (b). The image in the x-z plane is not affected by the motion, as shown in Figure 5(c).

[FIGURE 9 OMITTED]

Figures 6(a) to (c) demonstrate the corresponding image results produced by the proposed matched filtering function. From the comparison of Figures 5 and 6, the improvements due to the proposed correction are visible. The dynamic range of the images is 30 dB.

The next simulation is intended to consider a cubic target consisting of 8 scatterers distributed at the vertices of the cube as shown in Figure 7. All scatterers have the same radar cross section (RCS): 1 m2. Figure 8 shows the projection of the images onto the three main planes, with the matched filtering function that makes the stop-and-go approximation. The image results produced by the proposed matched filtering function are shown in Figure 9. The intersections of the dashed lines represent the real positions of the scatterers. The defocusing and range walk are clearly seen in Figure 8, which are the same as shown in Figure 5.

Figure 10 shows the comparison of the computational time (CPU time) between the Stolt interpolation-based method and the NUFFT-based method. It is evident that the NUFFT-based method is much faster than the Stolt interpolation-based method, especially with the increase of the sample points at the 2-D aperture.

[FIGURE 10 OMITTED]

5. CONCLUSIONS

A wavenumber domain algorithm with motion compensation to process the near-field 3-D imaging has been developed. In this algorithm, the effects of the array motion within the signal duration time have been described and processing solutions were given. The unwanted terms caused by the motion have been removed and the focusing performance is improved. To reduce the computational time of the image reconstruction process, we substitute the Stolt interpolation and FFT with a 3-D NUFFT.

The proposed algorithm has been validated with simulations. The algorithm can also be extended to process the 3-D imaging for cylindrical scanning geometry.

Received 24 November 2011, Accepted 31 December 2011, Scheduled 11 January 2012

REFERENCES

[1.] Sheen, D. M., D. L. McMakin, and T. E. Hall, "Three-dimensional millimeter-wave imaging for concealed weapon detection," IEEE Trans. on Microwave Theory and Techniques, Vol. 49, No. 9, 1581-1592, Sep. 2001.

[2.] Appleby, R. and R. N. Anderton, "Millimeter-wave and submillimeter-wave imaging for security and surveillance," Proceedings of the IEEE, Vol. 95, No. 8, 1683-1690, Aug. 2007.

[3.] Yeom, S., D. Lee, H. Lee, J. Son, and V. P. Gushin, "Distance estimation of concealed objects with stereoscopic passive millimeter-wave imaging," Progress In Electromagnetics Research, Vol. 115, 399-407, 2011.

[4.] Huang, Y., P. V. Brennan, D. Patrick, I. Weller, P. Roberts, and K. Hughes, "FMCW based MIMO imaging radar for maritime navigation," Progress In Electromagnetics Research, Vol. 115, 327-342, 2011.

[5.] Meta, A., P. Hoogeboom, and L. P. Ligthart, "Signal processing for FMCW SAR," IEEE Trans. on Geoscience and Remote Sensing, Vol. 45, No. 11, 3519-3532, Nov. 2007.

[6.] Wang, R., O. Loffeld, H. Nies, S. Knedlik, M. Hagelen, and H. Essen, "Focus FMCW SAR data using the wavenumber domain algorithm," IEEE Trans. on Geoscience and Remote Sensing, Vol. 48, No. 4, 2109-2118, Apr. 2010.

[7.] Lee, M. S., "Signal modeling and analysis of a planar phased-array FMCW radar with antenna switching," IEEE Antennas and Wireless Propagation Letters, Vol. 10, 179-182, 2011.

[8.] De Wit, J. J. M., A. Meta, and P. Hoogeboom, "Modified rangedoppler processing for FM-CW synthetic aperture radar," IEEE Geoscience and Remote Sensing Letters, Vol. 3, No. 1, 83-87,Jul. 2006.

[9.] Mittermayer, J., A. Moreira, and O. Loffeld, "Spotlight SAR data processing using the frequency scaling algorithm," IEEE Trans. on Geoscience and Remote Sensing, Vol. 37, No. 5, 2198-2214, Sep. 1999.

[10.] Jiang, Z. H., K. Huang-Fu, and J. W. Wan, "A chirp transform algorithm for processing squint mode FMCW SAR data," IEEE Geoscience and Remote Sensing Letters, Vol. 4, No. 3, 377-381, Jul. 2007.

[11.] Sheen, D. M., H. D. Collins, T. E. Hall, D. L. McMakin, R. P. Gribble, R. H. Severtsen, J. M. Prince, and L. D. Reid, "Real-time wideband holographic surveillance system," U.S. Patent 5557283, Sep. 17, 1996.

[12.] Sheen, D. M., D. L. McMakin, H. D. Collins, T. E. Hall, and R. H. Severtsen, "Concealed explosive detection on personnel using a wideband holographic millimeter-wave imaging system," Proceedings of SPIE, Vol. 2755, 503-513, 1996.

[13.] Sheen, D., D. McMakin, and T. Hall, "Near-field three-dimensional radar imaging techniques and applications," Applied Optics, Vol. 49, E83-E93, Jun. 2010.

[14.] Tan, W., W. Hong, Y. Wang, and Y. Wu, "A novel spherical-wave three-dimensional imaging algorithm for microwave cylindrical scanning geometries," Progress In Electromagnetics Research, Vol. 111, 43-70, 2011.

[15.] Subiza, I., E. Gimeno-Nieves, J. M. Lopez-Sanchez, and J. Fortuny-Guasch, "An approach to SAR imaging by means of non-uniform FFTs," IEEE International Geoscience and Remote Sensing Symposium, Vol. 6, 4089-4091, 2003.

[16.] Song, J., Q. H. Liu, P. Torrione, and L. Collins, "Two-dimensional and three-fimensional NUFFT migration method for landmine detection using ground-penetrating radar," IEEE Trans. on Geoscience and Remote Sensing, Vol. 44, No. 6, 1462-1469, Jun. 2006.

[17.] Dutt, A. and V. Rokhlin, "Fast Fourier transforms for nonequispaced data," SIAM Journal on Scientific Computing, Vol. 14, No. 6, 1368-1393, Nov. 1993.

[18.] Liu, Q. H. and N. Nguyen, "An accurate algorithm for nonuniform fast Fourier transforms (NUFFT)," IEEE Microwave and Guided Letters, Vol. 8, No. 1, 18-20, Jan. 1998.

[19.] Nguyen, N. and Q. H. Liu, "The regular Fourier matrices and nonuniform fast Fourier transforms," SIAM Journal on Scientific Computing, Vol. 21, No. 1, 283-293, Sep. 1999.

[20.] Liu, Q. H., X. M. Xu, B. Tian, and Z. Q. Zhang, "Applications of nonuniform fast transform algorithms in numerical solutions of differential and integral equations," IEEE Trans. on Geoscience and Remote Sensing, Vol. 38, No. 4, 1551-1560, Jul. 2000.

[21.] Fessler, J. A. and B. P. Sutton, "Nonuniform fast Fourier transforms using min-max interpolation," IEEE Trans. on Signal Processing, Vol. 51, No. 2, 560-574, Feb. 2003.

[22.] Greengard, L. and J. Y. Lee, "Accelerating the nonuniform fast Fourier transform," SIAM Review, Vol. 46, No. 3, 443-454, Jul. 2004.

[23.] Lee, J. Y. and L. Greengard, "The type 3 nonuniform FFT and its applications," Journal of Computational Physics, Vol. 206, No. 1, 1-5, Jun. 2005.

[24.] Fessler, J. A., "On NUFFT-based gridding for non-cartesian MRI," Journal of Magnetic Resonance, Vol. 188, 191-195, 2007.

[25.] Zhou, X., H. Sun, J. He, and X. Lu, "NUFFT-based iterative reconstruction algorithm for synthetic aperture imaging radiometers," IEEE Geoscience and Remote Sensing Letters, Vol. 6, No. 2, 273-276, Apr. 2009.

[26.] Li, S. Y., H. J. Sun, B. C. Zhu, and R. Liu, "Two-dimensional NUFFT-based algorithm for fast near-field imaging," IEEE Antennas and Wireless Propagation Letters, Vol. 9, 814-817, 2010.

[27.] Carrara, W. G., R. S. Goodman, and R. M. Majewski, Spotlight Synthetic Aperture Radar Signal Processing Algorithms, Artech House, Boston, MA, 1995.

[28.] Cumming, I. G. and F. H. Wong, Digital Processing of Synthetic Aperture Radar Data Algorithms and Implementation, Artech House, Norwood, MA, 2005.

[29.] Guo, D., H. Xu, and J. Li, "Extended wavenumber domain algorithm for highly squinted sliding spotlight SAR data processing," Progress In Electromagnetics Research, Vol. 114, 17-32, 2011.

S.Y. Li*, B.L. Ren, H.J. Sun, W.D. Hu, and X. Lv

School of Information and Electronics, Beijing Institute of Technology, Beijing 100081, China

* Corresponding author: Shiyong Li (lisy_0723@hotmail.com).

The increasing threat of terrorism has made personnel surveillance increasingly important. Conventional security systems have been effective for detecting metal targets. However, some modern threats, such as plastic or ceramic handguns and knives, as well as extremely dangerous liquid explosives, etc., cannot be detected with conventional security technology, such as metal detectors. X-ray systems can present an effective solution to the problem. However, potential health risks due to exposure to ionizing radiation make it less acceptable to public safety.

Thus, airports, the traveling public, and other secure locations have a need for a new personnel surveillance system that is capable of detecting concealed body-worn threats that are undetectable with conventional security technology. One of the choices is to use millimeter waves which are capable of penetrating common clothing barriers to form an image of a person as well as any concealed items with different reflectivity or emissivity [1-3]. Millimeter waves are nonionizing and, therefore, pose no known health hazard at moderate power levels.

The combination of frequency-modulated continuous-wave (FMCW) and MMW imaging techniques leads to compact, light-weight, cost-effective, low-power operation, and high-resolution imaging systems, which are especially suitable for security and detection application. FMCW signal has been widely used in imaging radar [4-7]. Typically, FMCW systems use long pulses. Therefore, conventional SAR imaging algorithms that use the stop-and-go approximation need to be modified for FMCW image processing. Several conventional algorithms, such as range-Doppler algorithm, frequency-scaling algorithm, and chirp-scaling algorithm, etc., have been modified to focus the FMCW SAR data [8-10]. Especially, an accurate received two-dimensional (2-D) signal model was proposed in [6], in which the effect of the variation of the instantaneous slant range on the transmitted and received signal was accurately represented, and a range-azimuth coupling term was formulated for the first time in the FMCW SAR community.

The 3-D MMW imaging techniques for concealed weapon detection were discussed in [1,11-14], etc., whereas the effects of the continuous motion of the array on the transmitted and echoed signal were not addressed. In this paper, we begin with an analysis of the 3D signal model which accurately includes the effects of the continuous motion of the array on the transmitted and received signal during the pulse time. The main effects of the motion on the image result turn out to be a range walk and out of focus. These effects can be corrected by a matched filter multiplication in the 3-D frequency domain.

After the corrections, the Stolt mapping can be implemented to compensate for the range curvature of all scatterers by an appropriate warping of the wavenumber domain backscattered data. For the 3D imaging, the Stolt mapping is implemented by a 1-D interpolation process in the range-frequency domain, which is repeated for every azimuth-frequency and elevation-frequency point. The complex image can then be reconstructed by the 3-D inverse FFT (IFFT). However, the repeated interpolation process reduces the computational efficiency of the algorithm. In [15,16], the Stolt interpolation and IFFT were substituted by a nonuniform fast Fourier transform (NUFFT). The NUFFT approaches have been developed to overcome the limitation of equally spaced sampling needed by FFT [17-23]. NUFFTs have wide applications in magnetic resonance image reconstruction [24], ground penetrating radar with the range migration algorithm [16], synthetic aperture imaging radiometers [25], and near-field imaging [26]. In this paper, the 3-D NUFFT is used to replace the Stolt interpolation and the followed 3-D IFFT to improve the computational efficiency.

This paper is organized as follows. In Section 2, an accurate backscattered signal model is derived with the corresponding spectrum analysis. The 3-D imaging procedure is given in Section 3. Section 4 shows the simulation results to verify the effectiveness of the motion compensation. The comparison of the computational time between the interpolation-FFT and NUFFT are also shown. Section 5 summarizes the conclusions.

2. FMCW THREE-DIMENSIONAL SIGNAL MODEL AND SPECTRUM ANALYSIS

This section derives an analytical model of the FMCW backscattered signal in the 3-D spatial-frequency domain. We consider the 3-D imaging geometry, as shown in Figure 1. The 1-D antenna array aligns with the x direction, and the array is scanned along the minus y direction.

For an FMCW transceiver, the transmitted signal can be expressed as [27]:

[S.sub.T] (t) = exp [j2[pi] ([f.sub.0]t + 1/2 K[t.sup.2])], (1)

where [f.sub.0] is the center frequency, t is the time variable varying within one cycle of signal transmitting, and K is the frequency sweep rate of the transmitted signal.

The signal is transmitted at an arbitrary time [rho]. And the corresponding instantaneous range between the antenna element and the target is R ([rho]). Let the time [[rho].sub.d] be the round-trip delay time. Then the backscattered signal is received at time t + Td with the corresponding instantaneous range R([rho] + [[rho].sub.d]). Thus, the round-trip delay time can be given by

[[rho].sub.d] R ([rho]) + R([rho] + [[rho].sub.d])/c (2)

where

R([rho]) = [square root of [(x'-x).sup.2] + [(y'-y).sup.2] + [(z'-z).sup.2]] = [square root of [(x'-x).sub.2] + ([y.sub.0]-[[v.sub.[rho]]-y).sup.2] + [(z'- z).sup.2]]0, R([rho] + [[rho].sub.d]) = [square root of [(x'-x).sup.2] + [(y'-y).sup.2] + [(z'-z).sup.2]] = [square root of [(x'-x).sup.2] + [[[y'.sub.0]-v ([rho] + [[rho].sub.d])- y].sup.2] + [(z'-z).sup.2],

and c is the speed of light. (x, y, z) is the location of a point target, and (x', y', z') is the position of the antenna element, both within the same coordinates, as shown in Figure 1. And y' = [y'.sub.0]-v[rho], where [y'.sub.0] is the initial position of the array, v is the scanning velocity, and t = [nT.sub.x] + [mT.sub.y] + t is the continuous time; [T.sub.x] is the period of signal transmitting along the array elements, assuming equal to the signal duration time of one cycle, thus t [member of] [0,[T.sub.x]), and [T.sub.y] is the period of sampling time along y direction; n = 0, 1, ..., N-1, and m = 0, 1,..., M-1, where N is the number of the array elements (corresponding to the number of samples along x direction), and M is the number of samples along y direction.

[FIGURE 1 OMITTED]

Due to the fact that the antenna array is near the target for the personnel surveillance system, (2) can be accurately approximated as

[[rho].sub.D] [approximately equal to] 2R([rho])/c

Omitting the time scaling influences on the envelope, the backscattered signal from a point target is given by

[S.sub.R] (x',y'; t) = [sigma] (x,y,z) [S.sub.T] (t-[[rho].sub.d]) (4)

where [sigma] (x, y, z) is the backscattering coefficient of the point target.

The dechirp-on-receive technology is generally used in the FMCW systems [27]. The received signal and a delayed version of the transmitted signal are mixed in order to reduce the required sampling rate. The intermediate frequency signal after mixing can be written as

[S.sub.IF] (x',y'; t) = [sigma](x, y, z)exp[-j2[pi][f.sub.0]([[rho].sub.d]- [[rho].sub.c])]exp[-j2[pi]K([[rho].sub.d]- [[rho].sub.c])([rho]-[[rho].sub.c])] exp [-j[pi]K ([[rho].sub.d]-[([[rho].sub.c]).sup.2] (5)

where [[rho].sub.c] is the delayed time of the transmitted signal. The last exponential term of (5) is known as the residual video phase (RVP), which is caused by the difference between any two displaced linear FM waveforms [27]. The compensation for RVP is conducted in the frequency domain, which is not the emphasis of this paper. Assuming the RVP has been removed in the following, i.e.,

[S.sub.IF] (x', y'; t) = [sigma] (x, y, z) exp [-j 2[pi][f.sub.0] ([[rho].sub.d]-[[rho].sub.c])] exp [-j2[pi]K ([[rho].sub.d]-[[rho].sub.c]) ([rho]-[t.sub.c])] (6)

Using the substitution/= K(t-[[rho].sub.c]), (6) can be rewritten as

[S.sub.IF] (x', y';f) = [sigma] (x, y, z) exp [-j2[pi] (f + [f.sub.0]) ([[rho].sub.d]-[[rho].sub.c])] (7)

Substituting t = [mT.sub.y] + [nT.sub.x] + t = [[rho].sub.m] + [[rho].sub.n] + t into (7) yields

[S.sub.IF] (x', y', [??]; [[rho].sub.m] [[rho].sub.n]; f) = [sigma]([??])exp[- j4[pi](f + [f.sub.0]) (T([rho])/c- [R.sub.c]/c = [sigma] ([??]) exp[-j4[pi](f + [f.sub.0])(R([[rho].sub.m] + [[rho].sub.n] + t/c-[R.sub.c]/c] (8)

where [??] = (x, y, z) and [R.sub.c] = [ct.sub.c]/2.

Then the 2-D Fourier transform of (8) with respect to the spatial variables x' and [y'.sub.m] ([y'.sub.m] = [v[rho].sub.m] = vmTy) is given by

[S.sub.IF] ([k'.sub.x],[k'.sub.y], f;[[rho].sub.n])

= 1/v [integral] [S.sub.IF](x', y', [??]; [[rho].sub.m], [[rho].sub.n];f) exp[- [jk.sub.x']x') exp(- [jk.sub.y'][y'.sub.m]) dx'[dy'.sub.m]

= 1/v[sigma]([??])[integral] exp[-j[PHI]([k.sub.x'], [k.sub.y'], f; x', [y.sub.m'], [[rho].sub.n])]dx'[d'.sub.ym] (9)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and [k.sub.r] = 2[pi](f + [f.sub.0])/c.

The integral in (9) can be solved by means of the principle of the stationary phase method [27, 28]. At the point of stationary phase, the first partial derivatives of the phase [PHI] ([k.sub.x]',[k.sub.y]', f;x',[y'.sub.m], [[rho].sub.n]) are equal to zeros, such as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (10)

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (11)

Solving (10) and (11) for the stationary phase points, i.e., [y'.sub.m0] and [x.sub.0]', respectively, yields

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (12)

and

[x.sub.0]' = x-[k.sub.x]'(z-z')/[square root of [4k.sup.2.sub.r]- [k.sup.2.sub.x']-[k.sup.2.sub. y']. (13)

Then the 3-D spectrum of the backscattered signal can be expressed as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (14) where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (15)

Based on the relation between frequency and time of FMCW signal, we have t = f/K + [2R.sub.c]/c. Consequently, the phase of (14) can also be expressed as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (15)

where [R.sub.0] =-z'.

Note that the first three terms of (15), i.e., [k.sub.x']x, [k.sub.y']y and [square root of [4k.sup.2.sub.r]-[k.sup.2.sub.x']-[k.sup.2.sub.y']z are linearly dependent on the target position x, y and z, respectively. Thus the image can be reconstructed by using inverse Fourier transform. The fourth term [k.sub.y'][y.sub.0]' is introduced by the initial position of the array. The term [k.sub.y']v[[rho].sub.n] represents the phase variation corresponding to the elevation-range walk caused by the motion of the nth antenna element. [k.sub.y']vf/K is a space-invariant term also caused by the motion of the array within the signal duration time. The terms [k.sub.y']v([2R.sub.c]/c) and [2k.sub.r][R.sub.c] refer to the constant azimuth and range shifts, respectively, and are introduced by the dechirp-on-receive approach. The last term [square root of [4k.sup.2.sub.r]-[k.sup.2.sub.x']- [k.sup.2.sub.y'] [R.sub.0] represents the constant range shift.

3. THREE-DIMENSIONAL IMAGING PROCEDURE

Based on the spectrum analysis of the 3-D FMCW backscattered signal in Section 2, we present a wavenumber domain algorithm for the image reconstruction. The wavenumber domain algorithm (or called range migration algorithm) [27-29] is an accurate algorithm to process the near-field data.

The processing steps of the algorithm are shown in Figure 2. Removing RVP is not the emphasis of this paper. So it is not addressed. As can be seen from Figure 2, the phase corrections are the key parts for the image reconstruction. To remove the unwanted phase in (15), a 3-D matched filter is given by

[G.sub.F] ([k.sub.x'],[k.sub.y'],[k.sub.r];[[rho].sub.n], [R.sub.ref]) = [G.sub.F1]([k.sub.y'],[k.sub.r]; [[rho].sub.n]) x [G.sub.F2] {[k.sub.x'], [k.sub.y'],[k.sub.r]; [R.sub.ref]), (16)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and

[G.sub.F2] ([k.sub.x'],[k.sub.y'], [k.sub.r]; [R.sub.ref]) = exp[j([square root of [4k.sup.2.sub.r]-[4k.sup.2.sub.y']- [k.sub.r]; [R.sub.ref]]],

where [R.sub.ref] is the reference range for focus processing, which is generally defined as the range from the scene center to the aperture center.

Note that the term [G.sub.F1] must be first multiplied exactly after the elevation Fourier transform and before the azimuth Fourier transform, as shown in Figure 2, due to the fact that the term [k.sub.y']V[[rho].sub.n] exists in the elevation-frequency domain and the azimuth-time domain. Then the term [G.sub.F2] is multiplied to correct the range curvature of all scatterers having minimum range [R.sub.ref]. Residual range curvature is still present on those scatterers with a minimum range that is different with [R.sub.ref] [27]. The residual range curvature can be removed by Stolt mapping, which will be shown in the next.

After matched filtering, the remaining signal becomes

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (17)

where

[[PHI].sub.F] ([k.sub.x'],[k.sub.y'], [k.sub.r]; [R.sub.ref]) = [k.sub.x']x- [k.sub.y']y + [square root of [4k.sup.2.sub.r]-[k.sup.2.sub.x']- [k.sup.2.sub.y']z] + [square root of [4k.sup.2.sub.r]-[k.sup.2.sub.x']- [k.sup.2.sub.x']]([R.sub.0]-[R.sub.ref]). (18)

The minus in "-[k.sub.y']y " is introduced by the scanning direction of the array. The distinction between the primed and unprimed coordinate systems can now be dropped since the coordinate systems coincide.

[FIGURE 2 OMITTED]

As can be seen from the phase [[PHI].sub.F], (17) suffices for the image reconstruction if the data are defined continuously in [k.sub.x], [k.sub.y], and [k.sub.r]. However, for practical imaging configurations, the data are discretely sampled at uniform intervals of position and frequency. After 2-D Fourier transform with respect to x and y, and matched filtering, a sampled version of [S.sub.F] ([k.sub.x], [k.sub.y], [k.sub.r]; [R.sub.ref]) is obtained, which are uniformly spaced in [k.sub.x], [k.sub.y] and [k.sub.r]. However, the samples are nonuniformly spaced in [square root of [4k.sup.2.sub.r]-[k.sup.2.sub.x]- [k.sup.2.sub.y]]. In order to reconstruct the image via 3-D IFFT in (17), the data [S.sub.F] will need to be resampled to uniformly spaced positions in [4k.sup.2.sub.r]-[k.sup.2.sub.x]-[k.sup.2.sub.y]. This is accomplished by Stolt interpolation, which converts [S.sub.F] into a 3-D linear phase grating for targets at all ranges [27]. The appropriate change of variables is [square root of [4k.sup.2.sub.r]-[k.sup.2.sub.x]-[k.sup.2.sub.y]] [right arrow] [k.sub.z] (19)

where [k.sub.z]s are uniformly spaced.

After Stolt interpolation, (18) becomes

[PHI]F ([k.sub.x], [k.sub.y], [k.sub.z]; [R.sub.ref]) = [k.sub.x]x-[k.sub.y]y + [k.sub.z]z + [k.sub.z] (R-[R.sub.ref]) (20)

The aforementioned description is focused on a point target. For a volume target, we have

[S.sub.F] ([k.sub.x],[k.sub.y],[k.sub.z]; [R.sub.ref]) = [integral] (1/v) f (x, y, z) exp [-j[k.sub.z] ([R.sub.0]- [R.sub.ref])]

exp [-j ([k.sub.x]x-[k.sub.y]y + [k.sub.z]z)] dxdydz (21)

where f (x, y, z) is the backscattering coefficient of a volume target.

Clearly, a 3-D IFFT can compress the signal in range, azimuthrange and elevation-range, i.e.,

f (x,y,z) = v [integral] [S.sub.F] ([k.sub.x], [k.sub.y], [k.sub.z] ; [R.sub.ref]) exp [j[k.sub.z] ([R.sub.0]- [R.sub.ref])] exp [j ([k.sub.x]x-[k.sub.y]y + [k.sub.z]z)] [dk.sub.x] [dk.sub.y] [dk.sub.z] (22)

The constant factor relative to 2n is ignored in (22).

However, the change of variables and the interpolation procedures will introduce approximation errors [16] and are time-consuming. In this paper, we substitute the Stolt interpolation and IFFT with the 3D inverse NUFFT to reconstruct the image. The NUFFT idea was first presented by Dutt and Rokhlin [17]; later a new approach was proposed by Liu and Nguyen, using the regular Fourier matrices [18, 19]. A minmax method for fast and accurate approximation of the NUFFT was developed by Fessler and Sutton [21], which is optimal in the min-max sense of minimizing the worst-case approximation error over all signals of unit norm.

Here, the NUFFT algorithm developed by Fessler and Sutton is implemented. For the sake of simplicity, we consider only the 1-D case:

[x.sub.n] = [M-1.summation over (m-0)] [y.sub.m] exp (j[[omega].sub.m]n, n = 0,1. (23)

where [y.sub.m] represents the nonuniform frequency domain sample, [[omega].sub.m] is the nonuniform digital frequency in radian, and [x.sub.n] is the uniform output signal in time domain.

The computation of (23) can be summarized as follows:

1) The "gridding" step: [X.sub.k] = [m-1.summation over (m=0)] [v.sub.mk][y.sub.m], where [V.sub.mk]s denote interpolation coefficients. This step requires O (JM) operations, where J is the number of the nearest neighbors to wm.

2) K-point inverse discrete Fourier transform (DFT) of [X.sub.k] scaled

by K: [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], 0 [less than or equal to] n [less than or equal to] N-1. The output

signal [[??].sub.n]s are the first N signal values. This step requires O (K [log.sub.2] N) operations by using the reduced FFT.

3) Scaling each [[??].sub.n] by [s.sup.*.sub.n] to get [x.sub.n], where the sns are the scaling factors, and "*" denotes conjugate. The overall operation count per NUFFT is O(JM) + O(K[log.sup.2] N). Typically, we choose K = 2N and J [less than or equal to] 10.

Therefore, the overall computational requirements are similar to an FFT but with a larger constant. For more details, the readers are referred to [21]. The computational complexity of the 3-D NUFFT is known as O ([N.sup.3] [log.sub.2] N).

The flow diagram of the image reconstruction process is demonstrated in Figure 2 in which the 3-D inverse NUFFT can also be replaced by 1-D inverse NUFFT and 2-D IFFT due to the fact that [k.sub.x]s and [k.sub.y]s are uniformly spaced.

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

4. SIMULATION RESULTS

In this section, simulations are performed to validate the performance of the proposed method. The simulation parameters are listed in Table 1.

For a first validation of the derivation, the frequency domain data in (14) is multiplied by the conventional and the proposed matched filters, respectively, for a target in the scene center.

Under the stop-and-go approximation, the conventional matched filter [27] is given by

[G.sub.F] ([k.sub.x],[k.sub.y],[k.sub.r]; [R.sub.ref]) = exp [j([square root of [4k.sup.2.sub.r],-[k.sup.2.sub.x]- [k.sup.2.sub.y], [R.sub.ref] + [k.sub.y] [y.sub.0]'-[2k.sub.r][R.sub.c], (24)

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

Figure 3 shows the phase of the signal after matched filtering with a reference function that makes the stop-and-go approximation, i.e., the phase of the remaining signal after (14) is multiplied by (24). The phase of the remaining signal in the elevation-frequency and range-frequency domain, i.e., in the [k.sub.y]-[k.sub.r] domain, is shown in Figure 3(a), and the phase in the elevation-frequency and azimuth-frequency domain, i.e., in the [k.sub.y]-[k.sub.x] domain, is illustrated in Figure 3(b). The phase in the azimuth-frequency and range-frequency domain is not demonstrated due to the fact that it is not affected by the motion of the array, as shown in Figure 1.

Figure 4 illustrates the phase results produced by filtering the signal (14) with the proposed function in (16), i.e., the phase of (17).

[FIGURE 8 OMITTED]

If the signal (14) is completely compensated for by the matched filter, the phase of the remaining signal, i.e., the phase of (17), should be constant and does not vary with [k.sub.x], [k.sub.y] and [k.sub.r] in the case that the simulated point target is in the scene center, i.e., (x, y, z) = (0, 0, 0) and [R.sub.0] = [R.sub.ref]. Clearly, only in Figure 4 does the phase in the support band is practically constant, indicating that the matched filter does match the received signal.

Figures 5(a) to (c) show the slices of the reconstructed 3-D image of a point target located in the scene center. These slices correspond to the x-y (with z = 0), y-z (with x = 0), and x-z (with y = 0) section planes, respectively, with the matched filtering function that makes the stop-and-go approximation. Clearly, the motion of the array has a significant effect on the azimuth-elevation (the x-y plane) image, as shown in Figure 5(a). It can be seen that the nulls of the image are raised which make the slice image, i.e., the 2-D impulse response (or called point spread function) departure from the ideal 2-D sinc function. In addition, a range walk exists in the elevation direction which is caused by the motion, as shown in Figures 5(a) and (b). The image in the x-z plane is not affected by the motion, as shown in Figure 5(c).

[FIGURE 9 OMITTED]

Figures 6(a) to (c) demonstrate the corresponding image results produced by the proposed matched filtering function. From the comparison of Figures 5 and 6, the improvements due to the proposed correction are visible. The dynamic range of the images is 30 dB.

The next simulation is intended to consider a cubic target consisting of 8 scatterers distributed at the vertices of the cube as shown in Figure 7. All scatterers have the same radar cross section (RCS): 1 m2. Figure 8 shows the projection of the images onto the three main planes, with the matched filtering function that makes the stop-and-go approximation. The image results produced by the proposed matched filtering function are shown in Figure 9. The intersections of the dashed lines represent the real positions of the scatterers. The defocusing and range walk are clearly seen in Figure 8, which are the same as shown in Figure 5.

Figure 10 shows the comparison of the computational time (CPU time) between the Stolt interpolation-based method and the NUFFT-based method. It is evident that the NUFFT-based method is much faster than the Stolt interpolation-based method, especially with the increase of the sample points at the 2-D aperture.

[FIGURE 10 OMITTED]

5. CONCLUSIONS

A wavenumber domain algorithm with motion compensation to process the near-field 3-D imaging has been developed. In this algorithm, the effects of the array motion within the signal duration time have been described and processing solutions were given. The unwanted terms caused by the motion have been removed and the focusing performance is improved. To reduce the computational time of the image reconstruction process, we substitute the Stolt interpolation and FFT with a 3-D NUFFT.

The proposed algorithm has been validated with simulations. The algorithm can also be extended to process the 3-D imaging for cylindrical scanning geometry.

Received 24 November 2011, Accepted 31 December 2011, Scheduled 11 January 2012

REFERENCES

[1.] Sheen, D. M., D. L. McMakin, and T. E. Hall, "Three-dimensional millimeter-wave imaging for concealed weapon detection," IEEE Trans. on Microwave Theory and Techniques, Vol. 49, No. 9, 1581-1592, Sep. 2001.

[2.] Appleby, R. and R. N. Anderton, "Millimeter-wave and submillimeter-wave imaging for security and surveillance," Proceedings of the IEEE, Vol. 95, No. 8, 1683-1690, Aug. 2007.

[3.] Yeom, S., D. Lee, H. Lee, J. Son, and V. P. Gushin, "Distance estimation of concealed objects with stereoscopic passive millimeter-wave imaging," Progress In Electromagnetics Research, Vol. 115, 399-407, 2011.

[4.] Huang, Y., P. V. Brennan, D. Patrick, I. Weller, P. Roberts, and K. Hughes, "FMCW based MIMO imaging radar for maritime navigation," Progress In Electromagnetics Research, Vol. 115, 327-342, 2011.

[5.] Meta, A., P. Hoogeboom, and L. P. Ligthart, "Signal processing for FMCW SAR," IEEE Trans. on Geoscience and Remote Sensing, Vol. 45, No. 11, 3519-3532, Nov. 2007.

[6.] Wang, R., O. Loffeld, H. Nies, S. Knedlik, M. Hagelen, and H. Essen, "Focus FMCW SAR data using the wavenumber domain algorithm," IEEE Trans. on Geoscience and Remote Sensing, Vol. 48, No. 4, 2109-2118, Apr. 2010.

[7.] Lee, M. S., "Signal modeling and analysis of a planar phased-array FMCW radar with antenna switching," IEEE Antennas and Wireless Propagation Letters, Vol. 10, 179-182, 2011.

[8.] De Wit, J. J. M., A. Meta, and P. Hoogeboom, "Modified rangedoppler processing for FM-CW synthetic aperture radar," IEEE Geoscience and Remote Sensing Letters, Vol. 3, No. 1, 83-87,Jul. 2006.

[9.] Mittermayer, J., A. Moreira, and O. Loffeld, "Spotlight SAR data processing using the frequency scaling algorithm," IEEE Trans. on Geoscience and Remote Sensing, Vol. 37, No. 5, 2198-2214, Sep. 1999.

[10.] Jiang, Z. H., K. Huang-Fu, and J. W. Wan, "A chirp transform algorithm for processing squint mode FMCW SAR data," IEEE Geoscience and Remote Sensing Letters, Vol. 4, No. 3, 377-381, Jul. 2007.

[11.] Sheen, D. M., H. D. Collins, T. E. Hall, D. L. McMakin, R. P. Gribble, R. H. Severtsen, J. M. Prince, and L. D. Reid, "Real-time wideband holographic surveillance system," U.S. Patent 5557283, Sep. 17, 1996.

[12.] Sheen, D. M., D. L. McMakin, H. D. Collins, T. E. Hall, and R. H. Severtsen, "Concealed explosive detection on personnel using a wideband holographic millimeter-wave imaging system," Proceedings of SPIE, Vol. 2755, 503-513, 1996.

[13.] Sheen, D., D. McMakin, and T. Hall, "Near-field three-dimensional radar imaging techniques and applications," Applied Optics, Vol. 49, E83-E93, Jun. 2010.

[14.] Tan, W., W. Hong, Y. Wang, and Y. Wu, "A novel spherical-wave three-dimensional imaging algorithm for microwave cylindrical scanning geometries," Progress In Electromagnetics Research, Vol. 111, 43-70, 2011.

[15.] Subiza, I., E. Gimeno-Nieves, J. M. Lopez-Sanchez, and J. Fortuny-Guasch, "An approach to SAR imaging by means of non-uniform FFTs," IEEE International Geoscience and Remote Sensing Symposium, Vol. 6, 4089-4091, 2003.

[16.] Song, J., Q. H. Liu, P. Torrione, and L. Collins, "Two-dimensional and three-fimensional NUFFT migration method for landmine detection using ground-penetrating radar," IEEE Trans. on Geoscience and Remote Sensing, Vol. 44, No. 6, 1462-1469, Jun. 2006.

[17.] Dutt, A. and V. Rokhlin, "Fast Fourier transforms for nonequispaced data," SIAM Journal on Scientific Computing, Vol. 14, No. 6, 1368-1393, Nov. 1993.

[18.] Liu, Q. H. and N. Nguyen, "An accurate algorithm for nonuniform fast Fourier transforms (NUFFT)," IEEE Microwave and Guided Letters, Vol. 8, No. 1, 18-20, Jan. 1998.

[19.] Nguyen, N. and Q. H. Liu, "The regular Fourier matrices and nonuniform fast Fourier transforms," SIAM Journal on Scientific Computing, Vol. 21, No. 1, 283-293, Sep. 1999.

[20.] Liu, Q. H., X. M. Xu, B. Tian, and Z. Q. Zhang, "Applications of nonuniform fast transform algorithms in numerical solutions of differential and integral equations," IEEE Trans. on Geoscience and Remote Sensing, Vol. 38, No. 4, 1551-1560, Jul. 2000.

[21.] Fessler, J. A. and B. P. Sutton, "Nonuniform fast Fourier transforms using min-max interpolation," IEEE Trans. on Signal Processing, Vol. 51, No. 2, 560-574, Feb. 2003.

[22.] Greengard, L. and J. Y. Lee, "Accelerating the nonuniform fast Fourier transform," SIAM Review, Vol. 46, No. 3, 443-454, Jul. 2004.

[23.] Lee, J. Y. and L. Greengard, "The type 3 nonuniform FFT and its applications," Journal of Computational Physics, Vol. 206, No. 1, 1-5, Jun. 2005.

[24.] Fessler, J. A., "On NUFFT-based gridding for non-cartesian MRI," Journal of Magnetic Resonance, Vol. 188, 191-195, 2007.

[25.] Zhou, X., H. Sun, J. He, and X. Lu, "NUFFT-based iterative reconstruction algorithm for synthetic aperture imaging radiometers," IEEE Geoscience and Remote Sensing Letters, Vol. 6, No. 2, 273-276, Apr. 2009.

[26.] Li, S. Y., H. J. Sun, B. C. Zhu, and R. Liu, "Two-dimensional NUFFT-based algorithm for fast near-field imaging," IEEE Antennas and Wireless Propagation Letters, Vol. 9, 814-817, 2010.

[27.] Carrara, W. G., R. S. Goodman, and R. M. Majewski, Spotlight Synthetic Aperture Radar Signal Processing Algorithms, Artech House, Boston, MA, 1995.

[28.] Cumming, I. G. and F. H. Wong, Digital Processing of Synthetic Aperture Radar Data Algorithms and Implementation, Artech House, Norwood, MA, 2005.

[29.] Guo, D., H. Xu, and J. Li, "Extended wavenumber domain algorithm for highly squinted sliding spotlight SAR data processing," Progress In Electromagnetics Research, Vol. 114, 17-32, 2011.

S.Y. Li*, B.L. Ren, H.J. Sun, W.D. Hu, and X. Lv

School of Information and Electronics, Beijing Institute of Technology, Beijing 100081, China

* Corresponding author: Shiyong Li (lisy_0723@hotmail.com).

Table 1. Simulation parameters. Parameter Value Carrier frequency 35 GHz Bandwidth 5 GHz Sampling interval: [DELTA]x' 0.007 m Number of array elements 64 Sampling interval: [DELTA]y' 0.007 m Number of samples along y direction 64 Scanning velocity 1 m/s Referenced range [R.sub.0] 1 m

Printer friendly Cite/link Email Feedback | |

Author: | Li, S.Y.; Ren, B.L.; Sun, H.J.; Hu, W.D.; Lv, X. |
---|---|

Publication: | Progress In Electromagnetics Research |

Article Type: | Report |

Geographic Code: | 9CHIN |

Date: | Feb 1, 2012 |

Words: | 5039 |

Previous Article: | Design of dual-band bandpass filters with controllable bandwidths using new map-ping function. |

Next Article: | A novel wideband antenna array with tightly coupled octagonal ring elements. |

Topics: |