DOA Tracking Based on Unscented Transform Multi-Bernoulli Filter in Impulse Noise Environment

Aiming at the problem of multiple-source direction of arrival (DOA) tracking in impulse noise, this paper models the impulse noise by using the symmetric α stable (SαS) distribution, and proposes a DOA tracking algorithm based on the Unscented Transform Multi-target Multi-Bernoulli (UT-MeMBer) filter framework. In order to overcome the problem of particle decay in particle filtering, UT is adopted to select a group of sigma points with different weights to make them close to the posterior probability density of the state. Since the α stable distribution does not have finite covariance, the Fractional Lower Order Moment (FLOM) matrix of the received array data is employed to replace the covariance matrix to formulate a MUSIC spatial spectra in the MeMBer filter. Further exponential weighting is used to enhance the weight of particles at high likelihood area and obtain a better resampling. Compared with the PASTD algorithm and the MeMBer DOA filter algorithm, the simulation results show that the proposed algorithm can more effectively solve the issue that the DOA and number of target are time-varying. In addition, we present the Sequential Monte Carlo (SMC) implementation of the UT-MeMBer algorithm.


Introduction
Multi-target Direction of Arrival (DOA) estimation is an essential issue in array processing and has a wide range of applications in source location, radar, sonar, and wireless communications [1,2]. Sparse representation and compressive sensing methods are used for DOA estimation of coprime array [3][4][5][6], while these methods are only applied in the case where the sources are stationary. In addition, difficulties also arise from the uncertainties of the source dynamics: the source may be moving or static. Thus, it is significant to extend the static DOA estimation algorithm to the dynamic DOA tracking algorithm.
The representative dynamic DOA tracking algorithms include the subspace tracking algorithm and the particle filter (PF) algorithm. The subspace tracking algorithm includes Projection Approximation Subspace Tracking (PAST) [7] and the Projection Approximation Subspace Tracking with Deflation (PASTD) [8]. In essence, these algorithms transform the determination of the eigensubspace into solving an unconstrained optimization problem, and combine the recursive least squares (RLS) theory to achieve effective tracking of the eigensubspace of time-varying sources. However, the RLS method is very sensitive to impulse noise, and the PAST algorithm's subspace tracking performance will degrade sharply in the impulse noise environment [9][10][11]. In an army of acoustic applications, such as underwater and room acoustic signal processing, the noise environment is non-Gaussian and is impulsive in nature [12,13]. Under investigation, it was found that α stable distribution (0 < α ≤ 2) is a Z(t) = A(θ)S(t) + N(t) (1) where N M×1 (t) = [n 1 (t), n 2 (t), . . . , n M (t)] T represents the impulsive noise vector which is not correlated with signals. Z M×1 (t) = [z 1 (t), z 2 (t), . . . , z M (t)] T is the measurement at time t, A M×P (t) = [a(θ 1 ), a(θ 2 ), . . . , a(θ P )] T is array manifold, S P×1 (t) = [s 1 (t), s 2 (t), . . . , s P (t)] T denotes the acoustic sources matrix, and a θ p = 1, e − j 2π λ d sin θ p , . . . , e − j 2π is the steering vector with λ denoting the wavelength of the carrier, and d is the array space.

α Stable Distribution
Most of the traditional research methods estimating the DOA are based on Gaussian noise models. In practical situations, such as radar echo and low-frequency atmospheric noise, they consist of impulse noise with a short duration and large amplitude. The performance of the algorithm will drop significantly when the Gaussian noise model is still modeled in an impulse noise environment. The α stable distribution is a good example of such a type with significant spike noise and a Gaussian distribution. The α stable distribution's probability function does not have the closed form, which can be conveniently described by its characteristic function as where = tan απ 2 , α 1 2 π log|t|, α = 1 (4) α is the characteristic exponent, whose size can affect the degree of impulse and the range is 0 < α ≤ 2. γ is a dispersion parameter whose mean is consistent with the variance of the Gaussian distribution. β is a symmetric parameter, and the distribution at β = 0 is a symmetric α stable (SαS) distribution. a is the positional parameter. When α = 2, β = 0, it is a Gaussian distribution model. When α = 1, β = 0, it is the Cauchy distribution model. When α = 1/2, β = −1, it is the Pearson distribution model. A crucial difference between the Gaussian distribution and the α stable distribution is that the latter does not have second-order statistics so that its covariance is inaccurate.

