Intrinsic synchrosqueezing analysis on micro ‐ Doppler features for small unmanned aerial vehicles identification with dual ‐ channel radar

The micro ‐ Doppler (m ‐ D) effect depends on the rotation of rotor blades in addition to the translation of the platform. Thus it is a characteristic for identifying small unmanned aerial vehicles (UAVs). However, compared with the Doppler signal induced by the translation of the platform, the m ‐ D signal is weak. In this article, a highly localised data ‐ association method, intrinsic synchrosqueezing analysis (ISA), is proposed for estimating m ‐ D characteristics from the returned signal of small UAVs with a dual ‐ channel radar. Employing synchrosqueezing transform on intrinsic mode functions derived from noise ‐ assisted multivariate empirical mode decomposition, the proposed ISA method separates the Doppler signal and enables denoising and sharpening time ‐ frequency representation of the m ‐ D signal. Simulation results confirm the theoretical analysis, showing the feasibility of estimating m ‐ D features in a noisy environment. Applications on field data illustrate brighter prospects for identifying small UAVs.


| INTRODUCTION
Because of their low price, ready availability, and suitable operational complexity, small unmanned aerial vehicles (UAVs) are rapidly proliferating [1,2]. Most commercial small UAVs are used for recreational purposes such as aerial photography, but they are also exploited for criminal and antisocial activities, such as intruding reconnaissance and the transportation of explosives. Emerging threats in the military or civilian field highlight the need to monitor small UAVs. To this end, research on reliable and robust measurements capable of identifying small UAVs is imperative.
Although various sensors such as laser sensors and visible light cameras are used in crowded urban and littoral environments, radar is more suitable because it can detect a target at night and in all weather conditions. Furthermore, dual-channel frequency modulated continuous wave (FMCW) radar can provide excellent measurement performance for the distance, angle, and speed of moving targets in portable and low-power applications. However, low altitude, slow speed, and small radar cross-section (RCS) are factors complicating the identification of small UAVs. Traditional detection techniques designed for FMCW radar, such as moving target identification, cannot acquire Doppler frequency when the UAV is in hover. Moreover, naturally present targets flying at a similar speed, such as birds [2], may potentially saturate signal processing for the Doppler frequency. Instead, the rotation of rotor blades captured by micro-Doppler (m-D) analysis contains unique kinematic features of small UAVs, which could be deployed for identification.
The m-D effect in radar signal processing, induced by rotation, vibration, or procession, was originally investigated by Chen [3]. Consequently, time-frequency (TF) distribution is introduced for an analysis of m-D features [4,5]. Nevertheless, conventional TF techniques, such as short time Fourier transform (STFT), share the same limitation, known as the uncertainty principle, stipulating that one cannot localise a signal with arbitrary precision in both time and frequency [6][7][8]. Moreover, owing to multiple components of the returned signal, the rotation parameters cannot be directly estimated from TF distribution with reasonably high precision.
To this end, synchrosqueezing transform (SST) was first investigated in Daubechies et al. [9], originally based on the continuous wavelet transform (CWT) [10,11], to enhance the localization property by rearranging the energies of TF coefficients in the data. SST can overcome the limitations of TF analysis by producing high-resolution ridges without spreading the energy of each term. Subsequently, the separation, reconstruction, and denoising of the multivariate signal are effectively proposed by multivariate synchrosqueezing wavelet transform (M-SWT) [12]. However, because CWT is restricted by the resolution and produces additional noise-wavelet coefficients that correspond to undesired frequency components, the performance of those methods degrades.
To mitigate these problems caused by CWT, the multivariate STFT-based SST algorithm (M-SST) [13] exploits various approximations of varying instantaneous frequencies for amplitude and phase with higher accuracy, which could obtain sharper representation along with better mode reconstruction. Nevertheless, because the returned signal of small UAVs is composed of complicated components, multiple energy ridges exist in the TF plane, which poses challenges for estimating m-D characteristics.
Theoretically, multivariate empirical mode decomposition (MEMD) provides localised TF representation by the adaptive decomposition of a multivariate signal into a finite series of intrinsic mode functions (IMFs) without making an assumption about its linearity and stationarity features [14]. Each IMF of the same channel exhibits detailed information about the channel data in the certain time scale or frequency band. However, IMFs with the same index contain segments of different scales from different channel signals, which is called the mode-mixing problem [1]. The mode-mixing problem occurs in two aspects. First, the number of IMFs for each channel is different, mainly because each channel is treated independently, so the IMFs of multiple channels fail to be arranged into the same index according to the frequency responses. Second, based on the spatial relationship between the radar and UAV, the platform would shelter some electromagnetic wave that should been reflected by the rotor blades. In this case, the rotation signal is discontinuous in certain periods. Theoretically, potential intermittency of the rotation signal would also induce the mode-mixing problem [14].
Within such a context, we aim to decompose the complicated returned signal and extract rotation features for identifying small UAVs with dual-channel radar. In this article, intrinsic synchrosqueezing analysis (ISA) is proposed to estimate m-D parameters from a bivariate returned signal, which could be extended to multivariate data. To reduce the modemixing problems, this article uses noise-assisted MEMD (NA-MEMD) [15] to decompose the returned signal of small UAVs. According to the phase spectrum of various IMFs after NA-MEMD, we select the corresponding IMFs within the rotation frequency band as the input signal for the STFTbased SST (SSST). Simulations demonstrate that the proposed ISA method could enhance the TF resolution and separate monocomponent m-D signal for the analysis of rotation characteristics with higher precision compared with other prevailing techniques. Field experiments with dual-channel continuous wave radar are performed to render this technique promising for identifying small UAVs. The main contribution of the work is twofold.
First, the input data of NA-MEMD consists of the returned signal and additional independent white Gaussian noise (WGN) in one or more separate channels to align IMFs in accordance with the quasidyadic filter bank structure, in which each component carries a frequency subband of the returned signal. In this case, the proposed ISA method could separate the Doppler signal and decrease the contamination of additive noise components.
Second, the phase spectrum is deployed to select IMFs corresponding to the rotation frequency band. Employing the SSST on the selected IMFs, the purely energy-concentrated monocomponent TF representation of rotation signal is acquired by preserving the high-frequency components with the rotation property. Therefore, m-D interference owing to the lower-frequency band induced by the vibration of small UAVs could be separated.
The rest of this article is organised as follows: Section II defines the signal model. Next, in Section III, we give the details of the proposed ISA method. Section IV compares prevailing theories with the proposed ISA method via the results of simulation and field experiments. Finally, in Section V, we provide concluding remarks.

