# Structure-Aware Bayesian Compressive Sensing for Near-Field Source Localization Based on Sensor-Angle Distributions

Authors: Si Qin; Yimin D. Zhang; Qisong Wu; Moeness G. Amin
Journal: International Journal of Antennas and Propagation
 (2015)

---

## Abstract

A novel technique for localization of narrowband near-field sources is presented. The technique utilizes the sensor-angle distribution (SAD) that treats the source range and direction-of-arrival (DOA) information as sensor-dependent phase progression. The SAD draws parallel to quadratic time-frequency distributions and, as such, is able to reveal the changes in the spatial frequency over sensor positions. For a moderate source range, the SAD signature is of a polynomial shape, thus simplifying the parameter estimation. Both uniform and sparse linear arrays are considered in this work. To exploit the sparsity and continuity of the SAD signature in the joint space and spatial frequency domain, a modified Bayesian compressive sensing algorithm is exploited to estimate the SAD signature. In this method, a spike-and-slab prior is used to statistically encourage sparsity of the SAD across each segmented SAD region, and a patterned prior is imposed to enforce the continuous structure of the SAD. The results are then mapped back to source range and DOA estimation for source localization. The effectiveness of the proposed technique is verified using simulation results with uniform and sparse linear arrays where the array sensors are located on a grid but with consecutive and missing positions.


---

## Body

## 1. Introduction

Localization of emitters is an important topic with many applications in radar, sonar, and wireless communications. The majority of research work in array signal processing has been focused on far-field signals; that is, the source signals are assumed to be located far from the receiver. Accordingly, the propagation waves assume a planar wavefront, making source localization solely based on the direction-of-arrival (DOA) of the sources. However, this characterization becomes invalid when sources are located in the near-field region. In this case, each source is characterized by both DOA and range that vary with sensor positions. To solve near-field source localization problems, subspace-based techniques, for example, MUSIC, were extended to a two-dimensional field [1]. However, this entails a high computational load required for both DOA and range searching.To decouple these two variables of the DOA and range for simplified computations, a more accurate quadratic approximation of the wavefront, in lieu of the linear approximation, can be made with respect to the sensor positions [2]. Consequently, the phase delay is no longer linear with the sensor positions but is decided by the DOA and range of the sources jointly. A number of methods have been proposed to jointly estimate the DOA and range information for near-field source localization. For instance, algorithms based on high-order statistics were proposed in [2–4]. The work in [2] exploits the rotational invariance in certain cumulant-domain signal subspace for range and bearing estimation. The technique in [3] applies a symmetric uniform linear array (ULA) and uses eigenvalues together with the corresponding eigenvectors of two properly designed matrices to jointly estimate signal parameters. As such, it avoids any spectral peak searching. Two orthogonally polarized sensors were utilized in [4] to improve performance. These methods, however, do not meet the desired performance and require certain array properties, for example, a ULA with an interelement spacing less than half a wavelength. In practical applications, such approaches are very costly due to receiver hardware and computational complexity.Sparse arrays [5–7] use fewer elements to achieve the same array aperture of fully populated arrays. They have similar mainbeam properties and, therefore, provide similar performance in terms of angular accuracy, resolution, and detection of targets close to interference directions, all being achieved with reduced size, weight, power consumption, and cost. An algorithm was proposed in [8] to cope with the near-field source localization problem using sparse arrays. However, it relies on the assumption that the sparse array configuration must be symmetric. In addition, similar to [2–4], the second-order Taylor approximation of the range between the source signals and the sensors is assumed to decouple the two variables. As a result, the performance is compromised when the sources are placed close to the array.In this paper, we propose an alternative approach using sparse arrays by mapping the source range and DOA information into sensor-dependent phase progression. The sensor-angle distribution (SAD) is used to represent the phase progression in the joint space and spatial frequency domain with a signature that resembles those of instantaneously narrowband frequency-modulated (FM) signals in the joint time-frequency domain [9, 10]. In this respect, the signal processing algorithms introduced for quadratic time-frequency analysis can be readily applied for the processing of the SAD for near-field source localization. It is important to note that, while the SAD was originally developed in the form of spatial Wigner distribution [9, 10], which is effective only for second-order (linear FM) spatial frequency signatures, we consider it in the general form within Cohen’s class with a proper time-frequency kernel for improved sensor-angle relationship [11]. The proposed technique permits the sensor-dependent phase progression as a high-order polynomial phase signal (PPS) to account for close sources, whereas the coefficients are estimated from the SAD signature with simple least squares fitting, which does not require parameter searching.Another important contribution of this paper is the consideration of the source localization problem in a sparse array configuration where the array sensors are located in a uniform linear grid but with missing positions. As such, conventional time-frequency analysis tools become infeasible due to the undesirable artifacts caused from missing sensor entries [12]. Instead, we use the emerging compressive sensing (CS) [13] techniques to reconstruct the SAD in the joint space and spatial frequency domain. The problem is considered in the Bayesian compressive sensing (BCS) framework that generally achieves desired performance and enables flexible utilization of various signature structures through the application of proper priors [14–17]. In particular, we use the structure-aware BCS, which was developed for effective reconstruction of sparse time-frequency distributions of FM signals, to exploit the sparsity and continuity of the spatial frequency signatures in the SAD domain. In doing so, the isolated or sporadic entries in the reconstructed SAD caused by the presence of missing sensors and measurement noise can be significantly suppressed, leading to more reliable SAD reconstruction and spatial frequency estimation. The obtained results are then mapped back to the estimation of source range and DOA information for effective source localization. Compared to the approach proposed in [8] using a sparse symmetric array, the proposed technique offers higher flexibility and robustness to array design and improves the source localization accuracy when the sources are close to the array and thus the phase exhibits a high-order polynomial phase relationship with the sensor positions.The rest of the paper is organized as follows. Section2 introduces the signal model of near-field localization based on a third-order approximation. The SAD and its sparse reconstruction are described in Section 3, and the structure-aware Bayesian compressive sensing technique applied for sparse SAD reconstruction and spatial frequency estimation is described in Section 4. Simulation results are provided in Section 5 to numerically evaluate and compare the performance. Such results reaffirm and demonstrate the effectiveness of the proposed technique. Section 6 concludes this paper.Notations. We use lower-case (upper-case) bold characters to denote vectors (matrices). In particular, I
N denotes the N
×
N identity matrix. (
·
)
∗ denotes complex conjugate, whereas (
·
)
T and (
·
)
H, respectively, denote the transpose and conjugate transpose of a matrix or vector. diag
(
x
) represents a diagonal matrix that uses the elements of x as its diagonal elements. ·
2
2 implies the Euclidean (l
2) norm of a vector. P
r
(
·
) expresses the probability density function (pdf) and p
(
x
∣
-
) denotes the conditional pdf of random variable x, given other parameters. Beta(x
,
y) computes the beta function for corresponding elements x and y. C
N
(
x
∣
a
,
b
) describes that random variable x follows a complex Gaussian distribution with mean a and variance b. Furthermore, δ
(
x
) is the Dirac delta function of x.


## 2. Signal Model