Multi-Target Bayesian Theory
Assume that the state of the sources at time k is x k = θ k , . θ k T , where θ k is the DOA and moves at a speed of . θ k rad/s. The state and number of sources are changing at time k + 1, which can be described by RFS. From [20], the sources state set in multiple sources tracking can be regarded as an RFS, namely where X k represents a set of sources at time k, and the element of the set may be one or more or null. Z k denotes the measurement set generated by all sources received time k, and the element is only one. Single-target Bayesian filtering can be extended to multi-target tracking by modeling the above source states and measured values. The single target posterior probability density function (pdf) p k|k (x k |Z 1:k ) is replaced by the joint multi-target posterior p k|k (X k |Z 1:k ). The Bayes joint filter recursion includes two stages: prediction and update. The prediction and update at time k in [24] are and where δ is the set integral and Z 1:k−1 represents all the measurement sets up to time k − 1. g(Z k |X k ) is a multi-target joint likelihood function and f k|k−1 (X k |X k−1 ) is a multi-target state transition probability density function. p k|k−1 (X k |Z 1:k−1 ) represents the multi-target joint prediction probability density and p k|k (X k |Z 1:k ) is the multi-target joint posterior probability density function.

Multi-Target Multi-Bernoulli Filter
A Bernoulli set X has a probability 1 − r of being a null set, and has a probability r of containing a single element x that is distributed via a pdf s(·). The probability of a Bernoulli RFS can be expressed in [21] as A Multi-Bernoulli RFS X can be considered as union of a fixed number of independent Bernoulli sets that have existence probability r ( j) ∈ (0, 1), j = 1, . . . J and the pdf s ( j) , such that where the j th Bernoulli set is described by its two parameters: the existence probability r ( j) and the pdf s X ( j) . So a Multi-Bernoulli RFS can be characterized by a posterior parameter set r where J k|k indicates the number of sources. Z k = z 1,k , z 2,k , . . . , z M,k T denotes the sensor measurement data and Z k ∈ Z, in which Z is the measurement space of the sensor. Target birth and survival are determined by birth probabilities p b,k (X k ) and survival probabilities p s,k (X k ), respectively. The source motion model is represented by the transition probability density f k|k−1 (X k |X k−1 ), and the prior probability of Multi-Bernoulli is described as According to Equation (7), the prediction part can be described as where J k|k−1 = J P,k|k−1 + J B,k , J P,k|k−1 = J k−1 . The number of Multi-Bernoulli parameter sets for survival sources and newborn sources are represented by J P,k|k−1 and J B,k , respectively. According to Equation (8), if the predicted Multi-Bernoulli parameter set can be expressed as r where g(Z k |X k ) denotes the likelihood function. If the covariance of the general sensor array at time k in Gaussian noise environment is R k , the likelihood function can be expressed as The frame of Formula (18) is not held in impulse noise, so we propose to replace the likelihood function with a spatial spectrum method.

Improved Algorithm for Likelihood Function
In the practical engineering application, to guarantee the real-time and effectiveness of the estimation, the observation matrix of the array is obtained with a limited number of snapshots. Assuming L observations at time k, the array covariance matrix is calculated asR k = X(t k )X(t k ) H /L. We assume that the noise vector N(t) is independent to the target signal and has a SαS distribution with a characteristic exponent of α. From [25], if the array observation matrix Z k at time k is obtained, the FLOM matrix is defined as where ψ i,j represents the (i, j)th element of Ψ k , and (·) * represents conjugate operation. The dimension of matrix Ψ k is M × M. In [25], the authors derived the form of the FLOM matrix as where R s and r represent the source and additive noise of the FLOM matrix, respectively. As can be seen from Equation (20), the (i, j)th FLOM matrix element is defined as Fractional moment p must satisfy 1 < p < α ≤ 2. The FLOM is used to replace the covariance matrix of the signal in impulse noise, and then the eigendecomposition is performed on Ψ k in the MUSIC algorithm to obtain the noise subspace U n . The form of the FLOM-MUSIC spatial spectrum estimation function is where C = [1, 0], and the CX k represents source azimuth information. a(·) is a space vector, and ζ ∈ R + represents an exponential weighting of the spatial spectrum. The response of the traditional MUSIC spatial spectral beamformer in an impulse noise environment is distorted, which can result in a significant degradation in the performance of the resampling step. After being weighted, the particles can be moved to the high likelihood region to the resampling performance.