| SIGNAL MODEL
m-D properties induced by the rotation of rotor blades are specific to a given target type. Besides the m-D signal, the returned signal contains a Doppler signal induced by the translation of the platform as well as the noise. Figure 1 [1] shows the geometry of the radar and UAV. Generally, the radar is stationary and located at Q. The origin O in the Cartesian coordinate system O-XYZ represents the bottom of the UAV; OZ is the perpendicular axis. Moreover, the stick on the top of the UAV indicates the rotor blade with the uniform density, which has the rotation angular velocity ω r around the centre O 1 .
For simplification, we suppose that the line of sight (LOS) of radar is on the YOZ plane and α(t) represents the angle between the LOS and Z-axis. The orientation vector of the LOS is expressed as [1]: Suppose the translation path of the UAV is along the LOS direction. We introduce the reference coordinate system O 0 -X 0 Y 0 Z 0 with the same translation velocity v. Generally, the detection range meets the far field condition; to this end, we could assume that αðtÞ ≐ α; L → ðtÞ ≐ L → are approximately constant in a small period t. Then, the orientation vector of edge points is defined: where the initial angle between the rotor blade and Y-axis is ϕ. At the instant time t, the edge point P moves to P 0 , whose location is given as P Therefore, the instantaneous range for QP 0 is written as [16]: ‖ represents distance QO 1 and φ is the initial phase. A 0 indicates the amplitude of the micromotion, which is the function of the radius of the turn and angle α.
Generally, the hover or other movements of small UAVs are always accompanied with undesired vibration. As a typical micromotion, the vibration also introduces the sinusoidal frequency modulation with a lower frequency than the rotation. Therefore, m-D data of the returned signal are composed of multiple components with different micromotion frequencies.
In this regard, if the radar transmits a continuous wave with a carrier frequency f 0 [17], the baseband of the returned signal is the summation of Doppler signal S D (t), m-D signal S m (t) and noise η(t), given as: Both S D (t) and S m (t) have similar translation components V(t) in the phase term, where c denotes the speed of the electromagnetic wave propagation. K is a constant. σ 1 represents the scattering coefficients of the platform. Instead, the reflectivity function of the kth micromotion scattering centre is given as σ k . Moreover, A k indicates the m-D amplitude for the corresponding micromotion frequency f k .
To obtain the representation for the returned signal of a dual-channel radar, we could construct a vector for the generalisation of a multivariate analytic signal with N channels: where S n (t) indicates the returned signal of the nth channel n = 1, L, N.

