# High Accuracy 2D-DOA Estimation for Conformal Array Using PARAFAC

Authors: Liangtian Wan; Weijian Si; Lutao Liu; Zuoxi Tian; Naixing Feng
Journal: International Journal of Antennas and Propagation
 (2014)

---

## Abstract

Due to the polarization diversity (PD) of element patterns caused by the varying curvature of the conformal carrier, the conventional direction-of-arrival (DOA) estimation algorithms could not be applied to the conformal array. In order to describe the PD of conformal array, the polarization parameter is considered in the snapshot data model. The paramount difficulty for DOA estimation is the coupling between the angle information and polarization parameter. Based on the characteristic of the cylindrical conformal array, the decoupling between the polarization parameter and DOA can be realized with a specially designed array structure. 2D-DOA estimation of the cylindrical conformal array is accomplished via parallel factor analysis (PARAFAC) theory. To avoid parameter pairing problem, the algorithm forms a PARAFAC model of the covariance matrices in the covariance domain. The proposed algorithm can also be generalized to other conformal array structures and nonuniform noise scenario. Cramer-Rao bound (CRB) is derived and simulation results with the cylindrical conformal array are presented to verify the performance of the proposed algorithm.


---

## Body

## 1. Introduction

Conformal arrays mounted on curved surfaces are commonly applied in various areas such as radar, sonar, and wireless communication [1]. The conformal antennas could fulfill specific aerodynamics, space-saving, elimination of random-induced bore-sight error, potential increase in available aperture, and so on [2, 3]. Most researches focus on the design of antenna configuration [4–6], the transform between the local and global coordinate [7–10], and the pattern synthesis of conformal array [11–15].The conventional direction-of-arrival (DOA) estimation algorithms, for example, the multiple signal classification (MUSIC) [16] and the estimation of signal parameters via rotational invariance techniques (ESPRIT) [17], are not applicable due to the varying curvature of the conformal array. The “shadow effect” of the conformal array is caused by the metallic shelter, leading to the condition that not all elements have the ability to receive the signal. As a result of the incomplete steering vector, common DOA estimation algorithms cannot be used for conformal array. Recently, DOA estimation algorithms for conformal array have been proposed for high resolution [18–22]. With the help of fourth-order cumulant and ESPRIT algorithm, a blind DOA estimation algorithm is proposed in [18]. Based on the mathematical technique of geometric algebra, 2D-DOA and polarization parameter estimations are completed by exploiting iterative ESPRIT algorithm [19]. The state space and propagator method are used for joint frequency and 2D-DOA estimations for the cylindrical conformal array [20]. However, the main drawback of the algorithms in [18–20] is the need of parameter pairing.Parallel factor (PARAFAC) analysis attracted the attention of researchers when it was originally introduced in array signal processing in 2000 [21, 22]. It has been widely used for low-rank decomposition of three and higher way array. In this paper, a high accuracy 2D-DOA estimation algorithm for the conformal array is proposed using PARAFAC theory. The proposed algorithm does not need parameter pairing and performs well when the signal-to-noise ratio (SNR) is low or in situations where high estimation accuracy is needed. Simulation results demonstrate the efficiency and accuracy of the proposed algorithm.The rest of this paper is structured as follows. Section2 introduces the snapshot data model. Section 3 elaborates the design of the conformal array structure. Section 4 contains the core contributions of this paper, in which the PARAFAC model is used for high accuracy 2D-DOA estimation. Section 5 extends the proposed algorithm to more general cases. Section 6 analyses the computational complexity of the algorithm. Section 7 derives the Cramer-Rao bound (CRB). Section 8 discusses the array design. Section 9 presents the simulation result. Conclusions are drawn in Section 10.


## 2. The Snapshot Data Model

In order to achieve accurate parameter estimation, the snapshot data model of the conformal array must be precisely established. The pattern of the elements is different due to the varying curvature of carriers [8]. The mathematic model of the conformal array constructed in [8, 18] is adopted in this paper. The far-field narrowband signal impinging on the conformal array with elevationθand azimuthφis shown in Figure 1(a). The steering vector takes the following form:
(1)
a
(
θ
,
φ
)
=
[
h
1
e
-
j
2
π
(
(
p
1
•
u
)
/
λ
)
,
h
2
e
-
j
2
π
(
(
p
2
•
u
)
/
λ
)
,
…
,
m
h
2
m
+
2
e
-
j
2
π
(
(
p
2
m
+
2
•
u
)
/
λ
)
]
T
,
(2)
h
i
=
(
g
i
θ
2
+
g
i
φ
2
)
1
/
2
(
k
θ
2
+
k
φ
2
)
1
/
2
cos
⁡
(
θ
i
g
k
)
=
|
g
i
|
|
p
l
|
cos
⁡
(
θ
i
g
k
)
=
g
i
•
p
l
=
g
i
θ
k
θ
+
g
i
φ
k
φ
,
where P
i represents the position vector of the ith element, u is the propagation vector, and λ is wavelength of the incident signal. As shown in Figure 1(b), u
θ and u
φ are orthogonal unit vectors and k
θ and k
φ are the polarization components of the incident signal. h
i denotes the ith element’s response of the incident signal, g
i is the pattern of the ith element, and p
l is the electric field direction of the incident signal. θ
i
g
k is the angle between vector g
i and p
l. (
•
)
T denotes the transpose of matrix (
•
).(a) The direction vectoru impinges on the array. (b) The ith sensors response.
(a)
(b)Assuming thatr incident signals impinge on the designed conformal array, the snapshot data model can be represented as
(3)
X
(
n
)
=
G
•
A
S
(
n
)
+
N
(
n
)
=
(
G
θ
•
A
θ
K
θ
+
G
φ
•
A
φ
K
φ
)
S
(
n
)
+
N
(
n
)
=
B
S
(
n
)
+
N
(
n
)
,
G
θ
=
[
g
θ
(
θ
1
,
φ
1
)
,
g
θ
(
θ
2
,
φ
2
)
,
…
,
g
θ
(
θ
r
,
φ
r
)
]
,
G
φ
=
[
g
φ
(
θ
1
,
φ
1
)
,
g
φ
(
θ
2
,
φ
2
)
,
…
,
g
φ
(
θ
r
,
φ
r
)
]
,
A
θ
=
[
a
θ
(
θ
1
,
φ
1
)
,
a
θ
(
θ
2
,
φ
2
)
,
…
,
a
θ
(
θ
r
,
φ
r
)
]
,
A
φ
=
[
a
φ
(
θ
1
,
φ
1
)
,
a
φ
(
θ
2
,
φ
2
)
,
…
,
a
φ
(
θ
r
,
φ
r
)
]
,
K
θ
=
diag
⁡
(
k
1
θ
,
k
2
θ
,
…
,
k
r
θ
)
,
K
φ
=
diag
⁡
(
k
1
φ
,
k
2
φ
,
…
,
k
r
φ
)
,
S
(
n
)
=
[
s
1
(
n
)
,
s
2
(
n
)
,
…
,
s
r
(
n
)
]
T
,
W
(
n
)
=
[
w
1
(
n
)
,
w
2
(
n
)
,
…
,
w
r
(
n
)
]
T
,
where G denotes the pattern matrix and A denotes the manifold matrix. K
θ and K
φ are the diagonal matrices whose diagonal entries are k
1
θ
,
k
2
θ
,
…
,
k
r
θ and k
1
φ
,
k
2
φ
,
…
,
k
r
φ, respectively.k
i
θandk
i
φare the components of ith signal’s polarization parameters along the orthogonal unit vector u
θ and u
φ, respectively. S
(
n
) denotes the signal vector. W
(
n
) is additive Gaussian white noise with the covariance matrix
(4)
E
{
W
(
n
)
W
(
n
)
H
}
=
Q
=
σ
2
I
.
(
•
)
H denotes conjugate transpose of matrix (
•
). Collecting N snapshots,
(5)
X
=
B
S
+
W
,
where S is the r
×
N signal matrix and W is the (
2
m
+
2
)
×
N noise matrix.The definition of the pattern is in the local coordinate, meaning that transformation must be performed from the global coordinate to the local coordinate [8]. In addition, the polarization parameters couple with the signal parameters; thus the decoupling is definitely necessary for the 2D-DOA estimation of the conformal array.


## 3. The Structure of the Array

