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.

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) [1113], 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 [1930]. 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, , and mean inverse, transpose, and Hermitian transpose of a matrix. , (./), 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 Figure 1. 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 aswhere nmk(t) is the noise received. As for CD sources, the angular signal density function can be written aswhere si (t) reflects time characteristics of the ith CD sources, (θ, φ; ui) is the deterministic angular signal density function (ASDF) of the ith source, and is the power of the source. Define signal vector s (t) as follows:

Assumed different CD sources are incoherent, we have

The generalized manifold coefficient of the (m, k)th sensor for the ith 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

The 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):

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 aswhere A1 is the generalized manifold matrix of subarray X1, which has the following expression:

Define subarray Xk 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 iswhere the noise vector can be expressed as

Take the first M − 1 sensors of subarray Xk to form subarray and the last M − 1 sensors to form subarray , and the generalized manifold matrix of subarray can be expressed as

Generalized manifold matrix of subarray 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:

Receive vectors of subarrays and can be written as

Receive vectors of URA can be written aswhere A is the generalized manifold matrix of URA, which have the following expression:

Receive noise vector of URA can be written aswhere

The covariance matrix of X(t) can be written as

The signal subspace W composed of eigenvectors corresponding to the largest q 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:where

Define as a vector selecting the first M-1 rows of W and as the 2nd to Mth rows of W. is a vector selecting [( –1)2M + 1]th to [( –1)2M + M − 1]th rows of W, 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), , , and have the following relations:

According to the rotating invariant relations described by equations (17) and (18), we have

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 of and eigenvalues of have the same eigenvector.

Eigenvalue can be obtained by eigendecomposition . Then, nominal elevation can be obtained aswhere is the ith eigenvalue of , and superscript of denotes the eigenvalue calculated by and . The corresponding eigenvector of is . Then, the eigenvalue of , which has the same position as , can be expressed aswhere 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:

It can be concluded that the nominal elevation and the nominal azimuth are the mean values obtained from the K/2 group of , , and .

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 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. Then, 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

Authors of [16] have proved that expression (38) converges to the signal subspace of the covariance matrix.

Define compressed signal vector as

Then, equation (37) can be expressed aswhere

Rw(t) can be regarded as the cross covariance matrix of the received signal vector X(t) and the compressed signal vector V(t). is the square root matrix of Rw(t), which have the following expression:

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.

The 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 , , and (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 RZ, which is , and calculate a pseudoinverse of , which is .

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 Figure 1 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. 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. 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.

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

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

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

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.

Appendix

Define as the small deviations of θi and as the small deviations of φi. Replace θ with () and φ with (); we can obtain the approximate representations of cos θ sin φ and cos φ utilizing the first term of the Taylor series expansions. Consider both sides of equation (8).

As , the following equation can be obtained.

Consider the right side of equation (9).

Noticing , we obtain the relation between and .

Similarly, we have

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors would like to thank the editorial board and reviewers for the improvement of this paper. This research was funded by the National Natural Science Foundation of China (grant nos. 51776222 and 61903378) and Field Foundation of China (grant no. 61400020109).