UT-MeMBer DOA Particle Filter Tracking Algorithm
In this section, we describe the particle filter implementation of the UT-MeMBer algorithm.
From [22], if the multi-target probability density parameter set at time k − 1 is r , then the spatial posterior probability density at time k − 1 and can be expressed as: where s j k−1|k−1 is the spatial posterior probability density, which can be approximated as the weighted k−1 represents the state of the i th particle, including angle and speed, i.e., x k−1 denotes the weight, usually satisfying According to (12), the spatial posterior probability density of the prediction step consists of two items and can be written as where N k|k−1 = N k−1 + N B,k and J k|k−1 = J p,k|k−1 + J B,k represent the number of predicted particles and predicted MeMBer parameter sets, respectively. All particles can be sampled from two parts: Among them, N B,k denotes the number of newborn particles at time k, x k|k−1 [13]. Particle filtering suffers from missing sample diversity, resulting in depletion of the sampled particles. In order to solve this problem, the surviving particles will be subjected to UT operations. A set of sigma points with different weights are selected by UT operation, and then the posterior probability density of the state is approximated by these sigma points. The weight is where p s and p b represent the survival probability of particles and the newborn probability of particles, respectively. N k−1 is the number of surviving particles sampled from the transition probability density f k|k−1 , and B is the number of newborn particles from the proposal probability density β k . If the prediction MeMBer parameter sets can be expressed as r then the update MeMBer parameter sets can be written as r j k|k , ω . The weight is where p D,k is the detection probability, and the likelihood function g Z k x (i,j) k|k−1 calculated by the MUSIC algorithm can be expressed as k|k−1 represents the azimuth angle information, ζ is the exponential weighting factor. U n represents the noise subspace obtained by the MUSIC algorithm. The steps of the UT-MeMBer DOA particle filter tracking algorithm are shown in Algorithm 1.
Algorithm 1 UT-MeMBer DOA particle filter tracking algorithm Input: Predict the existence probability: r k−1 denotes the existence probability of survival model, k−1 represents the existence probability of newborn model.

2.
Calculate the predicted state of surviving particles: -Select a weighted sample point of 2n x + 1 for each particle x Construct a newborn target weighted particle: Calculate the prediction weight ω (i,j) k|k−1 , i = 1, . . . , N k|k−1 according to (26). 5.
Unite weighted particle set:
Update existence probability:

8.
The updated weight is calculated by (27) and normalized ω

Resample
Step Algorithm 1 gives the pseudo-code of UT-MeMBer DOA particle filter tracking algorithm. The prediction is made in steps 1-5.
Step 6 calculates each predicted particle likelihood function which is replaced by the MUSIC spatial spectral function. The update existence probability is calculated in step 7.
Step 8 calculates the normalized weight. Particle resampling is performed in step 9. The particle set ω approximates the spatial probability density function s j k|k , and the estimation of updated source can be expressed as k .

