Micro-motion frequency estimation of radar targets with complicated translations

Int. J. Electron. Commun. (AEÜ) 69 (2015) 903–914

Contents lists available at ScienceDirect

International Journal of Electronics and Communications (AEÜ) journal homepage: www.elsevier.com/locate/aeue

Micro-motion frequency estimation of radar targets with complicated translations Wenpeng Zhang ∗ , Kangle Li, Weidong Jiang School of Electronic Science and Engineering, National University of Defense Technology, Changsha 410073, China

a r t i c l e

i n f o

Article history: Received 30 August 2014 Accepted 18 February 2015 Keywords: Micro-Doppler Micro-motion Parameter estimation Complicated translations Time-frequency squared difference sequences

a b s t r a c t Micro-Doppler (m-D) signatures induced by micro-motion dynamics, which are of great importance for target classification, have received increasing attention among the radar community. However, most of the existing m-D signature-extraction methods are based on the assumption that translations of radar targets have been completely compensated by other methods or they merely involve scenarios when radar targets undergo low order polynomial translations. These methods may become invalid for maneuvering targets with complicated translations. In this work, radar targets with micro-motion and complicated translations are considered. Inspired by the differential element method and the shift invariant difference operation, a micro-motion frequency estimation method based on piecewise translation compensation (PTC) and time-frequency squared difference sequences (TFSDS) is proposed. The PTC can compensate the translations of radar targets at the largest extent. After compensation, the TFSDS-based frequency estimator gives an estimation of micro-motion frequency; this estimator can remove the effect of residual translation and avoid the separation of the multicomponent radar echo. The combination of PTC and TFSDS enables us to extract the micro-motion frequency of radar targets with complicated translations without a prior knowledge. Experiments with synthetic and measured data confirm the effectiveness and good performance of the proposed method. © 2015 Elsevier GmbH. All rights reserved.

1. Introduction The radar target undergoing micro-motion dynamics imposes a periodic time-varying frequency modulation on radar signals. This is known as the micro-Doppler (m-D) effect [1]. The m-D features, which reflect the unique dynamic and structural characteristics of the target, serve as additional target features for target recognition and classification that are complementary to those made available by existing methods [1–3]. Extraction of m-D feature for maneuvering space targets is a hot topic in recent years, application can be found in the research of ballistic warhead, decoy, etc. [4]. For such maneuvering targets which are always with high velocity and acceleration, their translations shift the m-D and disturb its periodicity. Furthermore, high accelerations lead to a large dynamic frequency range and even Doppler ambiguity. The complicated instantaneous frequency (IF) laws of radar signals introduced by complicated translations also cause mD estimation bias. On the other hand, radar echoes of targets with

∗ Corresponding author. Tel.: +86 18609472312. E-mail addresses: [email protected] (W. Zhang), [email protected] (K. Li), [email protected] (W. Jiang). http://dx.doi.org/10.1016/j.aeue.2015.02.011 1434-8411/© 2015 Elsevier GmbH. All rights reserved.

micro-motion are nonstationary and multicomponent. As a result, when the acceleration is high enough, components of the radar echo become unresolvable in time-frequency (TF) domain because of their interferences. Therefore, it is essential to compensate the translation of the radar target for most applications to resolve the Doppler ambiguity and to improve the estimation precision of m-D. Numerous m-D signature-extraction methods were proposed in the past decades. As for non-rigid radar targets, the main body signatures and m-D signatures are overlapped in the TF domain. Typically, signal decomposition methods which decompose radar signals into components associated with different parts of the radar target, e.g., chirplet decomposition, wavelet decomposition, empirical mode decomposition (EMD), can be carried out to extract m-D signature [5–7]. Another category of methods for extracting m-D signatures of non-rigid targets is based on the statistics of time-frequency distribution (TFD) [8]. In [8], an L-statisticsbased method for m-D effects removal is proposed; it produces better focused images of the rigid body than the other TF-based approaches. The m-D signatures present periodicity on the TF plane, thus methods based on the periodicity of the TFD of radar signals can be developed to estimate the micro-motion period [9–12]. Taking advantage of this property, two approaches based on the

904

W. Zhang et al. / Int. J. Electron. Commun. (AEÜ) 69 (2015) 903–914

one-dimensional Fourier transform [9] and the two-dimensional Fourier transform [10] of the TFD are proposed to extract the micro-motion period (frequency), respectively. In [11], a similarity measure of the TFD is employed to estimate the micro-motion period, which can be applied to different types of m-D signals. In [12], a method based on the mixed m-D time-frequency data sequences is proposed to estimate micro-motion parameters of free rigid targets, including spin rate, precession rate, nutation angle and inertia ratio. The construction of mixed m-D TF data sequences plays a crucial role in the performance of the method. Recently, the authors show that the radar echo in the presence of translation presents circular periodicity in TF domain. In view of this, the circular correlation coefficients of the TFD are employed to characterize the circular periodicity of the TFD and to provide an estimate of the micro-motion period. But its performance for complicated translation remains to be tested [13]. The aforementioned methods are mostly based on the assumption that the target undergoes a low order (usually second order) polynomial translation or that the translation is compensated by other methods. As for translation compensation, the Doppler rate is used to remove the translation [14], but it merely compensates the translation for single scatterer-radar targets or it need to separate the multicomponent radar echo in the TF domain. In [15], a translation compensation method based on the minimization of the spectrum entropy of the radar echo is proposed; it can compensate high accelerationtranslations without the separation of the multicomponent radar echo if the translation of the target can be approximated by a second polynomial function. An EMD-based translation compensation method which decomposes the IFs of the radar echo into intrinsic mode functions (IMFs) associated with the micro-motion and the translation is proposed in [16]. The proposed method in [16] can compensate the radar echo with multi-scatterers and seems to be effective for complicated translations. However, it becomes invalid when there is frequency aliasing in the TF domain. These methods may not be effective for scenarios when the radar target undergoes a complicated translation, especially when there is no a prior knowledge about the translation. Aiming at m-D feature extraction of radar targets with complicated translations, a micro-motion frequency estimation method based on piecewise translation compensation (PTC) and timefrequency squared difference sequences (TFSDS) is proposed in this paper. The PTC can compensate the translation of the radar target at the largest extent. After that, a TFSDS-based frequency estimator gives an estimation of micro-motion frequency, which can remove the effect of residual translation and avoid the separation of the multicomponent time-varying and intersecting micro-Doppler components. The remainder of this paper is organized as follows. In Section 2, the general model of the m-D of the radar target with a translation is introduced. Next, the derivation of the estimation algorithm is presented in Section 3, including the methodologies of PTC and TFSDS. Experiments with synthetic and measured data of cone-shaped targets validate the proposed method in Section 4.

