Sparse Bayesian learning for wideband direction-of-arrival estimation via beamformer power outputs in a strong interference environment

Sparsity-based direction-of-arrival (DOA) estimation has attracted considerable attention in underwater source localization because it can be applied under low signal-to-noise ratio (SNR) conditions compared with traditional high-resolution methods.¹ This kind of method represents the DOA estimation as an underdetermined linear problem because only a few sources exist in ocean acoustics. The method can be roughly classified into two categories: the l_p-norm minimization method^1–3 and the sparse Bayesian learning (SBL) method.^4,5 The l_p-norm minimization method uses the l_p-norm to enforce the sparsity on the solution and solves a regularized optimization problem with an appropriate regularization parameter to estimate the DOAs. However, it is difficult to tune the regularization parameter satisfactorily. In contrast, the SBL is a probabilistic parameter estimation approach that assigns the suitable prior for the signal and enforces the sparsity by automatically estimating the model parameters from the received data, thus avoiding tuning the regularization parameter.

Wideband processing plays a fundamental role in underwater source localization. Das et al. extended the SBL into the wideband condition and proposed a wideband SBL (WSBL) method.⁶ The WSBL makes full use of the common sparsity profile across the frequency bins by assigning the Gaussian priors to the signals in all frequency bins with the same precision. The DOAs are then estimated from raw array outputs in a Bayesian framework. The idea is also used in Ref. 7. The main difference between Refs. 6 and 7 is that the DOAs in Ref. 7 are estimated from the covariance matrix; thus, a higher array SNR than raw array outputs is obtained.⁸ 

However, the performances of the sparsity-based methods are undermined by the existence of the nearby interferences with strong power when the power of the target is weak.^9,10 To solve DOA estimation problem in an interference environment, Ren et al.¹¹ proposed an eigenanalysis-based adaptive interference suppression (EAAIS) method. They defined a power ratio to identify the eigenvectors dominated by the interferences adaptively. The interference subspace was then removed from the covariance matrix, and conventional beamforming (CBF)¹² was used to estimate the DOAs. This method requires only a rough bearing range of the target but is affected by low resolution.

Different from the EAAIS-based CBF (EAAIS-CBF),¹¹ Yang et al.¹⁰ used the matrix filter with nulling (MFN)¹³ as a preprocessor to suppress the interferences. The MFN is a spatial filter, which facilitates the passage of the signals while suppressing the interferences. To achieve high resolution, Yang et al. extended the sparse spectrum fitting (SpSF)² to estimate the DOAs from the MFN outputs and introduced a method named MFN-based SpSF (SpSF-MFN).¹⁰ However, this method is afflicted by low computational efficiency because much time is needed to compute the MFN. The beamformer^12,14 is another kind of spatial filter. References 15 and 16 adopted the CBF¹² as a preprocessor. Operation in the beam domain offers several advantages over that in the element domain, including reduced computational complexity and low SNR resolution thresholds.¹⁷ In these methods, the beamformer power outputs (BPO) were expressed as a linear relation of the signal power vector, and the SpSF was used to estimate the DOAs from this relation. However, the CBF has limitations in interference suppression because of high sidelobes. If the interferences are not sufficiently suppressed, the residual interference after beamforming will affect the following DOA estimation. Moreover, these methods need to tune the regularization parameters.

This letter fuses the beamformer and the SBL, allowing the SBL method to be used in a strong interference environment. The minimum variance distortionless response with diagonal loading (MVDR-DL)¹⁴ is adopted because it produces deep nulls toward the directions of the interferences to sufficiently suppress them. The main contribution of this study lies in the establishment of the likelihood for the relation of the BPO based on the perturbation model for this relation. Inspired by Ref. 6, the Gaussian priors with the same precision are assigned to the signal power vectors in all frequency bins. Thus, a BPO-based SBL (SBL-BPO) method is provided by adopting the variational Bayesian inference (VBI)¹⁸ to estimate the DOAs via the BPO. The simulation and experimental results show that the SBL-BPO achieves high estimation precision and computational efficiency in a strong interference environment.

