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

Authors: Huaxin Yu; Xiaofeng Qiu; Xiaofei Zhang; Chenghua Wang; Gang Yang
Journal: International Journal of Antennas and Propagation
 (2015)

---

## Abstract

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 searching but also has better angle estimation performance than estimation of signal parameters via rotational invariance techniques (ESPRIT) algorithm. Furthermore, the proposed algorithm has close angle estimation performance to trilinear decomposition. 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. Numerical simulations present the effectiveness of our approach.


---

## Body

## 1. 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–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 two-dimensional DOA (2D-DOA). 2D-DOA estimation with rectangular array has received considerable attention in the field of array signal processing [7–12]. ESPRIT algorithms in [8–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 two-dimensional MUSIC (2D-MUSIC) algorithm [13] can be used for 2D-DOA estimation. However, the requirement of two-dimensional (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. 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, l
1-SVD method [16] employed l
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 l
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–22]. In the signal processing field, trilinear decomposition can be regarded as a generalization of ESPRIT and joint approximate diagonalization [19–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 low-rank 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. Section2 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.Notation. Bold lower (upper) case letters are adopted to represent vectors (matrices). (
·
)
T, (
·
)
H, (
·
)
∗, (
·
)
-
1, and (
·
)
+ denote transpose, conjugate transpose, conjugate, matrix inversion, and pseudoinverse operations, respectively. I
P stands for a P
×
P identity matrix. D
n
(
A
) denotes a diagonal matrix with the entries of the matrix A’s nth row on the main diagonal. The ith entry of a given column vector g is denoted by g
(
i
). ⊙, ⊗, and ⊕ denote Khatri-Rao product, Kronecker product, and Hadamard product, respectively.IfA is a p-by-q matrix and B is an m-by-n matrix, then the Kronecker product A
⊗
B is the m
p-by-n
q block matrix: (1)
A
⊗
B
=
a
11
B
a
12
B
⋯
a
1
q
B
a
21
B
a
22
B
⋯
a
2
q
B
⋮
⋮
⋱
⋮
a
p
1
B
a
p
2
B
⋯
a
p
q
B
,where a
i
j is the (
i
,
j
) element of the matrix A.IfA is an I-by-F matrix and B is a J-by-F matrix, then the Khatri-Rao product A
⊙
B is the I
J-by-F block matrix: (2)
A
⊙
B
=
a
1
⊗
b
1
,
…
,
a
F
⊗
b
F
,where a
f and b
f are the fth column of the matrices A and B, respectively.IfA is an I-by-J matrix and B is an I-by-J matrix, then the Hadamard product A
⊕
B is (3)
A
⊕
B
=
a
11
b
11
a
12
b
12
⋯
a
1
J
b
1
J
a
21
b
21
a
22
b
22
⋯
a
2
J
b
2
J
⋮
⋮
⋱
⋮
a
I
1
b
I
1
a
I
2
b
I
2
⋯
a
I
J
b
I
J
,where a
i
j and b
i
j are the (
i
,
j
) element of the matrices A and B, respectively.


## 2. Data Model

A rectangular array consisted ofM
×
N elements is shown in Figure 1, where the distance between two adjacent elements is d. 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 K noncoherent or independent sources, and the number of sources is preknown. θ
k and ϕ
k are the elevation angle and the azimuth angle of the kth source, respectively. We assume the sources impinge on the array with different DOAs.Figure 1
The structure of uniform rectangular array.The received signal of the first subarray in the rectangular array isx
1
(
t
)
=
A
x
s
(
t
)
+
n
1
(
t
), where A
x
=
[
a
x
(
ϕ
1
,
θ
1
)
,
a
x
(
ϕ
2
,
θ
2
)
,
…
,
a
x
(
ϕ
K
,
θ
K
)
]
∈
C
M
×
K with a
x
(
ϕ
k
,
θ
k
)
=
[
1
,
e
j
2
π
d
s
i
n
θ
k
c
o
s
ϕ
k
/
λ, …
,
e
j
2
π
(
M
-
1
)
d
s
i
n
θ
k
c
o
s
ϕ
k
/
λ
]
T and λ is the wavelength. n
1
(
t
) is the received additive white Gaussian noise of the first subarray. s
(
t
)
∈
C
K
×
1 is the source vector. The received signal of the nth subarray in the rectangular array is x
n
(
t
)
=
A
x
Φ
n
-
1
s
(
t
)
+
n
n
(
t
), where Φ
=
d
i
a
g
(
e
j
2
π
d
s
i
n
θ
1
s
i
n
ϕ
1
/
λ
,
…
,
e
j
2
π
d
s
i
n
θ
K
s
i
n
ϕ
K
/
λ
) and n
n
(
t
) is the received additive white Gaussian noise of the nth subarray. Therefore, the received signal of the rectangular array is [27](4)
x
t
=
x
1
t
x
2
t
⋮
x
N
t
=
A
x
A
x
Φ
⋮
A
x
Φ
N
-
1
s
t
+
n
1
t
n
2
t
⋮
n
N
t
.The signal x
(
t
)
∈
C
M
N
×
1 in (4) can also be denoted by(5)
x
t
=
A
y
⊙
A
x
s
t
+
n
t
,where A
y
=
[
a
y
(
ϕ
1
,
θ
1
)
,
a
y
(
ϕ
2
,
θ
2
)
,
…
,
a
y
(
ϕ
K
,
θ
K
)
] with a
y
(
ϕ
k
,
θ
k
)
=
[
1
,
e
j
2
π
d
s
i
n
θ
k
s
i
n
ϕ
k
/
λ, …
,
e
j
2
π
(
N
-
1
)
d
s
i
n
θ
k
s
i
n
ϕ
k
/
λ
]
T; n
(
t
)
=
[
n
1
(
t
)
T
,
n
2
(
t
)
T
,
…
,
n
N
(
t
)
T
]
T
∈
C
M
N
×
1. ⊙ denotes Khatri-Rao product.According to the definition of Khatri-Rao product, the signal in (5) can be rewritten as(6)
x
t
=
a
y
ϕ
1
,
θ
1
⊗
a
x
ϕ
1
,
θ
1
,
…
,
i
i
a
y
ϕ
K
,
θ
K
⊗
a
x
ϕ
K
,
θ
K
s
t
+
n
t
,where ⊗ denotes Kronecker product. We collect L samples and define X
=
[
x
(
1
)
,
x
(
2
)
,
…
,
x
(
L
)
]
∈
C
M
N
×
L, which can be expressed as(7)
X
=
A
y
⊙
A
x
S
T
+
N
=
X
1
X
2
⋮
X
N
=
A
x
D
1
A
y
A
x
D
2
A
y
⋮
A
x
D
N
A
y
S
T
+
N
1
N
2
⋮
N
N
,where S
=
[
s
(
1
)
,
s
(
2
)
,
…
,
s
(
L
)
]
T
∈
C
L
×
K is source matrix and N
=
[
n
(
1
)
,
n
(
2
)
,
…
,
n
(
L
)
] is the received additive white Gaussian noise matrix. N
n
∈
C
M
×
L (n
=
1
,
…
,
N) is the noise matrix. Thus, X
n
∈
C
M
×
L in (7) is denoted as(8)
X
n
=
A
x
D
n
A
y
S
T
+
N
n
,
n
=
1,2
,
…
,
N
.Equation (7) can also be denoted with the trilinear model [18, 28](9)
x
m
,
n
,
l
=
∑
k
=
1
K
A
x
m
,
k
A
y
n
,
k
S
l
,
k
+
n
m
,
n
,
l
,
m
=
1
,
…
,
M
,
n
=
1
,
…
,
N
,
l
=
1
,
…
,
L
,where A
x
(
m
,
k
) is the (
m
,
k
) element of the matrix A
x and similarly for the others. n
m
,
n
,
l is noise part. X
n
∈
C
M
×
L (n
=
1
,
…
,
N) 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
m
=
S
D
m
(
A
x
)
A
y
T
+
N
m
y, m
=
1
,
…
,
M, and Z
l
=
A
y
D
l
(
S
)
A
x
T
+
N
l
z, l
=
1
,
…
,
L, where N
m
y and N
l
z are noise matrices. Then, we form the matrices of Y
∈
C
M
L
×
N and Z
∈
C
L
N
×
M:(10)
Y
=
Y
1
Y
2
⋮
Y
M
=
A
x
⊙
S
A
y
T
+
N
y
,
(11)
Z
=
Z
1
Z
2
⋮
Z
L
=
S
⊙
A
y
A
x
T
+
N
z
,where (12)
N
y
=
N
1
y
N
2
y
⋮
N
M
y
,
N
z
=
N
1
z
N
2
z
⋮
N
L
z
.Figure 2
Trilinear model.


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


### 3.1. Trilinear Model Compression

We compress the three-way dataΧ
_
∈
C
M
×
N
×
L into a smaller three-way data Χ
_
′
∈
C
M
′
×
N
′
×
L
′, where M
′
<
M, N
′
<
N, and L
′
<
L. The trilinear model compression processing is shown in Figure 3. We define the compression matrices as U
∈
C
M
×
M
′, V
∈
C
N
×
N
′, and W
∈
C
L
×
L
′, 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, 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].Figure 3
The compression of trilinear model.Then compressX
∈
C
M
N
×
L in (7) to a smaller one as X
′
∈
C
M
′
N
′
×
L
′
(13)
X
′
=
V
T
⊗
U
T
X
W
T
=
V
T
⊗
U
T
A
y
⊙
A
x
S
T
W
T
+
V
T
⊗
U
T
N
W
T
.According to the property of Khatri-Rao product [23], we know(14)
V
T
⊗
U
T
A
y
⊙
A
x
=
V
T
A
y
⊙
U
T
A
x
.Define A
x
′
=
U
T
A
x, A
y
′
=
V
T
A
y, and S
′
=
W
T
S. Equation (11) is also denoted as(15)
X
′
=
A
y
′
⊙
A
x
′
S
′
T
+
N
′
,where N
′
=
(
V
T
⊗
U
T
)
N
W
T. 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(16)
Y
′
=
A
x
′
⊙
S
′
A
y
′
T
+
N
′
y
,
(17)
Z
′
=
S
′
⊙
A
y
′
A
x
′
T
+
N
′
z
,where N
′
y and N
′
z 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.


