# A Particle Filter Track-Before-Detect Algorithm Based on Hybrid Differential Evolution

^{*}

Next Article in Journal

Previous Article in Journal

Department of Information and Communication Engineering, Harbin Engineering University, Harbin 150001, China

Author to whom correspondence should be addressed.

Academic Editor: Stefano Mariani

Received: 19 August 2015 / Revised: 24 October 2015 / Accepted: 27 October 2015 / Published: 3 November 2015

In this paper, we address the problem of detecting and tracking targets with a low signal-to-noise ratio (SNR) by exploiting hybrid differential evolution (HDE) in the particle filter track-before-detect (PF-TBD) context. Firstly, we introduce the Bayesian PF-TBD method and its weaknesses. Secondly, the HDE algorithm is regarded as a novel particle updating strategy, which is proposed to optimize the performance of the PF-TBD algorithm. Thirdly, we combine the systematic resampling approach to enhance the performance of the proposed algorithm. Then, an improved PF-TBD algorithm based on the HDE method is proposed. Experiment results indicate that the proposed method has better performance in detecting and tracking than previous algorithms when the targets have a low SNR.

The traditional target detection approach is called track-after-detect (TAD); it involves the application of a detection threshold at every scan, which leads to a loss of information, as some sensor responses my fail to pass the detection threshold. Additionally, TAD methods cannot detect the targets with a low signal-to-noise (SNR) under 10 dB [1]. An alternative approach is called track-before-detect (TBD) [2], which jointly processes the raw sensor data (infrared sensor, radar, passive sonar, etc.) without a threshold and jointly declares the presence of a target and its track. This is acceptable if the SNR is low.

A large number of TBD algorithms have been developed in the past several decades including direct maximum likelihood [3], the Hough transform [4], dynamic programming [5,6], and so on. These methods generally require discretization of the state space and are very computationally intensive. Recently, the popularity of the particle filter (PF) method has increased due to its flexibility to handle cases where the dynamic and observation models are non-linear and/or non-Gaussian [7]. Particle filter track-before-detect (PF-TBD) was first proposed by Salmond [8], in parallel with the work by Boers [9]. These algorithms have since been extended in [10,11,12], which form a single target filter for two targets. The literature [10] presents a performance comparison of three particle filters using several different particle proposal densities designed for a single target. In [11], a new model for solving the problem of detecting a single target arrival and tracking its state in a TBD context is presented. The work in [12] shows a modeling setup and algorithm for a multiple target recursive Bayesian TBD application in a radar context. Phase information is used to improve the PF-TBD algorithm in the literature [13]. In [14], the PF is replaced by a cost-reference particle filter (CRPF). However, these proposals have some disadvantages in terms of computational complexity.

In this paper, an improved PF-TBD algorithm is proposed, which exploits a novel particle updating strategy that regards the particle updating process as an optimization problem using the hybrid differential evolution (HDE) algorithm [15,16] before the resampling process. It can be drawn from this strategy that we can reduce effectively the particle impoverishment and optimize the performance of weak target detection. Simulation results are presented to demonstrate the validity of the algorithm.

The remainder of this paper is organized as follows. The TBD system model is introduced in Section 2. In Section 3, first of all, we introduce the Bayesian PF-TBD approach and its weaknesses. Then, the HDE algorithm is given, and an improved PF-TBD algorithm is proposed. We evaluate the performance of the proposed algorithm in Section 4. At last, conclusions and pointers for future research are presented in Section 5.

In this section, we will describe the models that will be used in the TBD application, the target dynamics models and the measurement models.

The target state is represented by a five-dimensional state vector ${x}_{k}={[{x}_{k}\text{\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}}{\dot{x}}_{k}\text{\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}}{y}_{k}\text{\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}}{\dot{y}}_{k}\text{\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}}{I}_{k}]}^{T}$, where $({x}_{k},{y}_{k})$, $({\dot{x}}_{k},{\dot{y}}_{k})$ and ${I}_{k}$ represent the position and velocity in the x and y dimensions and the intensity of the target at time $k$, respectively. $k$ is the time index. A constant-velocity (CV) model is widely used as an example.

The discrete-time system dynamics of this model is of the form:
where ${w}_{k}$ is the zero mean Gaussian noise process with covariance ${Q}_{CV}$ and ${F}_{CV}$ is the process transition matrix, which are defined by:
where ${T}_{CV}$ is the sampling period, ${q}_{1}$ and ${q}_{2}$ are the power spectral density of the acceleration noise and the noise in the rate of change of target return intensity, respectively.

$${x}_{k}=F{x}_{k-1}+{w}_{k}$$

$${F}_{CV}=\left[\begin{array}{ccccc}1& T& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 1& T& 0\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right]$$

$${Q}_{CV}=\left[\begin{array}{ccccc}{{q}_{1}{T}_{CV}}^{3}/3& {{q}_{1}{T}_{CV}}^{2}/2& 0& 0& 0\\ {{q}_{1}{T}_{CV}}^{2}/2& {q}_{1}{T}_{CV}& 0& 0& 0\\ 0& 0& {{q}_{1}{T}_{CV}}^{3}/3& {{q}_{1}{T}_{CV}}^{2}/2& 0\\ 0& 0& {{q}_{1}{T}_{CV}}^{2}/2& {q}_{1}{T}_{CV}& 0\\ 0& 0& 0& 0& {q}_{2}{T}_{CV}\end{array}\right]$$