In this paper, $(·)T$⁠, $(·)H$⁠, and $(·)*$ denote the transpose, the conjugate transpose, and the conjugate, respectively. $vec{·}$ represents the vectorization operator. $Re[·]$ and $Im[·]$ are the real and the imaginary parts of a complex variable, respectively. $diag(X)$ and $diag(x)$ are the vector with diagonal elements of X being its elements and the diagonal matrix with x being its diagonal elements, respectively. $IM$ is an $M×M$ identity matrix. $⊗$⁠, $∘$⁠, and $⊙$ denote the Kronecker, the Hadamard, and the Khatri-Rao products, respectively.

Consider that an array with M elements receives the signals from K_S targets and K_D interferences at $θ⌢S=[θ⌢1,…,θ⌢KS]$ and $θ⌣D=[θ⌣1,…,θ⌣KD]$⁠, respectively. The array-received data are divided into N blocks, and an L-point discrete Fourier transformation (DFT) is used for each block. Thus, the array output $xl(n)∈ℂM×1$ at the lth frequency bin is modeled as follows:

where $xl(n)∈ℂM×1$⁠, $sl(n)∈ℂKS×1$⁠, $dl(n)∈ℂKD×1$⁠, and $el(n)∈ℂM×1$ represent the DFT coefficients of the received data, the desired signals, the interfering signals, and the additive noise at the nth block, respectively. $Al(θ⌢S)=[al(θ⌢1),…,al(θ⌢KS)]$ and $Al(θ⌣D)=[al(θ⌣1),…,al(θ⌣KD)]$⁠, where $al(θ)$ is the steering vector pointing the direction $θ$⁠.

Assume that the desired and the interfering signals are uncorrelated. Furthermore, assume that the noise is uniform white Gaussian noise of the power $σl$ and uncorrelated with the signals. The covariance matrix $Rl∈ℂM×M$ is expressed as

where $PlS=diag(plS)$ and $PlD=diag(plD)$⁠, in which $plS$ and $plD$ contain the powers of the desired and the interfering signals, respectively. The sample covariance matrix $R̂l=∑n=1Nxl(n)xlH(n)/N$ is always used to approximate the true covariance matrix. The relation between $R̂l$ and $Rl$ is formulated as

In this equation, $El$ is the error matrix that obeys the complex Gaussian distribution¹⁹ 

where $Ql=RlT⊗Rl/N$⁠.

Assume that the targets lie in a sector $ΘS=[ΘSL,ΘSR]$⁠, where $ΘSL$ and $ΘSR$ are the left and the right limitations of $ΘS$⁠, respectively. The beamforming matrix in $ΘS$ is recorded as

where $wl(ϕ)$ represents the weight vector of a beamformer with a pointing angle $ϕ$⁠, and $ϕ=[ϕ1,…,ϕKB]$ contains the pointing angles, where $ΘSL≤ϕ1<ϕ2<…<ϕKB≤ΘSR$⁠. Thus, the BPO is given by

With the substitution of Eq. (3) into Eq. (6), the following can be obtained:

where $Ãl(θ⌢S)=(WlHAl(θ⌢S))*∘(WlHAl(θ⌢S))$⁠, $Ãl(θ⌣D)=(WlHAl(θ⌣D))*∘(WlHAl(θ⌣D))$⁠, $ĩlB=diag(WlHWl)$⁠, and $ε̃lB=diag(WlHElWl)$ is the perturbation vector for the BPO model.

The MVDR-DL beamformer¹⁴ is used in this letter because it produces deep nulls toward the directions of the interferences, and the corresponding beam responses decrease as the interference powers increase to suppress them sufficiently. The weight vector is¹⁴ 

where $ϕ$ is a pointing angle and $ηl$ is the DL level. The value is chosen as $ηl=10σ̂l$⁠,¹⁴ where $σ̂l$ is the noise power estimate. It is worth mentioning that other adaptive beamformer options are also possible. After beamforming, the interference powers are largely reduced in the BPO and can be omitted. Hence, Eq. (7) is approximated as

With the sector $ΘS$ divided into the discretized grid $φ=[φ1,…,φKG]$⁠, Eq. (9) is rewritten as

where $Ãl=(WlHAl(φ))*∘(WlHAl(φ))$ and $pl$ is a zero-padded version of $plS$⁠. For any $n=1,…,KS$⁠, $(pl)m=(plS)n$ holds if $φm=θ⌢n$⁠, where $(pl)m$ and $(plS)n$ are the mth and the nth entries of $pl$ and $plS$⁠, respectively. Otherwise, $(pl)m=0$⁠.

