Multi-Target Tracking in Multi-Static Networks with Autonomous Underwater Vehicles Using a Robust Multi-Sensor Labeled Multi-Bernoulli Filter

Abstract

This paper proposes a centralized MTT method based on a state-of-the-art multi-sensor labeled multi-Bernoulli (LMB) filter in underwater multi-static networks with autonomous underwater vehicles (AUVs). The LMB filter can accurately extract the number of targets and trajectories from measurements affected by noise, missed detections, false alarms and port–starboard ambiguity. However, its complexity increases as the number of sensors increases. In addition, due to the time-varying underwater environment, AUV detection probabilities are time-varying, and their mismatches often lead to poor MTT performance. Consequently, we detail a robust multi-sensor LMB filter that estimates detection probabilities and multi-target states simultaneously in real time. Moreover, we derive an effective approximate form of the multi-sensor LMB filter using Kullback–Leibler divergence and develop an efficient belief propagation (BP) implementation of the multi-sensor LMB filter. Our method scales linearly with the number of AUVs, providing good scalability and low computational complexity. The proposed method demonstrates superior performance in underwater multi-AUV network MTT simulations.

1. Introduction

As a crucial aspect of ocean surveillance, multi-target tracking (MTT) aims to accurately estimate the states and number of targets by processing the collected measurements [1,2,3,4]. While passive MTT using passive sonar has been the preferred method for underwater MTT for many years [5,6], the increased stealth capabilities of targets and higher ambient noise levels have given a renewed impetus to active sonar target tracking [7,8,9]. Recent advances in marine technology have shown that multi-static sonar networks, which involve multiple entities transmitting signals and receiving echoes reflected from targets, have a significant advantage in large-range underwater MTT [10,11].

Technological developments have made small, relatively low-cost underwater intelligent systems, especially autonomous underwater vehicles (AUVs), a reality [12,13,14]. With their mobility and ability to reconfigure locations, AUVs have gained considerable attention as “listeners” in multi-static networks [15,16]. As shown in Figure 1, AUVs passively receive echoes while spatially separated from transmitters, ensuring stealth and endurance. In addition, multiple AUVs with multiple sources constituting different source/receiver (S/R) pairs can produce different source–target–receiver configurations, improving the MTT performance and potentially increasing the coverage significantly.

However, MTT can be complicated by noise, false alarms, missed detections, and potential errors in associating measurements with specific targets. The classic probabilistic methods for MTT are Multiple Hypothesis Tracking (MHT) [17] and Joint Probability Data Association (JPDA) [18]. MHT evaluates the posterior probability of all feasible data association hypotheses. It maintains a hypothesis tree for each track, selects the most probable association hypotheses, and returns a track estimate conditioned on the hypothesis. JPDA deals with the situation in which measurements fall within the intersection region of multiple tracks. It calculates the association probability between each measurement and every target, and tries to marginalize the association variables to compute the marginal distribution for each track. To provide a tool for the false track identification process, joint integrated probabilistic data association (JIPDA) [19,20] and joint integrated track splitting (JITS) [20,21] include the target existence as an additional random variable in estimation. Unlike probabilistic methods, the later developed random finite set (RFS)-based methods [22,23] adopt several fundamental concepts of Bayesian estimation principles and provide a top-down mechanism to achieve MTT. The RFS-based methods mainly include the probability hypothesis density (PHD) [24], cardinality PHD (CPHD) [25], multi-Bernoulli (MB) [26], generalized labeled multi-Bernoulli (GLMB) [27], and labeled multi-Bernoulli (LMB) [28]. As with JIPDA and JITS, the labeled (unlabeled) MB filters use the target existence probability to manage target tracks. However, they differ in that they provide the basis for Bayesian updates by describing the existence or non-existence of targets through the Bernoulli process in the RFS framework. It is also worth noting that the GLMB and LMB filters are particularly advanced as they can introduce labels for targets, making it easier to identify target trajectories. Therefore, this paper mainly investigates the underwater MTT within the RFS framework.

In multi-static systems, multiple S/R pairs can be viewed as numerous independent sensors, making MTT a multi-sensor MTT problem. Multi-sensor MTT based on RFS can be classified into two main categories: data-level methods and estimate-level methods [29]. Data-level methods are centralized and calculate the joint likelihood of all sensor measurements. This approach is optimal for multi-sensor MTT, but presents the challenge of the NP-hard multidimensional assignment problem. To address this issue, various centralized multi-sensor MTT filters have been proposed, including multi-sensor PHD/CPHD [30], multi-sensor multi-Bernoulli (MS-MB) [31], and multi-sensor GLMB filters [32]. In contrast, estimate-level multi-sensor MTT methods usually use geometric or linear arithmetic averaging to fuse the MTT densities obtained locally by each sensor [29]. Although estimate-level strategies for multi-sensor MTT reduce the computational complexity to some extent and can facilitate distributed MTT, it can be challenging to establish a relationship between the local results and the fused MTT results.

Additionally, the probability of detecting a target in a specific state has a substantial impact on MTT performance in real underwater MTT. Detection probability mismatches can negatively affect the tracking accuracy, resulting in target loss or misidentification errors. Many MTT algorithms rely on fixed detection probabilities based on experience or training data. However, predetermined parameters are insufficient in a changing underwater environment. As a result, researchers have proposed to predict the detection probability using mathematical acoustic propagation models [33,34,35]. In this way, the detection probability corresponds to a desired false alarm density with a signal-to-noise ratio (SNR), which should consider the source parameters, receiver parameters, source–target–receiver geometry, and environments such as multipath propagation and reverberation. Thus, the use of a high-fidelity acoustic propagation model to forecast detection probabilities can be time-consuming and requires more accurate ocean environment and target information.

Accordingly, several RFS filters have been suggested to address the detection probability mismatch. These filters enable the online adaptive estimation of the target detection probabilities. In previous publications [36,37,38], this issue has been discussed for single-sensor systems. In the multi-sensor case, the inverse gamma Gaussian mixture approach was employed in the MS-MB and MS-CPHD filters to estimate target detection features [39]. The hidden Markov process also handles the time-varying detection probabilities [40,41]. Although these filters provide reliable estimates of current target states, they do not produce target tracks. Therefore, a multi-sensor MTT filter that bootstraps the detection probability estimated from the CPHD filter into the GLMB filter to achieve MTT has been proposed [42,43]. Nevertheless, the complexity of the above filters increases exponentially as the sensor number grows; this is a crucial consideration when designing real-world MTT systems.

Here, we propose a centralized MTT method using a robust multi-sensor LMB (R-MS-LMB) filter, which approximates the more sophisticated GLMB filter and can learn and estimate detection probabilities while performing MTT. Notably, for the data association of multiple S/R pairs, we employ an efficient belief propagation (BP) implementation of a multi-sensor LMB filter, which scales linearly and performs better. Similar methods have been successfully applied to approximate the LMB filter [44,45] and estimate detection probabilities. In contrast to these previous works, we derive and present an efficient method of implementing the R-MS-LMB filter using the BP scheme within the RFS framework. The following are the most significant contributions of the paper.

(1) A robust multi-sensor LMB model is applied to account for time-varying detection probabilities. The probabilities of a target being detected by different S/R pairs are integrated into each target state and adjusted online as unknown variables rather than known parameters.

(2) Following the update step, as the posterior distribution is of a GLMB form [28], we present a novel derivation of LMB approximation through Kullback–Leibler divergence (KLD) minimization. As a result, we can reconstruct the GLMB based on the distribution of the association variables. We approximate the LMB by marginalizing the association variables. This novel derivation serves as a foundation for efficiently implementing a multi-sensor LMB filter.

(3) We construct a BP-based framework on a factor graph to implement the R-MS-LMB filter efficiently. The inner BP scheme calculates the marginal probability distribution of data association variables. The outer BP scheme is employed at each step to calculate the marginal posterior of the target states and the detection probabilities of the targets. With the BP-based LMB filter implementation, we can avoid pruning GLMB components.

The remainder of this paper is presented as follows. Section 2 introduces the background and objective. Section 3 describes the detection probability model and derives the novel LMB filter form. We present a BP-based implementation of an adaptive multi-sensor LMB filter in Section 4. We show the simulation results in Section 5 and conclude the work in Section 6.

2. Background and Objective

This paper considers a mobile multi-static network for MTT, as shown in Figure 1. The network comprises $N^R$ AUVs that tow sonar arrays as receivers and $N^S$ sources that act as transmitters. The sets of receivers and sources are denoted as $\mathcal{O}=\{1,\dots ,N^O\}$ and $\mathcal{S}=\{1,\dots ,N^S\}$, respectively. AUVs and different sources can be combined into separate S/R pairs. Therefore, we define the S/R pair index set as $\mathcal{V}\triangleq \{1,\dots ,\mathrm{V}\}$ with $\mathrm{V}=N^R{N}^O$, and introduce an indexing function $\pi :\mathcal{V}\to \mathcal{S}\times \mathcal{O}$, where $\pi (v)=(o,s)$ denotes the S/R pair v composed of AUV o and source s. Furthermore, we represent the state of AUV o at time t using a vector $y_t^o={[y_{t,x}^o,y_{t,y}^o,{\psi }_t^o]}^T$, where ${[y_{t,x}^o,y_{t,y}^o]}^T$ and ${\psi }_t^o$ are the two-dimensional horizontal position and heading, respectively. The position of source s is defined as $q_t^s={[q_{t,x}^s,q_{t,y}^s]}^T$.

2.1. Multi-Target State and LMB RFS

It is evident that the state and number of targets change over time; therefore, the conventional use of vectors to represent multi-target states is inadequate. An alternative approach is to use random finite sets (RFSs) to describe the multi-target states, where the number of elements and their order are random. As elaborated in [46], an RFS representation is generally more suitable than a vector representation for MTT. Readers are referred to [22,23] or more mathematical concepts and definitions related to RFS. In the following, the single-target state is denoted by lower-case letters (e.g., $x$, $x$) while multi-target states are indicated by upper-case letters (e.g., $X$, $\mathit{X}$). Bold-face symbols are used for labeled versions (e.g., $x$, $\mathit{X}$) to differentiate them from unlabeled versions (e.g., $x$, $X$), and blackboard bold (e.g., $\mathbb{X}$, $\mathbb{L}$, $\mathbb{Z}$, etc.) is used to represent the spaces. Moreover, $ℱ(\cdot )$ represents the set of all finite subsets of a space.