### 3.2. 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(18)
m
i
n
A
x
′
,
A
y
′
,
S
′
X
′
-
A
y
′
⊙
A
x
′
S
′
T
Fand LS update for the matrix S
′ is(19)
S
^
′
T
=
A
^
y
′
⊙
A
^
x
′
+
X
′
,where A
^
x
′ and A
^
y
′ are previously obtained estimates of A
x
′ and A
y
′, respectively. According to (16), LS fitting is(20)
m
i
n
A
x
′
,
A
y
′
,
S
′
Y
′
-
A
x
′
⊙
S
′
A
y
′
T
Fand LS update for A
y
′ is(21)
A
^
y
′
T
=
A
^
x
′
⊙
S
^
′
+
Y
′
,where A
^
x
′ and S
^
′ stand for the previously obtained estimates of A
x
′ and S
′. Similarly, according to (17), LS fitting is(22)
m
i
n
A
x
′
,
A
y
′
,
S
′
Z
′
-
S
′
⊙
A
y
′
A
x
′
T
F
,where Z
~
′ is the noisy compressed signal. LS update for A
x
′ is(23)
A
^
x
′
T
=
S
^
′
⊙
A
^
y
′
+
Z
′
,where A
^
y
′ and S
^
′ stand for the previously obtained estimates of A
y
′ and S
′, respectively.DefineE
=
X
-
[
A
^
y
′
⊙
A
^
x
′
]
S
^
′
T, where S
^
′, A
^
y
′, and A
^
x
′ present the estimates of S
′, A
y
′, and A
x
′, respectively. The sum of squared residuals (SSR) in trilinear model is defined as S
S
R
=
∑
j
=
1
L
∑
i
=
1
M
N
e
i
j, where e
i
j is the (
i
,
j
) element of the matrix E. According to (19), (21), and (23), the matrices S
′, A
y
′, and A
x
′ are updated with least squares repeatedly until SSR attain apriorthreshold. The we obtain the final estimates S
^
′, A
^
y
′, and A
^
x
′.Theorem 1 (see [22]).
ConsideringX
n
′
=
A
x
′
D
n
(
A
y
′
)
S
′
T, n
=
1
,
…
,
N
′, where A
x
′
∈
C
M
′
×
K, S
′
∈
C
L
′
×
K, and A
y
′
∈
C
N
′
×
K, if(24)
k
A
x
′
+
k
A
y
′
+
k
S
′
T
≥
2
K
+
2
,where k
A is the k-rank of the matrix A [18], then A
x
′, S
′, and A
y
′ are unique up to permutation and scaling of columns.For the different DOAs and independent sources, we havek
A
x
′
=
m
i
n
(
M
′
,
K
), k
A
y
′
=
m
i
n
(
N
′
,
K
), and k
S
′
T
=
m
i
n
(
L
′
,
K
) in the trilinear model in this paper, and then the inequality in (24) becomes(25)
min
⁡
M
′
,
K
+
min
⁡
N
′
,
K
+
min
⁡
L
′
,
K
≥
2
K
+
2
.WhenM
′
≥
K, N
′
≥
K, and L
′
≥
K, the identifiable condition is 1
≤
K
≤
m
i
n
(
M
′
,
N
′
).WhenM
′
≤
K, N
′
≤
K, and L
′
≥
K, the identifiable condition is m
a
x
(
M
′
,
N
′
)
≤
K
≤
M
′
+
N
′
-
2. Hence, the proposed algorithm is effective when K
≤
M
′
+
N
′
-
2 and the maximum number of sources that can be identified is M
′
+
N
′
-
2.After the trilinear decomposition, we obtain the estimates of the loading matrices(26)
A
^
x
′
=
A
x
′
Π
Δ
1
+
E
1
,
S
^
′
=
S
′
Π
Δ
2
+
E
2
,
A
^
y
′
=
A
y
′
Π
Δ
3
+
E
3
,where Π is a permutation matrix, and Δ
1, Δ
2, Δ
3 note for the diagonal scaling matrices satisfying Δ
1
Δ
2
Δ
3
=
I
K. E
1, E
2, and E
3 are estimation error matrices. After the trilinear decomposition, the estimates of A
x
′, A
y
′, 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.


