Robust Sparse Bayesian Learning Scheme for DOA Estimation with Non-Circular Sources

: In this paper, a robust DOA estimation scheme based on sparse Bayesian learning (SBL) for non-circular signals in impulse noise and mutual coupling (MC) is proposed. Firstly, the Toeplitz property of the MC matrix is used to eliminate the effect of array MC, and the array aperture is extended by using the properties of the non-circular signal. To eliminate the effect of impulse noise, the outlier part of the impulse noise is reconstructed together with the original signal in the signal matrix, and the DOA coarse estimation is obtained by balancing the accuracy and efﬁciency of parameter estimation using the alternating SBL update algorithm. Finally, a one-dimensional search is used in the vicinity of the searched spectral peaks to achieve a high-precision DOA estimation. The effectiveness and robustness of the algorithm for dealing with the above errors are demonstrated by extensive simulations.


Introduction
In recent years, with the development of array signal processing, direction-of-arrival (DOA) estimation has gradually become a research hotspot [1], and its development has a significant impact on the fields of radar, sonar, and mobile communications [1,2]. With the innovation and development of technology, DOA estimation techniques are widely used in target localization and navigation systems [3][4][5]. The DOA estimation technique originated with Pisarenko's discovery that DOA information could be obtained from the second order statistics of the receiver matrix. Soon after, many traditional direction finding algorithms were proposed. Subspace class DOA estimation algorithms represented by multiple signal classification (MUSIC) based methods [6,7] and rotational invariant techniques (ESPRT) based methods [8,9] were then proposed, which mainly exploit the relationship between the signal subspace and the noise subspace to achieve DOA estimation. Under ideal conditions, subspace class methods can perfectly approximate the Cramer-Rao lower bound (CRLB) [7]. However, the development of science and technology has made the electromagnetic wave environment more and more complicated. Non-ideal conditions, such as the correlation between signal sources, low signal-to-noise ratio, and small snapshots, make subspace class algorithms have great limitations in practical applications. Then, subspace fitting algorithms represented by weighted subspace fitting (WSF) [10,11] and maximum likelihood (ML) algorithm [12,13] are proposed. Although these types of algorithms solve the problem of poor performance of subspace algorithms under non-ideal conditions to a certain extent, it is restricted by computational efficiency, and it is difficult to exert its performance advantages.
At a time when the development of traditional high-resolution DOA estimation algorithms is limited by environmental or efficiency constraints, the rise of compressive sensing (CS) techniques offers a new solution to target vectoring. CS technology that has emerged in recent years has improved the large amount of redundancy and waste caused by traditional sampling rules (i.e., the Nyquist sampling theorem) [14]. It directly samples useful information, which greatly saves the cost of data sampling, and can recover the original signal from a small number of measured values through the numerical optimization algorithm. From then on, a large number of DOA estimation methods based on sparse representation (SR) have been proposed, such as sparse Bayesian learning (SBL) [15,16] and convex relaxation class methods [17][18][19]. These methods can still effectively and accurately estimate DOAs in non-ideal situations with small snapshots, low SNR, and high source correlation. Although convex optimization class algorithms indirectly solve the NP-hard problem of l 0 -parametric optimization, they suffer from high computational complexity and difficulty in ensuring global convergence, which is well solved by SBL-based algorithms [15,16]. However, the advantages of SBL-based vectoring algorithms often rely on ideal array manifolds and Gaussian noise environments.
In the actual array signal processing, compared with a separate antenna, the antenna array is not only affected by its own magnetic field, but also by the electromagnetic field of adjacent elements. If the elements of array are not fully calibrated, it will directly affect the array manifold. Nevertheless, the accuracy of DOA estimation depends on the perfect array manifold, and the existence of mutual coupling (MC) effect will make most highresolution DOA estimation methods lose estimation accuracy to a certain extent. To solve this problem, a number of relevant solutions have been proposed [20,21], and an iterative mutual coupling calibration method is proposed in [20], which is able to complete the array calibration under the condition that the prior information of the flow pattern matrix is unknown, but it has high computational burden associated with iteration. Ye et al. proposed a non-iterative algorithm in literature [21], which can achieve the Toeplitz structure of the MC matrix only without any calibration source. This approach was subsequently extended and applied to various algorithms for DOA estimation by Dai et al. [22,23]. On the other hand, non-Gaussian noise with significant spike characteristics, called impulse noise, exists in the real environment. Under conventional Gaussian noise models, these large spikes can be considered as outliers. To obtain ideal DOA information under impulsive noise, many SR-based methods model impulsive noise as a combination of Gaussian noise and outlier matrices [24][25][26]. Dai et al. then rewrite the signal matrix with outlier matrices to avoid the effect on DOA estimation in SBL iterations [24]. Zheng et al. refer to the signal model in [24] and propose a variational SBL iterative estimation algorithm. This algorithm identifies and corrects outliers during the iterative process [25]. This iterative estimation algorithm thus becomes an efficient means of dealing with impulse noise.
Whether it is the subspace based methods or SR based methods, the performance of DOA estimation mostly depends on the array aperture. Using the properties of the signal itself to extend the array aperture is an effective way to improve the performance of DOA estimation. It has been shown that the non-zero elliptic covariance property of non-circular signals can not only distinguish non-circular signals from circular signals in the source, but also bring better performance to the traditional DOA estimation methods. Due to this property, non-circular signals represented by BPSK and ASK are widely used in digital signal processing [27][28][29][30][31]. For example, Zheng et al. used the property of non-zero elliptic covariance of non-circular signals to construct partial joint sparsity of non-circular sources in the covariance vector and pseudo-covariance vector [31]. The non-circular sources are then identified by setting a set of latent variables, and finally DOA estimation under the condition of co-existence of non-circular and circular sources is achieved. On the other hand, array expansion is accomplished by exploiting the property of non-zero elliptic covariance of non-circular signals, which can effectively expand the array aperture. For example, Zhang et al. exploited the properties of non-circular signals to process the received data and extended the array aperture, which improved the efficiency of the algorithm to some extent [30].
However, the impulsive noise environment was considered in the literature [30] without the array MC. This paper proposes a robust algorithm based on SBL for noncircular signals in impulsive noise and MC. The proposed method uses the selection matrix to eliminate the influence of the mutual coupling of the array, and uses the characteristics of the non-circular signal to expand the array and expand the array aperture. Then, the outlier matrix of impulse noise is put into the signal model, and the alternating update SBL algorithm is used to estimate the parameters of signal and noise, respectively. This procedure solves the problem of too high complexity of using traditional SBL to deal with impulse noise, and better balances the accuracy and efficiency of DOA estimation. Finally, after the rough estimate of DOA is obtained, an improved one-dimensional search method is used according to the reconstruction result to obtain a high-precision DOA estimate. A large number of simulation experiments show that, under the influence of impulse noise and array MC, the proposed algorithm can mitigate the impact of off-grid errors to a certain extent and has good DOA estimation performance.
Notations: Lowercase bold letters and uppercase bold letters, respectively, denote vectors and matrices. Inverse, conjugate, transpose, and conjugate transpose operations on a matrix are denoted by (·) −1 , (·) * , (·) T , and (·) H , respectively. tr(·) denotes the trace operation on a matrix. The diagonalization and vectorization operations are represented by diag{·} and vec(·), respectively. E A B is to find the expectation of A with respect to B. R{·} is the real part operator.

