Bayesian Direction of Arrival Estimation with Prior Knowledge from Target Tracker

Abstract

The performance of traditional direction of arrival (DOA) estimation methods always deteriorates at a low signal-to-noise ratio (SNR) or without sufficient observations. This paper investigates the Bayesian DOA estimation problem aided by the prior knowledge from the target tracker. The Bayesian Cramér–Rao lower bounds (CRLB) and the expected CRLB are first derived to evaluate the theoretical performance of Bayesian DOA estimation. Based on the maximum a posterior (MAP) estimator in the Bayesian framework, two methods are proposed. One is a two-step grid search method for a single target DOA case. The other is a gradient-based iterative solution for multiple targets DOA case, which extends the traditional Newton method by incorporating the prior knowledge. We also propose a minimum mean square error (MMSE) estimator using a Monte Carlo method, which requires trading off accuracy against computational complexity. By comparing with the maximum likelihood (ML) estimators and the MUSIC algorithm, the proposed three Bayesian estimators improve the DOA estimation performance in low SNR or with limited snapshots. Moreover, the performance is not affected by the correlation between sources.

1. Introduction

Direction of arrival (DOA) estimation using sensor arrays is a fundamental theme in array signal processing with many important applications in radar, sonar, wireless sensor networks and wireless communications [1,2]. There are numerous DOA estimation methods in the research by using the array observation data. These methods can be broadly classified to several categories, such as traditional spectral estimation methods [3], maximum likelihood (ML) methods [4,5], subspace-based methods [6,7,8] and sparse methods [9,10,11,12].

The conventional (or Bartlett) beamformer and Capon’s beamformer are two typical spectral estimation methods. The spatial spectrum in the former suffers from the same resolution limitation as the periodogram. The latter was later proposed to improve the resolution capability of closely spaced sources, but its performance deteriorated in low signal-to-noise ratio (SNR) [3]. The ML method is a well known and frequently used approach in signal processing. In DOA estimation, two different ML methods, i.e., deterministic ML and stochastic ML, were proposed, and they are based on two different statistical assumptions about signal and noise [5]. Although they require a multidimensional grid search or gradient-based search to find the estimates, they always demonstrate the best performance among the subspace-based methods [6,7,8] and the sparse methods [9,10,11,12], especially when the SNR is not high or the observation is insufficient. The subspace-based methods are computationally attractive. The well-known MUSIC method is equivalent to a large sample realization of the ML method when the sources are uncorrelated, and has a super-resolution compared to beamforming [3]. However, these methods suffer from performance degradation for coherent (or highly correlated) sources if there is no pre-processing technique [6,7].

In the last decade, several high resolution methods have been proposed for DOA estimation by exploiting the techniques in sparse signal representation and compressed sensing [9,13,14,15,16]. The sparse DOA estimation methods are usually solved by convex optimization, such as basis pursuit denoising (BPDN), least absolute shrinkage and selection operator (LASSO) and $l_1$-SVD methods [9]. However, the defect with these methods is that the user parameters are hard to decide and can easily affect the performance. The optimization based on the $l_q,0<q<1$ optimization is a nonconvex relaxation and it is tighter than the previous $l_1$ convex relaxation [15]. Nevertheless, the existing focal underdetermined system solver (FOCUSS) for $l_q$ optimization is an iterative solution. Other important methods are sparse Bayesian learning (SBL) methods [17,18,19]. In the sparse Bayesian framework, the DOA estimation problems are transformed to find the hyperparameters of the assumed sparse prior from the observation data [18]. Although SBL can promote the sparsity of its solution by exploiting the sparse structure of the sparse vector (spatial sparsity), the choice of the prior distribution is still a problem [20,21]. In addition, the prior used in SBL is non-informative with respect to the sources. However, in some scenarios when the prior information is available, the existing SBL method based on the non-informative sparse prior cannot be applied directly. In general, the key issue of the mentioned sparse methods is the grid mismatch problem that affects DOA estimation significantly, although the off-grid sparse methods have been proposed to improve the performance [19,20,21,22,23].

More recently, the gridless sparse methods based on atomic norm were proposed for uniform linear array (ULA) to resolve the grid mismatch problem [24,25]. Because the atomic norm relies on semi-definite programming (SDP), the gridless methods are not efficient for large scale problems [26]. In fact, almost all the sparse methods are more computationally expensive compared to most conventional methods, especially for the gridless sparse methods [10].

In the concept of cognitive radar or sonar, the system performance can also be improved by employing feedback from the fusion center (such as target tracker) and using prior knowledge stored in memory [27,28,29,30,31,32]. Inspired by the cognitive radar, the knowledge-aided signal processing can be introduced in the DOA estimation problem to improve the performance as well. A typical localization or tracking system employs a two-step processing: each sensor first estimates the measurements (such as range, angle and Doppler) of targets, and then the fusion center collects all the measurements for target localization or tracking. Different from the traditional feed-forward processing, each sensor can also employ the feedback knowledge from the target tracker in the concept of “cognitive”, as shown in Figure 1. The prediction of future target positions from the tracker can be used to construct a prior knowledge of the targets. More importantly, this prior knowledge is informative. In this situation, the DOA estimation method can take advantage of the feedback informative prior and improve the estimation performance. Some other sources of prior information may be generated from the constraints on the target location. For example, the roads may provide a constraint on the ground targets location in a ground moving target indicator (GMTI) radar [33,34]. More generally, even for targets not constrained to roads, knowledge of the terrain may be used to generate prior probabilities on target locations [30,35].

Nearly all DOA estimation methods from the literature are only based on the array observation data without exploiting the prior knowledge in processing although Bayesian parameter estimation have been introduced in array signal processing for a long time [6,7,8,9,10,11,12,14,16,19,20,21,22,23,26,36]. The recent work [35] studied the mean square error performance of the Bayesian DOA estimation when the prior information about the target location is available. They showed that the Bayesian estimator outperforms the ML estimator at low SNR. The aim of the paper [35] is to construct an error bound to accurately predict the mean-squared-error of the maximum a posterior (MAP) estimator. Hence, the signal model used in [35] is for one source with single snapshot. In addition, the DOA estimation method is a grid-search MAP. In general, the DOA estimation scenario considered in [35] is simple and direct.

This paper focuses on the Bayesian DOA estimation problem using multiple snapshot observations with the aid of the feedback prior knowledge from the target tracker. Considering the unknown DOAs are assumed to be random variables in the Bayesian framework, two Cramér–Rao lower bounds (CRLBs) are provided to evaluate the performance in this case. One is the Bayesian CRLB (BCRLB) incorporating the prior information, the other is the expected CRLB (ECRLB). Based on the MAP estimator, a grid search method is proposed for a single target DOA estimation. For a general multiple DOAs case, a minimum mean square error (MMSE) estimator based on a Monte Carlo method and an iterative solution based on the MAP estimator are proposed to improve the performance of DOA estimation with the aid of the prior knowledge. The former need to trade off accuracy against computational complexity. The latter requires that the prior knowledge is informative. Simulations show that all the proposed three Bayesian estimators aided by the prior knowledge are better than the methods by only using the observation data.

The main contributions of this work are highlighted as follows:

• The Bayesian DOA estimation problem is formed by combining the multiple snapshots from observation and the prior knowledge from target tracker. By deriving the BCRLB and ECRLB, the relationship between the observation information and prior knowledge is analyzed.
• The MAP estimator is constructed and two methods are proposed. One is a two-step grid search method for a single target DOA case. The other is a gradient-based iterative solution for multiple targets DOA case, which extends the traditional Newton method by incorporating the prior knowledge.
• The MMSE estimator is constructed. Considering the multidimensional integration is difficult to calculate, the Monte Carlo method is proposed to estimate the integration.
• By comparing with the ML estimators and the MUSIC algorithm, the performance improvement achieved by the proposed three Bayesian estimators is demonstrated in different simulation settings.

