Two-Dimensional DOA Estimation for Coherently Distributed Sources with Symmetric Properties in Crossed Arrays

In this paper, a novel algorithm is proposed for the two-dimensional (2D) central direction-of-arrival (DOA) estimation of coherently distributed (CD) sources. Specifically, we focus on a centro-symmetric crossed array consisting of two uniform linear arrays (ULAs). Unlike the conventional low-complexity methods using the one-order Taylor series approximation to obtain the approximate rotational invariance relation, we first prove the symmetric property of angular signal distributed weight vectors of the CD source for an arbitrary centrosymmetric array, and then use this property to establish two generalized rotational invariance relations inside the array manifolds in the two ULAs. Making use of such relations, the central elevation and azimuth DOAs are obtained by employing a polynomial-root-based search-free approach, respectively. Finally, simple parameter matching is accomplished by searching for the minimums of the cost function of the estimated 2D angular parameters. When compared with the existing low-complexity methods, the proposed algorithm can greatly improve estimation accuracy without significant increment in computation complexity. Moreover, it performs independently of the deterministic angular distributed function. Simulation results are presented to illustrate the performance of the proposed algorithm.


Introduction
In recent decades, the problem of direction-of-arrival (DOA) estimation, which plays an important role in radar, sonar and wireless communication systems, has attracted a lot of attention. The most commonly considered system model in the DOA finding techniques is the point source model, where the signals are assumed to arrive at the array via a single path [1][2][3][4]. When dealing with a point source, conventional subspace-based algorithms, such as the multiple signal classification (MUSIC) algorithm [5,6] and the estimation of signal parameters via rotational invariance techniques (ESPRIT) algorithm [7,8], have high DOA estimation resolution. However, in many practical applications, the signals will reach the array through many rays reflected or scattered from the vicinity, which causes angular spreading. In these cases, directly applying the MUSIC and ESPRIT algorithms may lead to biased estimations. Therefore, some researchers have considered a more realistic signal model called the spatially distributed source model [9,10]. Depending on the correlation among different rays, distributed sources are classified into two types: coherently and incoherently distributed (CD and ID) sources. In this paper, we only consider the DOA estimation of the CD sources.
Many DOA estimation techniques for CD sources have been published. Since conventional subspace-based methods cannot be directly applied to a distributed source, some researchers sub-ULAs. Resorting to such relations, the central elevation and azimuth DOAs are estimated based on polynomial-root-based search-free method, respectively. Then simple parameter matching is accomplished by searching the minimums of the cost function of the estimated 2D angular parameters. The proposed algorithm does not require that the angular distribution functions of the multiple distributed sources are the same and known. In addition, it does not suffer additional errors induced by Taylor series approximation and high computational complexity brought about by spectrum-peak searching.
The rest of this paper is organized as follows: Section 2 presents the data model. In Section 3, we describe the proposed algorithm in detail. Some simulation results which illustrate the validity and performance of the proposed method are given in Section 4. Section 5 concludes the paper.
The following notations will be used throughout this paper. Superscript (·) * , (·) T , and (·) H represent the conjugate, transpose and conjugate transpose operations, respectively. The symbol ⊗ denotes the Schur-Hadamard product; E[·] stands for the mathematical expectation and det(·) is the matrix determinant; diag[·] is a diagonal matrix and the values in the brackets are the diagonal elements.

Signal Model
Let us consider the plane crossed array presented in Figure 1. The array is centered at the origin of the three-dimensional coordinate system with two ULAs directed along the y-axis and z-axis. The ULAs Y a and Z a are composed of M y and M z omni-directional antenna elements, respectively. The distance between adjacent sensors is d in the two ULAs. We assume that there are D narrowband CD sources impinging on the crossed array. The observation vectors of Y a and Z a at time t are given by [17][18][19][20][21].
Sensors 2017, 17, 1300 3 of 14 establish the generalized rotational invariance relations inside the GAMs for the two sub-ULAs. Resorting to such relations, the central elevation and azimuth DOAs are estimated based on polynomial-root-based search-free method, respectively. Then simple parameter matching is accomplished by searching the minimums of the cost function of the estimated 2D angular parameters. The proposed algorithm does not require that the angular distribution functions of the multiple distributed sources are the same and known. In addition, it does not suffer additional errors induced by Taylor series approximation and high computational complexity brought about by spectrum-peak searching. The rest of this paper is organized as follows: Section 2 presents the data model. In Section 3, we describe the proposed algorithm in detail. Some simulation results which illustrate the validity and performance of the proposed method are given in Section 4. Section 5 concludes the paper.
The following notations will be used throughout this paper. Superscript is a diagonal matrix and the values in the brackets are the diagonal elements.

