# Beam-Doppler Unitary ESPRIT for Multitarget DOA Estimation

Authors: Cai Wen; Yan Zhou; Mingliang Tao; Jianxin Wu; Jinye Peng
Journal: International Journal of Antennas and Propagation
 (2018)

---

## 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.


---

## Body

## 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. Section2 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.
•T, •H, and •−1 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; IK is a K×K identity matrix; Re• and Im• 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 Figure1, 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 P statistically independent targets are located in the same range cell. The DOA and Doppler frequency of the pth (p=1,…,P) target is θp and fp, respectively. After range matched filtering, the space-time data vector of one range cell is
(1)x0=As+n,where A=a1,a2,…,aP is a KN×P steering matrix; ap=asμp⊗atνp denotes the space-time steering vector of the pth target; asμp=e−jN−1μp/21,ejμp,…,ejN−1μpT and atνp=e−jK−1νp/21,ejνp,…,ejK−1νpT are the spatial steering vector and temporal steering vector of the pth target, respectively, where μp=2πdsinθp/λ is the normalized spatial frequency, νp=2πfpTr is the normalized Doppler frequency, d is the spacing of column subarrays, λ is the operating wavelength, and Tr is the pulse repetition period; s=s1,…,sPT is a vector that consists of P complex amplitudes of targets, where sp is the complex amplitude of the pth target; and n is the KN×1 complex Gaussian white noise with zero mean and covariance matrix σ2IKN, where σ2 denotes the variance of the white noise. It is worth noting that both the spatial steering vector asμp and temporal steering vector atνp are conjugate centrosymmetric.Figure 1
Sketch map of multitarget detection for early warning radar system.In order to track multiple targets within one beam (as shown in Figure1), 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) isM (M<K). Then, a total of L=K−M+1 snapshots could be acquired by forward time-smoothing. The lth (l=1,…,L) snapshot could be expressed as
(2)xl=Jlx0,where Jl=IN⊗Jt,l is the selection matrix for choosing the lth space-time subaperture, and Jt,l=0M×l−1,IM,0M×K−l−M+1 is the M×K selection matrix for choosing the lth temporal subaperture. By substituting (1) into (2), we obtain
(3)xl=JlAs+n=A¯sl+nl,where A¯=a¯1,a¯2,…,a¯P denotes the MN×P subaperture steering matrix; a¯p=asμp⊗a¯tνp denotes the space-time subaperture steering vector, where a¯tνp=e−jM−1νp/21,ejνp,…,ejM−1νpT is the temporal subaperture steering vector of the pth target; sl=diagejl−1−K−M/2v1,ejl−1−K−M/2v2,…,ejl−1−K−M/2vPT⋅s denotes the target complex amplitude vector of the lth snapshot; and nl=Jln is the noise vector of the lth snapshot. It is shown that the temporal subaperture steering vector a¯tνp 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 Ws=ws,1,ws,2,…,ws,N and Wt=wt,1,wt,2,…,wt,M, respectively. The nth (n=1,…,N) column of Ws is
(4)ws,n=e−jN−1/2n2π/N1,ejn2π/N,…,ejN−1n2π/NT,where the nth spatial beam is formed with the spatial frequency μ=n2π/N. Similarly, the mth (m=1,…,M) column of Wt is
(5)wt,m=e−jM−1/2m2π/M1,ejm2π/M,…,ejM−1m2π/MT,where the mth temporal beam (i.e., Doppler bin) is formed with the Doppler frequency ν=m2π/M. Then, the space-time beam transformation matrix could be expressed as
(6)W=Ws⊗Wt.Thelth snapshot in beam-Doppler domain can be obtained as
(7)yl=WHA¯sl+WHnl=Bsl+n¯l,where n¯l=WHnl stands for the noise vector in beam-Doppler domain, and B denotes the MN×P space-time beamspace steering matrix, and it can be expressed as
(8)B=WHA¯=Bs⊗Bt,where Bs=bsμ1,bsμ2,…,bsμP and Bt=btν1,btν2,…,btνP are the beamspace steering matrix and Doppler-domain steering matrix, respectively, while bsμp=bs,0μp,bs,1μp,…,bs,N−1μpT and btνp=bt,0νp,bt,1νp,…,bt,M−1νpT denote the beamspace steering vector and Doppler-domain steering vector of the pth target, respectively. The nth spatial beam response of bsμp can be described as
(9)bs,nμp=ws,nHasμp=sinN/2μp−n2π/Nsin1/2μp−n2π/N.For allP targets, the beamspace steering matrix Bs satisfies a rotational invariant relationship, that is,
(10)Γ1BsΩμ=Γ2Bs.The derivation of (10) is given in the appendix, where Ωμ=diagtanμ1/2,tanμ2/2,…,tanμP/2T is a real-valued diagonal matrix whose diagonal elements contain the desired DOA information. Analogously, the rotational invariance relationship for Doppler-domain steering matrix Bt is derived as
(11)Γ3BtΩν=Γ4Bt,where Γ3 and Γ4 are defined similar to Γ1 and Γ2 with N replaced by M, and Ων=diagtanν1/2,tanν2/2,…,tanνP/2T is a real-valued diagonal matrix whose diagonal elements contain the desired Doppler information. With some manipulation and the Kronecker product property (i.e., A⊗BC⊗D=AC⊗BD), the rotational invariant relation in space-time beamspace can be derived as
(12)Γμ1BΩμ=Γμ2B,Γν1BΩν=Γν2B,where Γμ1=Γ1⊗IM and Γμ2=Γ2⊗IM are beamspace selection matrices, and Γν1=IN⊗Γ3 and Γν2=IN⊗Γ4 are Doppler-domain selection matrices.Note that a real-valued matrixY=Rey1,…,yL,Imy1,…,yL could be constructed by the beam-Doppler data. The matrix Y has a property of doubling the sample number. Performing singular-value decomposition on Y, we have
(13)Y=UΣVH,where U=u1,u2,…,uMN and V=v1,v2,…,v2L are the MN×MN and 2L×2L left singular vector matrix and right singular vector matrix, respectively, and
(14)Σ=Σ10MN×2L−MN,where Σ1=diagλ1,λ2,…,λMN, and λ1≥λ2≥⋯≥λMN are the nonzero singular values of Y. Thus, the real-valued signal subspace E can be constructed by the P “largest” left singular vectors, that is,
(15)E=u1,u2,…,uP.In theory, the signal subspaceE could be spanned by the column of B [20], that is,
(16)E=BT,where T is a P×P real-valued nonsingular matrix. Substituting B=ET−1 into (12), one gets
(17)Γμ1EΦμ=Γμ2E,Γν1EΦν=Γν2E,where Φμ=T−1ΩμT and Φν=T−1ΩνT. 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 Φμ+jΦν, that is,
(18)Φμ+jΦν=T−1Ωμ+jΩνT.That is to say,tanμp/2 and tanνp/2 (p=1,…,P) could be estimated by the real part and imaginary part of Ωμ+jΩν, respectively. Finally, the DOA and Doppler frequency of the pth target could be estimated as follows:
(19)θ^p=sin−1tan−1Reγpλπd,f^p=tan−1ImγpπTr,where γp is the pth (p=1,…,P) eigenvalue of Φμ+jΦν.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 dataxl 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 dataxl into the space-time beamspace via (7) and obtain yl.Step 4.
Construct the real-valued data matrixY=Rey1,…,yL,Imy1,…,yL and compute the signal subspace E via the P “largest” singular vectors of Y.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 Γ1, Γ2, Γ3, and Γ4.Step 6.
Perform the EVD to the complex matrixΦμ+jΦν 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 isO4KsKtMNL, where Ks and Kt denote the number of spatial beams and the number of Doppler bins, respectively. Secondly, the computational complexity of real-valued SVD is O2LKsKt2. Thirdly, the computational complexity of twice TLS is O2KsKt2P2+2P3. Finally, the computational complexity of complex EVD is O4P3. 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 Ks<N and Kt<M 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.Table 1
Comparison of computational complexities.
Algorithm
Beam transformation
Subspace estimation
TSL
Parameter estimation
B-UESPRIT [20]
O4KsNK
O2KKs2
OKs2P2+2P3
OP3
U-JAFE [12]
None
O2LMN2
O2MN2P2+2P3
O4P3
Proposed
O4KsKtMNL
O2LKsKt2
O2KsKt2P2+2P3
O4P3


