Sparsity-Aware Noise Subspace Fitting for DOA Estimation

Abstract

We propose a sparsity-aware noise subspace fitting (SANSF) algorithm for direction-of-arrival (DOA) estimation using an array of sensors. The proposed SANSF algorithm is developed from the optimally weighted noise subspace fitting criterion. Our formulation leads to a convex linearly constrained quadratic programming (LCQP) problem that enjoys global convergence without the need of accurate initialization and can be easily solved by existing LCQP solvers. Combining the weighted quadratic objective function, the ${\ell }_1$ norm, and the non-negative constraints, the proposed SANSF algorithm can enhance the sparsity of the solution. Numerical results based on simulations, using real measured ultrasonic, and radar data, show that, compared to existing sparsity-aware methods, the proposed SANSF can provide enhanced resolution with a lower computational burden.

1. Introduction

Direction-of-arrival (DOA) estimation is one of the key problems in array processing [1,2]. Over the last 50 years, a variety of techniques have been developed to address this problem [2,3,4,5,6,7,8,9,10,11,12,13,14], in which high resolution, low complexity, and good statistical performance are desired. Especially, the maximum likelihood (ML) method was investigated because it can achieve the Cramér-Rao low bound (CRLB) for large samples [5,15,16,17,18,19]. Although the ML method has asymptotic statistical efficiency, its computational cost is extremely high and its global solution is not guaranteed [20,21].

The noise subspace fitting (NSF) methods and their variants are popular approximations of the ML method. One of the well-known NSF methods is multiple signal classification (MUSIC), which is computationally simple and asymptotically equivalent to the ML method for uncorrelated signals [6]. In the case of correlated sources, however, MUSIC is unsatisfactory from a statistical viewpoint [5]. To deal with correlated sources, an improved NSF strategy called the optimally weighted noise subspace fitting (OWNSF) was developed to achieve the asymptotic property of ML, and the corresponding algorithm is named method of direction estimation (MODE) [5,16,21]. For uniform linear arrays (ULA), MODE can be implemented by eigen-analysis with the computational cost of $O\left(\frac{1}{2}(M+1)MT\right)$ flops, where M and T denote the number of sensors and snapshots, respectively [22]. For general array geometries, however, the implementation of MODE requires iterative non-linear search procedures (e.g., the Gauss–Newton search technique), in which an accurate initialization is crucial to avoid local minima [5,16]. An accurate initialization may be difficult when there exist closely spaced sources [16]. Therefore, MODE suffers from initialization difficulties and its global solution may not be guaranteed for general array geometries.

Based on the sparsity of the spatial spectrum, some sparsity-aware methods have been proposed recently to avoid the requirement on accurate initializations [8,9,10,23,24,25,26]. Among them, a representative method is the global matched filter (GMF), which is robust according to correlated and even coherent sources [27,28]. GMF is a quadratic minimization with ${\ell }_1$ norm regularization, so it has an unique global minimum [27]. In addition, by weighting the ${\ell }_1$ norm penalty item, a weighted GMF (WGMF) was proposed in [11]. However, WGMF requires solving a second order cone programming (SOCP), for which the computational cost of each iteration is on the order of $O\left(K^3{P}^3\right)$ complex operations with the interior-point algorithm [29], where P and K are the length of the sparse signal and the number of sources, respectively. A hyperparameter-free estimator, called the likelihood-based estimation of sparse parameters (LIKES), was also developed from the ML criterion [30]. The computational cost per iteration of LIKES is on the order of $O\left(2PMT\right)$ complex operations (typically $M\ll P$ for sparse representation) [31]. From a computational point of view, LIKES is attractive when $MT$ is moderate. However, LIKES also suffers from a large computational cost problem when $MT$ is large. A gridless sparse method called gridless sparse iterative covariance-based estimation (GL-SPICE, GLS) can estimate DOA without grid selection, which is very practical to relieve the grid mismatch problem of spare recovery [32,33]. The literature [34] addressed the DOA estimation with a block sparse reconstruction method in the presence of unknown mutual coupling. But GLS [33] and the method presented in [34] can only deal with linear arrays.

In this paper, we propose a sparsity-aware noise subspace fitting (SANSF) algorithm by extending the OWNSF criterion [16] to achieve high-resolution DOA estimation for general array geometries. The extension of the OWNSF criterion under the sparse representation framework leads to a new quadratic objective function that contains a positive semidefinite Hessian matrix. The total power and the non-negative property of the spectral components are used as the constraints of the objective function. In this way, the proposed SANSF algorithm is converted into a convex linearly constrained quadratic programming (LCQP) problem without the requirement of an accurate initialization. Furthermore, the proposed SANSF algorithm can be efficiently solved by many existing LCQP solvers [35,36,37]. Compared to SOCP, LCQP is computationally more efficient [38]. Simulations and experimental results using real measured ultrasonic data and radar data demonstrate that, compared to the existing sparsity-aware methods, SANSF can achieve higher resolution with relatively lower computational cost.

The proposed SANSF algorithm has three attractive features: (1) Its convex quadratic objective function enjoys global convergence without an accurate initialization requirement; (2) the SANSF problem is an LCQP so that it can be easily solved by efficient solvers with low computational costs; (3) not only the weighted quadratic objective function but also the ${\ell }_1$ norm and the non-negative constraints play the role of encouraging sparseness, which would be very helpful to improve the resolution.

The remainder of this paper is organized as follows. In Section 2, the signal model is introduced. In Section 3, the proposed SANSF algorithm is formulated, and its advantages are analyzed. In Section 4, experimental results on simulated and real measured data are presented to illustrate the performance of the proposed algorithm. Concluding remarks are given in Section 5.

2. Signal Model

2.1. Array Processing Model