Following [27], the number of elements of an RFS X is represented by |X|, the inner product is represented as c, and the generalized Kronecker delta function is represented by ${\mathsf{\delta }}_Y(X)\triangleq \{\begin{array}{cc}1,& X=Y\\ 0,& \mathrm{otherwise}.\end{array}$. As a general rule, the indicator function is given by $1_Y(X)\triangleq \{\begin{array}{cc}1,& \mathrm{if} X\subseteq Y\\ 0,& \mathrm{otherwise}.\end{array}$.

Let the vector $x_t=[p_{t,x}, p_{t,y}, {\dot{p}}_{t,x}, {\dot{p}}_{t,y}]\in \mathbb{X}$ represent the single-target state, where $[p_{t,x}, p_{t,y}]$ and $\left[{\dot{p}}_{t,x}, {\dot{p}}_{t,y}\right]$ denote the locations and velocities, respectively. The multi-target states are defined as an RFS $X_t=\{x_{t,1},\dots ,x_{t,n}\}\in ℱ(\mathbb{X})$ at time t. A single-target state $x_t$ is augmented with a unique label $\ell \in {\mathbb{L}}_t$, i.e., $x_t=(x_t,\ell )\in \mathbb{X}\times {\mathbb{L}}_t$ abbreviated as $x_t^{\ell }$ for distinguishing target identities. Up to time t, all the targets have been labeled as the union of disjoint sets ${\mathbb{L}}_t={\mathbb{L}}_{t-1}\cup {\mathbb{L}}_t^B$, including the label space ${\mathbb{L}}_{t-1}$ prior to time t and the birth targets’ label space at time t. The corresponding labeled multi-target state can be represented by $\mathit{X}=\{(x,{\ell }_1),\dots (x,{\ell }_n)\}$, and the set of labels of $\mathit{X}$ is given by $ℒ(\mathit{X})=\{ℒ(\mathit{x}):\mathit{x}\in \mathit{X}\}$, where $ℒ: \mathbb{X}\times \mathbb{L}\to \mathbb{L}$ is a projection with $ℒ(x,\ell )=\ell$.

Considering that a previously detected potential target (PT) may exist or disappear, a detected PT is represented by a Bernoulli process incorporating a state probability density function (pdf) $p(x)$ and an existence probability r. Existence probability r represents the estimated probability of the presence of the detected PT. The state pdf $p(x)$ denotes the probability statistic of the target state conditional on the presence of the PT. Typically, a single labeled Bernoulli RFS ${\mathit{X}}^{\ell }$ corresponding to a Bernoulli component with a distinct label $\ell \in \mathbb{L}$ is specified by $\{r^{\ell },p^{\ell }(x)\}$, where $r^{\ell }$ describes its existence probability and $p^{\ell }(x)=p(x,\ell )$ is its respective state pdf. The single labeled Bernoulli RFS ${\mathit{X}}^{\ell }$ with parameter $\{r^{\ell },p^{\ell }(x)\}$ has the following density:

$$f^{\ell }({\mathit{X}}^{\ell })=\{\begin{array}{ll}1-r^{\ell }\hfill & {\mathit{X}}^{\ell }=\varnothing \hfill \\ r^{\ell }{p}^{\ell }(x)\hfill & {\mathit{X}}^{\ell }=\{(x,\ell )\}\hfill \\ 0\hfill & \mathrm{otherwise} \hfill \end{array}.$$

(1)

where ${\mathit{X}}^{\ell }=\varnothing$ implies that the PT $\ell \in \mathbb{L}$ does not exist with probability $r^{\ell }$, and ${\mathit{X}}^{\ell }=\{(x,\ell )\}$ means that the PT $\ell \in \mathbb{L}$ exists with probability $r^{\ell }$ and the state pdf is $p^{\ell }(x)$.

An LMB realization with the parameter ${\{(r^{\ell },p^{\ell }(x))\}}_{\ell \in \mathbb{L}}$ is a set $\mathit{X}={\cup }_{i=1}^n{\mathit{X}}^{{\ell }_i}=\{(x,{\ell }_1),\dots (x,{\ell }_n)\}$ given by a union of independent labeled Bernoulli RFSs ${\mathit{X}}^{{\ell }_i}=\{(x,{\ell }_i)\}$, ${\ell }_i={\ell }_1,\dots ,{\ell }_n$. The density of an LMB $\mathit{X}$ can be expressed as follows [28]:

$$f(X)=\Delta X{\displaystyle \prod _{{\ell }^{\prime }\in \mathbb{L}\backslash ℒ(X)}(1-r^{{\ell }^{\prime }})}{\displaystyle \prod _{\ell \in ℒ(X)}{1}_{\mathbb{L}}(\ell )r^{\ell }{p}^{\ell }(x)}.$$

(2)

where $\Delta \left(\mathit{X}\right)={\delta }_{\left|\mathit{X}\right|}(ℒ(\mathit{X})$ is one if $\mathit{X}$ has distinct labels; otherwise, it is zero.

2.2. Measurement Model

This study assumes frequency-modulated (FM) pulses in pulsed active sonar modes for their simplicity and excellent range resolution. FM pulses offer the advantage of being orthogonal in a specific domain [47], which allows us to ignore multi-source signal interference. Each autonomous underwater vehicle (AUV) can independently determine the source of the signal, causing multiple source/receiver (S/R) pairs to be treated as independent sensors. For the S/R pair v with $\pi (v)=(o,s)$, the signal emitted by source s may reach AUV o after being reflected by a PT. In the bistatic setup shown in Figure 2, AUVs equipped with rigid arrays have the same heading as the array, so the AUV o can extract a bearing-range measurement [48]:

$$z_t=\left[\begin{array}{l}\varphi \\ d\end{array}\right]=f(x_t,y_t^o,q_t^s)+{\epsilon }_t$$

(3)

where $f(x_t,y_t^o,q_t^s)=\left[\begin{array}{c}\arctan \left(p_{t,x}-y_{t,x}^a,p_{t,y}-y_{t,y}^a\right)-{\psi }_t^a\\ \sqrt{{\left(p_{t,x}-y_{t,x}^a\right)}^2+{\left(p_{t,y}-y_{t,y}^a\right)}^2}+\sqrt{{\left(p_{t,x}-q_{t,x}^s\right)}^2+{\left(p_{t,y}-q_{t,y}^s\right)}^2}\end{array}\right]$, and ${\epsilon }_t\sim \mathcal{N}\left(\left[\begin{array}{l}0\\ 0\end{array}\right],\left[\begin{array}{cc}{\sigma }_d^2& 0\\ 0& {\sigma }_{\varphi }^2\end{array}\right]\right)$, where ${\sigma }_{\varphi }$ and ${\sigma }_d$ are the additive noise to the bearing and range, respectively.

Due to the port–starboard ambiguity (PSA) problem caused by the AUV carrying a linear hydrophone array, as shown in Figure 2, each measurement recorded by the AUV with a different bearing but in the same range will result in a “ghost” measurement [49]. In this case, we assume that the AUV only stores the port measurement, while the respective starboard measurement is obtained by ${\iota }_v(z_t)={[-\varphi ,r]}^T$. Given that each S/R pair outputs at most one measurement for each PT, the single-target-originated likelihood of the v-th S/R pair with $\pi (v)=(o,s)$ is expressed as follows:

$$l_v(z|x)=\{\begin{array}{cc}{h}_v(z|x)& \mathrm{if} x \mathrm{on} \mathrm{the} \mathrm{port}\\ h_v({\iota }_v(z)|x)& \mathrm{if} x \mathrm{on} \mathrm{the} \mathrm{starb}\end{array}$$

(4)

where $h_v(z|x)=\mathcal{N}\left(f(x_t,y_t^o,q_t^s),\left[\begin{array}{cc}{\sigma }_d^2& 0\\ 0& {\sigma }_{\varphi }^2\end{array}\right]\right)$.

Under the S/R pair $v\in \mathcal{V}$, in the case of a multi-target state ${\mathit{X}}_t$, each state ${\mathit{x}}_t=(x_t,\ell )\in {\mathit{X}}_t$ is either detected with probability $P_{D_v}(x_t)$ and generates a measurement $z_t$ with likelihood $l_v(z_t|x_t)$ or is missed with a probability of $1-P_{D_v}(x_t)$. False alarm measurements (clutter) are modeled by a Poisson RFS with a clutter rate of ${\lambda }_v$ and probability density function (pdf) of $c(z)$. For a combination of target-originated measurements and clutter, the v-th S/R pair output $M_{t,v}$ measurements at time t are denoted as an RFS $Z_{t,v}=$$\{z_{t,v}^1,\dots ,z_{t,v}^{M_{t,v}}\}\in ℱ(\mathbb{Z})$. We denote ${ℳ}_{t,v}\triangleq \{1,\dots ,M_{t,v}\}$ by the measurement index set of S/R pair v.

The standard multi-target likelihood function for the S/R pair v is given by [27]

$$g_v(Z_{t,v}|X_t)\propto {\displaystyle \sum _{{\theta }_{t,v}\in {\Theta }_{t,v}}{1}_{{\Theta }_{t,v}(ℒ(X))}({\theta }_{t,v})}{\displaystyle {\prod }_{(x_t,\ell )\in X_t}{\kappa }_v^{{\theta }_{t,v}(\ell )}(x_t,\ell )}$$

(5)

where ${\kappa }_v^{{\theta }_{t,v}(\ell )}(x_t,\ell )=$$\{\begin{array}{cc}\frac{P_{D_v}(x_t,\ell )l_v(z_{t,v}^{{\theta }_{t,v}(\ell )}|x_t)}{{\lambda }_{v}c(z_t^{{\theta }_{t,v}(\ell )})}& if {\theta }_{t,v}(\ell )\in {ℳ}_{t,v}\\ 1-P_{D_v}(x_t,\ell )& if {\theta }_{t,v}(\ell )=0\end{array}$, and ${\Theta }_{t,v}$ indicates the space of admissible target–measurement association maps ${\theta }_{t,v}(\ell ):{\mathbb{L}}_t\to \{0,1,\dots ,M_{t,v}\}$, where ${\theta }_{t,v}(\ell )=m\in {ℳ}_{t,v}$ specifies that target $\ell$ generates $z_t^m\in Z_{t,v}$, and any misdetection is assigned a 0. The admissible association map means that it assigns each measurement to at most one target and each target to at most one measurement. ${\Theta }_{t,v}(ℒ(X))$ is the subset of ${\Theta }_{t,v}$ with $ℒ(X)\in {\mathbb{L}}_t$.

Let us define the set of measurements for all S/R pairs at time t as $Z_t\triangleq (Z_{t,1},\dots ,Z_{t,\mathrm{V}})$, and the set of measurements for all times as $Z_{1:t}\triangleq (Z_1,\dots ,Z_t)$. Based on the assumption that all S/R pairs are conditionally independent, we can express the multi-target likelihood of multiple S/R pairs as follows [32]:

$$g(Z_t|X_t)={\displaystyle {\prod }_{v=1}^{\mathrm{V}}{g}_v(Z_{t,v}|X_t)}={\displaystyle {\sum }_{{\theta }_t\in {\Theta }_t}{1}_{{\Theta }_t(ℒ(X))}({\theta }_t)}{\displaystyle {\prod }_{(x_t,\ell )\in X_t}{\kappa }^{{\theta }_t(\ell )}(x_t,\ell )}$$

(6)

Here, ${\theta }_t\triangleq ({\theta }_{t,1},\dots ,{\theta }_{t,V})$ represents the multi-sensor association map, and ${\Theta }_t\triangleq {\Theta }_{t,1}\times ,\cdots ,\times {\Theta }_{t,\mathrm{V}}$ refers to the association map space. The function ${\kappa }^{{\theta }_{t,v}(\ell )}(x_t,\ell )$ is given by ${\kappa }^{{\theta }_t(\ell )}(x_t,\ell )={\displaystyle {\prod }_{v=1}^{\mathrm{V}}{\kappa }_v^{{\theta }_{t,v}(\ell )}(x_t,\ell )}$, and the indicator function $1_{{\Theta }_t(ℒ(X))}({\theta }_t)$ indicates $1_{{\Theta }_t(ℒ(X))}({\theta }_t)={\displaystyle {\prod }_{v=1}^{\mathrm{V}}{1}_{{\Theta }_{t,\mathrm{V}}(ℒ(X))}({\theta }_{t,v})}$.

2.3. Multi-Sensor LMB Filter

As shown in [27,28], given that the multi-target posterior $f({\mathit{X}}_{t-1}|Z_{1:t-1})$ at time t-1 is the LMB form with the parameters ${\left\{(r_{t-1}^{\ell },p_{t-1}^{\ell }(x_t))\right\}}_{\ell \in {\mathbb{L}}_{t-1}}$ and the newbirthed PTs have a pdf with the parameters ${\left\{\left(r_{t|t-1}^{\ell },p_{t|t-1}^{\ell }(x_t)\right)\right\}}_{\ell \in {\mathbb{L}}_t^B}$, the predicted multi-target density $f({\mathit{X}}_{t|t-1}|Z_{1:t-1})$ is again of the LMB form with the parameters ${\left\{(r_{t|t-1}^{\ell },p_{t|t-1}^{\ell }(x_t))\right\}}_{\ell \in {\mathbb{L}}_{t-1}}\cup$${\left\{\left(r_{t|t-1}^{\ell },p_{t|t-1}^{\ell }(x_t)\right)\right\}}_{\ell \in {\mathbb{L}}_t^B}$ abbreviated as ${\left\{(r_{t|t-1}^{\ell },p_{t|t-1}^{\ell }(x_t))\right\}}_{\ell \in {\mathbb{L}}_t}$, where ${\mathbb{L}}_t={\mathbb{L}}_{t-1}\cup {\mathbb{L}}_t^B$ is subject to ${\mathbb{L}}_{t-1}\cap {\mathbb{L}}_t^B=\varnothing$. The target $\ell \in {\mathbb{L}}_{t-1}$ with the preceding state $x_{t-1}=({\mathit{x}}_{t-1},\ell )$ survives with probability $PS({\mathit{x}}_{t-1})$ and evolves via the single-target transition pdf $f^1(x_t|x_{t-1},\ell )=f^1(x_t|x_{t-1})$, so the $r_{t|t-1}^{\ell }$ and $p_{t|t-1}^{\ell }(x_t))$, $\forall \ell \in {\mathbb{L}}_{t-1}$, can be obtained by [28]

$$\{\begin{array}{l}{r}_{t|t-1}^{\ell }=r_{t-1}^{\ell }\langle p_{t-1}(\cdot ,\ell ),PS(\cdot )\rangle \\ p_{t|t-1}^{\ell }(x_t)\propto \langle p_{t-1}(\cdot ,\ell )PS(\cdot ),f^1(x_t|\cdot )\rangle \end{array}.$$

(7)

It is worth noting that the survival probability $PS({\mathit{x}}_{t-1})$ expresses uncertainty about whether the PT with the state ${\mathit{x}}_{t-1}$ still exists or has disappeared. It is equivalent to the prediction process of the existence probability in JIPDA if the Markov chain propagation formula only retains the transition from the existence of the target to itself. In addition, according to the Chapman–Kolmogorov equation, the LMB filter naturally incorporates the target birth process into the prediction [27]. The initialization of the Bernoulli density of the new target $\ell \in {\mathbb{L}}_t^B$ is not an exact science, so the state pdf and existence probability corresponding to the target $\ell \in {\mathbb{L}}_t^B$ can be obtained by one-point track initialization, two-point differencing, or other heuristic methods [20]. As the target birth process is another issue, a fixed-position birth process is used (see Section 5).

Similarly, given the prior density $f({\mathit{X}}_{t|t-1}|Z_{1:t-1})$ and the multi-sensor measurements $Z_t$ at time t, the updated multi-target density using the Bayes rule is of the form GLMB [32], which is expressed as

$$f\left(X_t|Z_{1:t}\right)=\Delta \left(X_t\right){\displaystyle \sum _{L\in ℱ\left({\mathbb{L}}_t\right)}{\displaystyle \sum _{{\theta }_t\in \Theta (L)}{w}^{L,{\theta }_t}}}{\delta }_L\left(ℒ\left(X_t\right)\right){\displaystyle \prod _{\left(x_t,\ell \right)\in X_t}{p}_t^{\ell ,{\theta }_t(\ell )}(x_t)}$$

(8)

where $w^{L,{\theta }_t}\propto {\displaystyle \prod _{{\ell }^{\prime }\in {\mathbb{L}}_t\backslash L}(1-r_{t|t-1}^{{\ell }^{\prime }})}{\displaystyle \prod _{\ell \in L}{r}_{t|t-1}^{\ell }{\varsigma }^{\ell ,{\theta }_t(\ell )}}$ with ${\varsigma }^{\ell ,{\theta }_t(\ell )}=\langle p_{t|t-1}^{\ell }(\cdot ),{\kappa }^{{\theta }_t(\ell )}(\cdot )\rangle$, and $p_t^{\ell ,{\theta }_t(\ell )}(x_t)=\frac{p_{t|t-1}^{\ell }(x_t){\kappa }^{{\theta }_t(\ell )}(x_t,\ell )}{{\varsigma }^{\ell ,{\theta }_t(\ell )}}$.

As discussed in [28,45], the multi-target density expressed in the LMB form can be recovered by the first moment approximation of the updated GLMB. The existence probability and the state pdf of the PT $\ell \in {\mathbb{L}}_t$ are obtained by

$$\{\begin{array}{l}{r}_t^{\ell }={\displaystyle \sum _{L\in ℱ\left({\mathbb{L}}_t\right)}{\displaystyle \sum _{{\theta }_t\in \Theta (L)}{1}_L(\ell )w^{L,{\theta }_t}}}\\ p_t^{\ell }(x_t)=\frac{1}{r_t^{\ell }}{\displaystyle \sum _{L\in ℱ\left({\mathbb{L}}_t\right)}{\displaystyle \sum _{{\theta }_t\in \Theta (L)}{1}_L(\ell )w^{L,{\theta }_t}}{p}_t^{\ell ,{\theta }_t(\ell )}(x_t)}\end{array}.$$

(9)

2.4. Inference Objective

From likelihood (6), we know that the target detection probability $P_{D_v}(x_t,\ell )$ is crucial in MTT algorithms but is usually unavailable a priori. Furthermore, approximating the LMB in the update step shown in (9) requires knowledge of all data associations ${\Theta }_t$, which can increase exponentially in multi-sensor scenarios. To this end, we wish to propose a method to overcome the challenges posed by unknown detection probabilities while effectively reducing the complexity of the multi-sensor LMB filter.

3. R-MS-LMB Model and the LMB Approximation

3.1. Robust Multi-Sensor LMB Model

Since the characteristics of the target detection probabilities may be unknown and time-varying, we assume that the target detection probabilities are variables to be estimated online, rather than known parameters. This section will explain how to estimate the multi-target state and target detection probabilities together.

To obtain a self-tuning R-MS-LMB tracking filter, we provide the basic idea that each target state is augmented by a variable that indicates the unknown detection probability. Therefore, the parameter $P_D(\cdot )$ can be estimated simultaneously in the update step. Based on the method in [42], the variable $e\in [0,1]$, representing the target detection probability, is augmented to the labeled state $\mathit{x}=(x,\ell )$. Notably, the detection probability depends on the ocean environment and the spatial distribution of the receiver and the sound source. Therefore, different S/R pairs may have different detection probabilities for the same PT. In light of this, we write the augmented state with label $\ell$ as $\tilde{\mathit{x}}=(x,\ell ,e^{\ell })$, where $e^{\ell }={[e_1^{\ell },\dots ,e_{\mathrm{V}}^{\ell }]}^T$ denotes the vector of the detection probability variables of V S/R pairs for PT $\ell$. The joint distribution of the augmented state $\tilde{\mathit{x}}=(x,\ell ,e^{\ell })$, $\forall \ell \in \mathbb{L}$, is defined as

$$\mu (x,\ell ,e^{\ell })=p^{\ell }\left(x\right){\displaystyle {\prod }_{v=1}^{\mathrm{V}}{\xi }^{\ell }\left(e_v^{\ell }\right)}$$

(10)

where ${\zeta }^{\ell }\left(e_v^{\ell }\right)$ is the statistic distribution of $e_v^{\ell }$. In the latter, we abbreviate $\mu (x,\ell ,e^{\ell })$ as ${\mu }^{\ell }\left(\tilde{x}\right)$ and denote $(x,e^{\ell })$ as $\tilde{x}$ for $\tilde{\mathit{x}}=(x,\ell ,e^{\ell })$. As a result of the augmented state $\tilde{x}=(x,\ell ,e^{\ell })$, the integral is as follows:

$${\displaystyle \int {\mu }^{\ell }\left(\tilde{x}\right)d}\tilde{x}={\displaystyle \int {\displaystyle {\int }_0^1\dots {\displaystyle {\int }_0^1{p}^{\ell }\left(x\right){\displaystyle {\prod }_{v=1}^{\mathrm{V}}{\xi }^{\ell }\left(e_v^{\ell }\right)}}}dx}d{e}_1^{\ell }\dots d{e}_{\mathrm{V}}^{\ell }.$$

(11)

At time $t-1$, the single PT with state ${\tilde{\mathit{x}}}_{t-1}=(x_{t-1},\ell ,e_{t-1}^{\ell })$ survives to the next time t with probability $PS({\tilde{\mathit{x}}}_{t-1})=PS({\mathit{x}}_{t-1})$ and moves with transition density $f(\cdot |{\tilde{\mathit{x}}}_{t-1})=f^1(\cdot |x){\displaystyle {\prod }_{v=1}^{\mathrm{V}}{f}^2(\cdot |e_{t-1,v}^{\ell })}$, where $f^1(\cdot |x_{t-1})$ denotes the single PT state x transition pdf defined in (7), and $f^2(\cdot |e_{t-1,v}^{\ell })$ denotes the transition probability pdf of $e_v^{\ell }$. Furthermore, the detection probability $P_{D_v}(x_t,\ell )$ for ${\tilde{\mathit{x}}}_t=(x_t,\ell ,e_t^{\ell })$ is given by $P_{D_v}(x_t,\ell )=e_{t,v}^{\ell }$.

Next, the density of the labeled Bernoulli RFS ${\tilde{\mathit{X}}}^{\ell }{}_t$ with the parameter $\{(r_t^{\ell },{\mu }^{\ell }\left(\tilde{x}\right))\}$ is given by

$$f^{\ell }\left({\tilde{\mathit{X}}}_t^{\ell }|Z_{1:t}\right)=\{\begin{array}{ll}1-r_t^{\ell }\hfill & {\tilde{\mathit{X}}}_t^{\ell }=\varnothing \hfill \\ r_t^{\ell }{p}_t^{\ell }(x){\displaystyle {\prod }_{v=1}^{\mathrm{V}}{\xi }_t^{\ell }\left(e_{t,v}^{\ell }\right)}\hfill & {\tilde{\mathit{X}}}_t^{\ell }=\{(x_t,\ell ,e_t^{\ell })\}\hfill \\ 0\hfill & \mathrm{Otherwise} \hfill \end{array}.$$

(12)

Similarly, an LMB realization augmented with the detection probability variables is represented by $\tilde{\mathit{X}}=\{(\tilde{x},{\ell }_1),\dots (\tilde{x},{\ell }_n)\}$ given by the union ${\tilde{\mathit{X}}}^{{\ell }_i}=\{(\tilde{x},{\ell }_i)\}$, i.e., $\tilde{\mathit{X}}={\cup }_{{\ell }_1}^{{\ell }_n}{\tilde{\mathit{X}}}^{{\ell }_i}$. Likewise, $ℒ(\tilde{\mathit{X}})$ is indicated by $ℒ(\tilde{\mathit{X}})=ℒ(\mathit{X})$. Replacing $\mathit{X}$ with $\tilde{\mathit{X}}$ in (8), we can rewrite the updated GLMB as follows:

$$f\left({\tilde{\mathit{X}}}_t|Z_{1:t}\right)=\Delta \left({\tilde{\mathit{X}}}_t\right){\displaystyle \sum _{L\in ℱ\left({\mathbb{L}}_t\right)}{\displaystyle \sum _{{\theta }_t\in \Theta (L)}{w}^{L,{\theta }_t}}}{\delta }_L\left(ℒ\left({\tilde{\mathit{X}}}_t\right)\right){\displaystyle \prod _{\left({\tilde{x}}_t,\ell \right)\in {\tilde{\mathit{X}}}_t}{\mu }_t^{\ell ,{\theta }_t(\ell )}({\tilde{x}}_t)}.$$

(13)

Here, ${\mu }_t^{\ell ,{\theta }_t(\ell )}({\tilde{x}}_t)=\frac{1}{{\varsigma }^{\ell ,{\theta }_t(\ell )}}{p}_{t|t-1}^{\ell }(x_t){\displaystyle {\prod }_{v=1}^{\mathrm{V}}{\xi }_{t|t-1}^{\ell }\left(e_{t|t-1,v}^{\ell }\right){\kappa }^{{\theta }_{t,v}(\ell )}({\tilde{x}}_t,\ell )}$ and $w^{L,{\theta }_t}\propto {\displaystyle \prod _{{\ell }^{\prime }\in {\mathbb{L}}_t\backslash L}(1-r_{t|t-1}^{{\ell }^{\prime }})}{\displaystyle \prod _{\ell \in L}{r}_{t|t-1}^{\ell }{\varsigma }^{\ell ,{\theta }_t(\ell )}}$, where the notations ${\varsigma }^{\ell ,{\theta }_t(\ell )}$, ${\kappa }^{{\theta }_t(\ell )}(x_t,\ell )$, and ${\kappa }_v^{{\theta }_{v,t}(\ell )}({\tilde{x}}_t,\ell )$ are redefined as ${\varsigma }^{\ell ,{\theta }_t(\ell )}={\displaystyle \int {\mu }_{t|t-1}^{\ell }(\tilde{x}),{\kappa }^{{\theta }_t(\ell )}(\tilde{x})d}\tilde{x}$, ${\kappa }^{{\theta }_t(\ell )}({\tilde{x}}_t,\ell )={\displaystyle {\prod }_{v=1}^{\mathrm{V}}{\kappa }_v^{{\theta }_{v,t}(\ell )}({\tilde{x}}_t,\ell )}$, and ${\kappa }_v^{{\theta }_{v,t}(\ell )}({\tilde{x}}_t,\ell )$$=\{\begin{array}{cc}\frac{e_{t|t-1,v}^{\ell }{l}_v(z_t|x)}{{\lambda }_{v}c(z_t^{{\theta }_{v,t}(\ell )})}& if {\theta }_{v,t}(\ell )\in {ℳ}_{v,t}\\ 1-e_{t|t-1,v}^{\ell }& if {\theta }_{v,t}(\ell )=0\end{array}$, respectively.

3.2. LMB Approximation

This section derives an approximation of LMB that differs from (9). Analyzing the relationship between measurements and PTs, rather than using the map ${\theta }_t$, we also define the PT–measurement associations related to the target-oriented variable ${\mathrm{a}}_{t,v}(\ell )\in \left\{-1,0,\dots ,M_{t,v}\right\}$ proposed in [45] for $\ell \in {\mathbb{L}}_t$ under the v-th S-R pair at time t:

$${\mathrm{a}}_{t,v}(\ell )=\{\begin{array}{cc}m\in \{1,\dots ,M_{t,v}\}& \mathrm{PT} \ell \mathrm{is} \mathrm{associated} \mathrm{with} z_{t,v}^m\in Z_{t,v},\\ 0& \mathrm{PT} \ell \mathrm{is} \mathrm{not} \mathrm{detected},\\ -1& \mathrm{PT} \ell \mathrm{does} \mathrm{not} \mathrm{exist}.\end{array}$$

(14)

As with an admissible map ${\theta }_{v,t}$, the admissible association variable ${\mathrm{a}}_{t,v}(\ell )$ should follow $\forall v=1,\dots ,V$, if ${\mathrm{a}}_{t,v}(\ell )\ne 0,{ \mathrm{a}}_{t,v}({\ell }^{\prime })\ne 0, \forall \ell \ne {\ell }^{\prime }\Rightarrow {\mathrm{a}}_{t,v}(\ell )\ne {\mathrm{a}}_{t,v}({\ell }^{\prime })$ (at most, one measurement can be assigned to a single PT, and no measurement can be assigned to more than one PT). Furthermore, we introduce association variables: $a_t(\ell )={[{\mathrm{a}}_{t,1}(\ell ),\dots ,{\mathrm{a}}_{t,V}(\ell )]}^T$ for a single PT $\ell \in {\mathbb{L}}_t$, and all PT cases $a_t={[a_t({\ell }_1),\dots ,a_t({\ell }_{\left|{\mathbb{L}}_t\right|})]}^T$. In addition, let us define by ${\mathcal{A}}_t$ and ${\mathcal{A}}_t^{\ell }$ the set of all admissible association variables $a_t$ and $a_t(\ell )$, respectively. To prevent pathological association variables, we let the association variable satisfy that $a_t(\ell )\in {\mathcal{A}}_t^{\ell }=\Lambda \uplus \{-{\mathbf{1}}_{\mathrm{V}}\}$ is a V-tuple that either pertains to $\Lambda$ when PT $\ell \in {\mathbb{L}}_t$ exists, or $\{-{\mathbf{1}}_{\mathrm{V}}\}$ if PT $\ell \in {\mathbb{L}}_t$ does not exist, where $\Lambda \triangleq \{0,\dots ,M_{t,1}\}\times ,\dots ,\times \{0,\dots ,M_{t,V}\}$ and ${\mathbf{1}}_{\mathrm{V}}$ denotes an identity vector with V elements.

The Kronecker delta function ${\delta }_L(ℒ({\tilde{\mathit{X}}}_t))$ with $L\in ℱ\left({\mathbb{L}}_t\right)$ in (13) serves to ensure that the labels $ℒ({\tilde{\mathit{X}}}_t)\in ℱ\left({\mathbb{L}}_t\right)$, so it can be expressed as an indicator function ${\displaystyle {\prod }_{({\tilde{x}}_t,\ell )\in {\tilde{\mathit{X}}}_t}{1}_{\mathbb{L}}(\ell )}$. We then rewrite the updated GLMB form (13) as

$$f\left({\tilde{\mathit{X}}}_t|Z_{1:t}\right)=\Delta \left({\tilde{\mathit{X}}}_t\right){\displaystyle \sum _{{\theta }_t\in {\Theta }_t}}{w}^{{\theta }_t}{\displaystyle \prod _{\left(x_t,\ell \right)\in {\tilde{\mathit{X}}}_t}{1}_{{\mathbb{L}}_t}(\ell ){\mu }_t^{\ell ,{\theta }_t(\ell )}({\tilde{x}}_t)}.$$

(15)

Next, we replace ${\theta }_t$ with $a_t$ to obtain the following:

$$f\left({\tilde{\mathit{X}}}_t|Z_{1:t}\right)=\Delta \left({\tilde{\mathit{X}}}_t\right){\displaystyle {\sum }_{a_t\in {\mathcal{A}}_t}\rho (a_t,{\tilde{\mathit{X}}}_t)w^{a_t}}{\displaystyle \prod _{\left(x_t,\ell \right)\in {\tilde{\mathit{X}}}_t}{1}_{{\mathbb{L}}_t}(\ell ){\mu }_t^{\ell ,{\theta }_t(\ell )}({\tilde{x}}_t)}.$$

(16)

Here, we define the function $\rho (a_t,{\tilde{\mathit{X}}}_t)$ to ensure that the associated variable $a_t$ is admissible, where $\rho (a_t,{\tilde{\mathit{X}}}_t)=1$ if $a_t\in {\mathcal{A}}_t$ is subject to $a_t(\ell )\in \Lambda$ for $\ell \in ℒ\left({\tilde{\mathit{X}}}_t\right)$ and $a_t(\ell )=-{\mathbf{1}}_{\mathrm{V}}$ for $\ell \in {\mathbb{L}}_t\backslash ℒ\left({\tilde{\mathit{X}}}_t\right)$, and $\rho (a_t,{\tilde{\mathit{X}}}_t)=0$ otherwise. As $a_t(\ell )\in \Lambda$ implies that $({\tilde{x}}_t,\ell )\in {\tilde{\mathit{X}}}_t$, ${\mu }_t^{\ell ,{\theta }_t(\ell )}({\tilde{x}}_t)$ in (16) is equivalent to ${\mu }_t^{\ell ,{\theta }_t(\ell )}({\tilde{x}}_t,\ell )$ in (13) with the substitution of $a_t(\ell )$ for ${\theta }_t(\ell )$ concerning $\ell \in ℒ\left({\tilde{\mathit{X}}}_t\right)$. Additionally, the weights $w^{a_t}$ can be represented by

$$w^{a_t}\propto {\displaystyle {\prod }_{\ell \in {\mathbb{L}}_t}{\tau }^{\ell ,a_t(\ell )}}$$

(17)

where ${\tau }^{\ell ,a_t(\ell )}=(1-r_{t|t-1}^{\ell })$ if $a_t(\ell )=-{\mathbf{1}}_{\mathrm{V}}$, and ${\tau }^{\ell ,a_t(\ell )}=r_{t|t-1}^{\ell }{\varsigma }^{\ell ,a_t(\ell )}$ if $a_t(\ell )\in \Lambda$. Compared with $w^{L,{\theta }_t}$ in (13), $w^{a_t}$ in (17) depends on ${\mathbb{L}}_t$ instead of $ℒ\left({\tilde{\mathit{X}}}_t\right)$, so the weight $w^{a_t}$ can generally be interpreted as the probability mass function (pmf) of $a_t$. Therefore, we can represent the pmf of $a_t\in {\mathcal{A}}_t$ by $g(a_t)=w^{a_t}$, and normalize $g(a_t)$ in that ${\displaystyle {\sum }_{a_t\in \mathcal{A}}g(a_t)}=1$.

More precisely, the $\rho (a_t,{\tilde{\mathit{X}}}_t)$ can be expanded as $\rho (a_t,{\tilde{\mathit{X}}}_t)={\displaystyle \prod _{\ell \in {\mathbb{L}}_t}^{{\uplus }_{\ell \in \mathbb{L}}{\tilde{\mathit{X}}}^{\ell }={\tilde{\mathit{X}}}_t}\rho \left(a_t(\ell ),{\tilde{\mathit{X}}}_t^{\ell }\right)}$. Here, we introduce

$$\rho \left(a_t(\ell ),{\tilde{\mathit{X}}}_t^{\ell }\right)=\{\begin{array}{cc}1& \begin{array}{l}\left(a_t(\ell )\in \Lambda \mathrm{and} {\tilde{\mathit{X}}}_t^{\ell }=\{(x_t,\ell ,e_t^{\ell })\}\right)\\ \mathrm{or} \left(a_t(\ell )=-{\mathbf{1}}_{\mathrm{V}} \mathrm{and} {\tilde{\mathit{X}}}_t^{\ell }=\varnothing \right)\end{array}\\ 0& \mathrm{Otherwise}\end{array}$$

(18)

Thus, introducing $g(a_t)$ and expanding $\rho (a_t,{\tilde{\mathit{X}}}_t)$ in (16) yields

$$f\left({\tilde{\mathit{X}}}_t|Z_{1:t}\right)=\Delta \left({\tilde{\mathit{X}}}_t\right){\displaystyle {\sum }_{a_t\in {\mathcal{A}}_t}g(a_t){\displaystyle \prod _{\ell \in {\mathbb{L}}_t}^{{\uplus }_{\ell \in \mathbb{L}}{\tilde{\mathit{X}}}^{\ell }={\tilde{\mathit{X}}}_t}\rho \left(a_t(\ell ),{\tilde{\mathit{X}}}_t^{\ell }\right){\phi }^{\ell ,a_t(\ell )}({\tilde{\mathit{X}}}_t^{\ell })}}$$

(19)

Here, we define ${\phi }^{\ell ,a_t(\ell )}({\tilde{\mathit{X}}}_t^{\ell })=\{\begin{array}{ll}1\hfill & {\tilde{\mathit{X}}}_t^{\ell }=\varnothing \hfill \\ p_t^{\ell ,a_t(\ell )}(x){\displaystyle {\prod }_{v=1}^{\mathrm{V}}{\xi }_t^{\ell ,a_t(\ell )}\left(e_v^{\ell }\right)}\hfill & {\tilde{\mathit{X}}}_t^{\ell }=\{(x_t,\ell ,e_t^{\ell })\}\hfill \\ 0\hfill & \mathrm{otherwise} \hfill \end{array}$ to indicate the conditional density in the existence of PT $\ell$ corresponding to the single labeled Bernoulli ${\tilde{\mathit{X}}}_t^{\ell }$.

Given the GLMB density $f\left({\tilde{\mathit{X}}}_t|Z_{1:t}\right)$ in (19), the best LMB approximation of the GLMB as $\stackrel{⌣}{f}\left({\tilde{\mathit{X}}}_t|Z_{1:t}\right)=\Delta ({\tilde{\mathit{X}}}_t){\displaystyle \prod _{\ell \in {\mathbb{L}}_t}^{{\uplus }_{\ell \in \mathbb{L}}{\tilde{\mathit{X}}}^{\ell }={\tilde{\mathit{X}}}_t}{\stackrel{⌣}{f}}^{\ell }\left({\tilde{\mathit{X}}}^{\ell }{}_t|Z_{1:t}\right)}$ using KLD minimization is derived:

$$\begin{array}{l}\mathrm{D}(f\Vert \stackrel{⌣}{f})={\displaystyle \int f({\tilde{\mathit{X}}}_t|Z_{1:t})}\log \frac{f({\tilde{\mathit{X}}}_t|Z_{1:t})}{\stackrel{⌣}{f}\left({\tilde{\mathit{X}}}_t|Z_{1:t}\right)}\delta {\tilde{\mathit{X}}}_t\\ ={\displaystyle \iint \dots }{\displaystyle \int f}\left({\tilde{\mathit{X}}}^{{\ell }_1}\uplus ,\dots ,\uplus {\tilde{\mathit{X}}}^{{\ell }_{\left|\mathbb{L}\right|}}\right)\log \frac{f\left({\tilde{\mathit{X}}}_t^{{\ell }_1}\uplus ,\dots ,\uplus {\tilde{\mathit{X}}}_t^{{\ell }_{\left|{\mathbb{L}}_t\right|}}\right)}{\stackrel{⌣}{f}\left({\tilde{\mathit{X}}}_t^{{\ell }_1}\uplus ,\dots ,\uplus {\tilde{\mathit{X}}}_t^{{\ell }_{\left|{\mathbb{L}}_t\right|}}\right)}\delta {\tilde{\mathit{X}}}_t^{{\ell }_1}\uplus ,\dots ,\uplus {\tilde{\mathit{X}}}_t^{{\ell }_{\left|{\mathbb{L}}_t\right|}}\\ =c-{\displaystyle {\sum }_{\ell \in {\mathbb{L}}_t}{\displaystyle \int {\displaystyle {\sum }_{a_t\in {\mathcal{A}}_t}p(a_t){\phi }^{\ell ,a_t(\ell )}({\tilde{\mathit{X}}}_t^{\ell }){\displaystyle \prod _{{\ell }^{\prime }\in {\mathbb{L}}_t}\rho \left(a_t({\ell }^{\prime }),{\tilde{\mathit{X}}}_t^{{\ell }^{\prime }}\right)}\log {\stackrel{⌣}{f}}^{\ell }\left({\tilde{\mathit{X}}}^{\ell }{}_t|Z_{1:t}\right)\delta {\tilde{\mathit{X}}}^{\ell }}}}\end{array}$$

(20)

where ${\stackrel{⌣}{f}}^{\ell }\left({\tilde{\mathit{X}}}^{\ell }{}_t|Z_{1:t}\right)$ is the approximated density of single labeled Bernoulli ${\tilde{\mathit{X}}}^{\ell }$ for PT $\ell \in {\mathbb{L}}_t$, and c is a constant that is independent of $\stackrel{⌣}{f}\left({\tilde{\mathit{X}}}_t|Z_{1:t}\right)$. Maximizing with respect to ${\stackrel{⌣}{f}}^{\ell }\left({\tilde{\mathit{X}}}^{\ell }{}_t|Z_{1:t}\right)$, we obtain

$${\stackrel{⌣}{f}}^{\ell }\left({\tilde{\mathit{X}}}^{\ell }{}_t|Z_{1:t}\right)={\displaystyle {\sum }_{a_t\in {\mathcal{A}}_t}g(a_t){\phi }^{\ell ,a_t(\ell )}({\tilde{\mathit{X}}}_t^{\ell }){\displaystyle \prod _{{\ell }^{\prime }\in {\mathbb{L}}_t}\rho \left(a_t(\ell ),{\tilde{\mathit{X}}}_t^{{\ell }^{\prime }}\right)}}.$$

(21)

Next, we approximate $g(a_t)$ based on its marginal product, i.e.,

$$g(a_t)\approx {\displaystyle {\prod }_{\ell \in {\mathbb{L}}_t}g(a_t(\ell ))}$$

(22)

where $g(a_t(\ell ))={\displaystyle {\sum }_{{\widehat{a}}_t\in {\mathcal{A}}_t:{\widehat{a}}_t(\ell )=a_t(\ell )}g(a_t)}$, and we can redefine (21) as

$${\stackrel{⌣}{f}}^{\ell }\left({\tilde{\mathit{X}}}^{\ell }{}_t|Z_{1:t}\right)={\displaystyle {\sum }_{a_t\in {\mathcal{A}}_t}{\phi }^{\ell ,a_t(\ell )}({\tilde{\mathit{X}}}_t^{\ell }){\displaystyle {\prod }_{{\ell }^{\prime }\in {\mathbb{L}}_t}g(a_t({\ell }^{\prime }))\rho \left(a_t({\ell }^{\prime }),{\tilde{\mathit{X}}}_t^{{\ell }^{\prime }}\right)}}.$$

(23)

Splitting ${\displaystyle {\sum }_{a_t\in {\mathcal{A}}_t}}$ as ${\displaystyle {\sum }_{a_t\in {\mathcal{A}}_t}=}$${\displaystyle {\sum }_{a_t(\sim \ell )\in {\mathcal{A}}_t^{\sim \ell }}{\displaystyle {\sum }_{a_{t(\ell )}\in {\mathcal{A}}_t^{\ell }}}}$ and exploiting the sum separability, i.e., ${\displaystyle {\sum }_i{\displaystyle {\sum }_j{g}_i{h}_j}}={\displaystyle {\sum }_i{g}_i}{\displaystyle {\sum }_j{h}_j}$, we obtain the following result:

$$\begin{array}{l}{\stackrel{⌣}{f}}^{\ell }\left({\tilde{\mathit{X}}}_t^{{\ell }_n}|Z_{1:t}\right)\\ ={\displaystyle {\sum }_{a_{t(\ell )}\in {\mathcal{A}}_t^{\ell }}\rho (a_t(\ell ),{\tilde{\mathit{X}}}^{\ell }){\phi }^{\ell ,a_t(\ell )}({\tilde{\mathit{X}}}^{\ell })}{\displaystyle {\sum }_{a_{t(\sim \ell )}\in {\mathcal{A}}_t^{\sim \ell }}{\displaystyle {\prod }_{{\ell }^{\prime }\in {\mathbb{L}}_t\backslash \ell }p(a_t({\ell }^{\prime }))}\rho (a_t({\ell }^{\prime }),{\tilde{\mathit{X}}}^{{\ell }^{\prime }})}\\ ={\displaystyle {\sum }_{a_{t(\ell )}\in {\mathcal{A}}_t^{\ell }}\rho \left(a_t(\ell ),{\tilde{\mathit{X}}}^{\ell }\right){\phi }^{\ell ,a_t(\ell )}({\tilde{\mathit{X}}}^{\ell })}\end{array}$$

(24)

where $a_t(\sim \ell )$ and ${\mathcal{A}}_t^{\sim \ell }$ are the vector $a_t$ with the component $a_t(\ell )$ removed and the set of all admissible association variables $a_t(\sim \ell )$, respectively.

Next, by comparing the expression ${\stackrel{⌣}{f}}^{\ell }\left({\tilde{\mathit{X}}}^{\ell }{}_t|Z_{1:t}\right)$ in (24) with the density of the label Bernoulli in (12) and evaluating $\rho \left(a_t(\ell ),{\tilde{\mathit{X}}}^{\ell }\right)$, we can obtain the approximated density of the single labeled Bernoulli ${\tilde{\mathit{X}}}_t^{\ell }$ for $\ell \in {\mathbb{L}}_t$ by

$${\stackrel{⌣}{f}}^{\ell }\left({\tilde{\mathit{X}}}_t^{\ell }|Z_{1:t}\right)=\{\begin{array}{ll}g\left(a_t(\ell )=-{\mathbf{1}}_{\mathrm{V}}\right)\hfill & {\tilde{\mathit{X}}}^{\ell }=\varnothing \hfill \\ \frac{1}{{\displaystyle {\sum }_{a_t(\ell )\in \Lambda }g\left(a_t(\ell )\right)}}{\displaystyle {\sum }_{a_t(\ell )\in \Lambda }g\left(a_t(\ell )\right){\mu }_t^{\ell ,a_t(\ell )}({\tilde{x}}_t)}\hfill & {\tilde{\mathit{X}}}^{\ell }=\{(x,\ell ,e^{\ell })\}\hfill \\ 0\hfill & \mathrm{otherwise} \hfill \end{array}.$$

(25)

Moreover, using (25), the updated existence probability and augmented state pdf for PT $\ell \in {\mathbb{L}}_t$ are given by

$$\{\begin{array}{l}{r}_t^{\ell }=1-g\left(a_t(\ell )=-{\mathbf{1}}_{\mathrm{V}}\right)={\displaystyle {\sum }_{a_t(\ell )\in \Lambda }g\left(a_t(\ell )\right)}\\ {\mu }_t^{\ell }({\tilde{x}}_t)=\frac{1}{r_t^{\ell }}{\displaystyle {\sum }_{a_t(\ell )\in \Lambda }g\left(a_t(\ell )\right){\mu }_t^{\ell ,a_t(\ell )}({\tilde{x}}_t)}\end{array}.$$

(26)

Finally, based on (26), the updated LMB can be parameterized by ${\left\{(r_t^{\ell },{\mu }_t^{\ell }({\tilde{x}}_t))\right\}}_{\ell \in {\mathbb{L}}_t}$. Therefore, the approximated posterior density of LMB after the update can be expressed as follows:

$$\begin{array}{l}f\left({\tilde{\mathit{X}}}_t\mid Z_{1:t}\right)\approx \stackrel{⌣}{f}\left({\tilde{\mathit{X}}}_t|Z_{1:t}\right)\\ =\Delta ({\tilde{\mathit{X}}}_t){\displaystyle {\prod }_{{\ell }^{\prime }\in {\mathbb{L}}_t\backslash ℒ({\tilde{\mathit{X}}}_t)}g\left(a_t({\ell }^{\prime })=-{\mathbf{1}}_{\mathrm{V}}\right)}{\displaystyle {\prod }_{\ell \in ℒ({\tilde{\mathit{X}}}_t)}{\displaystyle {\sum }_{a_t(\ell )\in \Lambda }g\left(a_t(\ell )\right){\mu }_t^{\ell ,a_t(\ell )}({\tilde{x}}_t)}}\end{array}$$

(27)

Above, we derive the LMB approximation form obtained from (26) and (27). Unlike the form in (9), the LMB approximation in (27) requires only performing marginalization over the association variable $a_t$ instead of $(L,{\theta }_t)$ with $L\in ℱ\left({\mathbb{L}}_t\right)$ and ${\theta }_t\in \Theta (L)$, which avoids the solution of the maps $\Theta (L)$ for all subsets $L\in ℱ\left({\mathbb{L}}_t\right)$. As a result, the approximation problem of LMB in the update step can be transformed into the issue of marginalization.

In addition, the approximate LMB density in (27) is very similar to JIPDA. The difference is that the LMB introduces labels to manage the target tracks and incorporates a Bayesian model for target birth and death. Another difference is the data association, where the existence of the target not detected and the non-existence of the target are separated into two events, while, in JIPDA, the two events are merged and interpreted as one event in which the target is not associated with any measurement.

4. The BP-Based Framework of the R-MS-LMB Filter

4.1. Joint Posterior Distributions

Recalling Theorem 1, the approximation of LMB becomes the problem of calculating the marginal posterior distribution quickly and efficiently at each time step. Thus, we are motivated to infer the multi-target states and detection probabilities by marginalizing the joint posterior distribution. Given the measurements of $Z_{1:t}$, we present the joint posterior distribution with association variables as follows:

$$f\left({\tilde{\mathit{X}}}_t,a_t|Z_{1:t}\right)=\Delta \left({\tilde{\mathit{X}}}_t\right)g(a_t){\displaystyle \prod _{\ell \in {\mathbb{L}}_t}^{{\uplus }_{\ell \in \mathbb{L}}{\tilde{\mathit{X}}}^{\ell }={\tilde{\mathit{X}}}_t}\rho \left(a_t(\ell ),{\tilde{\mathit{X}}}_t^{\ell }\right){\phi }^{\ell ,a_t(\ell )}({\tilde{\mathit{X}}}_t^{\ell })}.$$

(28)

Thus, the target-oriented association variables $a_{t,v}$ and $a_t$ have been introduced. Following [50,51], we describe the alternative PT–measurement associations in terms of the measurement-oriented association variables $b_{t,v}(m)$, $\forall m\in \{0,\dots ,M_{t,v}\}$ under the v-th S/R pair at time t:

$$b_{t,v}(m)=\{\begin{array}{cc}\ell \in {\mathbb{L}}_t& \begin{array}{l}\mathrm{Measurement} z_{t,v}^m \mathrm{is} \mathrm{associated}\\ \mathrm{with} \mathrm{PT} \ell \end{array}\\ 0& z_{t,v}^m \mathrm{is} \mathrm{not} \mathrm{from} \mathrm{any} \mathrm{PTs}\end{array}$$

(29)

Similarly, the stacked vector of association variables for all S/R pairs is defined as $b_t=(b_{t,1}\cdots b_{t,V})$. There is a redundant relationship between $a_t$ and $b_t$, which means that it is possible to derive $a_t$ from ${\mathsf{b}}_t$ and vice versa. The indicator function enforces the admissibility and consistency between $a_t$ and $b_t$ [41]

$$\Psi (a_t,b_t)\triangleq {\displaystyle {\prod }_{\ell \in {\mathbb{L}}_t}{\displaystyle {\prod }_{v=1}^{\mathrm{V}}{\displaystyle {\prod }_{m=1}^{M_{t,v}}\Psi \left(a_{t,v}(\ell ),b_{t,v}(m)\right)}}}$$

(30)

with $\Psi \left(a_{t,v}(\ell ),b_{t,v}(m)\right)\triangleq \{\begin{array}{cc}0& \begin{array}{l}\mathrm{if} a_{t,v}(\ell )=m \mathrm{and} b_{t,v}(m)\ne \ell \mathrm{or}\\ a_{t,v}(\ell )\ne m \mathrm{and} b_{t,v}(m)=\ell , \end{array}\\ 1& \mathrm{otherwise}\end{array}$.

The joint association pmf of $a_t$ and $b_t$ can be expressed by $p(a_t,b_t)=p(a_t)\Psi (a_t,b_t)$. It has been demonstrated in [50,52] that the joint association pmf of $a_t$ and redundant $b_t$ allows a fast parallel method to be developed for BP-based probabilistic data association. Introducing the joint pmf $p(a_t,b_t)$ into (28), we obtain the joint posterior distribution:

$$f\left({\tilde{\mathit{X}}}_t,a_t,b_t|Z_{1:t}\right)=\Delta \left({\tilde{\mathit{X}}}_t\right)g(a_t)\Psi (a_t,b_t){\displaystyle \prod _{\ell \in {\mathbb{L}}_t}^{{\uplus }_{\ell \in \mathbb{L}}{\tilde{\mathit{X}}}^{\ell }={\tilde{\mathit{X}}}_t}\rho \left(a_t(\ell ),{\tilde{\mathit{X}}}_t^{\ell }\right){\phi }^{\ell ,a_t(\ell )}({\tilde{X}}^{\ell })}.$$

(31)

Finally, the main challenge is to develop a computationally feasible way of calculating the marginal $f\left({\tilde{\mathit{X}}}_t,a_t,b_t|Z_{1:t}\right)$ of (31).

4.2. Factorization of the Joint Posterior Distribution

Joint posterior distribution $f\left({\tilde{\mathit{X}}}_t,a_t,b_t|Z_{1:t}\right)$ involves LMB, association variables, detection probabilities, and measurements at time t. This results in the exponential growth of the total computational amount with increasing PTs and measurements, so marginalization would be infeasible in practice. Fortunately, we can implement effective marginalization by implementing the BP scheme on a factor graph embedded in the joint posterior distribution’s factorization.

To derive this factorization, we replace $g(a_t)$ with ${\displaystyle {\prod }_{\ell \in \mathbb{L}}{\tau }^{\ell ,a_t(\ell )}}$ in (31), which yields

$$\begin{array}{l}f\left({\tilde{\mathit{X}}}_t,a_t,b_t|Z_{1:t}\right)=\Delta \left({\tilde{\mathit{X}}}_t\right)g(a_t)\Psi (a_t,b_t){\displaystyle \prod _{\ell \in {\mathbb{L}}_t}^{{\uplus }_{\ell \in \mathbb{L}}{\tilde{\mathit{X}}}^{\ell }={\tilde{\mathit{X}}}_t}\rho \left(a_t(\ell ),{\tilde{\mathit{X}}}_t^{\ell }\right){\phi }^{\ell ,a_t(\ell )}({\tilde{\mathit{X}}}_t^{\ell })}\\ =\Delta \left({\tilde{\mathit{X}}}_t\right)\Psi (a_t,b_t){\displaystyle \prod _{\ell \in {\mathbb{L}}_t}^{{\uplus }_{\ell \in \mathbb{L}}{\tilde{\mathit{X}}}^{\ell }={\tilde{\mathit{X}}}_t}{\tau }^{\ell ,a_t(\ell )}\rho \left(a_t(\ell ),{\tilde{\mathit{X}}}_t^{\ell }\right){\phi }^{\ell ,a_t(\ell )}({\tilde{\mathit{X}}}_t^{\ell })}\end{array}$$

(32)

Introducing ${\tau }^{\ell ,a_t(\ell )}$, $\Psi (a_t,b_t)$, $\rho \left(a_t(\ell ),{\tilde{\mathit{X}}}_t^{\ell }\right)$, and ${\phi }^{\ell ,a_t(\ell )}({\tilde{X}}^{\ell })$ into (32), we evaluate and group the terms to obtain the final factorization:

$$f\left({\tilde{\mathit{X}}}_t|Z_{1:t}\right)=\Delta \left({\tilde{\mathit{X}}}_t\right){\displaystyle \prod _{\ell \in {\mathbb{L}}_t}^{{\uplus }_{\ell \in \mathbb{L}}{\tilde{\mathit{X}}}^{\ell }={\tilde{\mathit{X}}}_t}f({\tilde{\mathit{X}}}_{t|t-1}^{\ell }){\displaystyle {\prod }_{v=1}^{\mathrm{V}}{v}_{t,v}({\tilde{\mathit{X}}}_{t|t-1}^{\ell },{\mathrm{a}}_{t,v}(\ell );Z_{t,v}){\displaystyle {\prod }_{m=1}^{M_{t,v}}\Psi \left(a_{t,v}(\ell ),b_{t,v}(m)\right)}}}$$

(33)

where $v_{t,v}({\tilde{\mathit{X}}}_{t|t-1}^{\ell },{\mathrm{a}}_{t,v}(\ell );Z_{t,v})=\{\begin{array}{cc}\frac{e_{t|t-1,v}^{\ell }{l}_v(z_{t,v}^m|x_{t|t-1})}{{\lambda }_{v}c(z_t^m)}& \mathrm{if} a_{t,v}(\ell )=m>0 \mathrm{and} {\tilde{\mathit{X}}}_{t|t-1}^{\ell }=\{({\tilde{x}}_{t|t-1},\ell )\}\\ 1-e_{t|t-1,v}^{\ell }& if a_{t,v}(\ell )=0 \mathrm{and} {\tilde{\mathit{X}}}_{t|t-1}^{\ell }=\{({\tilde{x}}_{t|t-1},\ell ),\\ 1& if a_t^s(\ell )=-1 \mathrm{and} {\tilde{\mathit{X}}}_{t|t-1}^{\ell }=\varnothing ,\\ 0& \mathrm{Otherwise}\end{array}$.

We show the factor graph corresponding to (33) at time t in Figure 3. In the factor graph, the nodes and edges constitute the factor graph, where the edges connect two nodes, the square nodes represent the functions of the variables, and the circle nodes represent the variables.

4.3. BP Scheme

We will detail the BP-based scheme for inferring the marginal densities of (33). BP based on a factor graph, also known as the sum–product algorithm [53], is an iterative message-passing algorithm that efficiently calculates the marginal posterior distribution of each expected variable. In BP, the calculated belief approximates the marginal posterior distribution. The message is the information transmitted from one node to another, and a node in the factor graph sends messages to its adjacent nodes. Readers can learn more about BP by referring to [54,55].

At the beginning of time step t, the prediction message ${\chi }^{\ell }$ with the parameter $\left\{(r_{t|t-1}^{\ell },{\mu }_{t|t-1}^{\ell }(\tilde{x}))\right\}$ for PT $\ell \in {\mathbb{L}}_t$ is given by

$${\chi }^{\ell }({\tilde{\mathit{X}}}_{t|t-1}^{\ell })=f_{t|t-1}^{\ell }\left({\tilde{\mathit{X}}}_{t|t-1}^{\ell }|Z_{1:t-1}\right)=\{\begin{array}{ll}1-r_{t|t-1}^{\ell }\hfill & {\tilde{\mathit{X}}}_{t|t-1}^{\ell }=\varnothing \hfill \\ r_{t|t-1}^{\ell }{p}_{t|t-1}^{\ell }(x){\displaystyle {\prod }_{v=1}^{\mathrm{V}}{\xi }_{t|t-1}^{\ell }\left(e_v^{\ell }\right)}\hfill & {\tilde{\mathit{X}}}_{t|t-1}^{\ell }=\{(x_{t|t-1},\ell ,e_{t|t-1}^{\ell })\}\hfill \\ 0\hfill & \mathrm{Otherwise} \hfill \end{array}$$

(34)

where $f_{t|t-1}^{\ell }\left({\tilde{\mathit{X}}}_{t|t-1}^{\ell }|Z_{1:t-1}\right)$ can be determined by the standard LMB filter prediction.

The message-passing iteration scheme is divided into the outer BP loop and inner BP loop, as shown in Figure 1. The inner BP loop focuses on the data association. Synchronously and in parallel for all S/R pairs, the outer BP executes the iterative message-passing scheme P times. Next, we also point out the corresponding messages in the factor graph of Figure 1.

Prior to the implementation of message passing, the message ${\partial }_v^{\ell }({\tilde{\mathit{X}}}_{t|t-1}^{\ell })$ with parameter $\left\{(r_{\ast }^{\ell },b_{\ast }^{\ell }(\tilde{x}))\right\}$ (sent from variable node ”${\tilde{\mathit{X}}}_{}^{\ell }$” to factor node “$v_v^{\ell }$”) is initialized with the prediction message ${\chi }^{\ell }$ with the parameter $\left\{(r_{t|t-1}^{\ell },{\mu }_{t|t-1}^{\ell }(\tilde{x}))\right\}$, $\forall \ell \in {\mathbb{L}}_t$, such that $r_{\ast }^{\ell }=r_{t|t-1}^{\ell }$, $b_{\ast }^{\ell }(\tilde{x})={\mu }_{t|t-1}^{\ell }(\tilde{x})$. The association message ${\beta }_v^{\ell }(a_{t,v}(\ell ))$ from node “$v_v^{\ell }$” to node “$a_v^{\ell }$” is then computed locally:

$${\beta }_v^{\ell }(a_{t,v}(\ell ))={\displaystyle \int }{v}_{t,v}({\tilde{X}}_{}^{\ell },{\mathrm{a}}_{t,v}(\ell );Z_{t,v}){\partial }_v^{\ell }({\tilde{\mathit{X}}}^{\ell })\delta {\tilde{\mathit{X}}}^{\ell }=\{\begin{array}{cc}{\displaystyle \int }{r}_{\ast }^{\ell }\frac{e_{*,v}^{\ell }{l}_v(z_t|x)}{{\lambda }_{v}c(z_t^{{\theta }_{v,t}(\ell )})}{b}_{\ast }^{\ell }(\tilde{x})d\tilde{x}& a_{t,v}(\ell )=m>0\\ r_{\ast }^{\ell }-r_{\ast }^{\ell }{\displaystyle \int e_{*,v}^{\ell }{b}_{\ast }^{\ell }(\tilde{x})d\tilde{x}}& a_{t,v}(\ell )=0\\ 1-r_{\ast }^{\ell }& a_{t,v}(\ell )=-1\end{array}.$$

(35)

The message ${\beta }_v^{\ell }(a_{t,v}(\ell ))$ in the first outer loop iteration is evaluated by marginalizing the predicted density of PT $\ell$. Next, the message ${\beta }_v^{\ell }(a_{t,v}(\ell ))$ is used to evaluate the messages ${\eta }_v^{\ell }(a_{t,v}(\ell ))$ by the inner BP loop. Descriptions of an efficient implementation of the inner BP loop with MATLAB code are available in [56].

Subsequently, the likelihood messages ${\gamma }_v^{\ell }({\tilde{\mathit{X}}}_{t|t-1}^{\ell })$ sent from the factor node “$v_v^{\ell }$” to the node “${\tilde{\mathit{X}}}_{}^{\ell }$” are as follows:

$$\begin{array}{l}{\gamma }_v^{\ell }({\tilde{\mathit{X}}}_{t|t-1}^{\ell })={\displaystyle \sum _{m=0}^{M_s}}{\eta }_v^{\ell }(m)v_{t,v}({\tilde{\mathit{X}}}_{t|t-1}^{\ell },m;Z_{t,v})+{\eta }_v^{\ell }(-1)v_{t,v}({\tilde{\mathit{X}}}_{t|t-1}^{\ell },-1;Z_{t,v})\\ =\{\begin{array}{cc}{\eta }_v^{\ell }(0)\cdot \left[1-e_{t|t-1,v}^{\ell }\right]+{\displaystyle \sum _{m=1}^{M_s}}{\eta }_v^{\ell }(m)\frac{e_{t|t-1,v}^{\ell }{l}_v(z_{t,v}^{\ell }|x_{t|t-1},\ell )}{{\lambda }^{v}c(z_{t,v}^m)}& {\tilde{\mathit{X}}}_{t|t-1}^{\ell }=\{({\tilde{x}}_{t|t-1},\ell )\}\\ {\eta }_v^{\ell }(-1)& {\tilde{\mathit{X}}}_{t|t-1}^{\ell }=\varnothing \end{array}\end{array}$$

(36)

Finally, at the last iteration of the outer loop, given the messages ${\gamma }_v^{\ell }({\tilde{\mathit{X}}}_{t|t-1}^{\ell })$ for PT $\ell \in {\mathbb{L}}_t$ from all the S/R pairs $v=1,\dots ,\mathrm{V}$, we approximate the marginal posterior pdf $f_t^{\ell }\left({\tilde{\mathit{X}}}_t^{\ell }|Z_{1:t}\right)$ by the belief $b({\tilde{\mathit{X}}}_t^{\ell })$ with the parameter $\{(r_t^{\ell },b_t^{\ell }(\tilde{x}))\}$:

$$f_t^{\ell }\left({\tilde{\mathit{X}}}_t^{\ell }|Z\right)=b({\tilde{\mathit{X}}}_t^{\ell })\propto {\chi }^{\ell }({\tilde{\mathit{X}}}_{t|t-1}^{\ell }){\displaystyle {\prod }_{v\in \mathcal{V}}{\gamma }_v^{\ell }({\tilde{\mathit{X}}}_{t|t-1}^{\ell })}=\{\begin{array}{cc}{r}_t^{\ell }{\mu }_t^{\ell }(\tilde{x})& {\tilde{\mathit{X}}}_t^{\ell }=\{({\tilde{x}}_t,\ell )\}\\ (1-r_t^{\ell })& {\tilde{\mathit{X}}}_t^{\ell }=\varnothing \end{array}$$

(37)

where

$$r_t^{\ell }=\frac{r_{t|t-1}{\displaystyle {\prod }_{v\in \mathcal{V}}(1-{\eta }_v^{\ell }(-1))}}{(1-r_{t|t-1}){\displaystyle {\prod }_{v\in \mathcal{V}}{\eta }_v^{\ell }(-1)}+r_{t|t-1}{\displaystyle {\prod }_{v\in \mathcal{V}}(1-{\eta }_v^{\ell }(-1))}}$$

(38)

$${\mu }_t^{\ell }({\tilde{x}}_t)\propto {\mu }_t^{\ell }({\tilde{x}}_{t|t-1}){\displaystyle {\prod }_{v\in \mathcal{V}}\left({\eta }_v^{\ell }(0)\cdot \left[1-e_{t|t-1,v}^{\ell }\right]+{\displaystyle \sum }_{m=1}^{M_s}{\eta }_v^{\ell }(m)\frac{e_{t|t-1,v}^{\ell }{l}_v(z_v^{\ell }|x_{t|t-1},\ell )}{{\lambda }^{v}c(z_v^m)}\right)}.$$

(39)

Furthermore, the messages ${\partial }^{s,\ell }({\tilde{X}}^{\ell })$ with $\left\{(r_{\ast }^{\ell },b_{\ast }^{\ell }(\tilde{x}))\right\}$ sent from “${\tilde{\mathit{X}}}^{\ell }$” to the factor node “$v_v^{\ell }$” for the next outer loop iteration are calculated by ${\partial }_v^{\ell }({\tilde{\mathit{X}}}_{t|t-1}^{\ell })\propto {\chi }^{\ell }({\tilde{\mathit{X}}}_{t|t-1}^{\ell }){\displaystyle {\prod }_{v^{\prime }\in \mathcal{V}\backslash v}{\gamma }_v^{\ell }({\tilde{\mathit{X}}}_{t|t-1}^{\ell })}$.

Unlike previous methods that rely on truncating data associations for LMB approximation, this approach does not discard any information and is less prone to information loss. It is worth noting that the main advantage of the BP method used in this work is its scalability. The complexity of performed operations (35) to (37) scales as $O(V|{\mathbb{L}}_t|{\displaystyle {\prod }_{v=1}^{\mathrm{V}}{M}_{t,v}})$, where $M_{t,v}$ progressively increases with the false alarms and the number of PTs $\left|{\mathbb{L}}_t\right|$. As a result, the computational complexity of the BP scheme scales linearly with $\left|{\mathbb{L}}_t\right|$, the number of measurements per S/R pair $M_{t,v}$, and the number of S/R pairs V. Therefore, the computational complexity of the density approximation for each PT $\ell \in {\mathbb{L}}_t$ (see (37)) scales only linearly with the number of S/R pairs.

Based on the message-passing approach, an online estimation of the detection probabilities for the PTs is performed. The result of the previous outer BP iteration is utilized in the subsequent iteration. This approach ensures that delayed estimation problems are avoided. Once the beliefs of the PTs are obtained, PT tracks with existence probabilities below a threshold $r_{\tau }$ can be pruned. According to the findings presented in [27], a commonly used method involves extracting all PT tracks and detection probabilities with existence probabilities exceeding an application-specific threshold:

$${\widehat{\mathit{X}}}_t=\left\{({\widehat{x}}_t,\ell ),\{{\widehat{e}}_{t,v}^{\ell }{\}}_v:r_t^{\ell }\ge r_{\delta }\right\}$$

(40)

where the state ${\widehat{x}}_t$ and detection probability ${\widehat{e}}_{t,v}^{\ell }$ estimation could be achieved by the minimum mean-squared-error (MMSE) estimator: ${\widehat{x}}_t={\displaystyle \int x_t{p}_t^{\ell }(x)}dx$ and $P_{D_v}(x_t,\ell )={\widehat{e}}_{t,v}^{\ell }=$${\displaystyle \int e_{t,v}^{\ell }{\xi }_{t,v}(e_{t,v}^{\ell })}d{e}_{t,v}^{\ell }$. In addition, the estimated number of targets is given by ${\widehat{\mathit{N}}}_t=\left|{\widehat{\mathit{X}}}_t\right|$.

At the same time, unlike JIPDA track management including a “confirmed track”, “terminated track”, and “tentative track”, LMB employs a pruning step to delete the PTs whose existence probabilities are less than a certain value and an estimation step to extract the target states whose existence probabilities are greater than a certain value.

4.4. Implementation

This section presents a closed-form implementation for the R-MS-LMB filter based on the Gaussian-Beta (GB) proposed by Mahler in [36]. The augmented state of a labeled PT $\ell \in {\mathbb{L}}_t$ is defined as

$${\mu }_t^{\ell }({\tilde{x}}_t)=\mathcal{N}(x_t;{\overline{x}}_{t,j},P_{t,j}){\displaystyle {\prod }_{v=1}^{\mathrm{V}}\mathrm{Beta}(e_{t,v};{\alpha }_{t,v},{\beta }_{t,v})}$$

(41)

where $\mathcal{N}(\cdot ;{\overline{x}}_{t,j},P_{t,j})$ represents the Gaussian density of the target state with the covariance $P_{t,j}$ and the mean ${\overline{x}}_{t,j}$. $\mathrm{Beta}(e_{t,v};{\alpha }_{t,v},{\beta }_{t,v})$ describes the PT $\ell \in {\mathbb{L}}_t$ detection probability of the S/R pair v, where $\alpha >1$ and $\beta >1$ are the given parameters, and the mean and variance can be obtained by $\overline{e}=\alpha /(\alpha +\beta )$ and ${\sigma }_{\mathrm{B}}^2=\frac{\alpha \beta }{{(\alpha +\beta )}^2(\alpha +\beta +1)}$, respectively.

For the GB model, we give the following assumptions.

(1)

The augmented state of the PT $\ell \in {\mathbb{L}}_t^B$ for newly birthed LMB is a GB form:

$${\mu }_t^{\ell }(\tilde{x})=\mathcal{N}(x;{\overline{x}}_t,P_t){\displaystyle {\prod }_{v=1}^{\mathrm{V}}\mathrm{Beta}(e_v;{\alpha }_{t,v},{\beta }_{t,v})}$$

(42)

where ${\overline{x}}_t$, $P_t$, ${\alpha }_{t,v}$, and ${\beta }_{t,v}$ are the given parameters.

(2)

The kinematic model for the single PT is an acceleration model [3]:

$$f^1(x_t|x_{t-1})=\mathcal{N}(x_t;F_t{x}_{t-1},Q_t).$$

(43)

Here, $F_t=\left[\begin{array}{cc}{I}_2& \mathrm{T}{I}_2\\ 0_2& I_2\end{array}\right]$ denotes the state transition matrix and $Q_t={\sigma }_v^2\left[\begin{array}{cc}{\mathrm{T}}^4/4I_2& {\mathrm{T}}^3/2I_2\\ {\mathrm{T}}^3/2I_2& {\mathrm{T}}^2{I}_2\end{array}\right]$ is the process noise covariance, where T is the time scan and ${\sigma }_v$ is the noise generated by the acceleration process.

(3) Given $\mathrm{Beta}(e_{t-1,v};{\alpha }_{t-1,v},{\beta }_{t-1,v})$, the prediction distribution is as follows [57]:

$$\mathrm{Beta}(e_{t|t-1,v};{\alpha }_{t|t-1,v},{\beta }_{t|t-1,v})={\displaystyle \int f^2(e_{t|t-1,}|e_{t-1})}\mathrm{Beta}(e_{t-1};{\alpha }_{t-1,v},{\beta }_{t-1,v})d{e}_{t-1}$$

(44)

where ${\alpha }_{v,t|t-1}={\tau }_v{\alpha }_{v,t-1}, {\beta }_{v,t|t-1}={\tau }_v{\beta }_{v,t-1},$ in which ${\tau }_v=\frac{{\alpha }_{v,t-1}{\beta }_{v,t-1}-{\sigma }_{v,t-1}^2{({\alpha }_{v,t-1}+{\beta }_{v,t-1})}^2}{{\sigma }_{\mathrm{B}\_\mathrm{dsi}}^2{({\alpha }_{v,t-1}+{\beta }_{v,t-1})}^3}$ if ${({\sigma }_{\mathrm{B},v,t-1})}^2\le {\sigma }_{\mathrm{B}\_\mathrm{dsi}}^2$, and ${\tau }_v=1$ otherwise. This prediction can mitigate the risk of the beta distribution reducing to a Dirac delta function and enhance the variability.

Non-linear measurement models can be analyzed using techniques such as UKF and EKF. Before the outer BP iteration, if the augmented target state density of ${\chi }^{\ell }({\tilde{\mathit{X}}}_{t|t-1}^{\ell })$ of (34) is a GB form, the augmented target state density of the message ${\partial }^{s,\ell }({\tilde{X}}^{\ell })$ is also a GB form. However, the message ${\gamma }_v^{\ell }({\tilde{\mathit{X}}}_{t|t-1}^{\ell })$ of (36) is a multi-modal likelihood function. Consequently, the augmented target state densities of $b({\tilde{\mathit{X}}}_t^{\ell })$ of (37) and ${\partial }_v^{\ell }({\tilde{\mathit{X}}}_{t|t-1}^{\ell })$ are represented by GB mixtures (GBM). After each outer loop iteration, we use first- and second-order moment matching to convert the GBM into a single BG.

In addition, the proposed R-MS-LMB multi-target tracking method is depicted in Figure 4.

5. Numerical Study

Simulations are implemented in this section to demonstrate the performance of the proposed underwater MTT method. In what follows, we describe two simulation scenarios and present the corresponding simulation results. Assuming that acoustic signal propagation losses are predominantly dependent on the range, we employ the Fermi function to simulate the true but unknown target detection probabilities [58]:

$${P_D}^{\prime }=\{\begin{array}{cc}\frac{1}{1+{10}^{(d_{\mathrm{equal}}/d_0-1)/b}}& d_{\mathrm{equal}}\ge d_b\\ 0.2& d_{\mathrm{equal}}<d_b\end{array}$$

(45)

where b, $d_0$, $d_{\mathrm{equal}}$, and $d_b$ are the given parameters, and $d_{\mathrm{equal}}=\sqrt{d_s{d}_r}$ is the equivalent monostatic range. To demonstrate poor detection at the end-fire of AUVs, let $P_D=\cos (\varphi ){P_D}^{\prime }$ if the measurement angle $\left|\theta \right|<{\theta }_{ef}$; otherwise, $P_D={P_D}^{\prime }$.

In simulations, the clutter pdf $c(z)$ is uniform on $[-180°,180°]$ and $[0, 60 \mathrm{km}]$. Each new birth target is modeled by $\left(r^{\ell },p^{\ell }(\tilde{x})\right)$, $\forall \ell \in {\mathbb{L}}_t^B$, where $r^{\ell }=0.01$ and $p^{\ell }(\tilde{x})=\mathcal{N}(x;\overline{x},P){\displaystyle {\prod }_{v=1}^{\mathrm{V}}\mathrm{Beta}(e_v;{\alpha }_v,{\beta }_v)}$ with ${\overline{e}}_v=0.5$, $P=\mathrm{diag}({200}^{^2},{200}^{^2},{10}^{^2},{10}^{^2})$, and $\overline{x}$ is located at the actual starting position of the target. To prevent overconfidence, we define the number of outer BP iterations as one and the number of inner iterations as four [51]. The simulation parameters are reported in

Table 1. Simulation parameters.

Parameter	Value	Specification
T	20 s	Time scan
${\sigma }_v$	10−2 m/s2	Process noise
${\sigma }_d$	100 m	Rang std. deviation
${\sigma }_{\theta }$	1°	Bearing std. deviation
$PS$	0.95	Prob. of survival
$r_{\tau }$	10−3	Deleting threshold of PTs
$r_{\sigma }$	10−3	Extracting threshold of PTs
b	0.5	Fermi function Para.b
$d_0$	20 km	Fermi function Para.R0
$d_b$	1.5 km	Fermi function Para.Rb
${\theta }_{ef}$	30°	End-fire angle
$\lambda$	10	False alarm rate

. All simulations are conducted in MATLAB with an Intel i7 processor in the 2.5 to 4.8 GHz frequency range.

5.1. First Scenario

In this scenario, we consider a mobile multi-static multi-sensor MTT network consisting of three AUVs, two fixed sound sources, and three targets. Figure 5 illustrates the positions of the sound source and the tracks of the AUVs and targets over 200 time steps. The simulated target detection probabilities follow (45).

Figure 6 presents detection probability maps for all S/R pairs at time t = 1. The figure clearly illustrates that the likelihood of detecting the target is contingent on the spatial configuration between the AUV, the target, and the source. Therefore, the detection parameters are time-varying as the targets and AUVs move. In a multi-static operation, multiple sources and receivers can eliminate each other’s acoustic detection blind areas and have high coverage. Moreover, utilizing mobile AUVs as receivers can enhance the ability to detect and track multiple targets.

To verify that the R-MS-LMB filer can estimate the target detection probabilities online, Figure 7 shows the mean estimates and associated standard deviations of the detection probabilities of the three targets for different S/R pairs. These results are from 100 independent Monte Carlo (MC) runs that vary in the S/R pair measurements. It can be seen how the detection probabilities change with time. The mean estimates are generally within the range of reality, and the behavior changes are roughly in line with the actual situation. Figure 6 highlights that the proposed method can flexibly adjust the detection probabilities.

5.2. Second Scenario

To assess the performance of the R-MS-LMB filter, we consider a challenging scenario over 200 time steps, as shown in Figure 8a. We compare the proposed R-MS-LMB filter with several other filters, including the robust multi-sensor GLMB (R-MS-GLMB) filter used in [42,43], the standard multi-sensor LMB (Std-LMB) filter with true detection probabilities for the R-MS-LMB filter, the MS-LMB filter with fixed high detection probabilities, named the High-MS-LMB filter, and the MS-LMB filter with fixed low detection probabilities, named the Low-MS-LMB filter. Additionally, the R-MS-GLMB filter samples 200 components. The High-MS-LMB and Low-MS-LMB filters presented in [32] are used with fixed detection probabilities of 0.9 and 0.5, respectively. We apply the generalized sub-pattern alignment OSPA (GOSPA) [59] metrics with $\mathrm{p}=2$ and $\mathrm{c}=2 \mathrm{km}$ to assess these filters’ general performance. The GOSPA metric with the parameter $\alpha =2$ can show the localization error (LE), the false-target error (FE), and the missed-target error (ME). All metrics are averaged over 100 independent MC runs throughout the simulation experiments, and UKF implementations are used in all filters.

The results of one MC run, shown in Figure 8b, illustrate the proposed R-MS-LMB filter’s excellent detection and estimation capabilities. Figure 9 shows the average GOSPA metric results from the overall MC runs. From Figure 9, the Std-LMB filter with accurate detection probabilities has the lowest GOSPA error, LE, ME, and FE. In addition, it should be noted that the R-GLMB and R-LMB filters can adaptively learn and estimate target detection probabilities online, which leads to better performance than filters that use fixed detection probabilities. However, the R-GLMB filter estimates the average detection probability based on the nested CPHD filters. This estimation may be lower than the actual detection probability for some targets. This may lead to false alarms and inferior FE performance compared to the Std-LMB and R-LMB filters. On the one hand, the detection probability of the High-MS-LMB filter is often higher than the actual detection probabilities of most targets, which results in target loss and a more significant ME. On the other hand, the High-MS-LMB filter uses a fixed lower detection probability, resulting in poorer FE error performance. The time average GOSPA and the average runtime are detailed in

Table 2. Time average GOSPA and the average runtime.

Filters	R-MS-LMB	R-MS-GLMB	Std-MS-LMB	Low-MS-LMB	High-MS-LMB
GOSPA (m)	307.81	353.87	271.81	845.82	425.00
LE (m)	175.65	201.87	169.35	152.46	220.27
ME (m)	82.76	100.73	81.21	753.96	83.33
FE (m)	38.50	38.49	30.83	52.13	164.34
Runtime (s)	21.34	88.26	20.62	50.49	28.97

. Compared to other filters, the R-MS-LMB filter performs more similarly to the Std-MS-LMB filter and has a manageable runtime, which reflects its superior capability for underwater MTT in challenging scenarios.

Finally, we analyze the proposed method’s complexity using the average runtime over 100 MC runs. Accordingly, we use a fixed but unknown detection probability of 0.9. To evaluate the relationship between the R-MS-LMB filter and the number of S/R pairs, we utilize the four sound sources presented in Scenario 2, with one, two, three, and four AUVs, resulting in 4, 8, 12, and 16 S/R pairs, respectively. As Figure 10 shows, the R-MS-LMB filter exhibits a low computational cost. It can also be observed that the R-MS-LMB filter runtime increases linearly with the number of S/R pairs. Based on simulations, the proposed R-MS-LMB filter can provide superior performance with a moderate computational cost and is appropriate for MTT in underwater multi-static networks involving multiple AUVs.

6. Conclusions

This paper proposes a robust multi-sensor LMB MTT method for a multi-static network with multi-AUV. In this work, we mainly focus on overcoming unknown target detection probabilities with a multi-sensor LMB filter and how to implement the filter efficiently. First, we augment the target state by treating the detection probabilities as variables rather than known parameters, allowing for the adaptation of unknown and time-varying detection probabilities. Then, we derive an approximate multi-sensor LMB filter by KL divergence and construct the factorization of the joint posteriori distribution on the basis of the derived form. Moreover, we detail an efficient implementation of the robust multi-sensor LMB filter using the belief propagation scheme. Simulation experiments show that the proposed method can simultaneously estimate the targets’ state and detection probabilities, resulting in better MTT performance than other methods. In addition, the results also show that the complexity of the proposed method scales linearly with the number of S/R pairs, demonstrating superior scalability and low computational complexity.

This work assumes that the locations of AUVs and sources are known. Future work will focus on simultaneous localization and tracking to use collected measurements to improve the location accuracy of AUVs and sources while performing MTT.