### 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
(20)CRBθ=σ22RePHDθHΠA⊥DθP−RePHDθHΠA⊥DfP⋅…⋅Re−1PHDfHΠA⊥DfP⋅RePHDfHΠA⊥DθP−1,where P=diags, ΠA⊥=IKN−AAHA−1AH, Dθ=∂a1/∂θ1,…,∂aP/∂θP, Df=∂a1/∂f1,…,∂aP/∂fP, and app=1,…,P stands for the pth column of steering matrix A.


## 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) isM (M<K). Then, a total of L=K−M+1 snapshots could be acquired by forward time-smoothing. The lth (l=1,…,L) snapshot could be expressed as
(2)xl=Jlx0,where Jl=IN⊗Jt,l is the selection matrix for choosing the lth space-time subaperture, and Jt,l=0M×l−1,IM,0M×K−l−M+1 is the M×K selection matrix for choosing the lth temporal subaperture. By substituting (1) into (2), we obtain
(3)xl=JlAs+n=A¯sl+nl,where A¯=a¯1,a¯2,…,a¯P denotes the MN×P subaperture steering matrix; a¯p=asμp⊗a¯tνp denotes the space-time subaperture steering vector, where a¯tνp=e−jM−1νp/21,ejνp,…,ejM−1νpT is the temporal subaperture steering vector of the pth target; sl=diagejl−1−K−M/2v1,ejl−1−K−M/2v2,…,ejl−1−K−M/2vPT⋅s denotes the target complex amplitude vector of the lth snapshot; and nl=Jln is the noise vector of the lth snapshot. It is shown that the temporal subaperture steering vector a¯tνp 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 Ws=ws,1,ws,2,…,ws,N and Wt=wt,1,wt,2,…,wt,M, respectively. The nth (n=1,…,N) column of Ws is
(4)ws,n=e−jN−1/2n2π/N1,ejn2π/N,…,ejN−1n2π/NT,where the nth spatial beam is formed with the spatial frequency μ=n2π/N. Similarly, the mth (m=1,…,M) column of Wt is
(5)wt,m=e−jM−1/2m2π/M1,ejm2π/M,…,ejM−1m2π/MT,where the mth temporal beam (i.e., Doppler bin) is formed with the Doppler frequency ν=m2π/M. Then, the space-time beam transformation matrix could be expressed as
(6)W=Ws⊗Wt.Thelth snapshot in beam-Doppler domain can be obtained as
(7)yl=WHA¯sl+WHnl=Bsl+n¯l,where n¯l=WHnl stands for the noise vector in beam-Doppler domain, and B denotes the MN×P space-time beamspace steering matrix, and it can be expressed as
(8)B=WHA¯=Bs⊗Bt,where Bs=bsμ1,bsμ2,…,bsμP and Bt=btν1,btν2,…,btνP are the beamspace steering matrix and Doppler-domain steering matrix, respectively, while bsμp=bs,0μp,bs,1μp,…,bs,N−1μpT and btνp=bt,0νp,bt,1νp,…,bt,M−1νpT denote the beamspace steering vector and Doppler-domain steering vector of the pth target, respectively. The nth spatial beam response of bsμp can be described as
(9)bs,nμp=ws,nHasμp=sinN/2μp−n2π/Nsin1/2μp−n2π/N.For allP targets, the beamspace steering matrix Bs satisfies a rotational invariant relationship, that is,
(10)Γ1BsΩμ=Γ2Bs.The derivation of (10) is given in the appendix, where Ωμ=diagtanμ1/2,tanμ2/2,…,tanμP/2T is a real-valued diagonal matrix whose diagonal elements contain the desired DOA information. Analogously, the rotational invariance relationship for Doppler-domain steering matrix Bt is derived as
(11)Γ3BtΩν=Γ4Bt,where Γ3 and Γ4 are defined similar to Γ1 and Γ2 with N replaced by M, and Ων=diagtanν1/2,tanν2/2,…,tanνP/2T is a real-valued diagonal matrix whose diagonal elements contain the desired Doppler information. With some manipulation and the Kronecker product property (i.e., A⊗BC⊗D=AC⊗BD), the rotational invariant relation in space-time beamspace can be derived as
(12)Γμ1BΩμ=Γμ2B,Γν1BΩν=Γν2B,where Γμ1=Γ1⊗IM and Γμ2=Γ2⊗IM are beamspace selection matrices, and Γν1=IN⊗Γ3 and Γν2=IN⊗Γ4 are Doppler-domain selection matrices.Note that a real-valued matrixY=Rey1,…,yL,Imy1,…,yL could be constructed by the beam-Doppler data. The matrix Y has a property of doubling the sample number. Performing singular-value decomposition on Y, we have
(13)Y=UΣVH,where U=u1,u2,…,uMN and V=v1,v2,…,v2L are the MN×MN and 2L×2L left singular vector matrix and right singular vector matrix, respectively, and
(14)Σ=Σ10MN×2L−MN,where Σ1=diagλ1,λ2,…,λMN, and λ1≥λ2≥⋯≥λMN are the nonzero singular values of Y. Thus, the real-valued signal subspace E can be constructed by the P “largest” left singular vectors, that is,
(15)E=u1,u2,…,uP.In theory, the signal subspaceE could be spanned by the column of B [20], that is,
(16)E=BT,where T is a P×P real-valued nonsingular matrix. Substituting B=ET−1 into (12), one gets
(17)Γμ1EΦμ=Γμ2E,Γν1EΦν=Γν2E,where Φμ=T−1ΩμT and Φν=T−1ΩνT. 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 Φμ+jΦν, that is,
(18)Φμ+jΦν=T−1Ωμ+jΩνT.That is to say,tanμp/2 and tanνp/2 (p=1,…,P) could be estimated by the real part and imaginary part of Ωμ+jΩν, respectively. Finally, the DOA and Doppler frequency of the pth target could be estimated as follows:
(19)θ^p=sin−1tan−1Reγpλπd,f^p=tan−1ImγpπTr,where γp is the pth (p=1,…,P) eigenvalue of Φμ+jΦν.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 dataxl 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 dataxl into the space-time beamspace via (7) and obtain yl.Step 4.
Construct the real-valued data matrixY=Rey1,…,yL,Imy1,…,yL and compute the signal subspace E via the P “largest” singular vectors of Y.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 Γ1, Γ2, Γ3, and Γ4.Step 6.
Perform the EVD to the complex matrixΦμ+jΦν 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 isO4KsKtMNL, where Ks and Kt denote the number of spatial beams and the number of Doppler bins, respectively. Secondly, the computational complexity of real-valued SVD is O2LKsKt2. Thirdly, the computational complexity of twice TLS is O2KsKt2P2+2P3. Finally, the computational complexity of complex EVD is O4P3. 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 Ks<N and Kt<M 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.Table 1
Comparison of computational complexities.
Algorithm
Beam transformation
Subspace estimation
TSL
Parameter estimation
B-UESPRIT [20]
O4KsNK
O2KKs2
OKs2P2+2P3
OP3
U-JAFE [12]
None
O2LMN2
O2MN2P2+2P3
O4P3
Proposed
O4KsKtMNL
O2LKsKt2
O2KsKt2P2+2P3
O4P3