Assume that there are M sensors with general array geometries and K far-field uncorrelated narrowband signals $\{s_k\left(t\right),k=1,2,\cdots ,K\}$, where the signals are impinging on the array from distinct directions $\{{\theta }_k,k=1,2,\cdots ,K\}$. The model of the received signal can be expressed as:

$$\mathbf{y}\left(t\right)=\mathbf{As}\left(t\right)+\mathbf{n}\left(t\right),t=1,2,\cdots ,T,$$

(1)

where $\mathbf{y}\left(t\right)$, $\mathbf{s}\left(t\right)$, $\mathbf{n}\left(t\right)$, and T denote the measurement vector, the signal vector, the noise vector, and the number of snapshots, respectively. The matrix $\mathbf{A}\in {\mathbb{C}}^{M\times K}$ is the array response matrix given by $\mathbf{A}=[\mathbf{a}\left({\theta }_1\right),\cdots ,\mathbf{a}\left({\theta }_K\right)]$, where $\mathbf{a}\left({\theta }_k\right)\in {\mathbb{C}}^{M\times 1}$ is the so-called steering vector depending on the array geometry. $\mathbf{n}\left(t\right)$ is an additive noise vector with zero-mean and covariance matrix ${\sigma }^2\mathbf{I}$. The noise $\mathbf{n}\left(t\right)$ is assumed to be uncorrelated to $\mathbf{s}\left(t\right)$.

The covariance matrix of $\mathbf{y}\left(t\right)$ can be expressed as

$$\begin{array}{ccc}\hfill \mathbf{R}& =& \mathrm{E}\left\{\mathbf{y}\left(t\right){\mathbf{y}}^{\mathrm{H}}\left(t\right)\right\}\hfill \\ & =& \mathbf{A}\mathrm{E}\left\{\mathbf{s}\left(t\right){\mathbf{s}}^{\mathrm{H}}\left(t\right)\right\}{\mathbf{A}}^{\mathrm{H}}+{\sigma }^2{\mathbf{I}}_M\hfill \\ & =& {\mathbf{AR}}_S{\mathbf{A}}^{\mathrm{H}}+{\sigma }^2{\mathbf{I}}_M,\hfill \end{array}$$

(2)

where ${\mathbf{R}}_S=\mathrm{E}\left\{\mathbf{s}\left(t\right){\mathbf{s}}^{\mathrm{H}}\left(t\right)\right\}$, ${\mathbf{I}}_M$ is a $M\times M$ identity matrix, and ${(·)}^{\mathrm{H}}$ denotes the conjugate transpose.

Under the assumption of uncorrelated signals, the eigen-decomposition of $\mathbf{R}$ can be written as

$$\begin{array}{c}\hfill \mathbf{R}=\sum _{m=1}^M{\lambda }_m{\mathbf{u}}_m{\mathbf{u}}_m^{\mathrm{H}}={\mathbf{U}}_S{\Lambda }_S{\mathbf{U}}_S^{\mathrm{H}}+{\sigma }^2{\mathbf{U}}_N{\mathbf{U}}_N^{\mathrm{H}},\end{array}$$

(3)