### 3.3. Angle Estimation via Sparse Recovery

Usea
^
x
k
′ and a
^
y
k
′ to denote the kth column of estimates A
^
x
′ and A
^
y
′, respectively. According to the compression matrices, we have(27a)
a
^
x
k
′
=
U
T
∂
x
k
a
x
k
+
n
x
k
,
k
=
1
,
…
,
K
,
(27b)
a
^
y
k
′
=
V
T
∂
y
k
a
y
k
+
n
y
k
,
k
=
1
,
…
,
K
,where a
x
k and a
y
k are the kth column of A
x, A
y, respectively. n
x
k and n
y
k are the corresponding noise, respectively. ∂
x
k and ∂
x
k are the scaling coefficients. Construct two Vandermonde matrices A
s
x
∈
C
M
×
P and A
s
y
∈
C
N
×
P (P
≫
M, P
≫
N) composed of steering vectors corresponding to each potential source location as its columns:(28a)
A
s
x
=
a
s
x
1
,
a
s
x
2
…
,
a
s
x
P
=
1
1
⋯
1
e
j
2
π
d
g
1
/
λ
e
j
2
π
d
g
2
/
λ
⋮
e
j
2
π
d
g
P
/
λ
⋮
⋮
⋮
⋮
e
j
2
π
M
-
1
d
g
1
/
λ
e
j
2
π
M
-
1
d
g
2
/
λ
⋯
e
j
2
π
M
-
1
d
g
P
/
λ
,
(28b)
A
s
y
=
a
s
y
1
,
a
s
y
2
…
,
a
s
y
P
=
1
1
⋯
1
e
j
2
π
d
g
1
/
λ
e
j
2
π
d
g
2
/
λ
⋮
e
j
2
π
d
g
P
/
λ
⋮
⋮
⋮
⋮
e
j
2
π
N
-
1
d
g
1
/
λ
e
j
2
π
N
-
1
d
g
2
/
λ
⋯
e
j
2
π
N
-
1
d
g
P
/
λ
,where g is a sampling vector and its pth elements is g
(
p
)
=
-
1
+
2
p
/
P, p
=
1,2
,
…
,
P. The matrices A
s
x and A
s
y can be regarded as the completed dictionaries. Then (27a)-(27b) can be expressed as(29a)
a
^
x
k
′
=
U
T
A
s
x
x
s
+
n
x
k
,
k
=
1
,
…
,
K
,
(29b)
a
^
y
k
′
=
V
T
A
y
x
y
s
+
n
y
k
,
k
=
1
,
…
,
K
,where x
s and y
s are sparse. The estimates of x
s and y
s can be obtained via l
0-norm constraint:(30a)
m
i
n
a
^
x
k
′
-
U
T
A
s
x
x
s
2
2
,
s
t
.
u
u
n
i
l
l
x
s
0
=
1
,
(30b)
m
i
n
a
^
y
k
′
-
V
T
A
s
y
y
s
2
2
,
s
t
.
u
u
n
i
l
l
y
s
0
=
1
,where ·
0 denotes the l
0-norm. x
s
0
=
1; that is to say, there is only one nonzero element in the vector x
s, similar to y
s
0
=
1. We can use the OMP recovery method [26] to find the nonzero element in x
s or y
s. 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
s and y
s, respectively, noted as p
x and p
y. According to the corresponding columns in A
s
x and A
s
y, we obtain g
(
p
x
) and g
(
p
y
), which are estimates of s
i
n
θ
k
c
o
s
ϕ
k and s
i
n
θ
k
s
i
n
ϕ
k. We define r
k
=
g
(
p
x
)
+
j
g
(
p
y
), and then the elevation angles and azimuth angles can be obtained via (31a)
θ
^
k
=
sin
-
1
⁡
a
b
s
r
k
,
k
=
1
,
…
,
K
,
(31b)
ϕ
^
k
=
a
n
g
l
e
r
k
,
k
=
1
,
…
,
K
,where a
b
s
(
·
) is the modulus value symbol and a
n
g
l
e
(
·
) is to get the angle of an imaginary number. As the columns of the estimated matrices A
^
x
′ and A
^
y
′ are automatically paired, then the estimated elevation angles and azimuth angles can be paired automatically.


