# DOA Tracking of Two-Dimensional Coherent Distribution Source Based on Fast Approximated Power Iteration

Authors: Tao Wu; Pengtao Zhang; Yiwen Li; Yangjun Gao; Chaoqi Fu; Bo Feng
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Aiming at two-dimensional (2D) coherent distributed (CD) sources, this paper has proposed a direction of arrival (DOA) tracking algorithm based on signal subspace updating under the uniform rectangular array (URA). First, based on the hypothesis of small angular spreads of distributed sources, the rotating invariant relations of the signal subspace of the receive vector of URA are derived. An ESPRIT-like method is constructed for DOA estimation using two adjacent parallel linear arrays of URA. Through the synthesis of estimation by multiple groups of parallel linear arrays within URA arrays, the DOA estimation method for 2D CD sources based on URA is obtained. Then, fast approximated power iteration (FAPI) subspace tracking algorithm is used to update the signal subspace. In this way, DOA tracking of 2D CD sources can be realized by DOA estimation through signal subspace updating. This algorithm has a low computational complexity and good real-time tracking performance. In addition, the algorithm can track multiple CD sources without knowing the angular signal distribution functions, which is robust to model errors.


---

## Body

## 1. Introduction

In the field of underwater acoustic and radar detection, since the position of target is constantly changing, it is necessary to obtain the target orientation information in real time. The traditional DOA tracking methods assume that the target is a point source. Considering stationary targets based on the point source model, some important achievements have been made in the field of mixed signal and sparse array and multiangle estimation recently. Based on the uniform rectangular array (URA), the two-dimensional (2D) direction of arrival (DOA) estimation of circular and noncircular mixed signals is proposed in [1]. In [2], authors have proposed a reduced dimension noncircular MUSIC algorithm for DOA and polarization estimation of circular and noncircular mixed signals in the virtual multiple input multiple output (MIMO) system. The authors of [3] have proposed two kinds of extended coprime arrays, namely, sliding extended coprime array and coprime array with displaced subarrays, which can alleviate mutual coupling and increased degree of freedom significantly. Based on a coprime array, authors of [4] have proposed an estimator where the Toeplitz matrix functions is utilized to resolve off-grid sources. The authors of [5] have realized the estimation of 2D DOA, receiving and transmitting angle simultaneously in the electromagnetic MIMO system. The authors of [6] have presented direction of departure and DOA estimation for the MIMO radar where a reduced dimension multiple signal classification algorithm requiring one-dimension search is derived. The authors of [7] have proposed an estimator for DOA and mutual coupling self-calibration in the uniform linear arrays-based bistatic MIMO radar where a two-step framework is proposed for DOA estimation, and mutual coupling coefficients are obtained via the least square method. In [8], authors have proposed an estimator for joint 2D DOD, 2D DOA, and polarization angles without pairing in the bistatic MIMO system with electromagnetic vector sensors.There are two main types of DOA tracking methods based on the point source model. One is a mature method based on the Kalman filter [9, 10], which has been used to track the parameters of DOA, target distance, and so on. The other is based on signal subspace updating methods, such as projection approximation subspace tracking (PAST) and PAST based on the deflation technique (PASTd) [11–13], orthonormal projection approximation subspace tracking (OPAST) [14, 15], and fast approximated power iteration (FAPI) [16, 17]. In the field of underwater acoustic detection, due to the existence of multipath propagation between the receiving array and target, especially with the reduction of the distance between the target and receiving array, the spatial scattering characteristics of the target cannot be ignored, and the DOA estimation or tracking algorithm based on the point source model is no longer established. Under this context, scholars put forward the distributed source models [18].According to the time correlation of scattering elements, distributed sources can be divided into incoherent and coherent sources. The incoherent distributed (ID) source assumes that different scattering elements of a target are incoherent, and the coherent distributed (CD) source assumes that different scattering elements of a target are coherent. According to the spatial distribution dimension, the distributed sources can be classified into one-dimensional (1D) and two-dimensional (2D) distributed sources. The 2D distributed sources assume that the target and receiving array are in a three-dimensional space, which is more general. For CD sources, the spatial distribution function describing the scattering elements of the target is a deterministic angular signal distribution function (ASDF), which is a probability density function depending on the distribution of scattering points in the space.As the distributed source model not only includes all the parameters of the point source but also adds the angular spread parameters, the complexity of parameter estimation increases significantly. For 2D ASDF function, it can be modeled as Gaussian, uniform, or other function forms. The parameters included are the nominal azimuth, nominal elevation, azimuth spreads, and elevation spreads. Nominal azimuth and nominal elevation can be collectively called DOA reflecting direction of arrival of the target center.Considering the targets are static, scholars have presented different algorithms to model and solve DOA estimation of distributed sources under diverse arrays [19–30]. Nevertheless, in the background of a moving target, the DOAs of distributed sources are time-varying. The research studies on DOA tracking of distributed sources are relatively fewer. Based on the covariance matrix renewed by the exponential window function, literature [31] has presented a DOA tracking method of 1D CD sources where the signal subspace is updated by FAPI; then, DOA is resolved by TLS-ESPRIT. Authors of [32] have presented a DOA tracking model based on the support vector machine (SVM) for 1D CD sources where vectorizations of the covariance matrix renewed by the exponential window function are regarded as the input of SVM, and DOAs are used as the output. In literature [33], an adaptive DOA tracking method for 1D ID sources is proposed where DOAs are estimated by the covariance fitting technique using the central moments of the angular power density; then, utilizing the constant acceleration model, the DOA tracking model is established in the framework of the Kalman filter.This paper presents a DOA tracking algorithm for 2D CD sources based on URA. First, the rotating invariant relations of signal subspace are derived by exploring the URA structure under the premise of small angular spreads of sources. Then, the signal subspace is updated by FAPI tracking algorithm in real time. Based on the rotation invariant relations of the signal subspace, DOAs of 2D CD sources can be calculated. This algorithm has a low computational complexity and good real-time tracking performance. In addition, the algorithm can track multiple distributed sources without knowing the ASDF of distributed sources and shows robustness to model errors. The main contributions of this article are listed below.(i)
The rotation invariance relations within the URA with respect to nominal elevation and nominal azimuth are deduced under small angular spreads assumption(ii)
The 2D DOA estimation algorithm of CD sources is derived based on URA, which does not need to know the ASDF form of sources and need not spectral peak search and can be applied for DOA estimation for point sources under the background of the nonmoving target(iii)
Based on the FAPI framework, a DOA tracking method is proposed. The algorithm has good tracking accuracy and can track multiple 2D CD sources in real time. It is also suitable for DOA tracking of point sources and has robustness in the case of mismatch between the point source and CD source model.The structure of the paper is as follows. Part 2 introduces the signal model and the receive vectors of arrays proposed. In part 3, the solution approach of DOA based on URA is detailed. Part 4 elaborates DOA tracking based on signal subspace iteration. Part 5 illustrates the simulation results. Part 6 draws the conclusion.Notations: scalar variables are denoted by italic letters, and vectors and matrices are denoted by bold letters.⋅−1, ⋅T, and ⋅H mean inverse, transpose, and Hermitian transpose of a matrix. E⋅, (./), and ⋅+ denote expectation, elementwise division, and pseudoinverse operation; diag⋅ is the diagonal matrix; arg⋅ is the phase of a complex number.


