MODIFIED CYLINDRICAL HOLOGRAPHIC ALGO- RITHM FOR THREE-DIMENSIONAL

19 Pages·2012·0.61 MB·English

Progress In Electromagnetics Research, Vol. 128, 519–537, 2012
MODIFIED CYLINDRICAL HOLOGRAPHIC ALGO-
RITHM FOR THREE-DIMENSIONAL MILLIMETER-
WAVE IMAGING
B. L. Ren, S. Y. Li*, H. J. Sun, W. D. Hu, and X. Lv
School of Information and Electronics, Beijing Institute of Technology,
Beijing 100081, China
Abstract—Millimeter-wave (MMW) imaging techniques have been
developedforthedetectionofconcealedweaponsandplasticexplosives
carriedonpersonnelatmajortransportationhubsandsecurelocations.
The combination of frequency-modulated continuous-wave (FMCW)
technologyandMMWimagingtechniquesleadstowideband,compact,
and cost-effective systems which are especially suitable for security
detection. Cylindrical three-dimensional (3-D) imaging technique,
with the ability of viewing multiple sides, is an extension of rectilinear
3-Dimaging techniqueonly viewing a single side. Due to the relatively
long signal sweep time, the conventional stop-and-go approximation of
the pulsed systems is not suitable for FMCW systems. Therefore, a
3-D backscattered signal model including the effects of the continuous
motion within the signal duration time is developed for cylindrical
imagingsystems. Then, aholographicimagereconstructionalgorithm,
with motion compensation, is presented and demonstrated by means
of numerical simulations.
1. INTRODUCTION
Personnel surveillance has been increasingly important due to the
increasing threat of terrorism. Some modern threats, such as plastic
or ceramic handguns and knives, as well as extremely dangerous liquid
explosives, have emerged in recent years, which make the conventional
security systems ineffective. Although some systems, such as X-ray
systems and terahertz systems, can present effective solutions to this
problem,X-raysystems[1]areunacceptableduetothepotentialhealth
risks, and terahertz systems [2] are immature and too expensive for
Received 16 April 2012, Accepted 30 May 2012, Scheduled 8 June 2012
* Correspondingauthor: ShiyongLi(lisy [email protected]).
520 Ren et al.
civilian applications. Thus, we need a personnel surveillance system
withtheabilityofdetectingalltypesofmodernthreatscost-effectively
and with the advantage of no harm to the person under surveillance.
One of the choices is to use millimeter-wave which can penetrate
common clothing barriers to form an image of a person as well as any
concealed items with different reflectivity [3–5]. Millimeter waves are
nonionizing and, therefore, pose no known health hazard at moderate
power levels.
The combination of FMCW and MMW imaging techniques leads
to wideband, compact, cost-effective, low-power operation, and high-
quality imaging systems, which are especially suitable for security
and detection application. FMCW signal has been widely used
in imaging radar [6–9]. The continuous motion of the antenna
while transmitting and receiving the signal for the FMCW systems
is no longer negligible. Therefore, the stop-and-go approximation
which is used in the conventional SAR imaging algorithms needs
to be modified for FMCW image processing. Several conventional
algorithms, such as wavenumber domain algorithm, frequency-scaling
algorithm, range-Doppler algorithm and chirp-scaling algorithm, etc.,
havebeenmodifiedtofocustheFMCWSARdata[8,10–12]. In[8], an
accurate two-dimensional (2-D) received signal model was presented,
in which a range-azimuth coupling term was formulated for the first
time in the FMCW SAR community.
The rectilinear and cylindrical 3-D imaging techniques for
concealed weapon detection were described along with several
application results in [3,13–17]. Wavenumber domain algorithms (or
called range migration algorithms, abbreviated to RMA) [1,13,20]
have been widely used in rectilinear 2-D and 3-D imaging system.
Cylindrical3-Dimagingtechniqueisanextensionofrectilinearimaging
technique. A cylindrical holographic algorithm was first developed
by Soumekh [18,19] and was used for security imaging by Sheen et
al. [15]. A 3-D RMA was developed for cylindrical geometries by
backpropagatingthebackscattereddataontoaplanaraperturein[21].
Cylindrical 3-D imaging technique has the ability of viewing
multiple sides [14], so we choose it for personal security detection.
The radius of the scan should be generally on the order of 1m,
obviously in the near-field of millimeter-wave imaging for personnel
surveillance. Thecylindricalalgorithmsmentionedabovearegenerally
used without considering the effects of the continuous motion within
the signal duration time. In this paper, we analyze the 3-D signal
model which takes the characteristics of FMCW systems into account.
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 compensated for
Progress In Electromagnetics Research, Vol. 128, 2012 521
Figure 1. Near-field 3-D imaging configuration.
by a matched filter multiplication in the 3-D wavenumber domain.
Then, a deconvolution of the Green’s function is performed. After
the corrections, the range curvature of all scatterers having the same
minimum range as the scatterer at the reference range is corrected
exactly . The residual curvature of the scatterers at other locations
is compensated for by implementing the bilinear mapping from the
polar format to rectangular format. The complex image can then be
reconstructed by the 3-D inverse fast Fourier transform (IFFT).
This paper is organized as follows. In Section 2, an accurate
backscattered signal model and the imaging algorithm are presented.
The 3-D imaging procedure is given in Section 3. Section 4
shows the simulation results to verify the effectiveness of the motion
compensation. Section 5 summarizes the conclusions.
2. FMCW THREE-DIMENSIONAL SIGNAL MODEL
AND IMAGING ALGORITHM
This section derives a holographic imaging algorithm based on the
analytical model of the FMCW backscattered signal in the 3-D
wavenumberdomainforcylindricalimagingsystems. The3-Dimaging
geometryisshowninFigure1. Theone-dimensional(1-D)arcantenna
array aligned with the azimuth direction is scanned along the minus y
direction.
The transmitted FM signal can be expressed as [22]:
(cid:183) (cid:181) (cid:182)(cid:184)
1
s (t) = exp j2π f t+ Kt2 , (1)
T 0
2
where f is the carried base frequency, t is the time variable varying
0
within one cycle of signal transmitting, and K is the frequency sweep
522 Ren et al.
rate of the transmitted signal.
Assuming that the signal is transmitted at an arbitrary time τ
and received at time τ +τ . The time τ is the round-trip delay time.
d d
The corresponding instantaneous ranges between the antenna element
and the target are R(τ) and R(τ +τ ), respectively. Consequently,
d
the round-trip delay time can be expressed as
R(τ)+R(τ +τ )
d
τ = , (2)
d
c
where
(cid:113)
R(τ) = (R cosθ(cid:48)−x)2+(y(cid:48)−y)2+(R sinθ(cid:48)−z)2
0 0
(cid:113)
= (R cosθ(cid:48)−x)2+(y(cid:48) −vτ −y)2+(R sinθ(cid:48)−z)2,
0 0 0
(cid:113)
R(τ+τ ) = (R cosθ(cid:48)−x)2+(y(cid:48)−y)2+(R sinθ(cid:48)−z)2
d 0 0
(cid:113)
= (R cosθ(cid:48)−x)2+[y(cid:48)−v(τ+τ )−y]2+(R sinθ(cid:48)−z)2,
0 0 d 0
y(cid:48) = y(cid:48) −vτ,
0
τ=nT +mT v+t, n=0,1,...N−1, m=0,1,...M−1, t∈(0,T )
θ(cid:48) y θ(cid:48)
and c is the speed of light, y(cid:48) represents the starting scanning position
0
of the arc antenna array, R is the radius of the arc antenna array, v
0
is the scanning velocity and τ is the continuous time. (x, y, z) is the
location of a point target in a 3-D Cartesian coordinate system. And
(R cosθ(cid:48), y(cid:48), R sinθ(cid:48))isthepositionoftheantennaelement, wherey(cid:48)
0 0
represents the scanning position of the arc antenna array, θ(cid:48) represents
the angle in the cylindrical coordinate system as shown in Figure 1 .
The origins of the two coordinate systems are the same. The period of
signaltransmittingalongthearcarrayelementsintheazimuthdomain
is T which is assumed to be equal to the pulse repetition interval. T
θ(cid:48) y
is the period of sampling time along y direction. The number of the
array elements is N corresponding to the number of samples along
the azimuth direction. M is the number of samples along elevation
direction.
Owing to the fact that the distance between the antenna array
and the target is very short for the personnel security detection, (2) is
accurately approximated as
2R(τ)
τ ≈ . (3)
d
c
Neglectingthetimescalinginfluencesontheenvelope,thereceived
signal from a single point target can be expressed as
(cid:161) (cid:162)
s θ(cid:48),y(cid:48);t = σ(x,y,z)s (t−τ ), (4)
R T d
Progress In Electromagnetics Research, Vol. 128, 2012 523
where σ(x, y, z) is the reflectivity of the single point target.
In dechirp-on-receive SAR systems, the received signal is mixed
with a reference signal in order to reduce the sampling requirements
anddatarates[22]. Thereferencesignalisaversionofthetransmitted
signal with delayed time τ . The dechirped signal can be written as
c
(cid:161) (cid:162)
s θ(cid:48),y(cid:48);t = σ(x,y,z)exp[−j2πf (τ −τ )]
IF 0 d c
(cid:104) (cid:105)
exp[−j2πK(τ −τ )(t−τ )]exp jπK(τ −τ )2 .(5)
d c c d c
The last exponential term of (5) is known as the residual video phase
(RVP), which is introduced by the operation of the dechirp-on-receive
technology [22]. The compensation for RVP is conducted in the
frequency domain, which has been introduced in [23]. Assuming the
RVP has been removed in the remainder of the derivation, i.e.,
(cid:161) (cid:162)
s θ(cid:48),y(cid:48);t
IF
= σ(x,y,z)exp[−j2πf (τ −τ )]exp[−j2πK(τ −τ )(t−τ )]. (6)
0 d c d c c
Substituting f = K(t−τ ) into (6), we obtain
c
(cid:161) (cid:162)
s θ(cid:48),y(cid:48);f = σ(x,y,z)exp[−j2π(f +f )(τ −τ )]. (7)
IF 0 d c
Applying the substitution of τ = mT +nT +t = τ +τ +t to (7):
y θ(cid:48) m n
(cid:161) (cid:162)
s θ(cid:48),f,y(cid:48);τ ,τ ,t
IF m n
(cid:183) (cid:181) (cid:182)(cid:184)
R(τ +τ +t) R
m n c
= σ(x,y,z)exp −j4π(f +f ) − , (8)
0
c c
where R = cτ /2.
c c
In order to obtain the reflectivity σ(x, y, z),we should first give
the Fourier analysis of the function s (θ(cid:48), f,y(cid:48); τ , τ , t) in the
IF m n
elevation frequency domain.
So performing 1-D Fourier transform to (8) with respect to the
spatial variables y(cid:48) (y(cid:48) = vτ = vmT ), yields
m m m y
(cid:90)
(cid:161) (cid:162) 1 (cid:161) (cid:162) (cid:161) (cid:162)
S θ(cid:48),f,k(cid:48);τ ,t = s θ(cid:48),f;y(cid:48) ,τ ,t exp −jk y(cid:48) dy(cid:48)
IF y n v IF m n y(cid:48) m m
(cid:90)
1 (cid:163) (cid:161) (cid:162)(cid:164)
= σ(x,y,z) exp−jΦ θ(cid:48),k ,k ;y(cid:48) ,τ ,t dy(cid:48) ,(9)
v r y(cid:48) m n m
where
(cid:161) (cid:162) 4π(f +f )
Φ θ(cid:48),k ,k ; y(cid:48) ,τ ,t = 0 [R(τ +τ +t)−R ]+k y(cid:48)
r y(cid:48) m n c m n c y(cid:48) m
= 2k [R(τ +τ +t)−R ]+k y(cid:48)
r m n c y(cid:48) m
and k = 2π(f +f )/c.
r 0
524 Ren et al.
The above integral can be evaluated by applying the principle of
the stationary phase method [19,22,24]. At the point of stationary
phase, the phase Φ(θ(cid:48), k , k ; y(cid:48) , τ , t) takes an extreme value, such
r y(cid:48) m n
as
(cid:161) (cid:162)(cid:175)
∂Φ θ(cid:48),kr,ky(cid:48);ym(cid:48) ,τn,t (cid:175)(cid:175)
(cid:175) = 0, (10)
∂y(cid:48) (cid:175)
m
y(cid:48) =y(cid:48)
(cid:161) (cid:162)(cid:175) m m0
∂2Φ θ(cid:48),kr,ky(cid:48);ym(cid:48) ,τn,t (cid:175)(cid:175)
(cid:175) (cid:54)= 0. (11)
∂2y(cid:48) (cid:175)
m
y(cid:48) =y(cid:48)
m m0
Solving (10) for the stationary phase point y(cid:48) , yields
m0
(cid:163) (cid:164)
(4k2−k2 ) y(cid:48) +v(τ +t)−(y(cid:48) −y) 2
r y(cid:48) m0 n 0
(cid:104) (cid:105)
(cid:161) (cid:162) (cid:161) (cid:162)
= k2 R cosθ(cid:48)−x 2+ R sinθ(cid:48)−z 2 .
y(cid:48) 0 0
After derivation, we can obtain
k R
y(cid:48) = −y+y(cid:48) −vτ − (cid:113) y(cid:48) xz −vt
m0 0 n 4k2−k2
r y(cid:48)
k R
= −y− (cid:113) y(cid:48) xz +y(cid:48) −vτ −vt, (12)
0 n
4k2−k2
r y(cid:48)
where (cid:112)
R = (R cosθ(cid:48)−x)2+(R sinθ(cid:48)−z)2. (13)
xz 0 0
If the point target is positioned at the center of the coordinates
and y(cid:48) ∈ (−L, L), then
(cid:195) (cid:33)
L L
k ∈ −2k (cid:112) ,2k (cid:112)
y(cid:48) r r
R2+L2 R2+L2
0 0
and 2L is the size of the elevation synthetic aperture [18].
Then (9) can be expressed as
(cid:161) (cid:162) 1 (cid:169) (cid:163) (cid:161) (cid:162)
S θ(cid:48),k ,k ;τ ,t = σ(x,y,z)exp −j −k y+k y(cid:48)−vτ −vt
IF r y(cid:48) n v y(cid:48) y(cid:48) 0 n
(cid:113) (cid:105)(cid:111)
−2k R +R 4k2−k2 . (14)
r c xz r y(cid:48)
Consequently, by the substitution of t = f/K +2R /c, (14) can
c
also be rewritten as
(cid:189) (cid:183)
(cid:161) (cid:162) 1 f
S θ(cid:48),k ,k ;τ = σ(x,y,z)exp −j −k y+k y(cid:48)−k vτ −k v
IF r y(cid:48) n y(cid:48) y(cid:48) y(cid:48) n y(cid:48)
v 0 K
(cid:184)(cid:190)
(cid:113)
2R
−k v c −2k R +R 4k2−k2 . (15)
y(cid:48) c r c xz r y(cid:48)
Progress In Electromagnetics Research, Vol. 128, 2012 525
The minus in “−k y” is introduced by the scanning direction of the
y(cid:48)
array.
For a volume target we can write
(cid:90) (cid:90) (cid:90) (cid:189) (cid:183)
(cid:161) (cid:162) 1
S k ,k ;θ(cid:48) = σ(x,y,z)×exp −j −k y+k y(cid:48)−k vτ
IF y(cid:48) r y(cid:48) y(cid:48) y(cid:48) n
v 0
x y z
(cid:184)(cid:190)
(cid:113)
f 2R
−k v −k v c−2k R +R 4k2−k2 dxdydz
y(cid:48) K y(cid:48) c r c xz r y(cid:48)
(cid:189) (cid:183) (cid:184)(cid:190) (cid:90)(cid:90)(cid:90)
f 2R 1
= exp −j k y(cid:48)−k vτ −k v −k v c−2k R × σ(x,y,z)
y(cid:48) y(cid:48) n y(cid:48) y(cid:48) r c
0 K c v
xy z
(cid:112)
exp(jk(cid:48)y)exp[−jk (R cosθ(cid:48)−x)2+(R sinθ(cid:48)−y)2]dxdydz, (16)
y xz 0 0
(cid:113)
where k = 4k2−k2 .
xz r y(cid:48)
The above model can be given by
(cid:189) (cid:183) (cid:184)(cid:190)
(cid:161) (cid:162) 1 f 2R
S k ,k ;θ(cid:48) = exp −j k y(cid:48)−k vτ −k v −k v c−2k R
IF y(cid:48) r y(cid:48) y(cid:48) n y(cid:48) y(cid:48) r c
v 0 K c
(cid:90) (cid:90)
× F (x,z,k )g∗(k ,x,z)dxdz, (17)
y y(cid:48) θ(cid:48) xz
x z
where (cid:90)
F (x,z,k ) = σ(x,y,z)exp(jk y)dy
y y(cid:48) y(cid:48)
y
and
(cid:104) (cid:112) (cid:105)
g (k ,x,z) = exp jk (R cosθ(cid:48)−x)2+(R sinθ(cid:48)−z)2
θ(cid:48) xz xz 0 0
.
Based on the generalized Parseval’s theorem [18,19], this model can
be expressed as
(cid:189) (cid:183) (cid:184)(cid:190)
(cid:161) (cid:162) 1 f 2R
S k ,k ;θ(cid:48) = exp −j k y(cid:48)−k vτ −k v −k v c−2k R
IF y(cid:48) r v y(cid:48) 0 y(cid:48) n y(cid:48) K y(cid:48) c r c
(cid:90) (cid:90)
× F k ,k ,k )G∗ (k ,k ,k )dk dk . (18)
( x z y(cid:48) θ(cid:48) xz x z x z
kx kz
Transforming the above double integral into the polar coordinate in
526 Ren et al.
the wavenumber domain, one can obtain
(cid:189) (cid:183) (cid:184)(cid:190)
(cid:161) (cid:162) 1 f 2R
S k ,k ;θ(cid:48) = exp −j k y(cid:48)−k vτ −k v −k v c−2k R
IF y(cid:48) r v y(cid:48) 0 y(cid:48) n y(cid:48) K y(cid:48) c r c
 
(cid:90) (cid:90)
 
× kxz Fkxz(kxz,φ,ky(cid:48))G∗θ(cid:48)kxz(kxz,θ(cid:48)−φ)dφdkxz.(19)
kxz φ
Using the Fourier properties of circular symmetric function [18,19], we
have
G (k ,φ) = exp[jk R cos(φ)].
θ(cid:48)kxz xz xz 0
Theterminsidethe[]of(19)representsaconvolutionintheθ(cid:48) domain,
such as
(cid:90)
(cid:161) (cid:162) (cid:161) (cid:162) (cid:161) (cid:162) (cid:161) (cid:162)
F k ,φ,k G∗ k ,θ(cid:48)−φ dφ =F k ,θ(cid:48),k ⊗G∗ k ,θ(cid:48) .
kxz xz y(cid:48) θ(cid:48)kxz xz kxz xz y(cid:48) kxz xz
φ
Based on the derivation in [18,19], yields
(cid:161) (cid:162)
F k ,θ(cid:48),k
kxz xz y(cid:40)(cid:48) (cid:169) (cid:161) (cid:162) (cid:163) (cid:161) (cid:162)(cid:164)(cid:170)(cid:41)
FFT S k ,k ;θ(cid:48) exp jΦ k ,k ;f
(θ(cid:48)) IF y(cid:48) r y(cid:48) r
= v∗IFFT (cid:163) (cid:164) , (20)
(θ(cid:48)) FFT G∗ (k ,θ(cid:48))
(θ(cid:48)) kxz xz
where
(cid:161) (cid:162) (cid:161) (cid:162) (cid:163) (cid:161) (cid:162)(cid:164)
S(cid:48) k ,k ;θ(cid:48) = S k ,k ;θ(cid:48) exp jΦ k ,k ;f
IF y(cid:48) r IF y(cid:48) r y(cid:48) r
and
(cid:161) (cid:162) f 2R
Φ k ,k ;f = k y(cid:48) −k vτ −k v −k v c −2k R . (21)
y(cid:48) r y(cid:48) y(cid:48) n y(cid:48) y(cid:48) r c
0 K c
Note that the first term k y(cid:48) in (21) is introduced by the initial
y(cid:48)
0
position of the array. The term k vτ represents the phase variation
y(cid:48) n
corresponding to the elevation-range walk caused by the motion of the
nthantennaelement. k vf/K isaspace-invarianttermalsocausedby
y(cid:48)
the motion of the array within the sweep time. The termsk v(2R /c)
y(cid:48) c
and 2k R refer to the constant shifts of the azimuth and range,
r c
respectively, and are introduced by the dechirp-on-receive technology.
After the reference function multiplication (RFM), we can obtain
F (k , θ(cid:48), k ) as shown in (20). The process of RFM is discussed
kxz xz y(cid:48)
in detail in Section 3.
The cylindrical samples of F (k , θ(cid:48), k ) are then converted to
kxz xz y(cid:48)
the samples of F (k ,k ,k )via an interpolation algorithm, where
x z y(cid:48)
k = k cosθ,
x xz
k = k sinθ,
z xz
k = k .
y y(cid:48)
Progress In Electromagnetics Research, Vol. 128, 2012 527
Figure 2. Flow diagram of the proposed image reconstruction
method.
The distinction between the primed and unprimed coordinate
systems can now be dropped since the coordinate systems coincide.
Since σ(x, y, z) and F(k , k , k ) form a Fourier transform pair, the
x y z
3-D inverse Fourier transform of this signal is the desired image, i.e.,
(cid:90)(cid:90)(cid:90)
σ(x,y,z)=v F(k ,k ,k )exp[j(k x−k y+k z)]dk dk dk . (22)
x y z x y z x y z
The constant factor relative to 2π is ignored in (22).
3. THREE-DIMENSIONAL IMAGING PROCEDURE
The computational procedure of the algorithm which is proposed in
Section 2 is summarized in Figure 2 and is outlined in detail in the
following steps.
Step 1: Remove RVP. Removing RVP is not the key parts in this
paper, so it is not addressed.
528 Ren et al.
Step 2:PerformthefastFouriertransform(FFT)ofs (θ(cid:48), y(cid:48); f)
IF
acquired in step 1 with respect to y(cid:48) for each frequency and azimuth
angle. As a result, the raw data in the y(cid:48) domain is transformed into
the k domain.
y(cid:48)
Step 3: RFM. Corrections of the phase are the emphasis of this
paper. The processing steps of the matched filtering can be split into
two parts.
1) Remove the unwanted phase in the wavenumber domain.
a) Considering the motion within the sweep, a 3-D matched filter is
given by
(cid:183) (cid:181) (cid:182)(cid:184)
(cid:161) (cid:162) f
G k ,k ; τ = exp j k y(cid:48) −k v −k vτ
F y(cid:48) r n y(cid:48) y(cid:48) y(cid:48) n
0 K
(cid:183) (cid:181) (cid:182)(cid:184)
2R
c
exp j −2k R −k v . (23)
r c y(cid:48)
c
b) Under the stop-and-go approximation, the conventional matched
filter [22] is shown as follows
(cid:183) (cid:181) (cid:182)(cid:184)
(cid:161) (cid:162) 2R
c
G k ,k ; τ = exp j −2k R −k v . (24)
F y(cid:48) r n r c y(cid:48)
c
Notethatthetermk vτ existsintheelevation-frequencydomain
y(cid:48) n
and the azimuth-time domain, therefore the term G should be first
F
multiplied exactly before the azimuth Fourier transform and after the
elevation Fourier transform, as shown in Figure 2.
2) Correct the range curvature of all scatterers at the reference
range R .
0
In this paper, R is considered as the reference range. The
0
correction progress should be calculated in the k domain in the
θ(cid:48)
proposed algorithm, such as the following steps.
a) Calculate the FFT of S(cid:48) (k , k ; θ(cid:48)) acquired in step 3 with
IF y(cid:48) r
respect to θ(cid:48), so the function S(cid:48) (k , k ; k ) is obtained.
IF y(cid:48) r θ(cid:48)
b) Obtain the FFT of the function G∗ (k , θ(cid:48)) with respect to θ(cid:48),
kxz xz
so as to yield the function G∗ (k , k ).
kxz xz θ(cid:48)
c) TakethecomplexproductofS(cid:48) (k , k ; k )and1/G∗ (k , k )
IF y(cid:48) r θ(cid:48) kxz xz θ(cid:48)
in the k domain.
θ(cid:48)
d) Perform the IFFT of the product with respect to k in order to
θ(cid:48)
obtain F (k , θ(cid:48), k ).
kxz xz y(cid:48)
In order to further alleviate the computational load of the
algorithm, the term 1/G∗ (k , k ) which is defined as the focusing
kxz xz θ(cid:48)
function[17]inthispapercanbecalculatedfirstandstoredinmemory.