For the simplicity of presentation, the signal model is presented for the single-source case, but the proposed technique can be straightforwardly extended to the multiple source problems and simulation results are provided for two source scenarios.Consider a near-field narrowband source impinging on a uniform or sparse linear array withK
=
2
Q
+
1 elements, as depicted in Figure 1, where Q
∈
Z
+. Note that we assume K to be an odd integer without loss of generality, but the problem can be similarly formulated for an even value of K. Examples of even values of K are provided in Section 5. Let the index of the reference element be 0, and the sensors are located at p
k
d, where p
k
∈
Z for k
=
-
Q
,
…
,
Q, and d is the unit interelement spacing. To avoid spatial ambiguity, d is usually taken as half-wavelength, denoted as d
=
λ
/
2. The source is located with a range r and azimuth angle θ with respect to a reference sensor. Note that the array may be asymmetric; that is, it is not necessary to assume p
-
k
=
-
p
k, for k
=
-
Q
,
…
,
Q.Figure 1
The geometry for near-field source localization.In a near-field environment, the range between the source and the sensors varies for each of theK sensors. It is straightforward to show that the range between the source signal and sensor k is given by (1)
r
k
=
r
2
+
p
k
2
d
2
-
2
r
p
k
d
cos
⁡
(
θ
)
.As such, the near-field steering vector is given by (2)
a
=
exp
⁡
j
2
π
r
-
Q
λ
,
…
,
exp
⁡
j
2
π
r
Q
λ
T
.For source signal s
(
t
), the data received at the sensor array is expressed as (3)
x
t
=
a
s
t
+
n
t
,
t
=
1
,
…
,
T
,where n
(
t
) is the additive Gaussian noise vector and T is the number of observed time samples.To understand the significance of the steering vector, consider the special case of fully populated array that is uniformly spaced with interelement spacingd; that is, p
k
=
k. Then, (1) can be expressed as (4)
r
k
=
r
2
+
k
2
d
2
-
2
k
d
r
cos
⁡
(
θ
)
.Assume that r is reasonably larger than Q
d. Then, using the Taylor series expansion, we can obtain the third-order approximation, which is given by (5)
r
k
≈
r
-
k
d
cos
⁡
θ
+
k
2
d
2
2
r
sin
2
θ
+
k
3
d
3
2
r
2
cos
⁡
θ
sin
2
θ
.It is indicated that the near-field steering vector given in (2) is a uniformly sampled FM signal with respect to k. In this respect, the signal processing techniques devised for nonstationary signal analysis and quadratic time-frequency signal representations can be readily applied for near-field source localization in the joint space and spatial frequency domain. The corresponding spatial frequency is obtained from the phase difference with respect to k, given as (6)
f
(
k
)
=
-
1
2
λ
r
k
+
1
-
r
k
-
1
=
b
0
+
b
1
k
+
b
2
k
2
,where (7)
b
0
=
d
λ
cos
⁡
θ
-
d
3
2
r
2
λ
cos
⁡
θ
sin
2
θ
,
b
1
=
-
d
2
r
λ
sin
2
θ
,
b
2
=
-
3
d
3
2
r
2
λ
cos
⁡
θ
sin
2
θ
.It is clear that the polynomial coefficientsb
0, b
1, and b
2 are related to the range and the DOA. As such, these parameters can be estimated once the sensor-dependent spatial frequency is obtained. Note that, for a large range (r
≫
1), f
(
k
) behaves as a chirp when the second-order term in (6) is ignored. Furthermore, when the source is located in the far field, as r approaches infinity, f
(
k
) approaches a constant, implying a single spatial frequency associated with a sensor-independent DOA. According to (7), from the estimated values of b
0, b
1, and b
2, the unknown range and DOA can be, respectively, estimated as (8)
θ
^
=
arccos
λ
d
b
0
-
b
2
3
,
r
^
=
-
d
2
b
1
λ
sin
2
θ
^
.Specifically, when d
=
λ
/
2, the above expressions become (9)
θ
^
=
arccos
2
b
0
-
b
2
3
,
r
^
=
-
d
2
b
1
sin
2
(
θ
^
)
.


## 3. Sensor-Angle Distribution and Sparse Reconstruction

Several methods are available for the estimation of linear FM and PPS parameters (see, e.g., [18–20]). Such methods, however, are effective only when the samples are uniformly applied (i.e., the array is uniform linear in the underlying problem). When the array is sparse, as considered in this paper for more effective array aperture utilization, these methods become difficult to apply. For illustrative purposes, we consider a sparse linear array structure, consisting of a set of uniform linear subarrays with different coprime interelement spacings, as an example. The coprime array structure was proposed in [21] for systematical sparse array configuration, which generally achieves O
(
N
2
/
2
) coarray sensors with N physical sensors. Therefore, the aperture and the number of degrees-of-freedom (DOFs) can be improved using correlation-aware techniques. Coprime arrays have found broad applications in DOA estimation and interference suppression [22–25].For the spatial frequency estimation, we exploit the SAD to describe the relationship between spatial frequency and the sensor positions. The SAD draws parallel to the time-frequency distribution. In essence, the SAD uses sensor positions in lieu of time samples whereas the spatial frequency replaces temporal frequency. With this analogy, SAD becomes a tool to reveal the local behavior of the spatial frequency along the space axis with respect to the array sensor positions. For sparse arrays, the SAD suffers from missing artifacts due to missing samples. We use time-frequency kernel and sparse reconstruction to achieve satisfactory SAD reconstruction and spatial frequency estimation. In the following, we describe the SAD definition in Section3.1 and discuss the effect of array sparsity in Section 3.2. The sparse reconstruction of SAD from a sparse array is addressed in Section 3.3, and the structure-aware Bayesian compressive sensing that implements the sparse reconstruction is detailed in Section 4.


### 3.1. Sensor-Angle Distribution

The SAD corresponding to theith sensor position and spatial frequency ξ at time instant t is expressed as (10)
D
~
i
,
ξ
;
t
=
∑
l
=
-
∞
∞
‍
∑
v
=
-
∞
∞
‍
ϕ
v
,
l
x
i
+
v
+
l
;
t
·
x
∗
(
i
+
v
-
l
;
t
)
e
-
j
4
π
ξ
l
,where x
(
i
;
t
) denotes the ith element of vector x
(
t
) and ϕ
(
v
,
l
) is a kernel function that characterizes the distribution and is a function of sensor position and sensor lag [26]. The above equation defines Cohen’s class of power distribution over the joint variables i and ξ. In this respect, all the standard kernel designs applied in the time-frequency analysis literature can be employed with the SAD. The distribution in (10) can be averaged over the T samples to reduce the effect of observation noise, yielding (11)
D
i
,
ξ
=
1
T
∑
t
=
1
T
‍
D
~
i
,
ξ
;
t
.Stacking D
(
i
,
ξ
) with respect to i and ξ forms an SAD matrix D.


### 3.2. Data-Dependent Time-Frequency Kernel

