A) General scheme

Scheme of the proposed classification system is illustrated in Fig. 2. The approach processes the radar signal sampled at two different sampling rates: 32 and 3.2 kHz to discriminate objects with wide (helicopters) and narrow (birds, planes) Doppler bandwidth. All information is contained in the branch of 32 kHz, however, to decrease the processing time we use a second sampling rate (3.2 kHz) with different parameters of feature extraction procedure. There is no practical need of the second ADC (3.2 kHz) as these information could be retrieved from the radar data sampled at 32 kHz followed by a down-sampler. Once the signal from radar is sampled, it passes through the m-D signature estimation, filtering, alignment, feature extraction, and classifier blocks, then the final decision (with confidence score) is made. Let us consider all of these blocks in detail.

Fig. 2. Block scheme of the proposed approach.

B) Estimation of m-D signature

The m-D signature of the moving target is estimated as the magnitude of the spectrogram of the Doppler signal:

(3)$$S\lpar w\comma \; t\rpar =\left\vert {\sum\limits_{m=0}^{M - 1} W\lpar m\rpar s\lpar m+\lpar t - 1\rpar \lpar M - L\rpar \rpar e^{ - j2\pi wm/M} } \right\vert \comma \;$$

where W(m) is the smoothing window function of length M, for instance, Hamming window can be used; w is the frequency index; L is the overlap of successive Fourier lengths, expressed in samples; the dimensionality of the spectrogram is $S \in {\rm \opf R}^{M \times Q}\comma \; Q=\left\lfloor {{\textstyle{{N - L} \over {M - L}}}} \right\rfloor $; N is a total number of samples (length of s).

C) Filtering

The noise ε together with clutter could be removed from the observation s by the “spectral subtraction” noise reduction procedure [Reference Vaseghi15]. The main idea of this approach is subtraction of an estimate of the average noise spectrum from the noisy signal spectrum. We assume that noise exhibits a static frequency profile with varied gain.

The average noise spectrum is estimated by a periodogram:

(4)$$\Upsilon \lpar w\rpar =\displaystyle{1 \over Q}\sum\limits_{t=1}^Q S^{\left( \epsilon \right) } \lpar w\comma \; t\rpar \comma \;$$

where S ^(ε) is STFT of noise only signal ε.

The gain of noise is estimated as a normalized length of the projection of the noisy signal spectrum onto noise spectrum:

(5)$$G\lpar t\rpar =\sum\limits_{w=1}^M \displaystyle{{\Upsilon \lpar w\rpar S\lpar w\comma \; t\rpar } \over {\Vert \Upsilon \Vert ^2 }}\comma \;$$

where $\Vert \Upsilon \Vert =\sqrt {\Upsilon \lpar 1\rpar ^2+\Upsilon \lpar 2\rpar ^2+\ldots+\Upsilon \lpar M\rpar ^2 } $ is a Euclidean norm of the noise pattern.

The noise can be removed from the signature as:

(6)$$U\lpar w\comma \; t\rpar =S\lpar w\comma \; t\rpar - G\lpar t\rpar \Upsilon \lpar w\rpar.$$

In order to make U non-negative:

(7)$$U\lpar{\cdot}\rpar =\left\{{\matrix{ {U\lpar{\cdot}\rpar \comma \; } & {{\rm if}\; U\lpar{\cdot}\rpar \geq \beta S\lpar{\cdot}\rpar } \cr {\beta S\lpar{\cdot}\rpar \comma \; } & {i{\rm f}\; U\lpar{\cdot}\rpar \lt \beta S\lpar{\cdot}\rpar } \cr } } \right.$$

The threshold is set to be at level β = 0.01.

D) Alignment of the m-D signatures

In order to extract behavior of m-D features the motion of the target must be compensated. This can be done by tracking the change of velocity of the target's body.

Let us assume that the target's body is a scatterer with the highest reflectivity coefficient. Therefore, the target's body will appear as the maximum at the STFT. We need to track this maximum and then compensate the motion. However, the contribution of the stationary clutter is not removed yet. For this purpose, we propose to apply a weighting function before estimating the maximum:

(8)$$h\lpar w\rpar =- \gamma e^{ - {\textstyle{{w^ 2 } \over {2\sigma ^2 }}}}\comma \;$$