Signal Model
Let us consider the plane crossed array presented in Figure 1. The array is centered at the origin of the three-dimensional coordinate system with two ULAs directed along the y-axis and z-axis. The The distance between adjacent sensors is d in the two ULAs. We assume that there are D narrowband CD sources impinging on the crossed array. The observation vectors of a Y and a Z at time t are given by [17][18][19][20][21].
where y(t) is the M y × 1 array output vector of the sub-array Y a ; z(t) is the M z × 1 array output vector of the sub-array Z a ; s i (θ, γ, t; µ i ) is the complex random angular signal density function of the i-th source. The vector µ i = (θ i , σ θ i , γ i , σ γ i ) determines the central azimuth DOA θ i , the azimuth angular spread σ θ i , the central elevation DOA γ i and the elevation angular spread σ γ i of the i-th sensor; n y (t) and n z (t) are Gaussian white noise with zero-mean and variance, while σ 2 n ; a y (θ, γ) and a z (θ, γ) are the two array manifold vectors at direction (θ, γ): where η = 2πd/λ, and λ is the wavelength of the coming signal. For a 2D CD source, the angular signal density function s i (θ, γ, t; µ i ) can be written as: where s i (t) is a random variable and ρ i (θ, γ; µ i ) is the deterministic angular distribution function. Define the GAM vectors of distributed source for subarray Y a and Z a as follows: For small angular extensions, we have the following closed forms [11,20]: where g y (µ i ) and g z (µ i ) are the ASDW vectors. The observation vectors in (1) and (2) can be written as: where s(t) = [s 1 (t), s 2 (t), . . . , s D (t)] T is a D × 1 signal vector, and B y (µ) and B z (µ) are the GAM matrices, which are composed of D GAM vectors: The total array output vector is expressed as:

The Proposed Algorithm
This section consists of four parts. Firstly, the symmetric property of an ASDW vector is identified in a centro-symmetric array. Then, by making use of the symmetric property of the ASDW vectors in the two sub-ULAs Z a and Y a , we establish two generalized rotational invariance relations for the GAM vectors. On the premise of such relations, the central elevation and azimuth DOAs are obtained by using a polynomial-root-based search-free approach, respectively. Afterwards, and when multiple CD sources exist, a simple parameter matching approach is addressed. Finally, we provide the algorithm's realization steps and an analysis of the computational complexity.

Symmetric Property of an ASDW Vector in a Centro-Symmetric Array
In this part, the symmetric property of an ASDW vector in a centro-symmetric array is derived in detail. Let us consider a centro-symmetric array consisting of M identical antenna elements centered at the coordinate origin, where the m-th sensor is placed at (x m , y m , z m ) for m = 1, 2 . . . , M. The array manifold vector in direction (θ, γ) is expressed as: If we define θ = θ i + θ and γ = γ i + γ, in which θ i , γ i are the central azimuth DOA and the central elevation DOA of the i-th source, and θ, γ are the corresponding random angular deviations , the GAM vector can be presented as: For small angular extensions, the m-th element of b(µ i ) can be written as [11]: Thus, the m-th element of ASDW vector is given by: where: We have (x m , y m , z m ) = −(x m+M/2 , y m+M/2 , z m+M/2 ) in the centrosymmetric array, thus ς m = −ς m+M/2 . Respecting the fact that ρ( θ, γ; µ i ) is an even function (see Appendix A), we can obtain the symmetric property of the ASDW vector such as: It is obvious that Y a and Z a are centrosymmetric arrays, thus we have:

Central Elevation DOA Estimation
For sub-array Z a , and owing to the symmetric property of the ASDW vector in (16), we can establish the following generalized rotational invariance relation of the GAM vector: where Π M Z is the M z × M z exchange matrix with ones on its anti-diagonal and zeros elsewhere. Ψ(γ i ) is an M z × M z diagonal matrix which is given by: If we define the complex variable k = e jη cos γ i , Ψ(γ i ) can be written as: According to the observation vector z(t) in (7), the covariance matrix of z(t) is expressed as: where R s = E s(t)s H (t) is the signal covariance matrix of the CD sources. The eigenvalue decomposition of R z is given by: where U z is the signal subspace matrix, whose columns are the eigenvectors corresponding to the D largest eigenvalues of R z . When R s is of full rank, the subspace spanned by the columns of U z is equal to the subspace spanned by the columns of B z (µ). At this point, there exists a unique non-singular According to the generalized rotational invariance relation in (17), we can formulate a matrix F z (k): Therefore, when Ψ(k i ) = Ψ(k) for i = 1, 2, . . . , D, the i-th column of F z (k) is a zero vector. Hence, F z (k) is rank deficient and the determinant of F H z (k)F z (k) is zero. The central elevation DOA estimationsγ i (i = 1, 2 . . . , D) can be obtained by finding the highest D peaks of the following spectrum function: However, estimator (23) involves computationally intensive spectral-peak searching. In order to reduce the complexity, we derive a polynomial-root-based search-free approach. The denominator of (23) can be written as the following polynomial: It is obvious that the central elevation DOAs can be obtained by rooting this polynomial. Note that the roots of (24) appear in conjugate reciprocal pairs, as in the conventional root-MUSIC [25]. To find the D central elevation DOAs, we select the D roots k i (i = 1, 2 . . . , D) inside the unit circle that maximize (23). Finally, the central elevation DOA estimates are obtained by: Sensors 2017, 17, 1300 7 of 13

Central Azimuth DOA Estimation
The method is similar to that of the central elevation DOA estimate, thus we simplify the process of deduction. For the centrosymmetric sub-array Y a , we also have the generalized shift invariance relation: where Ω(θ i , ϕ i ) is a M y × M y diagonal matrix which is given by: If we define the complex variable l = e jη sin θ i sin γ i , Ω(θ i , ϕ i ) can be written as: Let U y be the signal subspace matrix, whose columns are the eigenvectors corresponding to the D largest eigenvalues of R y = E y(t)y H (t) . Similarly, there exists a unique non-singular D × D matrix T 2 such that U y = B y (µ)T 2 . Thus, we introduce a matrix F y (l) which is expressed as: The central azimuth DOA estimationsθ i (i = 1, 2 . . . , D) can be obtained by rooting the following polynomial: Similarly, we use the roots inside the unit circle, and select the D roots l i (i = 1, 2 . . . , D) that maximize the spectral function such as: The values of sinθ i · sinγ i for i = 1, 2, . . . , D are obtained as:

The Parameter Matching Method
For only one CD source, the central elevation and azimuth DOAs can be estimated directly using (25) and (31). However, when multiple CD sources exist, the estimated elevation and azimuth DOAs are required to be matched. To perform the pair-matching procedure, we need to consider the GAM vector of the whole cross array such as: Let J be an (M y + M z ) × (M y + M z ) selection matrix which is defined as: Based on the symmetric property of the ASDW vector in (16), we have the following generalized rotational invariance relation: where Φ(θ i , γ i ) is an (M y + M z ) × (M y + M z ) diagonal matrix given by: Let U x be the signal subspace matrix of R x = E x(t)x H (t) . Similarly, there exists a unique non-singular D × D matrix T 3 such that U x = B x (µ)T 3 . We can also introduce a matrix F x (θ, γ) which is written as: Similarly, when Φ(θ i , γ i ) = Φ(θ, γ) for i = 1, 2, . . . , D, the i-th column of F x (θ, γ) is a zero vector. Therefore, the central elevation and azimuth DOA estimations can be matched by minimizing of the following cost function: If we pickγ i from the elevation DOA estimations {γ 1 ,γ 2 , . . .γ D }, there will be D pairs of 2D central DOAs forγ i , which is given by (θ i1,γi ), (θ i2,γi ), . . . (θ iD,γi ) . We then substitute the DOA estimations into (37) and calculate the function value f (θ, γ). If f (θ ij ,γ i ) for j = 1, 2, . . . D is the largest, then (θ ij ,γ i ) is the correct match.