Simulation Results
Since the traditional MUSIC algorithm cannot solve the multi-source tracking problem when target number is varying, this paper uses FLOM matrix to substitute the covariance matrix to obtain the corresponding MUSIC spatial spectrum, which can be as the particle likelihood function. We proposed a UT-MeMBer DOA tracking algorithm under RFS framework, which can be named as UT-MB-FLOM-MUSIC algorithm. The Generalized Signal to Noise Ratio (GSNR) is defined as where γ represents the noise dispersion parameter, and GSNR represents the ratio of signal intensity and noise dispersion. In the simulation, different characteristic indices α describe the degree of impact of different noises. In the following simulation experiments, the estimated performance is evaluated by the root mean square error (RMSE), which is defined as θ k (t) rad/s, the constant velocity (CV) model is given as where the transfer matrix F k and G are defined by respectively, where ∆T = 1s denotes the time step, and v k is a zero-mean real Gaussian process used to model the disturbance on the source velocity, i.e., v k ∼ N (0, Σ k ) with Σ k = 1. Experimental conditions are as follows: The number of array elements is M = 10, d = λ/2, the observation time is K = 50 s , L = 100, GSNR = 10 dB, MC = 100, and ξ = 5. The source survival probability p s,k (x k ) = 0.99, and the source detection probability p D,k (x k ) = 0.98. In the UT-MB-FLOM-MUSIC algorithm prediction step, we assume that there are six new sources at each time, i.e., J B,k = 6, all obeying a uniform distribution on [−π/2 , π/2] and each new source produces 300 particles, i.e., N B,k = 300. In the update step, the MUSIC spatial spectral function is used to replace the likelihood function and is exponentially weighted, which improves the feasibility of the algorithm.
In the impulse noise model, the noise is Gaussian noise when α = 2. The DOA estimation method  Figure 1a shows the RMSE of angles for four algorithms when running 100 MC at α = 2, GSNR = 10 dB, and Figure 1b shows two source trajectories for a single MC. It can be seen from Figure 1a that the UT-MB-FLOM-MUSIC algorithm proposed in this paper is obviously better than the traditional PASTD and has the highest accuracy when the number of targets is constant. It can be seen in Figure 1b that the algorithm can effectively track the target trajectory, while the PASTD algorithm deviates from the real trajectory at several times.

Scenario 2: The Number of Targets Is Time-Varying
Consider a linear multi-source scenario with three sources. The number of sources is timevarying due to births and deaths, the survival time of the four sources is 1-50 s, 10 Figure 3a shows the RMSE of angles for three algorithms for running 100 MC at α = 2, L = 100 and GSNR = 10 dB, and Figure 3b shows three sources trajectory for a single MC. It can be seen from Figure 3 that the likelihood function of the MUSIC spatial spectrum instead of the Multi-Bernoulli particle filter update stage can effectively estimate the target number and motion state, and also verify the feasibility of the literature [14] in the Gaussian noise environment. Although the error is large at We show the RMSE for tracking the multi-source motion when α = 1.3, GSNR = 10 dB, MC = 100, and L = 100 in Figure 2a. It can be seen from Figure 2a that the RMSE of the UT-MB-FLOM-MUSIC algorithm is smaller than that of the other three algorithms. The accuracy of the MB-MUSIC algorithm is significantly reduced in impulse noise, and the PAST algorithm is more accurate than MB-MUSIC. It can be seen from Figure 2b that the MB-MUSIC algorithm cannot effectively track the target trajectory in impulse noise, and the PASTD algorithm also has the problem of inaccurate target tracking. Based on the fact that the above target numbers are unchanged, we will analyze the target time-varying DOA tracking.  Figure 3a shows the RMSE of angles for three algorithms for running 100 MC at α = 2, L = 100 and GSNR = 10 dB, and Figure 3b shows three sources trajectory for a single MC. It can be seen from Figure 3 that the likelihood function of the MUSIC spatial spectrum instead of the Multi-Bernoulli particle filter update stage can effectively estimate the target number and motion state, and also verify the feasibility of the literature [14] in the Gaussian noise environment. Although the error is large at time 35, the overall error is below 2 degrees. It can also be seen from Figure 3a that the RMSE of the UT-MB-FLOM-MUSIC algorithm is also smaller than other algorithms even in the Gaussian noise environment.