| PROPOSED ISA METHOD
Based on the multivariate returned signal model of small UAVs, the core for extracting m-D characteristics is to separate the Doppler signal and the noise, as well as m-D interference. In this article, we propose the ISA method to estimate the rotation frequency for identifying small UAVs. This method is composed of three stages. First, the NA-MEMD algorithm is employed on the multivariate returned signal to discard Doppler information according to Doppler frequency bands of corresponding IMFs, but also to reduce the mode-mixing problem within the IMF. Second, the phase spectrum is developed to select IMFs exhibiting rotation. Finally, we use the SSST method to acquire the purely energy-concentrated monocomponent of TF representation to estimate the m-D frequency.

| Noise-assisted multivariate empirical mode decomposition algorithm
MEMD is a data-adaptive method of decomposing a multivariate signal into a finite set of multichannel IMFs according to the levels of their local oscillation or frequency in a physical domain. Each channel IMF exhibits detailed information in a certain frequency band, but IMFs with the same index of all channels contain different frequency scales derived from different channel signals. On the other hand, because the instantaneous frequency of m-D signal is time-variant, m-D information may be exhibited in different IMFs in certain periods. In this case, rotation features are barely distinguishable owing to the severe mode-mixing problem. NA-MEMD is an advanced extension of the MEMD algorithm. With the assistance of noise, multivariate input data, consisting of the original signal and added independent WGN [16], enable NA-MEMD to align IMFs in accordance with the quasidyadic filter bank structure, which could reduce the mode-mixing problem within the extracted IMFs [16]. The details of NA-MEMD are outlined in Algorithm 1 [18]. Algorithm 1 Noise-assisted multivariate empirical mode decomposition Input: The input multivariate (n channels) vector S n�d , the uncorrelated WGN time series (m channels) W m�d with the same length d.

Procedure:
1) Add the WGN signal W m�d to the input multivariate signal S n�d , obtaining multivariate signal X (n+m)�d with n + m channels.
2) Process the resulting signal X (n+m)�d using the MEMD algorithm to calculate multivariate IMFs I (n+m)�d�y , where y represents the number of IMFs for each channel data.
3) From the acquired n + m variate IMFs, discard m channels belonging to the WGN and select a set of n channels IMFs corresponding to the original signal. Output: The IMFsÎ n�d�y : In the case of detecting small UAVs with a single channel radar, because the RCS ratio of rotor blades to the platform is small, the literature indicates that the Doppler signal would be extracted in the first few IMFs by the empirical mode decomposition [19], whereas the m-D features are exhibited in subsequent IMFs. Heuristically, owing to the quasidyadic filter property, the extra noise channels of NA-MEMD serve as a reference in the TF plane and enable a more stable estimate of IMFs in the signal-channel space, giving more accurate TF analysis. Applying this method to the returned signal of the dual-channel radar, the behaviour of the Doppler signal and m-D components could be distinctly exhibited in different IMFs with higher precision and less unwanted mode-mixing.