### 3.4. 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
M
×
N
×
L, then compress the three-way matrix into a much smaller three-way matrix Χ
-
′
∈
C
M
′
×
N
′
×
L
′ via the compression matrices U
∈
C
M
×
M
′, V
∈
C
N
×
N
′, and W
∈
C
L
×
L
′.Step 2.
Perform trilinear decomposition through TALS algorithm for the compressed three-way matrix to obtain the estimation ofA
x
′, A
y
′, and S
′.Step 3.
Estimate the sparse vectors.Step 4.
Estimate 2D-DOA via (31a)-(31b).Remark A. Because the trilinear decomposition brings the same permutation ambiguity for the estimates A
x
′, A
y
′, 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, s
i
n
θ
k
c
o
s
ϕ
k (or s
i
n
θ
k
s
i
n
ϕ
k) 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 s
i
n
θ
k
c
o
s
ϕ
k and s
i
n
θ
k
s
i
n
ϕ
k, respectively. Afterwards, the elevation and azimuth angles are estimated through the estimates of s
i
n
θ
k
c
o
s
ϕ
k and s
i
n
θ
k
s
i
n
ϕ
k.Remark C. If the number of sources K 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 lower-complexity 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.


### 3.5. Complexity Analysis and CRB

The proposed algorithm has much lower computational cost than conventional trilinear decomposition-based method. The proposed algorithm requiresO
(
K
3
+
M
′
N
′
L
′
K
) operations for a iteration, while the trilinear decomposition algorithm needs O
(
K
3
+
M
N
L
K
) operations [28] for a iteration, where M
′
<
M, N
′
<
N, and L
′
<
L.We defineA
=
[
a
y
(
θ
1
)
⊗
a
x
(
θ
1
)
,
…
,
a
y
(
θ
K
)
⊗
a
x
(
θ
K
)
]. According to [37], we can derive the CRB(32)
C
R
B
=
σ
2
2
L
R
e
D
H
Π
A
⊥
D
⊕
P
^
s
T
-
1
,where L denotes the number of samples, a
k is the kth column ofA, and P
^
s
=
(
1
/
L
)
∑
t
=
1
L
s
t
s
H
t. σ
2 is the noise power. Π
A
⊥
=
I
M
N
-
A
(
A
H
A
)
-
1
A
H and(33)
D
=
∂
a
1
∂
θ
1
,
∂
a
2
∂
θ
2
,
…
,
∂
a
K
∂
θ
K
,
∂
a
1
∂
ϕ
1
,
∂
a
2
∂
ϕ
2
,
…
,
∂
a
K
∂
ϕ
K
.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 Figures6-7.
(3)
The proposed algorithm can achieve paired elevation angles and azimuth angles automatically.


