Time-Frequency Feature Extraction of HRRP Using AGR and NMF for SAR ATR

A new approach to classify synthetic aperture radar (SAR) targets based on high range resolution profiles (HRRPs) is presented. Features from each of the target HRRPs are extracted via the nonnegative matrix factorization (NMF) algorithm in time-frequency domain represented by adaptive Gaussian representation (AGR). Firstly, SAR target images have been converted into HRRPs. And the time-frequency matrix for each of HRRPs is obtained by using AGR. Secondly, the time-frequency feature vectors are extracted from the time-frequency matrix utilizing NMF. Finally, hidden Markov models (HMMs) are employed to characterize the timefrequency feature vectors corresponding to one target and are used to being the recognizer. To demonstrate the performance of the proposed approach, experiments are performed in the 10-target MSTAR public dataset. The results support the effectiveness of the proposed technique for SAR automatic target recognition (ATR).


Introduction
Automatic target recognition has become an important research topic for military applications.It is important to develop efficient and robust algorithms for automatic target recognition (ATR), especially for ground target identification in battle surveillance.In this paper, the problem of synthetic aperture radar (SAR) ATR is focused on, which is widely used in many surveillance tasks owing to its all-weather ability and other advantages [1].Different from passive vision system, SAR images are formed by reflections of a coherent source, which are difficult to be interpreted directly and their characteristics vary quickly and abruptly with small change in azimuth and depression.Due to the unique characteristics of SAR image formation process, such as specular reflection and multiple bounces, it is very difficult to extract effective features for ATR as used in the optical image.For SAR ATR, a broad class of feature extraction method is from twodimensional SAR images [2][3][4], while an alternative class of feature extraction method is from one-dimensional high range resolution profiles (HRRPs) obtained from SAR [5].The latter have more advantages when SAR target images are blurred due to target motion.However, due to the large footprint of the SAR illuminating beam, it is very difficult to separate target HRRPs directly from SAR raw echoes which include a large amount of ground clutter.In contrast, HRRPs can be converted from SAR images with ground clutter removed by target segmentation in image domain [5,6].
Features from HRRPs via Relax algorithm are investigated for target recognition [5].Power spectrum features extracted from HRRPs are employed for SAR ATR [6].HRRPs superresolution scattering centers features are also used to recognize SAR targets [7].In the HRRPs-based SAR ATR approach, the key is to extract robust and effective features of HRRPs.
In terms of HRRP feature extraction, it has been found that the exploitation of target time-frequency signatures can be effective for target discrimination.Kim et al. derived geometrical moments features in joint T-F domain [8].Raj et al. develop methods for the T-F analysis of human gait radar signals [9].T-F analysis techniques also successfully apply to other radar applications as shown in [10][11][12].
The most critical issue when using time-frequency features for radar target recognition is to reduce the dimension of the time-frequency plane data while preserving as much discriminative information as possible.Our focus in this paper concentrates on developing a new feature extraction technique for T-F domain based on nonnegative matrix factorization (NMF).Over the last decade, NMF has emerged as a useful feature extraction method in areas related to speech recognition and image processing [13][14][15].
In this paper, we propose a new SAR ATR strategy based on adaptive Gaussian representation (AGR) and NMF.First, we propose to construct the time-frequency matrix of HRRP by using AGR.Then, NMF is performed on the time-frequency matrix to obtain the coefficient matrix composed of spectral vectors and the base matrix composed of temporal vectors.The nonnegative constraints of NMF lead to part-based representation because they allow only additive combination, which can get distinct information of the time-frequency matrix.Finally, we extract several novel features from the spectral and temporal vectors in a way that they successfully represent joint time-frequency signatures of the HRRP.Experiments based on the proposed approach are performed with HMM classifier over 10-target MSTAR public datasets.

SAR Images to HRRPs
The MSTAR public SAR dataset is considered for the experiments to evaluate the proposed algorithm.The data typically consist of image chips.As discussed above, these image chips have to be converted into HRRPs by several filtering operations.The brief procedure of how MSTAR SAR image chips are converted to HRRPs is provided here as follows [5], which is summarized in Figure 1.
Consider a complex-valued SAR image I(, ), where  reflects the downrange dimension and  the cross range.A two-dimensional (2D) inverse FFT is taken off I(, ) to obtain the corresponding phase history data.Next, the deconvolution of the weighting and removal of the zeropadding is performed for the phase history data due to the operation of window weighting and zero-padding in SAR image formation.
Then, a 2D FFT is applied to produce a deconvolved and Nyquist-sampled image I  (, ).Note that both the target and surrounding clutter exist in I  (, ).So, to remove the clutter, a target segmentation procedure is taken to image I  (, ).Then, an inverse FFT is performed in the cross range dimension for all , from which each -dependent waveform, for a fixed , corresponds to a HRRP.

