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

. 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 al-gorithm can track multiple CD sources without knowing the angular signal distribution functions, which is robust to model errors.


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. e 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. e 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. e authors of [5] have realized the estimation of 2D DOA, receiving and transmitting angle simultaneously in the electromagnetic MIMO system. e 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. e 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.
ere 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. e other is based on signal subspace updating methods, such as projection approximation subspace tracking (PAST) and PAST based on the deflation technique (PASTd) [11][12][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. e 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. e 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. e 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][20][21][22][23][24][25][26][27][28][29][30]. Nevertheless, in the background of a moving target, the DOAs of distributed sources are time-varying. e 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.
is 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.
en, 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. is 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. e main contributions of this article are listed below.
(i) e rotation invariance relations within the URA with respect to nominal elevation and nominal azimuth are deduced under small angular spreads assumption (ii) e 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. e 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. e 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.

Array Structure and Signal Model
e structure of the array is shown in Figure 1. e URA is on the xoz plane. e distance between sensors along the xaxis and z-axis is d. ere 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 i th source, and φ i is the nominal elevation, 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 where n mk (t) is the noise received. As for CD sources, the angular signal density function can be written as where s i (t) reflects time characteristics of the i th CD sources, g i (θ, φ; u i ) is the deterministic angular signal density function (ASDF) of the i th source, and σ 2 i � ‖s i (t)‖ 2 is the power of the source. Define signal vector s (t) as follows: Assumed different CD sources are incoherent, we have e generalized manifold coefficient of the (m, k) th sensor for the i th source is defined as Reflecting the response to all CD sources, the generalized manifold vector of the (m, k) th sensor can be expressed as e signal received by the (m, k) th sensor can be expressed as 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):

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)(9)(10), an ESPRIT-like DOA estimation method for the 2D CD signal source is constructed between two adjacent parallel linear arrays of URA. rough 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 X 1 , and the received signals of subarray X 1 can be expressed as where A 1 is the generalized manifold matrix of subarray X 1 , which has the following expression: Define subarray X k as the linear array parallel to the xaxis and with the interval of (k − 1) d from subarray X 1 . It can be concluded that the received signal of subarray X k is where the noise vector can be expressed as Take the first M − 1 sensors of subarray X k to form subarray X 1 k and the last M − 1 sensors to form subarray X 2 k , and the generalized manifold matrix of subarray X 1 k can be expressed as Generalized manifold matrix of subarray X 2 k can be expressed as According to equations (8)-(10), generalized manifold matrices have the following rotating invariant relations: where rotational operations have the following expressions: Mathematical Problems in Engineering Φ x � diag e jπ cos θ 1 sin ϕ 1 , e jπ cos θ 2 sin ϕ 2 , . . . , e jπ cos θ q sin ϕ q , Φ z � diag e jπ cos ϕ 1 , e jπ cos ϕ 2 , . . . , e jπ cos ϕ q .
Receive vectors of subarrays X 1 k and X 2 k can be written as Receive vectors of URA can be written as where A is the generalized manifold matrix of URA, which have the following expression: Receive noise vector of URA can be written as where e covariance matrix of X(t) can be written as e signal subspace W composed of eigenvectors corresponding to the largest q eigenvalues can be obtained by eigendecomposition of R Z . e space of W is the same as that of A. erefore, a q × q dimensional nonsingular matrix F satisfies the following relation: where According to the rotating invariant relations described by equations (17) and (18), we have It is known that the rotation operators V x and V z are the diagonal matrices. e corresponding eigenvectors of the elements in the same position of the diagonal matrices V x and V z are the same. erefore, the eigenvalues e jπ cos ϕ i of W 1 2ς (W 1 2ς−1 ) + and eigenvalues e jπ cos θ i sin ϕ i of W 2 2ς−1 (W 1 2ς−1 ) + have the same eigenvector.
Eigenvalue e jπcos ϕ i can be obtained by eigendecomposition W 1 2ς (W 1 2ς−1 ) + . en, nominal elevation can be obtained as  Mathematical Problems in Engineering where μ ς i is the i th eigenvalue of W 1 2ς (W 1 2ς−1 ) + , and superscript ς of μ ς i denotes the eigenvalue calculated by W 1 2ς−1 and W 1 2ς . e corresponding eigenvector of μ ς i is ξ ς i . en, the eigenvalue ] ς i of W 2 2ς−1 (W 1 2ς−1 ) + , which has the same position as μ ς i , can be expressed as where 1 1 × 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: 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 W 1 2ς−1 , W 2 2ς−1 , and W 1 2ς .

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 data X(t) (t � 1, 2, . . ., N) can be recursively updated in form of an exponential window as follows: where a > 0 is the forgetting factor. Define W(t) is the signal subspace of the sample covariance matrix at time t. en, the power iteration method can be written follows: 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: Equation (37) can guarantee the orthogonality of W(t) in the iterative process, that is to say (38) Authors of [16] have proved that expression (38) converges to the signal subspace of the covariance matrix.
Define compressed signal vector as en, equation (37) can be expressed as where R w (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 R w (t), which have the following expression: e 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. is 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. e proposed 2D DOA fast tracking algorithm steps are listed as follows.
(1) Initialize W(0) and L(0) (2) Update snapshot data X(t) and obtain the signal subspace W(t) at time t by the FAPI algorithm described in Table 1 (3) Divide the subspace into K/2 groups of W 1 2ς−1 , W 2 2ς−1 , and W 1 2ς (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 azimuth It 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 R Z , which is O[(MK) 3 ], and calculate a pseudoinverse of W 1 2ζ−1 , which is O[(MK) 3 ]. { } N t�1 , which is N snapshots of receive vectors of URA. Output: W(t), which is the signal subspace of the covariance matrix of X(t). Tracking of θ by PASTd [11] Tracking of ϕ by PASTd [11] The ture trajectory Tracking of θ by PASTd [11] Tracking of ϕ by PASTd [11] The ture trajectory Tracking of θ by PASTd [11] Tracking of ϕ by PASTd [11] The ture trajectory (c) Figure 2: DOA tracking of a 2D CD source by PASTd for the point source.

Results and Discussion
In this section, three experiments are conducted to investigate the tracking performance of the proposed algorithm. e array model shown in Figure 1 is used. e array geometric parameters are M � K � 4, and d � 0.5 λ. e 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-tonoise ratio (SNR) on DOA tracking of single CD source and effectiveness of the method proposed. e initial DOA parameters are (30°, 45°) and transformed to (70°, 85°) during 2000 snapshots uniformly. Figure 2 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. e 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. Experiment 2 examines the effect of angular spread on DOA tracking of the single CD source. e changing process of DOA is the same as experiment 1. Figure 4 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. e 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. is 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. 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°. e ASDF of the first CD source obeys Gaussian distribution, while the second and the third CD source obey the uniform distribution. e first source is transformed from (30°, 45°) to (70°, 25°) in 2000 snapshots. e second and third sources are transformed from (40°, 45°) (50°, 45°) to (80°, 25°) (90°, 25°) in 2000 snapshots, respectively. Figure 5 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   Mathematical Problems in Engineering method proposed in this paper does not need to know the specific form of ASDF of sources during the DOA tracking process. e results show that the algorithm can effectively track the DOA of multiple distribution sources without ASDF information. 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. e transformation processes are in 2000 snapshots. Four kinds of transformation are investigated. e 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 Figure 6 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.

Conclusions
In this paper, a fast tracking algorithm based on subspace updating is proposed for moving 2D CD sources under URA. e 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. e 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.