Algorithm Implementation and Complexity Analysis
Now, we can summarize the proposed algorithm as follows: Step 1: Calculate the covariance matrix R z = E z(t)z H (t) . Through the eigen-decomposition of R z , obtain the signal subspace matrix U z .
Step 2: Construct the matrix F z (k) in (22), and root the polynomial in (24) to obtain the central elevation DOA estimationsγ i for i = 1, 2, . . . D. It is noted that the roots are inside a unit circle and maximize (23).
Step 3: Calculate the covariance matrix R y = E y(t)y H (t) . Through the eigen-decomposition of R y , obtain the signal subspace matrix U y .
Step 4: Construct the matrix F y (k) in (28), and root the polynomial in (29) to obtain sinθ i · sinγ i for i = 1, 2, . . . D. It is noted that the roots are inside a unit circle and maximize (30).
Step 5: Compute all the possible 2D DOAs (θ ij ,γ i ) for the elevation DOA estimationsγ i . Calculate the function values f (θ ij ,γ i ) for j = 1, 2, . . . D in (37). The largest one is the correct match.
Step 6: Repeat the process in Step 5 to match all the parameters.
Next, when the number of sensor elements M and the number of snapshots L change, we analyze the computational complexity of the proposed algorithm in comparison with the SOS algorithm in [18], the CC algorithm in [21] and Zheng's algorithm in [20]. The main computational cost of the proposed algorithm is mostly made of four operations: the estimation of the covariance matrix, the eigen-decomposition of the covariance matrix, the polynomials rooting, and the pair-matching procedure. Specifically, the cost involved by the estimation of covariance matrices R x , R y and R z is found to be in O(6M 2 L). The eigen-decomposition of the covariance matrices R x , R y and R z needs O(10M 3 ) multiplications. The computational complexity of rooting polynomials h z (k) and h y (k) is found to be in O(2M 3 ), and the pair-matching procedure adds O(D 5 + D 4 M) multiplications to the proposed algorithm. In above, the computational complexity of the proposed algorithm is O(12M 3 + 6M 2 L + D 5 + D 4 M). Moreover, the main computational complexity of the SOS algorithm, the CC algorithm and Zheng's algorithm is given in Table 1 (N denotes the number of searching points). Table 1. Comparison of different algorithms in computational complexity.

Algorithm Main Computational Complexity
Proposed algorithm When the number of searching points N is large, it is clear to see that the propose algorithm provides lower computational cost compared to the SOS algorithm. Although the computational complexity of the proposed algorithm is higher than Zheng's algorithm and the CC algorithm, it is not significant increment since the proposed algorithm does not require any spectrum searching. In addition, and unlike the SOS algorithm, Zheng's algorithm and the CC algorithm, the proposed algorithm does not use the Taylor series approximation to establish the rotational invariance relation, as this approximation may introduce additional errors.

Simulation Results and Performance Analysis
In the following experiments, noise is a complex Gaussian process with zero mean. The number of snapshots is 200. We use the root mean squared error (RMSE) to evaluate the estimation performance, where the RMSE of the central azimuth and elevation DOAs (RMSE(θ) and RMSE(γ)) are defined as: whereθ i and θ i are the estimated and true central azimuth DOA of the i-th source, respectively. Additionally,γ i and γ i are the estimated and true central elevation DOA of the i-th source, respectively. In the following simulations, the signal power of sources is assumed to be the same. In addition, signal-to-noise ratio (SNR) is defined as: where σ 2 s is the signal power of sources, while σ 2 n is the variance of noise.

Effect of Different Deterministic Angular Distributed Functions
In this part, we examine if the proposed algorithm works properly for different angular distributed functions. The sub-arrays Y a and Z a are both composed of M y = M z = 16 sensors. The distance between adjacent sensors is 0.5λ. The parameters of two CD sources are µ 1 = (20 • , 3 • , 20 • , 5 • ) and . Their deterministic angular distributed functions are Gaussian and uniform shaped, respectively. The SNR is 15 dB. For 30 independent trials, the central DOA estimations of CD sources are plotted in Figure 2. It can be seen that the proposed algorithm can give the correct DOA estimations for cases where different CD sources have different deterministic angular distributed functions, or unknown deterministic angular distributed functions.

Performance Comparison
In this part, we compare the estimation accuracy of the proposed algorithm with the SOS algorithm in [18], the CC algorithm in [21] and Zheng's algorithm in [20] with respect to SNR from 0 dB to 30 dB. The Cramér-Rao lower bound (CRLB) is also used as a benchmark [26]. The sub-arrays a Y and a Z of the crossed array in proposed algorithm are both composed of 16   y z M M sensors. The arrays in the SOS algorithm and the CC algorithm are both composed of 32 sensors. Since the number of antenna elements in Zheng's algorithm must be odd, we set it to 33. In these algorithms, the distance between adjacent sensors in a sub-array is 0 5λ . , the vertical distance between the two sub-arrays is 0 5λ . , and the radius of the UCA is 4 16 λ π / ( sin( / )) . The parameters of two Gaussianshaped CD sources are 1 Figure 3. It is clearly indicated that the estimation accuracy of the proposed algorithm is higher

