# Joint 2D Direction-of-Arrival and Range Estimation for Nonstationary Sources

Authors: Jian Chen; Hui Zhao; Xiaoying Sun; Guohong Liu
Journal: International Journal of Antennas and Propagation
 (2014)

---

## Abstract

Passive localization of nonstationary sources in the spherical coordinates (azimuth, elevation, and range) is considered, and a parallel factor analysis based method is addressed for the near-field parameter estimation problem. In this scheme, a parallel factor analysis model is firstly constructed by computing five time-frequency distribution matrices of the properly chosen observation data. In addition, the uniqueness of the constructed model is proved, and both the two-dimensional (2D) direction-of-arrival (DOA) and range can be jointly obtained via trilinear alternating least squares regression (TALS). The investigated algorithm is well suitable for near-field nonstationary source localization and does not require parameter-pairing or multidimensional search. Several simulation examples confirm the effectiveness of the proposed algorithm.


---

## Body

## 1. Introduction

Bearing estimation has been a strong interest in radar and sonar as well as communication. In the last three decades, various high-resolution algorithms for direction finding of multiple narrowband sources assume that the propagating waves are considered to be plane waves at the sensor array. However, when the sources are located in the Fresnel region [1] of the array aperture, the wavefronts emitted from these sources are spherical rather than planar at each sensor position and characterized by both the DOA and the range parameters; thus, the existing DOA estimation schemes, such as MUSIC and ESPRIT [2], would fail in estimating near-field localization parameters.By applying the Fresnel approximation to the near-field sources localization, the two-dimensional (2D) MUSIC method, the high-order ESPRIT method, and the path-following method were, respectively, proposed in [3–7] in order to cope with the problem of estimating azimuth and range. In recent years, several methods that can obtain azimuth, elevation, and range have been developed. For instance, in [8], an original 3D higher order statistics based localization algorithm has been presented. By translating the 1D uniform linear array of near-field into a virtual rectangular array of virtual far-field, Challa and Shamsunder [9] proposed a unitary-ESPRIT method. Aberd Meraim and Hua [10] used only the second-order statistics and proposed a higher resolution 3D near-field source localization; however, a parameter-pairing process leading to the poor performance in lower signal-to-noise ratio had to be taken into account. References [11, 12] also proposed fourth-order cumulant based algorithms to estimate 2D DOA and range, but they all coincided with [10]. All abovementioned algorithms rely on the assumption that impinging source signals are stationary; when nonstationary FM signals exist, they will show an unsatisfactory performance.While quadratic time-frequency distribution [13, 14] has been sought out and properly investigated into sensor and spatial signal processing, and its evaluation of the observation data across the sensor array yields spatial time-frequency distribution (STFDs), the main advantage of STFDs is that it can well handle signals of nonstationary waveforms that are highly localized in the time-frequency domain and effectively improve the robustness of localization methods by spreading the noise power into the whole time-frequency domain. The STFDs based algorithm to locate near-field nonstationary sources has been presented in [15] and showed a satisfactory parameters estimation accuracy; however, it required 2D search and only estimated 1D DOA and range.In this paper, by exploiting favorable characteristics of a uniform cross array, we present a joint 2D DOA and range estimation algorithm. We first compute five time-frequency matrices to construct a parallel factor (PARAFAC) analysis model. Then, we obtain three-dimensional (3D) near-field parameters via trilinear alternating least squares regression (TALS). Compared with the other methods, the main contribution for the proposed method can be summarized as follows:(
1
) we obtain 3D near-field sources parameters (elevations, azimuths, and ranges) of nonstationary signals rather than stationary waves; (
2
) we creatively incorporate STFDs with parallel factor analysis to well avoid both parameter pairing and multidimensional search.The rest of this paper is organized as follows. Section2 introduces the signal model of near-field localization based on cross array. Section 3 develops a joint estimation algorithm of three parameters in near-field. Section 4 shows simulation results. Section 5 presents the conclusion of the whole paper.


## 2. Near-Field Signal Model Based Cross Array

### 2.1. Near-Field Signal Model