Likelihood: It is difficult to directly obtain the distribution of $ε̃lB$ from the relation $ε̃lB=diag(WlHElWl)$⁠. Hence, the following transformation is performed:

where $J∈ℝKB×KB2$ is a matrix whose the (m, n)th element $[J]mn=1$ if $n=KB(m−1)+m$⁠; otherwise, $[J]mn=0$⁠. In this manner, the distribution of $ε̃lB$ can be easily obtained based on Eq. (4), which is given by

where $Q̃lB=[J(WlT⊗WlH)]Ql[J(WlT⊗WlH)]H=(WlHRlWl)*∘(WlHRlWl)/N$⁠. Thus, the likelihood function for $r̃lB$ is

Equation (13) can be converted to the following real-valued Gaussian distribution:²⁰ 

It is unnecessary to consider the distribution of $Im[r̃lB]$ because $r̃lB$ is defined in the real domain. Hence, only the distribution of $Re[r̃lB]$ is retained, that is,

where $Q¯lB=Q̃lB/2$⁠. $Re[·]$ is omitted for simplification.

Prior: A real Gaussian prior is assigned to $pl,l=1,…,L$⁠,

where $Γ=diag(γ)$⁠, and $γ=[γ1,…,γKG]T$ is the precision vector that controls the sparsity of $pl,l=1,…,L$⁠. Moreover, a real Gaussian prior is employed on $σl$⁠, as given by

where $γlN$ is the variance of $σl$⁠.

Based on Eqs. (15)–(17), the posterior distribution is expressed as follows:

where $Πl=1Lp(pl,σl,r̃lB;γ,γlN)=Πl=1Lp(r̃lB|pl,σl)p(pl;γ)p(σl;γlN)$ is the joint distribution of all unknown and observed quantities. The $pl$ and $σl$ are treated as the latent variables, whereas $γ$ and $γlN$ are the parameters.

The VBI¹⁸ is used to approximate Eq. (18). The approximated posterior distribution of the latent variables is expressed as

where $q(·)$ is the approximated posterior, and $Ω=[p1T,…,pLT,σ1,…,σL]T$⁠. The approximated posterior with respect to the individual variable is updated by¹⁸ 

where $Ωi$ is the ith entry of $Ω$⁠, $⟨·⟩q(·)$ represents the expectation with respect to $q(·)$⁠, and $q(Ω\Ωi)$ denotes the approximated posterior except $Ωi$⁠. The entries in $Ω$ are updated alternately while fixing the rest. The mean of $Ωi$⁠, recorded as $⟨Ωi⟩$⁠, provides the estimate of $Ωi$⁠. No information priors are assigned to $γ$ and $γlN$⁠. They are updated by maximizing the expected log-likelihood function.

According to Eq. (20), the approximated posterior with respect to $pl$ is obtained as

With the substitution of Eqs. (15) and (16) into Eq. (21), the approximated posterior distribution of $pl$ is as follows:

In this equation,

where $Γ̂$ is the estimate of $Γ$⁠. Similarly, $q̂(σl)$ also follows a real Gaussian distribution, whose covariance and mean are

respectively, where $γ̂lN$ is the estimate of $γlN$⁠.

The parameter $γ$ is updated by maximizing the expected log-likelihood function, i.e.,

After differentiating Eq. (25) with respect to $γp$ and setting it to zero, the estimate of $γp$ is expressed as

Furthermore, the parameter $γlN$ is updated similarly as given by

A summary of the SBL-BPO is provided in Table 1. The initial values of the iteration are set to $γ̂k(0)=L/∑l=1L|wlH(φk)R̂lwl(φk)|2, k=1,…,KG$⁠, $⟨σl⟩(0)=σ̂l$⁠, and $(γ̂lN)(0)=σ̂l2$⁠. The iteration is terminated when $‖γ̂(i)−γ̂(i−1)‖2/‖γ̂(i−1)‖2≤τ$ for a small tolerance $τ$ or a maximum number of iterations $Itermax$ is reached, where the superscript (i) denotes the ith iteration. Once the iteration converges, the estimated spectrum $PSBL−BPO$ is obtained as