We need to emphasize that the Bayesian DOA estimation problem has less been considered in the existing DOA estimation literature although Bayesian parameter estimation is a conventional method by incorporating the observation information and the prior information. It is shown that the DOA estimation performance can be improved in low SNR with the aid of the prior information by using the conventional routines in Bayesian framework to solve this Bayesian DOA estimation problem. On the other hand, the proposed three Bayesian DOA estimation methods expands the existing ML DOA estimation method when the prior knowledge is applied.

We start with a brief introduction to DOA estimation methods and the Bayesian DOA estimation aided by the informative prior knowledge in Section 1. In Section 2, the signal model and the prior knowledge are described. In Section 3, two different CRLBs, i.e., BCRLB and ECRLB are presented. In Section 4, three different methods are proposed by exploring two Bayesian estimators: the MAP estimator and MMSE estimator. Then, in Section 5, the performances of the proposed methods are explored using simulated experiments in a single-target DOA case and a multiple-targets DOA case, respectively. Finally, conclusions are made in Section 6.

In this paper, we use capital italic bold letters to represent matrices, and lowercase italic bold letters to represent vectors. The symbols $\mathit{I}$ and $\mathit{O}$ are the identity and zero matrices, and $\mathbf{1}$ and $\mathbf{0}$ are the all-one and zero vectors, with a subscript indicating the size when needed. For a given matrix $\mathit{A}$, ${\mathit{A}}^T$ and ${\mathit{A}}^H$ denote the transpose and conjugate transpose of $\mathit{A}$. ${\left|\right|\mathit{A}\left|\right|}_F$ denotes the Frobenius norm of $\mathit{A}$. For a given vector $\mathit{x}$, ${\left|\right|\mathit{x}\left|\right|}_2$ is the $l_2$ norm. The operations ⊙ and ⊗ are the Hadamard product and Kronecker product. $E[·]$ is the mathematical expectation.

2. Signal Model and Prior Knowledge

2.1. Signal Model for ULA

Let us consider a uniform linear array (ULA) with N omnidirectional sensors. We assume there are K far-field sources located at ${\theta }_1,{\theta }_2,\dots ,{\theta }_K$ impinging on the array sensors. The signal $s_k\left(t\right),k=1,2,\cdots ,K$ radiated from the kth source is narrowband. The output of the ith sensor is $x_i\left(t\right)$ and it is corrupted by an additive noise $n_i\left(t\right),i=1,2,\cdots ,N$. Let $\mathit{x}\left(t\right)={[x_1\left(t\right),x_2\left(t\right),\dots ,x_{{}_N}\left(t\right)]}^T$ and $\mathit{s}\left(t\right)={[s_1\left(t\right),s_2\left(t\right),\dots ,s_{{}_K}\left(t\right)]}^T$. Hence, after demodulation, the signal received by the array at time t is expressed as [7,13]

$$\mathit{x}\left(t\right)=\sum _{k=1}^K\mathit{a}\left({\theta }_k\right)s_k\left(t\right)+\mathit{n}\left(t\right)=\mathit{As}\left(t\right)+\mathit{n}\left(t\right)$$

(1)

for $1⩽t⩽T$ and $\mathit{n}\left(t\right)={[n_1\left(t\right),n_2\left(t\right),\dots ,n_{{}_N}\left(t\right)]}^T$ is the additive complex measurement noise vector assumed to be a zero-mean circular complex Gaussian random vector [3,36,37], i.e., $\mathit{n}\left(t\right)\sim \mathcal{CN}(\mathbf{0},{\mathit{Q}}_n)$. The noise covariance matrix ${\mathit{Q}}_n$ is assumed to be known and it can be estimated using conventional methods based on target-free observations. $\mathit{A}=[\mathit{a}\left({\theta }_1\right),\mathit{a}\left({\theta }_1\right),\cdots ,\mathit{a}\left({\theta }_K\right)]\in {\mathbb{C}}^{N\times K}$ is the array manifold matrix, whose $(i,k)$th element contains the delay and gain information from the kth source to the ith sensor. The kth column $\mathit{a}\left({\theta }_k\right)$ in $\mathit{A}$ is called a steering vector and it is

$$\mathit{a}\left({\theta }_k\right)={[1,e^{\mathrm{i}\frac{2\pi }{\lambda }dsin\left({\theta }_k\right)},\cdots ,e^{\mathrm{i}\frac{2\pi }{\lambda }(N-1)dsin\left({\theta }_k\right)}]}^T,$$

(2)

where ${\theta }_k\in (-{90}^{°},{90}^{°})$, $\lambda$ is the wavelength and d is the sensor spacing. All the complex signal samples $s_k\left(t\right)$ for $t=1,2,\cdots ,T$ and $k=1,2,\cdots ,K$ are also unknown parameters. They can be expressed as

$$s_k\left(t\right)=a_{kt}{e}^{\mathrm{i}{\varphi }_{kt}},$$

(3)

where $a_{kt}$ and ${\varphi }_{kt}$ are the amplitude and phase of the kth source in the time index t, respectively. The traditional DOA estimation methods only use the knowledge of array observations $\mathit{x}\left(t\right)$ to estimate the DOAs $\mathit{\theta }={[{\theta }_1,{\theta }_2,\cdots ,{\theta }_K]}^T$.

Note that although the signal model (1) is for ULA, the proposed methods are easily extended to any other linear arrays including sparse arrays. For a sparse array, only the expressions for the steering vectors $\mathit{a}\left({\theta }_k\right)$ have to be modified. The algorithms that are presented can be applied by replacing the array steering matrix for ULA with that for a sparse array since $\mathit{a}\left({\theta }_k\right)$ is a known function of ${\theta }_k$ for a given array.

2.2. Prior Knowledge from the Tracker

In addition to the observation data, the target tracker in the fusion center will naturally output an expected future position of the target. This can be used to provide a prior knowledge of the target. The block diagram of the target tracker prediction feedback system is shown in Figure 1. The mean and variance of the prior DOA can be obtained from the corresponding target position and estimation covariance matrix by some simple and direct manipulations.

Considering a 2D localization scenario shown in Figure 2, the predicted location of the kth target is assumed as ${\mathit{t}}_k^o={[t_{kx}^o,t_{ky}^o]}^T$ and the estimation covariance matrix is ${\mathit{Q}}_{tk}$. Without loss of generality, assuming the sensor array is placed in the origin of coordinates. Hence, the angle of arrival of the kth target with respect to the array is

$${\phi }_k^o=arctan(t_{ky}^o/t_{kx}^o),$$

(4)

where ${\phi }_k^o\in [0,2\pi )$ and $arctan(·)$ is the four-quadrant inverse tangent. The corresponding variance ${\sigma }_k^2$ can be approximated using first-order Taylor expansion as

$${\mathit{\sigma }}_k^2={\mathit{\beta }}^T{\mathit{Q}}_{tk}\mathit{\beta },\mathit{\beta }={\left[-\frac{t_{ky}^o}{\left|\right|{\mathit{t}}_k^o{\left|\right|}^2},\frac{t_{kx}^o}{\left|\right|{\mathit{t}}_k^o{\left|\right|}^2}\right]}^T.$$

(5)