Benefiting from the analogy of SAD to time-frequency distributions, we understand that the sparse placement of array sensors is deemed to produce artifacts that clutter the SAD, making the spatial frequency estimation difficult. Analysis of the effects of missing samples on time-frequency distribution is given in [12]. To mitigate effects of the induced artifacts and improve time-frequency signature estimation, the first step is to apply a time-frequency kernel [12, 17, 27]. Time-frequency kernel allows mitigation of cross-terms between different signal components as well as the effect of artifacts due to missing samples (dual of missing sensors in the underlying problem).Time-frequency kernels are generally represented as a multiplying weight (mask) matrix in the ambiguity domain. The ambiguity function corresponding to the underlying SAD is defined with respect to sensor position lagl and spatial frequency shift p as (12)
Q
~
p
,
l
;
t
=
∑
i
=
-
∞
∞
‍
∑
v
=
-
∞
∞
‍
x
i
+
l
;
t
x
∗
i
-
l
;
t
e
-
j
4
π
p
i
.The averaged AF is then obtained by averaging over the T time samples as (13)
Q
p
,
l
;
t
=
1
T
∑
t
=
1
T
‍
Q
~
p
,
l
;
t
.The time-frequency kernels are based on the fact that majority of autoterm distributions are located near the origin in the ambiguity domain, whereas the cross-terms and artifacts are distant from the time lag and frequency-shift (sensor lag and spatial frequency shift in the underlying problem) axes. This property has motivated researchers to propose educed-interference distribution kernels with low-pass filter characteristics to suppress cross-terms and preserve autoterms. Such kernels can be signal-independent or signal-dependent. In this paper, we use the adaptive optimal kernel (AOK) and modify it to specifically deal with missing positions. The AOK is formulated as an optimization problem in [28]: (14)
max
⁡
C
∑
r
‍
∑
ϕ
‍
r
Q
r
,
ϕ
C
r
,
ϕ
2
subject
to
C
(
r
,
ϕ
)
=
e
-
r
2
/
2
σ
2
(
ϕ
)
,
1
2
π
∑
‍
σ
2
(
ϕ
)
≤
α
,where C denotes a matrix with all C
(
r
,
ϕ
) entries. Note that the ambiguity function Q
(
r
,
ϕ
) and the kernel C
(
r
,
ϕ
) are defined in the polar coordinates with the radius r and angle ϕ. The objective of the AOK is to maximize the energy in the ambiguity domain in order to preserve the autoterms. At the same time, the Gaussian constraint in C
(
r
,
ϕ
) leads to cross-term reduction. The choice of ϕ influences the tradeoff between cross-term suppression and autoterm preservation.


### 3.3. Sparse Reconstruction of SAD from Sparse Array

As addressed above, a key step to achieve improved SAD reconstruction and spatial frequency estimation in the presence of missing sensors is through sparse reconstruction. In extending the techniques developed for sparse reconstruction of time-frequency distribution in [12, 17, 27] to the underlying sparse SAD reconstruction, we define the following discrete sensor autocorrelation function (SAF) for time instant t: (15)
A
~
i
,
l
;
t
=
∑
v
=
-
∞
∞
‍
ϕ
v
,
l
x
i
+
v
+
l
;
t
x
∗
i
+
v
-
l
;
t
.The averaged SAF (16)
A
(
i
,
l
)
=
1
T
∑
t
=
1
T
‍
A
~
(
i
,
l
;
t
)can be similarly defined using all available samples to reduce the effect of additive noise. The SAF matrix A is formed that contains A
(
i
,
l
) for all the entries of i and l. As such, at each sensor position, the columns of the averaged SAD matrix D and that of the averaged SAF matrix A are related by a one-dimensional (1D) Fourier transform along the lag dimension. Denote a
[
i
] as a column vector that contains all the elements in the ith column of A and d
[
i
] as a vector corresponding to the ith column of D. Then these two vectors are related by (17)
a
[
i
]
=
F
d
[
i
]
,
i
=
1
,
…
,
K
,where F is the inverse Fourier transform matrix. Equation (17) describes a linear model relating the sparse vector d
[
i
] to the compressed observation vector a
[
i
] using the inverse Fourier dictionary. As such, d
[
i
] can be obtained by solving a sparse optimization problem, which is commonly used in compressive sensing techniques. The positions of nonzero entries in d
[
i
] indicate the estimations of spatial frequency. Note that F can be designed as a wide matrix so as to achieve a higher spatial frequency resolution in the estimated d
[
i
].As demonstrated for time-frequency analysis, the 1D relationship between the instantaneous autocorrelation function (IAF, the dual of SAF) and the time-frequency distribution (the dual of SAD) is much more effective than the conventional two-dimensional (2D) approach between the ambiguity function and the time-frequency distribution [29, 30] and allows local time-frequency reconstruction for each time instant (the dual of sensor positions). In particular, such local reconstruction capability is fundamental to enable the development of the structure-aware Bayesian compressive sensing technique to be described in Section 4.


## 3.1. Sensor-Angle Distribution

The SAD corresponding to theith sensor position and spatial frequency ξ at time instant t is expressed as (10)
D
~
i
,
ξ
;
t
=
∑
l
=
-
∞
∞
‍
∑
v
=
-
∞
∞
‍
ϕ
v
,
l
x
i
+
v
+
l
;
t
·
x
∗
(
i
+
v
-
l
;
t
)
e
-
j
4
π
ξ
l
,where x
(
i
;
t
) denotes the ith element of vector x
(
t
) and ϕ
(
v
,
l
) is a kernel function that characterizes the distribution and is a function of sensor position and sensor lag [26]. The above equation defines Cohen’s class of power distribution over the joint variables i and ξ. In this respect, all the standard kernel designs applied in the time-frequency analysis literature can be employed with the SAD. The distribution in (10) can be averaged over the T samples to reduce the effect of observation noise, yielding (11)
D
i
,
ξ
=
1
T
∑
t
=
1
T
‍
D
~
i
,
ξ
;
t
.Stacking D
(
i
,
ξ
) with respect to i and ξ forms an SAD matrix D.


## 3.2. Data-Dependent Time-Frequency Kernel

Benefiting from the analogy of SAD to time-frequency distributions, we understand that the sparse placement of array sensors is deemed to produce artifacts that clutter the SAD, making the spatial frequency estimation difficult. Analysis of the effects of missing samples on time-frequency distribution is given in [12]. To mitigate effects of the induced artifacts and improve time-frequency signature estimation, the first step is to apply a time-frequency kernel [12, 17, 27]. Time-frequency kernel allows mitigation of cross-terms between different signal components as well as the effect of artifacts due to missing samples (dual of missing sensors in the underlying problem).Time-frequency kernels are generally represented as a multiplying weight (mask) matrix in the ambiguity domain. The ambiguity function corresponding to the underlying SAD is defined with respect to sensor position lagl and spatial frequency shift p as (12)
Q
~
p
,
l
;
t
=
∑
i
=
-
∞
∞
‍
∑
v
=
-
∞
∞
‍
x
i
+
l
;
t
x
∗
i
-
l
;
t
e
-
j
4
π
p
i
.The averaged AF is then obtained by averaging over the T time samples as (13)
Q
p
,
l
;
t
=
1
T
∑
t
=
1
T
‍
Q
~
p
,
l
;
t
.The time-frequency kernels are based on the fact that majority of autoterm distributions are located near the origin in the ambiguity domain, whereas the cross-terms and artifacts are distant from the time lag and frequency-shift (sensor lag and spatial frequency shift in the underlying problem) axes. This property has motivated researchers to propose educed-interference distribution kernels with low-pass filter characteristics to suppress cross-terms and preserve autoterms. Such kernels can be signal-independent or signal-dependent. In this paper, we use the adaptive optimal kernel (AOK) and modify it to specifically deal with missing positions. The AOK is formulated as an optimization problem in [28]: (14)
max
⁡
C
∑
r
‍
∑
ϕ
‍
r
Q
r
,
ϕ
C
r
,
ϕ
2
subject
to
C
(
r
,
ϕ
)
=
e
-
r
2
/
2
σ
2
(
ϕ
)
,
1
2
π
∑
‍
σ
2
(
ϕ
)
≤
α
,where C denotes a matrix with all C
(
r
,
ϕ
) entries. Note that the ambiguity function Q
(
r
,
ϕ
) and the kernel C
(
r
,
ϕ
) are defined in the polar coordinates with the radius r and angle ϕ. The objective of the AOK is to maximize the energy in the ambiguity domain in order to preserve the autoterms. At the same time, the Gaussian constraint in C
(
r
,
ϕ
) leads to cross-term reduction. The choice of ϕ influences the tradeoff between cross-term suppression and autoterm preservation.