## 2. Array Structure and Signal Model

The structure of the array is shown in Figure1. The URA is on the xoz plane. The distance between sensors along the x-axis and z-axis is d. There are M sensors along the x-axis and K sensors along the z-axis, and the array contains MK sensors totally. It is assumed that there are q 2D CD sources with wavelength λ in the far-field space with DOAs as (θi, φi) (i = 1, 2, …, q) impinging on the array. θi is the nominal azimuth of the ith source, and φi is the nominal elevation, θi∈[0, π], φi∈[0, π]. It is assumed that the noise of each sensor is additive Gaussian white noise, which is not related to the signal. Assuming that the origin sensor position is m = k = 1, the signal received by the (m, k)th sensor in URA can be expressed as(1)ymkt=∑i=1q∬ejπm−12dcosθsinϕ/λejπk−12dcosϕ/λsiθ,ϕ,tdθdϕ+nmkt,where nmk(t) is the noise received. As for CD sources, the angular signal density function can be written as(2)siθ,ϕ,t=sitgiθ,ϕ;ui,where si (t) reflects time characteristics of the ith CD sources, gi (θ, φ; ui) is the deterministic angular signal density function (ASDF) of the ith source, and σi2=sit2 is the power of the source. Define signal vector s (t) as follows:(3)st=s1t,s2t,…,sqtT.Figure 1
Uniform rectangular array structure.Assumed different CD sources are incoherent, we have(4)Λ=EststH=diagσ12,σ22,…,σq2.The generalized manifold coefficient of the (m, k)th sensor for the ith source is defined as(5)ηmkθi,ϕi=∬ej2πdm−1cosθsinϕ/λej2πdk−1cosϕ/λgiθ,ϕ;uidθdϕ.Reflecting the response to all CD sources, the generalized manifold vector of the (m, k)th sensor can be expressed as(6)amkθ,ϕ=ηmkθ1,ϕ1,ηmkθ2,ϕ2,…,ηmkθq,ϕq.The signal received by the (m, k)th sensor can be expressed as(7)ymkt=amkθ,ϕst+nmkt.Under the premise of small angular spreads of sources and the distance between sensors,d is set as the half of the wavelength λ, generalized manifold coefficients of the (m, k)th, the (m−1, k)th, and the (m, k−1)th, and the (m−1, k−1)th sensors have the rotating relations as follows (Appendix):(8)ηmkθi,ϕi≈ejπcosθisinϕiηm−1kθi,ϕi,(9)ηmkθi,ϕi≈ejπcosϕiηmk−1θi,ϕi,(10)ηmkθi,ϕi≈ejπcosθisinϕiejπcosϕiηm−1k−1θi,ϕi.