HRRP Time-Frequency
Representation by AGR where The adjustable parameters   ,   , and   for Gaussian basis functions and   for the coefficient can be obtained such that   () is most similar to ℎ  (): where ℎ  () is the remainder after the orthogonal projection of ℎ −1 () onto  −1 () and this iterative procedure is described as Since the projection integral in (3) is the Fourier transform of ℎ  () with the Gaussian window () = (  ) −0.25 exp[−(( −   ) 2 /2  )], the adjustable T-F center (  ,   ) and associated variance   can be obtained using FFT and the specific search procedure in (2).  ,   ,   , and   finally obtained give the solution of (3) and these four parameters completely describe one Gaussian T-F basis function at the th iteration.
After  max stages of AGR decomposition, the following relationships hold: (5) Therefore, the AGR iteration in (4) continues until the reconstruction error ‖ℎ  max +1 ()‖ 2 is sufficiently small; hence, the upper limit  max is determined.
After   ,   ,   , and   ,  = 0, 1, 2, . . .,  max , are obtained via AGR processing, the T-F matrix, which represents a signal energy distribution in the joint T-F plane, (, ), is given by The dimension of each (, ) is the same, which is determined by time domain sample number of each HRRP and frequency domain sample number in AGR.In this paper, each (, ) is one matrix of 100×100.Compared to the wellknown WVD technique, T-F matrix by AGR can give a joint T-F distribution that is nonnegative, cross-term interference free [8].As for comparison to the wavelet decomposition which form a complete set, AGR is not restricted to a regular sampling grid, and in addition, they are not subject to the number of data samples.In radar target electromagnetic feature extraction, the scattering mechanisms are usually too complicated, and consequently, for accurate representation of a radar signature, it is desirable to have the elementary functions on a flexible sampling grid as in AGR processing rather than the elementary functions on a regular grid as in the wavelet decomposition methods mentioned above.The advantage of AGR processing for radar applications has been well described in [8].In Figure 2, T-F matrices of several HRRPs from different targets using WVD are compared with T-F matrices utilizing AGR technique, from which it can be seen that AGR is superior to WVD especially in cross-term interference.

Time-Frequency Matrix Feature Extraction Using NMF
The next stage in this paper is to derive T-F features from the T-F matrix of HRRPs.There are several matrix decomposition techniques available, such as principal component analysis (PCA), independent component analysis (ICA), and NMF.Each of these techniques considers different sets of criteria with some properties; for example, PCA finds a set of orthogonal bases that minimizes the mean squared error of the reconstructed data; ICA decomposes a dataset into components that are as independent as possible; NMF decomposes a nonnegative matrix to its nonnegative components [16].
Compared to other matrix decomposition techniques, NMF can obtain components in the original matrix with a higher representation and localization property [17,18].Therefore, the features extracted from the HRRP T-F matrix by NMF will represent a HRRP signature with a better time and frequency localization.The basic principle of NMF for the HRRP T-F matrix is to find a locally optimal factorization of the matrix into two submatrices, of which the first one named base matrix represents the spectra of the scattering events in the HRRP and the second one named coefficient matrix represents the temporal characteristic of the scattering events in the HRRP.In the present paper, following the HRRP T-F matrix decomposed by NMF, T-F features are extracted from base matrix and coefficient matrix.
4.2.NMF Factorization Algorithm.Factorization in NMF approach is usually achieved by iterative minimization of cost functions.In this work, we choose the following function: The function (W, H) has turned out to yield perceptually good results at a reasonable computational cost, which can be minimized iteratively with the multiplicative update This factorization algorithm is the basis for several recent NMF-based techniques [19].The value of the parameter  is determined by the iteration number of NMF; here,  is equal to 10 according to the experiments.Figure 3 draws out vectors of base matrix and the coefficient matrix obtained by NMF corresponding to each AGR T-F matrix in Figure 2.