## 3.1. Trilinear Model Compression

We compress the three-way dataΧ
_
∈
C
M
×
N
×
L into a smaller three-way data Χ
_
′
∈
C
M
′
×
N
′
×
L
′, where M
′
<
M, N
′
<
N, and L
′
<
L. The trilinear model compression processing is shown in Figure 3. We define the compression matrices as U
∈
C
M
×
M
′, V
∈
C
N
×
N
′, and W
∈
C
L
×
L
′, 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, 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].Figure 3
The compression of trilinear model.Then compressX
∈
C
M
N
×
L in (7) to a smaller one as X
′
∈
C
M
′
N
′
×
L
′
(13)
X
′
=
V
T
⊗
U
T
X
W
T
=
V
T
⊗
U
T
A
y
⊙
A
x
S
T
W
T
+
V
T
⊗
U
T
N
W
T
.According to the property of Khatri-Rao product [23], we know(14)
V
T
⊗
U
T
A
y
⊙
A
x
=
V
T
A
y
⊙
U
T
A
x
.Define A
x
′
=
U
T
A
x, A
y
′
=
V
T
A
y, and S
′
=
W
T
S. Equation (11) is also denoted as(15)
X
′
=
A
y
′
⊙
A
x
′
S
′
T
+
N
′
,where N
′
=
(
V
T
⊗
U
T
)
N
W
T. 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(16)
Y
′
=
A
x
′
⊙
S
′
A
y
′
T
+
N
′
y
,
(17)
Z
′
=
S
′
⊙
A
y
′
A
x
′
T
+
N
′
z
,where N
′
y and N
′
z 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.