Performance Comparison
In this part, we compare the estimation accuracy of the proposed algorithm with the SOS algorithm in [18], the CC algorithm in [21] and Zheng's algorithm in [20] with respect to SNR from 0 dB to 30 dB. The Cramér-Rao lower bound (CRLB) is also used as a benchmark [26]. The sub-arrays Y a and Z a of the crossed array in proposed algorithm are both composed of M y = M z = 16 sensors. The arrays in the SOS algorithm and the CC algorithm are both composed of 32 sensors. Since the number of antenna elements in Zheng's algorithm must be odd, we set it to 33. In these algorithms, the distance between adjacent sensors in a sub-array is 0.5λ, the vertical distance between the two sub-arrays is 0.5λ, and the radius of the UCA is λ/(4 sin(π/16)). The parameters of two Gaussian-shaped CD sources are µ 1 = (20 • , 2 • , 60 • , 3 • ) and µ 2 = (15 • , 3 • , 70 • , 2 • ), respectively. Based on 500 Monte Carlo experiments, the RMSE curves of the central DOA estimations versus SNR are shown in Figure 3. It is clearly indicated that the estimation accuracy of the proposed algorithm is higher than the SOS algorithm, the CC algorithm and Zheng's algorithm, which arises from the fact that the proposed algorithm does not suffer additional errors brought about by Taylor series approximation.

Effect of Snapshots
In this part, we illustrate the influence of the number of snapshots on the performance of the proposed algorithm. The number of snapshots varies from 100 to 900. The SNR is fixed to 15 dB and the other parameters are the same as in Section 4.2. Based on 500 Monte Carlo experiments, the RMSE curves for different algorithms are shown in Figure 4, from which we can draw similar conclusion as in Section 4.2.

Effect of Snapshots
In this part, we illustrate the influence of the number of snapshots on the performance of the proposed algorithm. The number of snapshots varies from 100 to 900. The SNR is fixed to 15 dB and the other parameters are the same as in Section 4.2. Based on 500 Monte Carlo experiments, the RMSE curves for different algorithms are shown in Figure 4, from which we can draw similar conclusion as in Section 4.2.

Effect of Snapshots
In this part, we illustrate the influence of the number of snapshots on the performance of the proposed algorithm. The number of snapshots varies from 100 to 900. The SNR is fixed to 15 dB and the other parameters are the same as in Section 4.2. Based on 500 Monte Carlo experiments, the RMSE curves for different algorithms are shown in Figure 4, from which we can draw similar conclusion as in Section 4.2.

Effect of the Central Elevation and Azimuth DOAs
In the last example, we examine the performance of the proposed method versus the central elevation and azimuth DOAs for a Gaussian-shaped CD source with First, let us consider in Figure 5a  ), and our method has still satisfying estimation performance for approaching the lower bound.

Effect of the Central Elevation and Azimuth DOAs
In the last example, we examine the performance of the proposed method versus the central elevation and azimuth DOAs for a Gaussian-shaped CD source with σ θ i = σ γ i = 1 • . First, let us consider in Figure 5a the influence of the central azimuth DOA on performance. Assume that γ i = 30 • , SNR = 10 dB and the number of snapshots L = 200. As can be seen from the Figure 5a, the RMSE ofθ i estimated by the proposed method increases dramatically when the central azimuth DOA approaches the lower bound (θ i = −90 • ) or the upper bound (θ i = 90 • ), but our method can still estimate effectively the central azimuth DOA. Next, the influence of the central elevation DOA on performance is examined in Figure 5b. At this time, we assume that θ i = 30 • , SNR = 10 dB and L = 200. Similarly, the RMSE ofγ i estimated by the proposed method also increases dramatically when the central elevation DOA approaches the lower bound (γ i = 0 • ), and our method has still satisfying estimation performance for approaching the lower bound.

Conclusions
In this paper, we have presented a new approach for estimating the 2D central DOA of CD sources using a centrosymmetric crossed array. Instead of using the Taylor series approximation, we derive the symmetric property of an ASDW vector in a centrosymmetric array, based on which the generalized shift invariance relations inside the GAMs are established in the two sub-ULAs. Resorting to such relations, the central elevation and azimuth DOAs are estimated based on the

Conclusions
In this paper, we have presented a new approach for estimating the 2D central DOA of CD sources using a centrosymmetric crossed array. Instead of using the Taylor series approximation, we derive the symmetric property of an ASDW vector in a centrosymmetric array, based on which the generalized shift invariance relations inside the GAMs are established in the two sub-ULAs. Resorting to such relations, the central elevation and azimuth DOAs are estimated based on the polynomial-root-based search-free method, respectively. After that, the pair-matching method is presented. The proposed algorithm performs independently of the deterministic angular distributed function. Compared to the existing low-complexity algorithms, the proposed algorithm does not suffer additional errors brought by the Taylor series approximation, which allows it to achieve a higher estimation accuracy. Moreover, the proposed algorithm does not suffer high computational complexity brought by spectrum-peak searching.