The presence of a target in the data is denoted by a target existence variable ${E}_{k}$, a target existing in the data, ${E}_{k}=1$, otherwise ${E}_{k}=0$. The probability of target existence is modeled as a two-state Markov chain, the probability transition matrix of which is defined by:
where ${P}_{b}$ and ${P}_{d}$ denote the probability of target birth and death, respectively.

$$\mathrm{\Pi}=\left[\begin{array}{cc}1-{P}_{d}& {P}_{b}\\ {P}_{d}& 1-{P}_{b}\end{array}\right]$$

Generally, we use the results of the sensor signal processing, which are a series of images, and the signal is assumed to be two-dimensional images [8,9,10], each frame consisting of n cells in the x dimension and m cells in the $y$ dimension. The measurement at time $k$ is the intensity of each cell in the image, which is given by ${z}_{k}=\{{z}_{k}^{(i,j)}:i=1,\dots ,n;j=1,\dots ,m\}$. As in Chapter 11 of [7], for each cell, the measurement model can be defined as follows:
where ${h}_{k}^{(i,j)}({x}_{k})$ is a contribution from the intensity of the target in the cell $(i,j)$ and ${v}_{k}^{(i,j)}$ is the measurement noise in the cell, which is a known Gaussian distribution with zero mean, variance ${\sigma}^{2}$. In this work, we adopt the sensor point spread function, a point target of intensity ${I}_{k}$ at position $({x}_{k},{y}_{k})$; the contribution to cell $(i,j)$ can be written as:

$${z}_{k}^{(i,j)}=\{\begin{array}{cc}{h}_{k}^{(i,j)}({x}_{k})+{v}_{k}^{(i,j)}& {E}_{k}=1\\ {v}_{k}^{(i,j)}& {E}_{k}=0\end{array}$$

$${h}_{k}^{(i,j)}({x}_{k})=\frac{{\Delta}_{x}{\Delta}_{y}{I}_{k}}{2\pi {\mathrm{\Sigma}}^{2}}\mathrm{exp}\left(-\frac{{({x}_{k}-i\Delta x)}^{2}+{({y}_{k}-i{\Delta}_{y})}^{2}}{2{\mathrm{\Sigma}}^{2}}\right)$$

Here, the parameter $\mathrm{\Sigma}$ represents the extent of blurring and ${\Delta}_{x}$ and ${\Delta}_{y}$ denote the size of a cell in the $x$ and $y$ dimensions, respectively. Then, the set of complete measurements up to time $k$ is denoted by $Z=\{{z}_{1},\dots ,{z}_{k}\}$.

Under the assumptions of Gaussian measurement noise, which is independent of cell-to-cell, the likelihood function in the presence and absence of a target can be written as:

$$\text{p}({z}_{k}^{(i,j)}|{x}_{k},{E}_{k}=1)=\frac{1}{\sqrt{2\pi {\sigma}^{2}}}\mathrm{exp}\left(-\frac{{[{z}_{k}^{(i,j)}-{h}^{(i,j)}({x}_{k})]}^{2}}{2{\sigma}^{2}}\right)$$

$$\text{p}({z}_{k}^{(i,j)}|{E}_{k}=0)=\frac{1}{\sqrt{2\pi {\sigma}^{2}}}\mathrm{exp}\left(-\frac{{[{z}_{k}^{(i,j)}]}^{2}}{2{\sigma}^{2}}\right)$$

Then:

$$\text{p}({z}_{k}|{x}_{k},{E}_{k})={\displaystyle \prod _{i=1}^{n}{\displaystyle \prod _{j=1}^{m}\text{p}({z}_{k}^{(i,j)}|{x}_{k},{E}_{k})}}$$

The likelihood ratio for the cell is:
where:

$$l({z}_{k}|{x}_{k},{E}_{k})=\{\begin{array}{cc}\frac{p({z}_{k}|{x}_{k},{E}_{k}=1)}{p({z}_{k}|{E}_{k}=0)}& {E}_{k}=1\\ 1& {E}_{k}=0\end{array}$$

$$l({z}_{k}|{x}_{k},{E}_{k}=1)={\displaystyle \prod _{i=1}^{n}{\displaystyle \prod _{j=1}^{m}\mathrm{exp}\left(-\frac{{h}^{(i,j)}({x}_{k})[{h}^{(i,j)}({x}_{k})-2{z}_{k}^{(i,j)}]}{2{\sigma}^{2}}\right)}}$$

In this section, we give a derivation of the PF-based TBD algorithm and analyze the weaknesses of the PF-TBD, and then, an improved algorithm based on HDE is proposed.

The main problem of TBD can be formulated in the framework of recursive Bayesian estimation in different expansion forms. The target state is not defined when the target does not exist in the data, where ${E}_{k}=0$. Then, target detection and tracking can be implemented based on the posterior density $p({x}_{k},{E}_{k}|{Z}_{k})$ when the target does exist, where ${E}_{k}=1$.

$$\begin{array}{l}p({x}_{k},{E}_{k}=1|{Z}_{k})=\frac{p({z}_{k}|{x}_{k},{E}_{k}=1)p({x}_{k},{E}_{k}=1|{Z}_{k-1})}{p({z}_{k}|{Z}_{k-1})}\\ \text{\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}}\propto l({z}_{k}|{x}_{k},{E}_{k}=1)p({x}_{k},{E}_{k}=1|{Z}_{k-1})\end{array}$$

The predicted target state can be written in terms of the target and existence at the previous time, giving:
where ${p}_{b}({x}_{k})$ is the prior density of a target when it has appeared in the data between time $k-1$ and $k$. The algorithm proceeds by constructing two sets of particles. The first set, the birth particles, estimate $p\left({x}_{k}/{E}_{k}=1,{E}_{k-1}=0,{Z}_{k}\right)\text{\hspace{0.05em}}$, that is the case where the target did not exist in the data at time $k-1$, but does at time k. The second set, the continuing particles, estimates $p\left({x}_{k}/{E}_{k}=1,{E}_{k-1}=1,{Z}_{k}\right)\text{\hspace{0.05em}}$, which is the case where the target has continued to exist in the data k−1 to k. Both sets of particles require a proposal distribution. For the continuing particles, the target dynamics are used in Equation (4). For the birth particles, a uniform distribution can be used. Then, the former is approximated using interacting sets of particles, and the calculation of the latter is based on the particle weights in these two sets.

$$\begin{array}{l}p({x}_{k},{E}_{k}=1/{Z}_{k-1})={\displaystyle \int [p({x}_{k},{E}_{k}=1/{x}_{k-1},{E}_{k-1}=1)p({x}_{k-1},{E}_{k-1}=1/{Z}_{k-1})}\\ \text{\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}}+p({x}_{k},{E}_{k}=1,{E}_{k-1}=0/{Z}_{k-1})]d{x}_{k-1}\\ \text{\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}}={\displaystyle \int p({x}_{k}/{x}_{k-1},{E}_{k}=1,{E}_{k-1}=1)p({E}_{k}=1/{E}_{k-1}=1)}\\ \text{\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}}\cdot p({x}_{k-1},{E}_{k-1}=1/{Z}_{k-1})d{x}_{k-1}\\ \text{\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}}+\text{\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}}p({x}_{k}/{x}_{k-1},{E}_{k}=1,{E}_{k-1}=1,{Z}_{k-1})p({E}_{k}=1/{E}_{k-1}=0)\\ \text{\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}}={\displaystyle \int p({x}_{k}/{x}_{k-1},{E}_{k}=1,{E}_{k-1}=1)(1-{P}_{d})p({x}_{k-1},{E}_{k-1}=1/{Z}_{k-1})}d{x}_{k-1}\\ \text{\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}}+{p}_{b}({x}_{k}){P}_{b}\end{array}$$