We consider a near-field scenario ofM uncorrelated narrowband signals impinging on a cross array signed with the x- and y-axes (Figure 1), which consists of L
=
4
P
+
1 elements with interelement spacing d. Let the array center be the phase point; the signals received by the (
l
,
0
)th and the (
0
,
l
)th can be, respectively, expressed as
(1)
x
l
,
0
(
t
)
=
z
l
,
0
(
t
)
+
n
l
,
0
(
t
)
=
∑
m
=
1
M
s
m
(
t
)
e
j
(
γ
x
m
l
+
φ
x
m
l
2
)
+
n
l
,
0
(
t
)
,
(2)
x
0
,
l
(
t
)
=
z
0
,
l
(
t
)
+
n
0
,
l
(
t
)
=
∑
m
=
1
M
s
m
(
t
)
e
j
(
γ
j
m
l
+
φ
j
m
l
2
)
+
n
0
,
l
(
t
)
,
with n
(
t
) being additive noise, and
(3)
γ
x
m
=
-
2
π
d
λ
sin
α
m
cos
⁡
θ
m
,
ϕ
x
m
=
π
d
2
λ
r
m
(
1
-
si
n
2
α
m
cos
⁡
2
⁡
θ
m
)
,
γ
y
m
=
-
2
π
d
λ
sin
α
m
sin
θ
m
,
ϕ
y
m
=
π
d
2
λ
r
m
(
1
-
si
n
2
α
m
si
n
2
θ
m
)
,
where α
m, θ
m, and r
m indicate elevation, azimuth, and range of mth signal, respectively, and λ is wavelength of source signal.Figure 1
Sensor-source configuration for the near-field problem.Themth source signal with phase ψ
m
(
t
) can be modeled as
(4)
s
m
(
t
)
=
D
m
e
j
ψ
m
(
t
)
,
where m
=
1,2
,
…
,
M and D
m is the amplitude of the mth source signal.


### 2.2. Assumption of Signal Model

The main problem addressed in this paper is to jointly estimate the sets of parameters(
α
m
,
θ
m
,
r
m
); then the following assumptions are assumed to hold.(1)
The source signals
(
t
) is narrowband, independent, and nonstationary process.
(2)
The additive noise is spatially white Gaussian with zero-mean and independent from the source signals.
(3)
For unique estimation, we required
=
λ
/
4, M
≤
N
-
1, and -
γ
x
m
+
ϕ
x
m
≠
γ
x
m
+
ϕ
x
m as well as -
γ
y
m
+
ϕ
y
m
≠
γ
y
m
+
ϕ
y
m.


### 2.3. Parallel Factor Analysis

We need to introduce the following notation that will be used in the sequel.Definition 1 (see [16]).
Letx
k
,
n
,
p stand for the (
k
,
n
,
p
) element of a three-dimensional tensor X, if
(5)
x
k
,
n
,
p
=
∑
m
=
1
M
a
k
,
m
b
n
,
m
c
p
,
m
,
where a
k
,
m denotes the (
k
,
m
) element of matrix A, and similarly for the others. Equation (5) indicates x
k
,
n
,
p as a sum of triple products, which is variably known as the trilinear model, trilinear decomposition, triple product decomposition, canonical decomposition, or parallel factor (PARAFAC) analysis.Definition 2 (see [16]).
For a matrixB
∈
C
I
×
J, if all I
<
J columns of B are linearly independent but there exists a collection of I
+
1 linearly dependent columns of B, then it has Kruskal-rank (k-rank) k
B
=
I.Theorem 3 (see [16]).
Consider a three-dimensional tensorX as defined in (5), and M represents the common dimension; if
(6)
k
A
+
k
B
+
k
C
T
≥
2
M
+
2
,
then A, B
, and C are unique up to permutation and (complex) scaling of columns.


## 2.1. Near-Field Signal Model

We consider a near-field scenario ofM uncorrelated narrowband signals impinging on a cross array signed with the x- and y-axes (Figure 1), which consists of L
=
4
P
+
1 elements with interelement spacing d. Let the array center be the phase point; the signals received by the (
l
,
0
)th and the (
0
,
l
)th can be, respectively, expressed as
(1)
x
l
,
0
(
t
)
=
z
l
,
0
(
t
)
+
n
l
,
0
(
t
)
=
∑
m
=
1
M
s
m
(
t
)
e
j
(
γ
x
m
l
+
φ
x
m
l
2
)
+
n
l
,
0
(
t
)
,
(2)
x
0
,
l
(
t
)
=
z
0
,
l
(
t
)
+
n
0
,
l
(
t
)
=
∑
m
=
1
M
s
m
(
t
)
e
j
(
γ
j
m
l
+
φ
j
m
l
2
)
+
n
0
,
l
(
t
)
,
with n
(
t
) being additive noise, and
(3)
γ
x
m
=
-
2
π
d
λ
sin
α
m
cos
⁡
θ
m
,
ϕ
x
m
=
π
d
2
λ
r
m
(
1
-
si
n
2
α
m
cos
⁡
2
⁡
θ
m
)
,
γ
y
m
=
-
2
π
d
λ
sin
α
m
sin
θ
m
,
ϕ
y
m
=
π
d
2
λ
r
m
(
1
-
si
n
2
α
m
si
n
2
θ
m
)
,
where α
m, θ
m, and r
m indicate elevation, azimuth, and range of mth signal, respectively, and λ is wavelength of source signal.Figure 1
Sensor-source configuration for the near-field problem.Themth source signal with phase ψ
m
(
t
) can be modeled as
(4)
s
m
(
t
)
=
D
m
e
j
ψ
m
(
t
)
,
where m
=
1,2
,
…
,
M and D
m is the amplitude of the mth source signal.