Given the intersection angle between the array and the x-axis, the mean DOA of the kth target ${\theta }_k^o$ is obtained from ${\phi }_k^o$. Because the predicted target position error is usually modeled by Gaussian distribution, the prior knowledge of DOA that was constructed from the predicted target positions and position errors is approximated by Gaussian distribution for simplicity, i.e., ${\theta }_k\sim \mathcal{N}({\theta }_k^o,{\sigma }_k^2)$. Hence, the prior probability density function (PDF) for the DOA vector $\mathit{\theta }$ is

$$p\left(\mathit{\theta }\right)=\frac{1}{{\left(2\pi \right)}^{\frac{K}{2}}{\left|{\mathit{Q}}_{\theta }\right|}^{\frac{1}{2}}}exp\left[-\frac{1}{2}{(\mathit{\theta }-{\mathit{\theta }}^o)}^T{\mathit{Q}}_{\theta }^{-1}(\mathit{\theta }-{\mathit{\theta }}^o)\right],$$

(6)

where ${\mathit{\theta }}^o={[{\theta }_1^o,{\theta }_2^o,\cdots ,{\theta }_K^o]}^T$ and ${\mathit{Q}}_{\theta }=\mathrm{diag}\left([{\sigma }_1^2,{\sigma }_2^2,\cdots ,{\sigma }_K^2]\right)$. In this paper, we would like to obtain $\mathit{\theta }$ using observation of experimental data $\mathit{x}\left(t\right)$ together with the prior knowledge $p\left(\mathit{\theta }\right)$.

3. Cramér–Rao Lower Bounds

The CRLB provides a lower bound on the error covariance matrix for any unbiased estimator. This section evaluates the CRLB of the unknown DOA vector $\mathit{\theta }$ that has prior knowledge. Because of the existence of the unknown parameters $s_k\left(t\right)$, all the unknown parameters can be defined as

$$\mathit{\alpha }={[{\mathit{\theta }}^T,{\mathit{a}}^T,{\mathit{\varphi }}^T]}^T,$$

(7)

where $\mathit{a}={[a_{11},a_{12},\cdots ,a_{1T},a_{21},\cdots ,a_{KT}]}^T$ and $\mathit{\varphi }={[{\varphi }_{11},{\varphi }_{12},\cdots ,{\varphi }_{1T},{\varphi }_{21},\cdots ,{\varphi }_{KT}]}^T$. Note that $\mathit{\alpha }$ is a hybrid real vector including random parameter $\mathit{\theta }$ and the deterministic parameters $\mathit{a}$ and $\mathit{\varphi }$.

For the sake of simplicity, let us transform the observation Equation (1) to a vector form, which is more compact. By defining the observations $\mathit{x}\left(t\right)$ as a vector $\overline{\mathit{x}}={[{\mathit{x}}^T\left(1\right),{\mathit{x}}^T\left(2\right),\cdots ,{\mathit{x}}^T\left(T\right)]}^T$, (1) can be transformed into

$$\overline{\mathit{x}}=\overline{\mathit{A}}\overline{\mathit{s}}+\overline{\mathit{n}},$$

(8)

where $\overline{\mathit{A}}={\mathit{I}}_T\otimes \mathit{A}$, $\overline{\mathit{s}}={[{\mathit{s}}^T\left(1\right),{\mathit{s}}^T\left(2\right),\cdots ,{\mathit{s}}^T\left(T\right)]}^T$ and $\overline{\mathit{n}}={[{\mathit{n}}^T\left(1\right),{\mathit{n}}^T\left(2\right),\cdots ,{\mathit{n}}^T\left(T\right)]}^T$ with zero mean and covariance matrix ${\mathit{Q}}_{\overline{n}}={\mathit{I}}_T\otimes {\mathit{Q}}_n$. Hence, by using the complex Gaussian PDF, the likelihood function or the conditional probability distribution function of $\overline{\mathit{x}}$ given $\mathit{\theta }$ and $\overline{\mathit{s}}$ is [38]

$$p\left(\overline{\mathit{x}}\right|\mathit{\theta },\overline{\mathit{s}})=\frac{1}{{\pi }^{NT}\left|{\mathit{Q}}_{\overline{n}}\right|}exp\left[-{(\overline{\mathit{x}}-\overline{\mathit{A}}\overline{\mathit{s}})}^H{\mathit{Q}}_{\overline{n}}^{-1}(\overline{\mathit{x}}-\overline{\mathit{A}}\overline{\mathit{s}})\right]$$

(9)

Finally, the the joint PDF for $\overline{\mathit{x}}$ and $\mathit{\theta }$ is

$$p(\overline{\mathit{x}},\mathit{\theta }|\overline{\mathit{s}})=p\left(\overline{\mathit{x}}\right|\mathit{\theta },\overline{\mathit{s}})p\left(\mathit{\theta }\right).$$

(10)

3.1. Bayesian CRLB

In Bayesian estimation, the performance of any estimator can be bounded by the Bayesian CRLB (BCRLB). The Bayesian information matrix for $\mathit{\alpha }$ from the observations $\overline{\mathit{x}}$ and prior knowledge $p\left(\mathit{\theta }\right)$ is [39]

$${\mathit{J}}_b=E_{\overline{\mathit{x}},\mathit{\theta }}\left[\frac{\partial lnp(\overline{\mathit{x}},\mathit{\theta }|\overline{\mathit{s}})}{\partial \mathit{\alpha }}\frac{\partial lnp(\overline{\mathit{x}},\mathit{\theta }|\overline{\mathit{s}})}{\partial {\mathit{\alpha }}^T}\right],$$

(11)

whose inverse is the BCRLB of $\mathit{\alpha }$ and the upper left block $K\times K$ matrix gives the BCRLB for $\mathit{\theta }$. The expectation in (11) is taken jointly over $\overline{\mathit{x}}$ and $\mathit{\theta }$. After substituting (10) in (11), it yields

$${\mathit{J}}_b={\mathit{J}}_d+{\mathit{J}}_p,$$

(12)

where

$${\mathit{J}}_d=E_{\overline{\mathit{x}},\mathit{\theta }}\left[\frac{\partial lnp\left(\overline{\mathit{x}}\right|\mathit{\theta },\overline{\mathit{s}})}{\partial \mathit{\alpha }}\frac{\partial lnp\left(\overline{\mathit{x}}\right|\mathit{\theta },\overline{\mathit{s}})}{\partial {\mathit{\alpha }}^T}\right],$$

(13)

$${\mathit{J}}_p=E_{\mathit{\theta }}\left[\frac{\partial lnp\left(\mathit{\theta }\right)}{\partial \mathit{\alpha }}\frac{\partial lnp\left(\mathit{\theta }\right)}{\partial {\mathit{\alpha }}^T}\right].$$

(14)

Roughly speaking, ${\mathit{J}}_d$ represents the information from the observations and ${\mathit{J}}_p$ represents the information from prior. When the SNR is high and the number of observations is large, the BCRLB will be dominated by ${\mathit{J}}_d$, i.e., the observation information. On the other hand, ${\mathit{J}}_p$ will affect the BCRLB in the situations with limited number of observations or low SNR.

First let us focus on ${\mathit{J}}_d$. Note that ${\mathit{J}}_d$ is similar to the traditional Fisher information matrix for deterministic parameters. Let us define the Fisher information matrix ${\mathit{J}}_f$ as

$${\mathit{J}}_f=E_{\overline{\mathit{x}}}\left[\frac{\partial lnp\left(\overline{\mathit{x}}\right|\mathit{\theta },\overline{\mathit{s}})}{\partial \mathit{\alpha }}\frac{\partial lnp\left(\overline{\mathit{x}}\right|\mathit{\theta },\overline{\mathit{s}})}{\partial {\mathit{\alpha }}^T}\right],$$