Since the signal’s amplitude of the sensors received will decrease because of the “shadow effect,” the whole conformal array is divided into three subarrays, and each subarray covers12
0
°. For the single curvature and symmetry of the cylinder, each subarray possesses and remains the same structure as well as the parameter estimation mechanism. Thus, this paper only considers one subarray.As shown in Figure2, the distance between two sensors which are in the same intersecting surface is λ
/
4 and the distance between two adjacent intersecting surfaces is 5
λ.Figure 2
The structure of cylindrical conformal array.As shown in Figure2, 1
~
m constitute array 1, 2
~
m
+
1 constitute array 2, m
+
2
~
2
m
+
1 constitute array 3, and m
+
3
~
2
m
+
2 constitute array 4 (1
~
2
m
+
2 are sensors of the array). The distance vector between array 1 and array 2 is Δ
P
1 and d
1
=
|
Δ
P
1
|
=
λ
/
4. The distance vector between array 1 and array 3 is Δ
P
2 and d
2
=
|
Δ
P
2
|
=
λ
/
4 as shown in Figure 3. The sensors’ patterns of the same generatrix are identical. The sensors 1
~
m
+
1 possess the same pattern g
1, and the sensors m
+
2
~
2
m
+
2 possess the same pattern g
2. The decoupling between polarization parameter and angle information can be realized by taking full advantage of the characteristic of the array structure.(a) The distance vectorΔ
P
1. (b) The distance vector Δ
P
2.
(a)
(b)Assume thatX
1, X
2, X
3, and X
4 represent the data that array 1, array 2, array 3, and array 4 receive, respectively. The origin of the coordinate is the reference point. So, the received data can be expressed as
(6)
X
1
=
B
S
+
N
1
,
X
2
=
B
Ψ
1
S
+
N
2
,
X
3
=
B
Ψ
2
S
+
N
3
,
X
4
=
B
Ψ
1
Ψ
2
S
+
N
4
.
The covariance matrices among the received data are
(7)
R
1
=
E
{
X
1
X
1
H
}
=
B
R
s
B
H
+
Q
1
,
R
2
=
E
{
X
2
X
1
H
}
=
B
Ψ
1
R
s
B
H
+
Q
2
,
R
3
=
E
{
X
3
X
1
H
}
=
B
Ψ
2
R
s
B
H
+
Q
3
,
R
4
=
E
{
X
4
X
1
H
}
=
B
Ψ
1
Ψ
2
R
s
B
H
+
Q
4
,
where R
s
=
diag
⁡
{
s
1
2
,
…
,
s
r
2
} represents the signal covariance matrix and the noise covariance matrices are Q
1
–
Q
4, respectively.Consider(8)
Ψ
1
=
diag
⁡
[
exp
⁡
(
-
j
ω
11
)
,
exp
⁡
(
-
j
ω
12
)
,
…
,
exp
⁡
(
-
j
ω
1
r
)
]
,
(9)
Ψ
2
=
diag
⁡
[
h
3
(
θ
1
,
φ
1
)
h
1
(
θ
1
,
φ
1
)
exp
⁡
(
-
j
ω
21
)
,
m
m
m
i
h
3
(
θ
2
,
φ
2
)
h
1
(
θ
2
,
φ
2
)
exp
⁡
(
-
j
ω
22
)
,
…
,
m
m
m
i
h
3
(
θ
r
,
φ
r
)
h
1
(
θ
r
,
φ
r
)
exp
⁡
(
-
j
ω
2
r
)
]
,
(10)
ω
1
i
=
(
2
π
λ
)
Δ
P
1
•
u
i
=
(
2
π
d
1
λ
)
×
[
sin
(
θ
Δ
P
1
)
cos
⁡
(
φ
Δ
P
1
)
sin
(
θ
i
)
cos
⁡
(
φ
i
)
+
sin
(
θ
Δ
P
1
)
sin
(
φ
Δ
P
1
)
sin
(
θ
i
)
sin
(
φ
i
)
+
cos
⁡
(
θ
Δ
P
1
)
cos
⁡
(
θ
i
)
]
,
(11)
ω
2
i
=
(
2
π
λ
)
Δ
P
2
•
u
i
=
(
2
π
d
2
λ
)
×
[
sin
(
θ
Δ
P
2
)
cos
⁡
(
φ
Δ
P
2
)
sin
(
θ
i
)
cos
⁡
(
φ
i
)
+
sin
(
θ
Δ
P
2
)
sin
(
φ
Δ
P
2
)
sin
(
θ
i
)
sin
(
φ
i
)
+
cos
⁡
(
θ
Δ
P
2
)
cos
⁡
(
θ
i
)
]
,
where θ
Δ
P
i and φ
Δ
P
i represent the elevation and azimuth of the distance vector in the global coordinate. The array structure is shown in Figures 2 and 3 in which the distance vector Δ
P
1 is parallel to the Z-axis. The elevation and azimuth of Δ
P
1 in the global coordinate are
(12)
θ
Δ
P
1
=
0
,
φ
Δ
P
1
=
π
2
.
Δ
P
2 is parallel to the X-axis. The elevation and azimuth of Δ
P
2 in the global coordinate are
(13)
θ
Δ
P
2
=
π
2
,
φ
Δ
P
2
=
0
.In real applications, the ideal covariance matrix has to be replaced by the sample covariance matrix which is obtained by a finite number of snapshotsR
k  (
k
=
1
,
…
,
4
); that is,
(14)
R
k
=
1
N
∑
k
=
1
N
‍
X
(
n
)
X
(
n
)
H
=
1
N
X
X
H
.


## 4. The DOA Estimation Based on PARAFAC

Parallel factor (PARAFAC) analysis is a method of multiple data analysis, which has been first introduced in psychometrics. It has been used in many fields such as statistics, arithmetic complexity, and chemometrics. In this section, first, the PARAFAC theory is introduced briefly; second, based on the PARAFAC model, the multiple way array is constructed and the identifiability of inherently unresolvable source permutation and scaling ambiguities are analyzed; finally, the trilinear alternating least squares (TALS) regression algorithm is used to fit the PARAFAC model; thus, the unknown parameters in the PARAFAC model can be estimated.


### 4.1. Parallel Factor Analysis

In this subsection, the basic idea of PARAFAC theory is introduced. First, the trilinear decomposition of the element of the three-way array is elaborated. Second, based on the symmetrical characteristic of the trilinear decomposition, the three-way array is decomposed into three dimensions.Considering aC
×
D
×
E three-way array X
_ (the entry of the array denotes x
c
,
d
,
e), the trilinear decomposition of x
c
,
d
,
e can be expressed as [21]
(15)
x
c
,
d
,
e
=
∑
p
=
1
M
‍
s
¯
c
,
p
t
¯
d
,
p
u
¯
e
,
p
,
where c
=
1
,
…
,
C, d
=
1
,
…
,
D, and e
=
1
,
…
,
E. The three-way array X
_ is represented as the sum of M rank-1 factor. Here, the rank of X
_ is defined as the minimum number of rank-1 three-way components which are needed to decompose X
_. The vectors s
¯
p
∈
𝒞
C
×
1, t
¯
p
∈
𝒞
D
×
1, and u
¯
p
∈
𝒞
E
×
1 are called load vector, score vector, and factor profiles, respectively.AssumeM
=
3; an example is given as follows and thus we can understand the idea of trilinear decomposition intuitively. The decomposition of three-way array X
_ is shown in Figure 4; it can be seen that three dimensions can be decomposed.Figure 4
The decomposition of three-way array.The typical element of theC
×
P matrix S
¯ is defined as S
¯
(
c
,
p
)
∶
=
s
¯
c
,
p; the typical element of the D
×
P matrix T
¯ is defined as T
¯
(
c
,
p
)
∶
=
t
¯
c
,
p; the typical element of the E
×
P matrix U
¯ is defined as U
¯
(
c
,
p
)
∶
=
u
¯
c
,
p. In addition, a D
×
E matrix X
c, a C
×
E matrix X
d, and a C
×
D matrix X
e are defined. The typical element of each matrix is X
c
(
d
,
e
)
∶
=
X
d
(
c
,
e
)
∶
=
X
e
(
c
,
d
)
∶
=
x
c
,
d
,
e. The three-way matrix X
_ can be “sliced” along three different dimensions as shown in Figure 4.Consider(16)
X
c
=
T
¯
Λ
c
(
S
¯
)
U
¯
T
,
c
=
1
,
…
,
C
,
X
d
=
U
¯
Λ
d
(
T
¯
)
S
¯
T
,
d
=
1
,
…
,
D
,
X
e
=
S
¯
Λ
e
(
U
¯
)
T
¯
T
,
e
=
1
,
…
,
E
,
where Λ
c
(
S
¯
) is the diagonal matrix constructed by the cth row of the matrix S
¯. According to the definition of Khatri-Rao product, X
_
(
C
D
×
E
) can be written in another form as
(17)
X
_
(
C
D
×
E
)
=
|
X
c
=
1
X
c
=
2
⋮
X
c
=
C
|
=
|
T
¯
Λ
1
(
S
¯
)
T
¯
Λ
2
(
S
¯
)
⋮
T
¯
Λ
C
(
S
¯
)
|
U
¯
T
=
(
S
¯
⊙
T
¯
)
U
¯
T
,
where “⊙” is used to denote the Khatri-Rao (KR) product. |
•
| stands for the three-way array which is constructed by stacking the matrices into a column. The definition of Khatri-Rao (KR) product can be found in the appendix. The three-way array X
_ can also be expressed in two other forms:
(18)
X
_
(
D
E
×
C
)
=
(
T
¯
⊙
U
¯
)
S
¯
T
,
(19)
X
_
(
E
C
×
D
)
=
(
U
¯
⊙
S
¯
)
T
¯
T
.The reason of decomposition along three different dimensions is that the PARAFAC model can be solved by using TALS algorithm. The matricesT
¯, U
¯, and S
¯ can be updated alternately.


### 4.2. The Identifiability of the Model