Additionally, the estimate of the probability of existence can be calculated by particle weights in these two sets. State estimation can be obtained after resampling from these combined sets. There are three main resampling methods: multinomial resampling, residual resampling and systematic resampling [17].

To a certain extent, the resampling process is introduced to improve the degradation phenomenon of the particles. However, it also brought some problems. Particles that have higher weights will be reproduced many times. Additionally, those that have lower weights will fade away. Then, all particles will collapse to a small area, which will lead to the sample point sets not being sufficient after several iterations. The sample point sets are described the a posteriori probability density function. In the PF-TBD algorithm, if the particles collapse into a small area at time $k$, most of the particles, after the resampling process, are “the continuing particles”. This leads to these particles, which are produced by the state transition equation at time $k+1$, also collapsing into a small area. Moreover, if the area deviates from the true location, the target will be lost.

In this paper, to improve the particle impoverishment problem, produced by the resampling process, a novel particle updating strategy is proposed before the resampling process to overcome the problem. We regard the particle updating process as an optimization problem using the HDE algorithm to increase the diversity of particles and to avoid the particles collapsing into a small area.

The HDE algorithm based on the simulated annealing (SA) method is utilized for the search capability of the SA to enhance the convergence capability of differential evolution (DE) in later periods and to improve the robustness of the DE algorithm [15,16]. It uses the Metropolis criterion of the SA algorithm to replace the section step of the DE algorithm. It seems that the computation of the HDE step is more complicated than DE. However, the HDE algorithm can achieve better results than DE with much lower iterations [15,16]. The HDE algorithm saves time to make up for its complexity. Therefore, we choose the better optimized performance of HDE to obtain sampling particles that conform more to the true state distribution.

The HDE algorithm relies on the initial population generation, mutation, recombination and the new selection to probe the search space through iterative progress until the termination criteria are met.

Detailed steps are presented accordingly in the sequel.

Step 1: Creating the initial population.

The first step of HDE is to create the initial population samples (the number of generations is $\text{g}=0$) in n dimensional space as follows:
where $i=1,2,\dots ,NP;j=1,2,\dots ,n$, and $NP$ is the population size. ${x}_{ij}^{U}$ and ${x}_{ij}^{L}$ denote the upper and lower limits of the j-th variable in the population respectively. $ran{d}_{ij}(0,1)$ represents a uniformly-distributed random value within [0,1].

$${x}_{ij}(0)={x}_{ij}^{L}+ran{d}_{ij}(0,1)({x}_{ij}^{U}-{x}_{ij}^{L})$$

Step 2: Mutation operation.