(15)

and ${\mathit{J}}_d$ is related to ${\mathit{J}}_f$ by

$${\mathit{J}}_d=E_{\mathit{\theta }}\left[{\mathit{J}}_f\right].$$

(16)

Hence, ${\mathit{J}}_d$ can be obtained from ${\mathit{J}}_f$ by taking expectation of the components over $\mathit{\theta }$. However, this expectation is hard to derived analytically due to the involved expression of ${\mathit{J}}_f$, and therefore has to be evaluated numerically in general using a Monte Carlo simulation. Now, the focus is on deriving ${\mathit{J}}_f$. Considering the Fisher information matrix ${\mathit{J}}_f$ involves the estimation of a real-valued parameters-vector $\mathit{\alpha }$ from complex circular Gaussian vector $\overline{\mathit{x}}$, ${\mathit{J}}_f$ is given by [38]

$${\mathit{J}}_f=2\mathrm{Re}\left(\frac{\partial {\left(\overline{\mathit{A}}\overline{\mathit{s}}\right)}^T}{\partial \mathit{\alpha }}{\mathit{Q}}_{\overline{n}}^{-1}\frac{\partial \overline{\mathit{A}}\overline{\mathit{s}}}{\partial {\mathit{\alpha }}^T}\right),$$

(17)

where

$$\frac{\partial \overline{\mathit{A}}\overline{\mathit{s}}}{\partial {\mathit{\alpha }}^T}=\left[\frac{\partial \overline{\mathit{A}}\overline{\mathit{s}}}{\partial {\mathit{\theta }}^T},\frac{\partial \overline{\mathit{A}}\overline{\mathit{s}}}{\partial {\mathit{a}}^T},\frac{\partial \overline{\mathit{A}}\overline{\mathit{s}}}{\partial {\mathit{\varphi }}^T}\right],$$

(18)

and the corresponding elements are in Appendix A.

Next, we will focus on ${\mathit{J}}_p$. Putting (6) in (14) and taking the expectation over $\mathit{\theta }$ produce

$${\mathit{J}}_p=\left[\begin{array}{cc}{\mathit{Q}}_{\theta }^{-1}& \mathit{O}\\ \mathit{O}& {\mathit{O}}_{2KT}\end{array}\right].$$

(19)

3.2. Expected CRLB

The Fisher information matrix ${\mathit{J}}_f$ depends on the DOA vector $\mathit{\theta }$. In the traditional DOA estimation when $\mathit{\theta }$ is a deterministic parameter with true value equals to $\overline{\mathit{\theta }}$, ${\mathit{J}}_f$ is also a deterministic matrix by replacing $\mathit{\theta }$ by $\overline{\mathit{\theta }}$. When $\mathit{\theta }$ is a random parameter, a similar CRLB can be constructed directly by taking the expectation of ${\mathit{J}}_f^{-1}$

$${\mathit{C}}_e=E_{\mathit{\theta }}\left[{\mathit{J}}_f^{-1}\right],$$

(20)

which is also called the expected CRLB (ECRLB). We will gain some insights when using Jensen’s inequality to (20)

$${\mathit{C}}_e=E_{\mathit{\theta }}\left[{\mathit{J}}_f^{-1}\right]\ge {\left[E_{\mathit{\theta }}\left[{\mathit{J}}_f\right]\right]}^{-1}={\mathit{J}}_d^{-1},$$

(21)

Because BCRLB is approximately equal to ${\mathit{J}}_d^{-1}$ when there are enough observations, i.e., ${lim}_{T\to \infty }{\mathit{J}}_b^{-1}={\mathit{J}}_d^{-1}$, the ECRLB is greater than the BCRLB asymptotically [39]. When the SNR is high and the number of observations is large, ECRLB is a tight bound for the unbiased estimator [35]. Some other bounds, such as the Barankin bound, Ziv–Zakai bound, and Weiss–Weinstein bound, are proposed to capture the performance in the region of low SNR  [39,40]. These bounds are also called the large error bounds to better fit the threshold effect. The study of the lower bounds for the Bayesian DOA estimation problem is outside the scope of this work. Instead, the presented two bounds are intended as a bench mark to evaluate the performance of the proposed methods.

4. Bayesian DOA Estimation Methods with Prior Knowledge

In this section, we will concentrate on two principal Bayesian estimators for DOA estimation, i.e., MAP estimator and MMSE estimator. MAP estimator chooses ${\mathit{\theta }}_{\mathrm{map}}$ that maximize the posterior PDF. MMSE estimator ${\mathit{\theta }}_{\mathrm{mmse}}$ is defined to minimize the MSE (or Bayesian MSE) when averaged over all realization of $\mathit{\theta }$ and $\overline{\mathit{x}}$. The MMSE estimator is the mean of the posterior PDF or ${\mathit{\theta }}_{\mathrm{mmse}}=E\left(\mathit{\theta }\right|\overline{\mathit{x}})$ [38]. However, it is hard to find closed-form expressions for these two estimators due to the multidimensional maximization for the MAP estimator and the multidimensional integration for the MMSE estimator.

4.1. MAP Estimator Using Grid Search

Considering $p\left(\overline{\mathit{x}}\right|\mathit{\theta },\overline{\mathit{s}})$ in (9) includes the unknown parameters $\overline{\mathit{s}}$ (or $s_k\left(t\right)=a_{kt}{e}^{\mathrm{i}{\varphi }_{kt}}$) and they are hard to handle directly, they will be expressed by $\mathit{\theta }$ before applying Bayes’ theorem. Taking a partial derivative of $lnp\left(\overline{\mathit{x}}\right|\mathit{\theta },\overline{\mathit{s}})$ with respect to $\overline{\mathit{s}}$ and letting it go to zero gives

$$\frac{\partial lnp\left(\overline{\mathit{x}}\right|\mathit{\theta },\overline{\mathit{s}})}{\partial \overline{\mathit{s}}}=\mathit{0}.$$

(22)

Substituting (9) in (22) and solving for $\overline{\mathit{s}}$ yields

$$\overline{\mathit{s}}={\left({\overline{\mathit{A}}}^H{\mathit{Q}}_{\overline{n}}^{-1}\overline{\mathit{A}}\right)}^{-1}{\overline{\mathit{A}}}^H{\mathit{Q}}_{\overline{n}}^{-1}\overline{\mathit{x}}.$$

(23)

Putting (23) in (9) yields the conditional probability distribution function of $\overline{\mathit{x}}$ given $\mathit{\theta }$, expressed by $p\left(\overline{\mathit{x}}\right|\mathit{\theta })$.

$$p\left(\overline{\mathit{x}}\right|\mathit{\theta })=\frac{1}{{\pi }^{NT}\left|{\mathit{Q}}_{\overline{n}}\right|}exp\left[-{(\overline{\mathit{x}}-\Phi \overline{\mathit{x}})}^H{\mathit{Q}}_{\overline{n}}^{-1}(\overline{\mathit{x}}-\Phi \overline{\mathit{x}})\right],$$

(24)

where $\Phi =\overline{\mathit{A}}{\left({\overline{\mathit{A}}}^H{\mathit{Q}}_{\overline{n}}^{-1}\overline{\mathit{A}}\right)}^{-1}{\overline{\mathit{A}}}^H{\mathit{Q}}_{\overline{n}}^{-1}$. Note that $p\left(\overline{\mathit{x}}\right|\mathit{\theta })$ only contains the unknown $\mathit{\theta }$. By using the Bayes’ theorem, the posterior probability of the DOA vector $\mathit{\theta }$ given the observation $\overline{\mathit{x}}$ is