In this subsection, the PARAFAC model is constructed by the received data of the conformal array. The decomposition of the PARAFAC model must be unique; if not so, the PARAFAC model cannot be solved. Thus, some requirements of the matrices which construct the model have to be satisfied. Theorem2 guarantees the uniqueness decomposition of the model. If the PARAFAC model is constructed by the received data of the conformal array, then Theorem 2 can be used to judge the uniqueness decomposition of the model.On the basis of the PARAFAC theory [23], the m
×
m
×
4 three-way array of the cylindrical conformal array is constructed by (7)
(20)
|
R
(
:
,
:
,
1
)
R
(
:
,
:
,
2
)
R
(
:
,
:
,
3
)
R
(
:
,
:
,
4
)
|
=
|
R
1
R
2
R
3
R
4
|
=
|
B
R
s
B
H
B
Ψ
1
R
s
B
H
B
Ψ
2
R
s
B
H
B
Ψ
1
Ψ
2
R
s
B
H
|
+
Q
-
,
where R
i  (
i
=
1
,
…
,
4
) are the sample covariance matrices, respectively, and Q
- represents the noise in real observation. Let C
=
B
H; based on the definition of Khatri-Rao product, (20) can be transformed as
(21)
R
=
(
D
⊙
B
)
C
+
Q
-
,
(22)
R
X
=
(
C
T
⊙
D
)
B
T
+
Q
-
X
,
(23)
R
Y
=
(
B
⊙
C
T
)
D
T
+
Q
-
Y
,
(24)
D
=
[
Λ
-
1
(
R
s
)
Λ
-
1
(
Ψ
1
R
s
)
Λ
-
1
(
Ψ
2
R
s
)
Λ
-
1
(
Ψ
1
Ψ
2
R
s
)
]
,
where Λ
-
1
(
R
s
) is the row vector constructed by the diagonal entries of the diagonal matrix R
s. [
•
] stands for the matrix.Definition 1 (see [21]).
Considering a given matrixA
-
∈
𝒞
C
×
P, if and only if the matrix A
- contains at least r but not r
+
1 linearly independent columns, then the rank of matrix A
- is r
A
-
∶
=
Rank
(
A
-
)
=
r. If any k columns of matrix A
- are linearly independent, then the Kruskal rank (k-rank) is k
A
-
=
k [24]. Generally speaking, k
A
-
≤
r
A
-.In order to ensure the uniqueness of the PARAFAC model in (21)–(23), the sufficient condition must be met:
(25)
a
(
:
,
i
)
≠
a
(
:
,
j
)
.
For any two different incident signals, s
i
(
k
) and s
j
(
k
),
(26)
h
o
e
e
-
j
2
π
(
(
p
o
e
•
u
i
)
/
λ
)
≠
h
o
e
e
-
j
2
π
(
(
p
o
e
•
u
j
)
/
λ
)
ensures the uniqueness of the PARAFAC model, where
(27)
p
o
e
=
sin
(
θ
o
e
)
cos
⁡
(
φ
o
e
)
x
+
sin
(
θ
o
e
)
sin
(
φ
o
e
)
y
+
cos
⁡
(
θ
o
e
)
z
,
(28)
u
i
=
sin
(
θ
i
)
cos
⁡
(
φ
i
)
x
+
sin
(
θ
i
)
sin
(
φ
i
)
y
+
cos
⁡
(
θ
i
)
z
.
θ
o
e and φ
o
e are the elevation and azimuth of the sensor’s position in the global coordinate. θ
i and φ
i are the elevation and azimuth of the ith incident signal. Thus, (26) is equivalent to
(29)
sin
(
θ
o
e
)
cos
⁡
(
φ
o
e
)
ρ
1
+
sin
(
θ
o
e
)
sin
(
φ
o
e
)
ρ
2
+
cos
⁡
(
θ
o
e
)
ρ
3
≠
0
,
ρ
1
=
sin
(
θ
i
)
cos
⁡
(
φ
i
)
-
sin
(
θ
j
)
cos
⁡
(
φ
j
)
,
ρ
2
=
sin
(
θ
i
)
sin
(
φ
i
)
-
sin
(
θ
j
)
sin
(
φ
j
)
,
ρ
3
=
cos
⁡
(
θ
i
)
-
cos
⁡
(
θ
j
)
.Theorem 2.
The sufficient conditions of the PARAFAC model in (21)–(23) require that at least one of ρ
1 and ρ
2 and ρ
3 is not zero.Proof.
Thek-rank of matrices D, B, and C is k
D
=
min
⁡
(
4
,
r
), k
B
=
r, and k
C
T
=
r, respectively. The k-rank decomposition of the model is unique with k
D
+
k
B
+
k
C
T
≥
2
r
+
2 for k
≥
2 (see Theorem  1 in [21]).


### 4.3. The Trilinear Alternating Least Squares Regression Algorithm

The principle of the trilinear alternating least squares (TALS) regression algorithm is to fit the PARAFAC model in the noisy observation. The idea of TALS is very simple. In each step, only one matrix is updated. The updating algorithm is based on the estimated results of the last updating, and the remaining matrices are updated by least squares method. Repeat the step as mentioned above until the algorithm converges. The advantage of TALS algorithm is that the parameter adjustment is not necessary. A standard least square problem is solved in each step, and the performance of TALS is good [25]. Empirically, if the k-rank condition (Theorem 2) is satisfied, the TALS algorithm can converge to the global minimum [26]. The accelerating convergence technique of TALS algorithm can be found in [27, 28]. A brief introduction about TALS algorithm can be found in [21]. The application details of the TALS algorithm which is used to fit the PARAFAC model are elaborated as follows.On the basis of noisy observation, the problem (21) can be transformed into solving a least square problem
(30)
min
⁡
D
,
B
,
C
∥
R
-
(
D
⊙
B
)
C
∥
F
2
.
The principle of alternating least squares (ALS) can be used to fit the problem in (30). In the noiseless condition, ALS can be used to solve the matrices D, B, and C which constructed the three-way array R. The least square estimation of matrixCcan be expressed as
(31)
C
=
arg
min
⁡
C
∥
R
-
(
D
⊙
B
)
C
∥
F
2
.
Similarly, the matrices B and D can be expressed as
(32)
B
T
=
arg
min
⁡
B
∥
R
X
-
(
C
T
⊙
D
)
B
T
∥
F
2
,
D
T
=
arg
min
⁡
D
∥
R
Y
-
(
B
⊙
C
T
)
D
T
∥
F
2
.
In the iterative procedure, given matrices B and D, the matrix C can be represented as
(33)
C
=
(
D
⊙
B
)
†
R
.
The expression of matrices B
T and D
T is
(34)
B
T
=
(
C
T
⊙
D
)
†
R
X
,
D
T
=
(
B
⊙
C
T
)
†
R
Y
,
where (
•
)
† denotes the pseudoinverse of matrix (
•
).Now, the ALS algorithm steps can be summarized as follows:(1)
initializeB
(
0
)
∈
𝒞
M
×
P, D
(
0
)
∈
𝒞
4
×
P;
(2)
initializeε
>
0, k
=
0;
(3)
if∥
ρ
(
k
+
1
)
-
ρ
(
k
)
∥
/
ρ
(
k
)
>
ε, calculate matrices D, B, and C by (35)–(37). Update just one matrix at each time; then k
→
k
+
1;
(4)
else∥
ρ
(
k
+
1
)
-
ρ
(
k
)
∥
/
ρ
(
k
)
<
ε; the iteration is terminated.


### 4.4. The 2D-DOA Estimation Algorithm

The estimators of matricesD, B, and C are obtained by the TALS algorithm which is introduced in the last subsection. In this section, the 2D-DOA estimation can be obtained by the matrix D, which is shown as follows.The matrixD can be acquired by the TALS algorithm. ω
1
i and ω
2
i can be calculated by matrix D
(35)
ω
1
i
=
-
1
2
(
angle
⌈
D
2
i
D
1
i
⌉
+
angle
⌈
D
4
i
D
3
i
⌉
)
,
where D
j
i represents the jth row of D. ⌈
•
⌉stands for the absolute value operation. Because h
1 and h
3 are real numbers, D
3
i
/
D
1
i and D
4
i
/
D
2
i are squared to solve the ambiguity caused by the positive and negative values of h
1 and h
3.Consider(36)
ω
2
i
=
-
1
2
angle
(
[
h
3
(
θ
i
,
φ
i
)
h
1
(
θ
i
,
φ
i
)
exp
⁡
(
-
j
ω
2
i
)
]
2
)
=
-
1
2
angle
(
exp
⁡
(
-
j
2
ω
2
i
)
)
=
-
1
2
angle
(
[
D
3
i
D
1
i
]
2
)
=
-
1
2
angle
(
[
D
4
i
D
2
i
]
2
)
.
Then,
(37)
ω
2
i
=
-
1
4
|
angle
(
[
D
3
i
D
1
i
]
2
)
+
angle
(
[
D
4
i
D
2
i
]
2
)
|
.
Take (12), (13) into (10), (11), and the elevation θ
i and azimuth φ
i of ith incident signal can be expressed as
(38)
θ
i
=
arccos
(
λ
ω
1
i
2
π
d
2
)
=
arccos
(
2
ω
1
i
π
)
,
φ
i
=
arccos
(
λ
ω
2
i
2
π
d
2
sin
(
θ
i
)
)
=
arccos
(
2
ω
2
i
π
sin
(
θ
i
)
)
.Based on Theorem2, the estimators of D, B, and C have the same column permutation matrix; that is, the ith column of the steering matrix B corresponds to the ith of the matrix D. Thus, the elevation and azimuth pair with each other automatically.Combining PARAFAC theory with the TALS algorithm, the 2D-DOA estimation for the cylindrical conformal array can be summarized as follows:(1)
calculate the signals’ covariance matrices received by each subarray using (7);
(2)
construct the PARAFAC model by (20)–(23);
(3)
estimate the matrixD using the TALS algorithm;
(4)
Calculate (36) and (37) using the estimator of matrix D; then, obtain the estimators of ω
1
i and ω
2
i;
(5)
acquire the elevation and azimuth estimation from (40) and (41).


## 4.1. Parallel Factor Analysis

