Tracking Multiple Targets Using Bearing-Only Measurements in Underwater Noisy Environments

This article handles tracking multiple targets using bearing-only measurements in underwater noisy environments. For tracking multiple targets in underwater noisy environments, the Gaussian Mixture Probability Hypothesis Density (GM-PHD) filter provides good performance with its low computational load. Bearing-only measurements are passive and do not provide position information of a target. Note that the nonlinearity of the bearing-only measurements can be handled by Extended Kalman Filters (EKF) when applying the GM-PHD filter. However, range uncertainty of the target is large for bearing-only measurements. Thus, a single EKF leads to poor performance when it is applied in the GM-PHD. In this article, every bearing measurement gives birth to multiple target samples, which are distributed considering the feasible range of the passive sensor. Thereafter, every target sample is updated utilizing the measurement update step of the EKF. In this way, we run multiple EKFs associated to multiple target samples, instead of running a single EKF. To the best of our knowledge, our article is novel in tracking multiple targets in noisy environments, using the observer with bearing-only measurements. The effectiveness of the proposed GM-PHD is verified utilizing MATLAB simulations.


Introduction
Multiple Target Tracking (MTT) is applied as multiple targets are tracked with a camera [1] or sonar sensors [2,3]. The goal of MTT is to conjecture the target information from a series of measurement sets.
In this article, we consider a single observer, which has a passive sonar (receiver) array for measuring the bearing of sound emitted from a target. Array signal processing algorithms, such as the MUltiple SIgnal Classification (MUSIC) algorithm, have been used to estimate the bearing of sound generated from a target [4][5][6][7][8].
For underwater MTT, passive sonars are desirable, since they work in a stealthy manner and consume little energy. However, background noise (interference) may generate bearing measurements that do not originate from a true target. In addition, the birth and death of the target result in uncertain data association between target and bearing measurement. Therefore, data association in noisy environments is required for tackling underwater multiple target tracking problems.
For processing bearing-only measurements, one can apply the Extended Kalman filter (EKF), which is the nonlinear version of the Kalman filter linearization about an estimate of the current mean and covariance [11,14,[18][19][20]. In the EKF, the state transition and observation models do not have to be linear functions of the state but may instead be differentiable functions. The EKF has both low computational complexity and reasonable estimation accuracy; therefore, it has been widely used for tracking based on bearing-only measurements [11,14,19,20]. In our paper, we also use the EKF for processing bearingonly measurements.
For handling an observer tracking a target using bearing-only measurements, the Range-Parameterized Extended Kalman Filter (RPEKF) was developed [19,20]. The RPEKF is a Gaussian-sum filter with multiple EKFs, each initialized at an estimated target range. In this way, the RPEKF reduces the initial range estimation error. The RPEKF assumes that the maximum sensing range of the observer is known in advance. Inspired by the RPEKF [19,20], our paper also assumes that the maximum sensing range of an observer is known in advance.
The RPEKF is suitable for tracking a single target using bearing-only measurements, and the RPEKF does not consider noisy environments. Distinct from the RPEKF, our paper considers tracking multiple targets using bearing-only measurements in noisy environments.
For tracking multiple targets in noisy environments, the Multiple Hypotheses Tracking (MHT) [21] and Joint Probability Data Association (JPDA) [22,23] filters were developed. However, the computational burden of these methods increases as the number of targets or interference points increases.
Recently, the Random Finite Set (RFS) was developed for avoiding explicit associations between measurement and target [3,24,25]. The Probability Hypothesis Density (PHD) filter [3,24,25] was invented using the RFS. The PHD is computationally effective, since it avoids data association between target and measurement. The authors of [25] invented the Gaussian Mixture Probability Hypothesis Density (GM-PHD) filter, which is time-efficient for solving the data association problem. For tracking multiple targets that maneuver in noisy environments, ref. [26] developed a multiple model GM-PHD filter. As far as we know, no paper on MTT has considered tracking multiple targets using bearing-only measurements.
In our article, the observer measures the bearing of sound emitted from a target in real time. Note that the nonlinearity of the bearing-only measurements can be handled by the EKF [18] when applying the GM-PHD filter. However, range uncertainty of the target is large for bearing-only measurements. This results in poor performance when a single EKF is utilized in the GM-PHD.
Therefore, in this article, every bearing measurement gives birth to multiple target samples, which are distributed considering the feasible range of the passive sensor. Note that the feasible range of the passive sensor is available since we assume that the maximum sensing range of an observer is known in advance. Thereafter, every target sample is updated using the measurement update step of the EKF. In this way, we run multiple EKFs associated to multiple target samples instead of running a single EKF.
As far as we know, this paper is novel in tracking multiple targets in noisy environments, as the observer uses bearing-only measurements. The effectiveness of the proposed GM-PHD is verified utilizing MATLAB simulations.
The remainder of this article is organized as follows. Section 2 introduces the bearingonly target tracking as the preliminary information of this article. Section 3 addresses the GM-PHD filter for bearing-only target tracking. Section 4 discusses MATLAB simulations for demonstrating the performance of the proposed filter. Conclusions are addressed in Section 5.

Preliminary Information
We introduce the matrix notations utilized in this paper. Let C[i] define the i-th element in a column vector C. Let I n define the identity matrix with n rows and n columns. Let 0 n define the zero matrix with n rows and n columns. diag(a, b, c, . . .) defines a diagonal matrix having a, b, c, . . . as its diagonal elements in this order.

Bearing-Only Target Tracking
Let x k define the target state at sampling-stamp k. The target's location and velocity at sampling-stamp k are defined as Suppose that the maximum speed of a target is known a priori. Let v max define the target's maximum speed.
Let T define the sampling interval. Let [a x , a y ] indicate the target's acceleration, which is not known in advance. Considering a target with Constant Velocity (CV) motion, its process model is where the state transition matrix F is In (2), µ k defines the process noise with mean 0 and variance Q, where Here, diag(σ 2 a x , σ 2 a y ) is a diagonal matrix having σ 2 a x , σ 2 a y as its diagonal elements in this order. Here, σ a presents the standard deviation for a. In addition, (4) utilizes Let [x o k , y o k ] define the 2D position of the observer at sampling-stamp k. In addition, [ẋ o k ,ẏ o k ] defines the 2D velocity of the observer at sampling-stamp k. In bearing-only target tracking, observer maneuvering is required for target localization [27,28]. This implies that the observer changes its velocity [ẋ o k ,ẏ o k ] in order to obtain the observability on the target position.
Let r max define the maximum sensing range of the observer's sonar. In addition, the sonar's minimum sensing range is r min . Inspired by the RPEKF [19,20], we assume that both r min and r max are known in advance.
Once the relative distance between the target and the observer is within the interval [0, r max ], the bearing of the target can be measured by the observer. The bearing measurement at sampling-stamp k is Here, atan2(y, x) is the phase (angle) of a complex number x + jy. In addition, e k defines the measurement noise having a Gaussian distribution with mean 0 and standard deviation R.
Ignoring the noise term in (6), we obtain

GM-PHD Filter
We utilize the GM-PHD filter in [25]. The multitarget state space is defined as F(X ). At sampling-stamp k, target sample set X k and bearing measurement set Z k are defined as Here, N(k) defines the number of target samples and M(k) defines the number of bearing measurements at sampling-stamp k.
Every x k ∈ X k defines the target pose (position and velocity) addressed in (8). For a given state X k−1 at sampling-stamp k − 1, every x k−1 ∈ X k−1 survives at sampling-stamp k with probability P S (x k−1 ). Therefore, the survived target is represented using a RFS S k|k−1 (x k−1 ).
At sampling-stamp k, a new target appears by spontaneous birth or by spawning from an existing target at sampling-stamp k − 1. In addition, spontaneous birth sets Γ k and spawned target sets B k|k−1 (x k−1 ) are modeled as RFS at sampling-stamp k. Birth and spawn of a new target from bearing measurements are introduced in Section 3.
At sampling-stamp k, X k is In addition, at sampling-stamp k, a target is detected by bearing sensors with probability P D (x k ). Every target x k ∈ X k generates its associated bearing measurement b(x k ) at sampling-stamp k. See (6) for the equation of b(x k ). In addition, the bearing sensor measures an interference point at sampling-stamp k. An interference point is modeled as a RFS K k . Therefore, the RFS measurement set Z k is Let p k (X k |Z 1:k ) define the posterior density, f k|k−1 define the transition density, and g k (Z k |X) define the likelihood of every target. The posterior is calculated as Here, µ is an appropriate reference measure on F(X ) [29]. It is assumed that an interference is independent of target-originated measurements. In addition, it is assumed that an interference point and predicted RFS are distributed using a Poisson distribution.
The survival and detection probabilities do not depend on the state vector; hence, we utilize P S,k and P D,k instead of P S (x k ) and P D (x k ), respectively. Let N (x; m, P) denote the Gaussian density with mean m and covariance P. The birth intensity function of a newborn target at sampling-stamp k is where J γ,k represents the number of newborn Gaussian components, and w Here, λ k is the average number of interference points at sampling-stamp k. In addition, c k (z) defines the probability distribution of every interference point.
In the proposed GM-PHD, the prediction step is linear and the measurement update step is nonlinear. Thus, we only have to linearize the measurement update step. In our paper, the EKF is used for the measurement update step, since the EKF has both low computational complexity and reasonable estimation accuracy [11,14,19,20].
We perform the measurement update, once we obtain the newest measurement set Z k at sampling-stamp k. Since the measurement model in (6) is not linear, we can utilize the EKF in [18] for the measurement update step.
Under the GM-PHD [25,30], (15) is replaced by where Here, J k−1 is the number of components of the intensity at sampling-stamp k − 1. Moreover, we utilize In addition, we utilize Further, (16) is replaced by In (21), we use Here, one uses where In (21), we use Here, K is In addition, k|k−1 . In (21), we further use where I is the identity matrix. The detailed derivation of (21) is addressed in [25,30].
To reduce the computational load, Gaussian components must be pruned and merged [30]. Then, the posterior intensity is represented as Here, w (i) k is the weight of Gaussian distribution and J k is the number of components of the intensity at sampling-stamp k.

Birth or Spawn of New Target Samples Based on Bearing Measurements
Recall that at every sampling-stamp k, a new target sample is generated by spontaneous birth or by spawning from an existing target sample at sampling-stamp k − 1. Spontaneous birth sets Γ k and spawned target sets B k|k−1 (x k−1 ) are modeled as RFS at sampling-stamp k.
In this article, the observer measures the bearing of sound emitted from a target in real time. Bearing-only measurements are passive and do not provide position information of a target. Therefore, in this article, every bearing measurement gives birth to multiple target samples, which are distributed considering the feasible range of the passive sensor.
We update the mean m and covariance P of the newly born targets in Γ k as follows. Suppose that we generate B target samples using z k,g (g ∈ {1, 2, . . . , M(k)}), the g-th bearing measurement at sampling-stamp k. Let t b (b ∈ {1, 2, . . . , B}) define the b-th target sample generated on the bearing line at sampling-stamp k.
The feasible range of the true target is in the interval [r min , The range for the center of every subinterval is r b = r min + (b − 0.5) * sR (b ∈ {1, 2, . . . , B}). In addition, the velocity of every newly born target sample is zero. Under (7), the mean m of a newly born target sample where r b = r min + (b − 0.5) * sR.
We next derive the covariance P of t b . Under (6), the variance of the bearing error is R 2 . In addition, since every subinterval length is sR, the variance of the range error is sR 2 4 .
The error covariance related to the 2D position of t b (b ∈ {1, 2, . . . , B}) is where D is Under (30), we further obtain The covariance P of t b (b ∈ {1, 2, . . . , B}) is derived as In the case where the relative distance between a newly born target sample, e.g., t b (b ∈ {1, 2, . . . , B}), and an existing target sample, e.g., x k−1 , is less than a certain threshold α > 0, we say that the born target sample t b is spawned from the existing target sample x k−1 . This implies that we have t b ∈ B k|k−1 (x k−1 ).
We update the mean m and covariance P of the newly spawned target sample t b ∈ B k|k−1 (x k−1 ) as follows. We derive the covariance P of t b using (33). In addition, we use

MATLAB Simulations
We address MATLAB simulation results for verifying the effectiveness of the proposed GM-PHD filter. Target bearings are measured at every dt = 1 s. In addition, the simulation running time is 500 s.
Target detection probability is 0.98. The target's maximum speed is v max = 30 m/s. Considering the sonar sensing range, we set r max = 10, 000 m. In addition, the sonar's minimum sensing range is set to r min = 500 m. The standard deviation of bearing measurement noise is 1 degree.
In bearing-only target tracking, observer maneuvering is required for improving the tracking convergence [27,28]. Since the observer uses bearing-only measurements, the tracking convergence depends on the bearing rate of a target with respect to the observer. In other words, the geometry of a target with respect to the observer determines the tracking convergence of the target.
In the simulations, the observer maneuvers as follows. The observer moves with a constant speed v obs = 10 m/s, and it starts from the origin. Before 125 s, the observer moves with course 0 degree measured clockwise from the east direction. From 125 to 250 s, the observer moves with course 135 degrees measured clockwise from the east direction. From 250 to 375 s, the observer moves with course 0 degree measured clockwise from the east direction. After 375 s, the observer moves with course 90 degrees measured clockwise from the east direction.
See Figure 1 for the simulation scenario. In this figure, the observer's path is marked with yellow colors. Every target moves with a constant velocity. In Figure 1, we present the direction of target motion using symbols. A circle symbol indicates the start point of a target motion and a diamond symbol indicates the end point of a target motion. Since the sonar's maximum sensing range is 10 km, we deploy a target initially so that it can be detected by the observer. The position of each target at every second is marked with distinct colors. Target  For clear presentation of the observer's path, Figure 1 simulates the case where there is no interference point. In other words, Figure 1 shows the simulation result without interference points. In this figure, the estimated target position is marked with magenta triangles. Before the observer maneuver at 125 s, true target positions are not tracked by target estimations (magenta triangles). However, as time goes on, true target positions are tracked by target estimations.
Considering the simulation of the scenario in Figure 1, Figure 2 shows the simulation result with interference points. Recall that λ k (16) is the average number of interference points at sampling-stamp k. We set λ k = 1, which implies that, on average, one interference point is generated at each sampling-stamp. In Figure 2, the estimated target position is marked with magenta triangles. Due to interference points, many false targets are generated surrounding the observer. Before the observer maneuver at 125 s, true target positions are not tracked by target estimations (magenta triangles). However, as time goes on, true target positions are tracked by target estimations.
Considering the simulation of Figure 2, Figure 3 shows the bearing (in radians) change with respect to time (in seconds). In this figure, true target bearing measurements are marked with circles. In addition, interference points are marked with crosses. The estimated target bearing is marked with magenta circles. Before the observer maneuver at 125 s, true target bearings are not tracked by the estimated target bearing (magenta circles). However, as time goes on, true target bearings are tracked by estimated target bearings.   The estimated target bearing is marked with magenta circles. Before the observer maneuver at 125 s, true target bearings are not tracked by the estimated target bearing (magenta circles). However, as time goes on, true target bearings are tracked by estimated target bearings.

Monte-Carlo (MC) Simulations
A total of 100 Monte-Carlo (MC) simulations are utilized for verifying the outperformance of the proposed filter. Let M C = 100 denote the number of total MC simulations. We set λ k = 1, which implies that, on average, one interference point is generated at each sampling-stamp. We generate B = 50 target samples using every bearing measurement at every sampling-stamp. It takes 14 s to run a single MC simulation.
Tracking accuracy is checked with the cardinality (number of targets) and the Optimal Subpattern Assignment (OSPA) metric [31,32]. We derive the cardinality as the mean of the estimated number of targets under M C = 100 MC simulations. We also derive the mean of OSPA using M C = 100 MC simulations.
Recall that J k in (28) is the number of components of the intensity. The cardinality at sampling-stamp k is defined as The true cardinality at sampling-stamp k indicates the true number of targets detected at sampling-stamp k.
Let M(k) = {m 1 (k), m 2 (k), . . . , m J k (k)} denote the set of target position estimates at sample-step k. Recall that J k denotes the number of components of the intensity. Suppose that there are N(k) targets at sample-step k. Let Tgt(k) = {tgt 1 (k), tgt 2 (k), . . . , tgt N(k) (k)} denote the true target positions at sample-step k. The OSPA distance between the target estimates M(k) and the true target positions Tgt(k) is defined as Here, P N(k) denotes the set of permutations on {1, 2, . . . , N(k)}. The OSPA metric can be computed efficiently by using the Hungarian method for optimal point assignment [31,32].
The mean of OSPA at sample-step k is See Figure 4 for the OSPA ospa k and the cardinality card k of the proposed GM-PHD. Before 30 s pass, there are two targets (Targets 1 and 2). Target 3 is spawned from target 1 at 30 s. In bearing-only target tracking, observer maneuver is required for target localization [27,28]. The observer maneuvers at 125 s. See that the OSPA decreases gradually after the observer maneuver. In Figure 4, True plots the number of bearing measurements associated with true targets. In addition, Estimate is the mean of the estimated number of targets using 100 MC simulations. As time goes on, OSPA gradually approaches to zero and Estimate converges to True in cardinality. In other words, the cardinality error converges to almost zero as time goes on.  Note that the nonlinearity of the bearing-only measurements can be handled by the EKF when applying the GM-PHD filter. However, range uncertainty of the target is large for bearing-only measurements (r max = 10 km in MATLAB simulations). This results in poor performance when a single EKF is used in the GM-PHD. Setting B = 1 corresponds to the case where we apply a single EKF in the GM-PHD. As we increase B, we distribute more target samples considering the maximum sensing range of the bearing sensor. Consider the scenario in Figure 2. We next check the effect of changing B, the number of target samples generated using every bearing measurement at every sampling-stamp.
Considering the case where B = 1 is utilized, Figure 5 plots the OSPA and the cardinality of the proposed GM-PHD. Compared with the case where B = 50 is used (Figure 4), the estimation accuracy decreases considerably. This shows the effectiveness of the proposed sample distribution approach used in our paper. It takes 6 s to run a single MC simulation, as we utilize B = 1. Compared with the case where B = 50 is used (Figure 4), the computational load decreases.

The Effect of Changing B
Considering the case where B = 10 is utilized, Figure 6 plots the OSPA ospa k and the cardinality card k of the proposed GM-PHD. Compared with the case where B = 1 is used (Figure 5), the estimation accuracy increases considerably. It takes 6 s to run a single MC simulation. Compared with the case where B = 50 is used, the computational load decreases.
Considering the case where B = 20 is utilized, Figure 7 plots the OSPA and the cardinality of the proposed GM-PHD. Compared with the case where B = 1 is used ( Figure 5), the estimation accuracy increases considerably. It takes 7 s to run a single MC simulation. Compared with the case where B = 50 is used, the computational load decreases.  Table 1 shows the computational load (running time for one MC simulation) as B varies. Recall that setting B = 1 corresponds to the case where we apply a single EKF in the GM-PHD. Table 1 shows that as B increases, the computational load increases in general. However, considering both the computational load and the estimation accuracy ( Figure 6), setting B = 10 is desirable.

Conclusions
This article handles tracking multiple targets using bearing-only measurements in underwater noisy environments. We apply GM-PHD filters to solve our MTT problem based on bearing-only measurements. Since bearing-only measurements do not provide position information of a target, we make every bearing measurement give birth to B target samples, which are distributed considering the feasible range of the passive sensor. The effectiveness of the proposed GM-PHD is verified utilizing MATLAB simulations.
In the future, we will extend the proposed tracking filter so that one can handle the case where multiple observers are used for tracking multiple targets in noisy environments under bearing-only measurements. As we use multiple observers, false targets may appear as bearing lines of multiple observers' intersects. Thus, we require filters for removing false targets effectively.