For the SBL-BPO, the updating workloads for the signal and the noise parameters in each iteration are $O(LKBKG2)$ and $O(LKBKG)$⁠, respectively. The workload is then compared with that of the WSBL.⁶ Assume that the number of the grid points used in the WSBL is $K̃G$⁠. The workload is determined by signal parameter updating with $O(LMK̃G2)$ operation in each iteration. The SBL-BPO can achieve a lower computational workload than the WSBL mainly due to the following reasons. First, the SBL-BPO is defined in the beam domain whose computational workload decreases by setting a small $KB$⁠. $KB$ should be smaller than M to ensure that $Q¯lB$ is nonsingular, i.e., $KB<M$⁠. Second, the SBL-BPO can search for the DOAs only in the sector $ΘS$ to decrease the computational workload because the MVDR-DL is used to suppress the signals outside of $ΘS$⁠. For the WSBL, the bases of the out-of-sector signals should be included in the dictionary matrix. Otherwise, the powers of the out-of-sector signals will leak to $ΘS$⁠, affecting the DOA estimation in $ΘS$⁠. Hence, the searching area in the WSBL is larger than that in the SBL-BPO, i.e., $K̃G>KG$⁠.

The SBL-BPO is compared with the sparse reconstruction based on multiple beamspace measurements (MBM-SR),¹⁵ the SpSF-MFN,¹⁰ the EAAIS-CBF,¹¹ and the WSBL.⁶ A uniform linear array (ULA) with 32 sensors spaced at 4 m is considered. The directions of two targets and one interference are –10°, –7°, and 10°. The frequency range is [90, 180] Hz. The received data sampled at 2 kHz are divided into N blocks with the lengths of 1000. A 1000-point DFT is applied in each block. The sector $ΘS$ in the SBL-BPO and the MBM-SR is set to [ –16°, –2°] and divided with a step of 2° to obtain the beam pointing angles, i.e., $KB=8$⁠. The passband and the stopband in the MFN are set to [−16°, −2°] and [−90°, −23°]∪[5°, 90°], respectively. The stopband attenuation is set to −15 dB. The rough bearing range in the EAAIS-CBF is set to [−16°, −2°]. The SBL-BPO, the MBM-SR, and the SpSF-MFN search for the DOAs in $ΘS$ with a grid interval of 1°, i.e., $KG=15$⁠, for a fair comparison. The WSBL estimates the DOAs in the entire space with a grid interval of 1° to ensure the inclusion of the bases of the signals in the dictionary matrix. The tolerance $τ$ and the maximum number of iterations $Itermax$ are set to $10−3$ and 3000, respectively. All simulations are performed in matlab on a PC with an Intel Core i7–6820HQ CPU and 32 GB RAM.

The methods are first compared from the perspective of estimation precision, examined by the root mean square error (RMSE),

where $θ̂kw$ and $θk$ are the estimated DOA in the wth simulation and the true DOA, respectively. W is the sum of the Monte Carlo runs, which is set to 200.

Figure 1 shows the RMSE versus the signal-to-interference ratio (SIR) by fixing the SNR and the number of snapshots N to −10 dB and 30, respectively. The interference power increases with a decreased SIR. The performance of the SBL-BPO is hardly affected as the interference power increases due to the sufficient suppression of the interference by the MVDR-DL, decreasing its influence on the DOA estimation. By contrast, the performance of the WSBL degrades when $SIR<−20 dB$⁠. The performance of the SpSF-MFN is also unaffected by the interference because the MFN produces a deep null to suppress the interference sufficiently. For the MBM-SR, reconstructing the residual interference after beamforming on the corresponding basis is difficult due to the small entries of the basis outside of $ΘS$⁠, leading to the leakage of its energy to $ΘS$⁠. The DOA estimation is affected by this leakage seriously if the interference is not suppressed sufficiently. Hence, the performance of the MBM-SR degrades rapidly when $SIR<−15 dB$ due to the limitations of the CBF in interference suppression. The EAAIS-CBF cannot resolve two targets because of low resolution.

Figure 2 shows the RMSE versus the SNR by fixing the SIR and N to −30 dB and 30, respectively. The SBL-BPO and the SpSF-MFN employ the spatial filters to improve the signal-to-interference-and-noise ratio. Hence, their RMSEs are lower than those of other methods when $SNR<−2.5 dB$⁠. The performance of the WSBL is worse than that of the SBL-BPO due to the existence of the interference. The RMSE of the MBM-SR remains nearly unchanged because the estimation result is dominated by the leakage of the residual interference. Furthermore, the EAAIS-CBF fails due to the low resolution of the CBF.