## 3.2. 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(18)
m
i
n
A
x
′
,
A
y
′
,
S
′
X
′
-
A
y
′
⊙
A
x
′
S
′
T
Fand LS update for the matrix S
′ is(19)
S
^
′
T
=
A
^
y
′
⊙
A
^
x
′
+
X
′
,where A
^
x
′ and A
^
y
′ are previously obtained estimates of A
x
′ and A
y
′, respectively. According to (16), LS fitting is(20)
m
i
n
A
x
′
,
A
y
′
,
S
′
Y
′
-
A
x
′
⊙
S
′
A
y
′
T
Fand LS update for A
y
′ is(21)
A
^
y
′
T
=
A
^
x
′
⊙
S
^
′
+
Y
′
,where A
^
x
′ and S
^
′ stand for the previously obtained estimates of A
x
′ and S
′. Similarly, according to (17), LS fitting is(22)
m
i
n
A
x
′
,
A
y
′
,
S
′
Z
′
-
S
′
⊙
A
y
′
A
x
′
T
F
,where Z
~
′ is the noisy compressed signal. LS update for A
x
′ is(23)
A
^
x
′
T
=
S
^
′
⊙
A
^
y
′
+
Z
′
,where A
^
y
′ and S
^
′ stand for the previously obtained estimates of A
y
′ and S
′, respectively.DefineE
=
X
-
[
A
^
y
′
⊙
A
^
x
′
]
S
^
′
T, where S
^
′, A
^
y
′, and A
^
x
′ present the estimates of S
′, A
y
′, and A
x
′, respectively. The sum of squared residuals (SSR) in trilinear model is defined as S
S
R
=
∑
j
=
1
L
∑
i
=
1
M
N
e
i
j, where e
i
j is the (
i
,
j
) element of the matrix E. According to (19), (21), and (23), the matrices S
′, A
y
′, and A
x
′ are updated with least squares repeatedly until SSR attain apriorthreshold. The we obtain the final estimates S
^
′, A
^
y
′, and A
^
x
′.Theorem 1 (see [22]).
ConsideringX
n
′
=
A
x
′
D
n
(
A
y
′
)
S
′
T, n
=
1
,
…
,
N
′, where A
x
′
∈
C
M
′
×
K, S
′
∈
C
L
′
×
K, and A
y
′
∈
C
N
′
×
K, if(24)
k
A
x
′
+
k
A
y
′
+
k
S
′
T
≥
2
K
+
2
,where k
A is the k-rank of the matrix A [18], then A
x
′, S
′, and A
y
′ are unique up to permutation and scaling of columns.For the different DOAs and independent sources, we havek
A
x
′
=
m
i
n
(
M
′
,
K
), k
A
y
′
=
m
i
n
(
N
′
,
K
), and k
S
′
T
=
m
i
n
(
L
′
,
K
) in the trilinear model in this paper, and then the inequality in (24) becomes(25)
min
⁡
M
′
,
K
+
min
⁡
N
′
,
K
+
min
⁡
L
′
,
K
≥
2
K
+
2
.WhenM
′
≥
K, N
′
≥
K, and L
′
≥
K, the identifiable condition is 1
≤
K
≤
m
i
n
(
M
′
,
N
′
).WhenM
′
≤
K, N
′
≤
K, and L
′
≥
K, the identifiable condition is m
a
x
(
M
′
,
N
′
)
≤
K
≤
M
′
+
N
′
-
2. Hence, the proposed algorithm is effective when K
≤
M
′
+
N
′
-
2 and the maximum number of sources that can be identified is M
′
+
N
′
-
2.After the trilinear decomposition, we obtain the estimates of the loading matrices(26)
A
^
x
′
=
A
x
′
Π
Δ
1
+
E
1
,
S
^
′
=
S
′
Π
Δ
2
+
E
2
,
A
^
y
′
=
A
y
′
Π
Δ
3
+
E
3
,where Π is a permutation matrix, and Δ
1, Δ
2, Δ
3 note for the diagonal scaling matrices satisfying Δ
1
Δ
2
Δ
3
=
I
K. E
1, E
2, and E
3 are estimation error matrices. After the trilinear decomposition, the estimates of A
x
′, A
y
′, 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.


## 3.3. Angle Estimation via Sparse Recovery

