Identification-While-Scanning of a Multi-Aircraft Formation Based on Sparse Recovery for Narrowband Radar

It is known that the identification performance of a multi-aircraft formation (MAF) of narrowband radar mainly depends on the time on target (TOT). To realize the identification task in one rotated scan with limited TOT, the paper proposes a novel identification-while-scanning (IWS) method based on sparse recovery to maintain high rotating speed and super-resolution for MAF identification, simultaneously. First, a multiple chirp signal model is established for MAF in a single scan, where different aircraft may have different Doppler centers and Doppler rates. Second, based on the sparsity of MAF in the Doppler parameter space, a novel hierarchical basis pursuit (HBP) method is proposed to obtain satisfactory sparse recovery performance as well as high computational efficiency. Furthermore, the parameter estimation performance of the proposed IWS identification method is analyzed with respect to recovery condition, signal-to-noise ratio and TOT. It is shown that an MAF can be effectively identified via HBP with a TOT of only about one hundred microseconds for IWS applications. Finally, some numerical experiment results are provided to demonstrate the effectiveness of the proposed method based on both simulated and real measured data.


Introduction
Multi-aircraft formation (MAF) identification is an important task [1][2][3][4] for air search radar, where both the aircraft number and their motion information are of particular interest. For the narrowband radar, aircraft in the MAF can be densely distributed in a single beam and a single range cell of echoes. Although wideband waveforms have been introduced to distinguish between aircraft in difference range cells [4][5][6][7][8][9][10], a high A/D sampling frequency and huge data storage are required. Moreover, wideband waveforms may bring difficulties on target detection due to the energy being distributed along different range cells and the unknown range migration among different high range resolution profiles (HRRPs) [1]. Therefore, narrowband radar systems are still widely used in air search radars [1,11,12]. Though different aircraft in MAF may not be separated along the low range resolution profiles, different Doppler parameters, which are caused by different positions and motions of targets during a long time on target (TOT), provide a new approach for MAF identification [1,13,14]. In this regards, the third-order polynomial Fourier transform (PFT) [1] and many other time-frequency analysis (TFA) methods [13,14] have been proposed in a long TOT, e.g., one or several seconds. However, a long TOT cannot be obtained in many scenarios due to the following three reasons. First, the TOT is an important resource for radar, especially for wide-area search radar, and it needs to allocate the illumination time to detect targets from different directions. Second, modern radar is expected to have a quick response to environmental change, and the MAF information should be determined as soon as possible. Third, a dynamic frame-by-frame observation of MAF is expected in many applications to reflect the MAF time-varying motion. Specifically, for mechanically rotated or electronically scanned radar, it is hoped that an MAF can be identified in a single scan without decrease of the rotation speed as well as the data rate for target report. Therefore, the identification-while-scanning (IWS) in a short TOT is preferred for MAF identification in many practical applications.
The echoes of the MAF should be composed by different polynomial phase signals (PPS) with respect to different aircraft [1]. Furthermore, in this paper the echoes are derived as the composition of a multiple chirp signals (MCS) model in a short TOT of a single scan. Thus the IWS issue in a short TOT is equivalent to an estimation of multiple chirp components, in which the main problem is how to realize super-resolution of these components in a short TOT for IWS application. Fortunately, as the number of aircraft is limited, the MCS can be regarded as sparse in the whole parametric space. Therefore, in this paper the MAF identification is investigated by introducing sparse recovery (SR), which is recently discussed in radar area, especially for synthetic aperture radar (SAR) imaging [11,12,[15][16][17][18][19][20].
The SR theory indicates that if a signal is sparse on some basis, it can be exactly reconstructed from limited samples with high probability. SR is originally realized by l 0 norm minimization, which is an NP-hard problem and intractable to solve directly. Fortunately, if the measurement matrix satisfies the restricted isometry property (RIP) [21] with appropriate restricted isometric constants, the l 0 norm minimization can be approximated to the solvable l 1 norm minimization, which can be realized by convex optimization or greedy algorithms. Compared with greedy algorithms, such as orthogonal matching pursuit (OMP) [22,23], convex optimization algorithms like basis pursuit (BP) [24] have higher recovery accuracy but much more computational complexity. In the MAF identification issue, convex optimization is preferred since the target number needs to be estimated accurately. On the other hand, off-grid problem of SR [25,26] may be inevitable since the real Doppler parameters of MAF can be any continuous real values while the grid of the basis are pre-discretized. Although refining the grids can increase the recovery probability, it will also create a heavy computation burden and recovery efficiency loss. Therefore, to solve the off-grid problem and maintain a moderate computation complexity, hierarchical recovery has been proposed in the field of SR. In literatures [27,28], a hierarchical matching pursuit is proposed and applied to sparse coding in image classification. In [29,30], grid refinement is utilized for better SAR and ISAR imaging. In these former works, though, the effective grid setting rule is not well discussed. In this paper, a hierarchical basis pursuit (HBP) method is proposed. It is clearly shown that the MAF can be effectively identified via HBP in a short TOT of only one hundred microseconds for IWS applications.
The remainder of the paper is organized as follows. In Section 2, the MAF signal model is established based on narrowband signals in a short TOT. In Section 3, the MAF identification is discussed by SR and the HBP method for IWS application is proposed with detailed performance analysis. In Section 4, some results of simulations and real measured data are provided to demonstrate the effectiveness of the proposed method. In Section 5, some conclusions are drawn.

Multi-Aircraft Formation Signal Model in Short Coherent Integrated Time
We have proposed a two-dimensional MAF signal model in [1] in long coherent integrated time (LCIT), e.g., several seconds. In this paper, we focus on the model for IWS applications in a short TOT, e.g., from several tens to one hundred micro-seconds, in a single scan for coherent integration. The geometry of an MAF scenario is rewritten in Figure 1 with two aircraft for simplicity. Three Cartesian coordinates are introduced for the convenience of derivation. UOV is the radar coordinate and the radar is located at the coordinate origin, O. XOY is the reference coordinate, where the axis Y coincides with the axis V and X is parallel to U. |Oo| = R 0 and |·| is the absolute value operator. The third coordinate {x (k) o (k) y (k) , k = 1, ..., K} is the target coordinate, where o (k) is the k-th aircraft geometry center, K is the aircraft number of the MAF and the origin of the k-th targets' in XOY. The aircraft are assumed in the same beam and the same range cell from the constraint where θ is the azimuth beam width, ∆R is the range resolution of the narrowband radar. The radar light of sight (LOS) angle is θ (k) 0 , corresponding to the angle ∠UOo (k) . The radial direction is denoted as the axis of Y and the tangential direction is the axis of X.
Sensors 2016, 16,1972 3 of 14 geometry center, K is the aircraft number of the MAF and the origin of the k-th targets' coordination is located at The aircraft are assumed in the same beam and the same range cell from the constraint where  is the azimuth beam width, R  is the range resolution of the narrowband radar. The radar light of sight (LOS) angle is . The radial direction is denoted as the axis of Y and the tangential direction is the axis of X. Normally, the motions of aircraft in the MAF cannot change rapidly for the safety reason, thus only the first and second order motions needs to be considered, i.e., the radical velocity in radar coordinates and   T  is the transpose operator. The range from the scatterer P to the radar versus t can be expressed as Usually, closely spaced scatterers cannot be effectively identified in a limited illumination time, even in several seconds of LCIT, so the phase differences among scatterers in a same aircraft can be omitted, which means that each aircraft in MAF contributes one Doppler component. Therefore, the MAF signal model can be approximated to the combination of K components of four-order polynomial phase signal (PPS) as where   k l f is the first order phase coefficient, 1, 2,3, 4, l  which can be expressed as Normally, the motions of aircraft in the MAF cannot change rapidly for the safety reason, thus only the first and second order motions needs to be considered, i.e., the radical velocity v   x , should be considered.
Thus, an arbitrary scatterer P of the k-th target is located at x in radar coordinates and (·) T is the transpose operator. The range from the scatterer P to the radar versus t can be expressed as Usually, closely spaced scatterers cannot be effectively identified in a limited illumination time, even in several seconds of LCIT, so the phase differences among scatterers in a same aircraft can be omitted, which means that each aircraft in MAF contributes one Doppler component. Therefore, the MAF signal model can be approximated to the combination of K components of four-order polynomial phase signal (PPS) as where f (k) l is the first order phase coefficient, l = 1, 2, 3, 4, which can be expressed as where σ (k) is the backscattering amplitude of the k-th target, Ω (k) (t) is the equivalent rotation velocity generated by the k-th target's tangential motion. Usually, the rotation velocity Ω (k) (t) of each target can be assumed to be identical and denoted by Ω.
The linear phase coefficient f 1 causes the Doppler shift, while the quadratic term, cubic term and fourth-order term cause the Doppler dispersion. In the TOT, the Doppler frequency dispersion caused by the quadratic term, cubic term and fourth-order term are respectively denoted by Fds 2 , Fds 3 and Fds 4 as Fds (k) Normally, the velocity v (k) x,0 will be several hundred meters per second, the acceleration a (k) x < 10 m/s 2 , and the radial range between the MAF and the radar R 0 > 10 km. For the application of IWS, the TOT T is about 100 ms. Also, the wavelength of far-range search radar will be larger than 3 cm, the Doppler frequency dispersion by the cubic and fourth-order term meets where ρ f = 1/T represents the Doppler center resolution cell. Equations (11) and (12) indicate that the dispersion Fds 4 and Fds (k) 3 can be omitted in the short during time T as well as the cubic and fourth-order phase modulation. Moreover, the quadratic term Fds   3 . Thus, the MAF signal versus t for IWS application can be approximated as where the k-th aircraft is decided by its Doppler center f ak and Doppler rate f bk , jointly. Therefore, the narrowband coherent MAF signal in IWS application is equivalent to a model of K components of chirp signals, i.e., multiple chirp signal (MCS) model. The 2D geometry UOV established in is an observation plane composed of the radar LOS and the instant velocity vector of MAF at t = 0. Actually, the targets of MAF may move out of the UOV plane in a real 3D space. Fortunately, in a short TOT, e.g., one hundred microseconds for IWS usage, the 3D motion can well be assumed bounded in the plane of UOV and the high-order motions can be omitted. Thus the MAF signal still can be modeled as the combination of different chirp components with respect to different aircraft in a 3D space. That is, the dimension of the matrices and the complexity of the recovery in a 3D geometry will be as same as those in 2D geometry in a short TOT.
Besides, for safety reasons, the aircraft in a stable MAF may normally keep the identical motion, which is called rigid MAF formation. In that case, the MCS can be verified to have the same Chirp rate [31]. However, in some scenarios, the aircraft in a MAF may have different motions and change their relative positions in the flight, which is called non-rigid formation. Also, two different MAF formations may possibly appear in the same range cell of the same beam. Therefore, this paper mainly discusses the non-rigid MAF identification as in Equation (13), and the rigid one can be regarded as a special case of Equation (13).

IWS MAF Identification via Sparse Recovery Methods
It is shown in Equation (13) that MAF identification can be realized for IWS application by estimating the number of chirp components as well as their parameters. However, the TOT in one IWS single scan is too short to discriminate each chirp component for conventional methods. Therefore, sparse recovery approaches are introduced in this Section to solve the resolution problem. Moreover, it is shown that the number of chirp components as well as their parameters, i.e., Doppler centers and Doppler rates, can be estimated by the sparse recovery, simultaneously. Based on Equation (13), the dictionary of sparse recovery consisted of the chirp atoms with different Doppler centers and rates for the over complete representation. However, retrieving the chirp information is a two-dimensional parameter search problem with a considerable computational burden. To reduce the computational complexity and sustain a high accuracy for the convex optimization, a HBP recovery method is further proposed for MAF identification in this section.

Overview of Sparse Recovery Theory
The framework of SR indicates that with the assumption of sparsity, the N × 1 dimensional signal vector s can be expressed as where Ψ is a N × Z dimensional basis matrix, f is a Z × 1 dimensional reconstructed signal vector with sparsity of K, i.e., only K elements of f are nonzero and K << Z. From the compressed sensing theory, the s can be recovered from its limited observation y with length of M and M << N, which can be expressed as where y is the measurement vector and Φ is an M × N measurement matrix. Then we can obtain where the dictionary Θ = ΦΨ. According to Equations (14)- (16), the recovery of s is equivalent to the recovery of f. Thus the following will only focus on the recovery of f. If the measurement matrix Θ is chosen properly, f will be recovered from the locations of the nonzero entries of f.
The principle of BP is to find a representation of the signal whose coefficients have minimal l 1 norm, as represented by Equation (17) min||f|| 1 where y is the samples of the signal in a noise-free environment.
In the case of standard white Gaussian noise with deviation of σ, Basis Pursuit denoising (BPDN) is proposed, which is the solution of Assuming the dictionary is normalized, then the parameter λ is related to the noise power, usually set as λ C = σ 2log(C), with C is the cardinality of the dictionary. BPDN is implemented by searching the best solution to Equation (18) on the whole f axis, which requires a large amount of time for computation. To make BPDN feasible in practical application, a HBP method is proposed in this paper by hierarchically search to reduce the computation burden remarkably.

MAF Identification by HBP Method
At first we prove the sparsity of the formation signal represented in the transform domain. Rewrite the formation signal in Equation (13) in the discrete form as and s = (s [1] , s [2] , · · · , s [N]) T . Form the basis matrix Ψ in a row vector consisted by different Doppler centers and Doppler rates as where the subscript i in p i varies from 1 to P, and for each i the subscript j in q j traverses from 1 to Q. Thus there are P × Q combinations for p i and q j , corresponding to each columns in Ψ.
Assume that only the first M samples of the formation signal s[n] can be obtained in one scan, i.e., y = (s [1] , · · · , s [M]) T , which means the observation matrix Φ has the form as Then the dictionary Θ can be written as Thus, the MAF samples y is represented by The number of the nonzero elements of f in the basis matrix Ψ corresponds to the aircraft number, and the aircraft Doppler parameters can be estimated by the chirp basis corresponding to the nonzero elements in f. Specifically, the locations of the nonzero elements in f represent the aircraft's Doppler center and Doppler rate, jointly. Therefore, f is sparse with only a few non-zero entries, which can be solved by BPDN in Equation (18).
Since the computation for the optimization of Equation (18) largely depends on the size of the measurement matrix Θ, the atoms of Θ can be designed in two hierarchies with different grids.
Specifically, in the first hierarchy, the atoms in the measurement matrix Θ 1 , expressed as θ 1 (p i , q j ), should have a wider range to ensure that the signal components can totally be covered. To keep Θ 1 as small as possible, the grid among the atoms is allowed to be large to a great extent. While in the second hierarchy, atoms in the measurement matrix Θ 2 have a much smaller grid, and only need to cover the small range of the non-zero atoms based on the sparse recovery result of the first hierarchy. Compared to the traditional BP method, the HBP method can keep the same accuracy and reduce the computational complexity significantly.
An example is given to demonstrate the implementation of the HBP strategy. The signal is supposed to be constructed by two chirp components with Doppler center f a1 = 1.4 Hz, f a2 = 7.2 Hz and Doppler rate f b1 = 2 Hz/s, f b2 = 9 Hz/s respectively. The CPI is 100 ms. In the first hierarchy, the grids for the Doppler center and rate are set as 2 Hz and 10 Hz/s, as shown in Figure 2. In the second hierarchy, the dictionary is finer with more intensive grids of 0.1 Hz and 1 Hz/s. As shown in Figure 2 in the first hierarchy, the signal can be recovered by the items of the dictionary around the true value of the rough parameters. In the second hierarchy, the signal is recovered by the exact items as the Doppler parameters fall on the grid. The example above indicates that the tuning of the grid is important for the successfully recovery of the signal, which will be discussed in our future work.
respectively. The CPI is 100 ms. In the first hierarchy, the grids for the Doppler center and rate are set as 2 Hz and 10 Hz/s, as shown in Figure 2. In the second hierarchy, the dictionary is finer with more intensive grids of 0.1 Hz and 1 Hz/s. As shown in Figure 2 in the first hierarchy, the signal can be recovered by the items of the dictionary around the true value of the rough parameters. In the second hierarchy, the signal is recovered by the exact items as the Doppler parameters fall on the grid. The example above indicates that the tuning of the grid is important for the successfully recovery of the signal, which will be discussed in our future work.

The Discussion on Performance Analysis and the Grid Setting Rule for the HBP Method
To characterize the recovery ability of a given dictionary, the RIP and the coherence measure are two main metrics. Normally, the former is an NP-hard problem [19] to determine the RIP constants while the latter may be a simpler choice [32][33][34]. Although coherence measure is less accurate than RIP as the recovery condition, it has more explicit expression and is physically meaningful in our case of Equation (23). Therefore, in this section, the coherence is discussed and analyzed to test the recovery ability of the dictionary Θ. The coherence [17] of Θ is defined as the largest absolute inner product between any two columns θi, θj in Θ   where ,   is the inner product operator. According to Theorem 12 of chapter 1 in [35], a sparse vector f can be successfully recovered by CS framework given a length-N observed echo y and the dictionary Θ, as long as the aircraft number K satisfies According to the Welch bound [36], the coherence of a basis matrix is always in the range

The Discussion on Performance Analysis and the Grid Setting Rule for the HBP Method
To characterize the recovery ability of a given dictionary, the RIP and the coherence measure are two main metrics. Normally, the former is an NP-hard problem [19] to determine the RIP constants while the latter may be a simpler choice [32][33][34]. Although coherence measure is less accurate than RIP as the recovery condition, it has more explicit expression and is physically meaningful in our case of Equation (23). Therefore, in this section, the coherence is discussed and analyzed to test the recovery ability of the dictionary Θ. The coherence [17] of Θ is defined as the largest absolute inner product between any two columns θ i , θ j in Θ where ·, · is the inner product operator. According to Theorem 12 of chapter 1 in [35], a sparse vector f can be successfully recovered by CS framework given a length-N observed echo y and the dictionary Θ, as long as the aircraft number K satisfies According to the Welch bound [36], the coherence of a basis matrix is always in the range N − 1)), 1 . If N >> M, the lower bound for µ (Θ) approximates to 1/ √ M. According to Equation (25), the upper bound for the aircraft number K which can be uniquely reconstructed will be (1 + √ M)/2. The coherence of the atoms θ i , θ j in the measurement matrix Θ can be derived as where ∆p = p i i − p j i and ∆q = q i i − q j i . Notably, the coherence expression of Equation (26) has a similar expression to the ambiguity function [37], because they all reflect the correlation with the shifted parameters in a certain dimension.
The proposed sparse recovery method provides a super-resolution way to identify the aircraft. Since the theoretical resolution of the sparse recovery method depends on the dictionary's atoms, the resolution can reach as high as possible if the selected steps of ∆p and ∆q are small enough, ideally. However, the possible resolution is limited by the coherence of the measurement matrix Θ. That is, if the steps of ∆p and ∆q are too small, the MAF signal cannot be reconstructed successfully since the coherences for the atoms in the dictionary are too large. Therefore, it can be inferred that the resolution for sparse recovery methods can be defined as the minimum steps for ∆p and ∆q meeting Equation (25) where K = 2. That is It should be mentioned that the recovery capability defined by Equation (25) is an issue of probability, which gives a strict condition to recover the signal with high recovery probability. However, in most cases which cannot satisfy Equation (25), the sparse recovery can still succeed if the SNR is high enough. In brief, much higher resolution can be obtained for the proposed HBP method than those of Equation (27) to satisfy MAF identification in one IWS scan.
Accordingly, the choice of the hierarchical grids is vital to the performance of sparse recovery. From analysis above, we give the grid setting rule for the MAF as follows. The atoms in the first layer can be rough in a wider range, with the grids in the same magnitude of the resolutions. Then the second hierarchy will reconstruct the signal with higher accuracy. The atoms in the second layer are selected in the neighbors of the recovery result in the first hierarchy with grids as small as one-tenth of the resolutions. It should be pointed that although smaller grids may bring errors for recovery, signals can still be identified well in relatively high SNR conditions. Moreover, in the second layer, the number of atoms is greatly reduced which can improve the recovery performance to some extent.
At last, the computational complexity is discussed. As BP method consists of a solution of the convex optimization problems and can be recast as a linear programming [24,38] if the non-zero components are real, which has the polynomial computation, i.e., typically O(n 3 ), where n is the number of atoms in the dictionary for sparse recovery. Thus it may be computational expensive or infeasible in real applications when n is large. In contrast, the greedy algorithms have been proposed to reduce the computation complexity to O(nK 2 ) with lower recovery accuracy, where K is the number of non-zero entries [38] and K << n. According to the proposed flowchart of HBP in Section 3, the HBP method reduces the computational complexity of BP [24] remarkably from O(n 3 ) to O(n 3 1 ) + O(n 3 2 ), where O(n 3 1 ) is the computational complexity of the first hierarchy search and n 1 = n/C 1 is the atom number in the dictionary of the first hierarchy and C 1 is the grid distance times of first hierarchy grid on the original dense grid. O n 3 2 is the second hierarchy search complexity and n 2 = C 2 K is the atom number in dictionary of the second hierarch, where K is the sparsity, i.e., the aircraft number, and C 2 is the search atom numbers neighboring each of the K nonzero entries.

