SMC-CPHD Filter with Adaptive Survival Probability for Multiple Frequency Tracking

Abstract

We propose a sequential Monte Carlo-based cardinalized probability hypothesis density (SMC-CPHD) filter with adaptive survival probability for multiple frequency tracking to enhance the tracking performance. The survival probability of the particles in the filter is adjusted using the pre-designed exponential function related to the distribution of the estimated particle points. In order to ensure whether the proposed survival probability affects the stability of the filter, the error bounds in the prediction process are analyzed. Moreover, an inverse covariance intersection-based compensation method is added to enhance cardinality tracking performance by integrating two types of cardinality information from the CPHD filter and data clustering process. To evaluate the proposed method’s performance, MATLAB-based simulations are performed. As a result, the tracking performance of the multiple frequencies has been confirmed, and the accuracy of cardinality estimates are improved compared to the existing filters.

1. Introduction

Research on multiple target tracking filters has been actively conducted in various fields as well as in the target tracking field [1,2,3,4,5,6,7,8,9,10,11,12,13]. In particular, multiple target tracking filters based on a random finite set (RFS) can be used in many applications because the filters can track targets, even if the number of targets is unknown or if the number of targets changes over time [14,15,16,17,18,19,20,21]. Using these characteristics, we have recently conducted research to estimate and to finally mitigate multiple frequencies using the cardinalized probability hypothesis density (CPHD) filter among the RFS-based filters when a global navigation satellite system (GNSS) signal is received with multiple interference signals [22]. Studies have also been conducted to enhance the cardinality estimate performance of the CPHD filter, which results in increasing the performance of frequency estimates [22]. However, the performance is still not perfect, and additional methods are needed to enhance the estimate performance of the filter and the cardinality estimate performance in real situations, especially when clutter is present. As part of these studies, in this paper, the adaptive survival probability is applied to the sequential Monte Carlo-based CPHD (SMC-CPHD) filter.

In the initial version of the RFS filters, most of the filters assumed that the probability of survival is constant [7,8,10,23]. However, in a case with clutter in the measurements, fixing the probability of survival to a constant value (relatively large value) prevents the survival intensity from properly reflecting the situation and affects the estimate performance of the filter. Therefore, it is necessary to adaptively apply the probability of survival according to the measurements including clutter within the specified boundary. To address this problem, the previous paper also suggested adaptively setting the probability of survival to the GM-PHD filter [24]. Here, the authors of [24] defined a simple survival model to adaptively compute the survival probability of the state variable according to moving direction and position. In the paper [25], only the probability of detection was adaptively set. The proposed methods in previous papers were limited in the applicable systems or performed probability adjustment based on the observation area.

Based on the existing research results, we propose an adjustment technique of adaptive survival probability to set the appropriate value of survival probability for the multiple frequency tracking system, where the probability of survival is changed to be inversely proportional to the exponential function according to the Euclidean distance between the posterior state and the particles. In order to evaluate whether or not the modified survival probability affects the stability of the filter, the error bounds in the prediction process are analyzed to verify the stability of the filter. In addition to adaptive survival probability, an inverse covariance intersection (ICI)-based compensation technique is added to improve cardinality tracking performance. The ICI method [26,27] is one of the latest methods in the sensor fusion field. The ICI method provides more precise fusion results than those obtained by conventional covariance intersection (CI) methods [27]. The compensation method with ICI fused two types of cardinality information using the ICI method, which were the cardinality information estimated from the CPHD filter and the cardinality information additionally generated to compensate for the original cardinality estimate.

The paper is organized as follows. The processing steps of the SMC-CPHD filter are reviewed in Section 2. The probability of survival points in the SMC-CPHD filter is proposed in Section 3. In addition, the convergence of the mean square estimate is proved in case the proposed survival probability is applied in the existing SMC-CPHD filter. Subsequently, in Section 4, cardinality compensation with the ICI method is additionally applied in the proposed filter structure to enhance the estimate performance of the cardinality. In Section 5, various simulations are performed to evaluate the proposed method’s performance, and it is confirmed that the proposed method has better performance compared with conventional filtering methods. Section 6 summarizes the conclusions.

2. Processing Steps of the SMC-CPHD Filter

The CPHD filter is used for tracking multiple frequencies and is designed with an SMC structure [10,23]. This section presents a brief review of the SMC implementation in the CPHD filter with the target birth density driven by measurements [23].

Suppose at time $k$ there are $n_k$ targets with states $x_{k,1},\dots ,x_{k,n_k}$ in a state space, $X\subseteq {ℝ}^{n_x}$, and $m_k$ measurements, $z_{k,1},\dots ,z_{k,m_k}$, in an observation space, $Z\subseteq {ℝ}^{n_z}$. Then, a multiple target state and a multiple target observation are, respectively, defined as the finite sets: $X_k=\left\{x_{k,1},\dots ,x_{k,n_k}\right\}\in F\left(X\right)$ and $Z_k=\left\{z_{k,1},\dots ,z_{k,m_k}\right\}\in F\left(Z\right)$. The state vector, $x=\left(y,\beta \right)$, consists of the usual kinematic components denoted by $y$ and a label $\beta$, which distinguishes between a newborn target ($\beta =1$) and a persistent target ($\beta =0$) [23]. The PHD intensity function, $D_{k|k}\left(x|Z_{1:k}\right)\triangleq D_{k|k}\left(x\right)$, is referred to as ${\displaystyle \int {\delta }_x\left(x\right)f_{k|k}}\left(x|Z_{1:k}\right)\delta X$. ${\rho }_{k|k}\left(n|Z_{1:n}\right)\triangleq {\rho }_{k|k}\left(n\right)$ is the cardinality distribution at time $k$ and ${\rho }_{\Gamma ,k|k-1}\left(n\right)\triangleq {\rho }_{\Gamma }\left(n\right)$ is the cardinality distribution of new targets at time $k$. The CPHD filter jointly propagates $D_{k|k}\left(x\right)$ and ${\rho }_{k|k}\left(n\right)$. In this Algorithm 1, as the adaptive survival probability is not applied, the assumption is that the survival probability is constant, $p_S\left(x\right)=p_S$.