Scenario 2: The Number of Targets Is Time-Varying
Consider a linear multi-source scenario with three sources. The number of sources is timevarying due to births and deaths, the survival time of the four sources is 1-50 s, 10-50 s, 20-45 s, and the initial source states are Figure 3a shows the RMSE of angles for three algorithms for running 100 MC at α = 2, L = 100 and GSNR = 10 dB, and Figure 3b shows three sources trajectory for a single MC. It can be seen from Figure 3 that the likelihood function of the MUSIC spatial spectrum instead of the Multi-Bernoulli particle filter update stage can effectively estimate the target number and motion state, and also verify the feasibility of the literature [14] in the Gaussian noise environment. Although the error is large at time 35, the overall error is below 2 degrees. It can also be seen from Figure 3a that the RMSE of the UT-MB-FLOM-MUSIC algorithm is also smaller than other algorithms even in the Gaussian noise environment. Since Gaussian noise does not reflect true signal interference, the α stable distribution can reflect the impact of impulse noise. Figure 4 shows the RMSE and cardinality estimation error plots for three algorithms running 100 MC when the characteristic index α is different and the GSNR = 10 dB, L = 100. It can be seen from Figure 4a that, in α = 1.1~1.9 , the RMSE error of the three estimation algorithms first decreases, and finally tends to be flat. It also can be seen that the RMSE of the UT-MB-FLOM-MUSIC algorithm is significantly smaller than the MB-FLOM-MUSIC and MB-MUSIC algorithms when α = 1.1 or α = 1.2, so that the UT-MB-FLOM-MUSIC algorithm has a better effect when handling the impulse noise environment. Since the characteristic index is close to 2 when α = 1.8 or α = 1.9, Figure 4b shows that the cardinality estimation error of the three algorithms approaches 0. It also shows that it is feasible to use the MUSIC spatial spectrum as a substitute for the likelihood function when the noise environment is close to Gaussian noise while the MUSIC algorithm cannot effectively estimate the number of targets in an impulse noise environment. Since Gaussian noise does not reflect true signal interference, the α stable distribution can reflect the impact of impulse noise. Figure 4 shows the RMSE and cardinality estimation error plots for three algorithms running 100 MC when the characteristic index α is different and the GSNR = 10 dB, L = 100. It can be seen from Figure 4a that, in α = 1.1 ∼ 1.9, the RMSE error of the three estimation algorithms first decreases, and finally tends to be flat. It also can be seen that the RMSE of the UT-MB-FLOM-MUSIC algorithm is significantly smaller than the MB-FLOM-MUSIC and MB-MUSIC algorithms when α = 1.1 or α = 1.2, so that the UT-MB-FLOM-MUSIC algorithm has a better effect when handling the impulse noise environment. Since the characteristic index is close to 2 when α = 1.8 or α = 1.9, Figure 4b shows that the cardinality estimation error of the three algorithms approaches 0. It also shows that it is feasible to use the MUSIC spatial spectrum as a substitute for the likelihood function when the noise environment is close to Gaussian noise while the MUSIC algorithm cannot effectively estimate the number of targets in an impulse noise environment.
In Figure 5, we show the RMSE and cardinality estimation for tracking the multi-source motion when α = 1.3 and GSNR = 10 dB, MC = 100. It can be seen from Figure 5 that the RMSE of the UT-MB-FLOM-MUSIC algorithm is smaller than that of the other two algorithms. Although the RMSE will increase when the new target appears or disappears, it will decrease rapidly at the next time step. This phenomenon shows that the Multi-Bernoulli filter has a large recognition performance for the target and can quickly track the state of the target. Table 1  when handling the impulse noise environment. Since the characteristic index is close to 2 when α = 1.8 or α = 1.9, Figure 4b shows that the cardinality estimation error of the three algorithms approaches 0. It also shows that it is feasible to use the MUSIC spatial spectrum as a substitute for the likelihood function when the noise environment is close to Gaussian noise while the MUSIC algorithm cannot effectively estimate the number of targets in an impulse noise environment.  In Figure 5, we show the RMSE and cardinality estimation for tracking the multi-source motion when α = 1.3 and GSNR = 10 dB, MC = 100. It can be seen from Figure 5 that the RMSE of the UT-MB-FLOM-MUSIC algorithm is smaller than that of the other two algorithms. Although the RMSE will increase when the new target appears or disappears, it will decrease rapidly at the next time step. This phenomenon shows that the Multi-Bernoulli filter has a large recognition performance for the target and can quickly track the state of the target. Table 1   The operating environment includes an Intel (R) Core (TM) i5-8500 CPU @ 3.00 GHz processor and a 64-bit operating system MATLAB 2014. It can be seen from Table 1 that the UT-MB-FLOM-MUSIC algorithm RMSE is smaller than other algorithms when the running time is too long.    The operating environment includes an Intel (R) Core (TM) i5-8500 CPU @ 3.00 GHz processor and a 64-bit operating system MATLAB 2014. It can be seen from Table 1 that the UT-MB-FLOM-MUSIC algorithm RMSE is smaller than other algorithms when the running time is too long.  . let ε = 1 . It can be seen from Figure 6a that the MB-FLOM-MUSIC and UT-MB-FLOM-MUSIC algorithms have higher accuracy than the MB-MUSIC in an impulse noise environment, and the UT-MB-FLOM-MUSIC algorithm has higher accuracy under the high GSNR. It can be seen form Figure 6b that as the SNR increases, the PROC increases. And at the same GSNR, the performance of the MB-FLOM-MUSIC algorithm is more significant.