## 3. DOA Estimation under URA

In order to track the dynamic target in real time, DOA estimation for the 2D CD source based on URA should be derived based on the signal subspace, which is a basis of DOA tracking. In this part, according to the rotating relations (8–10), an ESPRIT-like DOA estimation method for the 2D CD signal source is constructed between two adjacent parallel linear arrays of URA. Through the synthesis solutions of multiple groups of parallel linear arrays within URA arrays, the 2D DOA estimation method based on URA is obtained.M sensors on the x-axis are defined as subarray X1, and the received signals of subarray X1 can be expressed as(11)x1t=A1st+nx1t,where A1 is the generalized manifold matrix of subarray X1, which has the following expression:(12)A1=a11θ,ϕa21θ,ϕ⋯aM1θ,ϕ=η11θ1,ϕ1η11θ2,ϕ2⋯η11θq,ϕqη21θ1,ϕ1η21θ2,ϕ2⋯η21θ2,ϕ2⋯⋯⋯⋯ηM1θ1,ϕ1ηM1θ1,ϕ1⋯ηM1θq,ϕq.Define subarrayXk as the linear array parallel to the x-axis and with the interval of (k − 1) d from subarray X1. It can be concluded that the received signal of subarray Xk is(13)xkt=Akst+nkt,where the noise vector can be expressed as(14)nkt=n1,k−1t,n2,k−1t,…,nM,k−1t.Take the firstM − 1 sensors of subarray Xk to form subarray Xk1 and the last M − 1 sensors to form subarray Xk2, and the generalized manifold matrix of subarray Xk1 can be expressed as(15)Ak1=IM−1×M−1|0M−1×1Ak.Generalized manifold matrix of subarrayXk2 can be expressed as(16)Ak2=0M−1×1|IM−1×M−1Ak.According to equations (8)–(10), generalized manifold matrices have the following rotating invariant relations:(17)Ak2≈Ak1Φx,(18)Ak1≈Ak−11Φz,where rotational operations have the following expressions:(19)Φx=diagejπcosθ1sinϕ1,ejπcosθ2sinϕ2,…,ejπcosθqsinϕq,Φz=diagejπcosϕ1,ejπcosϕ2,…,ejπcosϕq.Receive vectors of subarraysXk1 and Xk2 can be written as(20)xk1t≈Ak1st+nk1t,xk2t≈Ak2st+nk2t.Receive vectors of URA can be written as(21)Xt=x1tx2t⋯xKt=Ast+nXt,where A is the generalized manifold matrix of URA, which have the following expression:(22)A=A1A2⋯AK.Receive noise vector of URA can be written as(23)nXt=n1tT,n2tT,…,nKtTT,where(24)EnXtnXHt′=σn2IMKδt−t′,EnXtnXTt′=0,∀t,t′.The covariance matrix ofX(t) can be written as(25)RZ=EXtXtH=AΛAH+σn2IMK.The signal subspace W composed of eigenvectors corresponding to the largestq eigenvalues can be obtained by eigendecomposition of RZ. The space of W is the same as that of A. Therefore, a q × q dimensional nonsingular matrix F satisfies the following relation:(26)W=AF=A1A2⋯AKF=B1⋮Bς⋮BK/2F,where(27)B1=A11the mth row ofA1A21the2mth row ofA2=the mth row ofA1A12A21the2mth row ofA2,Bς=A2ς−11the mth row ofA2ς−1A2ς1the2ςth row ofA2ς=the1ςth row ofA2ς−1A2ς−12A2ς1the2ςth row ofA2ς.DefineW11 as a vector selecting the first M-1 rows of W and W12 as the 2nd to Mth rows of W. W2ς−11 is a vector selecting [(ς –1)2M + 1]th to [(ς –1)2M + M − 1]th rows of W, W2ς−12 is a vector selecting [(ς –1)2M + 2]th to [(ς –1)2M + M] of rows W, and ς = 1, 2, …, K/2; K is an even number. According to equation (26), W2ς−11, W2ς−12, and W2ς1 have the following relations:(28)W2ς−11=A2ς−11F,W2ς−12=A2ς−12F,W2ς1=A2ς1F.According to the rotating invariant relations described by equations (17) and (18), we have(29)W2ς−12W2ς−11+A2ς−11=A2ς−11ΦX,W2ς1W2ς−11+A2ς−11=A2ς−11ΦZ.It is known that the rotation operators Фx and Фz are the diagonal matrices. The corresponding eigenvectors of the elements in the same position of the diagonal matrices Фx and Фz are the same. Therefore, the eigenvalues ejπcosϕi of W2ς1W2ς−11+ and eigenvalues ejπcosθisinϕi of W2ς−12W2ς−11+ have the same eigenvector.Eigenvalueejπcosϕi can be obtained by eigendecomposition W2ς1W2ς−11+. Then, nominal elevation can be obtained as(30)ϕiς=arccosargμiςπ,i=1,2,…,q,where μiς is the ith eigenvalue of W2ς1W2ς−11+, and superscript ς of μiς denotes the eigenvalue calculated by W2ς−11 and W2ς1. The corresponding eigenvector of μiς is ξiς. Then, the eigenvalue νiς of W2ς−12W2ς−11+, which has the same position as μiς, can be expressed as(31)νiς=1q11×q⋅W2ς−12W2ς−11+ξiς.ξiς,where 11 ×q is the 1 × q dimensional vector with all the elements being 1; then, the nominal azimuth of the distribution source can be obtained as follows:(32)θiς=arccosargνiςπsinϕiς,i=1,2,…,q.It can be concluded that the nominal elevationϕi and the nominal azimuth θi are the mean values obtained from the K/2 group of W2ς−11, W2ς−12, and W2ς1.(33)ϕi=2K∑ς=1K/2ϕiς,(34)θi=2K∑ς=1K/2θiς.