NMF Feature Extraction.
In NMF feature extraction, we assume a linear signal model for HRRP, the time-frequency distribution of which can be expressed as linear combinations of spectra of several distinct scattering centers with different temporal characteristic.Thereby the coefficients are restricted to be nonnegative.One can interpret the columns of W as spectral components and the corresponding rows of H as their time-varying gains.
For spectral vectors of the basis matrix W, features are extracted utilizing several general spectral characteristic  parameters which include spectral centroid, spectral standard deviation normalized by the centroid, skewness, and kurtosis.For a given spectral vector w() of basis matrix W, spectral centroid  w , spectral standard deviation normalized by the centroid  w , skewness sk w , and kurtosis ku w are calculated as follows, respectively: where  is the dimension of w(), where  w is the standard deviation of w(), where (w) is the mean of w(), For feature extraction of temporal vectors, we calculate several time domain parameters including the root of the mean square (RMS)  h , standard deviation  h , skewness sk h , kurtosis ku h and shape factor sf h .The shape factor is defined as follows: where

𝑀 is the dimension of h(𝑖).
There are nine T-F features proposed from the significant T-F spectral and temporal components of a HRRP T-F matrix.When  = 10, both W matrix and H matrix are composed of ten vectors.

HMM Classifier
We use HMMs as the classifier for SAR targets, whose states represent the target orientation.Each of the states in the HMM corresponds to a special cover of target orientations and hence a special target feature varying with orientation.The parameters of the HMM specify a statistical characteristic of the target features.The number of HMM states needs to be chosen large enough to model the variation of the special range of target features but small enough to ensure that enough training features are available.The angular resolution corresponding to a state can be estimated from the ratio of the range resolution to the maximum target dimension.For MSTAR datasets, an angular resolution of 3 ∘ would thus lead 120 states in the HMM.
Let us assume that a target is partitioned into  distinct states, denoted by the set S = { 1 ,  2 , . . .,   }.As discussed above, each state corresponds an azimuthal partitioning.The state-transition probabilities of HMM are denoted by the matrix A = {  }, where   is the probability of transiting from state  to state .Further, the initial-state probabilities are denoted by the vector  = {  }, where   is the probability of sampling state  on the first measurement.For the sequence of HRR feature vectors, it is assumed that  represents the change in the target-sensor azimuthal orientation, upon continuous measurements.Let   represent the azimuthal angular range of state .We assume  <   , for all , implying that one of two continuous feature vectors may stay in the same state or transition in an adjacent state.This yields a tridiagonal state-transition matrix A. Based on the previous assumptions with regard to  and   , it can readily be demonstrated the following estimations for   : Moreover, it is assumed that initial target pose is uniformly distributed azimuthally, and we can have As discussed further below, ( 15) and ( 16) constitute initial estimates for A and , with these refined via Baum-Welch training algorithm.
In order to identify a target from sequential feature vectors, an HMM is designed for each target.A given set of sequential feature vectors of test target data is submitted to all HMMs and the likelihoods are computed.If the th HMM yields the largest likelihood, then we declare the sequential feature vectors are from the th target.For example, if we got the following observation sequence  = { 1 ,  2 ,  3 , . . .,   } that is associated with an unknown target , then the probability of the observation sequence  is given by summing the joint probability ( | S, )(S | ) over all possible state paths S: If the target   gives the maximum likelihood for the observation feature vector sequence , that is, then we declare the sequence  to be from the target   .
For the SAR ATR problem considered here, continuous HMMs are employed rather than discrete HMMs, because the latter has the issue of well-known distortion inherent in the quantization of feature vectors.

Database and Evaluation Methodology.
In this section, we evaluate the classification performance of the proposed approach using the MSTAR public database, which is a standard dataset for evaluating ATR algorithms, consisting of X-band SAR images with 1 ft × 1 ft resolution for multiple targets [20,21].These targets include several military vehicles and a few civilian vehicles, and they have a very similar shape.Some sample images are depicted in Figure 4 with visible light images.The data typically consist of image chips.For each target, images were acquired at several different angles over the full 360 ∘ aspect angles.These images have been converted into sequences of HRR waveform through several filtering operations as described in Section 2.
In the classification experiments, we analyze how the HMM parameters affect the identification performance and also explore how robust the proposed approach is to target variable.The training set includes total 10-target SAR images at 17degree depression angle and is used to train ten HMMs.