$$p\left(\mathit{\theta }\right|\overline{\mathit{x}})=\frac{p\left(\overline{\mathit{x}}\right|\mathit{\theta })p\left(\mathit{\theta }\right)}{\int p\left(\overline{\mathit{x}}\right|\mathit{\theta })p\left(\mathit{\theta }\right)d\mathit{\theta }}.$$

(25)

Because the denominator in (25) is $p\left(\overline{\mathit{x}}\right)=\int p\left(\overline{\mathit{x}}\right|\mathit{\theta })p\left(\mathit{\theta }\right)d\mathit{\theta }$ that is not related with $\mathit{\theta }$, maximizing the posterior probability $p\left(\mathit{\theta }\right|\overline{\mathit{x}})$ is equivalent to maximizing the numerator in (25). Hence,

$$\begin{array}{cc}\hfill {\mathit{\theta }}_{\mathrm{map}}& =argmaxp\left(\mathit{\theta }\right|\overline{\mathit{x}})=argmaxp\left(\overline{\mathit{x}}\right|\mathit{\theta })p\left(\mathit{\theta }\right)\hfill \\ \hfill & =argmax\left[lnp\left(\overline{\mathit{x}}\right|\mathit{\theta })+lnp\left(\mathit{\theta }\right)\right].\hfill \end{array}$$

(26)

For (26), using (6) and (24), and ignoring the constant terms produce the final MAP estimator

$$\begin{array}{cc}\hfill {\mathit{\theta }}_{\mathrm{map}}& =argmax\left[T_{\mathrm{ml}}-\frac{1}{2}{(\mathit{\theta }-{\mathit{\theta }}^o)}^T{\mathit{Q}}_{\theta }^{-1}(\mathit{\theta }-{\mathit{\theta }}^o)\right],\hfill \end{array}$$

(27)

where

$$T_{\mathrm{ml}}={\overline{\mathit{x}}}^H{\mathit{Q}}_{\overline{n}}^{-1}\overline{\mathit{A}}{\left({\overline{\mathit{A}}}^H{\mathit{Q}}_{\overline{n}}^{-1}\overline{\mathit{A}}\right)}^{-1}{\overline{\mathit{A}}}^H{\mathit{Q}}_{\overline{n}}^{-1}\overline{\mathit{x}}$$

(28)

is the maximum-likelihood test statistic. Hence, the ML estimation for DOA is

$${\mathit{\theta }}_{\mathrm{ml}}=argmax{T}_{\mathrm{ml}}.$$

(29)

Note that $T_{\mathrm{ml}}$ is a highly nonlinear and nonconvex function with respect to $\mathit{\theta }$ and it contains many local minimum. Hence, both the MAP estimator (27) and the ML estimator (29) require a multidimensional grid search method or a gradient-based methods (or other numerical iterative methods). Moreover, the MAP estimator in [35] is a special case of (27) for single snapshot. In this subsection, we focus on the grid search method, but the computational burden increases exponentially with the dimension K (number of targets). Nevertheless, it is still a good method for single target DOA estimation. A possible implementation of the grid search for single target is in the Algorithm 1.

Algorithm 1 MAP estimator for single target using grid search

Input: The observation vector $\overline{\mathit{x}}$ and noise covariance matrix $\overline{\mathit{Q}}$, the prior knowledge ${\theta }^o$ and ${\sigma }_{\theta }^2$ Output: An estimate ${\theta }_{\mathrm{map}}$     1: Create a rough grid of source locations ${\theta }_i^{\left(1\right)}$, $i=1,2,\cdots ,N_{\theta 1}$ with step size $d_{\theta }$     2: Compute the value using the statistic in (27) and find the maximum value ${\theta }_f^{\left(1\right)}$     3: Obtain a refined grid ${\theta }_i^{\left(2\right)}$, $i=1,2,\cdots ,N_{\theta 2}$ in the region $[{\theta }_f^{\left(1\right)}-d_{\theta }/2,{\theta }_f^{\left(1\right)}+d_{\theta }/2]$     4: Compute the value using the statistic in (27) and find the maximum value ${\theta }_{\mathrm{map}}$     5: return result

4.2. MMSE Estimator Using Monte Carlo Method

The MMSE estimator for DOA estimation is

$$\begin{array}{c}\hfill {\mathit{\theta }}_{\mathrm{mmse}}=\int p\left(\mathit{\theta }\right|\overline{\mathit{x}})\mathit{\theta }d\mathit{\theta }=\int \frac{\mathit{\theta }p\left(\overline{\mathit{x}}\right|\mathit{\theta })p\left(\mathit{\theta }\right)}{\int p\left(\overline{\mathit{x}}\right|\mathit{\theta })p\left(\mathit{\theta }\right)d\mathit{\theta }}d\mathit{\theta }.\end{array}$$

(30)

Because the denominator is not related with $\mathit{\theta }$, it can be expressed as

$$\begin{array}{c}\hfill {\mathit{\theta }}_{\mathrm{mmse}}=\frac{\int \mathit{\theta }p\left(\overline{\mathit{x}}\right|\mathit{\theta })p\left(\mathit{\theta }\right)d\mathit{\theta }}{\int p\left(\overline{\mathit{x}}\right|\mathit{\theta })p\left(\mathit{\theta }\right)d\mathit{\theta }}=\frac{E_{\mathit{\theta }}\left[\mathit{\theta }p\left(\overline{\mathit{x}}\right|\mathit{\theta })\right]}{E_{\mathit{\theta }}\left[p\left(\overline{\mathit{x}}\right|\mathit{\theta })\right]}.\end{array}$$

(31)

Both the numerator and denominator in (31) involve multidimensional integration, which makes it difficult to derive the closed-form expression. A numerical solution that resorts to the Monte Carlo method will be proposed. Supposing there are $N_m$ samples from the prior distribution $p\left(\mathit{\theta }\right)$, the expected values in (31) can be approximated by the finite sum of the samples, i.e.,

$${\mathit{\theta }}_{\mathrm{mmse}}\approx \frac{{\sum }_1^{N_m}{\mathit{\theta }}_{i}p\left(\overline{\mathit{x}}\right|{\mathit{\theta }}_i)}{{\sum }_1^{N_m}p\left(\overline{\mathit{x}}\right|{\mathit{\theta }}_i)},$$

(32)

where ${\mathit{\theta }}_i,i=1,2,\cdots ,N_m$ are generated from the known prior $p\left(\mathit{\theta }\right)$ and it is direct. The DOA estimation is more accurate when increasing $N_m$, but on the other hand, the computational complexity is increasing as well. Hence, there is a trade-off between the accuracy and the computational complexity. Putting (23) and (9) in (32) and rearranging yield

$${\mathit{\theta }}_{\mathrm{mmse}}\approx \frac{{\sum }_1^{N_m}{\mathit{\theta }}_{i}exp[-{\overline{\mathit{x}}}^H\mathit{R}\left({\mathit{\theta }}_i\right)\overline{\mathit{x}}]}{{\sum }_1^{N_m}exp[-{\overline{\mathit{x}}}^H\mathit{R}\left({\mathit{\theta }}_i\right)\overline{\mathit{x}}]},$$

(33)

where