In this subsection, the basic idea of PARAFAC theory is introduced. First, the trilinear decomposition of the element of the three-way array is elaborated. Second, based on the symmetrical characteristic of the trilinear decomposition, the three-way array is decomposed into three dimensions.Considering aC
×
D
×
E three-way array X
_ (the entry of the array denotes x
c
,
d
,
e), the trilinear decomposition of x
c
,
d
,
e can be expressed as [21]
(15)
x
c
,
d
,
e
=
∑
p
=
1
M
‍
s
¯
c
,
p
t
¯
d
,
p
u
¯
e
,
p
,
where c
=
1
,
…
,
C, d
=
1
,
…
,
D, and e
=
1
,
…
,
E. The three-way array X
_ is represented as the sum of M rank-1 factor. Here, the rank of X
_ is defined as the minimum number of rank-1 three-way components which are needed to decompose X
_. The vectors s
¯
p
∈
𝒞
C
×
1, t
¯
p
∈
𝒞
D
×
1, and u
¯
p
∈
𝒞
E
×
1 are called load vector, score vector, and factor profiles, respectively.AssumeM
=
3; an example is given as follows and thus we can understand the idea of trilinear decomposition intuitively. The decomposition of three-way array X
_ is shown in Figure 4; it can be seen that three dimensions can be decomposed.Figure 4
The decomposition of three-way array.The typical element of theC
×
P matrix S
¯ is defined as S
¯
(
c
,
p
)
∶
=
s
¯
c
,
p; the typical element of the D
×
P matrix T
¯ is defined as T
¯
(
c
,
p
)
∶
=
t
¯
c
,
p; the typical element of the E
×
P matrix U
¯ is defined as U
¯
(
c
,
p
)
∶
=
u
¯
c
,
p. In addition, a D
×
E matrix X
c, a C
×
E matrix X
d, and a C
×
D matrix X
e are defined. The typical element of each matrix is X
c
(
d
,
e
)
∶
=
X
d
(
c
,
e
)
∶
=
X
e
(
c
,
d
)
∶
=
x
c
,
d
,
e. The three-way matrix X
_ can be “sliced” along three different dimensions as shown in Figure 4.Consider(16)
X
c
=
T
¯
Λ
c
(
S
¯
)
U
¯
T
,
c
=
1
,
…
,
C
,
X
d
=
U
¯
Λ
d
(
T
¯
)
S
¯
T
,
d
=
1
,
…
,
D
,
X
e
=
S
¯
Λ
e
(
U
¯
)
T
¯
T
,
e
=
1
,
…
,
E
,
where Λ
c
(
S
¯
) is the diagonal matrix constructed by the cth row of the matrix S
¯. According to the definition of Khatri-Rao product, X
_
(
C
D
×
E
) can be written in another form as
(17)
X
_
(
C
D
×
E
)
=
|
X
c
=
1
X
c
=
2
⋮
X
c
=
C
|
=
|
T
¯
Λ
1
(
S
¯
)
T
¯
Λ
2
(
S
¯
)
⋮
T
¯
Λ
C
(
S
¯
)
|
U
¯
T
=
(
S
¯
⊙
T
¯
)
U
¯
T
,
where “⊙” is used to denote the Khatri-Rao (KR) product. |
•
| stands for the three-way array which is constructed by stacking the matrices into a column. The definition of Khatri-Rao (KR) product can be found in the appendix. The three-way array X
_ can also be expressed in two other forms:
(18)
X
_
(
D
E
×
C
)
=
(
T
¯
⊙
U
¯
)
S
¯
T
,
(19)
X
_
(
E
C
×
D
)
=
(
U
¯
⊙
S
¯
)
T
¯
T
.The reason of decomposition along three different dimensions is that the PARAFAC model can be solved by using TALS algorithm. The matricesT
¯, U
¯, and S
¯ can be updated alternately.


## 4.2. The Identifiability of the Model

In this subsection, the PARAFAC model is constructed by the received data of the conformal array. The decomposition of the PARAFAC model must be unique; if not so, the PARAFAC model cannot be solved. Thus, some requirements of the matrices which construct the model have to be satisfied. Theorem2 guarantees the uniqueness decomposition of the model. If the PARAFAC model is constructed by the received data of the conformal array, then Theorem 2 can be used to judge the uniqueness decomposition of the model.On the basis of the PARAFAC theory [23], the m
×
m
×
4 three-way array of the cylindrical conformal array is constructed by (7)
(20)
|
R
(
:
,
:
,
1
)
R
(
:
,
:
,
2
)
R
(
:
,
:
,
3
)
R
(
:
,
:
,
4
)
|
=
|
R
1
R
2
R
3
R
4
|
=
|
B
R
s
B
H
B
Ψ
1
R
s
B
H
B
Ψ
2
R
s
B
H
B
Ψ
1
Ψ
2
R
s
B
H
|
+
Q
-
,
where R
i  (
i
=
1
,
…
,
4
) are the sample covariance matrices, respectively, and Q
- represents the noise in real observation. Let C
=
B
H; based on the definition of Khatri-Rao product, (20) can be transformed as
(21)
R
=
(
D
⊙
B
)
C
+
Q
-
,
(22)
R
X
=
(
C
T
⊙
D
)
B
T
+
Q
-
X
,
(23)
R
Y
=
(
B
⊙
C
T
)
D
T
+
Q
-
Y
,
(24)
D
=
[
Λ
-
1
(
R
s
)
Λ
-
1
(
Ψ
1
R
s
)
Λ
-
1
(
Ψ
2
R
s
)
Λ
-
1
(
Ψ
1
Ψ
2
R
s
)
]
,
where Λ
-
1
(
R
s
) is the row vector constructed by the diagonal entries of the diagonal matrix R
s. [
•
] stands for the matrix.Definition 1 (see [21]).
Considering a given matrixA
-
∈
𝒞
C
×
P, if and only if the matrix A
- contains at least r but not r
+
1 linearly independent columns, then the rank of matrix A
- is r
A
-
∶
=
Rank
(
A
-
)
=
r. If any k columns of matrix A
- are linearly independent, then the Kruskal rank (k-rank) is k
A
-
=
k [24]. Generally speaking, k
A
-
≤
r
A
-.In order to ensure the uniqueness of the PARAFAC model in (21)–(23), the sufficient condition must be met:
(25)
a
(
:
,
i
)
≠
a
(
:
,
j
)
.
For any two different incident signals, s
i
(
k
) and s
j
(
k
),
(26)
h
o
e
e
-
j
2
π
(
(
p
o
e
•
u
i
)
/
λ
)
≠
h
o
e
e
-
j
2
π
(
(
p
o
e
•
u
j
)
/
λ
)
ensures the uniqueness of the PARAFAC model, where
(27)
p
o
e
=
sin
(
θ
o
e
)
cos
⁡
(
φ
o
e
)
x
+
sin
(
θ
o
e
)
sin
(
φ
o
e
)
y
+
cos
⁡
(
θ
o
e
)
z
,
(28)
u
i
=
sin
(
θ
i
)
cos
⁡
(
φ
i
)
x
+
sin
(
θ
i
)
sin
(
φ
i
)
y
+
cos
⁡
(
θ
i
)
z
.
θ
o
e and φ
o
e are the elevation and azimuth of the sensor’s position in the global coordinate. θ
i and φ
i are the elevation and azimuth of the ith incident signal. Thus, (26) is equivalent to
(29)
sin
(
θ
o
e
)
cos
⁡
(
φ
o
e
)
ρ
1
+
sin
(
θ
o
e
)
sin
(
φ
o
e
)
ρ
2
+
cos
⁡
(
θ
o
e
)
ρ
3
≠
0
,
ρ
1
=
sin
(
θ
i
)
cos
⁡
(
φ
i
)
-
sin
(
θ
j
)
cos
⁡
(
φ
j
)
,
ρ
2
=
sin
(
θ
i
)
sin
(
φ
i
)
-
sin
(
θ
j
)
sin
(
φ
j
)
,
ρ
3
=
cos
⁡
(
θ
i
)
-
cos
⁡
(
θ
j
)
.Theorem 2.
The sufficient conditions of the PARAFAC model in (21)–(23) require that at least one of ρ
1 and ρ
2 and ρ
3 is not zero.Proof.
Thek-rank of matrices D, B, and C is k
D
=
min
⁡
(
4
,
r
), k
B
=
r, and k
C
T
=
r, respectively. The k-rank decomposition of the model is unique with k
D
+
k
B
+
k
C
T
≥
2
r
+
2 for k
≥
2 (see Theorem  1 in [21]).


## 4.3. The Trilinear Alternating Least Squares Regression Algorithm

The principle of the trilinear alternating least squares (TALS) regression algorithm is to fit the PARAFAC model in the noisy observation. The idea of TALS is very simple. In each step, only one matrix is updated. The updating algorithm is based on the estimated results of the last updating, and the remaining matrices are updated by least squares method. Repeat the step as mentioned above until the algorithm converges. The advantage of TALS algorithm is that the parameter adjustment is not necessary. A standard least square problem is solved in each step, and the performance of TALS is good [25]. Empirically, if the k-rank condition (Theorem 2) is satisfied, the TALS algorithm can converge to the global minimum [26]. The accelerating convergence technique of TALS algorithm can be found in [27, 28]. A brief introduction about TALS algorithm can be found in [21]. The application details of the TALS algorithm which is used to fit the PARAFAC model are elaborated as follows.On the basis of noisy observation, the problem (21) can be transformed into solving a least square problem
(30)
min
⁡
D
,
B
,
C
∥
R
-
(
D
⊙
B
)
C
∥
F
2
.
The principle of alternating least squares (ALS) can be used to fit the problem in (30). In the noiseless condition, ALS can be used to solve the matrices D, B, and C which constructed the three-way array R. The least square estimation of matrixCcan be expressed as
(31)
C
=
arg
min
⁡
C
∥
R
-
(
D
⊙
B
)
C
∥
F
2
.
Similarly, the matrices B and D can be expressed as
(32)
B
T
=
arg
min
⁡
B
∥
R
X
-
(
C
T
⊙
D
)
B
T
∥
F
2
,
D
T
=
arg
min
⁡
D
∥
R
Y
-
(
B
⊙
C
T
)
D
T
∥
F
2
.
In the iterative procedure, given matrices B and D, the matrix C can be represented as
(33)
C
=
(
D
⊙
B
)
†
R
.
The expression of matrices B
T and D
T is
(34)
B
T
=
(
C
T
⊙
D
)
†
R
X
,
D
T
=
(
B
⊙
C
T
)
†
R
Y
,
where (
•
)
† denotes the pseudoinverse of matrix (
•
).Now, the ALS algorithm steps can be summarized as follows:(1)
initializeB
(
0
)
∈
𝒞
M
×
P, D
(
0
)
∈
𝒞
4
×
P;
(2)
initializeε
>
0, k
=
0;
(3)
if∥
ρ
(
k
+
1
)
-
ρ
(
k
)
∥
/
ρ
(
k
)
>
ε, calculate matrices D, B, and C by (35)–(37). Update just one matrix at each time; then k
→
k
+
1;
(4)
else∥
ρ
(
k
+
1
)
-
ρ
(
k
)
∥
/
ρ
(
k
)
<
ε; the iteration is terminated.


## 4.4. The 2D-DOA Estimation Algorithm