## 3.3. Sparse Reconstruction of SAD from Sparse Array

As addressed above, a key step to achieve improved SAD reconstruction and spatial frequency estimation in the presence of missing sensors is through sparse reconstruction. In extending the techniques developed for sparse reconstruction of time-frequency distribution in [12, 17, 27] to the underlying sparse SAD reconstruction, we define the following discrete sensor autocorrelation function (SAF) for time instant t: (15)
A
~
i
,
l
;
t
=
∑
v
=
-
∞
∞
‍
ϕ
v
,
l
x
i
+
v
+
l
;
t
x
∗
i
+
v
-
l
;
t
.The averaged SAF (16)
A
(
i
,
l
)
=
1
T
∑
t
=
1
T
‍
A
~
(
i
,
l
;
t
)can be similarly defined using all available samples to reduce the effect of additive noise. The SAF matrix A is formed that contains A
(
i
,
l
) for all the entries of i and l. As such, at each sensor position, the columns of the averaged SAD matrix D and that of the averaged SAF matrix A are related by a one-dimensional (1D) Fourier transform along the lag dimension. Denote a
[
i
] as a column vector that contains all the elements in the ith column of A and d
[
i
] as a vector corresponding to the ith column of D. Then these two vectors are related by (17)
a
[
i
]
=
F
d
[
i
]
,
i
=
1
,
…
,
K
,where F is the inverse Fourier transform matrix. Equation (17) describes a linear model relating the sparse vector d
[
i
] to the compressed observation vector a
[
i
] using the inverse Fourier dictionary. As such, d
[
i
] can be obtained by solving a sparse optimization problem, which is commonly used in compressive sensing techniques. The positions of nonzero entries in d
[
i
] indicate the estimations of spatial frequency. Note that F can be designed as a wide matrix so as to achieve a higher spatial frequency resolution in the estimated d
[
i
].As demonstrated for time-frequency analysis, the 1D relationship between the instantaneous autocorrelation function (IAF, the dual of SAF) and the time-frequency distribution (the dual of SAD) is much more effective than the conventional two-dimensional (2D) approach between the ambiguity function and the time-frequency distribution [29, 30] and allows local time-frequency reconstruction for each time instant (the dual of sensor positions). In particular, such local reconstruction capability is fundamental to enable the development of the structure-aware Bayesian compressive sensing technique to be described in Section 4.


## 4. Structure-Aware Bayesian Compressive Sensing

### 4.1. Continuous Structure Modeling with Spike-and-Slab Prior

A number of approaches, such as the OMP and LASSO [31, 32], can be effectively applied for the reconstruction of sparse SAD signatures. To exploit the continuous structure of the SAD signature, however, the BCS-based approaches are considered most effective due to their flexibility of employing proper priors for this purpose. We use the structure-aware BCS technique [17] to improve estimation performance.To encourage sparsity of the SAD signatures, we place a spike-and-slab prior to each element ofd
[
i
] [33]; that is, for its lth element, (18)
p
(
d
l
i
∣
π
l
i
,
β
0
i
)
=
1
-
π
l
i
δ
d
l
i
+
π
l
i
C
N
(
d
l
i
∣
0
,
β
l
i
-
1
)
,where π
l
[
i
] is the prior probability of a nonzero element and β
l
[
i
] is the precision (reciprocal of variance) of the complex Gaussian distribution, for the lth sensor position.To facilitate an analytical inference, which involves the delta function, the product of following two latent variables is introduced to follow the pdf in (18): (19)
d
l
[
i
]
=
z
l
[
i
]
η
l
[
i
]
,where (20)
η
l
[
i
]
~
C
N
η
l
[
i
]
∣
0
,
β
l
i
-
1
,
z
l
[
i
]
~
π
l
[
i
]
z
l
[
i
]
1
-
π
l
[
i
]
1
-
z
l
[
i
]is a binary variable following the Bernoulli distribution [17]. In addition, a Gaussian prior is placed on the additive noise as ϵ
l
[
i
]
~
C
N
(
ϵ
l
[
i
]
∣
0
,
α
-
1
), where the precision α is assumed to be independent of i and l. To acquire a trackable posterior distribution, we place Gamma priors, which are conjugate to the Gaussian distribution, on both β
l
[
i
] and α.To encourage continuity in the joint space and spatial frequency domain, we utilize data in the neighborhood sensor positions betweeni
-
1 and i
+
1, when the SAD signature in sensor position i is estimated. As shown in Figure 2, we categorize the relationship into three different patterns, termed asPattern 1 (“strong rejection”),Pattern 2 (“weak rejection”), andPattern 3 (“strong acceptance”), respectively. Define z
~
1
,
…
,
z
~
6 as illustrated in Figure 2. Based on the idea of continuity, we define two variables, κ
l
[
i
]
=
∑
n
=
1
6
‍
z
~
n and χ
l
[
i
]
=
∑
n
=
1
3
‍
z
~
n
z
~
7
-
n, which, respectively, represent the number of nonzero neighboring entries and the number of nonzero diagonal pairs. It should be noted that the vertical pixel pairs are excluded for being impractical in the underlying problem.Figure 2
Three patterns for 2D image. White squares denote entries with zeros amplitudes, blue shaded squares denote the entry with nonzero amplitudes, and red shaded squares denote the entry under test. (a) Pattern 1: strong rejection; (b) Pattern 2: examples for weak rejection withκ
l
[
i
]
=
1 case; (c) Pattern 3: examples for strong acceptance with χ
l
[
i
]
=
1 case.
(a)
(b)
(c)Pattern 1 includes two cases ofκ
l
[
i
]
=
0 or χ
l
[
i
]
=
3. In the former case, all neighbors are zero values, and thus the lth entry would be zero valued with a high probability. In the latter cases, all the neighbors are nonzero values. In both cases, e
0
<
f
0 is assumed in prior Beta(e
0,f
0) to encourage a small value of π
l
[
i
] so as to reject this entry.Pattern 2 illustrates those cases ofκ
l
[
i
]
>
0 or χ
l
[
i
]
≠
1. Figure 2(b) shows one of the examples with a single nonzero neighboring entry. In this case, the probability that the lth entry takes zero values is fair, and e
1
=
f
1 is used in the prior Beta(e
1,f
1) to exert noninformative prior on π
l
[
i
].Pattern 3 depicts those cases withχ
l
[
i
]
=
1. Figure 2(c) shows two examples of such cases. We confidently believe that the lth entry would take a nonzero value with a high probability when nonzero neighboring entries have such continuous and symmetric structure in the 2D space and spatial frequency domain. In this case, e
2
>
f
2 in the prior Beta(e
2,f
2) is assumed to encourage a large value of π
l
[
i
] to accept this entry.Note that the above patterns encourage instantaneously narrowband spatial frequency signatures and closely follow those developed for time-frequency signature reconstruction of FM signals [17]. These patterns are, however, different from those used to encourage merely spatial connectivity as described in [34, 35].