The function of mutation in HDE is to maintain the diversity of a population. A typical HDE mutation sample formulation is:
where $g$ represents the g-th generation, and ${h}_{ij}(g)$ are the mutated vector samples. $r1\ne r2\ne i$ and $r1$, $r2$ are randomly-selected integers within $NP$, $r1,r2\in \{1,2,\dots ,NP\}$. ${F}_{m}$ is the scaling factor.

$${h}_{ij}(g)={x}_{ij}(g)+{F}_{m}\cdot ({x}_{r1j}(g)-{x}_{r2j}(g))$$

Step 3: Crossover operation.

The basic crossover process is a discrete recombination, which employs a crossover constant $C\in [0,1]$ to determine if the newly-generated individual samples are to be recombined. The expression of the crossover process is given in Equation (16):
where ${v}_{ij}(g)$ are the trial vector samples.

$${v}_{ij}(g)=\{\begin{array}{cc}{h}_{ij}(g)& rand(0,1)\le C\\ {x}_{ij}(g)& rand(0,1)>C\end{array}$$

Step 4: New selection operation.

The HDE algorithm adopts the Metropolis criterion of the SA algorithm to select whether to accept the trial vector samples. We construct the parameter about fitness function values $\Delta f=f({v}_{ij}(g))-f({x}_{ij}(g))$, and then, we decide the trial vector samples by the Metropolis criterion.

Metropolis criterion: if $\Delta f<0$, the trial vector sample is accepted. If $\Delta f>0$, we calculate an acceptance probability $P(\Delta f)=\mathrm{exp}(-\frac{\Delta f}{{T}_{SA}})$ and generate a pseudo-random number $r$. If $P(\Delta f)\ge r$, the trial vector sample is accepted; otherwise, the point is discarded.

Step 5: Cool down operation.

In this step, the cool down operation of the SA algorithm is executed. We defined ${T}_{SA}={T}_{SA}\cdot \rho $; here, ${T}_{SA}$ is an annealing temperature, and $\rho \in (0,1)$ is an annealing control parameter.

When the new population is propagated, Step 2 to Step 5 is repeated until the pre-specified temperature $T0$ is reached.

Firstly, the particle updating process is optimized by the HDE algorithm. The importance-sampling particles are regarded as the initial population of the HDE algorithm, and the corresponding weights are treated as the fitness functions of the target vectors, respectively. The HDE resampling scheme recombines the particles by using an iterative process of mutation, crossover and the simulated annealing operator. Then, a new set of diverse particles is propagated. Secondly, the systematic resampling is chosen to be executed.

The proposed algorithm exploits a hybrid estimation technique, where the state is expanded to include the existence variable ${s}_{k}={[{x}_{k}^{T}\text{\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}}{E}_{k}]}^{T}$. Given a set of particles describing the joint posterior density at time $k-1$, ${\left\{{s}_{k-1}^{n},{w}_{k-1}^{n}\right\}}_{n=1}^{M}$ with uniform weights, $N$ is the particle number. Then, the proposed algorithm contains the following steps:

- Step 1:
- Calculate the target existence variable: $[{\left\{{E}_{\text{k}}^{n}\right\}}_{n=1}^{N}]=\text{\hspace{0.05em}}[{\left\{{E}_{\text{k-1}}^{n}\right\}}_{n=1}^{N}\uff0c\prod ]$.
- Step 2:
- Create a set of particles. These are two possible situations.
- (i)
- Set the ‘birth particle’ state as a sample from the proposal density ${x}_{\text{k}}^{n}~{q}_{b}({x}_{k}/{E}_{k}^{n}=1,{E}_{k-1}^{n}=0,{z}_{k})$. The target’s location $({x}_{k}^{n},{x}_{k}^{n})$ is uniform over the sensor field-of-view. Its velocity and intensity are assumed to be uniform as follows ${q}_{b}(({\dot{x}}_{k}^{n},{\dot{y}}_{k}^{n})/{z}_{k})=U[{v}_{\mathrm{min}},{v}_{\mathrm{max}}]$, ${q}_{b}({I}_{k}/{z}_{k})=U[{I}_{\mathrm{min}},{I}_{\mathrm{max}}]$, where ${v}_{\mathrm{min}}$ and ${v}_{\mathrm{max}}$ are the minimum and the maximum of the target velocity, respectively. ${I}_{\mathrm{min}}$ and ${I}_{\mathrm{max}}$ are the minimum and the maximum of the target intensity, respectively.
- (ii)
- Set the “survival particle” state sampled from the proposal ${x}_{\text{k}}^{n}~{q}_{b}({x}_{k}/{E}_{k}^{n}=1,{E}_{k-1}^{n}=1,{z}_{k})$. That is, the “survival particle” ${x}_{k}$ is sampled according to the target state transition Equation (1).
- (iii)
- For each particle, compute the un-normalized weight ${\tilde{\text{w}}}_{\text{k}}^{n}$ by using the likelihood ratio ${\tilde{\text{w}}}_{\text{k}}^{n}=l({z}_{k}/{x}_{k}^{n},{E}_{k}^{n}=1)$.

