Abstract

Direction of arrival (DOA) estimation algorithms based on sparse Bayesian inference (SBI) can effectively estimate coherent sources without recurring to extra decorrelation techniques, and their estimation performance is highly dependent on the selection of sparse prior. Specifically, the specified sparse prior is expected to concentrate its mass on the zero and distribute with heavy tails; otherwise, these algorithms may suffer from performance degradation. In this paper, we introduce a new sparse-encouraging prior, referred to as “Gauss-Exp-Chi²” prior, and develop an efficient DOA estimation algorithm for a mixture of uncorrelated and coherent sources under a hierarchical SBI framework. The Gauss-Exp-Chi² prior distribution exhibits a sharp peak at the origin and heavy tails, and this property makes it an appropriate prior to encourage sparse solutions. A three-layer hierarchical sparse Bayesian model is established. Then, by exploiting variational Bayesian approximation, the model parameters are estimated by alternately updating until Kullback-Leibler (KL) divergence between the true posterior and the variational approximation becomes zero. By constructing the source power spectra with the estimated model parameters, the number and locations of the highest peaks are extracted to obtain source number and DOA estimates. In addition, some implementation details for algorithm optimization are discussed and the Cramér-Rao bound (CRB) of DOA estimation is derived. Simulation results demonstrate the effectiveness and favorable performance of the proposed algorithm as compared with the state-of-the-art sparse Bayesian algorithms.

1. Introduction

Direction of arrival (DOA) estimation has been a crucial issue in various application areas involving radar, wireless communication, and navigation [1–3]. Multiple signal classification (MUSIC) [4] and estimation of signal parameter via rotational invariance technique (ESPRIT) [5], known as the two most classical subspace-based DOA estimation algorithms, have been proposed to resolve uncorrelated sources. However, in multipath propagation environments, sources from an identical target may undergo reflection from various surfaces, and thus, the received sources may be a mixture of uncorrelated and coherent sources. In such environments, the aforementioned subspace-based algorithms would suffer from serious performance deterioration or even fail [6].

In order to solve the aforementioned problem, several preprocessing techniques are developed for decorrelation. In the related studies, these preprocessing techniques are mainly classified into two categories: spatial smoothing (SS) [7] and matrix reconstruction (MR) [8]. Theoretically, the SS technique is implemented by dividing the whole array into multiple subarrays to combat the rank deficiency, while the MR technique is performed by rearranging the rank-deficient matrix into Hankel matrix, Toeplitz matrix, or other matrices to restore the rank. By combing the SS or MR techniques with subspace-based algorithms, several DOA estimation algorithms have been developed to handle the scenarios where uncorrelated and coherent sources coexist [9–13]. More specifically, the oblique projection spatial smoothing (OPSS) [9] and moduli property spatial smoothing (MPSS) [10] are the typical SS-based algorithms, and the oblique projection Toeplitz matrix reconstruction (OPT-MR) [11], symmetric non-Toeplitz matrix reconstruction (SNT-MR) [12], and real-valued Hankel matrix reconstruction (RVH-MR) [13] are the popular MR-based algorithms.

