Two-Dimensional Direction of Arrival (DOA) Estimation for Rectangular Array via Compressive Sensing Trilinear Model

. We investigate the topic of two-dimensional direction of arrival (2D-DOA) estimation for rectangular array. This paper links angle estimation problem to compressive sensing trilinear model and derives a compressive sensing trilinear model-based angle estimation algorithm which can obtain the paired 2D-DOA estimation. The proposed algorithm not only requires no spectral peak searchingbutalsohasbetterangleestimationperformancethanestimationofsignalparametersviarotationalinvariancetechniques (ESPRIT)algorithm.Furthermore,theproposedalgorithmhascloseangleestimationperformancetotrilineardecomposition.The proposedalgorithmcanberegardedasacombinationoftrilinearmodelandcompressivesensingtheory,anditbringsmuchlower computationalcomplexityandmuchsmallerdemandforstoragecapacity.Numericalsimulationspresenttheeffectivenessofour approach.


Introduction
Array signal processing has received a significant amount of attention during the last decades due to its wide application in radar, sonar, radio astronomy, and satellite communication [1].The direction of arrival (DOA) estimation of signals impinging on an antenna array is a fundamental problem in array signal processing, and many DOA estimation methods [2][3][4][5][6][7] have been proposed for its solution.They contain estimation of signal parameters via rotational invariance techniques (ESPRIT) algorithm [2,3], multiple signal classification (MUSIC) algorithm [4], Root-MUSIC [5], matrix pencil methods [6], and so on.Compared with linear arrays, uniform rectangular array can identify twodimensional DOA (2D-DOA).2D-DOA estimation with rectangular array has received considerable attention in the field of array signal processing [7][8][9][10][11][12].ESPRIT algorithms in [8][9][10][11] have exploited the invariance property for 2D-DOA estimation in uniform rectangular array.Parallel factor analysis (PARAFAC) in [12], which is also called trilinear decomposition method, was proposed for 2D-DOA estimation for uniform rectangular array, and it has better angle estimation performance than ESPRIT.MUSIC algorithm, as a subspace method, has good angle estimation performance and matches irregular arrays.It has been proved that twodimensional MUSIC (2D-MUSIC) algorithm [13] can be used for 2D-DOA estimation.However, the requirement of twodimensional (2D) spectrum searching renders much higher computational complexity.
Compressive sensing [14,15] has attracted a lot of attention recently, and it has been applied to image processing, machine learning, channel estimation, radar imaging, and penalized regression [16].According to the theory of compressive sensing, a signal that is sparse in some domain can be recovered via fewer samples than required by the Nyquist sampling theorem.The DOAs of sources form a sparse vector in the potential signal space, and, therefore, compressive sensing can be applied to DOA estimation.

International Journal of Antennas and Propagation
The superresolution property and ability of resolving coherent sources can be achieved when we apply it to the source location [17].Lots of the DOA estimation methods with compressive sensing just use one snapshot and are very sensitive to the noise.For multiple snapshots, ℓ 1 -SVD method [16] employed ℓ 1 norm to enforce sparsity and singular value decomposition to reduce complexity and sensitivity to noise, and sparse recovery for weighted subspace fitting in [17] improved the ℓ 1 -SVD method via the weight to the subspace.
Compared to matrix decomposition, trilinear decomposition has a distinctive and attractive feature: it is often unique [18][19][20][21][22].In the signal processing field, trilinear decomposition can be regarded as a generalization of ESPRIT and joint approximate diagonalization [19][20][21][22].The compressive sensing trilinear model-based algorithm discussed in this paper can be regarded as a combination of trilinear model and compressive sensing theory, which brings much lower computational complexity and much smaller demand for storage capacity.
The framework of compressive sensing for sparse lowrank tensor is proposed in [23] and used for signal detection and multiple-input-multiple-output radar in [24,25].In this paper, the problem of 2D-DOA estimation for rectangular array is linked to compressive sensing trilinear model.Exploiting this link, we derive a compressive sensing trilinear model-based 2D-DOA estimation algorithm for rectangular array.Firstly, we compress the received data to get a compressed trilinear model and then obtain the estimates of compressed direction matrices through performing trilinear decomposition for the compressed model.Finally, we formulate a sparse recovery problem through the estimated compressed direction matrices and apply the orthogonal matching pursuit (OMP) [26] to resolve it for 2D-DOA estimation.Due to compression, the proposed method has much lower computational complexity than conventional trilinear decomposition method [12] and 2D-MUSIC algorithm and requires much smaller storage capacity.We illustrate that the proposed algorithm has better angle estimation performance than ESPRIT algorithm.Furthermore, our algorithm can obtain paired elevation angles and azimuth angles automatically.We also derive the Cramer-Rao bound (CRB) for 2D-DOA estimation in rectangular array.Numerical simulations present the effectiveness of our approach.
The remainder of this paper is structured as follows.Section 2 presents the data model, and Section 3 proposes the compressed sensing trilinear model-based algorithm for 2D-DOA estimation in rectangular array.In Section 4, the simulation results are presented to verify improvement of the proposed algorithm, while the conclusions are drawn in Section 5.
If A is a -by- matrix and B is an -by- matrix, then the Kronecker product A ⊗ B is the -by- block matrix: where   is the (, ) element of the matrix A.
If A is an -by- matrix and B is a -by- matrix, then the Khatri-Rao product A ⊙ B is the -by- block matrix: where a  and b  are the th column of the matrices A and B, respectively.
If A is an -by- matrix and B is an -by- matrix, then the Hadamard product A ⊕ B is where   and   are the (, ) element of the matrices A and B, respectively.