- Step 3:
- Optimize the sample particles by HDE. The sampling particles are regarded as the initial population of the HDE algorithm, and the corresponding weights are regarded as the fitness functions. According to the HDE algorithm, the process is iterated until the optimal population is found or a pre-specified condition is reached, then we can get the optimal set of particles ${\left\{{{x}^{\prime}}_{\text{k}}^{n},{{\tilde{w}}^{\prime}}_{k}^{n}\right\}}_{n=1}^{N}$.
- Step 4:
- Normalize the weights of each particle: ${\text{w}}_{\text{k}}^{\prime n}={\tilde{\text{w}}}_{\text{k}}^{\prime n}/{\displaystyle \sum _{n=1}^{N}{\tilde{\text{w}}}_{\text{k}}^{\prime n}}=l({z}_{k}|{{x}^{\prime}}_{k}^{n},{E}_{k}^{n}=1)/{\displaystyle \sum _{n=1}^{N}\left(l({z}_{k}|{{x}^{\prime}}_{k}^{n},{E}_{k}^{n}=1)\right)}$.
- Step 5:
- Resample the particles: $[{\left\{{x}_{\text{k}}^{n},1/N\right\}}_{n=1}^{N}]=systematic-\text{re}sample[{\left\{{{x}^{\prime}}_{\text{k}}^{n},{{w}^{\prime}}_{k}^{n}\right\}}_{n=1}^{N}]$.
- Step 6:
- Output the state estimation: an estimate of the target state can be made from the set of particles resulting from the time algorithm above.

The posterior density of target existence at time k is ${\text{P}}_{k}\stackrel{\Delta}{=}P\left\{{E}_{k}=1/{Z}_{k}\right\}$, and then, the probability of existence is given by ${\widehat{\text{P}}}_{\text{k}}={\displaystyle \sum _{n=1}^{N}{E}_{k}^{n}}/N,\text{\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}}0\le {\widehat{\text{P}}}_{\text{k}}\le 1$. Additionally, the target state can be estimated by ${\widehat{\text{X}}}_{\text{k}}={\displaystyle \sum _{n=1}^{N}({x}_{k}^{n}{E}_{k}^{n}})/{\displaystyle \sum _{n=1}^{N}{E}_{k}^{n}}$.

We use the following simulated measurement data to test the performance of the proposed algorithm, which is widely used in the current study and some representative PF-TBD algorithms [2,3,4,5,6]. Meanwhile, the CV (constant velocity) and CT (coordinate turn) motion models are used for testing.

The simulated measurement data consists of a total of 30 frames, with the target appearing in the data at Frame 7 and disappearing at Frame 22. Each frame consists of an array of 20 × 20 pixels, $n=m=20$; each cell has unit dimensions ${\Delta}_{x}={\Delta}_{\text{y}}=1$. The probabilities of birth and death are set as ${P}_{b}={P}_{d}=0.05$. The standard deviation of the target spread function is $\sum =0.7$. The target initial state was given by ${x}_{0}={[4.2\text{\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}}0.45\text{\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}}7.2\text{\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}}0.25\text{\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}\hspace{0.05em}}20]}^{T}$. The birth particles proposal functions ${x}_{k}^{b}~U(0,n)$, ${y}_{k}^{b}~U(0,m)$, ${\dot{x}}_{k}^{b}~U(-1,1)$, ${\dot{y}}_{k}^{b}~U(-1,1)$ and ${I}_{k}^{b}~U(10,30)$. Set the observation noise variance ${\sigma}^{2}$ to get a different SNR, which is defined by $SNR=10\text{log}\left(\frac{{I}^{2}}{{\sigma}^{2}}\right)$. For example, if the SNR is 6 dB, the standard deviation of the background noise level in each pixel is $\sigma =10$ units. The HDE parameters were set according to [16]. ${F}_{m}=0.9$; $C=0.6$; the annealing initial temperature ${T}_{SA}=100$; the annealing parameter $\rho =0.9$; and the pre-specified temperature $T0=5$.

For the CV model, the transition matrix ${F}_{CV}$ is given by Equation (2) and ${Q}_{CV}$ by Equation (3), with ${q}_{1}=0.001$, ${q}_{2}=0.01$, ${T}_{CV}=1$ and the target detection threshold $\text{thr}=0.6$.

For the CT model, the system dynamics model function is given by Equation (1), and the covariances of noise process ${Q}_{CT}$ and the process transition matrix ${F}_{CT}$ are defined by:
where $w(k)$ is the turn rate, ${T}_{CT}$ is the sampling period and ${q}_{1}$ and ${q}_{2}$ are the power spectral density of the acceleration noise and the noise in the rate of change of target return intensity, respectively. Here, their values are set the same as in CV model. The turn angle is 4°, ${T}_{CT}=1$; and the target detection threshold $\text{thr}=0.7$.

$${F}_{CT}=\left[\begin{array}{ccccc}1& \frac{\mathrm{sin}\left(w(k){T}_{CT}\right)}{w(k)}& 0& -\frac{1-\mathrm{cos}\left(w(k){T}_{CT}\right)}{w(k)}& 0\\ 0& \mathrm{cos}\left(w(k){T}_{CT}\right)& 0& -\mathrm{sin}\left(w(k){T}_{CT}\right)& 0\\ 0& \frac{1-\mathrm{cos}\left(w(k){T}_{CT}\right)}{w(k)}& 1& \frac{\mathrm{sin}\left(w(k){T}_{CT}\right)}{w(k)}& 0\\ 0& \mathrm{sin}\left(w(k){T}_{CT}\right)& 0& \mathrm{cos}\left(w(k){T}_{CT}\right)& 0\\ 0& 0& 0& 0& 1\end{array}\right]$$