At time $k=0$, it is assumed that ${\rho }_{0|0}\left(n\right)=1$ if $n=0$, and ${\rho }_{0|0}\left(n\right)=0$ if $n=1,2,\dots ,n_{\max }$ ($n_{\max }$ denotes the targets’ maximum anticipated number), and $D_{0|0}\left(y,0\right)=D_{0|0}\left(y,1\right)=0$. The intensity function, $D_{k|k}\left(x\right)$, is approximated by using the random particles. At time $k-1$:

$$D_{k-1|k-1}\left(y,0\right)\approx {\displaystyle \sum _{n=1}^{N_{k-1}^p}{w}_{k-1,p}^{\left(n\right)}{\delta }_{y_{k-1,p}^{\left(n\right)}}\left(y\right)}$$

(1)

$$D_{k-1|k-1}\left(y,1\right)\approx {\displaystyle \sum _{n=1}^{N_{k-1}^b}{w}_{k-1,b}^{\left(n\right)}{\delta }_{y_{k-1,b}^{\left(n\right)}}\left(y\right)}$$

(2)

where ${\delta }_{y_0^{}}\left(y\right)$ represents the Dirac delta function [8] ${\left\{\left(w_{k-1,p}^{\left(n\right)},y_{k-1,p}^{\left(n\right)}\right)\right\}}_{n=1}^{N_{k-1}^p}$ and ${\left\{\left(w_{k-1,b}^{\left(n\right)},y_{k-1,b}^{\left(n\right)}\right)\right\}}_{n=1}^{N_{k-1}^b}$ denote the weighted particle sets for persistent targets and newborn targets. $N_{k-1}^p$ and $N_{k-1}^b$ define the number of persistent target particles and newborn target particles.

The predicted cardinality distribution, ${\rho }_{k|k-1}\left(n\right)$, can be expressed by a convolution [10]:

$${\rho }_{k|k-1}\left(n\right)={\displaystyle \sum _{j=0}^n{\rho }_{S,k|k-1}\left(j\right)}{\rho }_{\Gamma }\left(n-j\right)$$

(3)

$${\rho }_{S,k|k-1}\left(j\right)={\displaystyle \sum _{l=j}^{\infty }\frac{l!}{j!\left(l-j\right)!}{p}_S^j{\left(1-p_S\right)}^{l-j}{\rho }_{k-1|k-1}\left(l\right)}$$

(4)

For the predicted intensity function, the summation of $D_{k-1|k-1}\left(y,0\right)$ and $D_{k-1|k-1}\left(y,1\right)$ is preferentially conducted and can be presented by the union of two particle sets as follows:

$${\left\{\left(w_{k-1}^{\left(n\right)},y_{k-1}^{\left(n\right)}\right)\right\}}_{n=1}^{N_{k-1}^{}}={\left\{\left(w_{k-1,p}^{\left(n\right)},y_{k-1,p}^{\left(n\right)}\right)\right\}}_{n=1}^{N_{k-1}^p}\cup {\left\{\left(w_{k-1,b}^{\left(n\right)},y_{k-1,b}^{\left(n\right)}\right)\right\}}_{n=1}^{N_{k-1}^b}$$

(5)

The predicted intensity function, $D_{k|k-1}\left(y,0\right)$, can be approximately calculated by the particle set:

$$D_{k|k-1}\left(y,0\right)\approx {\displaystyle \sum _{n=1}^{N_{k-1}^{}}{w}_{k|k-1,p}^{\left(n\right)}{\delta }_{y_{k|k-1,p}^{\left(n\right)}}\left(y\right)}$$

(6)

$$y_{k|k-1,p}^{\left(n\right)}~q_k\left(\cdot |y_{k-1}^{\left(n\right)},Z_k\right)$$

(7)

$$w_{k|k-1,p}^{\left(n\right)}=\frac{p_S\left(y_{k-1}^{\left(n\right)}\right){\pi }_{k|k-1}\left(y_{k|k-1,p}^{\left(n\right)}|y_{k-1}^{\left(n\right)}\right)w_{k-1}^{\left(n\right)}}{q_k\left(y_{k|k-1,p}^{\left(n\right)}|y_{k-1}^{\left(n\right)},Z_k\right)}$$

(8)

For simplicity, the importance density is adopted as $q_k\left(\cdot |y_{k-1}^{\left(n\right)},Z_k\right)={\pi }_{k|k-1}\left(\cdot |y_{k-1}^{\left(n\right)}\right)$ [25]. For $\beta =1$, the predicted intensity function, $D_{k|k-1}\left(y,1\right)$, is approximated as:

$$D_{k|k-1}\left(y,1\right)\approx {\displaystyle \sum _{n=1}^{N_{k-1}^b}{w}_{k|k-1,b}^{\left(n\right)}{\delta }_{y_{k|k-1,b}^{\left(n\right)}}\left(y\right)}$$

(9)

$$y_{k|k-1,b}^{\left(n\right)}~b_k\left(\cdot |z\right)$$

(10)

$$w_{k|k-1,b}^{\left(n\right)}=\frac{v_{k|k-1}^b}{N_k^b}$$

(11)

where $b_k\left(\cdot |z\right)$ denotes the newborn target density and depends on $g_k\left(z|y\right)$ and prior knowledge. $N_k^b=M_b\cdot m_k$ defines the total number of newborn target particles and $M_b$ refers to the number of particles per newborn target.

The updated cardinality distribution, ${\rho }_{k|k}\left(n\right)$, is calculated by:

$${\rho }_{k|k}\left(n\right)=\frac{{\mathsf{\Upsilon }}_k^0\left[D_{k|k-1};Z_k\right]\left(n\right){\rho }_{k|k-1}\left(n\right)}{\langle {\mathsf{\Upsilon }}_k^0\left[D_{k|k-1};Z_k\right],{\rho }_{k|k-1}\rangle }$$