Usea
^
x
k
′ and a
^
y
k
′ to denote the kth column of estimates A
^
x
′ and A
^
y
′, respectively. According to the compression matrices, we have(27a)
a
^
x
k
′
=
U
T
∂
x
k
a
x
k
+
n
x
k
,
k
=
1
,
…
,
K
,
(27b)
a
^
y
k
′
=
V
T
∂
y
k
a
y
k
+
n
y
k
,
k
=
1
,
…
,
K
,where a
x
k and a
y
k are the kth column of A
x, A
y, respectively. n
x
k and n
y
k are the corresponding noise, respectively. ∂
x
k and ∂
x
k are the scaling coefficients. Construct two Vandermonde matrices A
s
x
∈
C
M
×
P and A
s
y
∈
C
N
×
P (P
≫
M, P
≫
N) composed of steering vectors corresponding to each potential source location as its columns:(28a)
A
s
x
=
a
s
x
1
,
a
s
x
2
…
,
a
s
x
P
=
1
1
⋯
1
e
j
2
π
d
g
1
/
λ
e
j
2
π
d
g
2
/
λ
⋮
e
j
2
π
d
g
P
/
λ
⋮
⋮
⋮
⋮
e
j
2
π
M
-
1
d
g
1
/
λ
e
j
2
π
M
-
1
d
g
2
/
λ
⋯
e
j
2
π
M
-
1
d
g
P
/
λ
,
(28b)
A
s
y
=
a
s
y
1
,
a
s
y
2
…
,
a
s
y
P
=
1
1
⋯
1
e
j
2
π
d
g
1
/
λ
e
j
2
π
d
g
2
/
λ
⋮
e
j
2
π
d
g
P
/
λ
⋮
⋮
⋮
⋮
e
j
2
π
N
-
1
d
g
1
/
λ
e
j
2
π
N
-
1
d
g
2
/
λ
⋯
e
j
2
π
N
-
1
d
g
P
/
λ
,where g is a sampling vector and its pth elements is g
(
p
)
=
-
1
+
2
p
/
P, p
=
1,2
,
…
,
P. The matrices A
s
x and A
s
y can be regarded as the completed dictionaries. Then (27a)-(27b) can be expressed as(29a)
a
^
x
k
′
=
U
T
A
s
x
x
s
+
n
x
k
,
k
=
1
,
…
,
K
,
(29b)
a
^
y
k
′
=
V
T
A
y
x
y
s
+
n
y
k
,
k
=
1
,
…
,
K
,where x
s and y
s are sparse. The estimates of x
s and y
s can be obtained via l
0-norm constraint:(30a)
m
i
n
a
^
x
k
′
-
U
T
A
s
x
x
s
2
2
,
s
t
.
u
u
n
i
l
l
x
s
0
=
1
,
(30b)
m
i
n
a
^
y
k
′
-
V
T
A
s
y
y
s
2
2
,
s
t
.
u
u
n
i
l
l
y
s
0
=
1
,where ·
0 denotes the l
0-norm. x
s
0
=
1; that is to say, there is only one nonzero element in the vector x
s, similar to y
s
0
=
1. We can use the OMP recovery method [26] to find the nonzero element in x
s or y
s. 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
s and y
s, respectively, noted as p
x and p
y. According to the corresponding columns in A
s
x and A
s
y, we obtain g
(
p
x
) and g
(
p
y
), which are estimates of s
i
n
θ
k
c
o
s
ϕ
k and s
i
n
θ
k
s
i
n
ϕ
k. We define r
k
=
g
(
p
x
)
+
j
g
(
p
y
), and then the elevation angles and azimuth angles can be obtained via (31a)
θ
^
k
=
sin
-
1
⁡
a
b
s
r
k
,
k
=
1
,
…
,
K
,
(31b)
ϕ
^
k
=
a
n
g
l
e
r
k
,
k
=
1
,
…
,
K
,where a
b
s
(
·
) is the modulus value symbol and a
n
g
l
e
(
·
) is to get the angle of an imaginary number. As the columns of the estimated matrices A
^
x
′ and A
^
y
′ are automatically paired, then the estimated elevation angles and azimuth angles can be paired automatically.


## 3.4. 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
M
×
N
×
L, then compress the three-way matrix into a much smaller three-way matrix Χ
-
′
∈
C
M
′
×
N
′
×
L
′ via the compression matrices U
∈
C
M
×
M
′, V
∈
C
N
×
N
′, and W
∈
C
L
×
L
′.Step 2.
Perform trilinear decomposition through TALS algorithm for the compressed three-way matrix to obtain the estimation ofA
x
′, A
y
′, and S
′.Step 3.
Estimate the sparse vectors.Step 4.
Estimate 2D-DOA via (31a)-(31b).Remark A. Because the trilinear decomposition brings the same permutation ambiguity for the estimates A
x
′, A
y
′, 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, s
i
n
θ
k
c
o
s
ϕ
k (or s
i
n
θ
k
s
i
n
ϕ
k) 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 s
i
n
θ
k
c
o
s
ϕ
k and s
i
n
θ
k
s
i
n
ϕ
k, respectively. Afterwards, the elevation and azimuth angles are estimated through the estimates of s
i
n
θ
k
c
o
s
ϕ
k and s
i
n
θ
k
s
i
n
ϕ
k.Remark C. If the number of sources K 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 lower-complexity 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.


## 3.5. Complexity Analysis and CRB