Data Model
We consider that, on a uniform linear array (ULA) composed of N identical sensors, K uncorrelated narrowband far-field signals are incident, where K ≤ N − 1. In addition, the distance between adjacent elements in a ULA is d = ν/2, where ν represents the wavelength of the incident signal. Therefore, for strictly non-circular signals, under conditions of array mutual coupling and coexistence of impulse noise, the data receiving matrix is rewritten as and T is the number of snapshots. In the data model shown in Equation (1), the presence of MC allows the array manifold under ideal conditions, which is represented as A = [a(θ 1 ), a(θ 2 ), . . . , a(θ K )] T to change intoÃ , where a(θ k ) = [1, e −jπ sin θ k , . . . , e −jπ(N−1) sin θ k ] is the steering vector for the k-th signal source, and band-symmetric Toeplitz matrice C is often used to describe MC effects in ULA, and it is indicated as where c is the MC coefficient. The non-circular signalS = [s(t 1 ),s(t 2 ), . . . ,s(t T )] is split and modeled as a noncircular phase matrix Ψ = diag[e jψ 1 /2 , e jψ 2 /2 , . . . , e jψ K /2 ] and a real-valued signal S, where ψ k is a non-circular phase and, when it is zero, Ψ is an identity matrix andS is a strictly non-circular signal satisfyingS * =S.
Similarly, the spike characteristics of impulse noise will cause errors in DOA estimation. To avoid the impact of its non-Gaussian part on the DOA estimation performance, the noise matrix is usually modeled as the sum of the outlier matrix H = [h(t 1 ), h(t 2 ), . . . , h(t T )] and the Gaussian noise matrix N = [n(t 1 ), n(t 2 ), . . . , n(t T )] according to the impulsive characteristics.