## 2.2. Assumption of Signal Model

The main problem addressed in this paper is to jointly estimate the sets of parameters(
α
m
,
θ
m
,
r
m
); then the following assumptions are assumed to hold.(1)
The source signals
(
t
) is narrowband, independent, and nonstationary process.
(2)
The additive noise is spatially white Gaussian with zero-mean and independent from the source signals.
(3)
For unique estimation, we required
=
λ
/
4, M
≤
N
-
1, and -
γ
x
m
+
ϕ
x
m
≠
γ
x
m
+
ϕ
x
m as well as -
γ
y
m
+
ϕ
y
m
≠
γ
y
m
+
ϕ
y
m.


## 2.3. Parallel Factor Analysis

We need to introduce the following notation that will be used in the sequel.Definition 1 (see [16]).
Letx
k
,
n
,
p stand for the (
k
,
n
,
p
) element of a three-dimensional tensor X, if
(5)
x
k
,
n
,
p
=
∑
m
=
1
M
a
k
,
m
b
n
,
m
c
p
,
m
,
where a
k
,
m denotes the (
k
,
m
) element of matrix A, and similarly for the others. Equation (5) indicates x
k
,
n
,
p as a sum of triple products, which is variably known as the trilinear model, trilinear decomposition, triple product decomposition, canonical decomposition, or parallel factor (PARAFAC) analysis.Definition 2 (see [16]).
For a matrixB
∈
C
I
×
J, if all I
<
J columns of B are linearly independent but there exists a collection of I
+
1 linearly dependent columns of B, then it has Kruskal-rank (k-rank) k
B
=
I.Theorem 3 (see [16]).
Consider a three-dimensional tensorX as defined in (5), and M represents the common dimension; if
(6)
k
A
+
k
B
+
k
C
T
≥
2
M
+
2
,
then A, B
, and C are unique up to permutation and (complex) scaling of columns.


## 3. PARAFAC Based 3D Near-Field Sources Localization

### 3.1. Computation of the Spatial Time-Frequency Distribution Matrices