2. The micro-Doppler model In optical region, a radar target can be modeled as a set of point scatterers and thus its radar echo is the combination of returns from these point scatterers. The motion of a scatterer is composed of the translation of the target (also refer to as macro-motion, which is identical for all scatterers) and its individual micro-motion. The typical micro-motions whose forms are sinusoidal functions with respect to the micro-motion frequency are employed during the derivation, while the form of translation is arbitrary and unknown.

Assuming that there are L scatterers on the target, the radial distance of the lth scatterer from radar associated with micro-motion can be expressed as: rM l (t) = Al sin(2fm t + ϕl ) + rl0

(1)

where Al is micro-motion amplitude, fm is micro-motion frequency, ϕl is initial phase, rl0 is the initial distance to radar. And total radial distance of the lth scatterer is rl (t) = rT (t) + rM l (t)

(2)

where rT (t) is the radial distance associated with the translation of the target. In general, the baseband signal of the returned radar signal can be expressed as s(t) =

L 

l (t) exp

 j4r (t)  l



l=1

= exp

L  j4r (t)   T



l (t) exp

 j4r

M l (t)





= sT (t)sM (t)

(3)

l=1

where sT (t) = exp

sM (t) =

 j4r (t)  T

L 

(4)

 l (t) exp

 j4r

M l (t)





(5)

l=1

sT (t) is the translation modulation part, sM (t) is the micromotion modulation part,  l (t) denotes the RCS of the lth scatterer,  is the wavelength of the carrier wave and L is the number of scatterers. The Doppler frequency of the lth scatterer is given by fl (t) =

2 drl (t) 2 drT (t) 2 drM l (t) 2 drT (t) = + =  dt  dt  dt  dt +

4fm Al cos(2fm t + ϕl ) 

(6)

where (2/)(drT (t)/dt) is the translation Doppler, (4fm Al /) cos(2fm t + ϕl ) is the m-D. The Doppler frequency of the lth scatterer can be simplified as follows: fl (t) = Al sin(2fm t + ϕl ) + ftrans (t)

(7)

where ftrans (t) = (2/)(drT (t)/dt) denotes the translation Doppler, fm is the micro-motion frequency of the target, and Al = 4fm Al / and ϕl = ϕl + /2 are the amplitude and initial phase of the m-D, respectively. As indicated in (3), (6) and (7), the radar echo of the radar target with micro-motions is a multicomponent amplitude and frequency modulated (AM-FM) signal. 3. The proposed estimation algorithm Because of the multicomponent and time-varying frequency modulation properties, radar echoes are overlapped in the TF domain, which makes it difficult to separate components associated with different scatterers. On the other hand, approaches based on the assumption that radar targets undergo a low order polynomial translation may become invalid for radar targets with complicated translations. Our approach contains two main steps: (1) PTC and TFD generation; (2) TFSDS-based micro-motion frequency estimation. At the first stage, the translation is compensated by PTC accompanied by a special TFD Generation strategy. The special TFD Generation strategy can reduce the end effect of TFD which results from the segmenting of the radar echo. After that, the TFSDS is utilized to

W. Zhang et al. / Int. J. Electron. Commun. (AEÜ) 69 (2015) 903–914

remove the residual translation and to estimate the micro-motion frequency of the radar target.

Searching the minimum entropy gives the optimal value of an . aˆ n = argminε (an ) an

3.1. Piecewise translation compensation and TFD generation

1 an t 2 2

1 ≤ n ≤ N,

(n − 1)T ≤ t ≤ nT

(8)

If the translation of the nth subsection is a strict second order polynomial function, then fR n (t) is equivalent to the actual translation of the nth subsection; the parameters in (6) have real physical meanings, which are known as initial range, initial velocity, and acceleration, respectively. Note that rn only generates a constant phase and does not affect the compensation effect of the nth subsection; therefore we set rn = 0. The nth subsection after compensation can be written as sn (t) = sn (t) exp

 −j4 · f

R n (t)



(9)



The parameters in (8) are calculated sequentially. The variation of translation makes the radar echo spread over a large frequency range. If the translation is compensated completely, the bandwidth of the radar echo will reduce to the bandwidth of the m-D. Therefore, the spectrum entropy which reflects the concentration of the spectrum of the radar echo can be used to evaluate the compensation effect. As the linear term vn t does not affect the spectrum entropy, we initially set vn = 0 and estimate an based on the minimization of the entropy of spectrum [15]. Considering the signal in frequency domain, the Fourier transform of sn (t) is Sn (f ) = F(sn )

(10)

The entropy of Sn (f) is defined as



∞

ε(an ) = −∞

|Sn (f )|2 ES ln df ES |Sn (f )|2

(11)

where