$${Q}_{CT}=\left[\begin{array}{ccccc}\frac{2\left(w\left(k\right){T}_{CT}-\mathrm{sin}\left(w\left(k\right){T}_{CT}\right)\right)\xb7{q}_{1}}{w{\left(k\right)}^{3}}& \frac{\left(1-\mathrm{cos}\left(w\left(k\right){T}_{CT}\right)\right)\xb7{q}_{1}}{w{\left(k\right)}^{2}}& 0& \frac{\left(w\left(k\right){T}_{CT}-\mathrm{sin}\left(w\left(k\right){T}_{CT}\right)\right)\xb7{q}_{1}}{w{\left(k\right)}^{2}}& 0\\ \frac{\left(1-\mathrm{cos}\left(w\left(k\right){T}_{CT}\right)\right)\xb7{q}_{1}}{w{\left(k\right)}^{2}}& {q}_{1}{T}_{CT}& \frac{-\left(w\left(k\right){T}_{CT}-\mathrm{sin}\left(w\left(k\right){T}_{CT}\right)\right)\xb7{q}_{1}}{w{\left(k\right)}^{2}}& 0& 0\\ 0& \frac{-\left(w\left(k\right){T}_{CT}-\mathrm{sin}\left(w\left(k\right){T}_{CT}\right)\right)\xb7{q}_{1}}{w{\left(k\right)}^{2}}& \frac{2\left(w\left(k\right){T}_{CT}-\mathrm{sin}\left(w\left(k\right){T}_{CT}\right)\right)\xb7{q}_{1}}{w{\left(k\right)}^{3}}& \frac{\left(1-\mathrm{cos}\left(w\left(k\right){T}_{CT}\right)\right)\xb7{q}_{1}}{w{\left(k\right)}^{2}}& 0\\ \frac{\left(w\left(k\right){T}_{CT}-\mathrm{sin}\left(w\left(k\right){T}_{CT}\right)\right)\xb7{q}_{1}}{w{\left(k\right)}^{2}}& 0& \frac{\left(1-\mathrm{cos}\left(w\left(k\right){T}_{CT}\right)\right)\xb7{q}_{1}}{w{\left(k\right)}^{2}}& {q}_{1}{T}_{CT}& 0\\ 0& 0& 0& 0& {q}_{2}{T}_{CT}\end{array}\right]$$

For convenience, in the experimental results, the proposed algorithm is marked as HDE PF-TBD; PF-TBD with the systematic resampling method is marked as SystematicR PF-TBD; and that using the multinomial resampling method is marked MultinomialR PF-TBD. We use exactly the same set of data and the same filter parameters. The number of particles used by each algorithm was set to 6000. The algorithm performance is gauged on an average over 100 Monte Carlo simulations.

Two sampled measurement images with SNR 6 dB are shown in Figure 1 in different frames: Frames 15 and 18. This shows that the target exists at Frames 15 and 18, but it is indistinguishable from the background noise. If using the TAD algorithm, this will inevitably lead to a loss of information at low SNR.

The detection performance of the proposed algorithm (HDE PF-TBD) is evaluated in simulations with the probability of existence. The average of the probability of detection with different SNRs of the three algorithms is compared in Table 1. Figure 2 shows the results in terms of the probability of existence at different frame times.

Algorithm | 9 dB | 6 dB | 3 dB | 1 dB |
---|---|---|---|---|

MultinomialR PF-TBD | 0.770 | 0.657 | 0.423 | 0.217 |

SystematicR PF-TBD | 0.774 | 0.682 | 0.498 | 0.300 |

HDE PF-TBD | 0.875 | 0.765 | 0.628 | 0.395 |

The target exists if the probability of existence is higher than 0.6, which a tested and sensible value and is plotted as a horizontal line in Figure 2. Moreover, the asterisk signs at the bottom of the figures indicate the time of the presence of the target. It can be seen from Figure 2 that for SNR = 1, 3, 6, 9 dB, three algorithms have target presence and absence declaration delays, and each algorithm follows the true existence of the target much more closely for the higher SNR cases, with performance degrading as SNR reduces. The HDE PF-TBD algorithm provides higher target detection confidence than the other two algorithms. Furthermore, the HDE PF-TBD algorithm performs very well among the four SNR values, in other words with less target presence and absence declaration delays than others. From Figure 2, we also can see that the probability of the existence of the proposed algorithm has a large increase when the SNR is 3 dB and 1 dB. For clarity, the data statistics of Table 1 have proven that the HDE PF-TBD has better detection sensitivity when the SNR is lower.

Moreover, it can be seen from Figure 3 that for SNR = 0, −1 dB, all three algorithms suffer a poor detection performance. However, compared to SystematicR PF-TBD and the MultinomialR PF-TBD, the HDE PF-TBD algorithm can detect the target at the end of the frames where the target exists. It is demonstrated that the proposed approach has a larger detectable and traceable SNR range and is capable of resolving large noise in the target. Meanwhile, it also proves that the improvements in the proposed algorithm can effectively enhance the detection sensitivity and that the proposed algorithm is more suitable for detecting weak targets with low SNR.