Data Model
A rectangular array consisted of  ×  elements is shown in Figure 1, where the distance between two adjacent elements is .We consider signals in the far field, in which case the signal sources are far away enough that the arriving waves are essentially planes over the array.We assume that the noise is independent of the sources.It is also assumed that there are  noncoherent or independent sources, and the number of sources is preknown.  and   are the elevation angle and the azimuth angle of the th source, respectively.We assume the sources impinge on the array with different DOAs.
The received signal of the first subarray in the rectangular array is x 1 () = A  s() + n 1 (), where A  = [a  ( 1 ,  1 ), a  ( 2 ,  2 ), . . ., a  (  ,   )] ∈ C × with a  (  ,   ) = [1,  2 sin   cos   / , . . .,  2(−1) sin   cos   / ]  and  is the wavelength.n 1 () is the received additive white Gaussian noise of the first subarray.s() ∈ C ×1 is the source vector.The received signal of the th subarray in the rectangular array is x  () = A  Φ −1 s() + n  (), where Φ = diag( 2 sin  1 sin  1 / , . . .,  2 sin   sin   / ) and n  () is the received additive white Gaussian noise of the th subarray.Therefore, the received signal of the rectangular array is [27] x The signal x() ∈ C ×1 in (4) can also be denoted by where According to the definition of Khatri-Rao product, the signal in ( 5) can be rewritten as where ⊗ denotes Kronecker product.We collect  samples and define X = [x(1), x(2), . . ., x()] ∈ C × , which can be expressed as ] where ), n(2), . . ., n()] is the received additive white Gaussian noise matrix.N  ∈ C × ( = 1, . . ., ) is the noise matrix.Thus, X  ∈ C × in ( 7) is denoted as Equation ( 7) can also be denoted with the trilinear model [18,28] where A  (, ) is the (, ) element of the matrix A  and similarly for the others. ,, is noise part.X  ∈ C × ( = 1, . . ., ) can be regarded as slicing the three-dimensional data in a series of slices, which is shown in Figure 2.There are two more matrix system rearrangements, in which we have Y  = S  (A  )A   + N   ,  = 1, . . ., , and Z  = A    (S)A   + N   ,  = 1, . . ., , where N   and N   are noise matrices.Then, we form the matrices of Y ∈ C × and where . . . ,

2D-DOA Estimation Based on Compressive Sensing Trilinear Model
We link the problem of 2D-DOA estimation for rectangular array to compressive sensing trilinear model and derive a compressive sensing trilinear model-based 2D-DOA estimation algorithm.Firstly, we compress the received data to get a compressed trilinear model and then obtain the estimates of compressed direction matrices through performing trilinear decomposition for the compressed model.Finally, we formulate the sparse recovery problem for 2D-DOA estimation.

Trilinear Model Compression.
We compress the threeway data Χ ∈ C ×× into a smaller three-way data Χ  ∈ C   ×  ×  , where   < ,   < , and   < .The trilinear model compression processing is shown in Figure 3.We define the compression matrices as U ∈ C ×  , V ∈ C ×  , and W ∈ C ×  , and the compression matrices U, V, and W can be generated randomly or obtained by Tucker3 decomposition [23,29].We can use the Tucker3 decomposition,