Unlike traditional subspace-based DOA estimation algorithms, the emerging sparse source reconstruction (SSR) algorithms [14–19], including matching pursuit (MP) algorithm [14], lp-norm optimization algorithms [15, 16], and sparse Bayesian inference (SBI) algorithms [17–19], provide a new perspective for DOA estimation. Since SSR-based algorithms realize DOA estimation via sparse source reconstruction, instead of calculating the covariance matrix, they can resolve the coherent sources directly without extra preprocessing techniques. Among these SSR-based algorithms, both the MP and the lp-norm optimization algorithms are restricted to solve an optimization problem and are relied on the point estimate, so that these algorithms ignore the uncertainty of the underlying source in the process of source reconstruction. By contrast, SBI-based algorithms specially consider the uncertainty of the underlying source and estimate the source via choosing an appropriate prior, which yield favorable reconstruction performance [20]. Moreover, they can provide good estimation performance in the case of low SNR or small number of snapshots [21]. In general, SBI-based algorithms first specify a sparsity-encouraging prior to the unknown source model and then the model parameters are estimated via Bayesian inference. Since the exact Bayesian inference is intractable, two mainstream approximation inference algorithms were presented to estimate the model parameters, in which one is evidence procedure [22] and the other is variational Bayesian approximation [23]. For the evidence procedure, some unknown hyperparameters with respect to the hierarchical prior model are estimated iteratively by maximizing the evidence. For the variational Bayesian approximation, the posterior distribution is approximated as the product of several tractable distributions, and the model parameters keep updating to minimize the Kullback-Leibler (KL) divergence between the true posterior and the variational approximation, which has attractive computational efficiency along with high estimation performance [24]. Note that these approximation inference algorithms operate under the premise that an appropriate sparse prior has been imposed on the source model for the purpose of encouraging sparse solutions. Many sparsity-encouraging prior models have been investigated in the SBI framework [25–27]. In [25], Bayesian compressed sensing (BCS) was proposed with a two-layer hierarchical Gaussian-inverse-gamma prior (or Student’s- prior), where the first layer is a Gaussian probability density function (pdf) and the second layer is a gamma pdf. Babacan et al. [26] proposed an equivalent Laplace prior that is generated by a Gaussian prior and an Exponential prior. In [27], a normal product (NP) prior was developed with two algorithms: NP-0 (using one-layer source model) and NP-1 (using two-layer source model). However, these priors are concentrated near the origin with relatively light tails, which may cause overshrinkage of the incident sources [28].

In this paper, we develop a new sparse-encouraging prior (called Gauss-Exp-Chi² prior) whose pdf distribution exhibits a sharp peak at the origin and heavy tails. With the proposed prior, the DOA estimation for a mixture of uncorrelated and coherent sources is performed under the hierarchical SBI framework using a uniform linear array (ULA). A three-layer hierarchical Bayesian model is established based on the Gauss-Exp-Chi² prior. Subsequently, according to the variational Bayesian approximation, the model parameters (including the mean and variance of sparse sources and hyperparameters) keep alternately updating until the KL divergence between the true posterior and the variational approximation tends to be zero. By exploiting the estimated model parameters, the source power spectra is constructed, from which the number and locations of the highest peaks are extracted to obtain source number and DOA estimates. Simulation results show that the proposed algorithm has superior estimation performance. Now we briefly summarize the contributions of this work as follows: (i)To encourage sparse solutions, we develop a new sparse-encouraging prior, called Gauss-Exp-Chi² prior, whose pdf distribution has a sharp peak at the origin and heavy tails.(ii)By constructing the source powers of all the potential directions in the angular space, both source number and DOA estimates are obtained.(iii)Variational approximations are adopted for the estimation of the hierarchical sparse Bayesian model parameters.(iv)Several implementation details for algorithm optimization including Woodbury matrix identity for dimension-reduction, pruning a basis function and the third kind Bessel function approximation are discussed, and the CRB of DOA estimation is derived.

The remainder of this paper is organized as follows. The DOA estimation model for mixed sources is formulated in Section 2. Section 3 presents the Gauss-Exp-Chi² prior and DOA estimation algorithm for a mixture of uncorrelated and coherent sources within the hierarchical SBI framework. The algorithm optimization and CRB of DOA estimation are discussed in Section 4. Section 5 presents the simulation results of the proposed algorithm. Conclusions are drawn in Section 6.

Notations. Vectors and matrices are denoted by lowercase and uppercase boldface letters, respectively. , , , and represent transpose, conjugate transpose, inverse, and the statistical expectation, respectively. is an identity matrix. denotes the Euclidean norm, and denotes a diagonal matrix. Additionally, represents the integral of from to .

2. Problem Formulation