The estimators of matricesD, B, and C are obtained by the TALS algorithm which is introduced in the last subsection. In this section, the 2D-DOA estimation can be obtained by the matrix D, which is shown as follows.The matrixD can be acquired by the TALS algorithm. ω
1
i and ω
2
i can be calculated by matrix D
(35)
ω
1
i
=
-
1
2
(
angle
⌈
D
2
i
D
1
i
⌉
+
angle
⌈
D
4
i
D
3
i
⌉
)
,
where D
j
i represents the jth row of D. ⌈
•
⌉stands for the absolute value operation. Because h
1 and h
3 are real numbers, D
3
i
/
D
1
i and D
4
i
/
D
2
i are squared to solve the ambiguity caused by the positive and negative values of h
1 and h
3.Consider(36)
ω
2
i
=
-
1
2
angle
(
[
h
3
(
θ
i
,
φ
i
)
h
1
(
θ
i
,
φ
i
)
exp
⁡
(
-
j
ω
2
i
)
]
2
)
=
-
1
2
angle
(
exp
⁡
(
-
j
2
ω
2
i
)
)
=
-
1
2
angle
(
[
D
3
i
D
1
i
]
2
)
=
-
1
2
angle
(
[
D
4
i
D
2
i
]
2
)
.
Then,
(37)
ω
2
i
=
-
1
4
|
angle
(
[
D
3
i
D
1
i
]
2
)
+
angle
(
[
D
4
i
D
2
i
]
2
)
|
.
Take (12), (13) into (10), (11), and the elevation θ
i and azimuth φ
i of ith incident signal can be expressed as
(38)
θ
i
=
arccos
(
λ
ω
1
i
2
π
d
2
)
=
arccos
(
2
ω
1
i
π
)
,
φ
i
=
arccos
(
λ
ω
2
i
2
π
d
2
sin
(
θ
i
)
)
=
arccos
(
2
ω
2
i
π
sin
(
θ
i
)
)
.Based on Theorem2, the estimators of D, B, and C have the same column permutation matrix; that is, the ith column of the steering matrix B corresponds to the ith of the matrix D. Thus, the elevation and azimuth pair with each other automatically.Combining PARAFAC theory with the TALS algorithm, the 2D-DOA estimation for the cylindrical conformal array can be summarized as follows:(1)
calculate the signals’ covariance matrices received by each subarray using (7);
(2)
construct the PARAFAC model by (20)–(23);
(3)
estimate the matrixD using the TALS algorithm;
(4)
Calculate (36) and (37) using the estimator of matrix D; then, obtain the estimators of ω
1
i and ω
2
i;
(5)
acquire the elevation and azimuth estimation from (40) and (41).


## 5. The Extendable Array Structure

The proposed algorithm can be extended to other array structures with little modification. First, some elements are added in the cylindrical conformal array, and then the elements can be arranged more flexibly. Second, the proposed algorithm is extended to conical conformal array.


### 5.1. The Extendable Cylindrical Conformal Array

First, the design of cylindrical conformal array is introduced. Second, the PARAFAC model is constructed by the received data. Finally, the 2D-DOA estimation is obtained.As shown in Figures2 and 3, the proposed algorithm is feasible when the distance vector Δ
P
2 between array 1 and array 3 is parallel to X-axis. However, the proposed algorithm must be modified to be used in general cases. Arrays 5 and array 6 are added so that the arrays can be designed more flexibly. The extendable cylindrical conformal array is shown in Figure 5. 2
m
+
3
~
3
m
+
2 is array 5, and 2
m
+
4
~
3
m
+
3 is array 6. The sensors 2
m
+
3
~
3
m
+
3 possess the same pattern g
3. The distance vector between array 1 and array 5 is Δ
P
3, and d
3
=
|
Δ
P
3
|. Under the design in Figure 5, the only restriction is that arrays are parallel to Z-axis.Figure 5
The extendable cylindrical conformal array structure.The received data of array 5 and array 6 are(39)
X
5
=
B
Ψ
3
S
+
N
5
,
(40)
X
6
=
B
Ψ
1
Ψ
3
S
+
N
6
,
where
(41)
Ψ
3
=
diag
⁡
[
h
5
(
θ
1
,
φ
1
)
h
1
(
θ
1
,
φ
1
)
exp
⁡
(
-
j
ω
31
)
,
m
m
h
5
(
θ
2
,
φ
2
)
h
1
(
θ
2
,
φ
2
)
exp
⁡
(
-
j
ω
32
)
,
…
,
m
m
h
5
(
θ
r
,
φ
r
)
h
1
(
θ
r
,
φ
r
)
exp
⁡
(
-
j
ω
3
r
)
]
,
(42)
ω
3
i
=
(
2
π
λ
)
Δ
P
3
•
u
i
=
(
2
π
d
2
λ
)
×
[
sin
(
θ
Δ
P
3
)
cos
⁡
(
φ
Δ
P
3
)
sin
(
θ
i
)
cos
⁡
(
φ
i
)
+
sin
(
θ
Δ
P
2
)
sin
(
φ
Δ
P
2
)
sin
(
θ
i
)
sin
(
φ
i
)
+
cos
⁡
(
θ
Δ
P
3
)
cos
⁡
(
θ
i
)
]
.
The cross-covariance matrices among array 1, array 5, and array 6 are
(43)
R
5
=
E
{
X
5
X
1
H
}
=
B
Ψ
3
R
s
B
H
+
Q
5
,
R
6
=
E
{
X
6
X
1
H
}
=
B
Ψ
1
Ψ
3
R
s
B
H
+
Q
6
,
whereQ
5andQ
6are noise covariance matrices. Then, (20) is modified into anm
×
m
×
6three-way array PARAFAC model
(44)
|
R
(
:
,
:
,
1
)
R
(
:
,
:
,
2
)
R
(
:
,
:
,
3
)
R
(
:
,
:
,
4
)
R
(
:
,
:
,
5
)
R
(
:
,
:
,
6
)
|
=
|
R
1
R
2
R
3
R
4
R
5
R
6
|
=
|
B
R
s
B
H
B
Ψ
1
R
s
B
H
B
Ψ
2
R
s
B
H
B
Ψ
1
Ψ
2
R
s
B
H
B
Ψ
3
R
s
B
H
B
Ψ
1
Ψ
3
R
s
B
H
|
+
Q
-
1
,
where Q
-
1 is the observed noise. The form of Khatri-Rao product is applied in (44)
(45)
R
-
=
(
D
-
⊙
B
)
C
+
Q
-
1
,
where
(46)
D
-
=
[
Λ
-
1
(
R
s
)
Λ
-
1
(
Ψ
1
R
s
)
Λ
-
1
(
Ψ
2
R
s
)
Λ
-
1
(
Ψ
1
Ψ
2
R
s
)
Λ
-
1
(
Ψ
3
R
s
)
Λ
-
1
(
Ψ
1
Ψ
3
R
s
)
]
.As long as Theorem2 holds, (45) is unique. ω
3
i can be obtained similarly with (37)
(47)
ω
3
i
=
-
1
4
|
angle
(
[
D
5
i
D
1
i
]
2
)
+
angle
(
[
D
6
i
D
2
i
]
2
)
|
after calculating matrix D by ALS algorithm.Assume thatΔ
p
11
=
sin
(
θ
Δ
P
1
)
cos
⁡
(
φ
Δ
P
1
), Δ
p
12
=
sin
(
θ
Δ
P
1
)
sin
(
φ
Δ
P
1
), and Δ
p
13
=
cos
⁡
(
θ
Δ
P
1
); similar to Δ
p
2
i and Δ
p
3
i  (
i
=
1,2
,
3
), γ
1
i
=
sin
(
θ
i
)
cos
⁡
(
φ
i
), γ
2
i
=
sin
(
θ
i
)
sin
(
φ
i
), and γ
3
i
=
cos
⁡
(
θ
i
); solving (10), (11), (42) (36), (37), and (24), we can obtain
(48)
λ
2
π
[
ω
1
i
d
1
ω
2
i
d
2
ω
3
i
d
3
]
=
[
Δ
p
11
Δ
p
12
Δ
p
13
Δ
p
21
Δ
p
22
Δ
p
23
Δ
p
31
Δ
p
32
Δ
p
33
]
[
γ
1
i
γ
2
i
γ
3
i
]
.
The solution of (48) is
(49)
[
γ
1
i
γ
2
i
γ
3
i
]
=
λ
2
π
[
Δ
p
11
Δ
p
12
Δ
p
13
Δ
p
21
Δ
p
22
Δ
p
23
Δ
p
31
Δ
p
32
Δ
p
33
]
-
1
[
ω
1
i
d
1
ω
2
i
d
2
ω
3
i
d
3
]
.
γ
1, γ
2, and γ
3 can be acquired by solving (49). Two methods can be used to obtainθ
i, φ
i. The first method can be described as
(50)
θ
i
=
arccos
γ
3
i
,
φ
i
=
arcsin
(
γ
1
i
cos
⁡
(
θ
i
)
)
or
φ
i
=
arcsin
(
γ
2
i
sin
(
θ
i
)
)
.
The second one is
(51)
φ
i
=
arctan
(
γ
2
i
γ
1
i
)
,
θ
i
=
arccos
(
γ
1
i
cos
⁡
(
φ
i
)
)
or
θ
i
=
arcsin
(
γ
2
i
sin
(
φ
i
)
)
.


### 5.2. The Extendable Conical Conformal Array

