Variational Bayesian cardinalized probability hypothesis density filter for robust underwater multi-target direction-of-arrival tracking with uncertain measurement noise

The direction-of-arrival (DOA) tracking of underwater targets is an important research topic in sonar signal processing. Considering that the underwater DOA tracking is a typical multi-target problem under unknown underwater environment with missing detection, false alarm, and uncertain measurement noise, a robust underwater multi-target DOA tracking method for uncertain measurement noise is proposed. First, a kinematic model of the multiple underwater targets and bearing angle measurement model with missing detection and false alarms are established. Then, the multi-target DOA tracking algorithm is derived by using the cardinalized probability hypothesis density (CPHD) filter, the performance of which largely depends on the accuracy of the parameter of measurement noise variance. In addition, the variational Bayesian approach is used to adaptively estimate the uncertain measurement of noise variance for each measurement of target in the real time of tracking. Thus, the robust underwater multi-target DOA tracking is carried out. Finally, comprehensive experimental validations and discussions are made to prove that the proposed algorithm can provide robust DOA tracking in the multi-target tracking scenario with uncertain measurement noise.


Introduction
The direction-of-arrival (DOA) estimation and tracking is an important research topic in sonar signal processing [1][2][3][4][5]. For the scenario of moving targets, the traditional DOA estimation methods cut the measurement of the output of sonar array signal into small periods in time to process, which ignores the kinematic characteristics of the targets [6][7][8]. The DOA tracking methods not only use the measurement information but also rely on the kinematic characteristics of the underwater targets [9][10][11][12][13][14][15][16][17][18]. Therefore, by not only depending on the current measurements but also utilizing the prior kinematic information of an unknown underwater target, the DOA tracking methods can provide more robust and accurate results than the traditional DOA estimation methods.
In view of the advantages of the DOA tracking methods, researchers have carried out a lot of work around DOA tracking. The DOA tracking methods first take the outputs of the sonar array signals or results of the traditional DOA estimation methods as the measurements. Then, considering the typical kinematics of the underwater targets, the kinematic model by the bearing angle is established. Then, based on the framework of the Bayesian filter theory [19], the bearing angle of the target can be recursively estimated from the current measurements. Depending on the bearing angle measurements obtained by the traditional DOA estimation methods, the Kalman filter (KF) is always utilized as the DOA tracking algorithm [9,10] for the linear relationship between the measurement and bearing angle of the target. Besides utilizing the bearing angle measurements, some research take the outputs of sonar array signals as measurements. However, under such a circumstance, the mathematical relationship between the bearing angle of the target and the measurement is non-linear. For such nonlinear problems, many tracking algorithms based on the non-linear Bayesian filter have proposed the extended Kalman filter (EKF) [11,12].The first-order Taylor series expansion of the measurement model was used to approximate the non-linear model to a linear model. The unscented Kalman filter (UKF) [20,21] uses a group of determined sigma points to linearize the non-linear model. The particle filter (PF) [22] generates and transforms a large number of arbitrary particles to represent the distribution of the target state. Thus, the PF can theoretically achieve accurate tracking in any non-linear tracking system with any distribution of the uncertainties, while the computational complexity of the PF stays significant.
Although the tracking techniques are widely utilized in the scenario of underwater target tracking, most of them are only applicable to the single-target tracking scenario [9][10][11][12][20][21][22]. However, the underwater DOA tracking issue is a typical multi-target tracking scenario where the multi-target tracking techniques should be considered and proposed. The methods for multi-target tracking problem can be divided into two categories: the traditional data association methods and random finite set (RFS)-based multi-target tracking methods. The data association methods establish the association between measurements and targets, so as to transform the multi-target tracking problem into a single-target tracking problem to use the abovementioned single-target tracking methods. However, these data association techniques have to match every measurement with its target which can make computation highly complex with large number of targets or false alarms caused by uncertain measurement noise [23,24]. From the beginning of this century, the RFS-based multi-target tracking methods have developed rapidly to overcome the drawbacks caused by the data association techniques. The RFS can be defined by a set with elements along with the number of the elements which are subjected to random distributions. The RFS-based multi-target tracking methods define the states and measurements of targets as RFSs such that the data association procedure can be avoided. As a result, by utilizing the RFS technique, the computational complexity of multi-target tracking can be hugely reduced, especially when the number of targets and false alarms is large. For the RFS-based multitarget tracking, Mahler first proposed the concept of "first-order moment filter," also known as probability hypothesis density (PHD) filter [25]. Since, the PHD filter has no closed-form solution in general, Vo et al. [26] proposed sequential Monte Carlo solution for the PHD filter (SMC-PHD), and Clark et al. [27] proposed the Gaussian mixture model implementation of the PHD filter (GM-PHD), which push the PHD filter from theory to application. The GM-PHD filter requires both the kinematic model of targets and the measurement model to be linear. The SMC-HD filter represents the distribution of the states of the targets by generating and transferring many of the particles, which is suitable for non-linear models. However, the calculation of the particles leads to high computation overload. Subsequently, Mahler [28,29] introduced the cardinality distribution to describe the number of targets on the basis of the PHD filter to make a more accurate estimation of the number of targets and proposed the CPHD filter. Similar to the PHD filter, the CPHD filter also has no closed-form solution in general and has to be implemented on the basis of the SMC or GM model [30]. Considering that the computational complexity of the RFS-based method is lower than that of the data association method, researchers have proposed many multi-target DOA tracking methods based on the PHD and CPHD filters [13][14][15][16]. These methods use the output of the array signal as measurements with the non-linear measurement model, thus implementing them based on SMC.
Due to high-dimension array signal-based measurements and the large number of particles for the SMC method, the existing RFS-based multi-target DOA tracking methods get a large computation overload. Moreover, in the real scenario of underwater tracking, the unknown ocean environment always results in uncertain measurement noise. Thus, besides considering the computational complexity, accomplishing robust and accurate tracking when the measurements are in a low SNR scenario is also important, especially in the underwater target tracking case. However, most of the existing DOA tracking methods assume the measurement noise to be a certain stochastic process which leads to the degradation of the tracking performance in the real scenario of uncertain measurement noise. To deal with the uncertain measurement noise, Zhang et al. [12] derived a robust singletarget DOA tracking method based on the EKF by estimating the measurement noise covariance matrix (MNCM) by using the improved Sage-Husa algorithm. By estimating the MNCM in a maximumlikelihood (ML) framework, the expectation-maximum adaptive Kalman filter (EM-AKF) was proposed [31,32]. Sarkka and Hartikainen [33,34] assumed the MNCM to be subject to an inverse Wishart distribution and iteratively estimated the MNCM by using variational Bayesian approach, namely, the variational Bayesian adaptive Kalman filter (VB-AKF). Huang et al. [37] assumed the mean square error matrix (MSEM) to also be subjected to an inverse Wishart distribution and jointly estimated it along with the MNCM to improve the performance of VB-AKF. In the existing works, the superior accuracy and stability of the VB-AKF has been demonstrated [33][34][35][36][37]. However, the MNCM estimation method has not been applied to robust multi-target DOA tracking yet. Thus, considering the uncertain underwater environment and its influences on the underwater target tracking missions, adopting variational Bayesian online estimation technique into the multitarget tracking scenario is inspiring and necessary.
In this article, the bearing angle estimates obtained by using the traditional DOA estimation method are regarded as the measurements, and the multi-target DOA tracking algorithm is derived by using the GM-CPHD filter. Different from most of the existing tracking algorithms, we considered the uncertain measurement noise caused by the unknown underwater environment, which always makes a certain tracking algorithm diverge. Thus, the variational Bayesian approach is utilized to estimate the covariance matrix of the measurement noise along with the states of targets in the framework of the GM-CPHD filter. In this way, the variational Bayesian GM-CPHD filter (VB-GMCPHD) for robust underwater multi-target DOA tracking is proposed for the scenario of robust tracking under uncertain measurement noise. Finally, the results of the experiment show the robustness and accuracy of the proposed method in real underwater multi-target DOA tracking scenario. The contributions of this article are summarized as follows: First, the multi-target DOA tracking algorithm is derived by using the GM-CPHD filter for the real underwater tracking scenario with missing detection and false alarm.
Second, the issue of uncertain measurement noise is addressed by using the variational Bayesian approach to estimate the measurement noise variance. Thus, the VB-GMCPHD for robust underwater multi-target DOA tracking is proposed.
Finally, the real experimental data is utilized to verify the superior accuracy and robustness of the proposed method in real underwater DOA tracking scenario.
The rest of this article is organized as follows: in Section 2, the problem of underwater multi-target DOA tracking with missing detection and false alarm is formulated. In Section 3, the GM-CPHD filter for DOA tracking is described. In Section 4, based on the variational Bayesian approach, an innovative multi-target DOA tracking method for uncertain measurement noise scenario, namely the VB-GMCPHD, is proposed. In Section 5, experimental validations are made and the results are proved, demonstrating the superior performance of the proposed VB-GMCPHD. Finally, the conclusions are drawn in Section 6.