Finally, the computational efficiency is compared considering the running time. The simulation settings are the same as those in Fig. 1, except that the SIR is fixed to −10 dB. The mean running time over 200 trials is shown in Table 2. The time usage of the proposed method is substantially smaller than that of the WSBL. The EAAIS-CBF achieves high computational efficiency because its workload mainly depends on eigen-decomposition, but its performance is affected by low resolution. The MBM-SR and the SpSF-MFN use the SpSF² to estimate the DOAs, which need to solve convex optimization problems in each frequency bin. Furthermore, the SpSF-MFN needs to address the other convex optimization problem to design the MFN. Hence, they are afflicted by lower computational efficiency than other competitors.

Experimental data, which have also been used in Ref. 10, were processed in this subsection. The experiment was conducted in the South China Sea with a maximum depth of approximately 150 m in September 2016. A ULA, with 32 sensors spaced at 4 m, was towed behind an experimental vessel at the depth of 50 m under the sea surface. According to an automatic identification system, a few vessels located at approximately −18°, −25°, and −29° (Ref. 10) were regarded as targets and recorded as Target1, Target2, and Target3. The radiation noises from the experimental vessel and other surface sources were regarded as interferences. The received data sampled at 2048 Hz are divided into 146 frames with 10 s of data, and the overlap between the adjacent frames is 80%. The data in each frame are divided into 39 blocks with 50% overlap, i.e., N = 39. A 1024-point DFT is applied in each block. The analyzed frequency ranges from 90 Hz to 180 Hz. The $ΘS$ in the SBL-BPO and the MBM-SR is set to [−30°, −12°]. The rough bearing range in the EAAIS-CBF is set to [−30°, −12°]. The passband and the stopband in the MFN are set to [−30°, −12°] and [−90°, −35°]∪[−7°, 90°], respectively. The SBL-BPO, the MBM-SR, the SpSF-MFN, and the EAAIS-CBF estimate the DOAs in $ΘS$ with a step of 1°. The WSBL estimates the DOAs in the entire space to ensure that all bases of the signals are included. Other parameters are set the same as those in the simulations. All results are calculated in matlab on a PC with an Intel Core i7–6820HQ CPU and 32 GB RAM.

Figure 3 shows the result of each method. The SBL-BPO and the SpSF-MFN can resolve three targets because the influence of the interferences is alleviated by the spatial filters. The performance of the MBM-SR is affected by the residual interferences, and the directions of Target2 and Target3 cannot be estimated in Fig. 3(b). The EAAIS-CBF is afflicted by low resolution, thus failing to resolve Target2 and Target3. The WSBL also has some difficulties in resolving Target2 and Target3, particularly from 80 s to 120 s. It is worth mentioning that the array deformation always exists in practical signal processing. The robustness of the SBL-BPO to the array deformation mainly depends on the robustness of the beamformer. The MVDR-DL is robust to the sensor position errors,¹⁴ which increases the robustness of the SBL-BPO to such array deformation.

Table 3 lists the mean running time of each method over frames. Similar to Table 2, the SBL-BPO requires less time than other sparsity-based methods because of the low computational workload. The EAAIS-CBF is computationally efficient but leaves a low resolution. Moreover, the MBM-SR and the SpSF-MFN require more time than other competitors because they need to solve convex optimization problems in each frequency bin. Therefore, the experimental results confirm that the SBL-BPO achieves high resolution and computational efficiency in a strong interference environment.

This letter fuses the beamformer and the SBL method to increase the robustness of the SBL in a strong interference environment. The MVDR-DL is used to suppress the interferences, such that the interference power is much reduced in the BPO. The perturbation for the relation of the BPO is deduced in this letter, thus establishing the likelihood for this relation. Based on this, the SBL-BPO is put forward by using the VBI to the established Bayesian model, avoiding tuning the regularization parameter. The simulation and experimental results confirm that the proposed method improves the performance of the SBL method in a strong interference environment. Furthermore, the proposed method also achieves higher computational efficiency than other sparsity-based methods.

This work was supported by the National Natural Science Foundation of China (Grants Nos. 11974286 and 11904274).