HMM Classification.
In the second stage of the HMM classification, the unknown target feature vector sequences are presented to the trained HMM models.For testing, we use a separate set of SAR images of 10-target SAR images at 15-degree depression angle.
Table 1 shows the confusion-matrix of classifying all test sequences of aperture 3 degrees from all the ten targets.We find that the HMM with the proposed NMF features yields an average classification rate of 87.27%.In Table 2, the confusionmatrix for 6-degree test sequences is shown.For the 6 degree angular extent, an average classification rate of 95.62% is obtained.The improvement of the classification performance is due to the increased number of state transitions that can occur in each test sequence.
Increasing the number of HMM states and Gaussian mixtures theoretically results in a complex model that is better able to model the target signature.However, if training data is not sufficient relatively, HMM acquired by training data will fail to generalize well to test data.Figure 5 shows the test set correct recognition rate for the same 10-target identification experiment using different number of HMM states  and Gaussian mixture components .When  = 60 the best performance is achieved with  = 2 and degrades for larger or smaller  due to overtraining or insufficiently modeling.So, in the remaining experiments we retained the initial configuration of  = 60 and  = 2.
For comparison, we compare the proposed approach with several state-of-the-art methods.The base line for the SAR ATR performance comparison is the template matching method [1].For our approach, it can be seen that the correct recognition rate of 95.6% is better than 91.8% by the template matching method.However, it is noted that our method in this paper does not require target pose estimation unlike template matching.
The obtained recognition accuracy of the proposed method is also competitive with the support vector machine (SVM).In [20], on 3-class SAR ATR tasks they employed SVM on SAR targets images to give the misclassification error of 9.01%, worse than our results of 4.4% from Table 2.For SAR ATR methods based on HRRPs, a comparison can be made between the correct recognition rate of 92% presented in the work of Liao et al. which used a Relax feature extraction approach with a 10-target identification experiment [5].Another spectral feature extraction approach achieved the correct recognition rate of 83.5% for SAR ATR [6].These clearly verify the superiority of the proposed method.

Robustness to Variant Depression Angles.
The robustness of a recognition algorithm to depression angle is important to the successful application of the algorithm for real scenarios, where the test target images may have been acquired from different depression angles.We will examine the invariance to depression angle for the approach.Some statistics of the variant depression angles dataset used in this experiment are summarized in Table 3.This is a subset of the MSTAR public database on 3 different targets (2S1, BRDM2, and ZSU234) at 4 different available depression angles (15 ∘ , 17 ∘ , 30 ∘ , and 45 ∘ ).The data collected at depression angle 17 ∘ are utilized for training and other data for testing.The classification results are summarized in Table 4.As can be seen from Table 4, the approach performs robust when there is a large change in depression angle (e.g., from 15 ∘ to 30 ∘ ).However, when the change is very large, the performances of the approach decrease due to the drastic change of the target signatures.

Conclusions
This paper has presented a novel feature extraction method for SAR ATR based on HRRP sequences, which achieves excellent performance on the MSTAR database.The method characterizes the target HRRP time-frequency signature by AGR and NMF.First, the AGR is applied to each HRRP to acquire the corresponding time-frequency matrix.Then, target time-frequency features are extracted by NMF.The use of the AGR can be accurately represented by one HRRP's complex electromagnetic signature in the time-frequency domain, which gets a time-frequency matrix.And the use of NMF technique can get remarkable features extracted from the time-frequency matrix.When performing classification using HMMs, the target poses are unknown.The approach was tested under MSTAR public release datasets.Through experiments, an average correct classification rate of 95% was achieved on a 10-target classification task.It was assessed how the performance is affected by the number of HMM states and the number of Gaussian mixture

Figure 1 :
Figure 1: The procedure from a SAR image to HRRPs.

Figure 2 :
Figure 2: Comparison between T-F matrices using WVD and T-F matrices using AGR.
Gaussian elementary functions   () with an adjustable T-F center (  ,   ) and a variance

Table 4 :
Recognition results with variant depression angles dataset.The robustness to variant depression angles was also tested.In summary, several experiments and the high classification accuracies achieved by the proposed technique clearly demonstrated the potential for radar target HRRP feature extraction and SAR ATR.