where γ = 128000/F _s, σ = 7·(3200/F _s)³ and F _s is sampling frequency. It should be noticed that the weighting function h is obtained in dB scale. The parameters of the weighting function are selected empirically to maximize the classification rate on the training data. The weighting function assumes that stationary clutter contributes to near zero Doppler frequencies only, with stationary amplitude, estimated on signal-free data; therefore, it has a form of constant function with amplitude 0 dB and a gap of −γ dB at zero Doppler.

As initial points to track the target's body motion, we take 10% of samples with the highest amplitude from the spectrogram. Then, unsupervised clustering is applied to reduce the number of initial points within the neighboring area. Then for each remaining initial point, the maximum is tracked within a local window with first increasing and then decreasing time index. In such a way, we obtain a number of possible tracks. Then the track corresponding to the maximum accumulated energy is selected as the target's body track V(t).

The tracked velocity V(t) is approximated by a polynomial of order 1 to reduce the number of outliers. In such a way the track is assumed to be a simple line, meaning that motion can contain only linear acceleration. The m-D signature is shifted then: $\hat U\lpar w\comma \; t\rpar =U\lpar w+\hat V\lpar t\rpar \comma \; t\rpar $. The last procedure is done by linear interpolation. Empirically, we found that for a short observation time (decision making interval is 0.5 s) linear track is a reasonable choice for considered types of UAVs.

Finally, the m-D signature is ready for feature extraction.

E) Feature extraction

The features are based on extraction of bases of the m-D signature. These bases are orthogonal to each other and contain essential information about the rotating parts of the target. After the alignment, we assume that the spectrogram can be viewed as a low rank matrix.

First, we need to compute the correlation between different frequency components. Correlation matrix is computed as:

(9)$$\Psi \lpar w_1\comma \; w_2 \rpar =\sum\limits_{t=1}^Q \hat U\lpar w_1\comma \; t\rpar \hat U\lpar w_2\comma \; t\rpar.$$

Second step is to estimate eigenpairs {v _r, λ_r}, where v _r is r-th eigenvector and λ_r corresponds r-th eigenvalue, such that λ₁ > λ₂ > …λ_r. Each eigenpair satisfies the following equality:

(10)$$\Psi v_r=\lambda _r v_r.$$

Eigenvectors are orthogonal and unique forming the basis functions of the signal's spectrum. Eigendecomposition of Ψ or singular value decomposition of $\hat U$ can be used for estimation of eigenpairs. The steps (9) and (10) are similar to principal component analysis (PCA) with only small difference [Reference Duda, Hart and Stork16]. For PCA, the data must be mean centered before calculating the correlation matrix. In our case, the first (and the most important) eigenvector corresponds to the mean vector of the spectrogram.

Next, the Fourier transform of the eigenvectors is computed to obtain features with strong “energy compaction” property, i.e. the features where most of the signal information is concentrated. Typically, the signal information is contained in just a few low-frequency components. Owing to this property, we can calculate only the I first (low-frequency) coefficients to represent the data:

(11)$$y_r=\bigcup\limits_{l=1}^I \left\{{\left\vert {\sum\limits_{p=1}^M v_r \lpar p\rpar e^{ - j2\pi pl/M} } \right\vert } \right\}.$$

Finally, the feature set is obtained by combining the first Z transformed eigenvectors and Z eigenvalues:

(12)$$F=\bigcup\limits_{r=1}^Z \lcub y_r\comma \; \lambda _r \rcub.$$

F) Classifier

In this paper, three different types of classifiers are considered in order to evaluate the robustness of the proposed features. The first two classifiers are the linear and non-linear support vector machines (SVMs) belonging to non-probabilistic classifiers, and the third one is the Naive Bayes classifier (NBC) belonging to probabilistic linear classifiers.

1) NAIVE BAYES CLASSIFIER

Naive Bayes classifier decides in favor of the class with the highest posterior probability, assuming that the features associated with each class are conditionally independent.

The decision is based on the maximum of class conditional probability:

(13)$$\hat C=\mathop {\max }\limits_C P\lpar C\vert F\rpar \comma \;$$

where F is a feature vector.

Considering the Bayes' theorem and assuming the conditionally independent properties of features, the decision can be expressed as:

(14)$$\hat C=\mathop {\max }\limits_C P\lpar C\rpar \prod\limits_{n=1}^N P\lpar F_n \vert C\rpar \comma \;$$