Multi-target tracking models
The sets of the states of the targets and the measurements are assumed to be RFSs. Assuming that n k targets exist within the detection range of the sonar at time step k and the state of the i n -th target is x in k , the set of the states of targets at time step k is expressed as X k x 1 k , x 2 k , ..., x n k k . Let θ k denote the bearing angle of the i n -th target (θ in k is the angle between the target and positive x-axis with respect to the positive counterclockwise.) and _ θ in k denote the change rate of θ in k , then the state of the i n -th target is expressed as where (·) T denotes the matrix transposition. Since the underwater targets are usually not maneuvering to save energy and keep concealed, the bearing angles of the targets are assumed to be subject to the constant velocity (CV) model. The CV model of the state of the i n -th target is expressed as where F k|k−1 and G k are the state transition matrix and the noise driving matrix, w k denotes the zero-mean Gaussian process noise with the covariance matrix Q k caused by the unknown underwater environment, and T is the interval between the adjacent time steps. The estimates of the bearing angles of the targets with error obtained by the traditional DOA estimation methods are taken as measurements. Considering the probability of missing detection and false alarm, the number of measurements and the states of the targets are always different in a multi-target tracking problem. It is assumed that m d k targets are detected and m f k false alarms exist, i.e., m k m d k + m f k measurements are obtained at time step k. The set of measurements at time step k are expressed as where H k denotes the measurement matrix and H k [1,0], v k is the error of the bearing angle estimate obtained by using traditional DOA estimation, which is subject to Gaussian distribution with zero mean and variance of σ 2 r,k . The false alarms z 1 k , z 2 k , ..., z m f k k are assumed to be subject to Poisson distribution with intensity of κ k .