∞

|Sn (f )|2 df

ES =

(12)

−∞

According

 nT

(n−1)T 

to

the

Parseval

law

ES =

 nT (n−1)T

|sn (t)|2 dt =

|sn (t)|2 dt, the entropy can be redefined as



∞

Sn (f )|2 ln |Sn (f )|2 df

ε (an ) = − −∞

(14)

Then vn can be estimated using (15).

As for the complicated translation, it is impossible to give an analytical formula of the translation without a prior knowledge or to approximate the translation by a low order polynomial function throughout the observation time. Inspired by the differential element method in mathematics, the radar echo is segmented into subsections and compensation is conducted on each subsection. When the time interval of each subsection is small enough, the translation of each subsection can be approximated by a second order polynomial function. A compensation function with second order polynomial phase is then multiplied to compensate the translation of each subsection, whose parameters are determined through minimizing the entropy of the spectrum [15]. The segmenting of the radar echo enables us to address complicated translations. In a word, we use a piecewise second order polynomial function to compensate the translation of the target; this makes the compensation achieved more easily and more effectively. The radar echo with an observation time T is equally segmented into N subsections, T = T/N is the interval of each subsection. A second order polynomial function is used to approximate the translation of nth subsection. fR n (t) = rn + vn t +

905

(13)

vˆ n =

fmax  2

(15)

where fmax = argmax|Sn (f )|an =ˆan

(16)

f

As a result of segmenting, the IF law of each scatterer is no more continuous, but it does not affect the estimation of micromotion frequency in our approach (discussed in Section 3.2). If we impose the continuity constraint on the compensation function, rn and vn should be calculated according to the continuity condition. As a consequence, the translation of each subsection cannot be compensated at the largest extent, which leads to a poor compensation performance in terms of the whole radar echo. The break of continuity enables us to compensate the translation as largely as possible. However, the discontinuity affects the precision of the TFD. In order to improve the precision of the TFD, a special TFD generation approach is adopted. After estimating parameters an and vn , we calculate the TFD of the nth subsection with extended data and results at additional instants are abandoned. This can avoid the influence of discontinuity and the end effect of TFD (results from segmenting). Fig. 1 is the demonstration of PTC and TFD generation for nth subsection. Line 1, Line 2 and Line 3 stand for the radar echo, the nth subsection, the extended nth subsection, respectively. Matrix 1, Matrix 2, Matrix 3 stand for the TFD of the radar echo, the TFD of the nth subsection, the TFD of the extended nth subsection, respectively. The length of lines and width of matrices stand for intervals of these signals. The details of the PTC and TFD generation procedure are presented as follows. Input: number of subsections N, extended ratio p, searching values of an , denoted as ˝. Output: TFD of the radar echo (t, f). Initial: n = 1, T = T/N Step 1: Extract the nth subsection sn (tn ) and the extended nth subsection sn e (tn e ); For sn (t), (n − 1)T ≤ tn ≤ nT; for sn e (tn e ), min {0, (n − 1)T − pT} ≤ tn e ≤ max {nT + pT, T}. Step 2: Estimate aˆ n for nth subsection sn (tn ); For each value an in ˝, calculate the entropy ε (an ) defined by (11), estimate aˆ n via (12). Step 3: Estimate vn via (13); Step 4: Compensate the extended  nth subsection sn e (t n e ); 2 v ˆ t + (1/2)ˆ a t n n e n n e sn e (tn e ) = sn e (tn e ) exp(−j4 · )   Step 5: Calculate TFD of sn e (tn e ), denoted as n e (tn e , ftn e , ftn e , f); Step 6: Find the corresponding TFD of nth subsection; (t, f ) = n e (tn e , f ) (n − 1)T ≤ t, tn e ≤ nT Step 7: n = n + 1, if n ≤ N, repeat Step 1 to Step 6; otherwise, stop the algorithm. Other issues of the PTC approach should be emphasized as well. The first is that there may be a residual high order translation after compensation as a result of the non-continuous compensating function. The second is the influence of the subsection number N. If N is too small, the translation of each subsection may not be approximated by a second order polynomial function. In that case, the proposed method may fail to address complicated translations. On

906

W. Zhang et al. / Int. J. Electron. Commun. (AEÜ) 69 (2015) 903–914

Fig. 1. Demonstration of PTC and TFD generation on nth subsection.

the other hand, if N is too large, the spectrum entropy fails to characterize the concentration of the Doppler spectrum of the radar echo (because the estimation of the spectrum entropy of a short signal is inaccurate.). One particular case is that when N equals the sample number of the radar echo (each subsection only has one sample), the spectrum entropy of each subsection will remain zero. In that case, the entropy of the compensated radar echo of each subsection with arbitrary compensating function would remain constant and the proposed method fails to compensate the translation. Therefore, the subsection number should be chosen properly. Relevant numeric experiments are presented in Section 4. 3.2. Time-frequency squared difference sequences-based frequency estimator

trans (t)

One remaining issue is the separation of components associated with different scatterers in TF domain. The break of continuity makes it impossible to achieve this. Here, a subtle method is used to avoid this problem. The squared value of (19) is |fl c (t) − fk c (t)|2 = A2lk sin2 (2fm t + ϕlk ) = A2lk

m(t) =



|fl c (t) − fk c (t)|2 =

l= / k

trans (t)

= ftrans (t) −

1 ≤ n ≤ N(n − 1)T ≤ t ≤ nT

Am =



A4lk +

l= / k

=

sin

(2fm t + ϕl ) + Ak

sin(2fm t + ϕk + )

= Alk sin(2fm t + ϕlk )

(19)

where



Alk =

A2 + A2 − 2Ak Al cos(ϕk − ϕl ) k l

(20)