## 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
(20)CRBθ=σ22RePHDθHΠA⊥DθP−RePHDθHΠA⊥DfP⋅…⋅Re−1PHDfHΠA⊥DfP⋅RePHDfHΠA⊥DθP−1,where P=diags, ΠA⊥=IKN−AAHA−1AH, Dθ=∂a1/∂θ1,…,∂aP/∂θP, Df=∂a1/∂f1,…,∂aP/∂fP, and app=1,…,P stands for the pth column of steering matrix A.


## 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
(21)Rn=∑l=1Ln¯ln¯lH.Rn 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 Rn, one gives Rn=UnΣnUnH, where Σn is a diagonal matrix that consists of eigenvalues of Rn and Un is the corresponding eigenvector unitary matrix. Then the whitened beamspace data can be expressed as [24]
(22)y¯l=UnΣn−1/2UnHyl.Assume that the signal subspace estimated by the whitened data isE¯. It is easy to derive that the relationship between the signal subspace with and without whitening is
(23)E=UnΣn1/2UnE¯.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.


## 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
(21)Rn=∑l=1Ln¯ln¯lH.Rn 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 Rn, one gives Rn=UnΣnUnH, where Σn is a diagonal matrix that consists of eigenvalues of Rn and Un is the corresponding eigenvector unitary matrix. Then the whitened beamspace data can be expressed as [24]
(22)y¯l=UnΣn−1/2UnHyl.Assume that the signal subspace estimated by the whitened data isE¯. It is easy to derive that the relationship between the signal subspace with and without whitening is
(23)E=UnΣn1/2UnE¯.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 ofd=λ/2. We assume that two far-field signals with equal power s1 and s2 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 s1 are 25.9° and 463 Hz, and those of s2 are 34.1° 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 Nm=500 Monte Carlo runs. The evaluation metric, that is, root-mean-square error (RMSE) of DOA, is defined as
(24)RMSEθ=1PNm∑p=1P∑l=1Nmθ^p,l−θp2,where θ^p,l 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. Section4.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 Ks=4 and Kt=3, 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].Figure 2
Effect of whitening on the DOA estimation error.