GM-CPHD filter for multi-target DOA tracking
The CPHD filter performs multi-target tracking by recursively calculating the PHD and the cardinality distribution to represent the distribution of the states and the number of the targets, respectively. The closed-form solution of the CPHD filter is given in the assumption of the linear Gaussian mixture (GM) model, which is called the GM-CPHD filter. Each component of the Gaussian mixture model represents the respective states of the targets [30].
The cardinality distribution p k−1 (n) and the PHD v k−1 (x) are assumed to be known, and v k−1 (x) is subject to the Gaussian mixture model as where w (i) k−1 denotes the weight, N(·; m, P) denotes the Gaussian distribution with the mean of m and covariance matrix of P, and m (i) k−1 and P (i) k−1 denote the estimate of the state of the target and the mean square error matrix (MSEM), respectively. Then, the process of the GM-CPHD filter at time step k is divided into prediction and update as follows: (1) Prediction Once the cardinality distribution p k−1 (n) at time step k − 1 is given, the predicted cardinality distribution is expressed as where is the cardinality distribution of birth targets, and p s,k is the probability of targets surviving.
Once the PHD v k−1 (x) at time step k − 1 is given, the predicted PHD is expressed as where v s,k|k−1 (x) and γ k (x) denote the PHD of surviving targets and birth targets, respectively. v s,k|k−1 (x) is expressed as where the predicted state m (j) s,k|k−1 and predicted MSEM P (j) s,k|k−1 of the j-th surviving target is given as follows: The PHD of birth targets γ k (x) is also subject to the Gaussian mixture model as follows: where J γ,k denotes the number of components of the Gaussian mixture model for γ k (x), and w (i) γ,k , m (i) γ,k , and P (i) γ,k denote the weight, state, and MSEM of the i-th birth target, respectively.
According to Equations 5, 6, and 9, the predicted PHD can be expressed as follows: where w (i) k|k−1 , m (i) k|k−1 , and P (i) k|k−1 denote the weight, predicted estimate of the state, and predicted MSEM of the i-th target, respectively.
(2) Update The cardinality distribution p k (n) and PHD v k (x) at time step k are obtained by using the measurement set Z k to update p k|k−1 (n) and v k|k−1 (x) as follows: where 〈α, β〉 denotes the inner product of α and β, i.e., 〈α, β〉 , and p D,k denotes the detection probability, and where |Z| denotes the number of the elements of the set Z, p K,k denotes the cardinality distribution of false alarm at time step k, e j (·) denotes elementary symmetric function of order j, and e j (Z) At last, the components with tiny weight are pruned away, the components with uniform distribution are combined, and the maximum number of components is limited [27]. Then, the updated PHD at time step k is expressed in the Gaussian mixture model as follows: where w (i) k , m (i) k , and P (i) k denote the weight, estimate of the state, and MSEM of the i-th target at time step k, respectively.
The n corresponding to the maximum value of the cardinality distribution p k (n) is the estimate of the number of targets. The m (i) k corresponding to the components with theN k largest weight in the PHD v k (x) is the estimate of the target states. The first element of the estimated state vector is the estimate of the bearing angle of the target. By substituting p k (n) and v k (x) into the next time step, the GM-CPHD filter can be used to recursively estimate the state of the target to carry out DOA tracking. 4 VB-GMPHD filter for robust multitarget DOA tracking with uncertain measurement noise The existing DOA tracking techniques usually assume that the measurement noise is a certain stochastic process, which means the MNCM is constant. However, in the real scenario of underwater tracking, the unknown underwater environment always results in uncertain measurement noise, which means the MNCM is time varying. Thus, the assumption of the existing DOA tracking technology on constant MNCM is inconsistent with the real scenario, which results in the decline of tracking performance. In order to improve the robustness of multi-target DOA tracking, variational Bayesian approach is used to estimate the MNCM in the real time of tracking. In this article, the measurement is a number instead of a vector, thus the MNCM reduces to the measurement noise variance.