1 − cos(4fm t + 2ϕlk ) 2 (23)



2A2lk A2pq cos(2ϕlk − 2ϕpq )

(24)

(l,k) = / (p,q),l = / k,p = / q



Alk sin 2ϕlk

l= / k



+ Alk cos 2ϕlk

 2

(25)

l= / k

(18)

fl c (t) − fk c (t) = Al sin(2fm t + ϕl ) − Ak sin(2fm t + ϕk )

A2lk

where

ϕm = arctan

The identical residual translation can be removed by the subtraction of Doppler frequencies of two scatterers, i.e.,

Al

1 − cos(2ωt + 2ϕlk ) 2 (22)

= Am sin(4fm t + ϕm ) + Bm

(17)

2 (ˆvn + aˆ n t), 

 l= / k

where fr trans (t) is the residual translation. fr

(21)

Al cos ϕl − Ak cos ϕk

An m-D signal which contains micro-motion frequency information, defined as TFSDS, is given as follows

As mentioned before, after PTC stage, there is still a residual complicated translation. As a result, the periodicity of m-D signatures is still not presented in TF domain, which means that methods based on the periodicity of the TFD still cannot work. Note that the residual translations of all scatterers remain identical after PTC stage. Therefore, the shift invariant difference operation is capable of removing the residual identical translation. The TFSDS which is made up of difference operation is utilized to remove the residual translation and extract the underlying periodicity of the TFD. After compensation, the Doppler frequency can be written as fl c (t) = A l sin(2fm t + ϕl ) + fr

Al sin ϕl − Ak sin ϕk

ϕlk = arctan

Bm =

1 2 Alk 2

(26)

l= / k

The TFSDS remains constant with any permutation order of scatterers; therefore the separation of scatterers is avoided. According to (23), m(t) is the sum of sinusoidal signals of frequency 2fm with different amplitudes, phases, and constant terms, which will turn into a sinusoidal signal or a constant value. Due to the different amplitudes and phases of scatterers, the probability of (23) to be a constant value is low and thus an effective sinusoidal signal which contains micro-motion frequency information is obtained.

W. Zhang et al. / Int. J. Electron. Commun. (AEÜ) 69 (2015) 903–914

After getting (23), the micro-motion frequency can be estimated as follows

where

1 fˆm = arg max|Fm (f )|2 2 f

ε(n) =

(27)



907

(vl (n) − vk (n))2 =

L 

l= / k

(L − 1)v2l (n) −

l=1

L 

(32)

 Fm (f ) =

m(t) exp(−j2ft) dt

(28)

3.3. Procedure of the micro-motion frequency estimation method

Using the linearity of the expectation operator and the i.i.d. property of the Guassian noise, we get E{ε(n)} = L(L − 1)v2

(33)

2

2

4

E{ε (n)} = L(L − 1)(4L − 5L + 3)v Before utilizing the TFSDS-based frequency estimator to estimate the micro-motion frequency, we need to estimate all scatterers’ Doppler frequencies. This can be achieved via maxima estimation in the TF plane. Errors may exist in the maxima estimation procedure in the presence of noise. Note that the micro-motion frequency is always low and the noise occupies the high frequencies; therefore, a denoising procedure is employed to suppress the noise. Fig. 2 is the process flow of proposed estimation algorithm. The details of the proposed estimation algorithm are presented as follows. Step 1: Perform PTC and TFD generation; Step 2: Estimate the Doppler frequencies via maxima estimation; Assume that the number of scatterer is known. The maxima in TF plane are estimated sequentially. Points with the strongest energy at each instant are extracted as the first component. After a component is extracted, the neighborhood of this component in TF plane is set zero before next maxima extraction stage. Repeat until the last component is extracted. Notice that the “component” in this case is not associated with a specific scatterer. Permutation may exist in the values of different components and thus a component may be the combination of different fractions (in terms of time) of different scatterers. The components of TFD are denoted as fˆl c (n), l = 1, 2, . . ., L. Step 3: Compute the m-D signal via (23); ˆ = m(n)



|fˆl c (n) − fˆl c (n)|2

3.4. Variance and bias expressions for the TFSDS-based frequency estimator In the case of noisy observations, the discrete-time version of (17) can be written as follows

 4f

m

fs n

+ ϕl



+ fr

trans (n) + vl (n)

(30)

where fs is the sampling frequency, n = 0, 1, . . ., N − 1, N is the sample number, vl (n), l = 1, 2, ..., L are zero-mean, white Gaussian stationary noise with variance v2 . Repeating the same deriving procedure and neglecting the cross-term of signal and noise, the discrete-time version of (23) is m(n) = Am sin

Bias (fˆm ) = 0

 4f

m

Fs n

+ ϕm



+ Bm + ε(n)

(31)

(35)

According to [18], the asymptotic variance of the discrete Fourier transform based frequency estimate at high SNR and for a constant amplitude and real sinusoidal signal is

  fˆ fs

Var

24

=

(2)

2

(A2 /E{ε2 (n)})N(N 2

− 1)

(36)

where A is the amplitude of sinusoidal signal, E {ε2 (n)} is the power of white noise. Considering the signal of (31), we have



Var

2fˆm fs



=

24 (2)

2

(A2m /E{ε2 (n)})N(N 2

− 1)

(37)

Therefore, the variance of the TFSDS-based frequency estimator can be expressed as follows Var(fˆm ) =

4(2)

2

24fs2 2 (Am /E{ε2 (n)})N(N 2

− 1)

(38)

It is difficult to study the relationship between the variance and the SNR of the input signal from (38). Define the SNR of the input signal (30) as follows SNR =