where P(C) is a prior probability of class C; P(F _n|C) is the likelihood function; F _n is the feature n.

Commonly, the prior probability is unknown and the equation is solved by maximizing the product of likelihood functions. Therefore, the decision rule is called the maximum likelihood (ML). Finally, the ML rule can be written as:

(15)$$\hat C=\mathop {\max }\limits_C \prod\limits_{n=1}^N P\lpar F_n \vert C\rpar.$$

According to the method of estimating P(F _n|C) the Naive Bayes classifier can be considered as normal or kernel-based. In the first case, the likelihood function is approximated from the Gaussian distribution. In the second case, the likelihood function is approximated by a kernel smoothing density function.

2) SUPPORT VECTOR MACHINE (SVM)

Support vector machines belong to the non-probabilistic kind of binary classifiers. In this paper, we use the version of a SVM with soft margins proposed in [Reference Cortes and Vapnik17]. Assuming the linearly separable classes, a hyperplane with maximal margins can be found to separate the classes with minimum error. The SVM can be extended to the multiple hypothesis classifier combining various SVMs using the one-against-one approach.

Assuming the training set of instance-label pairs $R_i={\rm \lcub }F^{\lpar i\rpar }\comma \; c_i {\mkern 1mu} \vert {\mkern 1mu} F^{\lpar i\rpar } \in {\rm \open R}^n\semicolon \; c_i \in {\rm \lcub }1\comma \; - 1{\rm \rcub \rcub }$ is given, the support vectors are found by the solution of the following optimization problem:

(16)$$\matrix{ &{\displaystyle{1 \over 2}w^T w+S\sum\limits_{i=1}^I \epsilon _i } \cr \mathop {{\rm subject}\; {\rm to}}\limits^{\mathop {\min }\limits_{w\comma b\comma \epsilon } } &\!\!\!{c_i \lpar w^T \phi \lpar F^{\lpar i\rpar } \rpar +b\rpar \geq 1 - \epsilon _i\comma \; } \cr &{\epsilon _i \geq 0\comma \; }}$$

where ε_i is the degree of misclassification of the data F ⁽ⁱ⁾; w and b describe the separation hyperplane; S is the penalty parameter of the error term; ϕ(·) is the function which maps the data F into a higher dimensional space.

Actually, the SVM is a linear classifier. However, not all data are linearly separable in the data space. Therefore, the function ϕ(·) is used for data mapping into a high dimensional space, where data are linearly separable with the minimum error (soft margin).

In this paper, we use a radial basis function: K(F ⁽ⁱ⁾, F ^(j)) = exp(−γ‖F ⁽ⁱ⁾ − F ^(j)‖Reference Koolhoven²), γ > 0. Thus, the function ϕ(·) is selected to satisfy K(F ⁽ⁱ⁾, F ^(j)) = ϕ(F ⁽ⁱ⁾)^Tϕ(F ^(j)).

A number of parameters are required to characterize the SVM classifier. The hyperplane that separates classes is described by a few support vectors (few feature vectors, F ⁽ⁱ⁾). These selected support vectors are based on the training set. The parameter γ of mapping function ϕ(·), and the penalty parameter S are estimated using cross-validation procedure for the training set using grid search. The considered parameters of S are equal to {exp(−5,…,5)}, and γ are {1/2exp(2σ)|σ ∈ {− 4,…,4}}. The parameters of the SVM are selected empirically to maximize the classification rate on the training data.

A) Data collection

The proposed ATR system is evaluated on real radar measurements. The radar data have been collected with a continuous-wave radar operating in X-band at radio frequency of 9.5 GHz. These radar measurements have been performed within the framework of D-RACE, the Dutch Radar Centre of Expertise, a strategic alliance of Thales Nederland B.V. and TNO. The measured data have been made available for the present study.

The UAVs considered for classification by the proposed technique are listed with their characteristics in Table 1. Different types of the UAVs are considered: two planes, a quadrocopter and three helicopters. For comparison, stationary rotating rotors with different number of blades as well as diverse species of birds are considered for classification. At least 30 s were recorded of each class. The data were recorded at open space with the distance to target less than 30 m. The aspect angle, distance and velocity were varying as within the real situation.