The discrete form of Cohen’s class of time-frequency distribution of a signalx
(
t
) can be expressed as
(7)
D
X
X
(
t
,
f
)
=
∑
k
=
-
∞
∞
∑
τ
=
-
∞
∞
φ
(
k
,
τ
)
x
(
t
+
k
+
τ
)
gggggggggg
·
x
*
(
t
+
k
-
τ
)
e
-
j
4
π
f
τ
,
where φ
(
k
,
τ
) is the time-frequency kernel and the superscript * denotes complex conjugate. Replacing x
(
t
) by the data snapshot x
-
l
-
1,0
(
t
) and x
-
l
,
0
(
t
), we obtain
(8)
D
-
l
-
1
,
-
l
x
(
t
,
f
)
=
∑
k
=
-
∞
∞
∑
τ
=
-
∞
∞
φ
(
k
,
τ
)
x
-
l
-
1
(
t
+
k
+
τ
)
g
g
g
g
g
g
g
g
g
·
x
-
l
,
0
*
(
t
+
k
-
τ
)
e
-
j
4
π
f
τ
.Substituting (1) into (8), D
-
l
-
1
,
-
l
x
(
t
,
f
) can be extended to the following form:
(9)
D
-
l
-
1
,
-
l
x
(
t
,
f
)
=
D
z
-
l
-
1
(
t
)
,
z
-
l
(
t
)
(
t
,
f
)
+
D
z
-
l
-
1
(
t
)
,
n
-
l
(
t
)
(
t
,
f
)
+
D
n
-
l
-
1
(
t
)
,
z
-
l
(
t
)
(
t
,
f
)
+
D
n
-
l
-
1
(
t
)
,
n
-
l
(
t
)
(
t
,
f
)
.Under the assumptions(
1
) and (
2
), it is obvious that
(10)
E
[
D
-
l
-
1
,
-
l
x
(
t
,
f
)
]
=
E
[
D
z
-
l
-
1
(
t
)
,
z
-
l
(
t
)
(
t
,
f
)
]
=
∑
m
=
1
M
E
[
D
s
m
(
t
,
f
)
]
e
j
(
-
γ
x
m
+
ϕ
x
m
)
e
j
2
l
ϕ
x
m
,
where D
s
m
(
t
,
f
) indicates the STFDs of source s
m
(
t
).Using a rectangular window of old lengthN, the pseudo Wigner-Ville distribution (PWVD) of s
m
(
t
) is given by
(11)
D
s
m
(
t
,
f
)
=
∑
τ
=
-
(
N
-
1
)
/
2
(
N
-
1
)
/
2
x
m
(
t
+
τ
)
x
m
*
(
t
-
τ
)
e
-
j
4
π
f
τ
.Assume that the third-order derivative of the phase can be negligible over the rectangular window lengthN, and f
m
=
d
ψ
m
(
t
)
/
(
2
π
d
t
), ψ
m
(
t
+
τ
)
-
ψ
m
(
t
+
τ
)
-
4
π
f
m
τ
=
0; then we obtain the approximated expression as
(12)
D
s
m
(
t
,
f
)
=
∑
τ
=
-
(
N
-
1
)
/
2
(
N
-
1
)
/
2
D
M
2
=
L
D
M
2
.We construct matrixD
1 with the (
k
,
q
)th element being given by
(13)
D
1
(
k
,
q
)
=
E
[
D
q
-
k
-
1
,
q
-
k
x
(
t
,
f
)
]
=
∑
m
=
1
M
E
[
D
s
m
(
t
,
f
)
]
e
j
(
-
γ
x
m
+
ϕ
x
m
)
e
j
2
(
k
-
q
)
ϕ
x
m
.On the other hand, following the same process described above, we can easily obtain(14)
D
2
(
k
,
q
)
=
E
[
D
l
+
1
,
l
x
(
t
,
f
)
]
=
E
[
D
k
-
q
+
1
,
k
-
q
x
(
t
,
f
)
]
=
∑
m
=
1
M
E
[
D
s
m
(
t
,
f
)
]
e
j
(
γ
x
m
+
ϕ
x
m
)
e
j
2
(
k
-
q
)
ϕ
x
m
,
D
3
(
k
,
q
)
=
E
[
D
-
l
-
1
,
-
l
y
(
t
,
f
)
]
=
E
[
D
q
-
k
-
1
,
q
-
k
y
(
t
,
f
)
]
=
∑
m
=
1
M
E
[
D
s
m
(
t
,
f
)
]
e
j
(
-
γ
y
m
+
ϕ
y
m
)
e
j
2
(
k
-
q
)
ϕ
y
m
,
D
4
(
k
,
q
)
=
E
[
D
l
+
1
,
l
y
(
t
,
f
)
]
=
E
[
D
k
-
q
+
1
,
k
-
q
y
(
t
,
f
)
]
=
∑
m
=
1
M
E
[
D
s
m
(
t
,
f
)
]
e
j
(
γ
y
m
+
ϕ
y
m
)
e
j
2
(
k
-
q
)
ϕ
y
m
,
D
5
(
k
,
q
)
=
E
[
D
l
,
l
+
1
y
(
t
,
f
)
]
=
E
[
D
q
-
k
,
q
-
k
+
1
y
(
t
,
f
)
]
=
∑
m
=
1
M
E
[
D
s
m
(
t
,
f
)
]
e
j
(
-
γ
y
m
-
ϕ
y
m
)
e
-
j
2
(
q
-
k
)
ϕ
y
m
.
And their matrices form becomes
(15)
D
1
(
t
,
f
)
=
A
D
s
(
t
,
f
)
Ω
1
*
Φ
1
A
H
,
D
2
(
t
,
f
)
=
A
D
s
(
t
,
f
)
Ω
1
Φ
1
A
H
,
D
3
(
t
,
f
)
=
A
D
s
(
t
,
f
)
Ω
2
*
Φ
2
A
H
,
D
4
(
t
,
f
)
=
A
D
s
(
t
,
f
)
Ω
2
Φ
2
A
H
,
D
5
(
t
,
f
)
=
A
D
s
(
t
,
f
)
Ω
2
*
Φ
2
*
A
H
,
where A
=
[
a
1
,
a
2
,
…
,
a
M
], a
m
=
[
1
,
e
j
2
ϕ
x
m
,
…
,
e
j
2
(
N
-
1
)
ϕ
x
m
], Ω
1
=
diag
⁡
(
e
j
γ
x
1
,
e
j
γ
x
2
,
…
,
e
j
γ
x
M
), Φ
1
=
diag
⁡
(
e
j
ϕ
x
1
,
e
j
ϕ
x
2
,
…
,
e
j
ϕ
x
M
), Ω
2
=
diag
⁡
(
e
j
γ
y
1
,
e
j
γ
y
2
,
…
,
e
j
γ
y
M
), and Φ
2
=
diag
⁡
(
e
j
ϕ
y
1
,
e
j
ϕ
y
2
,
…
,
e
j
ϕ
y
M
).


### 3.2. Construction of the Parallel Factor Analysis Model