where ${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_K\ge {\lambda }_{K+1}=\cdots ={\lambda }_M={\sigma }^2$ are the eigenvalues of $\mathbf{R}$, ${\Lambda }_S=\mathrm{diag}\{{\lambda }_1,\cdots ,{\lambda }_K\}$, ${\Lambda }_N=\mathrm{diag}\{{\lambda }_{K+1},\cdots ,{\lambda }_M\}$, ${\mathbf{U}}_S$ and ${\mathbf{U}}_N$ span the signal and noise subspaces, respectively, $\mathrm{diag}\{·\}$ denotes the diagonal matrix. Consistent estimates of the signal and noise subspaces can be obtained by the eigendecomposition of the sample covariance matrix

$$\begin{array}{ccc}\hfill \widehat{\mathbf{R}}& =& \frac{1}{T}\sum _{t=1}^{\mathrm{T}}\mathbf{y}\left(t\right){\mathbf{y}}^{\mathrm{H}}\left(t\right)\hfill \\ & =& \sum _{m=1}^M{\widehat{\lambda }}_m{\widehat{\mathbf{u}}}_m{\widehat{\mathbf{u}}}_m^{\mathrm{H}}\hfill \\ & =& {\widehat{\mathbf{U}}}_S{\widehat{\Lambda }}_S{\widehat{\mathbf{U}}}_S^{\mathrm{H}}+{\widehat{\mathbf{U}}}_N{\widehat{\Lambda }}_N{\widehat{\mathbf{U}}}_N^{\mathrm{H}},\hfill \end{array}$$

(4)

where ${\widehat{\mathbf{U}}}_S=[{\widehat{\mathbf{u}}}_1,\cdots ,{\widehat{\mathbf{u}}}_K]$, ${\widehat{\mathbf{U}}}_N=[{\widehat{\mathbf{u}}}_{K+1},\cdots ,{\widehat{\mathbf{u}}}_M]$, ${\widehat{\Lambda }}_S=\mathrm{diag}\{{\widehat{\lambda }}_1,\cdots ,{\widehat{\lambda }}_K\}$, ${\widehat{\Lambda }}_N=\mathrm{diag}\{{\widehat{\lambda }}_{K+1},\cdots ,{\widehat{\lambda }}_M\}$, and ${\widehat{\lambda }}_m$ and ${\widehat{\mathbf{u}}}_m$ are the estimates of ${\lambda }_m$ and ${\mathbf{u}}_m$, respectively.

2.2. OWNSF and Its Sparse Representation

Next, the OWNSF criterion and its sparse representation are reviewed. The OWNSF criterion can be expressed by minimizing the following function [5,16]:

$$\begin{array}{c}\hfill f=\parallel {\widehat{\mathbf{U}}}_N^{\mathrm{H}}\mathbf{A}{\mathbf{W}}_{\mathrm{opt}}^{1/2}{\parallel }_F^2\end{array}$$

(5)

where ${\mathbf{W}}_{\mathrm{opt}}={\mathbf{A}}^{†}{\mathbf{U}}_S{({\Lambda }_S-{\sigma }^2{\mathbf{I}}_M)}^2{\Lambda }_S^{-1}{\mathbf{U}}_S^{\mathrm{H}}{\mathbf{A}}^{†\mathrm{H}}$, ${\parallel ·\parallel }_F$, and ${(·)}^{†}$ denote the Frobenius norm and the pseudo-inverse, respectively. In practice, for general array geometries, an iterative non-linear search procedure [5,16] can be used to implement OWNSF after replacing ${\mathbf{W}}_{\mathrm{opt}}$ with its estimate ${\widehat{\mathbf{W}}}_{\mathrm{opt}}={\mathbf{A}}^{†}\left({\mathsf{\theta }}_0\right){\widehat{\mathbf{U}}}_S{({\widehat{\Lambda }}_S-{\widehat{\sigma }}^2{\mathbf{I}}_M)}^2{\widehat{\Lambda }}_S^{-1}{\widehat{\mathbf{U}}}_S^{\mathrm{H}}{\mathbf{A}}^{†\mathrm{H}}\left({\mathsf{\theta }}_0\right)$, where ${\mathsf{\theta }}_0$ is the initial estimate of the DOAs.

The following relationship was verified in [5,16]:

$${\mathbf{A}}^{†}{\mathbf{U}}_S{({\Lambda }_S-{\sigma }^2{\mathbf{I}}_M)}^2{\Lambda }_S^{-1}{\mathbf{U}}_S^{\mathrm{H}}{\mathbf{A}}^{†\mathrm{H}}={\mathbf{R}}_S{\mathbf{A}}^{\mathrm{H}}{\mathbf{R}}^{-1}\mathbf{A}{\mathbf{R}}_S.$$

(6)

Substituting Equation (6) into Equation (5), we have

$$\begin{array}{ccc}\hfill f& =& \parallel {\widehat{\mathbf{U}}}_N^{\mathrm{H}}\mathbf{A}{\mathbf{W}}_{\mathrm{opt}}^{1/2}{\parallel }_F^2\hfill \\ & =& \mathrm{Tr}\left\{{\mathbf{W}}_{\mathrm{opt}}{\mathbf{A}}^{\mathrm{H}}{\widehat{\mathbf{U}}}_N{\widehat{\mathbf{U}}}_N^{\mathrm{H}}\mathbf{A}\right\}\hfill \\ & =& \mathrm{Tr}\left\{{\mathbf{R}}_S{\mathbf{A}}^{\mathrm{H}}{\mathbf{R}}^{-1}\mathbf{A}{\mathbf{R}}_s{\mathbf{A}}^{\mathrm{H}}{\widehat{\mathbf{U}}}_N{\widehat{\mathbf{U}}}_N^{\mathrm{H}}\mathbf{A}\right\}\hfill \\ & =& \mathrm{Tr}\left\{{\mathbf{R}}^{-1}\mathbf{A}{\mathbf{R}}_S{\mathbf{A}}^{\mathrm{H}}{\widehat{\mathbf{U}}}_N{\widehat{\mathbf{U}}}_N^{\mathrm{H}}\mathbf{A}{\mathbf{R}}_S{\mathbf{A}}^{\mathrm{H}}\right\},\hfill \end{array}$$

(7)

where $\mathrm{Tr}\{·\}$ denotes the trace of a matrix.

Under the assumption of the uncorrelated signals (The proposed SANSF algorithm is robustness to this assumption. As can be seen from the experimental section, SNASF works well even if there are some correlated sources.), a sparse representation of $\mathbf{A}{\mathbf{R}}_S{\mathbf{A}}^{\mathrm{H}}$ can be expressed as

$$\begin{array}{ccc}\hfill \mathbf{A}{\mathbf{R}}_S{\mathbf{A}}^{\mathrm{H}}& =& \mathbf{B}\mathsf{\Gamma }{\mathbf{B}}^{\mathrm{H}},\hfill \end{array}$$

(8)

where

$$\begin{array}{ccc}\hfill \mathsf{\Gamma }& \triangleq & \left[\begin{array}{cccc}{\gamma }_1& 0& \cdots & \cdots \\ 0& {\gamma }_2& 0& \cdots \\ ⋮& 0& \ddots & ⋮\\ 0& \cdots & \cdots & {\gamma }_P\end{array}\right]\in {\mathbb{R}}^{P\times P},\hfill \end{array}$$

(9)

$${\gamma }_p\triangleq \mathrm{E}\left\{s_p\left(t\right)s_p^{*}\left(t\right)\right\},$$

(10)

and

$$\mathbf{B}\triangleq [{\mathbf{b}}_1,\cdots ,{\mathbf{b}}_p,\cdots ,{\mathbf{b}}_P]\in {\mathbb{C}}^{M\times P},$$

(11)

$P\gg M$, ${(·)}^{*}$ denotes the complex conjugate, ${\mathbf{b}}_p\in {\mathbb{C}}^{M\times 1}$ is the steering vector for the p-th angular grid point, $p=1,\cdots ,P$, and P is the number of angular grid points.

The definition of ${\gamma }_p$ in Equation (10) implies that ${\gamma }_p\ge 0$. ${\gamma }_p>0$ denotes that there exists a source at the p-th angular grid point with the power value ${\gamma }_p$, and ${\gamma }_p=0$ if and only there is no source at the p-th angular grid point. Hence, the spatial spectrum can be expressed by $\mathsf{\gamma }$, where $\mathsf{\gamma }={[{\gamma }_1,\cdots ,{\gamma }_P]}^{\mathrm{T}}$, with ${(·)}^{\mathrm{T}}$ denoting the transpose. The support of $\mathsf{\gamma }$ can be defined as $\Omega \triangleq \left\{p\right|{\gamma }_p>0\}\subseteq \{1,\cdots ,P\}$, which results in $\mathbf{A}={\mathbf{B}}_{\Omega }$, where ${\mathbf{B}}_{\Omega }$ denotes a matrix whose columns are indexed by the set $\Omega$.

Substituting Equation (8) into Equation (7), we have

$$\begin{array}{ccc}\hfill f& =& \mathrm{Tr}\left\{{\mathbf{R}}^{-1}\mathbf{B}\mathsf{\Gamma }{\mathbf{B}}^{\mathrm{H}}{\widehat{\mathbf{U}}}_N{\widehat{\mathbf{U}}}_N^{\mathrm{H}}\mathbf{B}\mathsf{\Gamma }{\mathbf{B}}^{\mathrm{H}}\right\}\hfill \\ & =& \mathrm{Tr}\left\{\mathsf{\Gamma }{\mathbf{B}}^{\mathrm{H}}{\mathbf{R}}^{-1}\mathbf{B}\mathsf{\Gamma }{\mathbf{B}}^{\mathrm{H}}{\widehat{\mathbf{U}}}_N{\widehat{\mathbf{U}}}_N^{\mathrm{H}}\mathbf{B}\right\}.\hfill \end{array}$$

(12)

According to the above derivation, it is clear that minimizing Equation (12) can be regarded as an extension of the OWNSF criterion under the sparse framework. The main difference between Equation (7) and Equation (12) is that $\mathbf{A}{\mathbf{R}}_S{\mathbf{A}}^{\mathrm{H}}$ (i.e., ${\mathbf{B}}_{\Omega }{\mathbf{R}}_S{\mathbf{B}}_{\Omega }^{\mathrm{H}}$) in Equation (7) is replaced with its sparse representation $\mathbf{B}\mathsf{\Gamma }{\mathbf{B}}^{\mathrm{H}}$ in Equation (12). Minimizing the OWNSF criterion in Equation (5) requires an initial estimate of $\Omega$ (i.e., ${\mathsf{\theta }}_0$) and the corresponding algorithm is called MODE [5,16].

3. The Proposed SANSF Algorithm for DOA Estimation

In this section, we formulate the proposed SANSF algorithm and discuss its performance.

3.1. Sparsity-Aware Noise Subspace Fitting

By replacing $\mathbf{R}$ with $\widehat{\mathbf{R}}$ in Equation (12) and then minimizing Equation (12) with respecting to the unknown spatial spectrum $\mathsf{\gamma }$, we propose the following quadratic objective function

$$\begin{array}{ccc}\hfill \widehat{\mathsf{\gamma }}& =& \underset{\mathsf{\gamma }}{\mathrm{argmin}}\mathrm{Tr}\left\{\mathsf{\Gamma }{\mathbf{B}}^{\mathrm{H}}{\widehat{\mathbf{R}}}^{-1}\mathbf{B}\mathsf{\Gamma }{\mathbf{B}}^{\mathrm{H}}{\widehat{\mathbf{U}}}_N{\widehat{\mathbf{U}}}_N^{\mathrm{H}}\mathbf{B}\right\}\hfill \\ & =& \underset{\mathsf{\gamma }}{\mathrm{argmin}}\mathrm{Tr}\left\{\mathsf{\Gamma }\mathbf{C}\mathsf{\Gamma }\mathbf{D}\right\}\hfill \\ & =& \underset{\mathsf{\gamma }}{\mathrm{argmin}}{\mathsf{\gamma }}^{\mathrm{T}}\mathbf{G}\mathsf{\gamma },\hfill \end{array}$$

(13)

where the last equality is derived from the diagonal characteristics of $\mathsf{\Gamma }$, $\mathbf{C}={\mathbf{B}}^{\mathrm{H}}{\widehat{\mathbf{R}}}^{-1}\mathbf{B}$, $\mathbf{D}={\mathbf{B}}^{\mathrm{H}}{\widehat{\mathbf{U}}}_N{\widehat{\mathbf{U}}}_N^{\mathrm{H}}\mathbf{B}$, $\mathbf{G}=\mathrm{Real}(\mathbf{C}\odot {\mathbf{D}}^{\mathrm{T}})$, ⊙ and $\mathrm{Real}(·)$ denote the Hadamard product and taking the real part operators, respectively. It is clear that $\mathbf{G}$ can be computed from the known array geometry and the measurement vector, so Equation (13) is a quadratic function of the spatial spectrum $\mathsf{\gamma }$.

From Equations (2), (3), and (8), we have

$$\mathrm{Tr}\left\{\mathbf{B}\mathsf{\Gamma }{\mathbf{B}}^{\mathrm{H}}\right\}=\mathrm{Tr}\left\{{\mathbf{U}}_S({\Lambda }_S-{\sigma }^2\mathbf{I}){\mathbf{U}}_S^{\mathrm{H}}\right\}.$$

(14)

Due to the diagonal characteristics of $\mathsf{\Gamma }$, Equation (14) can be rewritten as

$$M\sum _{p=1}^P{\gamma }_p=\mathrm{Tr}\{{\Lambda }_S-{\sigma }^2\mathbf{I}\}.$$

(15)

Namely, a constraint on $\mathsf{\gamma }$ can be expressed by

$${\parallel \mathsf{\gamma }\parallel }_1={\mathbf{1}}^{\mathrm{T}}\mathsf{\gamma }=\mathrm{Tr}\{{\Lambda }_S-{\sigma }^2\mathbf{I}\}/M,$$

(16)

where ${\parallel ·\parallel }_1$ denotes the ${\ell }_1$ norm and $\mathbf{1}$ is an all-one vector. The constraint in Equation (16) suggests the magnitude of the retrieved spectrum.

By combining Equation (13), Equation (16), and the non-negative property of the spectrum components, the proposed SANSF algorithm can be expressed as the following quadratic minimization with the linear constraints:

$$\begin{array}{c}\hfill \underset{\mathsf{\gamma }\ge \mathbf{0}}{\min }{\mathsf{\gamma }}^{\mathrm{T}}\mathbf{G}\mathsf{\gamma }\mathrm{s}.\mathrm{t}.{\mathbf{1}}^{\mathrm{T}}\mathsf{\gamma }=\mathrm{Tr}\{{\Lambda }_S-{\sigma }^2\mathbf{I}\}/M,\end{array}$$

(17)

where $\mathsf{\gamma }\ge \mathbf{0}$ means that $\{{\gamma }_p\ge 0,p=1,\cdots ,P\}$. To solve Equation (17) in practice, ${\Lambda }_S$ and ${\sigma }^2$ need to be replaced with correspondingly consistent estimates ${\widehat{\Lambda }}_S$ and ${\widehat{\sigma }}^2$, respectively, where ${\widehat{\sigma }}^2={\sum }_{\widehat{K}+1}^M{\widehat{\lambda }}_m/(M-\widehat{K})$, and $\widehat{K}$ is an estimate of K that can be obtained by the minimum description length (MDL) criterion [39].

The optimization problem in Equation (17) is an LCQP problem that can be solved by off-the-shelf toolbox [35,36,40]. Here, we employ the MATLAB’s quadprog solver to implement SANSF, in which the computational complexity is on the order of $O\left(8P^3\right)$ real operations [41].

3.2. About the Initialization

The major drawback of MODE for general array geometries is the requirement for an accurate initialization to guarantee the convergence to the true solution [5,16,21,42]. Although the initialization can typically be performed by MUSIC or alternating projection (AP) [20], they may not work well when there are closely spaced sources [16]. The inaccurate initialization will seriously degrade the performance of MODE [5,16]. In contrast, SANSF proposed in Equation (17) has global convergence without any requirement for accurate initialization. We explain it in detail below.

The convexity of Equation (17) depends on the convexity of Equation (13), since the constraint in Equation (17) is linear. The sufficient condition of the convexity of Equation (13) is that the Hessian matrix $\mathbf{G}=\mathrm{Real}(\mathbf{C}\odot {\mathbf{D}}^{\mathrm{T}})$ is positive semidefinite. We first analyze the property of $\mathbf{C}$:

$$\begin{array}{ccc}\hfill \mathbf{C}& =& {\mathbf{B}}^{\mathrm{H}}{\widehat{\mathbf{R}}}^{-1}\mathbf{B}\hfill \\ & =& {\mathbf{B}}^{\mathrm{H}}\widehat{\mathbf{U}}{\widehat{\Lambda }}^{-1/2}{\widehat{\Lambda }}^{-1/2}{\widehat{\mathbf{U}}}^{\mathrm{H}}\mathbf{B}\hfill \\ & =& {\mathbf{E}}^{\mathrm{H}}\mathbf{E},\hfill \end{array}$$

(18)

where $\widehat{\mathbf{U}}=[{\widehat{\mathbf{U}}}_S,{\widehat{\mathbf{U}}}_N]$, $\widehat{\Lambda }=\mathrm{diag}\{{\widehat{\lambda }}_1,\cdots ,{\widehat{\lambda }}_M\}$, and $\mathbf{E}={\widehat{\Lambda }}^{-1/2}{\widehat{\mathbf{U}}}^{\mathrm{H}}\mathbf{B}$. Because both $\widehat{\Lambda }\in {\mathbb{R}}^{M\times M}$ and $\widehat{\mathbf{U}}\in {\mathbb{C}}^{M\times M}$ are full rank matrices, we have

$$\begin{array}{c}\hfill \mathrm{rank}\left(\mathbf{E}\right)=\mathrm{rank}\left(\mathbf{B}\right)=M,\end{array}$$

(19)

where $\mathrm{rank}(·)$ denotes the rank of the matrix.

According to Equation (19), we can conclude that ${\mathbf{EE}}^{\mathrm{H}}$ is nonsingular and its eigenvalues are positive numbers, which means that ${\mathbf{E}}^{\mathrm{H}}\mathbf{E}$ has M positive eigenvalues and $P-M$ zero-value eigenvalues [43]. As a result, $\mathbf{C}$ is positive semidefinite. Similarly, $\mathbf{D}$ is also positive semidefinite. Theorem 5.2.1 in [44] indicates that, if $\mathbf{C}$ and $\mathbf{D}$ are positive semidefinite, $\mathbf{G}=\mathrm{Real}(\mathbf{C}\odot {\mathbf{D}}^{\mathrm{T}})$ is also positive semidefinite. Based on these observations, we can conclude that the objective function of Equation (17) is convex. This means that SANSF ensures global convergence without any requirement on an accurate initialization.

3.3. The Effect of Sparsity Enhancement

We define the Capon-like weight as

$$c_{i,j}={\left[\mathbf{C}\right]}_{i,j}={\mathbf{b}}_i^{\mathrm{H}}{\widehat{\mathbf{R}}}^{-1}{\mathbf{b}}_j,$$

(20)

and the MUSIC-like weight as

$$d_{i,j}={\left[\mathbf{D}\right]}_{i,j}={\mathbf{b}}_i^{\mathrm{H}}{\widehat{\mathbf{U}}}_N{\widehat{\mathbf{U}}}_N^{\mathrm{H}}{\mathbf{b}}_j,$$

(21)

where i and $j\in \{1,\cdots ,P\}$, respectively.

The combination of these two weights is defined as

$$g_{i,j}={\left[\mathbf{G}\right]}_{i,j}=\mathrm{Real}\left(c_{i,j}{d}_{i,j}^{*}\right).$$

(22)

Based on these definitions, we can conclude that the quadratic minimization in Equation (17) is a weighted ${\ell }_2$ minimization, whose weights come from the Hadamard product of the Capon-like and the MUSIC-like weighting matrices.

As can been deduced from Equation (20) and (21), if and only if $i=j$, $c_{i,j}$ and $d_{i,j}$ become the regular Capon weight $c_i$ and the regular MUSIC weight $d_i$, respectively, where

$$c_i={\mathbf{b}}_i^{\mathrm{H}}{\widehat{\mathbf{R}}}^{-1}{\mathbf{b}}_i,$$

(23)

and

$$d_i={\mathbf{b}}_i^{\mathrm{H}}{\widehat{\mathbf{U}}}_N{\widehat{\mathbf{U}}}_N^{\mathrm{H}}{\mathbf{b}}_i.$$

(24)

The Capon and MUSIC spectral estimations can be directly obtained from $1/\mathrm{Real}\left(c_i\right)$ [3] and $1/\mathrm{Real}\left(d_i\right)$ [6], respectively. Incidentally, similar to Capon and MUSIC [45], the proposed SANSF cannot deal with coherent sources, which is a disadvantage in comparison to the original OWNSF.

Next, we will show that the proposed weighting scheme is consistent with the principle of the weighted ${\ell }_p$ minimization method [7,46,47,48,49]. Namely, the entries whose indices are more likely to be outside of the support are assigned larger weights so that they are suppressed to near zero, while the entries whose indices are more likely to be inside of the support are assigned smaller weights so that they are preferentially recovered. Without loss of generality, the index set $\{1,\cdots ,P\}$ is divided into $\Omega$ and ${\Omega }^c$, i.e., $\Omega \cup {\Omega }^c=\{1,\cdots ,P\}$ and $\Omega \cap {\Omega }^c=\varnothing$. Rewriting the objective function of SANSF in Equation (17), we have

$$\begin{array}{ccc}\hfill {\mathsf{\gamma }}^{\mathrm{T}}\mathbf{G}\mathsf{\gamma }& =& \sum _{i=1}^P\sum _{j=1}^P{g}_{i,j}{\gamma }_i{\gamma }_j\hfill \\ & =& \sum _{i\in \Omega }\sum _{j\in \Omega }{g}_{i,j}{\gamma }_i{\gamma }_j+\sum _{i\in {\Omega }^c}\sum _{j\in \Omega }{g}_{i,j}{\gamma }_i{\gamma }_j+\sum _{i\in \Omega }\sum _{j\in {\Omega }^c}{g}_{i,j}{\gamma }_i{\gamma }_j+\sum _{i\in {\Omega }^c}\sum _{j\in {\Omega }^c}{g}_{i,j}{\gamma }_i{\gamma }_j.\hfill \end{array}$$

(25)

A proof in Appendix A shows that $\left\{\right|g_{i,j}|,i\in \Omega \mathrm{or}j\in \Omega \}$ are smaller than $\left\{\right|g_{i,j}|,i\in {\Omega }^c\mathrm{and}j\in {\Omega }^c\}$. It can be observed from Equation (25) that the set of smaller weights $\{g_{i,j},i\in \Omega \mathrm{or}j\in \Omega \}$ and the signal-related items $\{{\gamma }_i{\gamma }_j,i\in \Omega \mathrm{or}j\in \Omega \}$ are linked. Assigning the set of smaller weights $\{g_{i,j},i\in \Omega \mathrm{or}j\in \Omega \}$ to each signal-related item raises the priority of the spectrum peak $\{{\gamma }_i,i\in \Omega \}$ or $\{{\gamma }_j,j\in \Omega \}$ in the recovery process. The positions that more likely correspond to signal-unrelated entries, i.e., $\{{\gamma }_i,i\in {\Omega }^c\}$ and $\{{\gamma }_j,j\in {\Omega }^c\}$, are forced to near zero by weighting with larger values $\{g_{i,j},i\in {\Omega }^c\mathrm{and}j\in {\Omega }^c\}$. Therefore, $\mathbf{G}$ links smaller/larger weights with the signal-related/signal-unrelated items.

The constraints of SANSF contain two part, i.e., ${\parallel \mathsf{\gamma }\parallel }_1=\mathrm{Tr}\{{\Lambda }_S-{\sigma }^2\mathbf{I}\}/M$ and $\mathsf{\gamma }\ge \mathbf{0}$, which can also promote the sparsity of the solution. The sparsity enforcing property of the ${\ell }_1$ norm constraint is natural because it is the closest convex relaxation of the exact sparse ${\ell }_0$ norm constraint [50]. Moreover, theoretical results in [51,52] show that the non-negative constraint is very helpful for sparsity-promoting. Therefore, combining the weighted quadratic objective function and the ${\ell }_1$ norm and the non-negative constraints, SANSF has a remarkable sparsity-promoting property.

4. Numerical and Experimental Results

In this section, some numerical results based on simulated data and real measured ultrasonic data are presented to demonstrate the performance of the proposed SANSF algorithm. The results of several DOA estimation methods, including MUSIC [6], WGMF [11], LIKES [30], and MODE [5], are also provided for comparison purposes. Here, the number of the inner iterations of LIKES is 150, and the number of iterations of MODE is 30 and its initial estimates comes from MUSIC. SANSF is implemented in Matlab function quadprog with default options except OptimalityTolerance is set to $1e^{-10}$. In all simulations, we consider a linear array with $M=10$ sensors located at $[0.5,1.0,3.0,3.5,6.0,6.5,7.0,7.5,8.0,10.0]$ in half-wavelength units, unless specified otherwise.

4.1. Root-Mean-Squared-Errors (RMSEs)

In the first experiment, the root-mean-squared-errors (RMSEs) of the DOA estimates of these methods are compared. Three uncorrelated sources with the same amplitude are impinging on the array from $\{-34.2°,1.5°,35.8°\}$. For comparison of the RMSEs of the five algorithms and CRLB, the grid refinement method is used to achieve better precision [8,10,53]. The initial searching space of $[-90°,90°]$ is uniformly sampled by 1126 grid points with the grid spacing of $0.16°$. Following the coarse search, three successive iterations are used to further refine the grid spacing, where the local interval around the k-th peak of the spectrum, i.e., $[{\widehat{\theta }}_k^{(\tau -1)}-\frac{8°}{4^{\tau -1}},{\widehat{\theta }}_k^{(\tau -1)}+\frac{8°}{4^{\tau -1}}]$, is uniformly sampled by 401 grid points, ${\widehat{\theta }}_k^{(\tau -1)}$ is the estimate of the k-th source at step $\tau -1$ (e.g., ${\widehat{\theta }}_k^{\left(0\right)}$ is the coarse estimate of the k-th source), and $\tau =1,2,3$. As shown in Figure 1a, SANSF can reach CRLB with lower SNR than WGMF and LIKES. In addition, as can be seen from Figure 1b, with the increasing number of snapshots, MUSIC, MODE, LIKES, and SANSF can reach the CRLB, but a large number of snapshots is required for LIKES. From Figure 1b it is also noted that the proposed SANSF algorithm is slightly worse than WGMF and LIKES when T is not sufficiently large. The reason is that the Capon-like and the MUSIC-like matrices fail to provide reliable weighting values for small values of T [54,55]. Especially, the proposed SANSF algorithm does not work for $T<M$. Furthermore, in contrast to MUSIC, LIKES, and SANSF that exploit the grid refinement method, WGMF is not able to reach the CRLB when the number of snapshots is large enough with the same grid refinement scheme. This fact shows that WGMF requires more sophisticated operations.

4.2. Resolution

4.2.1. Probability of Resolution

In the second experiment, the spatial resolution of SANSF is investigated. The azimuth interval $[-90°,90°]$ is uniformly sampled with 1801 grids. Two closely spaced sources are located at $\pm \Delta /2$, where $\Delta$ denotes the source separation. Two widely used criteria of resolution are employed. The definition of Criterion 1 [56,57] is that two sources are successfully resolved provided that the estimate ${\widehat{\theta }}_k$ is located in the neighborhood of the true position $U({\theta }_k;\eta )$. Here, the radius of the neighborhood $\eta$ is set as $\eta =0.49\Delta$. As for Criterion 2 [58], a solution would be declared as the successful resolution if and only if $\left[\widehat{\mathsf{\gamma }}({\theta }_1-r:{\theta }_1+r)+\widehat{\mathsf{\gamma }}({\theta }_2-r:{\theta }_2+r)\right]/(4r+2)-\widehat{\mathsf{\gamma }}\left(0\right)>0$, where r denotes the searching radius of estimated peaks (in this paper, $r=3$) (Here, the definition of Criterion 2 is slightly modified by considering the searching radius of the estimated peaks. If $r=0$, it is exactly the same as the original definition in [58].). The estimate of K comes from the MDL criterion [39]. As can be seen from Figure 2 (In Figure 2b,d, the results of MODE are absent because MODE does not yield the spatial spectrum that is necessary for Criterion 2. MODE will not be present in Figures 4 and 5 for the same reason.), SANSF has higher probability of successful resolution than the other four methods. The reason is that the sparsity of the solution of SANSF is enhanced by embedding the weighting matrix $\mathbf{G}$ and forcing the ${\ell }_1$ norm and the non-negative constraints. It is also noted that MODE cannot yield satisfactory resolution performance in these experimental results because the initialization value from MUSIC may be inaccurate for closely spaced sources.

4.2.2. Histogram for 1D Closely Spaced Sources

In the third experiment, we compare the 1D resolution in the signal-dense situation. The histogram of DOA estimates over 500 Monte-Carlo trials are plotted in Figure 3, where the azimuth interval $[-90°,90°]$ is uniformly sampled with 1801 grids and five uncorrelated and closely spaced sources with the same amplitude exist on the grid points $[801,851,901,951,1001]$. As demonstrated in Figure 3, histograms of MUSIC, MODE, WGMF, and LIKES show only four peaks around the true positions of the sources and a number of spurious peaks. This indicates that they fail to deal with the signal-dense situation. It is noteworthy that the proposed SANSF algorithm clearly shows five peaks at the true positions of the sources.

4.2.3. Spatial Spectrum for 2D Closely Spaced Sources

In the fourth experiment, two top views of the 2D spatial spectrum are used to further demonstrate the resolution performance of SANSF in the signal-dense situation. Five uncorrelated sources with amplitudes $(15,15,15,15,15)$ from $(-6°,8°)$, $(6°,8°)$,$(0°,16°)$, $(-8°,22°)$, and $(8°,22°)$ are considered, where the number of snapshots is 800. The azimuth interval $[-15°,15°]$ and the elevation interval $[0°,30°]$ are uniformly sampled with 61 grids, respectively. In Figure 4, we use an L-shaped array consisting of six horizontal array elements and five vertical array elements with half-wavelength spacing. As shown in Figure 4, the proposed SANSF algorithm clearly shows five peaks with no visible artifacts, while other algorithms are failing to resolve all of them. In Figure 5, a uniform circular array (UCA) consisting of 16 elements with half-wavelength spacing is employed to obtain the spatial spectrum. Again, the resolution performance of the proposed SANSF algorithm outperforms other algorithms. Especially, Figure 4a and Figure 5a, and Figure 4c and Figure 5c, show that MUSIC and LIKES suffer from low-resolution and high-sidelobes. Although WGMF has higher resolution than MUSIC and LIKES, as can be seen from Figure 4b and Figure 5b it is apt to underestimate the sources for L-shaped array and may suffer from undesired artifacts for UCA in the signal-dense situation.

4.3. Robustness to the Assumption of Uncorrelated Sources

In the fifth experiment, we consider the correlated sources scenario, where three correlated sources have the same amplitude with the correlation coefficient matrix $[1,0.8,0.9;0.8,1,0.8;0.9,0.8,1]$. In Figure 6, the averaged spatial spectrum over 100 Monte-Carlo trials is plotted. It is obvious that the proposed SANSF algorithm can deal with the correlated sources and is robustness to the assumption of uncorrelated sources. A detailed explanation for the robustness can be found in the literature [45,59].

4.4. Computational Complexity

In the sixth experiment, the average running time yielded by the above methods (i.e., MUSIC, MODE, WGMF, LIKES, and SANSF) are compared, where all the experiments are implemented with Matlab 2018b on a PC with the Intel CORE i5-7200U CPU @2.50GHz. Figure 7 shows that, compared to WGMF and LIKES, the proposed SANSF has a lower computational burden, especially for the smaller number of grid points. It is also noted that, with increasing number of snapshots, the computational complexities of WGMF and SANSF only slightly increase. Although MUSIC has the lowest computational cost among these algorithms, its resolution is not competitive as shown in Figure 2, Figure 3, Figure 4 and Figure 5.

4.5. DOA Tracking Results with Real Measured Ultrasonic Data

The real measured ultrasonic data are from the University of Wyoming Source Tracking Array Testbed (UW STAT) [60]. Here, the data set Number 2, that has been used for evaluating DOA performance in [61,62], is employed. One stationary source and one moving source were present in the scene, the number of sensors was 6, the number of snapshots was 1533, the carrier frequency was 40 KHz, and the signal bandwidth was 200 Hz. The sensor array was an ULA with the interelement spacing $2.1$ times of the wavelength. The sliding window technique with the size of 150 snapshots and the forward–backward smoothing are used to obtain the source trajectories.

From Figure 8, it is obvious that not all methods are capable of dealing with the intricate scene that the sources are very closely spaced. However, one can observe that SANSF has better threshold performance than the other three methods. The algorithm performance would “break down” when these two sources are too closely spaced [61]. The performance degradation interval due to the sources becoming closely spaced is called a breakdown interval. The breakdown interval shows the ability of the method to resolve the closely spaced sources. The smaller breakdown interval is, the higher the resolution.

Table 1. Breakdown interval.

Algorithms	MUSIC	MODE-ULA	WGMF	LIKES	SANSF
Start	488	314	560	426	601
End	789	809	680	847	715

provides the breakdown interval of these methods. Compared with the other methods, SANSF has the narrowest breakdown interval which is about 114 snapshots. This further demonstrates that SANSF can yield superior resolution performance in practice.

4.6. DOA Eestimation Results with Real Radar Data

The cross-range imaging of radar can be regarded as the DOA estimation problem in the azimuth dimension [63]. Therefore, we can evaluate the DOA estimation resolution by the cross-range resolution. The real radar data were collected by the TI IWR1642BOOST evaluation board that is a frequency modulated continuous wave (FMCW) sensor. The sensor configuration is shown in

Table 2. Radar configuration paramerers.

Parameters	Tx Channel	Rx Channel	Start Frequency	Frequency Slope	ADC Samples	Sample Rate
Value	Tx1	Rx1,2,3,4	77 GHz	66.626 MHz/$\mathsf{\mu }$s	256	5000 ksps

, in which four receivers and one transmitter form a four-element ULA. As seen from Figure 9, two screws with a diameter of about 1.5 cm and a height of 11 cm were used as observation targets. The distance from the radar to the midpoint between the two targets is about 50 cm. Center-to-center cross-range spacings between these two screw targets are about 12.5 cm, 15.0 cm, 17.5 cm, and 20 cm, respectively. As a reference, the cross-range Rayleigh resolution limit [63] for our scenario is approximately 12.5 cm.

We obtained the cross-range (i.e., azimuth DOAs) of two screws after performing two FFT that were used to estimate the range and velocity information of screws. As can be seen from Figure 10, the proposed SANSF algorithm can deal with the real radar data with higher azimuth resolution compared to other algorithms. However, MODE, WGMF, and LIKES failed to distinguish two closely spaced targets when their range spacings are comparable to the Rayleigh resolution limit. MUSIC distinguishes two closely spaced targets but it has higher sidelobes. Obviously, all methods can clearly distinguish two screws when their cross-range spacing are large enough.

5. Conclusions

In this paper, a sparsity-aware SANSF algorithm has been introduced for DOA estimation. SANSF is derived from the OWNSF criterion and involves a convex quadratic minimization problem subject to the linear constraints. This allows us to solve SANSF by using LCQP solvers at an acceptable computational cost. A Hadamard product of the MUSIC-like and the Capon-like matrices is embedded into the objective function, which embodies the weighted ${\ell }_2$ minimization methodology. The weighting and the ${\ell }_1$ norm and the non-negative constraints can together enhance the sparsity of the solution. Simulations and experimental results on measured ultrasonic and radar data have demonstrated that the proposed SANSF algorithm has superior resolution performance over typical sparsity-driven methods for noncoherent sources.