The Proposed Method for NC Sources
As the proposed method requires the use of sparse recovery techniques to recover the signal, constructing a sparse-complete dictionary relies on a perfect array manifold. We use the banded complex symmetric Toeplitz structure of the MC matrix to construct the selection matrix J to remove the MC: where P is the number of MC coefficients, speaking of the selection matrix acting on the steering vector. The principle of eliminating mutual coupling is clearly expressed as where τ k = 1 + ∑ K k=1 2c k cos(kπsinθ k ). When the left and right ends of Equation (1) are respectively multiplied by the selection matrix J, the data receiving matrix after eliminating MC can be obtained whereǍ is the ideal array manifold, and it is composed of the middle N 1 line ofÃ with N 1 = N − 2P, which is illustrated in Figure 1.
when the angle is determined, this parameter τ k is a constant.
Uniform Linear Array (ULA) As shown in the array model in Figure 1, the decoupling process results in a loss of array aperture. However, the properties of non-circular signals can be exploited to extend the data reception matrix by constructing a virtual array. This allows for DOA estimation under more targets. At the same time, all possible source angles θˆk = [θ 1 , θ 2 , . . . , θK] are assumed to exist in the discrete space [−π/2, π/2], whereK is the number of discrete grids and satisfiesK N > K. The extended data reception matrix space domain discretization is denoted as where D is the permutation matrix with the anti-diagonal elements being 1, and the remain- The signal reconstruction process based on SBL techniques uses the observation matrix and prior probabilities to recover the signal matrix. The corresponding array manifolds are further obtained, and the final DOA estimates are achieved. Since the signal matrix is the object to be estimated in this technique, the method of reconstructing the outlier matrix into the signal matrix directly eliminates the detrimental effects on signal recovery. Therefore,S is the signal after re-modeling with the outlier matrix in Equation (6) whereS = Ṡ * Ṡ T , andṠ is the signal matrix after discretization of the null domain, with each of its non-zero rows corresponding to a real signal. Similarly,H = DȞ * Ȟ T ,N = DŇ * Ň T .

A Priori Setting of Parameters
According to the received data model described in Equation (6), we can obtain the prior probability of Z as where β = σ 2 is the power of Gaussian noiseN, which is usually unknown, so it will be assumed to satisfy the gamma distribution: where the hyperparameter ι is usually set to 1 and → 0 to obtain a broad hyperprior. Similarly, the signalS is assumed to satisfy a Gaussian distribution with variance λ = [λ 1 , λ 2 , . . . , λ 2K ], where λ i corresponds to the signal variance of the i-th row. Its prior probability is given by where Λ = diag(λ), and the hyperparameter λ i in the formula is further modeled as a gamma distribution where is usually set to = 0.01. The irregular characteristics of the outlier matrix result in sparse elements, so each element ofH needs to be assigned a Gaussian prior distribution with different variances: where ω n,l is used to represent the variance of the element in the n-th row and l-th column ofH, and [H] n,l =h n,l , [Ω] n,l = ω n,l Since Ω is unknown, it is modeled as a gamma distribution Gamma(ω n,l |1, ).
If the prior ofH andS were constructed in the form ofS, according to the Equation (6), it is expressed as where Γ l = diag([λ; ω l ]) and ω l represents the l-th column of Ω. According to the derivation process of SBL, the power Γ is updated in the last iteration. However, the sparsity of the elements of barbmH determines that its variance is three-dimensional, and each iteration is accompanied by two matrix inversion processes, which greatly affects the rate of the algorithm. In contrast to the joint estimation method proposed in Equation (20), this paper makes use of the alternating SBL update algorithm proposed in [24], which aims to repairS while updating the matrixH. This method reduces the dimensionality of the inverse matrix and thus the computing time to some extent. In this case, theS andH components ofS are estimated separately, and the changed data models are given separately among them, Q is an orthogonal matrix whose characteristics are used to reduce the complexity of subsequent matrix operations.