Considering the situation of limited samples, we build a parallel factor analysis model that uses the spatial time-frequency distribution as(16)
D
~
=
[
D
1
(
t
,
f
)
D
2
(
t
,
f
)
D
3
(
t
,
f
)
D
4
(
t
,
f
)
D
5
(
t
,
f
)
]
=
[
A
D
s
(
t
,
f
)
Ω
1
*
Φ
1
A
H
A
D
s
(
t
,
f
)
Ω
1
Φ
1
A
H
A
D
s
(
t
,
f
)
Ω
2
*
Φ
2
A
H
A
D
s
(
t
,
f
)
Ω
2
Φ
2
A
H
A
D
s
(
t
,
f
)
Ω
2
*
Φ
2
*
A
H
]
+
ν
1
.LettingC
=
A
*, the Khutrr-Rao product [16] for (16) shows
(17)
D
~
=
(
M
⊗
A
)
C
T
+
ν
1
,
where
(18)
M
=
[
g
-
1
(
D
s
(
t
,
f
)
Ω
1
*
Φ
1
)
g
-
1
(
D
s
(
t
,
f
)
Ω
1
Φ
1
)
g
-
1
(
D
s
(
t
,
f
)
Ω
2
*
Φ
2
)
g
-
1
(
D
s
(
t
,
f
)
Ω
2
Φ
2
)
g
-
1
(
D
s
(
t
,
f
)
Ω
2
*
Φ
2
*
)
]
,
g
-
1
(
D
s
(
t
,
f
)
Ω
1
*
Φ
1
) denoting a row vector consisting of diagonal matrix D
s
(
t
,
f
)
Ω
1
*
Φ
1.Similarly, (17) also yields
(19)
X
~
=
(
A
⊗
C
)
M
T
+
ν
2
,
Y
~
=
(
C
⊗
M
)
A
T
+
ν
3
.


### 3.3. Estimation of 2D Direction-of-Arrival and Range

As it stands,A and C are both Vander-monde matrices, and then they have Kruskal-rank (k-rank) k
A
=
k
C
T
=
M. On the other hand, the k-rank of M will be k
D
=
min
⁡
(
5
,
M
). When the condition that the number of signals being M
≥
2 holds, then A, C
, and M are unique up to permutation and scaling of columns. With trilinear alternating least squares regression, we obtain that
(20)
C
~
T
=
arg
min
⁡
C
T
∥
D
~
-
(
M
⊗
A
)
C
T
∥
F
2
,
M
~
T
=
arg
min
⁡
M
T
∥
X
~
-
(
A
⊗
C
)
M
T
∥
F
2
,
A
~
T
=
arg
min
⁡
A
T
∥
Y
~
-
(
C
⊗
M
)
A
T
∥
F
2
.Then using these estimates, we can get each pair(
γ
x
m
,
γ
y
m
,
φ
y
m
) as follows:
(21)
γ
~
x
m
=
1
2
arg
(
M
~
(
2
,
m
)
M
~
(
1
,
m
)
)
,
γ
~
y
m
=
1
2
arg
(
M
~
(
4
,
m
)
M
~
(
3
,
m
)
)
,
φ
~
y
m
=
1
2
arg
(
M
~
(
3
,
m
)
M
~
(
5
,
m
)
)
.Finally, the sources parameters can be estimated as(22)
α
~
m
=
a
sin
⁡
(
λ
2
π
d
(
γ
~
x
m
2
+
γ
~
y
m
2
)
1
/
2
)
,
θ
~
m
=
a
tan
⁡
(
γ
~
y
m
γ
~
x
m
)
,
r
~
m
=
π
d
2
(
1
-
si
n
2
α
~
m
si
n
2
θ
~
m
)
λ
ϕ
~
y
m
.


## 3.1. Computation of the Spatial Time-Frequency Distribution Matrices