### 4.2. Bayesian Inference

The Bayesian inference of the proposed algorithm is carried out by the Gibbs sampler [33]. The paired Gibbs sampler iteratively samples the following conditional pdf: (21)
p
z
l
[
i
]
,
η
l
[
i
]
∣
η
∖
l
[
i
]
,
z
∖
l
[
i
]
,
a
[
i
]
=
p
η
l
i
∣
z
l
[
i
]
,
η
∖
l
[
i
]
,
z
∖
l
[
i
]
,
a
[
i
]
×
p
z
l
[
i
]
∣
η
∖
l
[
i
]
,
z
∖
l
[
i
]
,
a
[
i
]
,where η
∖
l
[
i
] and z
∖
l
[
t
], respectively, denote η
l
[
i
] except the variable η
l
[
i
] and z
[
i
] except the variable z
l
[
i
].The probabilityp
z
l
[
i
]
=
1
∣
η
∖
l
[
i
]
,
z
∖
l
[
i
]
,
a
[
i
] is acquired analytically as (22)
p
z
l
[
i
]
=
1
∣
η
∖
l
[
i
]
,
z
∖
l
[
i
]
,
a
[
i
]
=
β
l
[
i
]
1
-
β
l
[
i
]
C
N
0
,
β
l
i
-
1
C
N
u
l
[
i
]
,
σ
l
[
i
]
,where u
l
[
i
] is derived as (23)
u
l
[
i
]
=
σ
l
i
-
1
α
F
l
H
a
∖
l
[
i
]
,with σ
l
[
i
]
=
α
F
l
H
F
l
+
β
l
[
i
]
-
1, a
∖
l
[
i
]
=
a
[
i
]
-
∑
k
≠
l
‍
F
l
z
k
[
i
]
η
k
[
i
], and F
l represents the lth column of the inverse Fourier transform matrix F. The conditional distribution of p
(
η
l
[
i
]
∣
z
l
[
i
]
=
1
,
η
∖
l
[
i
]
,
z
∖
l
[
i
]
,
a
[
i
]
) is expressed as (24)
p
η
l
[
i
]
∣
z
l
[
i
]
=
1
,
η
∖
l
[
i
]
,
z
∖
l
[
i
]
,
a
[
i
]
=
C
N
(
η
l
[
i
]
∣
u
l
[
i
]
,
σ
l
[
i
]
)
.For z
l
[
i
]
=
0, the value of η
l
[
i
] does not affect the result of d
l
[
i
]. Thus, we can draw the value of variable η
l
[
i
] from its prior. We adopt the procedure in [17] to update the mixing weight π
l
[
i
], signal precision β
l
[
i
], and noise precision α, in order to obtain z
l
[
i
] and η
l
[
i
]. Then, the sparse vector d
[
i
], whose positions of nonzero entries indicate the estimated nonzero spatial frequencies f
(
i
) in f
[
i
], can be obtained using (19).Based on (6), the following linear model can be obtained: (25)
f
=
B
b
,where b
=
[
b
0
,
b
1
,
b
2
]
T, f
=
[
f
(
1
)
,
f
(
2
)
,
…
,
f
K
′
]
T, and (26)
B
=
1
1
1
2
1
2
2
2
⋮
⋮
⋮
1
K
K
2
.As such, once the sparse spatial frequencies along the array axis are estimated, the polynomial phase coefficients inb can be obtained using least squares fitting as (27)
b
^
=
B
T
B
-
1
B
T
f
.Then, the range and DOA can be estimated based on (8) or (9).


## 4.1. Continuous Structure Modeling with Spike-and-Slab Prior

A number of approaches, such as the OMP and LASSO [31, 32], can be effectively applied for the reconstruction of sparse SAD signatures. To exploit the continuous structure of the SAD signature, however, the BCS-based approaches are considered most effective due to their flexibility of employing proper priors for this purpose. We use the structure-aware BCS technique [17] to improve estimation performance.To encourage sparsity of the SAD signatures, we place a spike-and-slab prior to each element ofd
[
i
] [33]; that is, for its lth element, (18)
p
(
d
l
i
∣
π
l
i
,
β
0
i
)
=
1
-
π
l
i
δ
d
l
i
+
π
l
i
C
N
(
d
l
i
∣
0
,
β
l
i
-
1
)
,where π
l
[
i
] is the prior probability of a nonzero element and β
l
[
i
] is the precision (reciprocal of variance) of the complex Gaussian distribution, for the lth sensor position.To facilitate an analytical inference, which involves the delta function, the product of following two latent variables is introduced to follow the pdf in (18): (19)
d
l
[
i
]
=
z
l
[
i
]
η
l
[
i
]
,where (20)
η
l
[
i
]
~
C
N
η
l
[
i
]
∣
0
,
β
l
i
-
1
,
z
l
[
i
]
~
π
l
[
i
]
z
l
[
i
]
1
-
π
l
[
i
]
1
-
z
l
[
i
]is a binary variable following the Bernoulli distribution [17]. In addition, a Gaussian prior is placed on the additive noise as ϵ
l
[
i
]
~
C
N
(
ϵ
l
[
i
]
∣
0
,
α
-
1
), where the precision α is assumed to be independent of i and l. To acquire a trackable posterior distribution, we place Gamma priors, which are conjugate to the Gaussian distribution, on both β
l
[
i
] and α.To encourage continuity in the joint space and spatial frequency domain, we utilize data in the neighborhood sensor positions betweeni
-
1 and i
+
1, when the SAD signature in sensor position i is estimated. As shown in Figure 2, we categorize the relationship into three different patterns, termed asPattern 1 (“strong rejection”),Pattern 2 (“weak rejection”), andPattern 3 (“strong acceptance”), respectively. Define z
~
1
,
…
,
z
~
6 as illustrated in Figure 2. Based on the idea of continuity, we define two variables, κ
l
[
i
]
=
∑
n
=
1
6
‍
z
~
n and χ
l
[
i
]
=
∑
n
=
1
3
‍
z
~
n
z
~
7
-
n, which, respectively, represent the number of nonzero neighboring entries and the number of nonzero diagonal pairs. It should be noted that the vertical pixel pairs are excluded for being impractical in the underlying problem.Figure 2
Three patterns for 2D image. White squares denote entries with zeros amplitudes, blue shaded squares denote the entry with nonzero amplitudes, and red shaded squares denote the entry under test. (a) Pattern 1: strong rejection; (b) Pattern 2: examples for weak rejection withκ
l
[
i
]
=
1 case; (c) Pattern 3: examples for strong acceptance with χ
l
[
i
]
=
1 case.
(a)
(b)
(c)Pattern 1 includes two cases ofκ
l
[
i
]
=
0 or χ
l
[
i
]
=
3. In the former case, all neighbors are zero values, and thus the lth entry would be zero valued with a high probability. In the latter cases, all the neighbors are nonzero values. In both cases, e
0
<
f
0 is assumed in prior Beta(e
0,f
0) to encourage a small value of π
l
[
i
] so as to reject this entry.Pattern 2 illustrates those cases ofκ
l
[
i
]
>
0 or χ
l
[
i
]
≠
1. Figure 2(b) shows one of the examples with a single nonzero neighboring entry. In this case, the probability that the lth entry takes zero values is fair, and e
1
=
f
1 is used in the prior Beta(e
1,f
1) to exert noninformative prior on π
l
[
i
].Pattern 3 depicts those cases withχ
l
[
i
]
=
1. Figure 2(c) shows two examples of such cases. We confidently believe that the lth entry would take a nonzero value with a high probability when nonzero neighboring entries have such continuous and symmetric structure in the 2D space and spatial frequency domain. In this case, e
2
>
f
2 in the prior Beta(e
2,f
2) is assumed to encourage a large value of π
l
[
i
] to accept this entry.Note that the above patterns encourage instantaneously narrowband spatial frequency signatures and closely follow those developed for time-frequency signature reconstruction of FM signals [17]. These patterns are, however, different from those used to encourage merely spatial connectivity as described in [34, 35].