## 4. DOA Tracking Based on FAPI

On the basis of obtaining the signal subspace, the DOA tracking of 2D CD signal sources can be realized by the above approach proposed. In this part, the signal subspace is updated in real time by FAPI algorithm under URA; then, DOA is solved by the estimation algorithm in the previous part. In this way, real-time tracking of DOA of 2D CD sources can be realized.For the tracking problem, the sample covariance matrix of snapshot dataX(t) (t = 1, 2, …, N) can be recursively updated in form of an exponential window as follows:(35)Rt=αRt−1+XtXtH,where a > 0 is the forgetting factor. Define W(t) is the signal subspace of the sample covariance matrix at time t. Then, the power iteration method can be written follows:(36)Wt=RtWt−1.In order to ensure that the modulus of each vector in the signal subspace is 1, multiply the right end of equation (36) by a coefficient matrix to obtain the following relationship:(37)Wt=RtWt−1Wt−1HR2t−1Wt−1−1/2.Equation (37) can guarantee the orthogonality of W(t) in the iterative process, that is to say(38)WtHWt=Iq.Authors of [16] have proved that expression (38) converges to the signal subspace of the covariance matrix.Define compressed signal vector as(39)Vt=Wt−1HXt.Then, equation (37) can be expressed as(40)WtΨt=RWt,where(41)Rwt=RtWt−1.Rw(t) can be regarded as the cross covariance matrix of the received signal vector X(t) and the compressed signal vector V(t). Ψt is the square root matrix of Rw(t), which have the following expression:(42)Ψt=RWtHRWt1/2.The power iterative method can be composed of the data compression process (39) and orthogonalization process (40). Although the power iterative method is very effective, the calculation complexity of the orthogonal process is high, which is not conducive to real-time application. This classical iterative subspace estimation method often cannot guarantee the good convergence and robust performance of the results.FAPI algorithm is proposed in [17], which are optimized and solved on the basis of the power iteration method. In order to track the signal subspace W(t), the detailed steps of FAPI algorithm are summarized in Table 1.Table 1
Pseudocode of the FAPI method.
Input:Xtt=1N, which is N snapshots of receive vectors of URA.Output: W(t), which is the signal subspace of the covariance matrix of X(t).Initialize: W0=Iq×q,0MK−q×qT,L0=Iq×q.Cost1.Vt=Wt−1HXtMKq2.gt=Lt−1Vtα+VtHLt−1HVt−1q23.γt2=XtHXt−VtHVtMK4.βt=1+γt2gt2q5.τt=γt2βt+βt1/2−116.ξt=1−τtgt217.X′t=Vtξt+gtτt2q8.ςt=τtηt−1Lt−1gt−gtX′tHgt2q29.Lt=α−1Lt−1−gtX′tHLt−1+ςtgtH2q210.e′t=Xtξt−Wt−1gtτtMKq + MK11.Wt=Wt−1+e′tgtHMKqThe proposed 2D DOA fast tracking algorithm steps are listed as follows.(1)
InitializeW(0) and L(0)(2)
Update snapshot dataX(t) and obtain the signal subspace W(t) at time t by the FAPI algorithm described in Table 1(3)
Divide the subspace intoK/2 groups of W2ς−11, W2ς−12, and W2ς1(4)
Calculate the nominal elevation from the eigenvalue obtained by eigendecomposition using equation (30) and calculate the azimuth angle using equations (31) and (32)(5)
Obtained the mean value by equations (33) and (34) as the estimated elevation and azimuthIt can be seen from table that the computational complexity required for each iteration of signal subspace updating is O (MKq), which is a very low value. After the signal subspace is obtained, the complexity of DOA estimation mainly lies in the eigendecomposition of RZ, which is OMK3, and calculate a pseudoinverse of W2ζ−11, which is OMK3.


