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.

1. 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 studied modifications to the MUSIC algorithm, which gave rise to the distributed signal parameter estimator [11], dispersed parametric estimator [12] and Vec-MUSIC estimator [13]. These three methods were established based on the assumption that the distribution shapes of multiply distributed sources are identical and known. In addition, the computational complexity is high because of the need for a two-dimensional (2D) search. The literature [14] used an unstructured model for the part of covariance matrix, where the 2D search problem was replaced by a successive one-dimensional (1D) search. With even lower computational complexity, a search-free algorithm called the spread root-MUSIC algorithm [15] was proposed to fit a two-ray model of the data. However, low-complexity was obtained in the special case where only one distributed source existed. In [16], the authors considered two identical and closely spaced sub-arrays. When the distance between the two sub-arrays was far shorter than the wavelength, an approximate rotational invariance relation between the two sub-arrays was obtained based on Taylor series approximation, and finally the central DOA of the CD source can be estimated by total least square estimation parameter via rotational invariance techniques (TLS-ESPRIT) using the generalized array manifold (GAM) model. All of the aforementioned works [11,12,13,14,15,16] were designed for 1D DOA and angular spread estimation of distributed sources. However, when the distributed source and the receiving sensor array are not in the same plane, it is reasonable to instead model the source as a 2D distributed source.

Since a 2D CD source is characterized by four parameters: the central azimuth direction, the azimuth angular extension, the central elevation direction and the elevation angular extension, the conventional optimum estimators will be computationally expensive owing to high-dimensional parameters [17]. Consequently, it is very necessary for 2D distributed sources to find some suboptimum algorithms to reduce the computational cost. To date, several low-complexity DOA estimation algorithms for 2D CD sources have been proposed. Specifically, the authors in [18] considered a pair of uniform circular arrays (UCAs). Preliminary estimations of central elevation DOAs were obtained using TLS-ESPRIT. Next, by using the estimated elevation DOAs, a sequential one-dimensional searching (SOS) method was proposed to estimate the central azimuth DOAs. In [19], using two parallel uniform linear arrays (ULAs), a low-complexity algorithm without searching was proposed for CD sources. Similarly, central elevation DOAs are obtained based on the approximate rotational invariance relation between the two ULAs. Instead of SOS, the quadric rotational invariance property (QRIP) of the GAM was used to estimate the central azimuth DOAs. Finally, a parameter matching approach was given to obtain the correct DOA estimation. In [20], the central elevation and azimuth DOAs were both estimated based on TLS-ESPRIT, which used two parallel ULAs, and the parameter matching method was also required. In order to avoid the parameter matching procedure, the literature [21] estimated the central elevation and central azimuth DOAs by applying the singular value decomposition method to the cross-correlation (CC) matrix of the received data in the double parallel ULAs. However, all the algorithms in [18,19,20,21] were all based on the special array geometry composed of two sub-arrays. The approximate rotational invariance relation between the two sub-arrays was obtained by using the one-order Taylor series approximation, which may introduce additional errors and affect the estimation accuracy.

In array processing, the crossed array is a commonly used 2D array geometry [22]. Compared to other 2D arrays such as the UCA and plane rectangular array, the crossed array can provide a larger aperture and hence offer better resolution for a given number of elements. Moreover, the crossed array consists of two intersecting ULAs working independently, thus the computational complexity is only double that of a single dimensional array. Several 2D estimation algorithms based on the crossed array have been proposed [23,24]. However, they are all based on the point source model. To the best of our knowledge, there have been few reports about the DOA estimation for CD sources in a crossed array.

In this paper, we consider a crossed array and divide it into two sub-ULAs. In particular, instead of using the Taylor series approximation, we prove the symmetric property of the angular signal distributed weight (ASDW) vector for an arbitrary centrosymmetric array, and use this property to 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 , , and represent the conjugate, transpose and conjugate transpose operations, respectively. The symbol denotes the Schur-Hadamard product; stands for the mathematical expectation and is the matrix determinant; is a diagonal matrix and the values in the brackets are the diagonal elements.

2. 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 and are composed of and omni-directional antenna elements, respectively. The distance between adjacent sensors is in the two ULAs. We assume that there are narrowband CD sources impinging on the crossed array. The observation vectors of and at time are given by [17,18,19,20,21]. where is the array output vector of the sub-array ; is the array output vector of the sub-array ; is the complex random angular signal density function of the i-th source. The vector determines the central azimuth DOA , the azimuth angular spread , the central elevation DOA and the elevation angular spread of the i-th sensor; and are Gaussian white noise with zero-mean and variance, while ; and are the two array manifold vectors at direction : where , and is the wavelength of the coming signal.

Figure 1

Geometry of the considered crossed array.

For a 2D CD source, the angular signal density function can be written as: where is a random variable and is the deterministic angular distribution function.

Define the GAM vectors of distributed source for subarray and as follows:

For small angular extensions, we have the following closed forms [11,20]: where and are the ASDW vectors. The observation vectors in (1) and (2) can be written as: where is a signal vector, and and are the GAM matrices, which are composed of D GAM vectors:

The total array output vector is expressed as:

3. 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 and , 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.