| Selection of intrinsic mode functions
As illustrated in section Ⅱ, the m-D signal is induced by the rotation of rotor blades and other micromotions of small UAVs, such as vibration. Although the mode-mixing problem of extracted IMFs is reduced by using the NA-MEMD method, m-D interference data, whose instantaneous frequency is time-variant and may be multiplied several times micromotion frequency, would mingle with rotation components in the corresponding IMFs. For the IMFs showing the Doppler property, the Fourier spectrum is employed to acquire the estimation of velocity or range directly. Nevertheless, the periodically sinusoidal modulation is the exclusive property for the instantaneous frequency (IF) of the m-D signal. Because the sinusoidal frequency modulation in the phase term is equivalent to the power series with infinite order terms, a series of multiple or difference frequency peaks of mixed m-D signal emerge in the Fourier spectrum. Therefore, the selection of the corresponding IMF with rotation features is intractable with Fourier transform. Analysing the signal model of the m-D signal, the phase term has the similar form as Doppler modulation, which is directly relevant to the spectrum of micromotion [1]. The similarity inspires us to discuss the phase spectrum for selecting IMFs related to the rotation signal.
Theoretically, the ith (0 < i ≤ y) IMF ofÎ n�d�y for each channel shares a similar frequency property [18]. In this case, M p,i is defined as the ith IMF of the pth channel for selecting IMFs with rotation. The phase spectrum is acquired as [2,20]: where PF(•) represents the phase Fourier transform and F(•) is the Fourier transform. The principle of phase Fourier transform is equivalent to the theory of the cyclic autocorrelation functions of the phase term in Zhao et al. [2]. The first term is the Fourier transform of phase terms. Moreover, the phase shift ambiguity is revised by PR.
The average rotation frequency of rotor blades is between 100 and 250 Hz. When the UAV is hovering, the rotation frequencies of four rotor blades are equal and relatively stable. However, if the small UAV is in other motions, the rotation frequencies of four rotor blades correspondingly shift to generate forward momentum. If the highest frequency peak of phase spectrum is within this bound, the IMF is regarded as the R-IMF, which exhibits rotation features. Although the mode-mixing problem within the IMF is reduced by the NA-MEMD algorithm, it also exists between the selected R-IMFs, owing to the periodically sinusoidal modulation of the m-D signal.

| The time-frequency representation with SIFT algorithm
The CWT-based SST (CSST) technique, given by a wavelet transform, could enhance the TF representation by reassigning coefficients using an estimate of IF. However, the m-D signal not composed of purely harmonic modes could not be analysed effectively by these techniques [20]. Because of the modemixing problem among R-IMFs, the TF distributions contain multiple TF ridges, which are infeasible for estimating rotation parameters. The main reason is that the CSST is limited to the processing scale variable of the CWT. As another significant variable, time shift should also be investigated, which would cause SST results to be blurred.
To mitigate these problems, we propose the SSST method to generate highly localised TF representations of R-IMFs. Compared with CSST, the SSST method could eliminate false ridges in the TF plane and extract the monocomponent harmonic ridge of the rotation signal.
If we select q IMFs with rotation fromÎ n�d�y , the R-IMF is given as: where CI r ,R 1 ≤ r ≤ R q denotes the rth IMF ofÎ n�d�y with size n � d.
To unify the energy of the R-IMF, the combined R-IMF is computed as: where SI j (t) indicates the jth (0 < j ≤ n) channel data of the combined R-IMF. Employing a window function to localise the signal in time, the STFT of each channel data of the combined R-IMF SI n�l is given by: The instantaneous frequency for each frequency index is calculated as: w s ðj; τ; ηÞ ¼ −j · STðj; τ; ηÞ −1 ∂STðj; τ; ηÞ ∂η ð10Þ Each point (j,τ,η) is reassigned to (j,w s (j,τ,η),η) for synchrosqueezing. Then, SST coefficients are carried out by inverting coefficients ST (j,τ,η) along the instantaneous frequency w s (j,τ,η): where δ(g) is the delta function.
To align SST coefficients of the TF distribution, the partitioning technique proposed in Ahrabian and Mandic [12] is used to separate the TF plane into 2 k frequency bands with equal width, depicted as [20]: where k = 0,1,L, K is the level of frequency bands, and l = 0,1, L,2 k −1 denotes the index for the k frequency band. As depicted in Ahrabian and Mandic [12] about the multivariate TF partitioning algorithm, a set of separated monocomponent signals could be distinguished by splitting a larger frequency band into smaller frequency subbands with the multivariate bandwidth in conjunction with frequency partitioning.
To calculate the multivariate bandwidth, a set of features should first be established. The joint spectrum of multivariate signal is given according to: where FTðwÞ ¼ F P w∈w k;l T s ðj; w; τÞ ! denotes the Fourier transform of the joint coefficients. E p (k,l) indicates the energy of joint spectrum [12]. Then, the joint bandwidth is obtained as [20]: Bðk; lÞ ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where � w ¼ 1 2π ∫ ∞ 0 wS k;l ðwÞdw represents the joint mean phase frequency [20].
Finally, those defined parameters could judge the instantaneous frequency and amplitude for each identified frequency band, after which the monocomponent TF representation of the m-D signal is achieved. Moreover, noise with a similar rotation frequency in the selected IMFs is eliminated according to the joint bandwidth of the SST coefficients.