Compression
Compressed trilinear model Trilinear model where tensor is decomposed into the core tensor, to obtain the compression matrices.The compression matrices should satisfy the restricted isometry property.And random Gaussian, Bernoulli, and partial Fourier matrices satisfy the restricted isometry property with number of measurements nearly linear in the sparsity level [30,31].
Then compress X ∈ C × in ( 7) to a smaller one as According to the property of Khatri-Rao product [23], we know Define A   = U  A  , A   = V  A  , and S  = W  S. Equation ( 11) is also denoted as where N  = (V  ⊗ U  )NW  .X  can be denoted by trilinear model.With respect to (10) and (11), we form the matrices of Y  and Z  according to the compressed data where N  and N  are the noise part.The compressed trilinear model may degrade the angle estimation performance.By trilinear model compression, the proposed method has much lower computational complexity than conventional trilinear decomposition method and requires much smaller storage capacity.Conventional compressive sensing is to compress the matrix, while our algorithm compresses the three-dimensional tensor.

Trilinear Decomposition.
Trilinear alternating least square (TALS) algorithm is an iterative method for estimating the parameters of a trilinear decomposition [18,28].We concisely show the basic idea of TALS: (1) update one matrix each time via LS, which is conditioned on previously obtained estimates of the remaining matrices; (2) proceed to update the other matrices; (3) repeat until convergence of the LS cost function [21,22].TALS algorithm is discussed as follows.
According to (15), least squares (LS) fitting is min and LS update for the matrix where Â  and Â  are previously obtained estimates of A   and A   , respectively.According to (16), LS fitting is min and LS update for A   is where Â  and Ŝ stand for the previously obtained estimates of A   and S  .Similarly, according to (17), LS fitting is min where Z is the noisy compressed signal.LS update for A   is where Theorem 1 (see [22]).Considering X   = A     (A   )S  ,  = 1, . . .,   , where A   ∈ C   × , S  ∈ C   × , and A   ∈ C   × , if where  A is the -rank of the matrix A [18], then A   , S  , and A   are unique up to permutation and scaling of columns.
When   ≤ ,   ≤ , and   ≥ , the identifiable condition is max(  ,   ) ≤  ≤   +   − 2. Hence, the proposed algorithm is effective when  ≤   +  −2 and the maximum number of sources that can be identified is   +   − 2.
After the trilinear decomposition, we obtain the estimates of the loading matrices where Π is a permutation matrix, and Δ 1 , Δ 2 , Δ 3 note for the diagonal scaling matrices satisfying Δ 1 Δ 2 Δ 3 = I  .E 1 , E 2 , and E 3 are estimation error matrices.After the trilinear decomposition, the estimates of A   , A   , and S  can be obtained.Scale ambiguity and permutation ambiguity are inherent to the trilinear decomposition problem.However, the scale ambiguity can be resolved easily by means of normalization, while the existence of permutation ambiguity is not considered for angle estimation.

Angle Estimation via Sparse Recovery. Use â󸀠
and â  to denote the th column of estimates Â  and Â  , respectively.According to the compression matrices, we have where a  and a  are the th column of A  , A  , respectively.n  and n  are the corresponding noise, respectively.  and   are the scaling coefficients.Construct two Vandermonde matrices A  ∈ C × and A  ∈ C × ( ≫ ,  ≫ ) composed of steering vectors corresponding to each potential source location as its columns: where g is a sampling vector and its th elements is g() = −1 + 2/,  = 1, 2, . . ., .The matrices A  and A  can be regarded as the completed dictionaries.Then (27a)-(27b) can be expressed as where x  and y  are sparse.The estimates of x  and y  can be obtained via  0 -norm constraint: where ‖ ⋅ ‖ 0 denotes the  0 -norm.‖x  ‖ 0 = 1; that is to say, there is only one nonzero element in the vector x  , similar to ‖y  ‖ 0 = 1.We can use the OMP recovery method [26] to find the nonzero element in x  or y  .The OMP algorithm tries to recover the signal by finding the strongest component in the measurement signal, removing it from the signal, and searching the dictionary again for the strongest atom that is presented in the residual signal [32].We extract the index of the maximum modulus of elements in x  and y  , respectively, noted as   and   .According to the corresponding columns in A  and A  , we obtain g(  ) and g(  ), which are estimates of sin   cos   and sin   sin   .We define International Journal of Antennas and Propagation   = g(  ) + g(  ), and then the elevation angles and azimuth angles can be obtained via θ = sin −1 (abs (  )) ,  = 1, . . ., , (31a) where abs(⋅) is the modulus value symbol and angle(⋅) is to get the angle of an imaginary number.As the columns of the estimated matrices Â  and Â  are automatically paired, then the estimated elevation angles and azimuth angles can be paired automatically.