(12)

$$\begin{array}{ll}{\mathsf{\Upsilon }}_k^u\left[D_{k|k-1},Z\right]\left(n\right)=& {\displaystyle \sum _{j=0}^{\min \left(\left|Z\right|,n\right)}\left(\left|Z\right|-j\right)!{\rho }_{\kappa ,k}\left(\left|Z\right|-j\right)P_{j+u}^n}\\ & \times \frac{{\langle 1-p_D(\cdot ),D_{k|k-1}\left(\cdot ,0\right)\rangle }^{u-\left(j+u\right)}}{{\langle 1,D_{k|k-1}\left(\cdot ,1\right)+D_{k|k-1}\left(\cdot ,0\right)\rangle }^n}\\ & \times e_j\left({\Xi }_k\left(D_{k|k-1},Z\right)\right)\end{array}$$

(13)

$${\Xi }_k\left(D_{k|k-1},Z\right)=\left\{\frac{\langle 1,{\kappa }_k\rangle }{{\kappa }_k\left(z\right)}\langle D_{k|k-1}\left(\cdot ,1\right)+p_D(\cdot )D_{k|k-1}\left(\cdot ,0\right),g\left(z|\cdot \right)\rangle :z\in Z\right\}$$

(14)

$$e_j\left(Z\right)={\displaystyle \sum _{W\subseteq Z,\left|W\right|=j}\left({\displaystyle \prod _{\zeta \in W}\zeta }\right)}$$

(15)

where $\langle g,f\rangle ={\displaystyle \int f\left(x\right)g\left(x\right)dx}$, ${\kappa }_k\left(z\right)$ defines the PHD of clutter, $p_{D,k}\left(x\right)\triangleq p_D\left(x\right)$ refers to the detection probability, and $g_k\left(z|x\right)$ represents the measurement likelihood at time $k$. $e_j\left(Z\right)$ denotes the elementary symmetric function of order $j$ for a finite set, $Z$ [8].

The updated intensity function, $D_{k|k}\left(y,0\right)$ and $D_{k|k}\left(y,1\right)$, can be respectively calculated by:

$$D_{k|k}\left(y,0\right)\approx {\displaystyle \sum _{n=1}^{N_{k-1}^{}}{w}_{k|k,p}^{\left(n\right)}{\delta }_{y_{k|k-1,p}^{\left(n\right)}}\left(y\right)}$$

(16)

$$w_{k|k,p}^{\left(n\right)}=\left(1-p_D\left(y_{k|k-1,p}^{\left(n\right)}\right)\right)w_{k|k-1,p}^{\left(n\right)}+{\displaystyle \sum _{z\in Z_k}\frac{p_D\left(y_{k|k-1,p}^{\left(n\right)}\right)g_k\left(z|y_{k|k-1,p}^{\left(n\right)}\right)w_{k|k-1,p}^{\left(n\right)}}{\theta \left(z\right)}}$$

(17)

$$D_{k|k}\left(y,1\right)\approx {\displaystyle \sum _{n=1}^{N_k^b}{w}_{k|k,b}^{\left(n\right)}{\delta }_{y_{k|k-1,b}^{\left(n\right)}}\left(y\right)}$$

(18)

$$w_{k|k,b}^{\left(n\right)}={\displaystyle \sum _{z\in Z_k}\frac{w_{k|k-1,b}^{\left(n\right)}}{\theta \left(z\right)}}$$

(19)

$$\theta \left(z\right)={\kappa }_k\left(z\right)+{\displaystyle \sum _{n=1}^{N_k^b}{w}_{k|k-1,b}^{\left(n\right)}}+{\displaystyle \sum _{n=1}^{N_{k-1}^{}}{p}_D\left(y_{k|k-1,p}^{\left(n\right)}\right)g_k\left(z|y_{k|k-1,p}^{\left(n\right)}\right)w_{k|k-1,p}^{\left(n\right)}}$$

(20)

For resampling steps, the estimated number of persistent targets, ${\widehat{v}}_k^p={\displaystyle {\sum }_{n=1}^{N_{k-1}}{w}_{k|k,p}^{\left(n\right)}}$, and the expected number of newborn targets, ${\widehat{v}}_k^b={\displaystyle {\sum }_{n=1}^{N_k^b}{w}_{k|k,b}^{\left(n\right)}}$, are respectively calculated. The number of particles, $N_k^p$, is selected by $\left[M_p{\widehat{v}}_k^p\right]$, where $[\cdot ]$ denotes the nearest integer and $M_p$ is the number of particles per persistent targets. Two intensity functions of (16) and (18) are respectively resampled $N_k^p$ and $N_k^b$ times to eliminate the particles with small weights and multiply the particles with large weights.

After resampling steps, the intensity function, $D_{k|k}\left(y,0\right)$ and $D_{k|k}\left(y,1\right)$, can be respectively approximated by:

$$D_{k|k}\left(y,0\right)\approx {\displaystyle \sum _{n=1}^{N_k^p}{w}_{k,p}^{\left(n\right)}{\delta }_{y_{k,p}^{\left(n\right)}}\left(y\right)}$$

(21)

$$D_{k|k}\left(y,1\right)\approx {\displaystyle \sum _{n=1}^{N_k^b}{w}_{k,b}^{\left(n\right)}{\delta }_{y_{k,b}^{\left(n\right)}}\left(y\right)}$$

(22)

where $w_{k,p}^{\left(n\right)}={\widehat{v}}_k^p/N_k^p$ and $w_{k,b}^{\left(n\right)}={\widehat{v}}_k^b/N_k^b$.

The whole process of the SMC-CPHD filter is presented in detail in [23], and the SMC-CPHD filter’s pseudocode at time $k$ is as shown below [23]:

Algorithm 1: SMC-CPHD filter

