#### 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^{2}” 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^{2} 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], l*p*-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 l*p*-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^{2} 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^{2} 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^{2} 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^{2} 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^{2} 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^{2} 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^{2} 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.

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^{2} 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^{2} 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^{2} 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^{2} prior is a sparsity-encouraging prior which has sharper peak and heavier tails than the existing priors.

**(a)**

**(b)**

##### 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

Then, the Fisher information matrix (FIM) [32] can be obtained as

It is noted that the following equality holds with and . Finally, by combing (32) and (33), we have

#### 5. Simulation Results

In this section, several simulations are presented to illustrate the performance of the proposed algorithm as compared with the BCS [25] and NP-1 [27] algorithms. Consider a 9-sensor ULA with sensor interspacing being a half of the carrier wavelength. In the proposed algorithm, hyperparameters and are initialized as and 0.1; the values of rate parameter , shape parameter _{,} and power threshold are, respectively, set to be , , and . Two hundred independent Monte Carlo trials are conducted for the following simulations, and the success rate and root mean squared error (RMSE) are two significant performance metrics for evaluating the DOA estimation performance. The success rate is defined as the rate between the number of successful estimates and the total number of Monte Carlo trials, where a successful estimate is declared if the estimation error is within a certain small Euclidean distance of the true directions. The RMSE is defined by
where is the estimates of in the th Monte Carlo trial.

In the first simulation, we evaluate the effectiveness of the proposed algorithm. Assume that two uncorrelated sources from and two coherent sources from with fading coefficients impinge on this ULA. The grid number, SNR, and number of snapshots are, respectively, set to be 361, 0 dB, and 200. Figure 3 plots the source power spectrum, which is a function of DOA.

As can be seen from Figure 3, the power peaks in the vicinity of the true DOAs can be distinguished from other peaks by their relatively strong powers and the biases between these peaks and the true DOAs are slight, which demonstrate the effectiveness of the proposed algorithm.

The second simulation compares the estimation performance of three algorithms, including the proposed algorithm, BCS, and NP-1 algorithms. Consider one uncorrelated source from and two coherent sources from with fading coefficients being . The grid number and SNR are, respectively, fixed at 361 and 0 dB. The success rate versus number of snapshots and SNR are depicted in Figures 4 and 5, whereas the RMSE versus number of snapshots and SNR are shown in Figures 6 and 7. The results from Figures 4 and 5 illustrate that all the three algorithms are able to estimate the DOA correctly under the condition that the SNR or number of snapshots are relatively high, but for low SNR or small number of snapshots, the proposed algorithm has a higher success rate than the BCS and NP-1 algorithms. Figures 6 and 7 show that the estimation accuracy of the three algorithms tends to improve with the increase of SNR or the number of snapshots, and the proposed algorithm yields the best estimation accuracy as compared to the other two algorithms. This is because the proposed Gauss-Exp-Chi^{2} prior is a three-layer sparsity-encouraging prior, which can help to improve the accuracy of source reconstruction.

The third simulation investigates the RMSE versus grid number. The simulation settings are the same as those of the second simulation, except that grid number in this simulation is ranged from 61 to 361. Figure 8 plots the RMSE curves versus grid number with fixed SNR 10 dB and number of snapshots 500.

It is observed that the RMSE curves of the three algorithms rapidly decrease as grid number increases from 61 to 181, and these curves tend to slowly decrease with the increase of grid number when grid number is beyond 181. In addition, the proposed algorithm has minimum estimation errors among the three algorithms. Note that the three algorithms recover sparse sources from a Bayesian perspective; thus, they are not restricted to the restricted isometry property. This implies that the grid number can continue to increase to further improve the estimation accuracy and the computational complexity would increase accordingly. Therefore, for the purpose of balancing estimation accuracy and computational complexity, the reasonable range of grid number is from 181 to 361.

In the last simulation, we test the angular resolution of the three algorithms. Consider two coherent sources from 10.1° and with the angular separation ranging from 3° to 10°. The fading coefficients are , and the grid number is set to be 361. Figure 9 depicts the RMSE versus angular separation with the fixed SNR 10 dB and number of snapshots 500. The results from Figure 9 show that both the BCS and NP-1 algorithms fail to work when the angular separation is less than 4°, whereas the proposed algorithm can still provide accurate DOA estimation as long as the angular separation is no less than 3°. Moreover, the proposed algorithm has smaller RMSE than the BCS and NP-1 algorithms at the same angular separation. Thus, we can say that the proposed algorithm has a higher angular resolution than the BCS and NP-1 algorithms.

#### 6. Conclusion

In this paper, we develop a DOA estimation algorithm for a mixture of uncorrelated and coherent sources using SBI. In the Bayesian framework, a Gauss-Exp-Chi^{2} prior is introduced to promote sparse solutions, and the corresponding three-layer hierarchical Bayesian model is established. Then, the model parameters are estimated via variational Bayesian approximations. Finally, we form the source power spectra by using the estimated model parameters, from which the number and the locations of the highest peaks are extracted to achieve source number and DOA estimates. Simulation results show that the proposed algorithm can effectively estimate the DOAs of mixed sources and has better estimation performance than the state-of-the-art BCS algorithm and NP-1 algorithm in terms of estimation accuracy, success rate, and angular resolution.

#### Appendix

#### A. The Derivation of the Probability Density Function for Gauss-Exp-Chi^{2} Prior

In this Appendix, we prove that (12) holds.

The expression on the right hand side of (12) can be rewritten as

Thus, we first calculate the integral over

By substituting (A3) into (A2), we have where . Let (A4) can be expressed as where . According to the definition of the confluent hypergeometric function, (A5) can be further simplified to

Thus, we have

#### Conflicts of Interest

The authors declare that they have no competing interests.

#### Authors’ Contributions

Pinjiao Zhao proposed the main idea and conceived the proposed approach. Pinjiao Zhao and Weijian Si discussed and designed the proposed algorithm. Pinjiao Zhao performed the experiments and wrote the paper. Guobing Hu and Liwei Wang reviewed and revised the manuscript. All authors read and approved the manuscript.

#### Acknowledgments

This work was financially supported by the High Level Talent Research Starting Project of Jinling Institute of Technology (Grant no. jit-b-201724), the National Natural Science Foundation of China (Grant no. 61671168), and the Natural Science Foundation of the Jiangsu Province (Project no. BK20161104) and the Six Talent Peaks Project of the Jiangsu Province (Project no. DZXX-022).