The discrete form of Cohen’s class of time-frequency distribution of a signalx
(
t
) can be expressed as
(7)
D
X
X
(
t
,
f
)
=
∑
k
=
-
∞
∞
∑
τ
=
-
∞
∞
φ
(
k
,
τ
)
x
(
t
+
k
+
τ
)
gggggggggg
·
x
*
(
t
+
k
-
τ
)
e
-
j
4
π
f
τ
,
where φ
(
k
,
τ
) is the time-frequency kernel and the superscript * denotes complex conjugate. Replacing x
(
t
) by the data snapshot x
-
l
-
1,0
(
t
) and x
-
l
,
0
(
t
), we obtain
(8)
D
-
l
-
1
,
-
l
x
(
t
,
f
)
=
∑
k
=
-
∞
∞
∑
τ
=
-
∞
∞
φ
(
k
,
τ
)
x
-
l
-
1
(
t
+
k
+
τ
)
g
g
g
g
g
g
g
g
g
·
x
-
l
,
0
*
(
t
+
k
-
τ
)
e
-
j
4
π
f
τ
.Substituting (1) into (8), D
-
l
-
1
,
-
l
x
(
t
,
f
) can be extended to the following form:
(9)
D
-
l
-
1
,
-
l
x
(
t
,
f
)
=
D
z
-
l
-
1
(
t
)
,
z
-
l
(
t
)
(
t
,
f
)
+
D
z
-
l
-
1
(
t
)
,
n
-
l
(
t
)
(
t
,
f
)
+
D
n
-
l
-
1
(
t
)
,
z
-
l
(
t
)
(
t
,
f
)
+
D
n
-
l
-
1
(
t
)
,
n
-
l
(
t
)
(
t
,
f
)
.Under the assumptions(
1
) and (
2
), it is obvious that
(10)
E
[
D
-
l
-
1
,
-
l
x
(
t
,
f
)
]
=
E
[
D
z
-
l
-
1
(
t
)
,
z
-
l
(
t
)
(
t
,
f
)
]
=
∑
m
=
1
M
E
[
D
s
m
(
t
,
f
)
]
e
j
(
-
γ
x
m
+
ϕ
x
m
)
e
j
2
l
ϕ
x
m
,
where D
s
m
(
t
,
f
) indicates the STFDs of source s
m
(
t
).Using a rectangular window of old lengthN, the pseudo Wigner-Ville distribution (PWVD) of s
m
(
t
) is given by
(11)
D
s
m
(
t
,
f
)
=
∑
τ
=
-
(
N
-
1
)
/
2
(
N
-
1
)
/
2
x
m
(
t
+
τ
)
x
m
*
(
t
-
τ
)
e
-
j
4
π
f
τ
.Assume that the third-order derivative of the phase can be negligible over the rectangular window lengthN, and f
m
=
d
ψ
m
(
t
)
/
(
2
π
d
t
), ψ
m
(
t
+
τ
)
-
ψ
m
(
t
+
τ
)
-
4
π
f
m
τ
=
0; then we obtain the approximated expression as
(12)
D
s
m
(
t
,
f
)
=
∑
τ
=
-
(
N
-
1
)
/
2
(
N
-
1
)
/
2
D
M
2
=
L
D
M
2
.We construct matrixD
1 with the (
k
,
q
)th element being given by
(13)
D
1
(
k
,
q
)
=
E
[
D
q
-
k
-
1
,
q
-
k
x
(
t
,
f
)
]
=
∑
m
=
1
M
E
[
D
s
m
(
t
,
f
)
]
e
j
(
-
γ
x
m
+
ϕ
x
m
)
e
j
2
(
k
-
q
)
ϕ
x
m
.On the other hand, following the same process described above, we can easily obtain(14)
D
2
(
k
,
q
)
=
E
[
D
l
+
1
,
l
x
(
t
,
f
)
]
=
E
[
D
k
-
q
+
1
,
k
-
q
x
(
t
,
f
)
]
=
∑
m
=
1
M
E
[
D
s
m
(
t
,
f
)
]
e
j
(
γ
x
m
+
ϕ
x
m
)
e
j
2
(
k
-
q
)
ϕ
x
m
,
D
3
(
k
,
q
)
=
E
[
D
-
l
-
1
,
-
l
y
(
t
,
f
)
]
=
E
[
D
q
-
k
-
1
,
q
-
k
y
(
t
,
f
)
]
=
∑
m
=
1
M
E
[
D
s
m
(
t
,
f
)
]
e
j
(
-
γ
y
m
+
ϕ
y
m
)
e
j
2
(
k
-
q
)
ϕ
y
m
,
D
4
(
k
,
q
)
=
E
[
D
l
+
1
,
l
y
(
t
,
f
)
]
=
E
[
D
k
-
q
+
1
,
k
-
q
y
(
t
,
f
)
]
=
∑
m
=
1
M
E
[
D
s
m
(
t
,
f
)
]
e
j
(
γ
y
m
+
ϕ
y
m
)
e
j
2
(
k
-
q
)
ϕ
y
m
,
D
5
(
k
,
q
)
=
E
[
D
l
,
l
+
1
y
(
t
,
f
)
]
=
E
[
D
q
-
k
,
q
-
k
+
1
y
(
t
,
f
)
]
=
∑
m
=
1
M
E
[
D
s
m
(
t
,
f
)
]
e
j
(
-
γ
y
m
-
ϕ
y
m
)
e
-
j
2
(
q
-
k
)
ϕ
y
m
.
And their matrices form becomes
(15)
D
1
(
t
,
f
)
=
A
D
s
(
t
,
f
)
Ω
1
*
Φ
1
A
H
,
D
2
(
t
,
f
)
=
A
D
s
(
t
,
f
)
Ω
1
Φ
1
A
H
,
D
3
(
t
,
f
)
=
A
D
s
(
t
,
f
)
Ω
2
*
Φ
2
A
H
,
D
4
(
t
,
f
)
=
A
D
s
(
t
,
f
)
Ω
2
Φ
2
A
H
,
D
5
(
t
,
f
)
=
A
D
s
(
t
,
f
)
Ω
2
*
Φ
2
*
A
H
,
where A
=
[
a
1
,
a
2
,
…
,
a
M
], a
m
=
[
1
,
e
j
2
ϕ
x
m
,
…
,
e
j
2
(
N
-
1
)
ϕ
x
m
], Ω
1
=
diag
⁡
(
e
j
γ
x
1
,
e
j
γ
x
2
,
…
,
e
j
γ
x
M
), Φ
1
=
diag
⁡
(
e
j
ϕ
x
1
,
e
j
ϕ
x
2
,
…
,
e
j
ϕ
x
M
), Ω
2
=
diag
⁡
(
e
j
γ
y
1
,
e
j
γ
y
2
,
…
,
e
j
γ
y
M
), and Φ
2
=
diag
⁡
(
e
j
ϕ
y
1
,
e
j
ϕ
y
2
,
…
,
e
j
ϕ
y
M
).