Inputs: cardinality distribution: ${\rho }_{k-1|k-1}\left(n\right),$

intensity function (case $\beta =0): D_{k-1|k-1}\left(y,0\right): {\left\{\left(w_{k-1,p}^{\left(n\right)},y_{k-1,p}^{\left(n\right)}\right)\right\}}_{n=1}^{N_{k-1}^p}$, Equation (1)

intensity function (case $\beta =1): D_{k-1|k-1}\left(y,1\right): {\left\{\left(w_{k-1,b}^{\left(n\right)},y_{k-1,b}^{\left(n\right)}\right)\right\}}_{n=1}^{N_{k-1}^b}$, Equation (2)

measurement set: $Z_k=\left\{z_{k,1},\dots ,z_{k,m_k}\right\} (m_k=\left|Z_k\right|)$

predicted cardinality distribution: ${\rho }_{k|k-1}\left(n\right)$, Equation (3)

input particle sets: ${\left\{\left(w_{k-1}^{\left(n\right)},y_{k-1}^{\left(n\right)}\right)\right\}}_{n=1}^{N_{k-1}^{}}$, Equation (5)

predicted intensity function (case $\beta =0): D_{k|k-1}\left(y,0\right)$, Equation (6)

for$n=1,\dots ,N_{k-1}$do

$y_{k|k-1,p}^{\left(n\right)}~{\pi }_{k|k-1}\left(\cdot |y_{k-1}^{\left(n\right)}\right)$, Equation (7)

$w_{k|k-1,p}^{\left(n\right)}=p_S\left(y_{k-1}^{\left(n\right)}\right)w_{k-1}^{\left(n\right)}$, Equation (8)

end for

number of newborn target particles: $N_k^b=M_b\cdot m_k$

predicted intensity function (case $\beta =1): D_{k|k-1}\left(y,1\right)$, Equation (9)

for$j=1,\dots ,m_k$do

for $l=1,\dots ,M_b$ do

$n=l+\left(j-1\right)M_b$

$y_{k|k-1,b}^{\left(n\right)}~b_k\left(\cdot |z_{k,j}^{}\right)$, Equation (10)

$w_{k|k-1,b}^{\left(n\right)}=v_{k|k-1}^b/N_k^b$, Equation (11)

end for

end for

updated cardinality distribution: ${\rho }_{k|k}\left(n\right)$, Equation (12)

for$\mathrm{all} \mathbf{z}\in Z_k$do

${\Xi }_k\left(D_{k|k-1},Z_k\backslash \left\{z\right\}\right)$, Equation (14)

$e_j\left({\Xi }_k\left(D_{k|k-1},Z_k\backslash \left\{z\right\}\right)\right)$, Equations (14) and (15)

${\mathsf{\Upsilon }}_k^u\left[D_{k|k-1},Z_k\backslash \left\{z\right\}\right]\left(n\right), \text{Equation} \left(13\right) (\mathrm{for} u=1,0)$

end for

updated intensity function: $D_{k|k}\left(y,0\right), D_{k|k}\left(y,1\right)$, Equations (16) and (18)

updated weights (case $\beta =0): w_{k|k,p}^{\left(n\right)}, n=1,\dots ,N_{k-1}^{}$, Equation (17)

updated weights (case $\beta =1) : w_{k|k,b}^{\left(n\right)}, n=1,\dots ,N_k^b$, Equation (19)

estimated number of persistent and newborn targets: ${\widehat{v}}_k^p={\displaystyle {\sum }_{n=1}^{N_{k-1}}{w}_{k|k,p}^{\left(n\right)}}, {\widehat{v}}_k^b={\displaystyle {\sum }_{n=1}^{N_k^b}{w}_{k|k,b}^{\left(n\right)}}$

number of particles for resampling: $N_k^p=\left[M_p{\widehat{v}}_k^p\right]$

$\mathrm{resample} N_k^p \mathrm{times} \mathrm{from} {\left\{\left(w_{k|k,p}^{\left(n\right)}/{\widehat{v}}_k^p,y_{k|k-1,p}^{\left(n\right)}\right)\right\}}_{n=1}^{N_{k-1}^{}} \mathrm{to} \mathrm{obtain} {\left\{\left(w_{k,p}^{\left(n\right)},y_{k,p}^{\left(n\right)}\right)\right\}}_{n=1}^{N_k^p}$

$\mathrm{resample} N_k^b \mathrm{times} \mathrm{from} {\left\{\left(w_{k|k,b}^{\left(n\right)}/{\widehat{v}}_k^b,y_{k|k-1,b}^{\left(n\right)}\right)\right\}}_{n=1}^{N_k^b} \mathrm{to} \mathrm{obtain} {\left\{\left(w_{k,b}^{\left(n\right)},y_{k,b}^{\left(n\right)}\right)\right\}}_{n=1}^{N_k^b}$

Outputs: updated cardinality distribution: ${\rho }_{k|k}\left(n\right),$

updated intensity function: $D_{k|k}\left(y,0\right): {\left\{\left(w_{k,p}^{\left(n\right)},y_{k,p}^{\left(n\right)}\right)\right\}}_{n=1}^{N_k^p}$, Equation (21)

3. SMC-CPHD Filter with Adaptive Survival Probability

The processing steps of the SMC-CPHD filter [23] were reviewed in the previous section. In this section, how to adaptively apply the survival probability is explained, and the result that changing the survival probability does not affect the stability of the filter is demonstrated through convergence analysis.

3.1. New Adaptive Survival Probability

As mentioned previously, the existing versions of the SMC-CPHD filter employ a predesigned value for the survival probability, which remains identical for all targets independent of distances between their posterior estimate and particles [7,8,10,23]. Here, we describe a simple survival model shown by $p_{S,k|k-1}\left(y_{k-1}^i\right)\triangleq p_S\left(y_{k-1}^i\right)$ to adaptively compute the survival probability of the state variable according to the distance from each particle, $i=1,2,\cdots ,N_{k-1}$.