## 4.2. Bayesian Inference

The Bayesian inference of the proposed algorithm is carried out by the Gibbs sampler [33]. The paired Gibbs sampler iteratively samples the following conditional pdf: (21)
p
z
l
[
i
]
,
η
l
[
i
]
∣
η
∖
l
[
i
]
,
z
∖
l
[
i
]
,
a
[
i
]
=
p
η
l
i
∣
z
l
[
i
]
,
η
∖
l
[
i
]
,
z
∖
l
[
i
]
,
a
[
i
]
×
p
z
l
[
i
]
∣
η
∖
l
[
i
]
,
z
∖
l
[
i
]
,
a
[
i
]
,where η
∖
l
[
i
] and z
∖
l
[
t
], respectively, denote η
l
[
i
] except the variable η
l
[
i
] and z
[
i
] except the variable z
l
[
i
].The probabilityp
z
l
[
i
]
=
1
∣
η
∖
l
[
i
]
,
z
∖
l
[
i
]
,
a
[
i
] is acquired analytically as (22)
p
z
l
[
i
]
=
1
∣
η
∖
l
[
i
]
,
z
∖
l
[
i
]
,
a
[
i
]
=
β
l
[
i
]
1
-
β
l
[
i
]
C
N
0
,
β
l
i
-
1
C
N
u
l
[
i
]
,
σ
l
[
i
]
,where u
l
[
i
] is derived as (23)
u
l
[
i
]
=
σ
l
i
-
1
α
F
l
H
a
∖
l
[
i
]
,with σ
l
[
i
]
=
α
F
l
H
F
l
+
β
l
[
i
]
-
1, a
∖
l
[
i
]
=
a
[
i
]
-
∑
k
≠
l
‍
F
l
z
k
[
i
]
η
k
[
i
], and F
l represents the lth column of the inverse Fourier transform matrix F. The conditional distribution of p
(
η
l
[
i
]
∣
z
l
[
i
]
=
1
,
η
∖
l
[
i
]
,
z
∖
l
[
i
]
,
a
[
i
]
) is expressed as (24)
p
η
l
[
i
]
∣
z
l
[
i
]
=
1
,
η
∖
l
[
i
]
,
z
∖
l
[
i
]
,
a
[
i
]
=
C
N
(
η
l
[
i
]
∣
u
l
[
i
]
,
σ
l
[
i
]
)
.For z
l
[
i
]
=
0, the value of η
l
[
i
] does not affect the result of d
l
[
i
]. Thus, we can draw the value of variable η
l
[
i
] from its prior. We adopt the procedure in [17] to update the mixing weight π
l
[
i
], signal precision β
l
[
i
], and noise precision α, in order to obtain z
l
[
i
] and η
l
[
i
]. Then, the sparse vector d
[
i
], whose positions of nonzero entries indicate the estimated nonzero spatial frequencies f
(
i
) in f
[
i
], can be obtained using (19).Based on (6), the following linear model can be obtained: (25)
f
=
B
b
,where b
=
[
b
0
,
b
1
,
b
2
]
T, f
=
[
f
(
1
)
,
f
(
2
)
,
…
,
f
K
′
]
T, and (26)
B
=
1
1
1
2
1
2
2
2
⋮
⋮
⋮
1
K
K
2
.As such, once the sparse spatial frequencies along the array axis are estimated, the polynomial phase coefficients inb can be obtained using least squares fitting as (27)
b
^
=
B
T
B
-
1
B
T
f
.Then, the range and DOA can be estimated based on (8) or (9).


## 5. Simulation Results

Simulation examples are presented to demonstrate the effectiveness of the proposed technique. For each sensor, the maximum number of iterations in the Gibbs sampling is 200, and the sampler with the maximum marginal likelihood in the last 20 samples is chosen to estimated
[
i
] for i
=
-
Q
,
…
,
Q. For all examples, the carrier frequency is considered to be 15 MHz, yielding a wavelength of λ
=
20 m. As such, the unit interelement spacing is set to d
=
λ
/
2
=
10 m. The source signals are assumed to be quadrature phase-shift keying (QPSK) modulated. In all simulations, we use a 32-element ULA as well as its subset described as a linear sparse array with missing sensors. For the latter, as shown in Figure 3, we consider a sparse array configuration consisting of three uniform linear subarrays with coprime interelement spacings of M
1
=
5, M
2
=
6, and M
3
=
7. The subarray sensors are, respectively, located at [
-
15
:
5
:
15
]
d, [
-
16
:
6
:
15
]
d, and [
-
15
:
7
:
15
]
d. Compared to the fully populated 32-sensor ULA, the sparse array consists only of 16 sensors, resulting in 50% of sensor reduction. The positions of the missing sensors are indicated by × in Figure 3. It should be noted that this configuration is not symmetric and, as such, the approach proposed in [8] becomes inapplicable.Figure 3
The sparse array configuration.
(a)
Subarray 1 with interelement spacing5
d
(b)
Subarray 2 with interelement spacing6
d
(c)
Subarray 3 with interelement spacing7
d
(d)
Coprime array configurationWe evaluate the source localization performance through Monte Carlo simulations. The average root mean square error (RMSE) of the estimated source location is used as the performance metrics, expressed as(28)
RMSE
=
1
M
∑
m
=
1
M
‍
Δ
x
m
2
+
Δ
y
m
2
,where (29)
Δ
x
(
m
)
=
r
^
m
cos
⁡
θ
^
m
-
r
cos
⁡
(
θ
)
,
Δ
y
(
m
)
=
r
^
m
sin
θ
^
m
-
r
sin
(
θ
)
,and θ
^
(
m
) and r
^
(
m
) are the estimates of θ and r, respectively, for the mth Monte Carlo trial, m
=
1
,
…
,
M. M
=
100 independent trials are used in all simulations.In the sequel, three examples are examined. The first and second examples consider the localization of a single near-field source but with different ranges, whereas the third example considers a two-source scenario.