Consider a total of far-field narrowband sources impinging on the uniform linear array (ULA) consisting of omnidirectional sensors with the interspacing between adjacent sensors being a half of the carrier wavelength , that is, . Then, the array output vector at the th snapshot can be expressed as where is the number of snapshots, is the steering vector corresponding to the direction of the th impinging source , and denotes the noise vector. Without loss of generality, the impinging source is a mixture of uncorrelated and coherent sources. Note that there exists power fading among the coherent source, and fading coefficients are used to describe the degree of power fading [5]. Specifically, consider far-field narrowband sources impinging on the ULA, in which the first sources are uncorrelated, and the last sources are coherent. Then, the uncorrelated source set is denoted as and the coherent source set is denoted as with are the fading coefficients.

Divide the entire angular space into sampling grids , where represents the grid number and generally satisfies . Assume that is the nearest grid to true direction ; thus, we have . Thus, can be rewritten in a sparse form where and . Due to the fact that has non-zero elements in elements, it is a sparse vector. In the case of multiple snapshots, the sparse sources at all the snapshots share the same support, and the array output matrix of snapshots can be represented by where , , and . The goal of this paper is to provide the DOA estimation under the coexistence of uncorrelated and coherent sources from a sparse Bayesian perspective.

3. DOA Estimation

In this section, a DOA estimation algorithm for a mixture of uncorrelated and coherent sources is proposed within the hierarchical SBI framework. A Gauss-Exp-Chi² prior is developed to encourage sparse solutions, and then the parameters of three-layer hierarchical Bayesian model are estimated via variational Bayesian approximations. By constructing source power spectra, the source number and DOA estimation are obtained.

3.1. Bayesian Model

In the Bayesian model, the pdf of a joint distribution with respect to all the unknown and observed quantities is required to be known. In this paper, the joint distribution can be expressed as where and are referred to as the hyperparameters; and are, respectively, referred to as the rate parameter and shape parameter. It is assumed that the components are the independently zero-mean stationary Gaussian noise with known variance . Thus, the pdf of the noise matrix is given by

Combining (3) and (5), the Gaussian likelihood model is obtained as follows:

From a Bayesian perspective, the pdf distribution of an assigned prior is appealing to exhibit a sharp peak at the origin and heavy tails, which favors strong shrinkage of noise sources and avoids overshrinkage of the interest sources. This property is generally considered as a desirable property for enforcing sparsity and variable selection [24]. Some typical sparse priors, such as Gaussian-inverse-gamma prior and Laplace prior, are widely used in the relevant research [25, 26], which, however, are concentrated near the origin with relatively light tails [28]. To alleviate this problem, we here develop a three-layer hierarchical prior, referred to as Gauss-Exp-Chi² prior, for . In the first layer of prior, we adopt a zero-mean Gaussian prior with being . In the second layer of prior, an exponential hyperprior is imposed on where denotes the exponential distribution. In the third layer of prior, a chi-square (Chi2) hyperprior is considered over , that is, where denotes the gamma function.

The first two layers (7) and (8) of the proposed prior result in a generative Laplace (Gaussian-Exponential) distribution [26], and the last layer (i.e., a Chi2 distribution) is embedded to obtain the proposed three-layer Gauss-Exp-Chi² prior. Thus, the proposed prior has more free parameters to control the degree of sparseness as compared with the Laplace prior [29].

Based on the above analysis, the directed graph of the sparse Bayesian model is shown in Figure 1, where arrows are used to point to the generative model. Note that the first five blocks from the left (corresponding to the variables , , , , and ) depict the three-layer hierarchical prior mentioned above.

Figure 1 

Directed graph of the sparse Bayesian model.

Since , it is obvious that most for are favored to zero, which leads to the fact that most are favored to zero. Thus, the proposed three-layer hierarchical Gauss-Exp-Chi² prior is a sparsity-encouraging prior, which can help to improve the performance of source reconstruction [30]. By combining the three layers of the proposed prior, the generated prior can be computed via where

Note that ; thus, can be computed by marginalizing over (see Appendix for detail). where is the confluent hypergeometric function defined by