Choice of prior distribution
In the multi-target DOA tracking scenario, when using the i m -th measurement to estimate the corresponding i n -th target, the one-Frontiers in Physics frontiersin.org step predicted probability density distribution (PDF) p(x in k |z im 1: k−1 ) and the likelihood PDF p(z im k |x in k ) are assumed to be subject to Gaussian distributions in the framework of the KF [35] as follows: where N(·; μ, Σ) denotes the PDF of the Gaussian distribution with mean μ and covariance matrix Σ, H k is the measurement matrix given in Equation 2, and σ 2 r,k denotes the measurement noise variance.x in k|k−1 and P in k|k−1 denote the predicted state and predicted MSEM of the i n -th target, respectively, which are expressed asx wherex in k−1|k−1 and P in k−1|k−1 are the estimates of the state and MSEM of the i n -th target at time step k-1, respectively.
In order to infer x in k along with σ 2 r,k , a conjugate prior distribution has to be selected for the fluctuant σ 2 r,k , since a conjugate distribution can guarantee the same functional forms of the prior distribution and posterior distribution. In the Bayesian theory, the inverse Wishart distribution is usually used as the conjugate prior for the covariance matrix of a Gaussian distribution with known mean [36]. Since σ 2 r,k is the variance of the Gaussian distribution, the prior distribution p(σ 2 r,k |z im 1: k−1 ) is selected as the inverse Wishart distribution, given by p σ 2 r,k z im where IW(·; λ, Ψ) denotes the PDF of the inverse Wishart distribution with degree of freedom (dof) λ and inverse scale matrix Ψ [37],û k|k−1 andÛ k|k−1 are the dof and inverse scale matrix of p(σ 2 r,k |z im 1: k−1 ), respectively. The posterior distribution p(σ 2 r,k−1 |z im 1: k−1 ) is also subject to an inverse Wishart distribution as follows: To guarantee that p(σ 2 r,k−1 |z im 1: k−1 ) is the inverse Wishart distribution given by Equation 29, the previous approximate posterior distribution is spread through a forgetting factor ρ ∈ (0, 1], which indicates the extent of time fluctuations of the MNCM. Then, the prior dofû k|k−1 and prior inverse scale matrix U k|k−1 are given as follows [34]: where r denotes the order of the MNCM σ 2 r,k .