The proposed algorithm has much lower computational cost than conventional trilinear decomposition-based method. The proposed algorithm requiresO
(
K
3
+
M
′
N
′
L
′
K
) operations for a iteration, while the trilinear decomposition algorithm needs O
(
K
3
+
M
N
L
K
) operations [28] for a iteration, where M
′
<
M, N
′
<
N, and L
′
<
L.We defineA
=
[
a
y
(
θ
1
)
⊗
a
x
(
θ
1
)
,
…
,
a
y
(
θ
K
)
⊗
a
x
(
θ
K
)
]. According to [37], we can derive the CRB(32)
C
R
B
=
σ
2
2
L
R
e
D
H
Π
A
⊥
D
⊕
P
^
s
T
-
1
,where L denotes the number of samples, a
k is the kth column ofA, and P
^
s
=
(
1
/
L
)
∑
t
=
1
L
s
t
s
H
t. σ
2 is the noise power. Π
A
⊥
=
I
M
N
-
A
(
A
H
A
)
-
1
A
H and(33)
D
=
∂
a
1
∂
θ
1
,
∂
a
2
∂
θ
2
,
…
,
∂
a
K
∂
θ
K
,
∂
a
1
∂
ϕ
1
,
∂
a
2
∂
ϕ
2
,
…
,
∂
a
K
∂
ϕ
K
.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 Figures6-7.
(3)
The proposed algorithm can achieve paired elevation angles and azimuth angles automatically.


## 4. Numerical Simulations

In the following simulations, we assume that the numerical simulation results converge when theS
S
R
≤
1
0
-
8. M, N, L, and K denote the number of antennas in x-axis, number of antennas in y-axis, samples, and sources, respectively. And we compress the parameters M, N, L to M
′, N
′, and L
′ (usually set M
=
16, N
=
20, L
=
100, and M
′
=
N
′
=
L
′
=
5 in numerical simulations). d
=
λ
/
2 is considered in the simulation. We present 1000 Monte Carlo simulations to assess the angle estimation performance of the proposed algorithm. Define root mean squared error (RMSE) as (34)
R
M
S
E
=
1
K
∑
k
=
1
K
1
1000
∑
l
=
1
1000
ϕ
^
k
,
l
-
ϕ
k
2
+
θ
^
k
,
l
-
θ
k
2
,where ϕ
k and θ
k denote the perfect elevation angle and azimuth angle of kth source, respectively. ϕ
^
k
,
l and θ
^
k
,
l are the estimates of ϕ
k and θ
k in the lth Monte Carlo trail. Assume that there are 3 noncoherent sources located at angles (
ϕ
1
,
θ
1
)
=
(
5
∘
,
1
0
∘
), (
ϕ
2
,
θ
2
)
=
(
1
5
∘
,
2
0
∘
), and (
ϕ
3
,
θ
3
)
=
(
4
0
∘
,
3
0
∘
).Figure4 presents the 2D-DOA estimation of the proposed algorithm for uniform rectangular array with N
=
20, M
=
16, L
=
100, K
=
3, and S
N
R
=
-
10 dB. Figure 5 depicts the 2D-DOA estimation performance with S
N
R
=
5 dB. Figures 4-5 illustrate that our algorithm is effective for 2D-DOA estimation.Figure 4
2D-DOA estimation of our algorithm inS
N
R
=
-
10 dB (N
=
20, M
=
16, L
=
100, and K
=
3).Figure 5
2D-DOA estimation of our algorithm inS
N
R
=
5 dB (N
=
20, M
=
16, L
=
100, and K
=
3).Figure 6
2D-DOA estimation performance comparison (N
=
20, M
=
16, L
=
100, and K
=
3).Figure 7
2D-DOA estimation performance comparison (N
=
16, M
=
16, L
=
200, and K
=
3).Figure6 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 N
=
20, M
=
16, L
=
100, and K
=
3, while Figure 7 depicts the 2D-DOA estimation performance comparison with N
=
16, M
=
16, L
=
200, and K
=
3. It is indicated that our algorithm has better angle estimation 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 M
′, N
′, and L
′.Figure8 depicts the 2D-DOA estimation performance of the proposed algorithm with different value of N (M
=
16, L
=
100, and K
=
3), while Figure 9 presents the 2D-DOA estimation performance of the proposed algorithm with different value of M. 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 8
Angle estimation performance of our algorithm with differentN (M
=
16, L
=
100, and K
=
3).Figure 9
Angle estimation performance of our algorithm with differentM (N
=
16, L
=
100, and K
=
3).Figure10 presents 2D-DOA estimation performance of the proposed algorithm with different value of L (N
=
20, M
=
16, and K
=
3). It illustrates that the angle estimation performance becomes better in collaboration with L increasing.Figure 10
Angle estimation performance of our algorithm with differentL (N
=
16, M
=
20, and K
=
3).


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



---

*Source: 297572-2015-05-19.xml*