Estimate SignalS
According to the a priori hypothesis process of Equation (8) to Equation (10), and based on an Equation (14) model, the prior probability of Z when only the signal is considered is given as The classical SBL-based DOA estimation algorithm [16] uses Bayesian formula to estimate the parameters of unknown signal matrixS based on known observation data Z and reasonable priori assumptions. However, the estimation ofS is related to the estimation of the parameters {λ, β}, which can be achieved by finding the corresponding parameter that makes P(λ, β|Z) maximal. This process can be simply expressed as The reason for the above transformation is that p(Z) is constant when the observation matrix Z is known, and the equation P(λ, β|Z)p(Z) = P(Z, λ, β) can be given by the conditional probability formula. Then, the posterior probability ofS can be obtained by using the Bayesian formula, i.e., where µ l = βΣA H z l , l = 1, 2, . . . , L, where C = (β −1 I 2N 1 + AΛA H ) [32]. According to the expectation maximisation (EM) algorithm, the following likelihood function is maximized to obtain the estimation of the corresponding parameters, which is For the updating of parameter λ, the term associated with β is ignored. Then, Equation (21) can be rewritten as E ln(p(S|λ)p(λ)) p((S|Z;λ,β)) = E ln(p(S|λ)) p((S|Z;λ,β)) + log(p(λ)) where ∆ l = µ l 2 2 + Σ. Similarly, when estimating the parameter β, the other irrelevant terms are ignored. Then, Equation (21) can be simplified to E ln(p(Z|S, β)p(β)) p((S|Z;λ,β)) = E ln(p(Z|S, β)) p((S|Z;λ,β)) + ln(p(β)) where Ξ = ∑ L l=1 z l − As l 2 F /L + tr(AΣA H ). Then, by taking the derivatives of Equations (22) and (23) with respect to λ and β, respectively, and making them zero, the updating rules for λ and β can be obtained as follows: After all, the modified signalS can be obtained as

Alternating Estimate Outlier MatrixH
WhenH is regarded as a signal, the prior probability of receiving data matrix is rewritten as SinceS in Equation (14) is replaced by the updatedS in equation Equation (26), we find that the noise variance has changed, which is denoted asβ here.
Similar to the procedure from Equation (21) to Equation (23), the posterior probability ofH with respect to data Z is deduced as where µ l =βΣ l Q H z l , l = 1, 2, . . . , L, Then, it can be easily obtained that the updating rule for parameters iŝ whereΞ = ∑ L l=1 z l − Qµ l 2 2 , and ∇ = (β new + (ω new n,l ) −1 ) −1 . Then, the updatedH can be expressed asH [µ 1 , µ 2 , . . . , µ L ]. (33) When the parameter exceeds the set threshold condition λ t − λ t−1 2 / λ t−1 2 < 10 −4 or exceeds the set maximum number of iterations T = 1500, the update is stopped. Then, the rough estimation of DOA is obtained by searching for the first K peaks ofS.
As can be seen from the spatial domain discretization model shown in Figure 1, the real signal source may not fall exactly on the divided grid points. The gap between true DOA and grid point is called the off-grid error. To achieve high-resolution DOA estimation, the influence of off-grid errors should be minimized. Thus, a fine one-dimensional search procedure around K spectral peaks is used to obtain more accurate DOA estimation [33], which is specifically expressed as whereh(θ) =ā(θ).
[·] ·k represents a matrix without the column related to θ k in matrix [·], and ·k = (β −1 I 2N 1 +Ā ·k Λ −kĀ H ·k ) −1 . Λ −k denotes the matrix after removing the row and column related to θ k in the matrix Λ. Generally, Θ k represents the peak search range of the k-th peak, and the search step is usually set to 0.01 or smaller.
Based on the above-mentioned derivations, the robust DOA estimation can be finally realized. The flow of the proposed algorithm is summarized as in Algorithm 1. The correction process ofS; 5: Reconstruct Z = Q(Z −H) according to Equation (14); 6: Estimate µ and Σ by Equations (19) and (20); 7: Update λ new and β new by Equations (24)  The update process ofH; 10: Reconstruct Z = Q(Z −ĀS) according to Equation (15); 11: Estimate µ and Σ by Equation (29) to Equation (30); 12: Updateβ new and Ω new by Equations (31) and (32); 13: ReviseH according to Equation (33); 14: end while 15: ObtainS ; 16: After the spatial spectrum search, a rough estimate of the K peaks by Equation (34) to obtain a fine estimate of the signal source.