In this section, the proposed algorithm is extended to conical conformal array. First, the design of conical conformal array is introduced. Second, the PARAFAC model is constructed by the received data. Finally, the 2D-DOA estimation is obtained. The analytic solution of the conical conformal array does not exist. Thus, the iterative numerical approximation method is used to obtain the optimal solution.Figure6 applies the proposed algorithm to the conical conformal array. 1
~
m is array 1, 2
~
m
+
1 is array 2, 1
~
2
m is array 3, and m
+
2
~
2
m
+
1 is array 4. The distance vector between array 1 and array 2 is Δ
P
1
¯
¯, and d
1
¯
¯
=
|
Δ
P
1
¯
¯
|. The distance vector between array 3 and array 4 is Δ
P
2
¯
¯, and d
2
¯
¯
=
|
Δ
P
2
¯
¯
|. 1
~
m
+
1 possess the same pattern g
1, and 1
~
2
m
+
1 possess the same pattern g
2. The decoupling between polarization parameter and DOA information can be completed by making full use of the geometry characteristic of the array structure.Figure 6
The extendable conical conformal array structure.The received data from array 1 to array 4 are(52)
X
1
=
B
1
S
+
N
1
,
(53)
X
2
=
B
1
Ψ
¯
1
S
+
N
2
,
(54)
X
3
=
B
2
S
+
N
3
,
(55)
X
4
=
B
2
Ψ
¯
2
S
+
N
4
,
where N
1
–
N
6 are the noise matrices received by each array. The cross-covariance matrices among the arrays are
(56)
R
1
=
E
{
X
3
X
1
H
}
=
B
2
R
s
B
1
H
+
Q
¯
¯
1
,
R
2
=
E
{
X
4
X
1
H
}
=
B
2
Ψ
¯
2
R
s
B
1
H
+
Q
¯
¯
2
,
R
3
=
E
{
X
3
X
2
H
}
=
B
2
Ψ
¯
1
H
R
s
B
1
H
+
Q
¯
¯
3
,
R
4
=
E
{
X
4
X
2
H
}
=
B
2
Ψ
¯
1
H
Ψ
¯
2
R
s
B
1
H
+
Q
¯
¯
4
,
where B
1 is the manifold matrix of array 1 and array 2 and B
2 is the manifold matrix of array 3 and array 4. Q
¯
¯
1
–
Q
¯
¯
4 are the noise covariance matrices received by the arrays. Ψ
¯
1 and Ψ
¯
2 possess the same forms as (8) and (9). According to PARAFAC theory, the m
×
m
×
4 three-way array for conical conformal array can be constructed by (56)
(57)
|
R
(
:
,
:
,
1
)
R
(
:
,
:
,
2
)
R
(
:
,
:
,
3
)
R
(
:
,
:
,
4
)
|
=
|
R
1
R
2
R
3
R
4
|
=
|
B
2
R
s
B
1
H
B
2
Ψ
¯
1
H
R
s
B
1
H
B
2
Ψ
¯
2
R
s
B
1
H
B
2
Ψ
¯
1
H
Ψ
¯
2
R
s
B
1
H
|
+
Q
-
2
and Q
-
2 is the observed noise. Let C
¯
¯
=
B
1
H, and (57) can be written in the following form of Khatri-Rao product:
(58)
R
¯
¯
=
(
D
¯
¯
⊙
B
2
)
C
¯
¯
+
Q
-
2
,
D
¯
¯
=
[
Λ
-
1
(
R
s
)
Λ
-
1
(
Ψ
¯
1
H
R
s
)
Λ
-
1
(
Ψ
¯
2
R
s
)
Λ
-
1
(
Ψ
¯
1
H
Ψ
¯
2
R
s
)
]
.
The matrix D can be obtained by the ALS algorithm if the uniqueness of (57) is guaranteed. By using the similar method, ω
¯
¯
1
i and ω
¯
¯
2
i can be represented as
(59)
ω
¯
¯
1
i
=
1
2
(
angle
⌈
D
¯
¯
2
i
D
¯
¯
1
i
⌉
+
angle
⌈
D
¯
¯
4
i
D
¯
¯
3
i
⌉
)
,
(60)
ω
¯
¯
2
i
=
-
1
2
(
angle
⌈
D
¯
¯
3
i
D
¯
¯
1
i
⌉
+
angle
⌈
D
¯
¯
4
i
D
¯
¯
2
i
⌉
)
,
(61)
ω
¯
¯
1
i
=
(
2
π
λ
i
)
Δ
P
1
¯
¯
•
u
i
=
(
2
π
d
1
¯
¯
λ
i
)
×
[
sin
(
θ
Δ
P
1
¯
¯
)
cos
⁡
(
φ
Δ
P
1
¯
¯
)
sin
(
θ
i
)
cos
⁡
(
φ
i
)
+
sin
(
θ
Δ
P
1
¯
¯
)
sin
(
φ
Δ
P
1
¯
¯
)
sin
(
θ
i
)
sin
(
φ
i
)
+
cos
⁡
(
θ
Δ
P
1
¯
¯
)
cos
⁡
(
θ
i
)
]
,
(62)
ω
¯
¯
2
i
=
(
2
π
λ
i
)
Δ
P
2
¯
¯
•
u
i
=
(
2
π
d
2
¯
¯
λ
i
)
×
[
sin
(
θ
Δ
P
2
¯
¯
)
cos
⁡
(
φ
Δ
P
2
¯
¯
)
sin
(
θ
i
)
cos
⁡
(
φ
i
)
m
+
sin
(
θ
Δ
P
2
¯
¯
)
sin
(
φ
Δ
P
2
¯
¯
)
sin
(
θ
i
)
sin
(
φ
i
)
m
+
cos
⁡
(
θ
Δ
P
2
¯
¯
)
cos
⁡
(
θ
i
)
]
.
The 2D-DOA estimation can be acquired by solving (61) and (62). The equations are nonlinear and the analytical solution does not exist. Thus, the iterative numerical approximation method is used to obtain the optimal solution. The nonlinear equations can be written as
(63)
f
1
(
x
1
,
x
2
)
=
0
,
f
2
(
x
1
,
x
2
)
=
0
,
where x
1
=
θ and x
2
=
φ. The equation with the same solution of (63) is expressed as
(64)
x
=
f
i
-
1
(
x
1
,
x
2
)
,
i
=
1,2
.
The iterative form can be represented as
(65)
x
k
+
1
=
F
i
(
x
1
k
,
x
2
k
)
,
i
=
1,2
.
The initial vector is selected as x
(
0
)
=
[
x
1
(
0
)
,
x
2
(
0
)
]
T. Calculating (65) until the sequence converges, that is, x
k
→
x
*, and regarding x
* as the iterative solution of (63), we can achieve the 2D-DOA estimation of the incident signal.Assuming that the noise is nonuniform, in theory through calculating the cross-covariance matrices among different arrays, nonuniform noise could be suppressed.


## 5.1. The Extendable Cylindrical Conformal Array

First, the design of cylindrical conformal array is introduced. Second, the PARAFAC model is constructed by the received data. Finally, the 2D-DOA estimation is obtained.As shown in Figures2 and 3, the proposed algorithm is feasible when the distance vector Δ
P
2 between array 1 and array 3 is parallel to X-axis. However, the proposed algorithm must be modified to be used in general cases. Arrays 5 and array 6 are added so that the arrays can be designed more flexibly. The extendable cylindrical conformal array is shown in Figure 5. 2
m
+
3
~
3
m
+
2 is array 5, and 2
m
+
4
~
3
m
+
3 is array 6. The sensors 2
m
+
3
~
3
m
+
3 possess the same pattern g
3. The distance vector between array 1 and array 5 is Δ
P
3, and d
3
=
|
Δ
P
3
|. Under the design in Figure 5, the only restriction is that arrays are parallel to Z-axis.Figure 5
The extendable cylindrical conformal array structure.The received data of array 5 and array 6 are(39)
X
5
=
B
Ψ
3
S
+
N
5
,
(40)
X
6
=
B
Ψ
1
Ψ
3
S
+
N
6
,
where
(41)
Ψ
3
=
diag
⁡
[
h
5
(
θ
1
,
φ
1
)
h
1
(
θ
1
,
φ
1
)
exp
⁡
(
-
j
ω
31
)
,
m
m
h
5
(
θ
2
,
φ
2
)
h
1
(
θ
2
,
φ
2
)
exp
⁡
(
-
j
ω
32
)
,
…
,
m
m
h
5
(
θ
r
,
φ
r
)
h
1
(
θ
r
,
φ
r
)
exp
⁡
(
-
j
ω
3
r
)
]
,
(42)
ω
3
i
=
(
2
π
λ
)
Δ
P
3
•
u
i
=
(
2
π
d
2
λ
)
×
[
sin
(
θ
Δ
P
3
)
cos
⁡
(
φ
Δ
P
3
)
sin
(
θ
i
)
cos
⁡
(
φ
i
)
+
sin
(
θ
Δ
P
2
)
sin
(
φ
Δ
P
2
)
sin
(
θ
i
)
sin
(
φ
i
)
+
cos
⁡
(
θ
Δ
P
3
)
cos
⁡
(
θ
i
)
]
.
The cross-covariance matrices among array 1, array 5, and array 6 are
(43)
R
5
=
E
{
X
5
X
1
H
}
=
B
Ψ
3
R
s
B
H
+
Q
5
,
R
6
=
E
{
X
6
X
1
H
}
=
B
Ψ
1
Ψ
3
R
s
B
H
+
Q
6
,
whereQ
5andQ
6are noise covariance matrices. Then, (20) is modified into anm
×
m
×
6three-way array PARAFAC model
(44)
|
R
(
:
,
:
,
1
)
R
(
:
,
:
,
2
)
R
(
:
,
:
,
3
)
R
(
:
,
:
,
4
)
R
(
:
,
:
,
5
)
R
(
:
,
:
,
6
)
|
=
|
R
1
R
2
R
3
R
4
R
5
R
6
|
=
|
B
R
s
B
H
B
Ψ
1
R
s
B
H
B
Ψ
2
R
s
B
H
B
Ψ
1
Ψ
2
R
s
B
H
B
Ψ
3
R
s
B
H
B
Ψ
1
Ψ
3
R
s
B
H
|
+
Q
-
1
,
where Q
-
1 is the observed noise. The form of Khatri-Rao product is applied in (44)
(45)
R
-
=
(
D
-
⊙
B
)
C
+
Q
-
1
,
where
(46)
D
-
=
[
Λ
-
1
(
R
s
)
Λ
-
1
(
Ψ
1
R
s
)
Λ
-
1
(
Ψ
2
R
s
)
Λ
-
1
(
Ψ
1
Ψ
2
R
s
)
Λ
-
1
(
Ψ
3
R
s
)
Λ
-
1
(
Ψ
1
Ψ
3
R
s
)
]
.As long as Theorem2 holds, (45) is unique. ω
3
i can be obtained similarly with (37)
(47)
ω
3
i
=
-
1
4
|
angle
(
[
D
5
i
D
1
i
]
2
)
+
angle
(
[
D
6
i
D
2
i
]
2
)
|
after calculating matrix D by ALS algorithm.Assume thatΔ
p
11
=
sin
(
θ
Δ
P
1
)
cos
⁡
(
φ
Δ
P
1
), Δ
p
12
=
sin
(
θ
Δ
P
1
)
sin
(
φ
Δ
P
1
), and Δ
p
13
=
cos
⁡
(
θ
Δ
P
1
); similar to Δ
p
2
i and Δ
p
3
i  (
i
=
1,2
,
3
), γ
1
i
=
sin
(
θ
i
)
cos
⁡
(
φ
i
), γ
2
i
=
sin
(
θ
i
)
sin
(
φ
i
), and γ
3
i
=
cos
⁡
(
θ
i
); solving (10), (11), (42) (36), (37), and (24), we can obtain
(48)
λ
2
π
[
ω
1
i
d
1
ω
2
i
d
2
ω
3
i
d
3
]
=
[
Δ
p
11
Δ
p
12
Δ
p
13
Δ
p
21
Δ
p
22
Δ
p
23
Δ
p
31
Δ
p
32
Δ
p
33
]
[
γ
1
i
γ
2
i
γ
3
i
]
.
The solution of (48) is
(49)
[
γ
1
i
γ
2
i
γ
3
i
]
=
λ
2
π
[
Δ
p
11
Δ
p
12
Δ
p
13
Δ
p
21
Δ
p
22
Δ
p
23
Δ
p
31
Δ
p
32
Δ
p
33
]
-
1
[
ω
1
i
d
1
ω
2
i
d
2
ω
3
i
d
3
]
.
γ
1, γ
2, and γ
3 can be acquired by solving (49). Two methods can be used to obtainθ
i, φ
i. The first method can be described as
(50)
θ
i
=
arccos
γ
3
i
,
φ
i
=
arcsin
(
γ
1
i
cos
⁡
(
θ
i
)
)
or
φ
i
=
arcsin
(
γ
2
i
sin
(
θ
i
)
)
.
The second one is
(51)
φ
i
=
arctan
(
γ
2
i
γ
1
i
)
,
θ
i
=
arccos
(
γ
1
i
cos
⁡
(
φ
i
)
)
or
θ
i
=
arcsin
(
γ
2
i
sin
(
φ
i
)
)
.