$$\mathit{R}\left({\mathit{\theta }}_i\right)={\mathit{Q}}_{\overline{n}}^{-1}-{\mathit{Q}}_{\overline{n}}^{-1}\overline{\mathit{A}}{\left({\overline{\mathit{A}}}^H{\mathit{Q}}_{\overline{n}}^{-1}\overline{\mathit{A}}\right)}^{-1}{\overline{\mathit{A}}}^H{\mathit{Q}}_{\overline{n}}^{-1}$$

(34)

when substituting ${\mathit{\theta }}_i$. Note that when the observation noise is Gaussian, the MMSE estimator is equivalent to the MAP estimator and the only difference is the calculation method. The method is summarized in Algorithm 2.

Algorithm 2 MMSE estimator for multiple targets using the Monte Carlo method

Input: The observation vector $\overline{\mathit{x}}$ and noise covariance matrix $\overline{\mathit{Q}}$, the prior knowledge ${\mathit{\theta }}^o$ and ${\mathit{Q}}_{\theta }$ Output: An estimate ${\mathit{\theta }}_{\mathrm{mmse}}$     1: Generate $N_m$ samples ${\mathit{\theta }}_i$ from the prior Gaussian distribution $\mathcal{N}({\mathit{\theta }}^o,{\mathit{Q}}_{\theta })$     2: Compute $\mathit{R}\left({\mathit{\theta }}_i\right)$ for each ${\mathit{\theta }}_i$ using (34) and ${\mathit{\theta }}_{\mathrm{mmse}}$ using (33)     3: return result

4.3. Iterative Solution Based on the MAP Estimator

4.3.1. Gradient-Based Iterative Solution

When the array observation noise $\mathit{n}\left(t\right)$ is homoscedastic, i.e., ${\mathit{Q}}_n={\sigma }_n^2\mathit{I}$, the likelihood function (9) can be expressed as

$$p\left(\overline{\mathit{x}}\right|\mathit{\theta },\overline{\mathit{s}})=\frac{1}{{\pi }^{NT}{\sigma }_n^{2NT}}exp\left[-\frac{1}{{\sigma }_n^2}\left|\right|\mathit{X}-{\mathit{AS}\left|\right|}_F^2\right],$$

(35)

where $\mathit{X}=\left[\mathit{x}\right(1),\mathit{x}(2),\cdots ,\mathit{x}(T\left)\right]$ and $\mathit{S}=\left[\mathit{s}\right(1),\mathit{s}(2),\cdots ,\mathit{s}(T\left)\right]$. Note that $p\left(\overline{\mathit{x}}\right|\mathit{\theta },\overline{\mathit{s}})$ includes the unknown parameters $\mathit{a}$ and $\mathit{\varphi }$ through $\mathit{S}$. The similar technique in (22) will be applied to express $\mathit{S}$ by $\mathit{A}$ and $\mathit{X}$.

$$\mathit{S}={\mathit{A}}^{†}\mathit{X}={\left({\mathit{A}}^H\mathit{A}\right)}^{-1}{\mathit{A}}^H\mathit{X}.$$

(36)

Substituting (36) in (35) yields the conditional probability distribution function of $\overline{\mathit{x}}$ given $\mathit{\theta }$

$$p\left(\overline{\mathit{x}}\right|\mathit{\theta })=\frac{1}{{\pi }^{NT}{\sigma }_n^{2NT}}exp\left[-\frac{1}{{\sigma }_n^2}\left|\right|\mathit{X}-{\Psi \mathit{X}\left|\right|}_F^2\right],$$

(37)

where $\Psi =\mathit{A}{\left({\mathit{A}}^H\mathit{A}\right)}^{-1}{\mathit{A}}^H$. The MAP estimator in (26) is equivalent to minimizing the following negative log-likelihood function:

$$\begin{array}{c}\hfill {\mathit{\theta }}_{\mathrm{map}}=argmin\left[-lnp\left(\overline{\mathit{x}}\right|\mathit{\theta })-lnp\left(\mathit{\theta }\right)\right].\end{array}$$

(38)

Putting (37) and (6) in (38) and ignoring the constant terms give the following cost function for minimization:

$$\begin{array}{c}\hfill J_{\mathrm{map}}=\frac{1}{{\sigma }_n^2}\left|\right|\mathit{X}-{\Psi \mathit{X}\left|\right|}_F^2+\frac{1}{2}{(\mathit{\theta }-{\mathit{\theta }}^o)}^T{\mathit{Q}}_{\theta }^{-1}(\mathit{\theta }-{\mathit{\theta }}^o),\end{array}$$

(39)

where $J_{\mathrm{ml}}\doteq \frac{1}{{\sigma }_n^2}\left|\right|\mathit{X}-{\Psi \mathit{X}\left|\right|}_F^2$ is corresponding to the cost function in the deterministic maximum likelihood (DML) estimation [3]. $J_{\mathrm{ml}}$ can be further expressed as

$$J_{\mathrm{ml}}=\frac{1}{{\sigma }_n^2}\left|\right|\mathit{X}-\mathit{A}{\mathit{A}}^{†}\mathit{X}{\left|\right|}_F^2=\frac{T}{{\sigma }_n^2}\mathrm{tr}\left({\mathbb{P}}_{\mathit{A}}^{\perp }\mathit{R}\right),$$

(40)

where

$$\begin{array}{c}\hfill {\mathbb{P}}_{\mathit{A}}^{\perp }=\mathit{I}-\mathit{A}{\mathit{A}}^{†},\\ \hfill \mathit{R}=\frac{1}{T}\mathit{X}{\mathit{X}}^H.\end{array}$$

(41)

Note that the second term in (39) is a quadratic function related to $\mathit{\theta }$, but the first term $J_{\mathrm{ml}}$ is nonlinear and nonconvex. It can be approximated by using second order Taylor expansion. The second order approximation of $J_{\mathrm{ml}}\left(\mathit{\theta }\right)$ at ${\mathit{\theta }}^{\left(j\right)}$ is

$$J_{\mathrm{ml}}\left(\mathit{\theta }\right)\simeq J_{\mathrm{ml}}\left({\mathit{\theta }}^{\left(j\right)}\right)+\nabla J_{\mathrm{ml}}^{\left(j\right)T}\Delta \mathit{\theta }+\frac{1}{2}\Delta {\mathit{\theta }}^T{\mathit{H}}^{\left(j\right)}\Delta \mathit{\theta },$$

(42)

where $\Delta \mathit{\theta }=\mathit{\theta }-{\mathit{\theta }}^{\left(j\right)}$. $\nabla J_{\mathrm{ml}}\left({\mathit{\theta }}^{\left(j\right)}\right)$ and ${\mathit{H}}^{\left(j\right)}$ are the gradient vector and the Hessian matrix evaluated at ${\mathit{\theta }}^{\left(j\right)}$, respectively. They are expressed as [5]

$$\nabla J_{\mathrm{ml}}^{\left(j\right)}={\left.-\frac{2T}{{\sigma }_n^2}\mathrm{Re}\left[\mathrm{diag}\left({\mathit{A}}^{†}\mathit{R}{\mathbb{P}}_{\mathit{A}}^{\perp }\mathit{D}\right)\right]\right|}_{\mathit{\theta }={\mathit{\theta }}^{\left(j\right)}},$$

(43)

$${\mathit{H}}^{\left(j\right)}={\left.\frac{2T}{{\sigma }_n^2}\mathrm{Re}\left[\left({\mathit{D}}^H{\mathbb{P}}_{\mathit{A}}^{\perp }\mathit{D}\right)\odot {\left({\mathit{A}}^{†}\mathit{R}{\mathit{A}}^{†H}\right)}^T\right]\right|}_{\mathit{\theta }={\mathit{\theta }}^{\left(j\right)}},$$