Analysis of the Maximum Identified Target Number
Under the signal model shown in Section 3 of this paper, the analysis of the maximum identified number of sources for the algorithm is as follows: First, the left and right sides of Equation (6) where Z contains 2N 1 L measurements and the dimension of S after the null domain discretization is 2K * L (where 2KL positions are large values compared to other positions). SinceS in the above equation is reconstructed fromS andH, the outlier matrixH provides E outliers toS-in summary, when satisfying Equation (36) has a consistent solution. Therefore, the maximum number of identified sources is K ≤ N 1 − (E/L) [24]. In theory, the array expansion shown in Equation (6) in Section 3 constructs a virtual array that doubles the array aperture. However, since the matrix Ξ generated by the decoupling process shown in Equation (5) reconstructs the signal matrixS, the proposed scheme can only identify K ≤ N 1 − (E/L) sources.

Convergence Analysis
The technique used in this paper is a DOA estimation method based on alternating updates of the SBL. The core of the method is that an alternating iterative strategy is used to update the signal power µ l and the noise power µ l . As with the conventional EM-based SBL estimation algorithm, there is good convergence. The specific analysis is presented in [34,35]. The simulation results below verify this property, where the error is set to The physical meaning of this error representation is the difference in signal power estimated from two adjacent iterations. Based on the convergence analysis diagram of the algorithm shown in Figure 2, it can be seen that the error curve tends to decrease asymptotically in one direction as the number of iterations increases. In addition, the trend continues to hold when the signal-to-noise ratio changes. What is clear is that all three curves shown in the graph flatten out and approach zero when the number of iterations is greater than 1000. This phenomenon indicates that the difference in power estimates at this point is small and can be considered constant. The iterative convergence of the algorithm is also demonstrated from a simulation perspective.