## 5. Results and Discussion

In this section, three experiments are conducted to investigate the tracking performance of the proposed algorithm. The array model shown in Figure1 is used. The array geometric parameters are M = K = 4, and d = 0.5λ. The forgetting factor α = 0.95. It is assumed that the angular spread parameters remain unchanged during the movement of distributed sources.Experiment 1 investigates the influence of the signal-to-noise ratio (SNR) on DOA tracking of single CD source and effectiveness of the method proposed. The initial DOA parameters are (30°, 45°) and transformed to (70°, 85°) during 2000 snapshots uniformly. Figure2 shows the tracking curves by PASTd [11] with SNR = 0 dB, 10 dB, and 20 dB. Figure 3 shows the tracking curves of the algorithm proposed when the angular spread is 2°. It can be seen that with the increase of SNR, the tracking effect of both algorithms is improved. The experimental results show that the proposed algorithm has a good tracking performance for the single 2D CD source. It is worth noting that even with the increase of SNR, both the estimated nominal elevation and nominal azimuth have a stable trend error with the passage of time in Figure 2, that is, the estimated nominal elevation is less than the true value and the estimated nominal azimuth is greater than the real value. Since PASTd [11] is for DOA tracking of the point source, when the target model is a distributed source, there is a large error due to model mismatch.Figure 2
DOA tracking of a 2D CD source by PASTd for the point source.
(a)(b)(c)Figure 3
DOA tracking of a 2D CD source by the method proposed.
(a)(b)(c)Experiment 2 examines the effect of angular spread on DOA tracking of the single CD source. The changing process of DOA is the same as experiment 1. Figure4 shows tracking curves with angular spreads set as 2°, 5°, and 10° when SNR = 10 dB. It can be seen that the tracking effect decreases with the increase of the angular spreads. The derivation of the rotating invariant relations in the DOA estimation stage is based on the assumption of the small angular spreads, so the increase of the angular spread will bring errors in the DOA estimation step. As a matter of fact, whether it is Gaussian or a uniform distributed source or any other distributed sources, when the angular spreads are 0°, the parameters describing the CD source only include DOA parameters containing azimuth and elevation, which means the distributed source can be considered as a point source. This experiment also shows that the method proposed in this paper is not only suitable for DOA tracking of the CD source but also suitable for DOA tracking of the point source. In the case of mismatch between the point source and CD source model, it has good robustness.Figure 4
DOA tracking of a 2D CD source with different angular spreads.
(a)(b)(c)Experiment 3 investigates DOA tracking effectiveness for multiple sources. Assuming that three CD sources are incident on the array, all the angular spreads of the sources are 2°. The ASDF of the first CD source obeys Gaussian distribution, while the second and the third CD source obey the uniform distribution. The first source is transformed from (30°, 45°) to (70°, 25°) in 2000 snapshots. The second and third sources are transformed from (40°, 45°) (50°, 45°) to (80°, 25°) (90°, 25°) in 2000 snapshots, respectively. Figure5 shows tracking curves of three sources based on the method proposed. From the derivation process in parts 3 and 4 and experimental results, we can see that the method proposed in this paper does not need to know the specific form of ASDF of sources during the DOA tracking process. The results show that the algorithm can effectively track the DOA of multiple distribution sources without ASDF information.Figure 5
DOA tracking of multiple CD sources by the method proposed.Experiment 4 investigates the tracking effect of the proposed algorithm under different change rates when the distributed source changes uniformly. All the CD sources are Gaussian type, the angular spreads are all 2°, and SNR = 10 dB. The transformation processes are in 2000 snapshots. Four kinds of transformation are investigated. The first case is that the source changes from (30°, 45°) to (55°, 70°); the second case changes from (30°, 45°) to (70°, 85°); the third case changes from (30°, 45°) to (85°, 100°); and the fourth case changes from (30°, 45°) to (100°, 115°). It can be seen from Figure6 that with the increase of the change rates, the tracking accuracy will slightly deteriorate; the FAPI algorithm is based on the updating of the sample covariance matrices in the form of an exponential window, which is the approximation of the covariance matrices of moving targets. But on the whole, the tracking effect is still satisfactory.Figure 6
DOA tracking under different change rates.
(a)(b)(c)(d)


## 6. Conclusions

In this paper, a fast tracking algorithm based on subspace updating is proposed for moving 2D CD sources under URA. The signal subspace is updated by fast approximated power iteration algorithm; the nominal azimuth and nominal elevation are calculated by an ESPRIT-like algorithm through the estimation synthesis of multiple groups of parallel linear arrays within URA arrays. The algorithm proposed in this paper has low complexity and is suitable for real-time tracking. Experimental results show that the algorithm proposed in this paper has high tracking accuracy and robustness in case of mismatch between the point source and CD source model and can track multiple distribution sources without knowing the angular signal distribution functions.



---

*Source: 3219516-2020-11-16.xml*