(44)

where

$$\mathit{D}=\left[\frac{\partial \mathit{a}\left({\theta }_1\right)}{\partial {\theta }_1},\frac{\partial \mathit{a}\left({\theta }_2\right)}{\partial {\theta }_2},\cdots ,\frac{\partial \mathit{a}\left({\theta }_K\right)}{\partial {\theta }_K}\right],$$

(45)

and $\frac{\partial \mathit{a}\left({\theta }_k\right)}{\partial {\theta }_k}$ is in Appendix A. After substituting $\mathit{\theta }=\Delta \mathit{\theta }+{\mathit{\theta }}^{\left(j\right)}$ and (42) in (39), minimizing $J_{\mathrm{map}}$ with respect to $\Delta \mathit{\theta }$ yields the update

$$\Delta {\mathit{\theta }}^{\left(j\right)}=-{\left({\mathit{H}}^{\left(j\right)}+{\mathit{Q}}_{\theta }^{-1}\right)}^{-1}\left(\nabla J_{\mathrm{ml}}^{\left(j\right)}+{\mathit{Q}}_{\theta }^{-1}({\mathit{\theta }}^{\left(j\right)}-{\mathit{\theta }}^o)\right).$$

(46)

Therefore, the estimate is iteratively calculated as

$${\mathit{\theta }}^{(j+1)}={\mathit{\theta }}^{\left(j\right)}+{\mu }_j\Delta {\mathit{\theta }}^{\left(j\right)}.$$

(47)

where ${\mu }_j$ is the step length that can be set as a constant or selected in each step using the Wolfe conditions. The iterative solution (47) can be regarded as a variant of the Newton method by incorporating the prior knowledge. We can gain some insight from (43) to (46). When the number of observations T is large enough and the observation noise variance ${\sigma }^2$ is not large (SNR is high), the update will be dominated by $\nabla J_{\mathrm{ml}}^{\left(j\right)}$ and ${\mathit{H}}^{\left(j\right)}$. The iterative solution will approximate to the traditional Newton method that is derived from the DML. On the other hand, when T is insufficient or the SNR is low, the prior knowledge will affect the solution. The analysis here is consistent with that in the Bayesian CRLB. Finally, the iterative solution for the MAP estimator is summarized in the Algorithm 3.

Algorithm 3 Iterative solution based on the MAP estimator

Input: The observation matrix $\mathit{X}$ and the noise variance ${\sigma }_n^2$, the prior knowledge ${\mathit{\theta }}^o$ and ${\mathit{Q}}_{\theta }$, initial guess ${\mathit{\theta }}^{\left(0\right)}$ and the error tolerance $ϵ$ Output: An estimate ${\mathit{\theta }}_{\mathrm{map}}$     1: Initialize $\nabla J_{\mathrm{ml}}^{\left(j\right)}$ and ${\mathit{H}}^{\left(j\right)}$ using (43) and (44), iteration count $j=0$     2: whilej is less than a given number of iterations do     3:      Compute $\Delta {\mathit{\theta }}^{\left(j\right)}$ using (46)     4:      Update ${\mathit{\theta }}^{(j+1)}$ using (47)     5:      if $\left|\right|\Delta {\mathit{\theta }}^{\left(j\right)}{\left|\right|}_2$ is less than the error tolerance $ϵ$ then     6:          return ${\mathit{\theta }}^{(j+1)}$     7:      else     8:          Increment j, and compute $\nabla J_{\mathrm{ml}}^{\left(j\right)}$ and ${\mathit{H}}^{\left(j\right)}$ using (43) and (44)     9: return result

4.3.2. Initialization for DOA

The initial guess for DOA can be obtained in two ways. One way is from the the prior knowledge by letting ${\mathit{\theta }}^{\left(0\right)}={\mathit{\theta }}^o$ when ${\sigma }_k^2$ is not large. The MUSIC algorithm (or any other DOA estimation method) could also be used as an initial guess for uncorrelated signals when the observations are enough with high SNR. In this paper, the prior knowledge is informative (${\sigma }_k^2$ is not large) and it is enough to provide a close initial guess.

4.4. Complexity Analysis

The complexities of the proposed estimators are summarized in

Table 1. Complexity of the estimators.

Estimator	Complexity *
MAP-grid search	$O\left(N_{\theta }^K(K^3{T}^3+K^2{T}^{3}N)\right)$
MMSE	$O\left(N_m(K^3{T}^3+K^2{T}^{3}N)\right)$
MAP-iterative	$O\left(M(K^3+N^3+N^2{T}^2)\right)$

* $N_{\theta }$ is the number of grid, $N_m$ is the number of samples in MMSE and M is the number of iterations in MAP-iterative.

. All the algorithms depend on the number of targets, the number of sensors and the number of snapshots. The complexity of MAP-grid search increases exponentially with the number of targets. Algorithm 1 is for single target and the complexity is $O\left((N_{\theta 1}+N_{\theta 2})\left(T^{3}N\right)\right)$.

5. Simulation

The performance of the proposed Bayesian estimator, i.e., MAP estimator using grid search, MMSE estimator using the Monte Carlo method and the iterative solution based on MAP will be evaluated in this section. For comparison, the MUSIC method and the corresponding ML estimators (grid search-based ML and iterative solution based ML) are also included in the simulation. Note that the ML estimators are always regarded as the benchmark of the methods including the mentioned subspace-based methods and the sparse methods. The existing Bayesian estimation method [18,23] in DOA estimation is sparse Bayesian learning or sparse Bayesian inference, and they also reply on the array observation data. Hence the performance of these methods is also bounded by the ML estimators.

The root mean square error (RMSE) of DOAs is defined as

$$\mathrm{RMSE}\left(\mathit{\theta }\right)=\sqrt{\frac{1}{LK}\left(\sum _{l=1}^L{\parallel {\widehat{\mathit{\theta }}}_{\left(l\right)}-{\mathit{\theta }}_{\left(l\right)}\parallel }_2^2\right)},$$

(48)

where L is the number of trails, ${\mathit{\theta }}_{\left(l\right)}$ and ${\widehat{\mathit{\theta }}}_{\left(l\right)}$ are the true random sample and estimated DOA in the lth trail, respectively. In the simulation part, because $\mathit{\theta }$ is a random vector with PDF $p\left(\mathit{\theta }\right)$, the true DOA values ${\mathit{\theta }}_{\left(l\right)}$ are generated randomly from the given distribution in each trial. The RMSE are calculated from ${\widehat{\mathit{\theta }}}_{\left(l\right)}$ (denote estimated DOA in the lth trail) by averaging over L trails. The DOA errors are evaluated in degree (Deg.) in the results.

We consider a uniform linear array of $N=10$ sensors separated by half a wavelength of the actual narrowband source signals. In the MMSE estimator, the number of samples is $N_m=$ 2000. The error tolerance in the MAP-based iterative solution is $ϵ={10}^{-6}$.

5.1. Single Target DOA Estimation

Because the grid search MAP estimator is proposed for single target, the DOA estimation accuracy for a single target will be evaluated at first. The prior knowledge of the target is ${\theta }^o={30}^{°}$ and ${\sigma }_{\theta }=2^{°}$ unless specified otherwise. The number of snapshots is 30.

Figure 3 directly displays the histograms of the estimation errors for each method obtained from 10,000 Monte Carlo simulations. The SNR is $-15$ dB. The errors of three Bayesian estimators are in the region $[-{10}^{°},{10}^{°}]$, but the errors of the ML-grid search and the MUSIC method are in a wider region $[-{90}^{°},{50}^{°}]$ with numerous outliers. Hence, all three Bayesian estimators outperform the estimators without prior knowledge in this scenario with low SNR, and the histograms of them are similar. The two ML estimators are better than the MUSIC method. However, the ML-iterative solution is better than the ML-grid search because the initial guess in the ML-iterative solution is ${\theta }^o$, which is close to the true value.