| RESULTS AND DISCUSSION
In this section, the proposed ISA method is compared with other prevailing approaches to show its superiority for the concentration of the m-D ridge in TF representation for both synthetic and real-world data.

| Synthetic signal
According to the separation of the multivariate returned signal of small UAVs illustrated in section Ⅲ, Doppler information is always distilled in the first few IMFs when applying the NA-MEMD method, whereas the following IMFs show the property of the m-D signal. In this case, we omit the impact of the Doppler signal and investigate the separation of the mixed m-D signal. In the signal model of Section Ⅱ, we assume that the primary m-D interference is induced by the vibration of the platform. Generally, the vibration frequency is less than the rotation frequency of the rotor blades. Therefore, the signal of interest is the rotation signal with a high micromotion frequency, whereas the m-D signal with a lower micromotion frequency is m-D interference. In the simulation, synthetic m-D data consist of a bivariate sinusoidal oscillation containing the signal of interest and low-frequency interference: where x m1 (t) and x m2 (t) are assumed to be the sum of the wanted signal with a micromotion frequency of 15 Hz and the m-D interference with the lower frequency of 2 Hz. Furthermore, the second m-D signal x m2 (t) contains the constant frequency deviation of 0.1 Hz between two channels. n 1 (t) and n 2 (t) denote the white Gaussian noise with a signal-to-noise ratio of 10 dB of two channels. First, the M-SWT algorithm [12] is applied to the mixed signal, as shown in Figure 2a. The ridge of the strongest energy does not indicate the frequency of the wanted signal. Similarly, the result of the M-SST in Figure 2b does not exhibit the separated m-D signals.
With the NA-MEMD method, the fourth to seventh IMFs are selected according to the frequency band of the phase spectrum. The TF distributions of CSST and the proposed ISA method are compared in Figure 3a,b, respectively. The result of the CSST method exhibits fluctuant ridges of the m-D components, but the energy spreads out around the TF representations. Instead, in the proposed ISA method, the energy of the high-frequency m-D signal is enhanced and the fluctuations in the ridge are removed.

| Application to the identification of small unmanned aerial vehicles
The experiment of identifying small UAVs is carried on the roof of a building with a dual-channel radar system, as shown in Figure 4a. The carrier frequency is 24 GHz and the transmitted power is 1.5 W. As shown in Figure 4b, the small UAV loaded with four rotor blades with a turning radius of 0.05 m, is receding from the radar with a velocity of v ≈ 1 m/s. The average rotation frequency is about 130 Hz with the top gear of the controller. The initial distance from the UAV to the radar is R 0 = 20 m.
The univariate returned signal and the STFT are depicted in Figure 5a and 5b, respectively. Because the m-D signal has a lower energy than the Doppler signal and noise, m-D features cannot be extracted directly from the STFT. Moreover, multiple Next, Figure 6a,b show TF distributions of M-SWT and M-SST corresponding to the two channel signals. The results of the two methods are distinctly superior to the STFT in both temporal and frequency resolutions. As expected, M-SST leads to a relatively sharper and purer TF representation compared with that in M-SWT. Unfortunately, the ridge with the strongest energy displays the features of frequency multiplication in the results of the two methods, which deteriorates the estimation of rotation frequency.
With the ability to decompose a noisy multivariate signal and address the mode-mixing problem, NA-MEMD could extract information from high-frequency noise as well as Doppler information induced by the translation in the first few IMFs. In this experiment, the phase spectra of the third to the sixth IMFs imply that they contain the information about the rotation signal. Figure 7a displays the TF representation with the CSST method. Because the coefficients of the CWT yield an overlap between bandwidths of the input signal, multiple ridges of low-frequency components contaminate the TF distribution of the rotation signal. Instead, as shown in Figure 7b for the proposed ISA method, the contamination of multiple false alarms is eliminated. The monocomponent ridge of the wanted m-D signal is extracted with frequencyf r ¼ 132:1Hz, where the estimate error is 1.6%.
Next, the proposed ISA method is applied to the returned signal of the motion of vertical ascent. In this case, the translation velocity in the vertical direction is about v ≈ 2 m/s. The power setting of the remote control is set to the second gear, in which the average frequency of the rotor blades is about 110 Hz. Because the UAV rises almost directly over the radar, the platform would shelter most of the electromagnetic wave that should been reflected by the rotor blades. Compared with the first experiment for horizontal movement, the rotation signal is weaker and more discontinuous in this experiment. Figure 8a,b show the univariate returned signal and the STFT, respectively. The periodic m-D features are not uniform and cannot be extracted in the TF distribution.
Moreover, the results of the M-SWT and M-SST methods exhibit severe contamination, as shown in Figure 9a,b, respectively. The concentrated TF distribution is in the lowfrequency band, but the frequency curve with the maximum energy is submerged in this area.
Because the IMFs are generated according to various frequency bands of the signal, R-IMF contains only the frequency band related to the rotation frequency. As displayed in Figure 10a, the CSST method yields the blurred TF distributions. Instead, the monocomponent ridge in Figure 10b indicates that the ISA method could eliminate the contamination of false alarms for the estimation of rotation frequency. The frequency of the ridge isf r ¼ 117:6Hz, where the estimate error is 6.9%.