## 5.2. The Extendable Conical Conformal Array

In this section, the proposed algorithm is extended to conical conformal array. First, the design of conical conformal array is introduced. Second, the PARAFAC model is constructed by the received data. Finally, the 2D-DOA estimation is obtained. The analytic solution of the conical conformal array does not exist. Thus, the iterative numerical approximation method is used to obtain the optimal solution.Figure6 applies the proposed algorithm to the conical conformal array. 1
~
m is array 1, 2
~
m
+
1 is array 2, 1
~
2
m is array 3, and m
+
2
~
2
m
+
1 is array 4. The distance vector between array 1 and array 2 is Δ
P
1
¯
¯, and d
1
¯
¯
=
|
Δ
P
1
¯
¯
|. The distance vector between array 3 and array 4 is Δ
P
2
¯
¯, and d
2
¯
¯
=
|
Δ
P
2
¯
¯
|. 1
~
m
+
1 possess the same pattern g
1, and 1
~
2
m
+
1 possess the same pattern g
2. The decoupling between polarization parameter and DOA information can be completed by making full use of the geometry characteristic of the array structure.Figure 6
The extendable conical conformal array structure.The received data from array 1 to array 4 are(52)
X
1
=
B
1
S
+
N
1
,
(53)
X
2
=
B
1
Ψ
¯
1
S
+
N
2
,
(54)
X
3
=
B
2
S
+
N
3
,
(55)
X
4
=
B
2
Ψ
¯
2
S
+
N
4
,
where N
1
–
N
6 are the noise matrices received by each array. The cross-covariance matrices among the arrays are
(56)
R
1
=
E
{
X
3
X
1
H
}
=
B
2
R
s
B
1
H
+
Q
¯
¯
1
,
R
2
=
E
{
X
4
X
1
H
}
=
B
2
Ψ
¯
2
R
s
B
1
H
+
Q
¯
¯
2
,
R
3
=
E
{
X
3
X
2
H
}
=
B
2
Ψ
¯
1
H
R
s
B
1
H
+
Q
¯
¯
3
,
R
4
=
E
{
X
4
X
2
H
}
=
B
2
Ψ
¯
1
H
Ψ
¯
2
R
s
B
1
H
+
Q
¯
¯
4
,
where B
1 is the manifold matrix of array 1 and array 2 and B
2 is the manifold matrix of array 3 and array 4. Q
¯
¯
1
–
Q
¯
¯
4 are the noise covariance matrices received by the arrays. Ψ
¯
1 and Ψ
¯
2 possess the same forms as (8) and (9). According to PARAFAC theory, the m
×
m
×
4 three-way array for conical conformal array can be constructed by (56)
(57)
|
R
(
:
,
:
,
1
)
R
(
:
,
:
,
2
)
R
(
:
,
:
,
3
)
R
(
:
,
:
,
4
)
|
=
|
R
1
R
2
R
3
R
4
|
=
|
B
2
R
s
B
1
H
B
2
Ψ
¯
1
H
R
s
B
1
H
B
2
Ψ
¯
2
R
s
B
1
H
B
2
Ψ
¯
1
H
Ψ
¯
2
R
s
B
1
H
|
+
Q
-
2
and Q
-
2 is the observed noise. Let C
¯
¯
=
B
1
H, and (57) can be written in the following form of Khatri-Rao product:
(58)
R
¯
¯
=
(
D
¯
¯
⊙
B
2
)
C
¯
¯
+
Q
-
2
,
D
¯
¯
=
[
Λ
-
1
(
R
s
)
Λ
-
1
(
Ψ
¯
1
H
R
s
)
Λ
-
1
(
Ψ
¯
2
R
s
)
Λ
-
1
(
Ψ
¯
1
H
Ψ
¯
2
R
s
)
]
.
The matrix D can be obtained by the ALS algorithm if the uniqueness of (57) is guaranteed. By using the similar method, ω
¯
¯
1
i and ω
¯
¯
2
i can be represented as
(59)
ω
¯
¯
1
i
=
1
2
(
angle
⌈
D
¯
¯
2
i
D
¯
¯
1
i
⌉
+
angle
⌈
D
¯
¯
4
i
D
¯
¯
3
i
⌉
)
,
(60)
ω
¯
¯
2
i
=
-
1
2
(
angle
⌈
D
¯
¯
3
i
D
¯
¯
1
i
⌉
+
angle
⌈
D
¯
¯
4
i
D
¯
¯
2
i
⌉
)
,
(61)
ω
¯
¯
1
i
=
(
2
π
λ
i
)
Δ
P
1
¯
¯
•
u
i
=
(
2
π
d
1
¯
¯
λ
i
)
×
[
sin
(
θ
Δ
P
1
¯
¯
)
cos
⁡
(
φ
Δ
P
1
¯
¯
)
sin
(
θ
i
)
cos
⁡
(
φ
i
)
+
sin
(
θ
Δ
P
1
¯
¯
)
sin
(
φ
Δ
P
1
¯
¯
)
sin
(
θ
i
)
sin
(
φ
i
)
+
cos
⁡
(
θ
Δ
P
1
¯
¯
)
cos
⁡
(
θ
i
)
]
,
(62)
ω
¯
¯
2
i
=
(
2
π
λ
i
)
Δ
P
2
¯
¯
•
u
i
=
(
2
π
d
2
¯
¯
λ
i
)
×
[
sin
(
θ
Δ
P
2
¯
¯
)
cos
⁡
(
φ
Δ
P
2
¯
¯
)
sin
(
θ
i
)
cos
⁡
(
φ
i
)
m
+
sin
(
θ
Δ
P
2
¯
¯
)
sin
(
φ
Δ
P
2
¯
¯
)
sin
(
θ
i
)
sin
(
φ
i
)
m
+
cos
⁡
(
θ
Δ
P
2
¯
¯
)
cos
⁡
(
θ
i
)
]
.
The 2D-DOA estimation can be acquired by solving (61) and (62). The equations are nonlinear and the analytical solution does not exist. Thus, the iterative numerical approximation method is used to obtain the optimal solution. The nonlinear equations can be written as
(63)
f
1
(
x
1
,
x
2
)
=
0
,
f
2
(
x
1
,
x
2
)
=
0
,
where x
1
=
θ and x
2
=
φ. The equation with the same solution of (63) is expressed as
(64)
x
=
f
i
-
1
(
x
1
,
x
2
)
,
i
=
1,2
.
The iterative form can be represented as
(65)
x
k
+
1
=
F
i
(
x
1
k
,
x
2
k
)
,
i
=
1,2
.
The initial vector is selected as x
(
0
)
=
[
x
1
(
0
)
,
x
2
(
0
)
]
T. Calculating (65) until the sequence converges, that is, x
k
→
x
*, and regarding x
* as the iterative solution of (63), we can achieve the 2D-DOA estimation of the incident signal.Assuming that the noise is nonuniform, in theory through calculating the cross-covariance matrices among different arrays, nonuniform noise could be suppressed.


## 6. The Computational Complexity Analysis

We only focus on the major part, which is the number of multiplications involved in calculating covariance matrices and the ALS algorithm. As mentioned above,N, r, and 2
m represent the snapshot, source number, and sensor number, respectively. The ESPRIT algorithm which is similar to the algorithm proposed in [18] needs to calculate the eigendecomposition of covariance matrices and parameter pairing (the covariance matrix is used instead of the four-order cumulant in [18]). The eigendecomposition of three 2
m
×
2
m matrices requires about O
(
24
m
3
). The algorithm in this paper uses COMFAC algorithm to fit an m
×
m
×
4 three-way array. The computational complexity per iteration is O
(
r
3
)
+
O
(
4
m
2
r
). For the simulation in Section 9, the proposed algorithm converges in only two iterations. Therefore, the total computation complexity of the proposed algorithm is O
[
K
(
r
3
+
4
m
2
r
)
] (K stands for the number of iterations). Since the initialization takes up a large proportion of computation time in each iteration, the proposed algorithm suffers from a higher computational complexity.