Variational approximations of posterior PDFs
According to the variational Bayesian approximation, the joint posterior PDF of the state of the i n -th target x in k and MNCM σ 2 r,k is approximated to where q(x in k ) and q(σ 2 r,k ) are the approximate posterior PDFs of x in k and σ 2 r,k , respectively [38,39]. The variational Bayesian approximation is formed by minimizing the Kullback-Leibler divergence (KLD) between the true joint distribution p(x in k , σ 2 r,k | z im 1: k ) and the approximate distribution q(x in k )q(σ 2 r,k ), i.e., where KLD(q(x) p(x)) denotes the KLD between q(x) and p(x) [38,39], and The optimal solution of Equation 34 satisfies the following equations [39]: where E x in k [·] and E σ 2 r,k [·] denote the expectation with regard to x in k and σ 2 r,k , respectively, and c x and c R denote the constants with respect to x in k and σ 2 r,k , respectively. Since the variational parameters of q(x in k ) and q(σ 2 r,k ) are coupled, a fixed point iteration process is applied to solve Equations 36, 37, i.e., the approximate posterior PDF q(x in k ) is updated to q (n+1) (x in k ) at the n + 1-th iteration using the posterior PDF q (n) (σ 2 r,k ), and q(σ 2 r,k ) is updated to q (n+1) (σ 2 r,k ) using the posterior q (n) (x in k ). According to Equations 25,26,29, the joint PDF is expressed as The posterior q (n+1) (x in k |z im 1: k−1 ) is updated according to the extended Kalman filter equations as where the mean vectorx (n+1) k|k and the covariance matrixP (n+1) k|k are given as follows: x n+1 (2) Update of σ 2 r,k According to Equation 38, log q (n) (σ 2 r,k ) is given by where the dofû (n+1) k and the inverse scale matrixÛ (n+1) k are given as follows:û Then, according to Equation 38, log q (n) (x in k ) is given by (49) The modified one-step predicted PDF p (n+1) (z k |x k ) at the n + 1-th iteration is defined as where the modified MNCM (σ (n+1) r,k ) 2 is formulated aŝ Finally, after N fixed-point iterations, the variational approximations of the posterior PDFs for the i n -th target are given as follows:

Algorithm of VB-GMCPHD filter
According to the above derivation, the pseudo-code of the variational Bayesian GM-CPHD filter (VB-GMCPHD) for DOA tracking at one time step is given in Algorithm 1. 3. Calculate the components of survive targets:

Initialize GM components
Add the components of birth Express the predicted GM components as Update: 6. Update the cardinality distribution p k (n) by using Equation 11; 7. Update GM components of targets: End End 9. J k|k J k|k−1 + J k|k−1 m k , Prune, merge, and limit the J k|k components, and new J k components are obtained; 10. Express the updated GM components as w (i) k , m (i) k , P (i)

Experimental setup
The open experimental data set SWellEx-96 [40] is used to verify the tracking performance of the proposed VB-GMCPHD filter for robust DOA tracking. The experiment was performed at the United States Marine Physical Laboratory from 10 to 18 May 1996, approximately 12 km from the tip of Point Loma near San Diego, California. The data of the north horizontal linear array of the S59 event of the experiment is used in this section. 900 s of experimental data is used to test the proposed method in this section. Before the experiment, a CTD is used to obtain the sound velocity in the experimental area. The experimental ship tows a continuous sound source at a depth of 54 m at a speed of 5 knots and sails north, which is named target 1 in this article. In addition, an uncooperative ship sails from northwest to southeast with continuous radiating noise, which is named target 2 in this article. The radar system of an experimental ship records the distance and bearing of the uncooperative ship and derives the latitude and longitude of the uncooperative ship. The horizontal linear array is placed on the seafloor at a depth of 213 m and continuously records the received acoustic signal at a sampling frequency of 3,276.8 Hz. The horizontal array consists of 32 elements, among which 27 elements provide effective data. The positions of the effective elements are shown in Figure 1, where the position of the first element is taken as the origin of coordinates. Figure 1 shows that the elements are actually not arranged in a straight line but in a slight curve. Therefore, the port and starboard ambiguity problem of a linear array is avoided.
To evaluate the difference between the real set and estimated set of the states of the targets, an evaluation has to be chosen. Because the number of elements of the real set and estimated set of the states of the targets is different, the root-mean-square type of evaluation (such as the root mean square error) is not usable. In this article, optimal sub-pattern assignment (OSPA) error is selected to evaluate the multi-target tracking performance, which is defined as follows [41]: .., y n are RFSs, Π denotes the set made up with m elements from 1, 2, ..., n { } , and d (c) (x, y) min d(x, y), c , d(x, y) is the Euclidean distance between x and y, while c and p denote the truncation parameter and order, and are set to 5 and 1, respectively. The smaller the OSPA error means the higher the precision.