### 5.1. Example I: Single Source with a Moderate Range

The first example considers the localization of a single near-field source with a range of 1000 m and DOA of 30° with respect to the reference sensor. The corresponding source coordinate is[
866.03
,
500
] m. In this case, the spatial frequency is approximately linear to the sensor position. The input signal-to-noise ratio (SNR) is 0 dB and 100 time samples are observed. In order to obtain a smooth figure, we oversample the spatial frequency on a multiple of 128. Figure 4 shows the SAD results for a ULA from both Wigner-Ville distribution (WVD) and sparse reconstruction, and the corresponding results for the sparse array are illustrated in Figure 5.Figure 4
Near-field source location estimation using the ULA (r
1
=
1000 m, θ
1
=
3
0
∘).
(a)
WVD
(b)
Sparse reconstructionFigure 5
Near-field source location estimation using the sparse array (r
1
=
1000 m, θ
1
=
3
0
°).
(a)
WVD
(b)
Sparse reconstruction
(c)
AOK distribution
(d)
Sparse reconstructionFigure4(a) shows the WVD of the SAD. Because there is no missing sensor, the signature can be estimated using a number of conventional methods as well as sparse reconstruction techniques, with respect to the SAF. The corresponding results obtained from the OMP and the structure-aware BCS are depicted in Figure 4(b). It is apparent that the signature is well approximated as a chirp and there is little difference noticed between the quadratic and third-order approximations. As such, in this case, the use of third-order approximation with the proposed BCS technique achieves similar performance. The PPS coefficients estimated from the proposed method are, respectively, b
=
[
0.4321
,
-
0.0013
,
-
3
×
1
0
-
6
]
T, and the yielding source position estimate is [
844.81
,
492.07
] m, which corresponds to a 22.6 m error to the true source position.The corresponding 16-element sparse array scenario is presented in Figure5, in which the conventional method developed in [8] is not effective. In addition, the SAD is contaminated with artifacts that may obscure the true feature estimation due to missing sensors, as shown in Figure 5(a). There are erroneous spatial frequency estimates for sparse reconstruction depicted in Figure 5(b). For the SAD technique, the modified AOK can be applied to suppress the effect of missing sensors. The results are shown in Figures 5(c) and 5(d). It is evident that a desired performance is achieved.Figure6 compares the RMSE performance as a function of the input SNR obtained from 100 independent trials. It is evident that the source localization performance is improved as the input SNR increases. The WVD in the ULA case presents the best performance since the WVD is known to be optimal in the analysis of linear FM signals [36, 37]. However, the performance of WVD degrades when the spatial frequency is nonlinear to the sensors positions or when some sensors are missing. The sparse array with the proposed method yields a similar performance but suffers from a gap of performance degradation. Such performance degradation partly stems from the signal power loss with the fewer number of sensors.Figure 6
RMSE versus SNR (r
1
=
1000 m, θ
1
=
3
0
∘).


### 5.2. Example II: Single-Source with a Close Range

In this example, the source is located at a range of 350 m and DOA of3
0
° with respect to the reference sensor. The corresponding coordinate is [
303.1
,
175
] m. As a consequence, the signature behaves as a PPS with respect to the sensor positions. The WVD and the estimated results obtained from the ULA and sparse array are shown in Figures 7 and 8, respectively. As indicated in Figure 7, the estimated spatial frequency is approximated as a third-order PPS with respect to the sensor positions. In this case, the WVD exhibits some cross-terms. The third-order PPS model yields more accurate results as compared to the quadratic approximation in this case to reflect the closer source range to the array. Under the third-order PPS approximation, the coefficient b
=
[
0.4277
,
-
0.004
,
-
0.0001
]
T is obtained, resulting in the source position [
286.14
,
173.16
] m, a 17.06 m shift from the true source position. For the sparse array, the effect of artifacts is again confirmed in Figures 8(a) and 8(b). This effect is significantly suppressed after the AOK is applied, as shown in Figures 8(c) and 8(d).Figure 7
Near-field source location estimation using the ULA (r
1
=
350 m, θ
1
=
3
0
∘).
(a)
WVD
(b)
Sparse reconstructionFigure 8
Near-field source location estimation using the sparse array (r
1
=
350 m, θ
1
=
3
0
∘).
(a)
WVD
(b)
Sparse reconstruction
(c)
AOK processing
(d)
Sparse reconstructionIn Figure9, we compare the RMSE performance of the second example. It is shown that the AOK-based results achieve better performance because of the effective cross-term and artifact suppression. In addition, there appears a floor in the obtained RMSE performance due to the spatial frequency resolution and bias. For the sparse array, a similar performance is attained to that of the ULA when the input SNR is higher than 0 dB and outperforms the WVD-based ULA results.Figure 9
RMSE versus SNR (r
1
=
350 m, θ
1
=
3
0
°).


### 5.3. Example III: Two-Source Scenario

Finally, we consider a two-source case in Figure10. The first source is located at r
=
300 m and θ
=
3
0
° with respect to the reference sensor, whereas the second source is located at r
=
1000 m and θ
=
6
0
∘. For the 32-element ULA, as shown in Figure 10(a), the SAD in this case is contaminated by the cross-terms between the two sources, making the spatial frequency estimation difficult. Such cross-term effect can be substantially reduced by exploiting the AOK, as shown in Figure 10(b). The sparse reconstruction results of OMP obtained from the AOK are depicted in Figure 10(c), where the spatial frequencies are not correctly identified towards both edges. Improved results are obtained from the proposed technique exploiting the structure-aware BCS, as depicted in Figure 10(d), due to the enforced spatial frequency continuity in the sparse SAD reconstruction. The latter evidently shows better performance for the multiple near-field source case. When the 16-element spare array is used, as shown in Figure 11(a), the SAD suffers from both cross-terms and artifacts due to missing sensor positions, yielding highly obscured SAD signatures. Both cross-terms and artifacts are significantly reduced when the AOK is used, as depicted in Figure 11(b). When it is combined with the proposed structure-aware BCS for sparse spatial frequency reconstruction, improved spatial frequency estimation, as shown in Figure 11(d), is achieved. For comparison, the OMP fails to identify one of the sources in the reconstructed spatial frequencies, as depicted in Figure 11(c).Figure 10
Near-field source location estimation using the ULA (r
1
=
300 m, θ
1
=
3
0
° and r
2
=
1000 m, θ
2
=
6
0
°).
(a)
WVD
(b)
AOK distribution
(c)
OMP
(d)
The proposed methodFigure 11
Near-field source location estimation using the sparse array (r
1
=
300 m, θ
1
=
3
0
° and r
2
=
1000 m, θ
2
=
6
0
°).
(a)
WVD
(b)
AOK distribution
(c)
OMP
(d)
The proposed method


## 5.1. Example I: Single Source with a Moderate Range