In the proposed algorithm, the survival model includes a set of the boundary with their corresponding posterior estimate. The probability illustrates the survival probability with respect to the posterior estimate, which is calculated by:

$$p_S\left(y_{k-1}^i\right)=\{\begin{array}{l}{p}_{S,\max }\\ \left(if{d}_{k-1}^i=0\right)\\ p_{S,\max }\exp \left(-d_{k-1}^i/d_{\max }\right)\\ \left(if0<d_{k-1}^i<d_{\max }and{d}_M^2\le T_g\right)\\ p_{S,\min }=0.3679\times p_{S,\max }\\ \left(if{d}_{k-1}^i=d_{\max }or{d}_M^2>T_g\right)\end{array}$$

(23)

where $p_S$ represents the adaptive survival probability at time $k$, $p_{S,max}$ represents the maximum value of $p_S$, and $\exp (\cdot )$ is the exponential function. $p_{S,\min }$ indicates the minimum value of $p_S$ and is equal to $0.3679\times p_{S,max}$ because $\exp (-1)=0.3679$ in the case of $d_{k-1}^i=d_{\max }$, with $d_{\max }$ being the maximum distance of a particle to the prior. $d_{k-1}^i=\left|y_{k-1}^i-{\widehat{y}}_{k-1}^{}\right|$ denotes the Euclidean distance between the posterior estimate, ${\widehat{y}}_{k-1}^{}$, and the $i$-th of each particle $y_{k-1}^i$ at time $k$, $d_{max}$ is the maximum value of the distance, $d_{k-1}^i$, and ${\widehat{y}}_{k-1}^{}$ is the observation location. $d_M^2\triangleq {\left({\epsilon }_k^{\left(ij\right)}\right)}^T{\left(P_{z,k|k-1}^{\left(j\right)}\right)}^{-1}{\epsilon }_k^{\left(ij\right)}$ is the Mahalanobis distance associated with ${\epsilon }_k^{\left(ij\right)}$ and $d_M^2\le T_g$ represents the validation gate area, which is selected to reduce false but apparently possible target measurements in the presence of the clutter or noise [28,29,30,31,32]. The other related notations are as follows: ${\epsilon }_k^{\left(ij\right)}\triangleq z_k^{\left(i\right)}-h\left(y_{k|k-1}^{\left(j\right)}\right)$ denotes the innovation vector associated with the $i$-th measurement, $z_k^{\left(i\right)}\in Z_k$, and with respect to that, the $j$-th predicted observation, $P_{\epsilon ,k|k-1}^{\left(j\right)}=P_{z,k|k-1}^{\left(j\right)}=H{P}_{y,k|k-1}^{\left(j\right)}{H}^T+R$ is the residual covariance matrix of ${\epsilon }_k^{\left(ij\right)}$, $R$ defines the measurement noise covariance matrix, $h\left(y_{k|k-1}^{\left(j\right)}\right)$ refers to the measurement model, and $T_g$ is the threshold of the $n_y$-dimensional validation gate and a design parameter.

Figure 1 shows the adaptive survival probability based on (23).

As shown in Figure 1, the closer the distance, $d_{k-1}^i$, the higher the probability of survival. As mentioned previously and shown in (23), the adaptive survival probability, $p_S$, does not consist of constant values but varying values according to the distance within a validation gate area, and $p_{S,\min }$ is not zero but a value determining by $p_{S,\max }$. The parameter $p_{S,\max }$ is a design parameter and is set to 1 in our simulations. In the setting of the value of $p_{S,\max }$, the value of 0.98 is often used for the fixed survival probability [23], while the value of one is used for consideration of the attenuation of the survival probability according to the distance [24].

However, an improvement of the proposed method is still needed. When we apply the calculated value $d_{\max }$ to set the survival probability, the sample impoverishment can occur in the worst case because the survival probability is directly affected by the distance. Therefore, an additional method such as distance leveling for obtaining enough samples may be necessary. Nevertheless, the proposed method is meaningful because the effect of resampling in the SMC-CPHD filter can be present; $p_S$ has an effect on the weights of the particles because $p_S$ is multiplied to the weights of the particles. Simulations for verifying the performance improvement of the proposed algorithm are shown in Section 5.

3.2. Convergence Analysis

The bounds for the mean square error (MSE) [33,34] of the SMC-PHD filter depend on some considerations related to the designed parameters of the model of the target system. In order to verify the convergence of the proposed algorithm, those analysis results are extended and applied to the proposed algorithm. In particular, the effect of the character that changed the survival probability in Section 3.1 on the convergence of the filter is analyzed. Some of the assumptions and considerations used in the convergence analysis are given as follows.

If ${\theta }^N$ defines a sequence of measures that rely on the number of particles, $N$, then ${\theta }^N$ converges to $\theta$ when $\forall \phi \in B\left(R^d\right)$.

$$\underset{N\to \infty }{\lim }E\left[{\left(\langle {\theta }^N,\phi \rangle -\langle \theta ,\phi \rangle \right)}^2\right]=0$$

(24)

where $B\left(R^d\right)$ denotes the set of bounded Borel measurable function on $R^d$ and $d$ refers to the dimension of the space. When the measure in the inner product, $\langle .,.\rangle$, is discrete, it is the summation inner product:

$$\langle D_{k|k}^N,\phi \rangle ={\displaystyle \sum _{i=1}^N{w}_k^{\left(i\right)}}\phi \left(x_k^{\left(i\right)}\right)$$

(25)

In order to prove convergence of the proposed algorithm, ${D^{\prime }}_{k|k-1}$ is defined as the density propagated from the previous time step, and the boundary of $\left|\langle D_{k|k}^N,\phi \rangle -\langle {D^{\prime }}_{k|k-1},\phi \rangle \right|$ is analyzed by the triangle inequality:

$$\begin{array}{l}\left|\langle D_{k|k-1}^N,\phi \rangle -\langle {D^{\prime }}_{k|k-1},\phi \rangle \right|\le \left|\langle D_{k|k-1}^N,\phi \rangle -\langle D_{k-1|k-1}^N,{\pi }_{k|k-1}\phi \rangle \right|\\ +\left|\langle D_{k-1|k-1}^N,{\pi }_{k|k-1}\phi \rangle -\langle D_{k-1|k-1}^{},{\pi }_{k|k-1}\phi \rangle \right|\end{array}$$

(26)

where ${\pi }_{k|k-1}$ refers to the transition density and is bounded by the system transition, $f_{k|k-1}$, and the PHD of spawned targets, $s_{k|k-1}$, as follows:

$${\pi }_{k|k-1}\left(y,y_{k-1}\right)=p_S\left(y_{k-1}\right)f_{k|k-1}\left(y|y_{k-1}\right)+s_{k|k-1}\left(y|y_{k-1}\right)$$

(27)

Let ${\varsigma }_{k-1}$ be the $\sigma$-algebra generated by the particles $\left\{y_{k-1}^{\left(i\right)}\right\}$. Then, the first term of the right side of (26) can be expected as:

$$E\left[{\left(\langle D_{k|k-1}^N,\phi \rangle -\langle D_{k-1|k-1}^N,{\pi }_{k|k-1}\phi \rangle \right)}^2|{\varsigma }_{k-1}\right]=E\left[{\langle D_{k|k-1}^N,\phi \rangle }^2\right]-{\langle D_{k-1|k-1}^N,{\pi }_{k|k-1}\phi \rangle }^2$$

(28)

The above equation can be expressed in sum using the independence of each particle as follows:

$$\begin{array}{l}E\left[{\langle D_{k|k-1}^N,\phi \rangle }^2\right]-{\langle D_{k-1|k-1}^N,{\pi }_{k|k-1}\phi \rangle }^2\\ ={\left(\frac{\widehat{T}}{N}\right)}^2{\displaystyle \sum _{i=1}^N\left(E\left[{\left(\phi \left({\tilde{y}}_k^{\left(i\right)}\right)p_s\left(y_{k-1}\right)\right)}^2|{\varsigma }_{k-1}\right]-\left({\pi }_{k|k-1}\phi \right){\left(y_{k-1}^{\left(i\right)}\right)}^2\right)}\end{array}$$

(29)

where $\widehat{T}$ denotes the estimated number of targets.

According to (29), the equation’s boundary can be expressed as:

$$\left|E\left[{\langle D_{k|k-1}^N,\phi \rangle }^2\right]-{\langle D_{k-1|k-1}^N,{\pi }_{k|k-1}\phi \rangle }^2\right|\le \frac{{\widehat{T}}^2}{N}{\Vert \phi \Vert }^2\left({\Vert p_s\left(y_{k-1}\right)\Vert }^2+{\Vert {\pi }_{k|k-1}\Vert }^2\right)$$

(30)

Using the inequality of Minkowski, (30) can be rearranged as:

$$\begin{array}{l}E{\left[{\left(\langle D_{k|k-1}^N,\phi \rangle -\langle {D^{\prime }}_{k|k-1},\phi \rangle \right)}^2\right]}^{1/2}\\ \le E{\left[{\left(\langle D_{k|k-1}^N,\phi \rangle -\langle D_{k-1|k-1}^N,{\pi }_{k|k-1}\phi \rangle \right)}^2\right]}^{1/2}\\ +E{\left[{\left(\langle D_{k-1|k-1}^N,{\pi }_{k|k-1}\phi \rangle -\langle D_{k-1|k-1}^{},{\pi }_{k|k-1}\phi \rangle \right)}^2\right]}^{1/2}\end{array}$$

(31)

If there are no newly created targets at time $k$, the boundary of (31) can be finally expressed as:

$$E{\left[{\left(\langle D_{k|k-1}^N,\phi \rangle -\langle {D^{\prime }}_{k|k-1},\phi \rangle \right)}^2\right]}^{1/2}\le \frac{1}{\sqrt{N}}\Vert \phi \Vert {\left(\widehat{T}{\Vert p_s\left(x_{k-1}\right)\Vert }^2+{\Vert {\pi }_{k|k-1}\Vert }^2\right)}^{1/2}+\sqrt{c_{k-1|k-1}}$$

(32)

In (32), $c_{k-1|k-1}$ defines the MSE’s boundary at the previous time step, $k-1$, and can be obtained from the assumption that:

$$E\left[{\left(\langle D_{k-1|k-1}^N,\phi \rangle -\langle D_{k-1|k-1}^{},\phi \rangle \right)}^2\right]\le c_{k-1|k-1}\frac{{\Vert \phi \Vert }^2}{N}$$

(33)

Finally, if $\Vert {\pi }_{k|k-1}\Vert \le 1+T_{k|k-1}$, where $T_{k|k-1}$ denotes the number of spawned targets at time step $k$:

$$E\left[{\left(\langle D_{k|k-1}^N,\phi \rangle -\langle {D^{\prime }}_{k|k-1},\phi \rangle \right)}^2\right]\le \frac{{\Vert \phi \Vert }^2}{N}{\left({\left(\widehat{T}{\Vert p_s\left(y_{k-1}\right)\Vert }^2+{\Vert 1+T_{k|k-1}\Vert }^2\right)}^{1/2}+\sqrt{c_{k-1|k-1}}\right)}^2$$

(34)

In the case of the proposed algorithm, we confirmed that the MSE’s convergence depends on only the number of targets, $\widehat{T}$, including the spawned targets, $T_{k|k-1}$, the particles’ number, $N$, and the maximum boundary of the survival probability. Thus, if the maximum range of the survival probability is specified ($p_{s,\max }$), as in the previous section, the SMC-CPHD filter will converge as $\frac{{\Vert \phi \Vert }^2}{N}{\left({\left(\widehat{T}{\Vert p_{s,\max }\Vert }^2+{\Vert 1+T_{k|k-1}\Vert }^2\right)}^{1/2}+\sqrt{c_{k-1|k-1}}\right)}^2$.

4. Cardinality Compensation with ICI

Information fusion is mainly applied to combine information using various sensors in the sensor fusion field. Typical methods, such as CI, cross covariance, sequential CI, and ellipsoidal intersection methods [26,27,35,36,37,38], have been consistently studied. They are classified as ellipsoidal methods and the ICI method [26,27], which is one of the latest data fusion methods. A data fusion method based on the ICI is proposed to obtain less conservative but still consistent estimates. A major advantage of the ICI method is that its fusion results are more accurate than those of the CI method [27]. Therefore, the ICI method is used for combining two pieces of cardinality information.

According to fusion methods, the cardinality, $C_f$, can be determined as:

$$C_f=\alpha C_1+\beta C_2$$

(35)

for any $\gamma \in \left[0,1\right]$. The covariance of the fused cardinality and gains in (35) are set as $\alpha =P_{C_f}\times \left({\left(P_{C_1}\right)}^{-1}-\gamma {\left(\gamma P_{C_1}+\left(1-\gamma \right)P_{C_2}\right)}^{-1}\right)$ and $\beta =P_{C_f}\times \left({\left(P_{C_1}\right)}^{-1}-\left(1-\gamma \right){\left(\gamma P_{C_1}+\left(1-\gamma \right)P_{C_2}\right)}^{-1}\right)$. $C_1$ is the estimated cardinality using the CPHD filter, and $P_{C_1}$ is the variance of $C_1$. $C_2$ is the additionally generated cardinality by iterative clustering to compensate for the original cardinality estimate, $C_1$. The iterative clustering algorithm [39] repeatedly performs fuzzy c-means (FCM) clustering using the resampled particles obtained from the CPHD filter. The clustering number that maximizes the clustering evaluation index among the candidates of the clustering number is selected as cardinality information, $C_2$. $P_{C_2}$ defines the quality information and represents the variance of $C_2$. Finally, the covariance of the fused cardinality can be calculated by:

$${\left(P_{C_f}\right)}^{-1}={\left(P_{C_1}\right)}^{-1}+{\left(P_{C_2}\right)}^{-1}-{\left(\gamma P_{C_1}+\left(1-\gamma \right)P_{C_2}\right)}^{-1}$$

(36)

The compensation method using the ICI method does not always guarantee good performance before fusing. A certain condition is required to continuously obtain the well-fused results and is analyzed by an inequality equation. Because the cardinality is a scalar value, Equation (36) can be expressed as:

$$\frac{1}{P_{C_f}}-\frac{1}{P_{C_1}}=\frac{1}{P_{C_2}}-\frac{1}{\gamma P_{C_1}+\left(1-\gamma \right)P_{C_2}}$$

(37)

If the difference in Equation (37) is positive, the fused result using the ICI method always has better estimated results compared with those of the CPHD filter. Therefore, to always make the above equation positive, the $\gamma P_{C_1}+\left(1-\gamma \right)P_{C_2}$ term guarantees the positive condition or is equal to zero. The term is multiplied in both sides of (37) and is then rewritten as follows:

$$\begin{array}{l}\left(\frac{1}{P_{C_f}}-\frac{1}{P_{C_1}}\right)\times \left(\gamma P_{C_1}+\left(1-\gamma \right)P_{C_2}\right)\\ =\left(\frac{1}{P_{C_2}}-\frac{1}{\gamma P_{C_1}+\left(1-\gamma \right)P_{C_2}}\right)\times \left(\gamma P_{C_1}+\left(1-\gamma \right)P_{C_2}\right)=\gamma \left(\frac{P_{C_1}}{P_{C_2}}-1\right)\end{array}$$

(38)

where $\gamma \in \left[0,1\right]$. In addition, if $P_{C_1}\ge P_{C_2}$, $\gamma \left(\frac{P_{C_1}}{P_{C_2}}-1\right)\ge 0$. This means that the fused result with the ICI method has a lower estimated error compared with the original CPHD filter. Finally, the estimated cardinality is only used in the $P_{C_1}\ge P_{C_2}$ case when performing the ICI-based cardinality compensation process.

5. Simulations

Simulations were performed to verify the performance of the proposed algorithm. This paper aims to track multiple frequencies using the CPHD filter to remove GNSS interference when GNSS signals are received with multiple GNSS interference signals.

In each simulation, GNSS data were generated by a SatGen simulator. The frequencies of received GNSS interference signals were designed as the combination of a constant frequency and a linear changed frequency in the time-frequency domain. We simulated two cases, and the true values of trajectory including the cardinality of each case are shown in Figure 2 and Figure 3. In Figure 2a and Figure 3a, the solid red line indicates the true trajectory, and the gray crosses show the frequency measurement, including the clutter of the CPHD filter. In both cases, the clutter conditions were set differently, as shown in Figure 2c, to confirm how resistant the proposed method is to the clutter. The jamming to signal ratio was set to 30 dB, the power of the received GNSS signal was set to −130 dBm, and the overall noise floor was set to −114 dB/MHz. In addition, the chirp rate was set to $3.65\times {10}^{12}\mathrm{Hz}/\mathrm{s}$ [40,41]. The initial digital frequency was set to 0.25 in case 1 and was respectively set to 0 and 0.5 in case 2. The intermediate frequency was 9.548 MHz, and the sampling frequency was 38.192 MHz.

The A_CPHD was compared with the following methods: The first method is “C_CPHD”, which is a conventional CPHD filter with k-means clustering. The second method is “FCM_CPHD”, which is a conventional CPHD filter applying FCM clustering so that it is robust in the presence of heavy noise and clutter. The last method is “CC_CPHD”, which is the FCM_CPHD applying cardinality compensation [22] with the ICI method and iterative clustering, which was recently proposed.

In the first and second cases, the tracking and cardinality estimate results from each method are presented as shown in Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11. The top panel of each figure indicates the frequency tracking result, and the bottom panel of each figure shows the cardinality estimate result. In each top graph, the true trajectory is indicated by the solid red line, the measurements with the clutter are shown by the gray crosses, and the black dots present the estimation results. In each bottom graph, the estimation result of cardinality is shown in the solid black line, and the true value of cardinality is presented by the solid black line. As shown in Figure 2c, the distribution of clutter varies with time. The clutter had the Poisson distribution property. Here, the Poisson average rate value ($\lambda$) in the frequency measurement was changed over a total of 300 simulation runs.

In Figure 4, the C_CPHD was used. As shown in this figure, tracking loss partially occurs, and the cardinality error was also severe. The results of FCM_CPHD are shown in Figure 5. Through these results, we confirmed that the tracking performance of the FCM_CPHD was a little better than that of the C_CPHD due to the robustness of the clustering method when there is a lot of noise in the measurement.

However, there was still tracking loss, and estimated errors of cardinality still remained even though FCM clustering was used. Therefore, we proposed the cardinality compensation method to improve the tracking performance. According to Figure 6, the performance of the frequency tracking and the cardinality estimate were ameliorated when the CC_CPHD was used.

However, as shown in Figure 6b, the true value of cardinality is 2 when the time index is 165, but the estimated value is 3. We eventually aim to eliminate the frequency of the interference signal by using the estimation result as the notch frequency. Therefore, the problem arises of increasing the dimension of the notch filter if the number of the frequency to be removed is estimated to be more than the true number. In addition, as shown in Figure 6a, a performance improvement of the algorithm is needed to prevent this situation because this is a result of estimating the wrong frequencies, resulting in the elimination of frequency components that should not be removed. Therefore, we have proposed the adaptive survival probability as earlier described, and the results are shown in Figure 7.

The results in Figure 7 show that the estimated results of both frequency and cardinality represent close to the true value. The cardinality value differs from the true value when the time index is from 101 to 109, as seen in Figure 7b, because the distance between the two frequencies at that time is too close, meaning that the estimate accuracy is reduced.

To evaluate the tracking performance, an optimal sub-pattern assignment (OSPA) [42], which is a general performance index used in the multiple target tracking field, was used as a performance index.

Table 1. The performance of tracking (case 1).

Algorithm	OSPA	Computation Time [sec.]
C_CPHD	16.2388	2.6815
FCM_CPHD	12.7585	2.7087
CC_CPHD (without adaptation)	6.6667	4.7724
A_CPHD (Proposed)	6.5437	4.7032

shows the results from Figure 4, Figure 5, Figure 6 and Figure 7 by the OSPA index, and the proposed method with adaptive survival probability (A_CPHD) has a small OSPA index compared to the existing methods. This means that the estimate performance of frequency and cardinality was improved. The improvement of the computational efficiency was also found through a smaller calculation time than the CC_CPHD without adaptive survival probability. These results can be considered as the effect of adaptive gating due to applying the probability by the conditions.

The second simulation assumes that more complex multiple frequencies are input than the first simulation case, as shown in Figure 3. The clutter conditions are the same as in case 1. Figure 8, Figure 9, Figure 10 and Figure 11 show the results of estimating frequency and cardinality for case 2 using the C_CPHD, FCM_CPHD, CC_CPHD, and A_CPHD methods, respectively.

Table 2. The performance of tracking (case 2).

Algorithm	OSPA	Computation Time [sec.]
C_CPHD	17.2284	3.5796
FCM_CPHD	16.9158	3.7380
CC_CPHD (without adaptation)	7.1634	9.8984
A_CPHD (Proposed)	6.9683	9.4285

shows the results in Figure 8, Figure 9, Figure 10 and Figure 11 by the OSPA performance index. As shown in Figure 8, Figure 9, Figure 10 and Figure 11 and Table 2, the performance of the proposed algorithm for case 2, like case 1, is superior to that of the existing methods. However, in the second scenario, the OSPA of the FCM_CPHD is not significantly improved compared with the results of C_CPHD, unlike the first scenario, which results from a strong tendency to unite in one when clustering by FCM in cases where the frequencies are dense. In addition, for CC_CPHD and A_CPHD using the iterative clustering method, the calculation time in the first scenario is twice as much as C_CPHD and FCM_CPHD, while the second scenario is tripled. This is because the number of targets (frequencies) to estimate increased, and the number of candidates for iterative clustering consequently increased, which results from increasing the calculation time.

From the simulation results for both cases, the combination of FCM and cardinality compensation methods (CC_CPHD) shows a significant improvement in performance over existing methods. However, depending on the scenarios and the degree of the clutter, the estimated cardinality may be more than the actual cardinality, as in Figure 6b and Figure 10b. From the perspective of ultimately eliminating frequencies, as mentioned briefly previously, the mitigation performance can be degraded because the number of frequencies to be removed is increased; thus it is hard to respond to this situation. However, the addition of the proposed adaptive survival probability shows more robust performance in the clutter or scenarios.

Consequently, through the above simulation results, we confirmed that the proposed adaptive survival probability is effective in improving performance compared to constant survival probability under any clutter conditions. This was because survival probability affected the weight of the particles, resulting in an adaptive gating effect. However, a review of the resampling methods may be necessary because sample impoverishment may occur if the adaptive boundary is incorrectly set.

6. Conclusions

We proposed an SMC-CPHD filter with adaptive survival probability for multiple frequency tracking to enhance the tracking performance in the conventional SMC-CPHD filter. The survival probability of the particles in the filter was adjusted using the pre-designed exponential function. In order to confirm whether the proposed survival probability affects the stability of the filter, the error bounds in the prediction process were analyzed. Moreover, the ICI based compensation technique was added to improve cardinality tracking performance. To evaluate the proposed method’s performance, MATLAB-based simulations were performed. As a result, we confirmed that the tracking performance of the multiple frequencies and the accuracy of cardinality estimates were improved compared to the existing filters.