To measure the estimation accuracy of the proposed algorithm, we introduce the root mean squared error ($RMSE$) in the position and its mean $\stackrel{\u2014}{R}\stackrel{\u2014}{M}\stackrel{\u2014}{S}\stackrel{\u2014}{E}$ for different SNRs. Additionally, the $RMSE$ for M simulations with different frames at time $k$ is defined as:
where ${x}_{k,m}$ and ${y}_{k,m}$ are the true values of the target state and ${\widehat{x}}_{k,m}$ and ${\widehat{y}}_{k,m}$ are defined as the estimation of the target state. We only consider the frames in which the target exists, and the period length is L. The $\stackrel{\u2014}{R}\stackrel{\u2014}{M}\stackrel{\u2014}{S}\stackrel{\u2014}{E}$ can be shown as:

$$RMSE={\left[\frac{1}{M}{\displaystyle \sum _{m=1}^{M}\left({\left({\widehat{x}}_{k,m}-{x}_{k,m}\right)}^{2}+{\left({\widehat{y}}_{k,m}-{y}_{k,m}\right)}^{2}\right)}\right]}^{1/2}$$

$$\stackrel{\u2014}{R}\stackrel{\u2014}{M}\stackrel{\u2014}{S}\stackrel{\u2014}{E}=\frac{1}{L}{\displaystyle \sum _{i=1}^{L}RMS{E}_{i}}$$

The results in terms of mean RMSE for the different SNR values of the three algorithms are compared in Table 2. Figure 4 shows the results in terms of RMSE in position with corresponding frames for different SNR values.

Algorithm | 9 dB | 6 dB | 3 dB | 1 dB |
---|---|---|---|---|

MultinomialR PF-TBD | 2.1757 | 2.4715 | 2.7528 | 3.0814 |

SystematicR PF-TBD | 1.9771 | 2.3427 | 2.6222 | 3.0622 |

HDE PF-TBD | 1.3810 | 1.5971 | 1.8221 | 2.3595 |

As expected, the tracking accuracy of each algorithm degrades as the SNR reduces. Meanwhile, the proposed algorithm shows a better performance on the low SNR target than the other two algorithms in Figure 4. It is can be seen that the proposed algorithm has much lower RMSE when the target began to appear. It is demonstrated that the HDE PF-TBD algorithm can more accurately estimate the initial location of the target appearance. Moreover, from Table 2, it can be clearly seen that the estimation accuracy of the proposed algorithm is always superior to that of the other two algorithms. Although the RMSE of the proposed method approximates the other two algorithms at Frames 16 and 17 when the SNR is 1 dB, the tracking accuracy of the HDE PF-TBD algorithm is still higher than the other two algorithms in the whole frame range where the target exists. Therefore, the proposed algorithm is capable of tracking the target with low SNR, and it outperforms the SystematicR PF-TBD and the MultinomialR PF-TBD.

Similarly, we also analyze detection and tracking performance under the CT model. The average of the detection probability of detection with different SNRs of the three algorithms for 100 Monte Carlo simulations (MC) is compared in Table 3. Figure 5 shows the results in terms of the probability of existence in different frame times for 100 MC.

Algorithm | 9 dB | 6 dB | 3 dB | 1 dB |
---|---|---|---|---|

MultinomialR PF-TBD | 0.747 | 0.621 | 0.378 | 0.275 |

SystematicR PF-TBD | 0.768 | 0.648 | 0.399 | 0.284 |

HDE PF-TBD | 0.827 | 0.743 | 0.501 | 0.367 |

The target exists if the probability of existence is higher than 0.7, which is plotted as a horizontal line in Figure 5. It can be seen from Figure 5 that, for SNR =1, 3, 6, 9 dB, the three algorithms also have target presence and absence declaration delays, and each algorithm follows the true existence of the target much more closely for the higher SNR cases, with performance degrading as SNR reduces under the CT model. For the CT model, the HDE PF-TBD algorithm provides higher target detection confidence than the other two algorithms. Moreover, the HDE PFTBD algorithm has less target presence delays than the others. Meanwhile, we also can see that the probability of existence of the proposed algorithm has a large increase when the SNR is 3 dB and 1 dB. For clarity, the data statistics of Table 3 have proven that the HDE PF-TBD has better detection sensitivity when the SNR is lower.

Algorithm | 9 dB | 6 dB | 3 dB | 1 dB |
---|---|---|---|---|

MultinomialR PF-TBD | 1.910 | 2.1556 | 2.7044 | 2.909 |

SystematicR PF-TBD | 1.745 | 1.8739 | 2.6258 | 2.81 |

HDE PF-TBD | 1.433 | 1.4477 | 2.3836 | 2.3541 |

The results in terms of mean RMSE for the different SNR values of the three algorithms for the CT model are compared in Table 4. Figure 6 shows the results in terms of RMSE in position with corresponding frames for different SNR values for the CT model.

For the CT model, it can be clearly seen from Table 4 and Figure 6 that the proposed algorithm also shows a better performance on the low SNR target than the other two algorithms. Although the RMSE of the proposed method is a little higher than the other two algorithms at Frame 10 when SNR is 9 dB, the tracking accuracy of the HDE PF-TBD algorithm is still higher than other two algorithms in the whole frame range that the target exists. Therefore, as can be seen from the two aspects of detection and tracking, the proposed algorithm is also capable of tracking the target with the CT motion model and low SNR, and it outperforms SystematicR PF-TBD and MultinomialR PF-TBD.