3.1. 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 for . The array manifold vector in direction is expressed as:

If we define and , in which , 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 can be written as [11]:

Thus, the m-th element of ASDW vector is given by: where:

We have in the centrosymmetric array, thus . Respecting the fact that is an even function (see Appendix A), we can obtain the symmetric property of the ASDW vector such as:

It is obvious that and are centrosymmetric arrays, thus we have:

3.2. Derivation

3.2.1. Central Elevation DOA Estimation

For sub-array , 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 is the exchange matrix with ones on its anti-diagonal and zeros elsewhere. is an diagonal matrix which is given by:

If we define the complex variable , can be written as:

According to the observation vector in (7), the covariance matrix of is expressed as: where is the signal covariance matrix of the CD sources. The eigenvalue decomposition of is given by: where is the signal subspace matrix, whose columns are the eigenvectors corresponding to the largest eigenvalues of . When is of full rank, the subspace spanned by the columns of is equal to the subspace spanned by the columns of . At this point, there exists a unique non-singular matrix such that . According to the generalized rotational invariance relation in (17), we can formulate a matrix :

Therefore, when for , the i-th column of is a zero vector. Hence, is rank deficient and the determinant of is zero. The central elevation DOA estimations can be obtained by finding the highest 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 central elevation DOAs, we select the roots inside the unit circle that maximize (23). Finally, the central elevation DOA estimates are obtained by:

3.2.2. 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 , we also have the generalized shift invariance relation: where is a diagonal matrix which is given by:

If we define the complex variable , can be written as: .

Let be the signal subspace matrix, whose columns are the eigenvectors corresponding to the largest eigenvalues of . Similarly, there exists a unique non-singular matrix such that . Thus, we introduce a matrix which is expressed as:

The central azimuth DOA estimations can be obtained by rooting the following polynomial:

Similarly, we use the roots inside the unit circle, and select the roots that maximize the spectral function such as:

The values of for are obtained as:

3.2.3. 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 be an 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 is an diagonal matrix given by:

Let be the signal subspace matrix of . Similarly, there exists a unique non-singular matrix such that . We can also introduce a matrix which is written as:

Similarly, when for , the i-th column of 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 from the elevation DOA estimations , there will be pairs of 2D central DOAs for , which is given by . We then substitute the DOA estimations into (37) and calculate the function value . If for is the largest, then is the correct match.

3.2.4. Algorithm Implementation and Complexity Analysis

Now, we can summarize the proposed algorithm as follows:

Calculate the covariance matrix . Through the eigen-decomposition of , obtain the signal subspace matrix .

Construct the matrix in (22), and root the polynomial in (24) to obtain the central elevation DOA estimations for . It is noted that the roots are inside a unit circle and maximize (23).

Calculate the covariance matrix . Through the eigen-decomposition of , obtain the signal subspace matrix .

Construct the matrix in (28), and root the polynomial in (29) to obtain for . It is noted that the roots are inside a unit circle and maximize (30).

Compute all the possible 2D DOAs for the elevation DOA estimations . Calculate the function values for in (37). The largest one is the correct match.

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 , and is found to be in . The eigen-decomposition of the covariance matrices , and needs multiplications. The computational complexity of rooting polynomials and is found to be in , and the pair-matching procedure adds multiplications to the proposed algorithm. In above, the computational complexity of the proposed algorithm is . 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	O(12M^3+6M^2L+D^5+D^4M)
SOS algorithm	O(8M^3+4M^2L+N(D^3+3D^2M))
CC algorithm	O(M^3+M^2L+D^3)
Zheng’s algorithm	O((2M+1)^3+(2M+1)^2L+D^3)

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.

4. 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 ( and ) are defined as: where and are the estimated and true central azimuth DOA of the i-th source, respectively. Additionally, and 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 is the signal power of sources, while is the variance of noise.

4.1. 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 and are both composed of sensors. The distance between adjacent sensors is . The parameters of two CD sources are 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.

Figure 2

The 2D central DOA estimation results of the proposed algorithm (30 trials).

4.2. 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 and of the crossed array in proposed algorithm are both composed of 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 , the vertical distance between the two sub-arrays is , and the radius of the UCA is . The parameters of two Gaussian-shaped CD sources are and , 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.

Figure 3

(a) RMSE curves of the central azimuth DOA estimations versus SNR; (b) RMSE curves of the central elevation DOA estimations versus SNR.

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

Figure 4

(a) RMSE curves of the central azimuth DOA estimations versus the number of snapshots; (b) RMSE curves of the central elevation DOA estimations versus the number of snapshots.

4.4. 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 the influence of the central azimuth DOA on performance. Assume that , SNR = 10 dB and the number of snapshots . As can be seen from the Figure 5a, the RMSE of estimated by the proposed method increases dramatically when the central azimuth DOA approaches the lower bound () or the upper bound (), 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 , SNR = 10 dB and . Similarly, the RMSE of estimated by the proposed method also increases dramatically when the central elevation DOA approaches the lower bound (), and our method has still satisfying estimation performance for approaching the lower bound.

Figure 5

(a) RMSE curves versus the central azimuth DOA; (b) RMSE curves versus the central elevation DOA.

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