Table 1. UAVs used for experimental verification of the proposed ATR system. Type ‘P’ is for Plane, ‘H’ is for Helicopter, ‘Q’ for quadrocopter, ‘B’ for birds, ‘S’ for stationary rotors.

For classification we formulate two problems:

11-classes problem where each class corresponds to different model of the UAV, and the goal is to determine the correct type as well as the model.

5-classes problem where each class consists of the data from the same type of aerial object, and the goal is to determine only the type of the flying object.

The radar signal is divided onto a number of non-overlapping segments of fixed length, by processing these segments a decision about the class label is made. The length of one segment is set to 0.5 s. To compute the STFT we use a sliding Hamming window of length M = 128 samples = 4 ms with overlapping of L = 90 corresponding to 0.9 × 4 = 3.6 msFootnote ¹. Finally, the features extracted for the two different sample rates are concatenated to one feature vector.

Examples of m-D signatures of the considered classes are illustrated in Fig. 3. The classes can be separated even visually. The features we extract are based on eigenvectors of the signature and the first two of them are also shown in Fig. 3.

Fig. 3. Examples of filtered and aligned micro-Doppler (m-D) signatures of 5 classes (first row); eigenvectors corresponding to the largest eigenvalue (second row) and the second largest eigenvalue (third row) of ψ for corresponding signatures. These eigenvectors are basis functions of the m-D signatures.

B) Performance analysis

To classify the type of UAV, we will use the machine learning approach. All available segments are divided into training and testing sets without overlapping. The decision about class membership is made by a predefined system of rules. The process of defining the system of rules is called learning; it is performed on the training set. The number of correctly classified segments determines the probability of correct classification as ratio to the total number of segments from the testing set.

The K = 10 cross-validation technique is applied to obtain robust classification rates. The initial data under analysis are split into K subsets of the same length, and K−1 subsets are used as a training dataset, and the remaining one as a testing dataset. The cross-validation process is repeated K−1 times (K−1-folds) with each of the K subsets used as a testing dataset. The K results from the folds are averaged to evaluate a single estimation. The cross-validation is applied in such a way that training and testing sets contain data from different trials (experiments).

First, we want to study the influence of the feature extraction parameters on the classification result. Therefore, the number of eigenpairs in (12) was varied as Z = 1,2,…,14 and the number of calculated Fourier coefficients in (11) was varied as I = 1,2,…,40. Total number of 560 classification results were obtained using ten-fold cross-validation for all possible pairs of Z and I. Then, the classification rates were averaged for all values of Z and I, to show the dependence on these variables in Fig. 4. Analyzing the dependences, we conclude the following:

• – With increasing the number of Fourier coefficients, SVM classifiers show increase in the performance;

• – For Naive Bayes classifier, the performance reaches the pick with 16 Fourier coefficients, with larger number the performance degrades due to the “curse of dimensionality” phenomenon;

• – For non-linear classifiers, such as non-linear SVM, with increasing the number of eigenpairs, the performance degrades; the best result is obtained for only one eigenpair extracted.

• – For linear classifiers, such as linear SVM and Naive Bayes classifier, the best performance is reached when several eigenpairs are extracted. Such for linear SVM, four eigenpairs provide the best performance, for Naive Bayes classifier five eigenpairs are selected.

According to the analysis we selected the following parameters:

• – Non-linear SVM: Z = 1, I = 29;

• – Linear SVM: Z = 4, I = 33;

• – Naive Bayes classifier: Z = 5, I = 16;

Computed probabilities of correct classification for different classifiers are shown in Table 2. The results are shown for different sampling frequencies of the signal: 3.2, 32 kHz and when two of them are applied. The results in Table 2 allow us to conclude the following:

• – For all classifiers the best results are obtained when two sampling frequencies are used;

• – The benefit of processing the signal sampled at two sampling frequencies, compared to the best result out of one sampling frequency, is: SVM linear 3.56%; SVM non-linear 5.85%; Naive Bayes 0.14%;

• – The best result is obtained with SVM nonlinear classifier: 92.3%.

Confusion matrix for non-linear SVM with applied ten-fold cross-validation is shown in Table 3 to estimate the intra class distribution of classification rates. It can be seen that all targets, except classes 1, 2 and 7, 8 are classified with probability of correct classification higher than 95%. The first two classes have more errors because blades are small and their m-D signatures are similar. The results for classes 5–8 show the capability to distinguish different number of rotors. The class #7, three stationary rotors, shows misclassification at level of 15% with class of four stationary rotors and helicopters class.