Simulation Results
To evaluate the performance of DOA estimation, we directly verify the effectiveness and robustness of the proposed DOA estimation method. Numerous simulation experiments are performed in this section. First of all, the physical meaning of some parameters and environmental settings is explained.
The impulse noise considered in this paper is described by the α-stable distribution with probability density function (PDF), which is where α is the characteristic index. The smaller the value of α, the more obvious the characteristic of impulse noise. When α = 2, it satisfies the Gaussian property. γ is the dispersion parameter, which is against the variance of Gaussian noise and the location parameter ς. In symmetric α-stable distribution, the symmetric parameters are both zero; then, the PDF is rewritten as Since SNR is an invalid description for α < 2, the generalized signal-to-noise ratio (GSNR) is used instead of SNR, which is defined by where δ 2 s denotes signal power, and γ α represents noise power. In order to illustrate the superiority of the proposed algorithm, the RIM-SBLfast algorithm in [30], and the IM-SBLfast algorithm proposed in [24] are introduced to compare with the proposed algorithm. During simulations, unless otherwise specified, the uniform linear array (ULA) is composed of N = 10 array elements, and assume that there are K = 2 uncorrelated signal sources, whose direction is θ = [−10.56 • 9.7 • ]. The mutual coupling coefficient is set to c = [1, 0.9174 + j * 0.3577], the number of snapshots is set to L = 20, GSNR = 10dB, α = 1.2, and the grid interval is set to r = 1 • . In addition, the root mean square error (RMSE) is used to evaluate the performance of different algorithms, which is denoted by where θ k is the true angle of the k-th signal, and ϑ m,k is the corresponding estimated angle in the m-th Monte-Carlo experiment. M represents the total number of Monte Carlo experiments, which is M = 200 in this paper. Simulation results demonstrate the performance of the proposed algorithm from two aspects. Figure 3 illustrates the effectiveness of the proposed algorithm and Figures 4-7 illustrate the robustness of the algorithm. As shown in Figure 3, this results shows that the proposed algorithm is effective to eliminate the mutual coupling and off-grid errors, and the performance of the proposed algorithm is stable under different conditions. Figure 3a shows that the proposed algorithm has relatively accurate DOA estimation when the GSNR changes from 0 dB to 10 dB. Moreover, in the process of increasing the GSNR, the main lobe of the signal spectrum is sharper and the side lobe is lower because the influence of noise is weaker and weaker.  Figure 3b shows that, under the different alpha conditions, our proposed algorithm is effective for the identification and processing of impulse noise. Figure 3c is an experiment conducted under different grid intervals. These results show that the proposed algorithm can reduce the off grid error to a certain extent. However, when the grid interval is set too large, the efficiency of the proposed algorithm will be improved, but some estimation accuracy will be lost.
As shown in Figure 4, the RMSE of different algorithms in terms of GSNR are compared. In order to obtain the validity of the results, the value range of GSNR is set as from −10 dB to 20 dB with a step of 5 dB. The results show that the influence of noise gradually decreases with the increase of GSNR, hence the three curves show an obvious downward trend. In addition, the RMSE of our algorithm has always been smaller than RIM-SBLfast and IM-SBLfast. The advantage of our proposed algorithm is more obvious when the GSNR is higher. When GSNR > 15 dB because the signal is dominant and our algorithm can effectively handle the mutual coupling problem, the performance is also better than the comparison algorithms. Figure 5 shows the performance comparison of different algorithms in terms of characteristic index α. In this simulation, the characteristic index α ∈ [1.1 2], and GSNR = 10 dB. The results show that, when α is smaller, the pulse characteristics are more obvious, and our proposed algorithm has better performance with the increasing of α. However, on the whole, the performance of our algorithm is better than the comparison algorithms in the whole range of α.  Last but not least, in sparse recovery-based DOA estimation algorithms, the high resolution of the DOA estimate often relies on how well the algorithm can attenuate the effects of off-grid errors. Figures 6 and 7 illustrate the performance of the proposed algorithm from different perspectives.   Figure 6 shows the RMSE comparison of the proposed algorithm in terms of grid interval when the GSNR is 0 dB and 10 dB, respectively. Obviously, the proposed algorithm has a smaller RMSE when GSNR = 10 dB because the influence of impulse noise is significantly reduced in this case. However, with the increasing grid interval, DOA estimation performance of the proposed algorithm is obviously reduced. Hence, it is reasonable to set grid interval to 1 • or 2 • .
On the other hand, Figure 7 is the performance comparison of different algorithms in terms of grid interval. According to Figure 7, we can find that the estimation accuracy of these three algorithms decreases with the increase of grid interval (the RMSE of the RIM-SBLfast algorithm decreases when grid = 3 • compared to grid = 2 • , which may be because of the fact that the grid falls just near the true signal). At the same time, it can be seen that our algorithm has the best performance among the other algorithms. This result in Figure 7 effectively illustrates the robustness of the proposed algorithm against off-grid error.

Conclusions
In this paper, a robust DOA estimation scheme based on sparse Bayesian learning has been proposed for non-circular sources in impulse noise and mutual coupling. The proposed algorithm can effectively reduce the influence of array mutual coupling, and can extend the array aperture by using the non-circular property. On the other hand, the proposed algorithm can achieve accurate DOA estimation under impulsive noise while balancing the accuracy and efficiency of DOA estimation. Finally, the robustness and superiority of the proposed algorithm have been demonstrated through numerous simulations.