The Procedures of the Proposed Algorithm.
Till now, we have achieved the proposal for the compressive sensing trilinear model-based 2D-DOA estimation for rectangular array.We show major steps of the proposed algorithm as follows.
Step 1. Form the three-way matrix Χ ∈ C ×× , then compress the three-way matrix into a much smaller threeway matrix Step 2. Perform trilinear decomposition through TALS algorithm for the compressed three-way matrix to obtain the estimation of A   , A   , and S  .
Step 3. Estimate the sparse vectors.
Remark A. Because the trilinear decomposition brings the same permutation ambiguity for the estimates A   , A   , and S  , the estimated elevation angles and azimuth angles are paired automatically.
Remark B. The conventional compressive sensing method formulates an angle sampling grid for sparse recovery to estimate angles.When it is applied to 2D-DOA estimation, both elevation and azimuth angles must be sampled, and it results in a two-dimensional sampling problem which brings much heavier cost for sparse signal recovery.In this paper, sin   cos   (or sin   sin   ) is bundled into a single variable in the range of −1 to 1.The bundled variable is sampled for sparse recovery to obtain the estimates of sin   cos   and sin   sin   , respectively.Afterwards, the elevation and azimuth angles are estimated through the estimates of sin   cos   and sin   sin   .
Remark C. If the number of sources  is unknown, it can be estimated by performing singular value decomposition for received data matrix X in (7) and finding the number of largest singular values [33].We also use some lowercomplexity algorithm in [34] for estimating the number of the sources.
Remark D. When the coherent sources impinge on the array, we can use the parallel profiles with linear dependencies (PARALIND) model [35,36], which is a generalization of PARAFAC suitable for solving problems with linear dependent factors, to resolve coherent DOA estimation problem.

Complexity Analysis and CRB.
The proposed algorithm has much lower computational cost than conventional trilinear decomposition-based method.The proposed algorithm requires ( 3 +       ) operations for a iteration, while the trilinear decomposition algorithm needs ( 3 + ) operations [28] for a iteration, where   < ,   < , and   < .
The advantages of the proposed algorithm can be presented as follows.
(1) The proposed algorithm can be regarded as a combination of trilinear model and compressive sensing theory, and it brings much lower computational complexity and much smaller demand for storage capacity.
(2) The proposed algorithm has better 2D-DOA estimation performance than ESPRIT algorithm and close angle estimation performance to TALS algorithm, which will be proved by Figures 6-7.
(3) The proposed algorithm can achieve paired elevation angles and azimuth angles automatically.

Numerical Simulations
In the following simulations, we assume that the numerical simulation results converge when the SSR ≤ 10    Figure 6 shows the 2D-DOA estimation performance comparison of the proposed algorithm, the ESPRIT algorithm, the TALS algorithm, and the CRB for the uniform rectangular array with  = 20,  = 16,  = 100, and  = 3, while Figure 7  performance than the ESPRIT algorithm and close angle estimation to TALS algorithm.The angle estimation performance of the proposed algorithm will be further improved through increasing the compressed parameters   ,   , and   .Figure 8 depicts the 2D-DOA estimation performance of the proposed algorithm with different value of  ( = 16,  = 100, and  = 3), while Figure 9 presents the 2D-DOA estimation performance of the proposed algorithm with different value of .It is clearly shown that the angle estimation performance of our algorithm is gradually improved with the number of antennas increasing.Multiple   antennas improve the angle estimation performance because of diversity gain.Figure 10 presents 2D-DOA estimation performance of the proposed algorithm with different value of  ( = 20,  = 16, and  = 3).It illustrates that the angle estimation performance becomes better in collaboration with  increasing.

Conclusions
In this paper, we have addressed the 2D-DOA estimation problem for rectangular array and have derived a compressive sensing trilinear model-based 2D-DOA estimation algorithm, which can obtain the automatically paired 2D-DOA estimate.The proposed algorithm has better angle estimation performance than ESPRIT algorithm and close angle estimation performance to conventional trilinear decomposition method.Furthermore, the proposed algorithm has lower computational complexity and smaller demand for storage capacity than conventional trilinear decomposition method.
The proposed algorithm can be regarded as a combination of trilinear model and compressive sensing theory, and it brings much lower computational complexity and much smaller demand for storage capacity.

Figure 1 :
Figure 1: The structure of uniform rectangular array.

Figure 3 :
Figure 3: The compression of trilinear model.
Figure 5 depicts the 2D-DOA estimation performance with SNR = 5dB.Figures 4-5 illustrate that our algorithm is effective for 2D-DOA estimation.
(34) θ, −   ) 2 ,(34)where   and   denote the perfect elevation angle and azimuth angle of th source, respectively.φ, and θ, are International Journal of Antennas and Propagation