Abstract
High-resolution direction of arrival (DOA) estimation is a critical issue for mainbeam multitarget tracking in ground-based or airborne early warning radar system. A beam-Doppler unitary ESPRIT (BD-UESPRIT) algorithm is proposed to deal with this problem. Firstly, multiple snapshots without spatial aperture loss are obtained by using the technique of time-smoothing. Then the conjugate centrosymmetric discrete Fourier transform (DFT) matrix is used to transform the extracted data into beam-Doppler domain. Finally, the rotational invariance property of the space-time beam is exploited to estimate DOA of the target. The DOA estimation accuracy is improved greatly because the proposed algorithm takes full advantage of temporal information of the signal. Furthermore, the computational complexity of the presented algorithm is reduced dramatically, because the degree of freedom after beam transformation is very small and most of the operations are implemented in real-number domain. Numerical examples are given to verify the effectiveness of the proposed algorithm.
1. Introduction
Mainbeam multitarget separation is a key problem for multitarget tracking in ground-based or airborne early warning radar system. When multiple targets fly in a formation, it is common that some targets fall into the same range cell or adjacent range cells. In this case, it is difficult to separate these targets in range domain only using pulse compression. Meanwhile, it is also difficult to distinguish them in Doppler domain, since the radial velocities of these targets are usually very close; meanwhile, the coherent accumulation time is usually very short in an early warning radar system. In this paper, we consider resolving closely spaced targets in spatial domain through high-accuracy DOA measurement.
Traditional angle measurement algorithms such as the classical monopulse technique [1] becomes invalid in the presence of multiple targets. In order to resolve multiple targets, a typical way is to improve radar resolution by narrowing antenna beam width, widening bandwidth, or extending dwell time. However, these solutions require additional hardware cost. Another possible solution is to achieve superresolution by employing array signal processing techniques. Classical spatial superresolution methods include maximum likelihood (ML) [2], multiple signal classification (MUSIC) [3], estimation of signal parameters via rotational invariance technique (ESPRIT) [4], and their variants [5–7]. In these methods, only spatial information is used, and the temporal pulses are treated as snapshots. Their performance degrades significantly in the presence of closely spaced targets, especially in the condition of low signal-to-noise ratio (SNR) and small number of snapshots. In fact, for early warning radar, the temporal information such as the difference of Doppler frequencies between two targets is extremely useful for DOA estimation. Essentially, the problem of spatial parameter estimation could be transformed into the problem of space-time parameter estimation if we make use of the temporal information of the signal. In general, the spatial distance between two targets will be amplified in space-time plane. Therefore, the accuracy of spatial parameter estimation, especially in the case of closely spaced targets, could be improved if the temporal information is used.
For the space-time signal model, the problem is to estimate multitarget DOA with a single snapshot. The ML method could be used to joint estimate of DOA and Doppler frequency [8]. However, the ML method requires a 2D parameter search process, which usually results in heavy computational load. In order to improve efficiency, the reduced-dimension ML method was proposed in [9]. In this method, the problem of DOA and Doppler frequency estimation is decoupled into two sequential 1D optimization problems. Nevertheless, this method needs a 1D parameter search. Moreover, it is easy to converge to the local optimal solution. The subspace-based superresolution methods combined with aperture smoothing could be practically extended for space-time parameter estimation, such as the 2D MUSIC method [10]. However, the 2D MUSIC needs a 2D parameter search yet. The joint angle and frequency estimation (JAFE) method based on ESPRIT is proposed to estimate DOA and carrier frequency of multicarrier signal for the communication system [11]. The JAFE method does not need parameter search, and thus the computational complexity of JAFE is much lower than that of search-based subspace method. The real value-based unitary JAFE (U-JAFE) is proposed in [12] with lower computational load than that of JAFE. However, in the JAFE and U-JAFE methods, the 2D parameters are estimated in element-pulse domain, which means the signal subspace is also estimated via eigenvalue decomposition (EVD) or singular-value decomposition (SVD) in the element-pulse domain. In the case of large temporal degree of freedom (DOF) or spatial DOF, the computational complexity of JAFE and U-JAFE is still very high. Recently, compressive sensing [13–15] (CS) is widely studied and has been applied in the array signal processing community [16–18]. Although CS could be extended to estimate the space-time parameter and a reasonable performance could be obtained in the case of small snapshots [19], it needs to solve a high-dimensional optimization problem with high computational complexity.
In this paper, the problem of estimating DOAs of multiple targets within mainbeam is considered. In order to improve the DOA estimation accuracy, the temporal information of these targets is exploited. Considering the computational complexities of existing methods are very high in original element-pulse domain, a beam-Doppler unitary ESPRIT algorithm is proposed, which aims to reduce the computational load with improved DOA estimation accuracy. Firstly, multiple space-time snapshots without array aperture loss are extracted from a single sample via time-smoothing. Then, a few space-time beams are formed pointing towards the targets according to the coarse-resolution information provided by the target detection stage. Evidently, the parameter could be estimated in the low-dimensional space-time beamspace with a reduced computation load. Further, the conjugate centrosymmetric DFT matrix is utilized for beamforming with a rotational invariance property being retained. Finally, the real value-based ESPRIT is carried out for estimating DOA with further reduced computational complexity. The proposed method is quite different from the beamspace ESPRIT method [20], in which only spatial information is utilized for DOA estimation. It is also different from the JAFE [11] and U-JAFE [12] methods, which estimate DOA in element-pulse domain, while our method is performed in a low-dimensional beam-Doppler domain.
The remainder of this paper is organized as follows. Section 2 develops the data model for an early warning radar system. Section 3 proposes a beam-Doppler unitary ESPRIT method for joint DOA and Doppler frequency estimation with some analyses of computational complexity and Cramer-Rao bound (CRB), and Section 4 gives some discussion of critical aspects for the proposed method. In Section 5, simulation results are presented to illustrate the performance improvement introduced by the proposed algorithm, while the conclusions are given in Section 6.
Notation. , , and denote transpose, conjugate-transpose, and inverse operations, respectively; represents the Kronecker product; diag(v) stands for diagonal matrix whose diagonal element is a vector v; is a identity matrix; and denotes the real part and imaginary part of a complex number, respectively.
2. Problem Formulation
In this paper, we consider a narrowband pulse-Doppler (PD) early warning radar system operating with a rectangular planar phased array. The sketch map of target detection for early warning radar is depicted in Figure 1, where the transmitting beam is formed by the whole antenna arrays and the receiver is composed of N uniformly spaced column-synthesized subarrays. The pulse number for coherent accumulation is K, assuming that statistically independent targets are located in the same range cell. The DOA and Doppler frequency of the th () target is and , respectively. After range matched filtering, the space-time data vector of one range cell is where is a steering matrix; denotes the space-time steering vector of the pth target; and are the spatial steering vector and temporal steering vector of the pth target, respectively, where is the normalized spatial frequency, is the normalized Doppler frequency, d is the spacing of column subarrays, is the operating wavelength, and is the pulse repetition period; is a vector that consists of complex amplitudes of targets, where is the complex amplitude of the pth target; and is the complex Gaussian white noise with zero mean and covariance matrix , where denotes the variance of the white noise. It is worth noting that both the spatial steering vector and temporal steering vector are conjugate centrosymmetric.
In order to track multiple targets within one beam (as shown in Figure 1), it is required that the radar system should separate these targets firstly. In this paper, we consider resolving closely spaced targets within one beam through high-accuracy DOA measurement. The DOA information is also useful for subsequent target tracking process. From the space-time signal model in (1), it is shown that the DOA estimation could be casted as the problem of multidimensional parameter estimation with a single snapshot. The ML method [9] or JAFE method [11] could be used to deal with this problem. However, they require a parameter search or high-dimensional eigenvalue (or singular-value) decomposition, which leads to a heavy computational load. In the next section, we propose a beam-Doppler unitary ESPRIT algorithm with reduced computational cost and increased DOA estimation accuracy.
3. Proposed Beam-Doppler Unitary ESPRIT Method
For the early warning radar system, the target parameter estimation is usually carried out after target detection. And the target detection is usually performed in the range-Doppler plane with PD processing [1]. Although the resolution of conventional PD processing is unable to break the Rayleigh limit [1], the coarse-resolution location information of targets is available. Hence, a few space-time beams could be formed to cover the region of interest. Then, the parameter estimation could be implemented in the low-dimensional space-time beamspace with a dramatically reduced computational complexity.
3.1. Time-Smoothing
The ESPRIT method needs multiple snapshots to estimate the signal subspace. However, there is only one snapshot for the space-time signal model presented in this paper. In order to cure this problem, the time-smoothing method is utilized to extract snapshots without spatial aperture loss. In general, for the early warning radar system, the temporal aperture is much easier to increase than is spatial aperture by increasing the pulse number. Assume that the length of temporal subaperture (also called time-smoothing factor) is M (). Then, a total of snapshots could be acquired by forward time-smoothing. The th () snapshot could be expressed as where is the selection matrix for choosing the lth space-time subaperture, and is the selection matrix for choosing the lth temporal subaperture. By substituting (1) into (2), we obtain where denotes the subaperture steering matrix; denotes the space-time subaperture steering vector, where is the temporal subaperture steering vector of the pth target; denotes the target complex amplitude vector of the lth snapshot; and is the noise vector of the lth snapshot. It is shown that the temporal subaperture steering vector is still conjugate centrosymmetric. Subsequently, the smoothed data will be transformed into beam-Doppler domain, and unitary ESPRIT will be implemented for estimation of space-time parameters.
3.2. Unitary ESPRIT in Beam-Doppler Domain
To retain the rotational invariance structure in the space-time beamspace, the conjugate centrosymmetric DFT matrix is applied as the transformation matrix [20]. Assume that the spatial beam transformation matrix and temporal beam transformation matrix are and , respectively. The th () column of is where the nth spatial beam is formed with the spatial frequency . Similarly, the th () column of is where the mth temporal beam (i.e., Doppler bin) is formed with the Doppler frequency . Then, the space-time beam transformation matrix could be expressed as
The lth snapshot in beam-Doppler domain can be obtained as where stands for the noise vector in beam-Doppler domain, and denotes the space-time beamspace steering matrix, and it can be expressed as where and are the beamspace steering matrix and Doppler-domain steering matrix, respectively, while and denote the beamspace steering vector and Doppler-domain steering vector of the pth target, respectively. The nth spatial beam response of can be described as
For all targets, the beamspace steering matrix satisfies a rotational invariant relationship, that is,
The derivation of (10) is given in the appendix, where is a real-valued diagonal matrix whose diagonal elements contain the desired DOA information. Analogously, the rotational invariance relationship for Doppler-domain steering matrix is derived as where and are defined similar to and with N replaced by M, and is a real-valued diagonal matrix whose diagonal elements contain the desired Doppler information. With some manipulation and the Kronecker product property (i.e., ), the rotational invariant relation in space-time beamspace can be derived as where and are beamspace selection matrices, and and are Doppler-domain selection matrices.
Note that a real-valued matrix could be constructed by the beam-Doppler data. The matrix has a property of doubling the sample number. Performing singular-value decomposition on , we have where and are the and left singular vector matrix and right singular vector matrix, respectively, and where , and are the nonzero singular values of . Thus, the real-valued signal subspace can be constructed by the “largest” left singular vectors, that is,
In theory, the signal subspace could be spanned by the column of [20], that is, where is a real-valued nonsingular matrix. Substituting into (12), one gets where and . Therefore, and could be obtained by solving (17) via least square (LS) or total least square (TLS) algorithm, respectively. It should be noted that and are real-valued matrices. Therefore, automatic pairing of the estimated spatial frequency and Doppler frequency can be realized by decomposition of complex matrix , that is,
That is to say, and () could be estimated by the real part and imaginary part of , respectively. Finally, the DOA and Doppler frequency of the pth target could be estimated as follows: where is the th () eigenvalue of .
In summary, the main steps of the proposed method for DOA and Doppler frequency estimation are as follows:
Step 1. Perform the time-smoothing to extract the subaperture data via (2).
Step 2. Construct the beam transformation matrix using the prior location information provided by target detection via (6).
Step 3. Transform the smoothed data into the space-time beamspace via (7) and obtain .
Step 4. Construct the real-valued data matrix and compute the signal subspace via the “largest” singular vectors of .
Step 5. Solve (17) by using TLS to estimate and , where the selection matrices with dimension reduction are constructed by the appropriate subblocks [20] of , , , and .
Step 6. Perform the EVD to the complex matrix and estimate the DOA and Doppler frequency via (19).
3.3. Computational Complexity Analysis
The algorithm proposed in this paper uses complex value operations except for the beam transformation and the final space-time parameter estimation, and the rest of them are real-valued operations. Here we mainly consider the computational complexity of multiplication, since the number of clocks consumed by adding operations could be ignored compared with that of multiplication. In addition, the number of clocks consumed by complex multiplication is about 4 times that of real multiplication [21]. Therefore, the complexity of real multiplication is adopted for computational cost measurement of various methods.
The computational cost of the proposed method is decomposed as follows: Firstly, the computational complexity of beam transformation is , where and denote the number of spatial beams and the number of Doppler bins, respectively. Secondly, the computational complexity of real-valued SVD is . Thirdly, the computational complexity of twice TLS is . Finally, the computational complexity of complex EVD is . A comparison of computational complexities of the beamspace unitary ESPRIT (B-UESPRIT) [20], the U-JAFE method [12], and the proposed method is given in Table 1. It is shown that the computational cost of the proposed method is much less than that of U-JAFE when and but slightly larger than that of B-UESPRIT due to the utilization of temporal information. A much more detailed comparison will be given in Section 5.
3.4. CRB for the Space-Time Data Model
Similar to the derivation in [11], the CRB for DOA estimation under the space-time signal model (1) could be derived as where , , , , and stands for the th column of steering matrix .
4. Performance Discussion
4.1. Influence of Source Number Detection on Target Tracking
Generally, the ESPRIT-like methods assume that the number of sources is known as a prior. In practice, the number of sources should be estimated before DOA estimation. Various methods have been proposed for source number detection, such as the information theoretic criteria-based approaches [22, 23]. It is worth noting that no matter which method is adopted, the estimate of source number may be biased with a certain probability (especially in the case of small samples and low SNR). For the target tracking in early warning radar, the overestimated source number will introduce false tracking points. These false tracking points are generated by noise with randomness and thus would make it difficult to form a stable and reliable track. On the other hand, the underestimated source number will lead incomplete track. In general, as long as the radar system does not successively lose tracking points, the effect on target tracking is weak and could be omitted.
4.2. Whitening as the Preprocessing Stage
Due to the fact that the successive subapertures of time-smoothing overlapped each other, there exists a correlation between the noise components of adjacent snapshots. This means that the noise is colored after time-smoothing. Although the parameters of the proposed method are estimated in space-time beamspace, the beam transformation is a linear transformation. Therefore, the noise in the beamspace is still colored. The estimated signal subspace is biased in the presence of color noise. In this context, the proposed algorithm could be preprocessed with a whitening filter. Considering the noise component in (7), and its noise covariance matrix is given by
is usually unknown in practice. However, it can be estimated by the pure noise data obtained from the clear zone of range-Doppler map. Performing the EVD on , one gives , where is a diagonal matrix that consists of eigenvalues of and is the corresponding eigenvector unitary matrix. Then the whitened beamspace data can be expressed as [24]
Assume that the signal subspace estimated by the whitened data is . It is easy to derive that the relationship between the signal subspace with and without whitening is
Substitute (23) into (17), and the proposed algorithm can still be implemented in the usual way. We will analyze the impact on the DOA estimation with and without whitening in the simulation stage.
4.3. Optimum Time-Smoothing Factor
In this paper, the technique of time-smoothing is utilized for acquiring snapshots. As discussed in the previous section, larger time-smoothing factor (i.e., the length of temporal subaperture) will lead to fewer snapshots, and vice versa. Both of them have significant influence on the performance of parameter estimation. On the one hand, the estimation accuracy of signal subspace and space-time parameter will be improved with the increase of snapshots. On the other hand, the estimation accuracy of space-time parameter, particularly the temporal parameter, will decrease with the shortened temporal subaperture. Therefore, there exists a tradeoff between the time-smoothing factor and parameter estimation performance. The error of DOA estimation is actually a function of time-smoothing factor. The optimum time-smoothing factor for the JAFE method is about one half of the total temporal aperture [11]. However, the proposed method is performed in beam-Doppler domain, whereas the JAFE method is performed in element-pulse domain. Consequently, the conclusion presented for JAFE is not suitable for our method. We have found that the optimum time-smoothing factor for our method is about 0.4K via Monte Carlo simulations.
5. Simulation Examples
In the simulation example, we consider an 8-element uniform linear array with element spacing of . We assume that two far-field signals with equal power and within the mainbeam are impinging on the antenna array. The source signals are narrowband signals. The pulse repetition frequency is 1000 Hz, and the number of coherent pulses is 32. The DOA and Doppler frequency of are and 463 Hz, and those of are and 537 Hz, respectively. Both the angular spacing and the Doppler spacing of these two targets are two-thirds of Rayleigh limit. All simulation results are based on Monte Carlo runs. The evaluation metric, that is, root-mean-square error (RMSE) of DOA, is defined as where denotes estimated DOA of the pth target at the lth Monte Carlo experiment.
5.1. Investigating the Effectiveness of Prewhitening
Results comparing the performances of the proposed algorithm implemented with and without prewhitening (viz. Section 4.2) are shown in Figure 2. (For clarity, only behaviors corresponding to the first source are shown. Similar situations are observed for the second source as well.) In the simulation, the time-smoothing factor is chosen to be M = 23, and the number of spatial beams and Doppler bins are chosen to be and , respectively. From the results, it is shown that the whitening process has a minor effect on the DOA estimation. This phenomenon is consistent with the observations in [11].
5.2. Influence of the Number of Spatial Beams and Doppler Bins on DOA Estimation Error
First, the effect of the number of Doppler bins on DOA estimation accuracy is presented. In this simulation, the time-smoothing factor and the number of spatial beams are chosen to be and , respectively. The signal-to-noise ratio (SNR) is set to be 20 dB. The RMSE of DOA with different numbers of Doppler bins is shown in Figure 3. It is shown that the DOA estimation accuracy improves with the increased Doppler bins, but the estimation accuracy becomes stable when the Doppler bins are more than one. There are two main reasons holding up this phenomenon: The first reason is that temporal DOFs required for enlarging the distance of sources in space-time domain are at least two, and the second reason is that the temporal information is mainly contained in a few Doppler bins covering targets. Considering the tradeoff between computational complexity and estimation accuracy, the number of Doppler bins is chosen to be in our method.
Next, we investigate the effect of the number of spatial beams on DOA estimation accuracy. The RMSE of DOA with different number of spatial beams is shown in Figure 4. In this simulation, the time-smoothing factor and the number of Doppler bins are chosen to be and , respectively. The SNR is set as 20 dB. From Figure 4, it is shown that the DOA estimation error decreases with the increase of spatial beams, but the estimation error decreases more slowly when the spatial beams are more than 2. This phenomenon is attributed to the fact that the spatial DOFs for resolving 2 targets are at least 3 since the ESPRIT-like method consumes one DOF to construct the rotational invariance relationship [4], and the energy (or information) of targets is mainly included in the transformed beams around the target angle, that is, the mainbeam and its neighboring beams (assume that the mainbeam is the th beam, and then the target information is mainly included in the ()th, th, and ()th beams).
5.3. Choice of Optimum Time-Smoothing Factor
In this simulation, the SNR is set to be 30 dB. The number of Doppler bins and spatial beams are chosen to be and , respectively. Figure 5 plots the variation of time-smoothing factor with the parameter estimation errors. It is shown that the DOA estimation error is minimum for . It should be noted that time-smoothing not only can improve the estimation accuracy but also can provide robustness against rank loss when there exist multiple signals with the same DOA.
5.4. Performance Comparison with the ESPRIT-like Algorithms
In this simulation, the ESPRIT-like algorithms such as B-UESPRIT [20] and U-JAFE [11] are chosen for performance comparison. The proposed method is denoted as BD-UESPRIT for short, and the number of Doppler bins, spatial beams, and time-smoothing factor is chosen to be , , and , respectively. The B-UESPRIT method is performed only in space domain, and all the temporal pulses are treated as snapshots. The number of spatial beams for B-UESPRIT is chosen as . The time-smoothing factor for the U-JAFE method is set to be the same as that for the proposed BD-UESPRIT method.
The variation of DOA estimation errors with SNR is shown in Figure 6. The CRB for DOA estimation is also plotted. It is shown that the proposed BD-UESPRIT algorithm has a much higher DOA estimation accuracy than has B-UESPRIT. This is because the proposed method takes full advantage of the temporal information of the target. Another observation is that the proposed method has a lower estimation error than has U-JAFE in the condition of small SNR. This is intuitive since the target energy is concentrated in a few beams after beam transformation.
5.5. Computational Load Comparison with the ESPRIT-like Algorithms
In this subsection, computational complexities of the method mentioned in Table 1 are analyzed. The number of targets, Doppler bins, spatial beams, and time-smoothing factor are set to be , , , and , respectively. Figure 7 depicts the complexity of ESPRIT-like algorithms against different DOFs, where the number of pulses is 2 times that of sensors (i.e., ). It is shown that the computational complexity of U-JAFE increases sharply with the increased number of sensors, whereas complexities of B-UESPRIT and BD-UESPRIT increases very slowly. The reason is that the computational complexities of B-UESPRIT and BD-UESPRIT algorithm are mainly influenced by the number of beams formed, which is unchanged with different DOFs. The computational complexity of BD-UESPRIT is slightly larger than that of B-UESPRIT due to the extra temporal DOF, but it is acceptable with the improvement of estimation accuracy.
6. Conclusions
In this paper, the unitary ESPRIT in beam-Doppler domain is developed for high-resolution DOA estimation with a single snapshot. Firstly, the technique of time-smoothing is applied to acquire snapshots. Then the conjugate centrosymmetric DFT matrix is used to transform the extracted data into beam-Doppler space. Finally, the invariance property of the spatial beam and the temporal beam is exploited, respectively, to calculate the DOA and the Doppler frequency of the target.
Experimental comparison results with the B-UESPRIT and U-JAFE methods demonstrate the superiority of the proposed method. Compared with the B-UESPRIT algorithm, the presented algorithm takes full advantage of the temporal information of the target. Therefore, the estimation accuracy of DOA is improved significantly with a slight increase in computational cost. Furthermore, the proposed method has a better DOA estimation accuracy than U-JAFE has in the case of small SNR. Moreover, the computational load of the proposed algorithm is much lower than that of U-JAFE. Several critical parameters affecting the performance are also discussed, including the choice of beam number and optimum time-smoothing factor.
Although the mathematical development of the proposed algorithm is based on the uniform linear array, it can be extended to the nonuniform array [25] case with the array interpolation technique [26].
Appendix
Derivation of (10)
From (9), it is shown that the beamspace steering vector is a real-valued vector and that the numerator of is the negative of that of (by replacing with in (9), we can get ). Then the first successive components of are related as
With the trigonometric transformation of (10), one obtains
It is worth noting that the spatial frequency of the 0th beam and the ()th beam are and , respectively. Therefore, these two beams are physically adjacent to each other. Further, we can derive
The first and last elements of are related by setting in (11), that is,
According to (A2) and (A4), all equations () lead to an identical relationship for , that is, where and are the matrices, which are defined as follows:
Thus, for all targets, the beamspace steering matrix satisfies the rotational invariant relationship as given in (10).
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported in part by the National Nature Science Foundation of China under Grants 61471285 and 61371233. This work was also supported by Postdoctoral Innovation Talent Support Program under Grant BX201700199.