Scenario 3: The Number of Targets Is Time-Varying and Maneuvering
Consider a nonlinear multi-source scenario with three sources. The number of sources is timevarying due to births and deaths, and the survival time of the three sources is 1-50 s, 10 where ω = rad 0.25 and other experimental conditions are the same as scenario 1. Figure 8 shows the maneuvering target trajectory of three algorithms running one MC when α = 1.3, L = 100, and GSNR = 10 dB. It can be clearly seen from Figure 8 where ω = 0.25 rad and other experimental conditions are the same as scenario 1. Figure 8 shows the maneuvering target trajectory of three algorithms running one MC when α = 1.3, L = 100, and GSNR = 10 dB. It can be clearly seen from Figure 8    In Figure 9, we show the RMSE and cardinality estimation for tracking the multi-source motion when α = 1.3 and GSNR = 10 dB, MC = 100. It can be seen from Figure 9a that the RMSE of the UT- In Figure 9, we show the RMSE and cardinality estimation for tracking the multi-source motion when α = 1.3 and GSNR = 10 dB, MC = 100. It can be seen from Figure 9a that the RMSE of the UT-MB-FLOM-MUSIC algorithm is smaller than that of the other two algorithms. As can be seen from  Table 2 shows the RMSE and computing performance of the MB-MUSIC algorithm, MB-FLOM-MUSIC algorithm and the UT-MB-FLOM-MUSIC algorithm. Compared with the results in Table 1, the RMSE and running time of the three algorithms are increased when the target is maneuvered. The RMSE of UT-MB-FLOM-MUSIC algorithm is smaller than other two algorithms when the running time is long.

Conclusions
A DOA tracking algorithm based on the UT-MeMBer particle filter in an impulse noise environment is proposed in this paper. Since the FLOM matrix is used instead of the covariance matrix, the spatial spectrum based on FLOM can well reflect the real DOA in impulse noise environment. For the persistent surviving particles, the sigma point is selected by UT to approximate  Table 2 shows the RMSE and computing performance of the MB-MUSIC algorithm, MB-FLOM-MUSIC algorithm and the UT-MB-FLOM-MUSIC algorithm. Compared with the results in Table 1, the RMSE and running time of the three algorithms are increased when the target is maneuvered. The RMSE of UT-MB-FLOM-MUSIC algorithm is smaller than other two algorithms when the running time is long.

Conclusions
A DOA tracking algorithm based on the UT-MeMBer particle filter in an impulse noise environment is proposed in this paper. Since the FLOM matrix is used instead of the covariance matrix, the spatial spectrum based on FLOM can well reflect the real DOA in impulse noise environment. For the persistent surviving particles, the sigma point is selected by UT to approximate the posterior density of the state to improve the performance of the particle. Then, the MUSIC spatial spectral function of the FLOM matrix is used to represent the likelihood function of the particle. And the weighting of the likelihood function can further increase the weight of the particles in the high likelihood region. The results show that the UT-MB-FLOM-MUSIC algorithm is more effective than the PASTD, MB-MUSIC, and MB-FLOM-MUSIC algorithms in an impulse noise environment. The advantage of this algorithm is that the MeMBer filter can operate the array data more directly, and can effectively track the target number of time-varying DOA. The shortcoming of this algorithm is that it takes a long time. Our future work will focus on how to improve the efficiency of the algorithm, maneuvering target tracking in other noisy environments, etc.