Step 4: Smooth the extracted signal obtained by Step 3 to suppress the noise; Step 5: Compute the spectrum of the smoothed extracted signal; Step 6: Estimate the maximum position of the Fourier spectrum and divide it by two to get an estimation of the micro-motion frequency (Fig. 2).

(34)

ε(n) is a i.i.d. noise. According to [17], the discrete Fourier transform based frequency estimate at high SNR for i.i.d. noise is unbiased.

(29)

l= / k

fl c (n) = Al sin

2vl (n)vk (n)

1≤l
max{A2 } l

(39)

v2

Using (34) and (39), (38) becomes Var(fˆm ) =

24L(L − 1)(4L2 − 5L + 3)fs2 2

4(2)2 (SNR)2 (Am / max{A2 }) N(N 2 − 1) l

(40)

While the variance of discrete Fourier transform based frequency estimator for single sinusoidal signal with the largest amplitude is Vars (fˆm ) =

24fs2 2

(2) (SNR)N(N 2 − 1)

(41)

It is easy to prove that

⎧ Var(fˆm ) > Vars (fˆm ) SNR <

⎪ ⎪ ⎨ ⎪ ⎪ ⎩

Var(fˆm ) = Vars (fˆm )

SNR =

Var(fˆm ) < Vars (fˆm )

SNR >

(42)

where

=