## 3.2. Construction of the Parallel Factor Analysis Model

Considering the situation of limited samples, we build a parallel factor analysis model that uses the spatial time-frequency distribution as(16)
D
~
=
[
D
1
(
t
,
f
)
D
2
(
t
,
f
)
D
3
(
t
,
f
)
D
4
(
t
,
f
)
D
5
(
t
,
f
)
]
=
[
A
D
s
(
t
,
f
)
Ω
1
*
Φ
1
A
H
A
D
s
(
t
,
f
)
Ω
1
Φ
1
A
H
A
D
s
(
t
,
f
)
Ω
2
*
Φ
2
A
H
A
D
s
(
t
,
f
)
Ω
2
Φ
2
A
H
A
D
s
(
t
,
f
)
Ω
2
*
Φ
2
*
A
H
]
+
ν
1
.LettingC
=
A
*, the Khutrr-Rao product [16] for (16) shows
(17)
D
~
=
(
M
⊗
A
)
C
T
+
ν
1
,
where
(18)
M
=
[
g
-
1
(
D
s
(
t
,
f
)
Ω
1
*
Φ
1
)
g
-
1
(
D
s
(
t
,
f
)
Ω
1
Φ
1
)
g
-
1
(
D
s
(
t
,
f
)
Ω
2
*
Φ
2
)
g
-
1
(
D
s
(
t
,
f
)
Ω
2
Φ
2
)
g
-
1
(
D
s
(
t
,
f
)
Ω
2
*
Φ
2
*
)
]
,
g
-
1
(
D
s
(
t
,
f
)
Ω
1
*
Φ
1
) denoting a row vector consisting of diagonal matrix D
s
(
t
,
f
)
Ω
1
*
Φ
1.Similarly, (17) also yields
(19)
X
~
=
(
A
⊗
C
)
M
T
+
ν
2
,
Y
~
=
(
C
⊗
M
)
A
T
+
ν
3
.


## 3.3. Estimation of 2D Direction-of-Arrival and Range

As it stands,A and C are both Vander-monde matrices, and then they have Kruskal-rank (k-rank) k
A
=
k
C
T
=
M. On the other hand, the k-rank of M will be k
D
=
min
⁡
(
5
,
M
). When the condition that the number of signals being M
≥
2 holds, then A, C
, and M are unique up to permutation and scaling of columns. With trilinear alternating least squares regression, we obtain that
(20)
C
~
T
=
arg
min
⁡
C
T
∥
D
~
-
(
M
⊗
A
)
C
T
∥
F
2
,
M
~
T
=
arg
min
⁡
M
T
∥
X
~
-
(
A
⊗
C
)
M
T
∥
F
2
,
A
~
T
=
arg
min
⁡
A
T
∥
Y
~
-
(
C
⊗
M
)
A
T
∥
F
2
.Then using these estimates, we can get each pair(
γ
x
m
,
γ
y
m
,
φ
y
m
) as follows:
(21)
γ
~
x
m
=
1
2
arg
(
M
~
(
2
,
m
)
M
~
(
1
,
m
)
)
,
γ
~
y
m
=
1
2
arg
(
M
~
(
4
,
m
)
M
~
(
3
,
m
)
)
,
φ
~
y
m
=
1
2
arg
(
M
~
(
3
,
m
)
M
~
(
5
,
m
)
)
.Finally, the sources parameters can be estimated as(22)
α
~
m
=
a
sin
⁡
(
λ
2
π
d
(
γ
~
x
m
2
+
γ
~
y
m
2
)
1
/
2
)
,
θ
~
m
=
a
tan
⁡
(
γ
~
y
m
γ
~
x
m
)
,
r
~
m
=
π
d
2
(
1
-
si
n
2
α
~
m
si
n
2
θ
~
m
)
λ
ϕ
~
y
m
.