### 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 beM=0.5K and Ks=4, 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 Kt=3 in our method.Figure 3
Effect of the number of Doppler bins on the DOA estimation error.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 Figure4. In this simulation, the time-smoothing factor and the number of Doppler bins are chosen to be M=0.5K and Kt=3, 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 mth beam, and then the target information is mainly included in the (m−1)th, mth, and (m+1)th beams).Figure 4
Effect of the number of spatial beams on the DOA estimation error.


### 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 beKt=3 and Ks=4, 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 M≈0.4K. 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.Figure 5
Parameterized DOA estimation error as a function of time-smoothing factorM.


### 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 Kt=3, Ks=4, and M=0.4K, 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 Ks=4. 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 Figure6. 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.Figure 6
Variation of the DOA estimation errors with SNR for various methods.


### 5.5. Computational Load Comparison with the ESPRIT-like Algorithms

In this subsection, computational complexities of the method mentioned in Table1 are analyzed. The number of targets, Doppler bins, spatial beams, and time-smoothing factor are set to be P=3, Kt=3, Ks=4, and M=0.4K, 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., K=2N). 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.Figure 7
Complexity comparison against number of sensorsN (number of pulses K=2N).


## 5.1. Investigating the Effectiveness of Prewhitening

Results comparing the performances of the proposed algorithm implemented with and without prewhitening (viz. Section4.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 Ks=4 and Kt=3, 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].Figure 2
Effect of whitening on the DOA estimation error.


## 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 beM=0.5K and Ks=4, 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 Kt=3 in our method.Figure 3
Effect of the number of Doppler bins on the DOA estimation error.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 Figure4. In this simulation, the time-smoothing factor and the number of Doppler bins are chosen to be M=0.5K and Kt=3, 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 mth beam, and then the target information is mainly included in the (m−1)th, mth, and (m+1)th beams).Figure 4
Effect of the number of spatial beams on the DOA estimation error.


## 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 beKt=3 and Ks=4, 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 M≈0.4K. 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.Figure 5
Parameterized DOA estimation error as a function of time-smoothing factorM.


## 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 Kt=3, Ks=4, and M=0.4K, 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 Ks=4. 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 Figure6. 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.Figure 6
Variation of the DOA estimation errors with SNR for various methods.


## 5.5. Computational Load Comparison with the ESPRIT-like Algorithms

In this subsection, computational complexities of the method mentioned in Table1 are analyzed. The number of targets, Doppler bins, spatial beams, and time-smoothing factor are set to be P=3, Kt=3, Ks=4, and M=0.4K, 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., K=2N). 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.Figure 7
Complexity comparison against number of sensorsN (number of pulses K=2N).


## 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].



---

*Source: 3568286-2018-03-26.xml*