In this paper, a novel particle filter track-before-detect algorithm based on hybrid differential evolution was developed. Firstly, we exploit a hybrid differential evolution as a novel particle updating strategy for overcoming the shortcoming of the particle impoverishment. Secondly, in order to reduce the influence of the particle degeneracy, the systematic resampling method is chosen for the proposed method. Finally, the proposed algorithm combines a more effective particle initiation method, which only places the particles in the highest intensity bins of the data. The proposed algorithm not only can make full use of the data information, but also decreases the particle degeneracy and enhances the diversity of particles. The detection and tracking of the performance advantage of the proposed algorithm have been demonstrated by two cases of a low signal-to-noise point target against a background of Gaussian noise with the CV and CT motion models, respectively. What is more, the proposed method is much more efficient than previous methods. Future work will continue to enhance the detection accuracy of the proposed algorithm and to explore the implementation of multiple targets within this framework.

This work is supported by the National Natural Science Foundation of China (Grant No. 61172159) and the Fundamental Research Funds for the Central Universities (HEUCFT1101).

The work presented here was carried out in collaboration between both authors. Chaozhu Zhang and Lin Li defined the research theme. Lin Li designed the methods and experiments, carried out the laboratory experiments, analyzed the data, interpreted the results and wrote the paper. Both authors have contributed to, seen and approved the manuscript.

The authors declare no conflict of interest.

- Hadzagic, M.; Michalska, H.; Lefebvre, E. Track-before-detect methods in tracking low-observable targets: A survey. Sens. Trans. Mag.
**2005**, 54, 374–380. [Google Scholar] - Rutten, M.; Gordon, N.; Maskell, S. Recursive track-before-detect with target amplitude fluctuations. IEE Proc. Radar Sonar Navig.
**2005**, 152, 345–352. [Google Scholar] [CrossRef] - Pohlig, S.C. An algorithm for detection of moving optical targets. IEEE Trans. Aerosp. Electron. Syst.
**1989**, 25, 56–63. [Google Scholar] [CrossRef] - Carlson, B.D.; Evans, E.D.; Wilson, S.L. Search radar detection and track with the hough transform. IEEE Trans. Aerosp. Electron. Syst.
**1994**, 30, 102–108. [Google Scholar] [CrossRef] - Huang, D.; Xue, A.; Guo, Y. Penalty dynamic programming algorithm for dim targets detection in sensor systems. Sensors
**2012**, 12, 5028–5046. [Google Scholar] [CrossRef] [PubMed] - Grossi, E.; Lops, M.; Venturino, L. A Novel dynamic programming algorithm for track-before-detect in radar systems. IEEE Trans. Signal. Process.
**2013**, 61, 2608–2619. [Google Scholar] [CrossRef] - Ristic, B.; Arulampalam, S.; Gordon, N. Beyond the Kalman. Filter-Particle Filters for Tracking Applications; Artech House: London, UK, 2004. [Google Scholar]
- Salmond, D.J.; Birch, H. A Particle Filter for Track-Before-Detect. In Proceedings of the American Control Conference, Arlington, VA, USA, 25–27 June 2001; pp. 3755–3760.
- Boers, Y.; Driessen, H. Particle Filter Based Detection for Tracking. In Proceedings of the American Control Conference, Arlington, VA, USA, 25–27 June 2001; pp. 4393–4397.
- Davey, S.J.; Rutten, M.G. A Comparison of Three Algorithms for Tracking Dim Targets. In Proceedings of the Information, Decision and Control (IDC’07), Queensland, Australia, 12–14 February 2007; pp. 342–347.
- Davey, S.; Rutten, M.; Cheung, B. Using phase to improve track before detect. IEEE Trans. Aerosp. Electron. Syst.
**2012**, 48, 832–849. [Google Scholar] [CrossRef] - Boers, Y.; Driessen, J. Multitarget particle filter track before detect application. IEE Proc. Radar Sonar Navig.
**2004**, 151, 351–357. [Google Scholar] [CrossRef] - Lepouter, A.; Rabaste, O. A Particle Filter for Target Arrival Detection and Tracking in Track-Before-Detect. In Proceedings of the Workshop on Sensor Data Fusion: Trends, Solutions, Applications (SDF’12), Bonn, Germany, 4–6 September 2012; pp. 13–18.
- Lu, J.; Shui, P.-L. Track-before-detect method based on cost-reference particle filter in non-linear dynamic systems with unknown statistics. IET Signal. Process.
**2014**, 8, 85–94. [Google Scholar] [CrossRef] - Hu, Z.B.; Xiong, S.W. Study of hybrid differential evolution based on simulated annealing. Comput. Eng. Des.
**2007**, 28, 1989–1992. [Google Scholar] - Zhang, C.; Li, L. Hybrid differential evolution particle filter for nonlinear filtering. Appl. Comput. Electromagn. Soc.
**2014**, 12, 1133–1139. [Google Scholar] - Douc, R.; Cappe, O. Comparison of Resampling Schemes for Particle Filtering. In Proceedings of the Image and Signal Processing and Analysis, Zagreb, Croatia, 15–17 September 2005; pp. 64–69.

© 2015 by the authors; licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution license (http://creativecommons.org/licenses/by/4.0/).