## 4. Computer Simulation Results

In this section, we explicit several simulation results to evaluate the performance of proposed method. For all examples, a symmetrical cross array with a number of 17 elements and interelements spacing of0.25
λ, where λ is the wavelength of the narrowband source signals. The noise used in this section is zero-mean, Gaussian distributed, and temporally white, and the root mean square error (RMSE) is defined as (23). All the following presented results are obtained by averaging the results of 200 independent Monte Carlo simulations. Consider
(23)
RMSE
(
θ
)
=
1
N
(
θ
-
θ
~
)
2
,
where N denotes the number of independent Monte Carlo simulations.In the first example, we examine the performance of the elevation, azimuth, and range estimations accuracy versus the SNR. The snapshot number is set at 512. Two linear frequency-modulated signals arrival at the sensor array with start and end frequencies(
0.1,0
) and (
0.5,0.4
), their 3D near-field parameters locate at (
2
5
∘
,
1
5
∘
,
0.3
λ
) and (
1
0
∘
,
3
0
∘
,
0.4
λ
). For the comparison, the fourth-order cumulant based method [9] is also displayed, in which the source signals are non-Gaussian and stationary process. When the SNR varies from 0 dB to 25 dB, the RMSEs of the 2D direction-of-arrival and range estimations using the proposed method and the fourth-order cumulant based method are shown in Figure 2. From Figure 2, we can see that the proposed method outperforms the fourth-order cumulant based method in elevation and azimuth as well as range estimation for all available SNRs. In addition, the RMSE of range estimations for the first source that is closer to the array is less than the second one. This phenomenon is in well agreement with the theoretical analysis that the sources closer to sensor array would hold a smaller standard deviation than the ones far away from the array.The RMSEs of the elevation, azimuth, and range estimation using the proposed method and the fourth-order cumulant based method versus SNRs.
(a)
Elevation
(b)
Azimuth
(c)
RangeIn the second example, the proposed method is used to deal with the situation that two near-field FM signals are impinging on the sensor array shown in Figure1. The elevation, azimuth, and range of the near-field sources are located at (
3
5
∘
,
4
0
∘
,
1
/
6
λ
) and (
2
0
∘
,
6
0
∘
,
0.4
λ
). Moreover, the snapshot number and SNR are set at 512 and 10 dB. Table 1 establishes the mean and variance of the elevation, azimuth, and range estimations using the proposed method. From Table 1, we can see that the proposed method shows a satisfactory performance in localizing the 3D near-field nonstationary sources.Table 1
The mean and variance of the proposed method for the second example.
True
Mean
Variance
Source 1
Elevation (°)
35
34.1036
0.0854
Azimuth (°)
40
38.8652
0.0514
Range
1/6λ
0.1956λ
8.7246
×
10
-
7
Source 2
Elevation (°)
20
19.8274
0.0963
Azimuth (°)
60
59.3610
0.1827
Range
0.4λ
0.4449λ
7.6641
×
10
-
6In the last example, we consider the situation when far-field and near-field nonstationary sources are incoming on the sensor array mentioned above, and they are located at(
3
5
∘
,
4
0
∘
,
1
/
6
λ
) and (
2
0
∘
,
6
0
∘
,
inf
⁡
), respectively. The snapshot number is fixed at 512. When the SNR varies from 0 dB to 25 dB, the RMSEs of the 2D direction-of-arrival and range estimations using the proposed method are shown in Figure 3. In addition, the results of the fourth-order cumulant based method are also displayed in the same figure for comparison. From Figure 3, it can be seen that the proposed method still performs a satisfactory estimation accuracy for the case that far-field and near-field sources exist simultaneously. However, the fourth-order cumulant based method shows a poor performance in the same situation.The RMSEs of the elevation, azimuth, and range estimation using the proposed method and the fourth-order cumulant based method versus SNRs.
(a)
Elevation
(b)
Azimuth
(c)
Range


## 5. Conclusion

We have developed a spatial time-frequency distribution based algorithm for 3D near-field nonstationary source localization problems. Additionally, with parallel factor analysis technique, there is no parameter pairing or multidimensional searching. Finally, the computer simulation results indicate that using spatial time-frequency distribution and parallel factor together significantly solves the problem of the joint estimation of elevation, azimuth, and range of nonstationary signals. However, the spatial time-frequency averaging methods may lead to the additional computation load.



---

*Source: 849039-2014-07-17.xml*