| Quality of representation
To evaluate the performance of these techniques regarding the TF concentration, we employ the approach in Pham and Meignen [6] to measure the energy concentration by computing the normalized energy associated with the first coefficients with the largest amplitude: the faster the growth of this energy toward 1, the sharper the representation [6]. Supposing the aligned elements along different dimensions are arranged in ascending order, the normalized energy is written as NE = A/B. A is acquired as the cumulative sums for the square of the elements of the aligned SST coefficients, and B represents the sums for the square of the elements of the aligned SST coefficients. Moreover, M denotes the number of coefficients, which indicates the locations of the elements in the aligned SST coefficients. For example, when M = 1, it denotes the Kth element of the aligned SST coefficients. Therefore, M should be less than the number of rows of the SST coefficients.
In Figure 11, the normalized energy of the mixed m-D signal x m (t) under the signal-to-noise ratio of 5 dB is investigated corresponding to the reassignment of the previous four methods with respect to the number of coefficients. The first remark is that both curves of M-SST and M-SWT rise to 1 with the same coefficients, which means they have a similar energy concentration of TF distribution for the mixed signal. In comparison, the energy of the proposed ISA method increases to 1 with the least number of coefficients, because the selected IMFs after the NA-MEMD method contain fewer nonlinear frequency modulations related to the signal of interest. On the other hand, although the frequency band of selected IMFs is narrowed for the CSST method, the degrees of expended energy between these frequency bands are higher than that of the M-SWT method. Therefore, the curve of CSST exhibits the slowest increase in energy to reach a stable value. On the other hand, the error rates for the estimated frequency of four methods with simulated data are compared in Table 1. The m-D frequency is estimated by calculating the average frequency of all sampling bins with the strongest energy. In general, the M-SWT and M-SST methods acquire the frequency property of the returned signal, so the ridges would spread over the TF plane. In this case, the estimation results of the M-SWT and M-SST methods are unacceptable. Instead, the proposed ISA method successfully improves the precision of the estimation. The superiority for the estimation accuracy is effective for identifying small UAVs.

| CONCLUSION
In this article, the proposed ISA method is developed to extract m-D features to identify small UAVs with a dual-channel radar. The multivariate returned signal is applied to the NA-MEMD method to separate the Doppler signal induced by the translation in the first few IMFs, whereas subsequent IMFs exhibit features similar to the rotation signal. Then, the ISA method is deployed on multivariate IMFs exhibiting rotation to separate the rotation signal from other interferences with similar frequency bands. Compared with other existing multivariate reassigned transforms, such as M-SWT, M-SST and CSST, the proposed ISA method successfully enhances the concentration of TF distribution for the monocomponent rotation signal, but it also reduces the TF energy for the Doppler signal and multifrequency false alarms. The effectiveness is validated by synthetic data. Experimental results on small UAVs indicate that the proposed ISA method has brighter prospects for identifying small UAVs.