L(L − 1)(4L2 − 5L + 3) }) 4(Am / max{A2 l

2

(43)

The TFSDS-based frequency estimator outperforms the discrete Fourier transform based frequency estimator for single sinusoidal signal at higher SNR. This property can be viewed as the result of

908

W. Zhang et al. / Int. J. Electron. Commun. (AEÜ) 69 (2015) 903–914

Fig. 2. Process flow of estimation algorithm.

the exploitation of more data, which uses all the sinusoidal components. 4. Experimental results and discussion 4.1. Performance of TFSDS-based frequency estimator Use the following multi-sinusoidal signals to simulate the m-D xi (n) = Ai sin

 2f

m

fs n

+ ϕi



+ ft (n) i = 1, 2, 3

(44)

where ft (n) is a random trend, which is used to simulate the complicated translation ft (n) =

k=n 

w(k)

(45)

k=0

w(n) is a random process. A1 = A2 = 100, A3 = 10, ϕ1 = 0, ϕ2 = , ϕ3 = /3. The sampling frequency is fs = 1024 Hz, the number of observation is N = 1024. The noise-free original sinusoidal signals and the corresponding signals with random trend are depicted in Fig. 3(a) and (b), respectively. As it can be seen from Fig. 3(a) and (b), the random trend in Fig. 3(b) disturb the periodicity of the sinusoidal signals.

The extracted sinusoidal signal generated by TFSDS is depicted in Fig. 3(c), which shows that the proposed method can recover sinusoidal signal (with double frequency) from multiple sinusoidal signals with an identical trend (Fig. 3(b)). The performance of the TFSDS-based frequency estimator under different SNRs is illustrated in Fig. 4. Monte Carlo simulations for 1000 realizations are run for each SNR. Additive white Guassian noise is added to each xi (n) independently. The SNR is computed by (39). The bias of the estimator depending on the SNR is illustrated in Fig. 4(a). The TFSDS-based frequency estimator provides unbiased estimate for SNR ≥ −1 dB. The estimate variance and the theoretical variance of the proposed estimator, the theoretical variance of discrete Fourier transform based frequency estimator for single sinusoidal signal are illustrated in Fig. 4(b). As it can be seen from Fig. 4(b), the estimate variance reachs at values near the theoretical variance for SNR ≥ 1 dB, which is less than −31.28 dB (7.4e−4 Hz2 ). In this case, the TFSDS-based frequency estimator outperforms the discrete Fourier transform based frequency estimator for SNR > 6.03 dB in terms of variance.

4.2. Simulation of a coning target A cone target with coning motion and complicated translation in our experiment is demonstrated in Fig. 5. The reference coordinate

Fig. 3. (a) Original sinusoidal signals, (b) sinusoidal signals with identical random trend, and (c) extracted sinusoidal signal.

W. Zhang et al. / Int. J. Electron. Commun. (AEÜ) 69 (2015) 903–914

909

Fig. 4. Bias and variance of estimator depended on SNR level.

modulation part is generated by using (4). The radial distance of the complicated translation is given by rT (t) = −9.08t + 30t 2 + sin(2 · 0.1t + 0.2 · ) + sin(2 · 0.5t + 0.4). The radar echo is generated as the product of micro-motion modulation part and translation modulation part. The simulated discrete-time radar echo with noise is given by y(n) = s(n) + w(n),

(47)

where s(n) is the sampled version of s(t); and the sampling frequency is equal to the PRF. w(n), l = 1, 2, ..., L are zero-mean, white 2 . The signal to noise ratio Gaussian stationary noise with variance w (SNR) of the noisy radar echo is calculated as

Fig. 5. Target model.

system O − XYZ and the target local coordinate system O − xyz are used to describe the motion of the target. There are three scatterers on the target, with one scatterer on the top and two scatterers on the bottom edge; their positions in O − xyz are P1 = (0, 0, 1.6) m, P2 = (0, − 0.2, − 0.4) m, and P3 = (0, 0.2, 0.4) m, while their RCS are  1 =  2 =  3 = 1. The coning frequency and coning angle of the target are fm = 1.5 Hz and  = 15◦ , respectively. The azimuth and elevation angle of the target center with respect to the radar in the reference coordinate system are ˛ = − 50◦ and = 270◦ , respectively. The radar operates at 10 GHz with a pulse repetition frequency (PRF) of 1024 Hz. Fig. 6 shows the experimental setup. Following the similar derivation in [1], the radial distances of scatterers from radar associated with micro-motion can be calculated as follows: r(t) = sin ˛(y sin  + z cos ) + cos ωt cos ˛(x cos + y sin cos  − z sin sin ) + sin ωt cos ˛(x sin − y cos cos  + z cos sin )

n = 0, 1, ..., N − 1

(46)

where (x, y, z) are the coordinates of the scatterers. According to the radial distances of scatterers and their RCS, we can generate the micro-motion modulation part by using (5), while the translation

N−1

SNR = 10 log 10 Experiment 1.

n=0

|s(n)|2

2 w

(48)

Validation of PTC and TFD generation.

A noisy radar echo with a SNR of 5 dB is considered. In order to demonstrate the performance of the PTC and TFD generation methods, we carried out experiments with and without continuity constraint as well as experiments with and without the proposed TFD generation approach. The hamming window with length 65 is adopted in the calculation of STFT. The search interval of an is [−40, 40] m/s2 with resolution of 0.1 m/s2 . The PTC is performed with subsection number N = 20. Fig. 7 shows the experimental results. The STFT of the radar echo without translation is plotted in Fig. 7(a); there are three sinusoidal curves in the TF plane. Fig. 7(b) shows the STFT of the radar echo, we can see that the translation of the target induces Doppler ambiguity and disturbs the periodicity of m-D. Besides this, the interferences between different components are serious, which makes the components of the radar echo inseparable at some instants, e.g., instant near 0.1 s, 0.6 s. Fig. 7(c) and (d) are the STFT of the radar echo after PTC with and without continuity constraint, respectively. As we can see from Fig. 7(c), the

Fig. 6. Experimental setup.

910

W. Zhang et al. / Int. J. Electron. Commun. (AEÜ) 69 (2015) 903–914

Fig. 7. (a) STFT of radar echo without translation, (b) STFT of radar echo, (c) STFT of radar echo after PTC with continuity constraint, N = 20, (d) STFT of radar echo after PTC without continuity constraint, N = 20, and (e) STFT of radar echo after PTC without continuity constraint, the STFT is calculated via our TFD generation method, N = 20. SNR = 5 dB.

PTC with continuity constraint compensates part of the translation of the target; but there is still a residual translation and it causes Doppler ambiguity in the TF domain. As depicted in Fig. 7(d), the PTC without continuity constraint can compensate more of the translation than the previous one with continuity constraint; although there is also a residual translation, it obtains a better compensation result without Doppler ambiguity. Fig. 7(e) is the STFT of the radar echo without continuity constraint which is calculated via the proposed TFD generation approach. Comparing the results in Fig. 7(d) and (e), it is clear that the STFT calculated via the proposed TFD generation method has better TF concentration and more accurate frequency estimation of the components. In the following experiments, the PTC stands for PTC without continuity constraint and the STFT is calculated via the proposed TFD generation method if without special declarations. Fig. 8 shows the result of the PTC and TFD generation approach under different subsection numbers. Fig. 8(a) is the STFT of the radar echo after PTC, Fig. 8(b) is the residual translation frequency under

different subsection numbers, N ∈ [1, 400] with an interval of 2. The residual frequency is calculated as follows:

   M−1  rf 2 (i) 

RMSEf =

rf (i) =

i=0

⎧ ⎪ ⎨ fs − |fr ⎪ ⎩ |fr

(49)

M trans (i)%fs |

trans (i)%fs |

|fr

trans (i)%fs |

≥

fs 2

|fr

trans (i)%fs |

<

fs 2

(50)

where i is time instant, M is the length of the discretized radar echo, % stands for the modulo operation, i.e. the operation which finds the remainder after division. From Fig. 8(b), we can see that the residual translation frequency between interval [17, 127] is significantly small; it reaches its minimum value 103 Hz at N = 61. This result indicates that the PTC works

W. Zhang et al. / Int. J. Electron. Commun. (AEÜ) 69 (2015) 903–914

911

Fig. 8. (a) STFT of radar echo after PTC, N = 60 and (b) residual translation frequency under different subsection numbers. SNR = 5 dB.

Fig. 9. (a) Auto correlation of radar echo, (b) STFT of radar echo after PTC, (c) Doppler frequencies of three scatterers, (d) extracted m-D signal before smoothing, (e) extracted m-D signal after smoothing and (f) spectrum of m-D signal. SNR = 5 dB.

912

W. Zhang et al. / Int. J. Electron. Commun. (AEÜ) 69 (2015) 903–914

Fig. 10. (a) STFT of radar echo, (b) STFT of radar echo after PTC, (c) spectrum of m-D signal. SNR = −2 dB, (d) RMSE of estimation error in different SNR cases.

effectively in a suitable section number interval. Too small or too large subsection number would lead to poor compensation performance. This result is consistent with the discussion in the end of Section 3.1. The micro-motion estimate errors with the same subsection number interval will be shown in the following experiment. Experiment 2. mation.

Performance of the micro-motion frequency esti-

In this experiment, the performance of the micro-motion frequency estimation is evaluated. The median filter with a window of 50 is used to suppress the noise. Fig. 9(a) is the auto correction function of the radar echo, indicating that methods utilizing auto correlation to estimate the micro-motion frequency of the target with complicated translations are invalid. Fig. 9(b) depicts the STFT of the radar echo after PTC. After compensation, the m-D effect is clear although the translation is not completely compensated. The Doppler frequency of each scatterer is not continuous after PTC as mentioned before. Fig. 9(c) demonstrates the Doppler frequencies of three scatterers extracted by the maxima estimation procedure. In this experiment, the maxima estimation is performed for three times to extract the three components of the radar echo. Points of the same color are extracted at the same maxima estimation stage. Notice that the Doppler frequencies near the intersecting region are not accurate and there are permutations among different components. But this does not affect the proposed micro-motion frequency estimation approach. Fig. 9(d) and (e) are the extracted m-D signals calculated based on the TFSDS before and after smoothing, respectively. As seen in Fig. 9(e), noise is partly suppressed after smoothing and a clear sinusoidal m-D signal is obtained. Fig. 9(f) shows the Fourier spectrum of the extracted m-D signal which reaches its peak at 3 Hz (2fm ), and thus the micro-motion frequency is estimated as 1.5 Hz, which agrees with the actual micro-motion frequency.

Fig. 11. Relative error of estimate under different subsection numbers. SNR = 5 dB.

The smoothing process is taken to suppress noise; thus the proposed method has the potential to work in low SNR case. With all values of other parameters remaining unchanged except that the SNR is set to −2 dB, we illustrate the results of the proposed method in Fig. 10(a)–(c). The Fourier spectrum of the extracted m-D signal has a maximum position of 3 Hz (2fm ), which is in accordance with the actual micro-motion frequency. Monte Carlo simulation under different SNRs is conducted to evaluate the performance of the proposed method. The SNR changes from −5 to 18 dB with an interval of 2 dB. The simulation times of each SNR are 100. Fig. 10(d) shows the experimental result. As depicted in Fig. 10(d), the RMSE of estimate decreases dramatically from SNR = −5 dB to −1 dB; after −1 dB, it decreases at a slow rate and tends to reach 0. The estimate error of SNR larger than −1 dB is less than 0.048 Hz in terms of RMSE. In the following experiment, the relative error of estimate is used to characterize the effectiveness of the proposed method under different subsection numbers. Fig. 11 shows the result. From Fig. 11, we can see that the proposed method is effective with a

W. Zhang et al. / Int. J. Electron. Commun. (AEÜ) 69 (2015) 903–914

913

Doppler ambiguity in the TF plane. This performance can be viewed as the result of the residual translation-removing capability of the TFSDS-based frequency estimator as well as the property of Fourier transform. However, it should be emphasized that the translation of the radar target should be compensated to get a TFD free of Doppler ambiguity and to get an accurate estimate of instantaneous frequencies of the scatterers. Otherwise, the statistical performance of micro-motion frequency estimation would be poor. 4.3. Experiment with measured radar data

Fig. 12. RCS vs. aspects.

large subsection number interval. In this case, the interval is [5, 341]. Compared to Fig. 8(b), this range is much larger than [17, 127], which means that the residual translation has little impact on the micro-motion frequency estimation as long as there is not

Chen et al. proposed that the motion of the object can be considered as a sequence of snapshots taken at each time instant, the variation of the EM backscattering field as if the object were static can be estimated [19]. In this experiment, this method is used to generate the radar echo. A cone-shaped target is adopted in this experiment. Fig. 12 shows the RCS vs. aspect angles. The target was measured with stepped frequency signals at 10 GHz in a microwave anechoic chamber. In this experiment, we use the same coordinate systems with Fig. 5. The experimental setup is the same with

Fig. 13. (a) Auto correlation of radar echo, (b) STFT of radar echo without translation, (c) STFT of radar echo, (d) STFT of radar echo after PTC, (e) extracted m-D signal after smoothing, and (f) spectrum of m-D signal. SNR = 4 dB.

914

W. Zhang et al. / Int. J. Electron. Commun. (AEÜ) 69 (2015) 903–914

Section 4.2 except the generation of the micro-motion modulation part. We assume that in real scenarios, the target aspects are merely determined by micro-motion dynamics; this would be suitable for far-field target. In short observation time, the RCS variation caused by translation is small and thus it is neglected. We first calculate the target aspects with respect to time during its motion, which is given by [1] ˇ(t) = arccos[cos  sin ˛ + sin  cos ˛ cos(ωt − )]

(51)

After that the micro-motion modulation part can be realized through (51) and Fig. 12 by converting target aspects into RCS. The micro-motion parameters are fm = 0.5 Hz,  = 15◦ ; the radial distance associated with translation is rT (t) = − t + 0.5t2 + 0.15t3 ; the PRF is 200 Hz; ˛ = − 50◦ and = 270◦ ; the rest of micro-motion parameters is the same with Section 4.2. Fig. 13(a) shows the STFT of the radar echo without translation (i.e., the micro-motion modulation part), as we can see from it, the m-D of one scatterer is non-sinusoidal, which means that the scatterer is sliding-type. Fig. 13(b) shows the auto correction function of the radar echo, indicating that methods utilizing auto correlation to estimate the micro-motion frequency of the target with complicated translations are invalid. Fig. 13(c) and (d) are the STFT of the radar echo before and after PTC, respectively. In the PTC stage, after segmenting the radar echo into N = 4 subsections, we search the possible acceleration value in [−100, 100] m/s2 with a resolution of 0.1 m/s2 . Before compensation, the m-D is severely deformed by translation; after compensation, although the translation is not completely removed, the m-D signature is clear enough for further exploitation. Fig. 13(e) is the extracted m-D signal (with double coning frequency) after smoothing. In this experiment, the type of micro-motion is not standard sinusoidal function as it can see from Fig. 13(a) but the TFSDS is effective to extract the periodicity of m-D signatures as plotted in Fig. 13(e). The Fourier spectrum of the extracted signal is depicted in Fig. 13(f). The estimate of the conning frequency is fˆm = 1.023/2 = 0.51125 Hz with a relative error of 2.23%. This experiment indicates that the proposed may have the potential to extract the micro-motion frequency of the radar target for some non-sinusoidal micro-motion forms, but further experiment and theoretical analysis are required to test this. 5. Conclusions A PTC and TFSDS-based method is proposed to estimate the micro-motion frequency of radar targets with complicated translations. The translation of the target is approximated by a piecewise second order polynomial function and this approximation makes the compensation of the complicated translation more easily and more effectively. The shift-invariant difference operation in the TFSDS enables us to remove residual translation and it avoids the separation of the multicomponent time-varying and intersecting micro-Doppler components as well. As a consequence, the proposed method in this work can address the complicated translations. Experimental results confirm the effectiveness of the proposed method. The effect of the subsection number on the proposed method is analyzed, but it is difficult to give a deterministic number interval which it is effective. Further research on how to completely compensate complicated translations based on the micro-motion frequency information should be carried out. Due to the lack of radar data, we only conducted experiments with cone-shaped target. Experiments and performance for other types of radar targets with micro-motions should be tested in the future.

Acknowledgement This work was supported in part by the National Natural Science Foundation of China under Grants 61302147, 61201334 and 61422114. References [1] Chen VC, Li F, Ho SS, Wechsler H. Micro-Doppler effect in radar: phenomenon, model and simulation study. IEEE Trans Aerosp Electron Syst 2006;42:2–21. [2] Molchanov P, Egiazarian K, Astola J, Totsky A, Leshchenko S, Jarabo-Amores MP. Classification of aircraft using micro-Doppler bicoherence-based features. IEEE Trans Aerosp Electron Syst 2014;50:1455–67. [3] Lei P, Wang J, Guo P, Cai D. Automatic classification of radar targets with micromotions using entropy segmentation and time-frequency features. AEU Int J Electron Commun 2011;65(10):806–13. [4] Pan X, Wang W, Liu J, Ma L, Feng D, Wang G. Modulation effect and inverse synthetic aperture radar imaging of rotationally symmetric ballistic targets with precession. IET Radar Sonar Navig 2013;7:950–8. [5] Thayaparan T, Abrol S, Riseborough E. Micro-Doppler feature extraction of experimental helicopter data using wavelet and time-frequency analysis. In: International conference on radar systems. 2004. [6] Li J, Ling H. Application of adaptive chirplet representation for ISAR feature extraction from targets with rotating parts. IET Proc Radar Sonar Navig 2003;150:284–91. [7] Cai C, Liu W, Fu JS, Lu L. Empirical mode decomposition of micro-Doppler signature. In: IEEE international radar conference. 2005. p. 895–9. [8] Stankovic L, Thayaparan T, Dakovic M, Popovic-Bugarin V. Micro-Doppler removal in the radar imaging analysis. IEEE Trans Aerosp and Electron Syst 2013;49(2):1234–50. [9] Molchanov P, Totsky A. On micro-Doppler period estimation. In: 2013 19th international conference on control systems and computer science. 2013. p. 325–30. [10] Liu L, Mclernon D, Ghogho M, Hu W, Huang J. Ballistic missile detection via micro-Doppler frequency estimation from radar return. Digital Signal Process 2012;22:87–95. [11] Dakovic M, Brajovic M, Thayaparan T, Stankovic L. An algorithm for microDoppler period estimation. In: 20th telecommunications forum TELFOR. 2012. p. 851–4. [12] Lei P, Sun J, Wang J, Hong W. Micromotion parameter estimation of free rigid targets based on radar micro-Doppler. IEEE Trans Geosci Remote Sens 2012;50(10):3776–86. [13] Zhang W, Li K, Jiang W. Parameter estimation of radar targets with macromotion and micro-motion based on circular correlation coefficients. IEEE Signal Process Lett 2015;22(5):633–7. [14] Yang YC, Tong NN, Feng CQ, He SS. Micro-Doppler extraction of rotating targets based on Doppler rate. In: 2011 IEEE international conference on signal processing, communications and computation (ICSPCC). 2011. p. 1–4. [15] Yu W, Wang J. Phase adjustment for extraction of micro-motion information of ballistic targets. In: 2012 5th international congress on image and signal processing (CISP). 2012. p. 1837–40. [16] Zhao Y, Wang G, Bai Y, Zhang Q, Li R. A new method of translational compensation for space spinning target based on EMD. In: 2013 IEEE international conference on signal processing, communications and computing (iCSPCC). 2013. p. 1–5. [17] Barkat B, Boashash B. Instantaneous frequency estimation of polynomial FM signals using the peak of the PWVD: statistical performance in the presence of additive Gaussian noise. IEEE Trans Signal Process 1999;47:2480–90. [18] Rife DC, Boorstyn RR. Single tone parameter estimation from discrete-time observations. IEEE Trans Inf Theory 1974;IT-20:591–8. [19] Chen VC, Lin C, Pala WP. Time-varying Doppler analysis of electromagnetic backscattering from rotating object. In: IEEE radar conference. 2006. p. 807–12. Wenpeng Zhang was born in Guangdong, China, in 1989. He received the B.S. degree in information engineering from the National University of Defense Technology (NUDT), Changsha, China, in 2012. His is currently working toward the master degree in information and communication engineering from NUDT. His research interests include radar signal processing, time-frequency analysis, and signal decomposition. Kangle Li was born in Anhui, China, in 1983. She received the B.S. and Ph.D. degrees in information and communication engineering from the National University of Defense Technology (NUDT), Changsha, China, in 2006 and 2010, respectively. She is currently a lecturer with School of Electronic Science and Engineering, NUDT. Her research interests include radar signal processing, time-frequency analysis, and automation target recognition. Weidong Jiang was born in Chongqing, China, in 1968. He received the B.S. and Ph.D. degrees in information and communication engineering from the National University of Defense Technology (NUDT), Changsha, China, in 1991 and 2001, respectively. He is currently a Professor with School of Electronic Science and Engineering, NUDT. His research interests include radar signal processing, radar systems, and automation target recognition.