The reason for choosing the Gauss-Exp-Chi² prior is twofold: its pdf distribution has a sharp spike at the origin and heavy tails; it forms in a conjugate manner since the exponential and Chi2 distributions belonging to the exponential distribution family are chosen as the 2nd layer and the 3rd layer of this prior, which significantly simplifies the form of posterior distributions [24]. The tail to spike behavior for Gauss-Exp-Chi² prior, Laplace prior, Student’s- prior, and Gaussian prior is illustrated in Figure 2, wherein the parameter selection of these priors is according to the standard derivation [25, 26]. It is observed that the proposed Gauss-Exp-Chi² prior is a sparsity-encouraging prior which has sharper peak and heavier tails than the existing priors.

(a)

(b)

(a)
(b)

Figure 2 

Four kinds of pdf curves with the standard derivation: (a) compares four distributions at origin and (b) compares four distributions in tail.

3.2. Variational Bayesian Approximations

It is well known that Bayesian inference operates on the basis of the posterior distribution . However, this posterior distribution is intractable for the reason that cannot be calculated analytically. In order to approximate this posterior distribution, we resort to mean-field variational Bayesian approximation [23], whose task is to seek several analytically approximate distributions that minimize the KL divergence between the posterior and its approximation distribution. According to [23], can be factorized as the product of three variational distributions , , and ; thus, we have

Denote a set , the distribution of each variable is expressed as , where represents the subset excluding . More specifically, the posterior approximation of , , and in this paper are, respectively, given by

(1) Update of . Keeping the terms of only depend on , we have with . By gathering together the similar terms, we have which implies that is the product of multiple Gaussian distributions with the mean and covariance being

(2) Update of . Similarly, according to (16), is derived as

Denoting and , (22) can be further simplified to

The posterior distribution of can be approximately equivalent to the product of a series of generalized inverse Gaussian (GIG) distributions, that is, , and thus, we have where is the third kind Bessel function with order . The case of in (24) gives the calculation of , which can be used in (25). The case of gives used for the evaluation of in (21).

(3) Update of . For , we have

Thus, it can be known that is a gamma distribution, that is, with the following mean

The estimation of model parameters is implemented by alternately updating the mean , variance , and hyperparameters and to minimize the KL divergence between the true posterior and the variational approximation, and the main steps are summarized in Algorithm 1.

3.3. Source Power and DOA Estimation

The DOAs of the impinging sources are estimated via the following two steps: (1) form spatial spectra using the estimated source powers of all the potential directions and (2) extract the number and the corresponding locations of the peaks beyond a power threshold from the spectra. Let and be the final estimates of and (the mean and variance with respect to ) and consider row by row, then we have . As outlined in [18], let be the source power of direction ; then we have

Therefore, the source powers of all the potential directions in the angular space are estimated by calculating . Suppose that there are peaks exceeding the power threshold and the corresponding grid indices are . Then, the estimated source number and DOAs will be and .

4. Discussions

4.1. Algorithm Optimization

For alleviating the computational complexity and speeding up the update efficiency, some implementation details are adopted for algorithm optimization in this subsection.

4.1.1. Woodbury Matrix Identity for Dimension-Reduction

At each iteration, it is inevitable to calculate a matrix inversion when updating according to (21), which requires computations. To reduce the computational complexity, we adopt the Woodbury matrix identity, and then (21) can be rewritten as with an matrix . Thus, the computational complexity is reduced to .

4.1.2. Pruning a Basis Function

In order to speed up the update efficiency, and the corresponding would be pruned from the original sparse Bayesian model when is smaller than a certain threshold. This procedure is referred to as “pruning a basis function”, and similar approaches have been used in [29, 30].

4.1.3. The Third Kind Bessel Function Approximation

It is known that (24) involves the third kind Bessel function, and this integral process is relatively complicated. According to [31], the third kind Bessel function of can be approximated as , if , and , if . Then, (24) can be simplified to

4.2. Cramér-Rao Bound (CRB)

To evaluate the DOA estimation performance, the CRB of the DOA estimation for a mixture of uncorrelated and coherent sources is derived in this subsection. Consider the sparse form of the array output vector mentioned in (2), the log-likelihood function of is firstly construct as

Let and be the real part and imaginary part of , that is, and . Then, the partial derivatives of (29) with respect to , , , and are, respectively, given as where and . Subsequently, define