Numeric Experiments of the Proposed HBP Method
In this section, the results of both the simulation data and the real measured data are provided to demonstrate the effectiveness of the proposed HBP method for IWS applications. The system parameters are listed as follows. The signal carrier frequency f c = 1 GHz, the pulse repetition frequency (PRF) is 300 Hz, transmitting signal bandwidth B = 1 MHz and the sampling frequency f s = 2 MHz. For a typical narrow-band coherent radar, assume that the antenna beam width is 4 • with rotation speed 6 revolutions per minute (RPM), then the short TOT for coherent integration will be about 100 ms for the MAF identification in IWS applications.

MAF Performance Analysis
The aircraft number of the MAF is two and the slant range between the MAF and radar is 12 km.
The locations in the reference coordinate are X  Table 1, from which the Doppler centers and Doppler rates for the two aircraft can be calculated from the parameters in Table 1.

Motion Parameters Aircraft1 Aircraft2
Radical Tangential Radical Tangential  Figure 3 shows the identification correct probability of the aircraft number versus SNR with 100 Monte Carlo samples, which illustrates that the correct probability reaches more than 90% when the SNR is higher than 30 dB. Figure 4a,b show the recovery accuracy for the Doppler center and the Doppler rate with 100 times of Monte Carlo simulation. The estimation accuracy of the proposed HBP method is compared with the CRLBs for each parameter. In Figure 4, both accuracies of the Doppler centers and Doppler rates all approaches to the CRLBs when SNR increases gradually, which shows that the proposed HBP method can obtain high accuracies.