## 7. CRB of Joint Polarization and DOA Estimation

Cramer-Rao bound (CRB) gives a lower bound of unbiased parameter estimation. In this section, the multiparameter estimation CRB is derived. For the sake of simplicity, the signal covariance matrixR
s is assumed to be known, and the noise is normalized as one. 4
r parameters are contained in the covariance matrix R, that is, r elevation parameters, r azimuth parameters, and 2
r polarization parameters, respectively. In this paper, the polarization parameters of the incident signals are assumed to be known. Only 2
r parameters remain to be estimated. The vector parameters to be estimated are represented as
(66)
v
T
=
[
θ
1
,
φ
1
,
θ
2
,
φ
2
,
…
,
θ
r
,
φ
r
]
.
The CRB of joint polarization parameters and angle parameters is defined as
(67)
E
[
(
v
^
-
v
)
(
v
^
-
v
)
T
]
≥
CRB
,
CRB
=
F
-
1
.
The 2
r
×
2
r Fisher information matrix (FIM) for the parameter v is given by
(68)
F
=
[
F
θ
θ
F
θ
φ
F
φ
θ
F
φ
φ
]
,
where F
θ
θ is the block matrix of elevation estimator and F
φ
φ is the block matrix of azimuth estimator. The entry of ith row and jth column of FIM is represented as
(69)
F
i
j
=
N
trace
[
R
-
1
∂
R
∂
v
i
R
-
1
∂
R
∂
v
j
]
=
2
N
Re
{
trace
[
D
i
R
s
B
H
R
-
1
B
R
s
D
j
H
R
-
1
]
m
m
m
m
m
+
trace
[
D
i
R
s
B
H
R
-
1
B
R
s
D
j
R
-
1
]
}
,
D
i
=
∂
B
∂
v
i
,
where trace
[
•
] is the trace of matrix [
•
] and ∂
R
/
∂
v
i is the partial derivative of matrix R.


## 8. Discussion

The proposed algorithm is able to operate on very irregular shapes, which is not following a circular geometry. However, some requirements have to be satisfied. The design of the array is discussed in two categories: the algorithm has an analytic solution; the algorithm does not have an analytic solution.(
1)  The Algorithm Has an Analytic Solution. Based on the design of array in Figure 2 and (12), (13), it can be seen that the distance vector Δ
P
1 is perpendicular to Δ
P
2. This is one requirement for obtaining the analytic solution. Another requirement is that the distance vector Δ
P
1 or Δ
P
2 is parallel or perpendicular to the coordinate axis. Then, the analytic solution can be obtained. For example, the design of the spherical conformal array is given as follows.The structure of spherical conformal array is shown in Figure7. The elements are arranged on the surface of the spherical conformal array. 1
~
m
-
1 constitute array 1, 2
~
m constitute array 2, m
+
1
~
2
m
-
1 constitute array 3, and m
+
2
~
2
m constitute array 4. The elements of the array are divided into two subarrays. Subarray 1 consists of array 1 and array 2, and subarray 2 consists of array 3 and array 4. The distance vector between array 1 and array 2 is Δ
P
1. The distance vector between array 3 and array 4 is Δ
P
2. As shown in Figures 7 and 8, the distance vector Δ
P
1is perpendicular to Δ
P
2. Also, the distance vector Δ
P
1 is parallel to Y-axis.Figure 7
The structure of spherical conformal array.The schematic diagram of the distance vector.
(a)
The distance vectorΔ
P
1
(b)
The distance vectorΔ
P
2(
2
)  The Algorithm Does Not Have an Analytic Solution. The only requirement is that Δ
P
1 is not parallel to Δ
P
2 as shown in Figure 2. In other words, the subarray 1 is not parallel to subarray 2 as shown in Figure 6. According to (61) and (62), the analytic solution of the proposed algorithm does not exist. Thus, the iterative numerical approximation method is used to obtain the optimal solution.For very irregular shapes array, if we can find two subarrays as mentioned above and the two distance vectors are not parallel to each other, then the proposed algorithm can be used for DOA estimation.


## 9. Simulation Results

In this section, we demonstrate the performance of the proposed algorithm via numerical simulation. Taking a cylindrical conformal array as an example, we use the array structure in Figure2 for simulation. The number of sensors is 16; that is, m
=
8. Without loss of generality, k
1
θ
=
0.5, k
1
φ
=
0.5; k
2
θ
=
0.3, k
2
φ
=
0.7. The sensor patterns are g
i
θ
=
sin
(
θ
j
-
φ
j
) and g
i
φ
=
cos
⁡
(
θ
j
-
φ
j
) in the global coordinate. θ
j and φ
j are the elevation and azimuth of the ith incident signal in the jth local coordinate, respectively. The details regarding how the pattern is transformed from the local coordinate to the global coordinate can be found in [8, 20]. 500 independent trials are considered. Root mean square error (RMSE) which indicates the performance of the proposed algorithm is defined as
(70)
RMSE
=
1
500
∑
t
=
1
500
‍
[
(
θ
^
i
,
t
-
θ
)
2
+
(
φ
^
i
,
t
-
φ
)
2
]
,
where θ
^
i
,
t and φ
^
i
,
t are the elevation and azimuth estimators of the ith incident signal in the tth trial. The ESPRIT algorithm proposed in [18] is simulated in the same scenarios (the ESPRIT algorithm is used for the covariance matrix rather than the four-order cumulant).For the following experiments, the COMFAC algorithm is used to fit them
×
m
×
4 three-way array. The initialization and fitting of the COMFAC are done in compressed space. The Tucker3 three-way model is used in data compression [29].In the first experiment, we show how the success rate and RMSE changed with different SNR. Here, a successful experiment is defined as the experiment with estimation error of less than 1 degree. The elevation and azimuth angles of two far-field narrow incident signals are(
9
5
∘
,
5
0
∘
) and (
10
0
∘
,
6
0
∘
), respectively. 500 snapshots are used in this experiment. The success rate and RMSE are shown in Figures 9(a) and 9(b), respectively. Figure 9(a) shows that the success rate of the proposed PARAFAC algorithm is higher than that of the ESPRIT algorithm when the SNR is low. In addition, the RMSE of the proposed algorithm is much smaller than that of the ESPRIT algorithm at high SNR. As the SNR increases, the RMSE of the proposed algorithm approximates to the CRB. Four covariance matrices are used to estimate DOA in the proposed algorithm. However, only two covariance matrices are used for the ESPRIT algorithm. The proposed algorithm utilizes more data information than the ESPRIT algorithm, which leads to better performance.The estimation performance versus SNR.
(a)
The success rate versus SNR
(b)
The RMSE versus SNRWe plot the curves of success rate and RMSE versus number of snapshots in Figures10(a) and 10(b). SNR is 0 dB and 10 dB in Figures 10(a) and 10(b), respectively. Other simulation conditions are the same as those in the first experiment. Figure 10(a) shows that the success rate of the proposed algorithm is higher than that of the ESPRIT algorithm. In particular, the success rate of source 2 estimated by the ESPRIT algorithm is lower than others. As shown in Figure 10(b), the RMSE of the proposed algorithm is much smaller than that of the ESPRIT algorithm at the same number of snapshots. The RMSE approaches to the CRB as the number of snapshots increases. The RMSE of source 2 with ESPRIT algorithm is relatively larger which is mainly caused by the fact that the ESPRIT algorithm only uses the autocovariance matrices of the array. However, the proposed algorithm achieves higher accuracy by using both autocovariance and cross-covariance matrices of the array.The estimation performance versus number of snapshots.
(a)
The success rate versus number of snapshots
(b)
The RMSE versus number of snapshotsFigure11 shows the resolution probability versus SNR for both the proposed PARAFAC algorithm and the ESPRIT algorithm when the sources are closely spaced (6 degrees separation). The elevation and azimuth of the incident signals are (
10
0
∘
,
5
4
∘
) and (
10
0
∘
,
6
0
∘
), respectively. Other simulation conditions are same as those in the first experiment. The two incident signals are resolved if |
θ
^
1
-
θ
1
|, |
θ
^
2
-
θ
2
| are smaller than |
θ
1
-
θ
2
|
/
2. Meanwhile, |
φ
^
1
-
φ
1
|, |
φ
^
2
-
φ
2
| are smaller than |
φ
1
-
φ
2
|
/
2. θ
^
i and θ
i represent the estimated and real elevations for ith incident signal, respectively. φ
^
i and φ
i represent the estimated and real azimuths for ith incident signal, respectively [30]. Figure 11 shows that the resolution of ESPRIT algorithm outperforms the proposed algorithm in relatively low SNR (−2 dB to 2 dB). When SNR is greater than 2 dB, the proposed algorithm performs better. According to (25) and Theorem 2, to ensure the uniqueness of the PARAFAC model, a
(
:
,
i
)
≠
a
(
:
,
j
) must be guaranteed. In other words, the two incident signals could not be too close. When the SNR is low, the effect of noise is fairly remarkable. All the above reasons result in the relatively poor resolution performance of the proposed algorithm in low SNR. Nevertheless, it is a good compromise for having excellent estimation accuracy and robustness to noise.Figure 11
The resolution probability versus SNR.


## 10. Conclusion

In this paper, a novel high accuracy 2D-DOA estimation algorithm for the conformal array is proposed. The ordinary DOA estimation algorithm cannot be used on conformal array because of the polarization diversity of the varying curvature. To avoid parameter pairing problem, the algorithm forms a PARAFAC model of covariance matrices in the covariance domain to estimate the 2D-DOA. The uniqueness condition of the PARAFAC model is derived. It can be generalized to other conformal array structure, that is, conical conformal array, and so forth. Computer simulation verifies the effectiveness of the proposed algorithm in terms of accuracy and robustness to noise.



---

*Source: 394707-2014-01-16.xml*