Table 2. Probability of correct classification obtained for 11-classes problem using ten-fold cross validation.

Table 3. Confusion matrix (in percentages) for linear SVM classifier, 11-classes problem.

To show the flexibility of the proposed ATR scheme to classify different types of the UAVs, the following classes were removed from the training set: 1, 5, 6, 7, 9, 10. In this way for each type of the UAVs (plane, birds, stationery rotors, quadrocopter, and helicopter) only a single class is used for training. The removed classes are then used to test the system. The results of classification by linear SVM classifier are listed in the Table 4. The class of stationary rotating rotor, and the class of helicopters are classified correctly with 100%; the data from class 1 are classified correctly as a class of plane with probability of 86%. We can claim that proposed ATR scheme is robust to variations inside the class, and therefore with high probability, other UAVs with similar flying concept will be classified correctly.

Table 4. Classification of UAVs hidden from training procedure, 11-classes problem.

For ATR system, it is important to reject an unknown target which class significantly differs from the ones used for training. To reject the unknown target a confidence threshold can be set. By confidence threshold, we mean a minimum value of the class conditional probability (posterior probability) for which the decision is assumed to be true, otherwise the target is marked as unknown. Even for non-probabilistic classifiers, such as SVM, it is possible to estimate the approximated posterior probability based on distances to the hyperplanes [Reference Wu, Lin and Weng18]. If we set the confidence threshold to 0.85 then with high probability targets from classes not used for training will be rejected. To demonstrate this Table 5 is shown. The following classes were removed from training set: [Reference Toscano1, Reference Koolhoven2, Reference Clemente, Balleri, Woodbridge and Soraghan5, Reference Chen6, Reference de Wit, Harmanny and Prémel-Cabic7, Reference Chen and Ling9, Reference Yoon, Kim and Kim10]. As classes 1 and 2 are removed, the classifier has no data to train the class of planes, and these classes are assumed to be from class of unknown target. According to Table 3, 83% of the data of the plane's class is classified as unknown; also, 16% of data from quadrocopter's class and 9% of the class of stationary rotating rotor are misclassified as unknown target (Table 5).

Table 5. Classification of the unseen for training classes, 5-classes problem.

For all experiments above we used half of a second of the observation time to classify the target. Next, we show what a performance can be achieved when shorter interval is used for decision making. The training set is constructed in the same manner as before (using 0.5 s as segment's duration); however, for training set, the features are extracted from shorter segments. This shows the robustness of the feature extraction procedure to the duration of the segment.

Figure 5 demonstrates the dependence of the classification rates on the dwell time for 11-classes problem. The following conclusions can be pointed out:

• – Both SVM classifiers demonstrate stable result with dwell time longer than 0.25 s;

• – Naive Bayes classifier provides a decrease of 4% with dwell time of 0.25 s;

• – With dwell time shorter than 0.25 s, the performance drops significantly with decrease of the dwell time. At dwell time 0.05 s the following result is achieved: linear SVM 75%, non-linear SVM 82%, Naive Bayes classifier 65%.

Noise can corrupt the signal and degrade the classification performance. We will add white Gaussian noise to experiments from test dataset and examine the performance of the system. To do this, we assume that collected radar signal contains no noise and the power of the raw data will be used to calculate Signal to Noise Ratios (SNRs). Additive white Gaussian noise is generated as variable ε(n) in (2). Figure 6 demonstrates the dependence of classification rate on SNR, for 11-classes problem. We can conclude the following:

• – The proposed approach together with SVM classifier demonstrates robustness to the SNR when it is higher than 30 dB;

• – Performance drops to half with SNR = 18 dB, and reaches probability of random guess at 0 dB;

• – Naive Bayes classifier together with the proposed approach is more sensitive to the SNR; drop to the half of the performance is obtained with 25 dB and than rapidly decreases.

Fig. 4. Probability of correct classification for 11-classes problem depending on: (a) the number of Fourier coefficients in (11); (b) depending on the number of eigenpairs in (12).

Fig. 5. Classification rate for 11-class problem depending on the dwell time.

Fig. 6. Classification rate depending on SNR, 11-classes problem.