Numeric Experiments of the Proposed HBP Method
In this section, the results of both the simulation data and the real measured data are provided to demonstrate the effectiveness of the proposed HBP method for IWS applications. The system parameters are listed as follows. The signal carrier frequency

MAF Performance Analysis
The aircraft number of the MAF is two and the slant range between the MAF and radar is 12 km. The locations in the reference coordinate are    Table 1, from which the Doppler centers and Doppler rates for the two aircraft can be calculated from the parameters in Table 1.

Motion Parameters Aircraft1 Aircraft2
Radical Tangential Radical Tangential  Figure 3 shows the identification correct probability of the aircraft number versus SNR with 100 Monte Carlo samples, which illustrates that the correct probability reaches more than 90% when the SNR is higher than 30 dB. Figure 4a,b show the recovery accuracy for the Doppler center and the Doppler rate with 100 times of Monte Carlo simulation. The estimation accuracy of the proposed HBP method is compared with the CRLBs for each parameter. In Figure 4, both accuracies of the Doppler centers and Doppler rates all approaches to the CRLBs when SNR increases gradually, which shows that the proposed HBP method can obtain high accuracies.     With the increase of the TOT, the MAF identification will be more accurate. Figure 5 gives the correct identification probability of the aircraft number with different TOTs. When the TOT added up to 200 ms, the SNR requirement for 92% identification rate was reduced by nearly 15 dB. The identification results for different aircraft numbers are shown in Figure 6a,b. For the TOT of 100 ms, the method can effectively identify as many as six targets with correct identification probability more than 80% when SNR is higher than 60 dB. For the TOT of 200 ms, MAF number can be correctly estimated with a probability of over 95% when SNR is higher than 30 dB.   With the increase of the TOT, the MAF identification will be more accurate. Figure 5 gives the correct identification probability of the aircraft number with different TOTs. When the TOT added up to 200 ms, the SNR requirement for 92% identification rate was reduced by nearly 15 dB.  With the increase of the TOT, the MAF identification will be more accurate. Figure 5 gives the correct identification probability of the aircraft number with different TOTs. When the TOT added up to 200 ms, the SNR requirement for 92% identification rate was reduced by nearly 15 dB. The identification results for different aircraft numbers are shown in Figure 6a,b. For the TOT of 100 ms, the method can effectively identify as many as six targets with correct identification probability more than 80% when SNR is higher than 60 dB. For the TOT of 200 ms, MAF number can be correctly estimated with a probability of over 95% when SNR is higher than 30 dB.   The identification results for different aircraft numbers are shown in Figure 6a,b. For the TOT of 100 ms, the method can effectively identify as many as six targets with correct identification probability more than 80% when SNR is higher than 60 dB. For the TOT of 200 ms, MAF number can be correctly estimated with a probability of over 95% when SNR is higher than 30 dB.  With the increase of the TOT, the MAF identification will be more accurate. Figure 5 gives the correct identification probability of the aircraft number with different TOTs. When the TOT added up to 200 ms, the SNR requirement for 92% identification rate was reduced by nearly 15 dB. The identification results for different aircraft numbers are shown in Figure 6a,b. For the TOT of 100 ms, the method can effectively identify as many as six targets with correct identification probability more than 80% when SNR is higher than 60 dB. For the TOT of 200 ms, MAF number can be correctly estimated with a probability of over 95% when SNR is higher than 30 dB.