Experimental results
The minimum variance distortionless response (MVDR) was used to obtain the measurement of the bearing angles of the targets at 48 Hz-52 Hz every 10 s, and the bearing angle scans from 0 to 360°e very 1°. According to the CTD data, the sound velocity was set to 1,493 m/s. The bearing time recording (BTR) obtained by the MVDR is given by the background 3D color diagram of Figure 2A. The true bearing angles of the two targets are given by the black line. The bearing angles corresponding to the peaks of the spectrum are extracted at each moment, and the obtained measurements of bearing angles of targets are shown by red dots in Figure 2A. Figure 2A shows trajectories of three targets. The bearing angle of target 1 moves from 135°to 50°, the bearing angle of target 2 moves from 285°to 275°, and the bearing angle of an unknown uncooperative target stays near 320°, which is named target 3 in this article. Figure 2A also shows some error, missing detection, and false alarm in the measurements.
The proposed robust multi-target DOA tracking method is used to process the bearing angle measurements, and the results of KF-JPDA [42], GM-PHD filter [27], and GM-CPHD filter [30] are also given for comparison. The process noise variance of the CV model is set to 2.5 × 10 −4 . The detection probability and false alarm probability are set to 0.9 and 0.1, respectively. The measurement noise variance σ 2 r,k of KF-JPDA, GM-PHD filter, and GM-CPHD filter is set to 25. The parameterŝ u 0 ,Û 0 , and ρ of the VB-GMCPHD are set to 12, 12 σ 2 r,k , and 0.95, respectively, and the number of VB iterations N is set to 5.
The tracking results of KF-JPDA, GM-PHD filter, GM-CPHD filter, and VB-GMCPHD filter are shown in Figures 2B-E, and the OSPA errors are shown in Figure 2F. The average OSPA errors tracking step are shown in Table 1; Figures 2B-E show that all multitarget DOA tracking methods carried out stable tracking of the three targets, and the tracking trajectories are consistent with the real trajectories. Figure 2F and Table 1 show that KF-JPDA, GM-PHD filter, GM-CPHD filter, and the proposed VB-GMCPHD filter significantly reduce OSPA errors on the basis of measurements, and the OSPA errors of the VB-GMCPHD filter are slightly smaller than those of the KF-JPDA, GM-PHD filter, and GM-CPHD filter.
Since the power of measurement noise hardly varies in the experiment, the performance of the proposed VB-GMCPHD filter for the robust multi-target DOA tracking improves slightly when compared to the other methods. In order to further test the robustness of the proposed robust underwater multi-target DOA tracking method in the scenario of uncertain measurement noise, a period of high-power noise was added to the experimental data from  Figure 3F, and the average OSPA errors per tracking step are shown in Table 2. Figure 3 shows that KF-JPDA, GM-PHD filter, and GM-CPHD filter provide stable tracking of the bearing angles of the three targets Frontiers in Physics frontiersin.org 09 when the measurement noise variance is stationary from 0 to 500 s. However, when the measurement noise variance increases from 500 s to 600 s, the assumption of fixed measurement noise variance of KF-JPDA, GM-PHD filter, and GM-CPHD filter becomes inconsistent with the real increasing measurement noise variance, which results in inaccurate Kalman filter gain. Therefore, the KF-JPDA, GM-PHD filter, and GM-CPHD filter carry out fluctuating tracking trajectory or even miss tracking, which results in an increasing of the OSPA error. The OSPA error of the proposed VB-GMCPHD filter for robust multitarget DOA tracking is significantly less than that of the KF-JPDA, GM-PHD filter, and GM-CPHD filter. The reason is that the VB-GMCPHD filter estimates the measurement noise variance in the real time of tracking so the Kalman filter gain calculated by using the estimate of measurement noise variance is more accurate. Therefore, the increasing measurement noise variance hardly affects the performance of the VB-GMCPHD filter, and the robust DOA tracking with uncertain measurement noise is carried out. Table 2 proves the superiority of the proposed VB-GMCPHD filter in robustness again.

Conclusion
The robust multi-target underwater DOA tracking problem in uncertain measurement noise scenario is studied in this article. In order to solve the underwater multi-target DOA tracking problem, a multitarget DOA tracking technique is derived on the basis of the GM-CPHD filter, which is combined the kinematics of the target, the bearing angle measurements, and the multi-target tracking scheme at the same time. Then, considering the uncertain measurement noise results from unknown underwater environment, the online measurement estimator is designed on the basis of the variational Bayesian approach to estimate the measurement noise variance along with the states of the targets as an integrated part during the multi-target DOA tracking procedure. Thus, a robust underwater multi-target DOA tracking method with uncertain measurement noise, named VB-GMCPHD filter, is proposed. From the experiment results validated by real sea trial data, the accuracy and robustness of the proposed VB-GMCPHD is verified and comprehensive discussions are made. From the experimental results and discussions, the proposed VB-GMCPHD can be regarded as an alternative method to accomplish DOA tracking missions, especially when the underwater environment is uncertain.

Data availability statement
Publicly available data sets were analyzed in this study. These data can be found at http://swellex96.ucsd.edu/.  Frontiers in Physics frontiersin.org 10