Figure 4 shows the RMSE of the proposed estimators at different SNR. When the SNR is high, BCRLB is slightly less than the ECRLB and there is a tiny gap between them. This gap becomes wider when the SNR is extremely low. When the SNR is large ($\mathrm{SNR}\ge -$5 dB), almost all the estimators can reach the CRLB except for the MUSIC method when $\mathrm{SNR}=$ 20 dB. The accuracy is limited in very high SNR because the search grid in MUSIC is not fine enough. In addition, the MUSIC method and the ML-grid search suffer from a severe threshold effect when the SNR is less than $-5$ dB, and the corresponding RMSEs deviate from the CRLB dramatically. The ML-iterative is not worse due to the good initial guess from the prior knowledge. On the contrary, the proposed three Bayesian estimators have better performance in staying with the ECRLB when the SNR becomes lower. Therefore, the threshold effect can be delayed after using the prior knowledge.

Figure 5 illustrates the RMSE of all the estimators when increasing the number of snapshots from 5 to 50. The SNR is set to $-10$ dB. All the Bayesian estimators can reach the CRLB in the entire range that are irrespective of the snapshots number. However, the performance of the MUSIC and the ML-grid search deteriorate when the number of snapshots is limited. As the number increases, their performance becomes better and is approaching the CRLB performance. The results in Figure 4 and Figure 5 are consistent with the previous analysis that the performance of the traditional DOA estimation methods using the observation data will only be close to the performance of the Bayesian estimators when the SNR is high and the number of snapshots is large.

The effect of the prior knowledge (i.e., ${\theta }^o$ and ${\sigma }_{\theta }$) on the DOA estimation performance is shown in Figure 6 and Figure 7, where SNR $=0$ dB. Figure 6 shows the RMSE of the proposed estimators in different angles from $0^{°}$ to ${80}^{°}$. The performance in the negative angle part is similar. The CRLB increases slowly as the angle increases. The proposed three Bayesian estimators can reach the CRLB in the entire angle range, but the other three estimators that use the observation data only deviate from the CRLB when the angles are in the region around ${80}^{°}$. Figure 7 shows the RMSE of the proposed estimators when the prior standard deviation ${\sigma }_{\theta }$ varies from $0.5^{°}$ to $3^{°}$. The value of ${\sigma }_{\theta }$ has less of an effect on either the CRLB or the RMSE of any estimator. Figure 8 and Figure 9 show the results when reducing the SNR to $-10$ dB and repeating Figure 6 and Figure 7. They show that the MUSIC method and the two ML estimators will deviate from the CRLB, yet the proposed Bayesian estimators still reach the CRLB accuracy. The results validate again that the prior knowledge will affect the DOA estimation performance in the situations with limited observations or low SNR.

5.2. Multiple Targets DOA Estimation

This subsection intends to present the performance of DOA estimation for multiple targets. Thus, the grid search estimators (ML and MAP) are excluded from the estimators in the simulation due to the time-consuming multidimensional search. We consider the case of two adjacent narrowband signals from $\mathit{\theta }$ with prior knowledge ${\mathit{\theta }}^o={[0,30]}^T$ and ${\mathit{Q}}_{\theta }=\mathrm{diag}\left([4,1]\right)$ impinge onto the array simultaneously. Both signals have equal amplitude. Hence, the SNRs are also equal.

Figure 10 illustrates the histograms of the estimation errors obtained from 10,000 Monte Carlo simulations for two uncorrelated sources. The SNR is $-15$ dB. The errors of three Bayesian estimators are in a narrow region close to zero, but the errors of the ML-iterative solution and the MUSIC method are in a wider region. In addition, the error distribution of the MUSIC method is asymmetrical, which means that it is a biased estimator in this situation. The results validate again that the Bayesian estimators outperform the estimators without prior knowledge at low SNRs.

Figure 11 shows the RMSE of the proposed methods for two uncorrelated sources. The behaviors of the estimators are consistent with those in Figure 4. The results validate again that the proposed Bayesian DOA estimators outperforms the other methods, especially when the SNR is low. However, the MMSE estimator cannot reach the CRLB when the SNR is high ($\mathrm{SNR}\ge 10\mathrm{dB}$). The reason is that the number of the samples $N_m$ is insufficient.

Figure 12 shows the RMSE of the proposed methods for two correlated sources. It is well known that the MUSIC method would fail to produce peaks at the DOA locations of correlated sources if there is no pre-processing technique, such as forward–backward (FB) averaging and spatial smoothing, to handle the rank deficiency of the source covariance matrix. Because there is no pre-processing step, the performance of the MUSIC method in Figure 12 is very poor. However, the other estimators do not experience such problems when faced with correlated signals. Moreover, the proposed Bayesian estimators have the best performance in staying with the CRLB.

Figure 13 further demonstrates the performance of the proposed methods for two sources when varying their correlation coefficients from 0 to 1 where SNR $=-5$ dB. The iterative adaptive approach for amplitude and phase estimation (IAA-APES) in [41] can work well with few snapshots (even one), uncorrelated, partially correlated, and coherent sources. Hence, it is included in the simulation for comparison as well. It shows that only the MUSIC method suffers from performance degradation for highly correlated signals when the correlation coefficient is greater than $0.5$. For highly correlated signals, the source covariance matrix is close to singular and the performance is prone to be affected by noise. In contrast, all the other estimators and the CRLBs are not affected by the signal correlation.

6. Conclusions

This paper investigates the Bayesian DOA estimation problem by combining multiple snapshots from observation and prior knowledge from the target tracker. We started with the signal model for ULA and expressed the prior knowledge by Gaussian distribution. We then derived the BCRLB and ECRLB for Bayesian DOA estimation. Analysis shows that the prior knowledge will improve the estimation performance in the situations with a limited number of observations or at low SNR. Based on the MAP estimator, a grid search method was proposed for single DOA estimation. Considering the computational burden of the MAP-grid search increases exponentially with the number of targets. Hence for a general multiple-targets case, a MMSE estimator based on Monte Carlo method and an iterative solution based on the MAP estimator were proposed for Bayesian DOA estimation. The former needs to trade off accuracy against computational complexity. The latter requires that the prior knowledge is informative. Simulations show that all three proposed Bayesian estimators aided by the prior knowledge are better than the methods (the ML estimators and the MUSIC method) only using the observation data.

The proposed Bayesian DOA estimation methods can be united with the tracking approaches, such as the Kalman filter or extended Kalman filter, by assuming dynamic equations for the targets. In fact, the existing Kalman filters are implemented in the measurement data level (such as angle and range), but the Bayesian DOA estimation is implemented in the signal level (array snapshots).

In this paper, the proposed DOA estimation is carried out in the conventional Bayesian framework, i.e., the MAP method and MMSE method. In general, the computationally complexity is high due to grid search or Monte Carlo computation. The existing sparse DOA estimation methods show good performance in harsh scenarios such as insufficient array observations and coherent imping signals. Future work can extend these sparse methods to Bayesian DOA estimation by incorporating the prior knowledge. In particular, the weighted atomic norm approach in the gridless sparse methods can take advantage of the prior knowledge to enhance the estimation performance [26,42,43]. We are currently conducting studies of the Bayesian DOA estimation problem by using sparse estimation methods.