Comparison with the PFT Method
In this part, HBP method is compared with the second-order PFT method proposed in [1], with the same MAF parameter of Table 1 with T = 0.1 s. The searching result for the second-order PFT method in f 2 − f 1 plane is give in Figure 7, where the two aircraft cannot be identified due to the short TOT as well as the low parameter resolution. The result shows that the PFT method failed to estimate the MAF information. In fact, only when the coherent time T > 0.8 s can the parameters be well estimated by PFT. The results indicate that the proposed HBP is a super-resolution method while PFT is not suitable for IWS applications.

Comparison with the PFT Method
In this part, HBP method is compared with the second-order PFT method proposed in [1], with the same MAF parameter of Table 1 with T = 0.1 s. The searching result for the second-order PFT method in 2 1 f f  plane is give in Figure 7, where the two aircraft cannot be identified due to the short TOT as well as the low parameter resolution. The result shows that the PFT method failed to estimate the MAF information. In fact, only when the coherent time T > 0.8 s can the parameters be well estimated by PFT. The results indicate that the proposed HBP is a super-resolution method while PFT is not suitable for IWS applications.

MAF Identification Based on Real Measured Data
The real measured data are obtained by an experimental radar system with parameters listed in Table 2. In the experiments, the location and velocity information can be obtained by the record of track for each aircraft, thus the real Doppler information can be calculated to evaluate the performance of the HBP method.   Table 3.