The first example considers the localization of a single near-field source with a range of 1000 m and DOA of 30° with respect to the reference sensor. The corresponding source coordinate is[
866.03
,
500
] m. In this case, the spatial frequency is approximately linear to the sensor position. The input signal-to-noise ratio (SNR) is 0 dB and 100 time samples are observed. In order to obtain a smooth figure, we oversample the spatial frequency on a multiple of 128. Figure 4 shows the SAD results for a ULA from both Wigner-Ville distribution (WVD) and sparse reconstruction, and the corresponding results for the sparse array are illustrated in Figure 5.Figure 4
Near-field source location estimation using the ULA (r
1
=
1000 m, θ
1
=
3
0
∘).
(a)
WVD
(b)
Sparse reconstructionFigure 5
Near-field source location estimation using the sparse array (r
1
=
1000 m, θ
1
=
3
0
°).
(a)
WVD
(b)
Sparse reconstruction
(c)
AOK distribution
(d)
Sparse reconstructionFigure4(a) shows the WVD of the SAD. Because there is no missing sensor, the signature can be estimated using a number of conventional methods as well as sparse reconstruction techniques, with respect to the SAF. The corresponding results obtained from the OMP and the structure-aware BCS are depicted in Figure 4(b). It is apparent that the signature is well approximated as a chirp and there is little difference noticed between the quadratic and third-order approximations. As such, in this case, the use of third-order approximation with the proposed BCS technique achieves similar performance. The PPS coefficients estimated from the proposed method are, respectively, b
=
[
0.4321
,
-
0.0013
,
-
3
×
1
0
-
6
]
T, and the yielding source position estimate is [
844.81
,
492.07
] m, which corresponds to a 22.6 m error to the true source position.The corresponding 16-element sparse array scenario is presented in Figure5, in which the conventional method developed in [8] is not effective. In addition, the SAD is contaminated with artifacts that may obscure the true feature estimation due to missing sensors, as shown in Figure 5(a). There are erroneous spatial frequency estimates for sparse reconstruction depicted in Figure 5(b). For the SAD technique, the modified AOK can be applied to suppress the effect of missing sensors. The results are shown in Figures 5(c) and 5(d). It is evident that a desired performance is achieved.Figure6 compares the RMSE performance as a function of the input SNR obtained from 100 independent trials. It is evident that the source localization performance is improved as the input SNR increases. The WVD in the ULA case presents the best performance since the WVD is known to be optimal in the analysis of linear FM signals [36, 37]. However, the performance of WVD degrades when the spatial frequency is nonlinear to the sensors positions or when some sensors are missing. The sparse array with the proposed method yields a similar performance but suffers from a gap of performance degradation. Such performance degradation partly stems from the signal power loss with the fewer number of sensors.Figure 6
RMSE versus SNR (r
1
=
1000 m, θ
1
=
3
0
∘).


## 5.2. Example II: Single-Source with a Close Range

In this example, the source is located at a range of 350 m and DOA of3
0
° with respect to the reference sensor. The corresponding coordinate is [
303.1
,
175
] m. As a consequence, the signature behaves as a PPS with respect to the sensor positions. The WVD and the estimated results obtained from the ULA and sparse array are shown in Figures 7 and 8, respectively. As indicated in Figure 7, the estimated spatial frequency is approximated as a third-order PPS with respect to the sensor positions. In this case, the WVD exhibits some cross-terms. The third-order PPS model yields more accurate results as compared to the quadratic approximation in this case to reflect the closer source range to the array. Under the third-order PPS approximation, the coefficient b
=
[
0.4277
,
-
0.004
,
-
0.0001
]
T is obtained, resulting in the source position [
286.14
,
173.16
] m, a 17.06 m shift from the true source position. For the sparse array, the effect of artifacts is again confirmed in Figures 8(a) and 8(b). This effect is significantly suppressed after the AOK is applied, as shown in Figures 8(c) and 8(d).Figure 7
Near-field source location estimation using the ULA (r
1
=
350 m, θ
1
=
3
0
∘).
(a)
WVD
(b)
Sparse reconstructionFigure 8
Near-field source location estimation using the sparse array (r
1
=
350 m, θ
1
=
3
0
∘).
(a)
WVD
(b)
Sparse reconstruction
(c)
AOK processing
(d)
Sparse reconstructionIn Figure9, we compare the RMSE performance of the second example. It is shown that the AOK-based results achieve better performance because of the effective cross-term and artifact suppression. In addition, there appears a floor in the obtained RMSE performance due to the spatial frequency resolution and bias. For the sparse array, a similar performance is attained to that of the ULA when the input SNR is higher than 0 dB and outperforms the WVD-based ULA results.Figure 9
RMSE versus SNR (r
1
=
350 m, θ
1
=
3
0
°).


## 5.3. Example III: Two-Source Scenario

Finally, we consider a two-source case in Figure10. The first source is located at r
=
300 m and θ
=
3
0
° with respect to the reference sensor, whereas the second source is located at r
=
1000 m and θ
=
6
0
∘. For the 32-element ULA, as shown in Figure 10(a), the SAD in this case is contaminated by the cross-terms between the two sources, making the spatial frequency estimation difficult. Such cross-term effect can be substantially reduced by exploiting the AOK, as shown in Figure 10(b). The sparse reconstruction results of OMP obtained from the AOK are depicted in Figure 10(c), where the spatial frequencies are not correctly identified towards both edges. Improved results are obtained from the proposed technique exploiting the structure-aware BCS, as depicted in Figure 10(d), due to the enforced spatial frequency continuity in the sparse SAD reconstruction. The latter evidently shows better performance for the multiple near-field source case. When the 16-element spare array is used, as shown in Figure 11(a), the SAD suffers from both cross-terms and artifacts due to missing sensor positions, yielding highly obscured SAD signatures. Both cross-terms and artifacts are significantly reduced when the AOK is used, as depicted in Figure 11(b). When it is combined with the proposed structure-aware BCS for sparse spatial frequency reconstruction, improved spatial frequency estimation, as shown in Figure 11(d), is achieved. For comparison, the OMP fails to identify one of the sources in the reconstructed spatial frequencies, as depicted in Figure 11(c).Figure 10
Near-field source location estimation using the ULA (r
1
=
300 m, θ
1
=
3
0
° and r
2
=
1000 m, θ
2
=
6
0
°).
(a)
WVD
(b)
AOK distribution
(c)
OMP
(d)
The proposed methodFigure 11
Near-field source location estimation using the sparse array (r
1
=
300 m, θ
1
=
3
0
° and r
2
=
1000 m, θ
2
=
6
0
°).
(a)
WVD
(b)
AOK distribution
(c)
OMP
(d)
The proposed method


## 6. Conclusion

In this paper, we proposed a simple near-field source localization technique using uniform and sparse linear arrays. We used the sensor-angle distribution to characterize the sensor-dependent phase progression as a function of the source range and its direction. Sparse reconstruction techniques were used to estimate the sensor-dependent spatial frequency signature for source localization. Since the degree of the polynomial phase signal viewed across sensors is unknown and depends on emitter proximity to the array, parameterized dictionary cannot be used. The effectiveness of the proposed technique was verified using simulation results.



---

*Source: 783467-2015-04-23.xml*