MAF Identification Based on Real Measured Data
The real measured data are obtained by an experimental radar system with parameters listed in Table 2. In the experiments, the location and velocity information can be obtained by the record of track for each aircraft, thus the real Doppler information can be calculated to evaluate the performance of the HBP method.  Figure 8a. The identification result by HBP in one CPI is shown in the Doppler center-Doppler rate plane of Figure 8b, with the parameter estimation result in Table 3.

Comparison with the PFT Method
In this part, HBP method is compared with the second-order PFT method proposed in [1], with the same MAF parameter of Table 1 with T = 0.1 s. The searching result for the second-order PFT method in 2 1 f f  plane is give in Figure 7, where the two aircraft cannot be identified due to the short TOT as well as the low parameter resolution. The result shows that the PFT method failed to estimate the MAF information. In fact, only when the coherent time T > 0.8 s can the parameters be well estimated by PFT. The results indicate that the proposed HBP is a super-resolution method while PFT is not suitable for IWS applications.

MAF Identification Based on Real Measured Data
The real measured data are obtained by an experimental radar system with parameters listed in Table 2. In the experiments, the location and velocity information can be obtained by the record of track for each aircraft, thus the real Doppler information can be calculated to evaluate the performance of the HBP method.  Figure 8a. The identification result by HBP in one CPI is shown in the Doppler center-Doppler rate plane of Figure 8b, with the parameter estimation result in Table 3.  Case 2. The data of MAF is formed by four aircraft. The frequency spectrum and recovery result are given as Figure 9a,b, where both aircraft number and their Doppler parameters are all well identified as seen in Table 4. The data of MAF is formed by four aircraft. The frequency spectrum and recovery result are given as Figure 9a,b, where both aircraft number and their Doppler parameters are all well identified as seen in Table 4.  From the above experiment results, the effectiveness of the HBP for MAF identification is verified and the estimation accuracy for the Doppler center and Doppler rate is about 1 Hz and 2 Hz/s, respectively, which agrees with the simulation results in Section 4.1.

Conclusions
In this paper, MAF identification is discussed for IWS application for narrowband radar, based on a proposed MCS signal model where each aircraft can be represented by a chirp signal with a different Doppler center and Doppler rate. To realize super-resolution discrimination in a short TOT, sparse recovery is introduced for IWS applications in this paper. Furthermore, block-sparse MAF Doppler distribution is exploited to form an HBP method, which can reduce the computation complexity of the convex optimization algorithms and sustain high reconstruction accuracy, simultaneously. Thus the aircraft number as well as their Doppler parameters can be identified via the proposed HBP for each aircraft in MAF. Furthermore, the recovery condition, accuracy, resolution and computational complexity are discussed for the proposed method to show its high performance, respectively. Finally, some numerical experiment results of both simulated data and real measured data are provided to demonstrate the effectiveness of the proposed HBP method.  From the above experiment results, the effectiveness of the HBP for MAF identification is verified and the estimation accuracy for the Doppler center and Doppler rate is about 1 Hz and 2 Hz/s, respectively, which agrees with the simulation results in Section 4.1.

Conclusions
In this paper, MAF identification is discussed for IWS application for narrowband radar, based on a proposed MCS signal model where each aircraft can be represented by a chirp signal with a different Doppler center and Doppler rate. To realize super-resolution discrimination in a short TOT, sparse recovery is introduced for IWS applications in this paper. Furthermore, block-sparse MAF Doppler distribution is exploited to form an HBP method, which can reduce the computation complexity of the convex optimization algorithms and sustain high reconstruction accuracy, simultaneously. Thus the aircraft number as well as their Doppler parameters can be identified via the proposed HBP for each aircraft in MAF. Furthermore, the recovery condition, accuracy, resolution and computational complexity are discussed for the proposed method to show its high performance, respectively. Finally, some numerical experiment results of both simulated data and real measured data are provided to demonstrate the effectiveness of the proposed HBP method.