Cost-Reference Particle Filter for Cognitive Radar Tracking Systems with Unknown Statistics

Abstract

A novel robust particle filtering algorithm is proposed for updating both the waveform and noise parameter for tracking accuracy simultaneously and adaptively. The approach is a significant step for cognitive radar towards more robust tracking in random dynamic systems with unknown statistics. Meanwhile, as an intelligent sensor, it would be most desirable for cognitive radar to develop the application of a traditional filter to be adaptive and to expand the adaptation to a wider scope. In this paper, after analysis of the Bayesian bounds and the corresponding cost function design, we propose the cognitive radar tracking method based on a particle filter by completely reconstructing the propagation and the update process with a cognitive structure. Moreover, we develop the cost-reference particle filter based on optimizing the cost function design according to the complicated system or environment with unknown statistics. With this method, the update of the estimation cost and variance arrives at the approximate optimization, and the estimation error can be more adjacent to corresponding low bounds. Simulations about the tracking implementation in unknown noise are utilized to demonstrate the superiority of the proposed algorithm to the existing methods in traditional radar.

1. Introduction

1.1. Problem Statement

Cognitive radar can basically be defined as: an intelligent radar system of hardware and software in which the transmit and receive parameters (i.e., pulse length, pulse repetition frequency (PRF), modulation, power, frequency, and polarization) are selected, in real-time, and use adaptation between the information extracted from the sensor/processor and the design and transmission of subsequent waveforms, in response to the observed scene to optimize the performance of a given application. The problem of target tracking in cognitive radar system has received considerable attention. It is well known that most of the state-space dynamical systems are nonlinear or non-Gaussian. For example, the multi-state transition of a drone from hovering to maneuvering is often nonlinear, and the measurement noise is mainly flicker noise or heavy-tailed noise [1].

Despite the performance decline of extended Kalman filter (EKF) and unscented Kalman filter (UKF) in highly nonlinear problems, and that general analytical solutions are intractable in nonlinear or non-Gaussian systems, solutions continue to emerge from different viewpoints. However, most of the methods rely on the assumption that the noise has known statistics, or they require accurate mathematical representation of the dynamics of the system evolution; otherwise, it is almost impossible to directly approximate the true distribution. In practice, sometimes the assumption is in accordance with the actual situation, but sometimes not. The same situation occurs in cognitive radar systems. Without explicit mathematical models or a priori information, how cognitive radar can still exert its advantages is a concern of this article.

We briefly review various existing methods to cognitive radar tracking problem that involve particle filter (PF) in some relevant manner.

1.2. Related Works

For handling arbitrary nonlinear models and arbitrary noise distributions, methods based on Monte Carlo (MC) methodology have emerged in Bayesian estimation. MC is a simulation-based method aimed at estimating the posteriori pdf of the state given the observations. Markov Chain Monte Carlo (MCMC) and Sequential MC (SMC) are two main tools in it, and they can sample from high dimensional probability distributions. The performance of MCMC would be unreliable when the proposal distributions that are used to explore the space are poorly chosen and/or if highly correlated variables are updated independently. While, PF became popular for it is particularly suitable for real-time estimation.

The Sequential Importance Sampling (SIS) algorithm was used for the first time to solve the problem of nonlinear filtering [2]. The formal establishment of PF was attributed to the proposal of resampling technology [3,4]. Meanwhile, the idea of Sampling Importance Resampling (SIR) was discovered and developed [5]. PF is allowed to express a complete and precise state posterior distribution, so any statistical data such as mean, variance, and modulus can be easily calculated, and theoretically, the accuracy is higher than that of EKF and UKF in nonlinear systems [6]. Hence, despite its computational pressure, it is still quite attractive to us. During the development of PF, there are several inevitable drawbacks need to be addressed: (1) How to approximate the optimal proposal distribution further; (2) How to overcome the problems of weight degeneracy, sample impoverishment, etc., to make the resampling more effective; (3) How to make the algorithm efficient and online.

Therefore, more variants of PF were produced [7]. Efficient importance sampling techniques were studied in [8,9] for the first problem, but these algorithms require that the posterior distribution of the states are assumed as a priori known and can be approximated by a Gaussian distribution. An auxiliary variable particle filter was used to deal with the second problem [10], but the filtering performance degrades when the state noise is strong, and because the likelihood function and weight value need to be calculated twice for each particle, the calculation amount increases. Adaptive PF (APF) can release the computation burden by adjusting the number of particles dynamically [11,12,13,14,15]. This method chooses a small particle number if the density is focused on a small part of the state-space, and chooses a large number if the state uncertainty is high [16]. KLD-APF is based on Kullback–Leibler (KL) information or KL distance (KLD) sampling, but it ignores any mismatch between the true and the proposal distribution [17]. When the mismatch happens, the addition of particle number will only increase the computational load. On the other hand, reducing the particle number may aggravate the sample impoverishment and further weaken the effect of resampling.

In addition to these drawbacks, PF would still be invalid when some issues are encountered. PF requires that the conditional pdf of the observed variable can be estimated, otherwise the weight of the particles cannot be calculated. Thus, particle MCMC (PMCMC) was proposed by combing SMC with MCMC as an efficient approach [18]. Based on PMCMC, the SMC² algorithm is motivated to tackle the intractable problem of probable increments in state-space models [19]. Similar to the SMC² scheme, nested particle filters for online parameter estimation were proposed but in a purely recursive manner, in order to address the problem of approximating the posterior probability distribution of the fixed parameters of a state-space dynamical system [20]. On robustness, cooperative parallel particle filters were designed for the dual purpose of Bayesian inference and on-Line Model Selection, with the online adaptation for the particle number [21]. Particle learning was proved in [22] to be outperforming existing PF alternatives and a competitor to MCMC. In the presence of model uncertainty where discrete data are encountered, a new SMC method was proposed for the filtering and prediction of time-varying signals [23]. A similar method with better predictive powers was proposed in [24], wherein the resampling step could be dynamically adjusted and the predictive powers could be updated sequentially as more data were observed. Moreover, what deserves attention is that Bayes is not the only choice, as a neural filter based on GRNN is also outperforming numerical filters in state vector estimation during dynamic changes of target movement parameters [25].

Inspired by the previous research, we might find a solution independent of the prior information through updating the proposal distribution and its stepwise approximation to the true distribution by iterative methods. The authors of [26] designed a cost function in PF to iteratively update the state and variance when the prior information is unknown. Naturally, it has been applied to many areas such as target tracking, autonomous vehicle positioning, sensor network, orthogonal frequency division multiplexing (OFDM) systems, wireless local area network (WLAN), and tilt estimation [27,28,29]. There are some improvements to the algorithm. In [30], a particle selection algorithm was proposed and analyzed for implementation with parallel computing devices and to circumvent the main drawback of the conventional resampling techniques. Authors of [31] melded random measures of two or more cost-reference particle filters to obtain a fused random measure that combines the information from the individual cost-reference particle filters. However, there is no further study on the optimization of the cost function, especially when the tracking structure can be adaptive. Cognitive radar designs the optimal estimator for Bayesian framework in [32], but the majority of the existing Bayesian tracking methods for cognitive radar applications are based on the Kalman filter for linear systems [33,34], and Kalman-like solutions, e.g., cubature KF (CKF) and continuous-discrete (CD)-CKF for nonlinear problems [35]. Few researches use PF. In [36], a particle filter combined with probabilistic data association is used as a tracker. In [37], a cognitive structure was designed as only a part of PF, namely, a parallel structure of PF, while EKF was used by waveform selection to adapt the particle number and reduce the computation cost.

1.3. Contributions and Organization of the Paper

When we specify PF as the nonlinear tracking method in a cognitive radar, we consider using the cognitive structure to expand the state-space of PF or its variants to another dimension where the waveform parameters can be changed from fixed to dynamic, and we can optimize the cost function design and corresponding lower bound of the estimation error. In turn, using the particle filter tracking method based on the optimal cost function can expand the scope of the application of cognitive radar. The main contributions of this paper will be presented by the following items:

(1) The cognitive radar tracking method based on PF is proposed and the mathematical model of which it is derived for the first time. Push the cognitive radar tracking framework from existing the Kalman-like to the developing SMC-like. Not only the data process but also the waveform design is in the PF iteration, that is, a fully cognitive PF.

(2) Refine the idea of cost-reference PF with cognitive framework. When the estimation of the parameters (i.e., signal-to-noise ratio (SNR)) in cognitive radar mismatches the real situation, a novel cognitive cost-reference particle filter algorithm is proposed to bring about robustness to the existing cognitive radar tracking methods and intelligence for the current adaptive method. The convergence is proofed mathematically.

(3) The Cramér–Rao Lower Bound (CRLB) of the proposed cognitive PF is derived, and the corresponding cost function is designed.

The rest of this paper is organized as follows: the standard PF, the mathematical model of cognitive radar, and the interface of PF to cognitive algorithm are presented, and the CRLB of the cognitive PF is derived in Section 2. The principle of the cost-reference particle filter is restated, the proposed cognitive scheme is presented along with the implementation steps, and the convergence of the solution is proofed in Section 3. Section 4 shows the dynamic model, the maneuvering target tracking in unknown non-Gaussian noise, and the simulation results. Section 5 presents the conclusions.

2. PF for Cognitive Radar Tracking

In this section, the cognitive tracking method based on PF is formed to address the non-linear and non-Gaussian state estimation problem in cognitive radar.

2.1. PF and Cognitive Radar Model

2.1.1. Standard PF

The basic idea is to generate a set of random samples in the state-space according to the empirical conditioning distribution of the system state vector. They are called particles, the weight and position of which are continuously adjusted according to the measurement, and the initial empirical conditioning distribution is modified. The MC estimate of the integral can be written as:

$$E\left[f\left(x\right)\right]={\displaystyle \int f\left(x\right)p\left(x\right)dx}={\displaystyle \int f\left(x\right)\left[p\left(x\right)/q\left(x\right)\right]q\left(x\right)dx}.$$

(1)

Let $\omega \left(x_{}^{\left(i\right)}\right)=p\left(x_{}^{\left(i\right)}\right)/q\left(x_{}^{\left(i\right)}\right)$ denote the importance weights. Equation (1) is calculated by generating $N\gg 1$ independent samples $\left\{x_{}^{\left(i\right)}:i=1,\dots ,N\right\}$ distributed according to $q\left(x\right)$. The weighted sum is formed and normalized as:

$$E_N\left[f\left(x\right)\right]=\frac{1}{N}{\displaystyle \sum _{i=1}^{N}f\left(x_{}^{\left(i\right)}\right)\omega \left(x_{}^{\left(i\right)}\right)},$$

(2)

$$E\left[f\left(x\right)\right]=\frac{\frac{1}{N}{\displaystyle \sum _{i=1}^{N}f\left(x_{}^{\left(i\right)}\right)\omega \left(x_{}^{\left(i\right)}\right)}}{\frac{1}{N}{\displaystyle \sum _{j=1}^N\omega \left(x_{}^{\left(j\right)}\right)}}={\displaystyle \sum _{i=1}^{N}f\left(x_{}^{\left(i\right)}\right)\tilde{\omega }\left(x_{}^{\left(i\right)}\right)}.$$

(3)

The weight update equation can be shown to be:

$$\omega \left(x_k^{\left(i\right)}\right)=\omega \left(x_{k-1}^{\left(i\right)}\right)\frac{p\left(z_k^{}|x_k^{\left(i\right)}\right)p\left(x_k^{\left(i\right)}|x_{k-1}^{\left(i\right)}\right)}{q\left(x_k^{\left(i\right)}|x_{k-1}^{\left(i\right)}\right)}.$$

(4)

Substitution of $q\left(x_k^{\left(i\right)}|x_{k-1}^{\left(i\right)}\right)=p\left(x_k^{\left(i\right)}|x_{k-1}^{\left(i\right)}\right)$ into Equation (4) yields:

$$\omega \left(x_k^{\left(i\right)}\right)=\omega \left(x_{k-1}^{\left(i\right)}\right)p\left(z_k^{}|x_k^{\left(i\right)}\right).$$

(5)

When $N\to \infty$, the posterior probability density can be approximated by the following:

$$p\left(x_k^{}|z_{1:k}^{}\right)\approx {\displaystyle \sum _{i=1}^N\omega \left(x_k^{\left(i\right)}\right)\delta \left(x_k^{}-x_k^{\left(i\right)}\right)},$$

(6)

where $\delta (\cdot )$ denotes the Dirac delta measure, with defining properties: (1) $\delta \left(x_k^{}-x_k^{\left(i\right)}\right)=0$ for $x_k^{}\ne x_k^{\left(i\right)}$, (2) $\delta \left(x_k^{}-x_k^{\left(i\right)}\right)=+\infty$ for $x_k^{}=x_k^{\left(i\right)}$; and (3) ${\displaystyle {\int }_{-\infty }^{\infty }\delta \left(x_k^{}-x_k^{\left(i\right)}\right)d}{x}_k^{}=1$. More details about the SIR-PF can refer to [38]. It has been termed standard PF (SPF) in the simulation below. The accuracy approximates the optimal estimation.

2.1.2. Cognitive Radar System Model

The radar transmitter through the emission exploiting waveform selection stimulates the background with the goal to obtain a response from it such as a radar echo. The mentioned response is perceived by radar receiver, which plays the equivalent role of the human senses. The detector performs a low-level processing of received sensor data. The scene analyzer is conceived by exploiting previously acquired information to estimate the statistical characteristics of the operating environment. The tracker acts as a high-level processor.

The functions manager receives requests for radar activities (e.g., search, tracking) from the keyboard according to a specific mode of working or operator requirements. The functional interaction is with the radar hardware, the signal processor, the data processor, and the external peripherals. Starting from the state estimate ${\widehat{\mathbf{x}}}_{k-1}^{}$ and state error covariance at time step k − 1, the transmitter uses waveform parameter ${\mathbf{\theta }}_{k-1}$ to illuminate the environment. The receiver performs a prediction and obtains ${\widehat{x}}_k^{\left(i\right)}$ (i.e., Equation (49)) with the previous value and current measurement. Then, ${\widehat{x}}_k^{\left(i\right)}$ is utilized to evaluate the CRLB, and thus, the measurement error covariance matrix $R_k\left(\mathbf{\theta }\right)$. For each waveform of the library, the predicted error covariance matrix $P_k$ is computed for evaluation (i.e., Equation (5) in Algorithm 1). At this point, a cost function (i.e., Equation (18)) is to be formulated, and the waveform at the next transmission is chosen as the one that minimizes the weighted mean square error (MSE) on the target state, namely minimizing the cost function by the optimal waveform design or parameters adjustment. The transmitted signal then provides a new optimized measurement ${\mathbf{z}}_k$ to feed the PF processor and track the target. The obtained target state estimate ${\widehat{x}}_k^{}$ is finally exploited to predict the next target state ${\widehat{x}}_{k+1}^{\left(i\right)}$, and optimizes the new transmitted waveform according to the closed-loop paradigm. This constructs the PAC in cognitive radar [39].

Summarizing, the block scheme of a cognitive tracking system architecture and the information-flow in the mathematical model is displayed in Figure 1. Notice that, if the blocks in red are cancelled, the rest blocks in blue degrade to the traditional radar system.

2.2. PF-Based Cognitive Tracking Algorithm

2.2.1. Bayesian Bounds for Cognitive Radar

The recursive Bayesian estimators, such as the EKF or PF, work in the non-linear system employing the feed-forward processing chain, and do not consider the feedback information from the receiver to the transmission of subsequent waveforms. So far, the corresponding recursive BCRB on the MSE matrix of the state vector was derived and presented in [40]. It can be obtained by calculating the Fisher information matrix (FIM) [41], the inverse of which is the unbiased estimator of the CRLB of the state estimation error covariance [42]. In this paper, the PF based on the cognitive structure characterized by feedback is analyzed, and the estimation is biased. Bayesian FIM (BIM) for this case can be developed and computed iteratively using this form [43]:

$$B_k^{\uparrow }\left({\mathbf{\theta }}_k|{\mathbf{z}}_{1:k-1};{\mathsf{\Theta }}_{k-1}\right)=B_k^{-}\left({\mathbf{\theta }}_k|{\mathbf{z}}_{1:k-1};{\mathsf{\Theta }}_{k-1}\right)+J_k^{-}\left({\mathbf{\theta }}_k|{\mathbf{z}}_{1:k-1};{\mathsf{\Theta }}_{k-1}\right),$$

(7)

where $B_k^{-}\left({\mathbf{\theta }}_k|{\mathbf{z}}_{1:k-1};{\mathsf{\Theta }}_{k-1}\right)$ represents the contribution of the prior information is prior, called the prior term, and is denoted as $J_P$; the $J_k^{-}\left({\mathbf{\theta }}_k|{\mathbf{z}}_{1:k-1}^{};{\mathsf{\Theta }}_{k-1}\right)$ represents the contribution of the data, called data term and denoted as $J_D$.

$$J_D=E_k\left\{J_F\left(x_k^{}\right)\right\}=E_k\left\{-E_{{\mathbf{z}}_k^{}|x_k^{};{\mathbf{\theta }}_k}\left\{\frac{{\partial }^2\ln f\left({\mathbf{z}}_k^{}|x_k^{};{\mathbf{\theta }}_k\right)}{\partial x_k^2}\right\}\right\},$$

(8)

$$J_P=-E_{x_k^{}}\left\{\frac{{\partial }^2\ln f\left(x_k^{}\right)}{\partial x_k^2}\right\}.$$

(9)

The lower bound is defined as the inverse of the BIM, as the form:

$$BCRB={\left[J_D+J_P\right]}^{-1}.$$

(10)

The lower bound of CR would be smaller than the general BCRB, which can be proved as:

$${\left[J_D\left({\mathbf{\theta }}_k\right)+J_P\left({\mathbf{\theta }}_k\right)\right]}^{-1}\le {\left[J_D+J_P\right]}^{-1}=BCRB.$$

(11)

2.2.2. Cost Function Design

As we know, CRLB is the possible MMSE of an unbiased estimation. It can be proved that minimizing the trace of the covariance is equal to the MMSE of the expectation of tracking state [38]; thus, we choose the minimum trace with respect to the waveform parameters as the processor cost function modeled by:

$${\mathrm{C}}_{\mathsf{\Theta }}\left({\mathbf{\theta }}_k|{\mathbf{z}}_{1:k-1}^{};{\mathsf{\Theta }}_{k-1}\right)=E_{x_k z_k|I_k {\mathbf{\theta }}_k}\left[{\left(x_k-{\widehat{x}}_k\left({\mathbf{z}}_k^{}\right)\right)}^T\left(x_k-{\widehat{x}}_k\left({\mathbf{z}}_k^{}\right)\right)\right]\approx \mathrm{tr}\left\{B_k^{\uparrow }{\left({\mathbf{\theta }}_k|{\mathbf{z}}_{1:k-1}^{};{\mathsf{\Theta }}_{k-1}\right)}^{-1}\right\},$$

(12)

where $\mathrm{tr}(\cdot )$ is an operator that extracts the trace of $B_k^{\uparrow }{\left({\mathbf{\theta }}_k|{\mathbf{z}}_{1:k-1}^{};{\mathsf{\Theta }}_{k-1}\right)}^{-1}$, and the waveform parameters are denoted as $\mathsf{\Theta }\equiv \left\{{\mathbf{\theta }}_1,{\mathbf{\theta }}_2,\dots ,{\mathbf{\theta }}_k\right\}$. ${\widehat{x}}_{k+1}\left(I_k,x_{k+1},z_{k+1};{\mathbf{\theta }}_k\right)$ is the posterior predictive value of the state estimation given waveform parameter. Since ${\widehat{x}}_{k+1}$ is related to ${\mathbf{z}}_{k+1}$, which is a nonlinear mapping of $x_{k+1}$, we use approximation rather than integration to solve this expectation. ${\mathrm{C}}_{\mathsf{\Theta }}\left({\mathbf{\theta }}_k|{\mathbf{z}}_{1:k-1}^{};{\mathsf{\Theta }}_{k-1}\right)$ is approximated by $\mathrm{tr}\left\{B_k^{\uparrow }{\left({\mathbf{\theta }}_k|{\mathbf{z}}_{1:k-1}^{};{\mathsf{\Theta }}_{k-1}\right)}^{-1}\right\}$ [33,35]. The model of $P_{k+1|k+1}$ in KF can refer to [44], while in PF, $P_{k+1|k+1}$ is modeled by a different form [38]. More details of the relation between $P_{k+1|k+1}$ and $R_k\left({\mathbf{\theta }}_{k-1}\right)$ are referred to in [35]. We derive the model of measurement noise covariance based on waveform selection as [45]:

$$R_k\left({\mathbf{\theta }}_{k-1}\right)=\mathsf{\Gamma }{B}_k^{\uparrow }{\left({\mathbf{\theta }}_k|{\mathbf{z}}_{1:k-1}^{};{\mathsf{\Theta }}_{k-1}\right)}^{-1}\mathsf{\Gamma }=\mathsf{\Gamma }{U}^{-1}\left({\mathbf{\theta }}_{k-1}\right)\mathsf{\Gamma }/\eta ,$$

(13)

$$\mathsf{\Gamma }\triangleq \mathrm{diag}\left[\frac{c}{2},\frac{c}{4\pi f_c}\right],$$

(14)

where $\eta =2E_R/N_0$ is the SNR, and $N_0$ denotes the spectral density of the complex noise envelope $\tilde{n}\left(k\right)$. $U\left({\mathbf{\theta }}_{k-1}\right)$ is a scaled version of FIM. $\mathsf{\Gamma }$ denotes a symmetric matrix defined as this form. $c$ is the speed of waveform propagation, and $f_c$ is the carrier frequency.

When the LFM signal with Gaussian amplitude modulation is selected as the transmit waveform, with the duration of Gaussian pulse λ and the chirp rate b chosen to form the parameter vector ${\mathbf{\theta }}_{k-1}=\left[\lambda ,b\right]$, the measurement noise covariance matrix can be modeled by $R_{k,r\dot{r}}\left({\mathbf{\theta }}_{k-1}\right)$ [41]:

$$R_k\left({\mathbf{\theta }}_{k-1}\right)=R_{k,r\dot{r}}\left({\mathbf{\theta }}_{k-1}\right)=\left[\begin{array}{cc}\frac{c^2{\lambda }^2}{2\eta }& -\frac{c^{2}b{\lambda }^2}{2\pi f_c\eta }\\ -\frac{c^{2}b{\lambda }^2}{2\pi f_c\eta }& \frac{c^2}{{\left(2\pi f_c\right)}^2\eta }\left(\frac{1}{2{\lambda }^2}+2b^2{\lambda }^2\right)\end{array}\right].$$

(15)

The sensor cost function is defined as:

$${\mathrm{C}}_{\mathsf{\Theta }}\left({\mathbf{\theta }}_k\right)=\{\begin{array}{ll}0& {\mathbf{\theta }}_k\in \mathsf{\Theta }\\ \infty & \mathrm{otherwise}\end{array}.$$

(16)

The loss function is defined with respective to the cost function of the processor and sensor, modeled by:

$$L_{\mathcal{C},\mathsf{\Theta }}\left({\mathbf{\theta }}_k|{\mathbf{z}}_{1:k-1}^{};{\mathsf{\Theta }}_{k-1}\right)\equiv {\mathrm{C}}_{\mathsf{\Theta }}\left({\mathbf{\theta }}_k|{\mathbf{z}}_{1:k-1}^{};{\mathsf{\Theta }}_{k-1}\right)+{\mathrm{C}}_{\mathsf{\Theta }}\left({\mathbf{\theta }}_k\right).$$

(17)

The choice of the optimal parameter is equivalent to the optimization problem given by:

$${\mathbf{\theta }}_k=\underset{{\mathbf{\theta }}_k^{}\in \mathrm{P}}{\arg \min }{\mathrm{C}}_{\mathsf{\Theta }}\left(\mathbf{\theta }|{\mathbf{z}}_{1:k-1}^{};{\mathsf{\Theta }}_{k-1}\right) \mathrm{s}.\mathrm{t}. \mathbf{\theta }\in \mathsf{\Theta }.$$

(18)

2.2.3. Algorithm and Numerical Simulation

We have restated the principle of the cognitive radar system, and have given the interface from PF to it to induce the human cognition mechanism into the PF algorithm. As with the statements before, SPF is available to be the tracking method of the cognitive radar for a tracking target with a nonlinear dynamic and non-Gaussian distribution. The basic algorithm steps of SPF along with the interface to cognitive tracking framework is given as follows Algorithm 1:

Algorithm 1. Cognitive Radar Tracking Recursion Based on PF

Initialization

1.    ${\widehat{x}}_0^{\left(i\right)}\sim p\left(x_0\right)$, ${\mathsf{\omega }}_k^{(i)}=1/N_s$, $i=1,\dots ,N_s$

Controller Optimization

2.    ${\widehat{x}}_k^{\left(i\right)}\sim q\left(x_k^{(i)}|{\widehat{x}}_{k-1}^{(i)},{\mathbf{z}}_k\right)$, the mean $E_{p_k\left(x_k^{}|{\widehat{x}}_{k-1}^{\left(i\right)}\right)}\left[{\widehat{x}}_k^{}\right]=f_x\left({\widehat{x}}_{k-1}^{\left(i\right)}\right)$, $i=1,\dots ,N_s$

3.    ${\mathsf{\omega }}_k^{(i)}\left(\mathbf{\theta }\right)\propto p\left({\mathbf{z}}_k|{\widehat{x}}_k^{(i)}\right)p\left({\widehat{x}}_k^{(i)}|x_{k-1}^{(i)}\right)/q\left(x_k^{(i)}|x_{0:k}^{(i)},{\mathbf{z}}_{1:k}\right)$$, i=1,\dots ,N_s$$, \mathrm{normalize} {\tilde{\mathsf{\omega }}}_k^{(i)}\left(\mathbf{\theta }\right)={\mathsf{\omega }}_k\left(x_{0:k}^{\left(i\right)}\right)/{\displaystyle \sum _{i=1}^{N_s}{\mathsf{\omega }}_k\left(x_{0:k}^{\left(i\right)}\right)}$

4.    $\mathrm{If} N_{eff}=1/{\displaystyle \sum _{i=1}^{N_s}{\left({\tilde{\mathsf{\omega }}}_k^{\left(i\right)}\right)}^2}<N_{th}\left(\mathrm{empirical}\right)$$, \mathrm{then} \left[{\left\{{\widehat{x}}_k^{\left(i\right)},{\widehat{\mathsf{\omega }}}_k^{(i)}\left(\mathbf{\theta }\right)\right\}}_{i=1}^{N_s}\right]=\mathrm{RESAMPLE}\left[{\left\{x_k^{\left(i\right)},{\tilde{\mathsf{\omega }}}_k^{(i)}\left(\mathbf{\theta }\right)\right\}}_{i=1}^{N_s}\right]$

5.    ${\mathbf{P}}_k\left(\mathbf{\theta }\right)={\displaystyle \sum _{i=1}^{N_s}{\mathsf{\omega }}_k^{(i)}\left(\mathbf{\theta }\right)\left({\widehat{x}}_k^{\left(i\right)}-{\widehat{x}}_k^{}\right){\left({\widehat{x}}_k^{\left(i\right)}-{\widehat{x}}_k^{}\right)}^T}$

6.    ${\mathbf{\theta }}_k^{*}=\arg \underset{{\mathbf{\theta }}_k^{}\in {P}}{\min }\left[\mathrm{Tr}\left({\overline{\mathbf{P}}}_{k+1|k+1}^{}\left({\mathbf{\theta }}_k\right)\right)\right]$$. \mathrm{Select} \mathrm{the} \mathrm{optimal} {\mathbf{\theta }}_k={\mathbf{\theta }}_k^{*}$

Motion Update and Measurement

7.    $f^{-}\left({\widehat{x}}_{k-1}^{(i)}\right)={\displaystyle \int q\left(x_k^{(i)}|x_{k-1}^{(i)};{\mathbf{\theta }}_k\right)f\left({\widehat{x}}_{k-1}^{(i)}\right)d{x}_{k-1}^{}}$$, {\widehat{x}}_{k-1}^{\left(i\right)}={\mu }_k^{-}=E_k^{-}\left[{\widehat{x}}_k^{}\right]$

8.    ${\widehat{y}}_k=\arg \underset{y}{\min }\ln f\left({\mathbf{z}}_k^{}|y;{\mathbf{\theta }}_k\right)$

Information Update and State Estimation

9.    $\left\{{\widehat{x}}_k^{\left(i\right)},{\mathsf{\omega }}_k^{(i)}\right\}=\left\{x_k^{\left(i\right)},{\tilde{\mathsf{\omega }}}_k^{(i)};{\mathbf{\theta }}_k\right\}$$, i=1,\dots ,N_s$$, {\widehat{x}}_k\approx {\displaystyle \sum _{i=1}^{N_s}{\mathsf{\omega }}_k^{\left(i\right)}{\widehat{x}}_k^{\left(i\right)}}$

As an illustration of the CRLB, consider the scenario that the target is assumed to move with a constant velocity along a straight line in the Cartesian plane. The target kinematic state (position, velocity) ${\left[x_{1,k}^{},x_{2,k}^{},x_{3,k}^{},x_{4,k}^{}\right]}^T$ is estimated from the noise-corrupted measurements ${\mathbf{z}}_k={\left[r_k, {\beta }_k\right]}^T+w_k$ with $r_k={\left(x_{1,k}^2+x_{2,k}^2\right)}^{1/2}$, ${\beta }_k=\arctan \left(x_{2,k}^{}/x_{1,k}^{}\right)$. Thus, the measurement equation is nonlinear. The measurement noise $w_k$ is assumed to be white, zero-mean Gaussian, with covariance $R=diag\left[{\sigma }_r^2,{\sigma }_{\beta }^2\right]$, and with the sampling interval T = 0.1 s. The covariance of the process noise is $Q={\sigma }_x^2\cdot diag\left[T_0,T_0\right]$, wherein $T_0=[0.5T^4 0.5T^3; 0.5T^3 T^2]$. The initial true state of the target is ${\left[4000, 80, 12,000, -20\right]}^T$, and the given initial state is set as ${\left[5000, 90, 15,000, -20\right]}^T$, with an estimation error covariance of $P_0=1000\mathrm{I}$.

Figure 2 shows the results in terms of posterior Cramér–Rao bound (PCRB) and root mean square errors (RMSEs) in the range estimation. As observed in Figure 2a, when we perform the proposed CPF for target tracking in cognitive radar systems, the MSE is initially larger than PCRB, due to the initialization that is not exactly matched to PCRB. Very soon though, CPF demonstrates a fast convergence. Then, MSE agrees with the theoretical curve for PCRB, presenting a considerably well overall performance. Figure 2b compare the PCRB curves obtained using a fixed waveform in a traditional radar and dynamic waveform in cognitive radar. Obviously, the PCRB in the cognitive radar provides a valid lower bound for all time steps, and the convergence is faster.

As we know, the CRLB does not care about the estimation method, but only reflects the best effect of using the existing available information to estimate the parameters. Thus, it can be seen from the figures that after the cognitive system utilizes the extra environmental information through the closed-loop structure, it is possible to lower the previous bounds further.

3. Cognitive Tracking Problem Model in Unknown Environment

In this section, we develop the cost-reference PF (CRPF) for a cognitive robust filter. This means that we extend the classic cognitive radar scenarios to an unknown environment. Specifically, the general description of the CRPF approach is listed first, then the details on fusing the cost with cognitive method are specified.

3.1. Cost-Reference Particle Filtering Approach

As we all know, SPF works well when the mathematical forms of the probability distributions of the noise are assumed explicit, but when the mathematical representation of the dynamics of the system evolution is unknown a priori, or the assumptions of probabilistic models cannot be achieved, SPF may fail and the estimation results may be inaccurate even if they have the cognitive structure.

Unlike basic PF, which only gives the state vectors in sequential algorithm, CRPF gives the state samples and associated costs, namely the WPS as the form:

$${\mathsf{\Xi }}_k={\left\{x_k^{\left(i\right)},{\mathcal{C}}_k^{\left(i\right)}\right\}}_{i=1}^M,$$

(19)

where ${\mathcal{C}}_k^{\left(i\right)}=\mathcal{C}\left(x_{0:k}^{\left(i\right)}|y_{1:k}^{},\lambda \right)$ is the cost of the particle $x_k^{\left(i\right)}$. The cost function can be denoted by a recursive additive structure:

$$\mathcal{C}\left(x_{0:k}^{}|{\mathbf{z}}_{1:k}^{},\lambda \right)=\lambda \mathcal{C}\left(x_{0:k-1}^{}|{\mathbf{z}}_{1:k-1}^{},\lambda \right)+\Delta \mathcal{C}\left(x_k^{}|{\mathbf{z}}_k^{}\right),$$

(20)

where $\lambda$ is a forgetting factor, which is used to avoid attributing an excessive weight to old observations. $\Delta \mathcal{C}\left(x_k^{}|{\mathbf{z}}_k^{}\right)$ is the incremental cost function, the prediction of which can be obtained as:

$$\Delta \mathcal{C}\left(f_x\left(x_{k-1}^{}\right)|{\mathbf{z}}_k^{}\right)=\mathcal{R}\left(x_{k-1}^{}|{\mathbf{z}}_k^{}\right),$$

(21)

where the one-step risk function is introduced, and the risk of particle $i$ can be computed as:

$$\mathcal{R}:{\mathbb{R}}^{L_x}\times {\mathbb{R}}^{L_y}\to \mathbb{R}, x_{k-1}^{},{\mathbf{z}}_k^{}\to \mathcal{R}\left(x_{k-1}^{}|{\mathbf{z}}_k^{}\right),$$

(22)

$${\mathcal{R}}_{k+1}^{\left(i\right)}=\lambda {\mathcal{C}}_k^{\left(i\right)}+\mathcal{R}\left(x_k^{\left(i\right)}|{\mathbf{z}}_{k+1}^{}\right).$$

(23)

We assign the particles in ${\mathsf{\Xi }}_{k+1}$ the PMF, which is, according to the cost, defined as:

$${\pi }_{k+1}^{\left(i\right)}\propto \mu \left({\mathcal{C}}_{k+1}^{\left(i\right)}\right),$$

(24)

where a monotonically decreasing function can be used as $\mu (\cdot )$, which is selected to guarantee an adequate discrimination of low-cost particles from higher ones. Then, we would obtain the minimum cost by maximum the function below:

$${\tilde{x}}_{0:k+1}^{\min }=x_{k+1}^{\left(i_0\right)},i_0=\arg \underset{i}{\max }\left\{{\pi }_{k+1}^{\left(i\right)}\right\}.$$

(25)

Notice that the mean value here is equal to the minimum cost estimate, and it even has slight advantages over the latter one. More details can be referred to in [26], and the presentation can be found in Algorithm 2.

The overall procedure of the CRPF algorithm is summarized in Algorithm 2.

Algorithm 2. CRPF algorithm for target tracking problem

Initialization

1.    $x_0^{\left(i\right)}\sim p_0\left(x_0\right)$$, {\mathcal{C}}_0^{\left(i\right)}=0$$, {\mathbf{\sigma }}_0^{2,\left(i\right)}$$, i=1,\dots ,N_s$$, \mathrm{the} \mathrm{weighted}\text{-}\mathrm{particle} \mathrm{set} {\mathsf{\Xi }}_0={\left\{x_0^{\left(i\right)},{\mathcal{C}}_0^{\left(i\right)}\right\}}_{i=1}^{N_s}$

PMF Update

2.    ${\mathcal{R}}_k^{\left(i\right)}=\lambda {\mathcal{C}}_{k-1}^{\left(i\right)}+{\Vert {\mathbf{z}}_k-f_y\left(f_x\left(x_{k-1}^{\left(i\right)}\right)\right)\Vert }^q$$, q=1,2$$, {\widehat{\mathbf{\pi }}}_k^{\left(i\right)}\propto \mu \left({\mathcal{R}}_k^{\left(i\right)}\right)=\frac{1}{{\left({\mathcal{R}}_k^{\left(i\right)}-\min {\left\{{\mathcal{R}}_k^{\left(i\right)}\right\}}_{i=1}^{N_s}+\delta \right)}^{\beta }}$$, i=1,\dots ,N_s$

3.    ${\left\{{\widehat{x}}_{k-1}^{\left(i\right)},{\widehat{\mathcal{C}}}_{k-1}^{\left(i\right)}\right\}}_{i=1}^{N_s}=\mathrm{RESAMPLE}\left[{\left\{{\widehat{\mathbf{\pi }}}_k^{\left(i\right)}\right\}}_{i=1}^{N_s}\right]$

Particle Propagation and Variance Update

4.    $E_{p_k\left(x_k^{}|{\widehat{x}}_{k-1}^{\left(i\right)}\right)}\left[x_k^{}\right]=f_x\left({\widehat{x}}_{k-1}^{\left(i\right)}\right)$$, Co{v}_{p_k\left(x_k^{}|{\widehat{x}}_{k-1}^{\left(i\right)}\right)}\left[x_k^{}\right]={\mathbf{\sigma }}_k^{2,\left(i\right)}{\mathit{I}}_{\left[x\right]}$$, x_k^{\left(i\right)}\sim p_k\left(x_k^{}|{\widehat{x}}_{k-1}^{\left(i\right)}\right).$

5.    ${\mathbf{\sigma }}_k^{2,\left(i\right)}=\{\begin{array}{ll}{\sigma }_{k-1}^{2,\left(i\right)}& \mathrm{t}\le 10\\ \frac{k-1}{k}{\sigma }_{k-1}^{2,\left(i\right)}+\frac{{\Vert x_k^{\left(i\right)}-g\left({\widehat{x}}_{k-1}^{\left(i\right)}\right)\Vert }^2}{k\times \dim \left[x\right]}& \mathrm{t}>10\end{array}$$, i=1,\dots ,N_s$

State Estimation

6.    ${\mathcal{C}}_k^{\left(i\right)}=\lambda {\mathcal{C}}_{k-1}^{\left(i\right)}+{\Vert {\mathbf{z}}_k-f_y\left(x_k^{\left(i\right)}\right)\Vert }^q$$, {\widehat{\mathbf{\pi }}}_k^{\left(i\right)}\propto {\mu }_2\left({\mathcal{C}}_k^{\left(i\right)}\right)=\frac{1}{{\left({\mathcal{C}}_k^{\left(i\right)}-\min {\left\{{\mathcal{C}}_k^{\left(i\right)}\right\}}_{i=1}^{N_s}+\delta \right)}^{\beta }}$$, \delta ,\beta >0$$, \mathrm{normalized} \mathrm{to} {\mathbf{\pi }}_k^{\left(i\right)},$

7.    ${\widehat{x}}_k=x_k^{mean}={\displaystyle \sum _{i=1}^{N_s}{\mathbf{\pi }}_k^{\left(i\right)}{x}_k^{\left(i\right)}}$$, i=1,\dots ,N_s$

3.2. Cognitive Cost-Reference Particle Filtering

Although the CRPF approach can be used in the dynamic system to deal with the uncertainty problem, by controlling the risk and cost and adapting the variance of the noise, it is also limited by the system performance, e.g., the waveform parameter, process noise, SNR, etc. If the iteration can be performed in another dimension, namely in the varying information entropy caused by the varying parameters, through perceiving the exterior environment, the estimation error may be decreased further. To address this challenge, we propose a cognitive cost-reference particle filter algorithm to perform the robust adaptive filter for target tracking.

3.2.1. Cost Function Design

Considering the more complicated feedback in cognitive CRPF (CCRPF), it seems not easy to achieve the closed-form solution of the CRLB even iteratively [46]. Thus, we use the MC approach to approximate the BIM.

Breaking through the limitation of the traditional CRPF method in Equation (25) that only minimizes the cost of particles in a one-time step, we are dedicated to designing a new cost function based on parameter adaption in cognitive radar framework. The processor cost function is defined with a recursive additive structure as follows:

$$C\left(x_{0:k}^{}|{\mathbf{z}}_{1:k}^{},\lambda ;{\mathbf{\theta }}_k\right)=\lambda {\widehat{C}}_k^{}+\Delta C_{k+1}^{}\left({\mathbf{\theta }}_k\right)={\displaystyle \sum _{i=0}^k{\lambda }^{\left(k-i\right)}\Delta C\left(x_k^{}|{\mathbf{z}}_k^{};{\mathbf{\theta }}_k\right)}.$$

(26)

Lemma 1.

Let ${\mu }_{k,M}$ be the unbiased estimation of mean value ${\mu }_k$, let ${\mathbf{\sigma }}_k^{\left(i\right)}\left({\mathbf{\theta }}_k\right)$ be the variance, and let $\Delta C\left(x_k^{}|y_k^{};{\mathbf{\theta }}_k\right)$ be a cost function.

If the two following conditions are met:

(1) The ${\widehat{\mathbf{\sigma }}}_k^{\left(i\right)}\left({\mathbf{\theta }}_k\right)$ is asymptotically to the ${\mathbf{\sigma }}_k^{\left(i\right)}\left({\mathbf{\theta }}_k\right)$, that is

$$\underset{k\to \infty }{\lim }\left|{\mathbf{\sigma }}_k^{\left(i\right)}\left({\mathbf{\theta }}_k\right)-{\widehat{\mathbf{\sigma }}}_k^{\left(i\right)}\left({\mathbf{\theta }}_k\right)\right|=0\hspace{1em}(i.p.)$$

(27)

where i.p. stands for “in probability”.

(2) The mean incremental cost $\overline{\Delta C_k}={\displaystyle \sum _{i=1}^M{\omega }_k^{\left(i\right)}\Delta C\left(x_k^{\left(i\right)}|y_k^{};{\mathbf{\theta }}_k\right)}$ converges to the minimal incremental cost

$$\underset{M\to \infty }{\lim }\left|\Delta C\left(x_k^{opt}|y_k^{};{\mathbf{\theta }}_k\right)-\overline{\Delta C_k}\right|=0\hspace{1em}(i.p.),$$

(28)

and then ${\mu }_{k,M}$ is asymptotically to ${\mu }_k$, so namely the error satisfies

$$\underset{k\to \infty }{\lim }\left|{\mu }_{k,M}-{\mu }_k\right|=0\hspace{1em}(i.p.).$$

(29)

See Appendix A for a proof.

It can be seen from Equation (26) that the cost function changes with the waveform parameters. According to Lemma 1, we can achieve the asymptotically optimal sequence of state vectors by using the MC approach. The loss function is defined with respect to the cost function of the processor and sensor, modeled by:

$$L_{\mathcal{C},\mathbf{\Theta }}\left({\mathbf{\theta }}_k|{\mathbf{z}}_{1:k-1}^{};{\mathsf{\Theta }}_{k-1}\right)\equiv L\left\{C_{k+1}^{}\left({\mathbf{\theta }}_k|{\mathbf{z}}_{1:k-1}^{};{\mathsf{\Theta }}_{k-1}\right),T_{\mathsf{\Theta }}\left({\mathbf{\theta }}_k\right)\right\},$$

(30)

where ${\mathcal{P}}^k$ is the waveform library at time k. The loss function should be modeled to balance the predicted conditional Bayes risk/cost and measurement cost. Thus, we define it to be the sum of them, that is,

$$L_{\mathcal{C},\mathbf{\Theta }}\left({\mathbf{\theta }}_k|{\mathbf{z}}_{1:k-1}^{};{\mathsf{\Theta }}_{k-1}\right)\equiv C_{k+1}^{}\left({\mathbf{\theta }}_k|{\mathbf{z}}_{1:k-1}^{};{\mathsf{\Theta }}_{k-1}\right)+T_{\mathbf{\Theta }}\left({\mathbf{\theta }}_k\right).$$

(31)

The next ${\mathbf{\theta }}_k^{*}$ is chosen to minimize the loss function, and the choice of the optimal parameter is equivalent to the optimization problem given by:

$${\mathbf{\theta }}_k^{*}=\underset{{\mathbf{\theta }}_k^{}\in \mathrm{P}}{\arg \min }{C}_{k+1}^{}\left({\mathbf{\theta }}_k|{\mathbf{z}}_{1:k-1}^{};{\mathbf{\Theta }}_{k-1}\right) \mathrm{s.t.} {\mathbf{\theta }}_k\in \mathbf{\Theta }.$$

(32)

Standard PF uses the statistical reference, which gives the particles the equal weighting after resampling. CRPF uses the cost reference, which preserves the particle cost after resampling to shift the random grid representation of the cost function toward its local minima. The proposed method is to shift the cost function toward another dimension, namely the vector space effected by the waveform library, to further minimize the cost signal estimates.

3.2.2. Sequential Algorithm and Design Issue

According to Equation (21) and the definition of the risk function, the risk function $\mathcal{R}\left(x_k^{}|{\mathbf{z}}_{k+1}^{};{\mathbf{\theta }}_k\right)$ is a prediction of the cost increment $\Delta C\left(x_{k+1}^{}|{\mathbf{z}}_{k+1};{\mathbf{\theta }}_k\right)$, so $\mathcal{R}(\cdot )$ and $\Delta C(\cdot )$ are related closely, which can be chosen simply as:

$$\Delta C\left(x_k^{}|{\mathbf{z}}_k;{\mathbf{\theta }}_k\right)={\Vert {\mathbf{z}}_k-f_y\left(x_k^{};{\mathbf{\theta }}_k\right)\Vert }^q,$$

(33)

$$\mathcal{R}\left(x_k^{}|{\mathbf{z}}_{k+1}^{};{\mathbf{\theta }}_k\right)={\Vert {\mathbf{z}}_{k+1}-f_y\left(f_x\left(x_k^{};{\mathbf{\theta }}_k\right)\right)\Vert }^q.$$

(34)

After resampling, particles are drawn according to the propagation density, in which we add the waveform parameter ${\mathbf{\theta }}_k$ because of the cognitive framework, denoted as:

$$x_{k+1}^{\left(i\right)}\sim p_{k+1}\left(x_{k+1}^{}|{\widehat{x}}_k^{\left(i\right)};{\mathbf{\theta }}_k\right).$$

(35)

The expected $E[\cdot ]$ with respect to the PDF $p(\cdot )$ is used as the constraint of Equation (35), denoted as:

$$E_{p_{k+1}\left(x_{k+1}^{}|x_k^{};{\mathbf{\theta }}_k\right)}\left[x_{k+1}^{}\right]=f_x\left(x_k^{};{\mathbf{\theta }}_k\right).$$

(36)

Therefore, the particles propagate randomly according to the dynamic model, which changes with ${\mathbf{\theta }}_k$, that is, the conditional pdf $p_{k+1}\left(x_{k+1}^{}|x_k^{}\right)$ still can be arbitrary, but the range of the constraint will be changed by the parameters in another dimension. A zero-mean Gaussian pdf with adaptive variance is used for the particle propagation, denoted as:

$$x_{k+1}^{\left(i\right)}\sim \mathcal{N}\left(f_x\left({\widehat{x}}_{k-1}^{\left(i\right)};{\mathbf{\theta }}_k\right),{\sigma }_k^{2,\left(i\right)}{I}_{L_x}\right),$$

(37)

where ${\sigma }_k^{2,\left(i\right)}$ is the variance, and $I_{L_x}$ is the $L_x\times L_x$ identity function. One choice of adaptive selection for ${\sigma }_k^{2,\left(i\right)}$ is to compute it by the recursive model.

$${\sigma }_k^{2,\left(i\right)}\left({\mathbf{\theta }}_k\right)=\frac{t-1}{t}{\sigma }_{k-1}^{2,\left(i\right)}+\frac{{\Vert x_k^{\left(i\right)}-f_x\left({\widehat{x}}_{k-1}^{\left(i\right)};{\mathbf{\theta }}_k\right)\Vert }^2}{k{L}_x}.$$

(38)

In terms of the cost function:

$$C_{k+1}^{\left(i\right)}\left({\mathbf{\theta }}_k\right)=\lambda {\widehat{C}}_k^{\left(i\right)}+\Delta C_{k+1}^{\left(i\right)}\left({\mathbf{\theta }}_k\right),$$

(39)

where $\Delta C_{k+1}^{\left(i\right)}$ is also affected by the newly updated ${\mathbf{\theta }}_k$, shown as:

$$\Delta C_{k+1}^{\left(i\right)}\left({\mathbf{\theta }}_k\right)=\Delta C\left(x_{k+1}^{\left(i\right)}|{\mathbf{z}}_{k+1}^{};{\mathbf{\theta }}_k\right).$$

(40)

The next transmitted waveform is selected. The function $\mu$ in Equation (41) is selected to assign large probability masses to lower-cost particles. When we add ${\mathbf{\theta }}_k$ to the cost function, the probability masses will be re-assigned as follows:

$${\widehat{\mathbf{\pi }}}_k^{\left(i\right)}\propto {\mu }_2\left\{C_k^{\left(i\right)}\left({\mathbf{\theta }}_k\right)\right\}=\frac{1}{{\left(C_k^{\left(i\right)}\left({\mathbf{\theta }}_k\right)-{\min }_p\left\{C_k^{\left(p\right)}\left({\mathbf{\theta }}_k\right)\right\}+\delta \right)}^{\beta }},$$

(41)

where $0<\delta <1$, $\beta >1$. The current state estimation and its prediction covariance matrix is updated by:

$$i_0=\arg \max {\widehat{\mathbf{\pi }}}_k^{\left(i\right)},{\widehat{x}}_k^{\min }={\widehat{x}}_k^{\left(i_0\right)}.$$

(42)

We can obtain the optimal sequence of state vectors and the pointwise solution, denoted as:

$${\widehat{x}}_k^{mean}={\displaystyle {\sum }_{i=1}^{N_s}{\widehat{\mathbf{\pi }}}_k^{\left(i\right)}{\widehat{x}}_k^{\left(i\right)}}.$$

(43)

The cognitive structure is combined with the cost-reference function couple, considering the internal and external mechanisms to form the complete CCRPF algorithm. The implementation details are given in Algorithm 3. This can be reduced to the CRPF if using a fixed waveform.

Algorithm 3. CCRPF algorithm for maneuvering target tracking problem.

Initialization, $(\mathrm{for} i=1,\dots ,N_s),$

1.    $\mathrm{generate} x_0^{\left(i\right)}\sim p_0\left(x_0\right)$$, \mathrm{assign} \mathrm{the} \cos \mathrm{t} C_0^{\left(i\right)}=0$$, \mathrm{and} \mathrm{initialize} {\mathbf{\sigma }}_0^{2,\left(i\right)}.$

PMF Update

2.    Start with the initial waveform parameter $\mathbf{\theta }$. For each $\mathbf{\theta }$ compute:

3.      $\left\{{\widehat{x}}_{k-1}^{\left(i\right)},{\mathsf{\omega }}_{k-1}^{(i)}\right\}=\left\{x_{k-1}^{\left(i\right)},{\mathsf{\omega }}_{k-1}^{(i)}\right\}$$, i=1,\dots ,N_s$$, {\mathbf{z}}_k=f_y\left(x_k^{},w_k;\mathbf{\theta }\right)$

4.      ${\mathcal{R}}_k^{\left(i\right)}\left(\mathbf{\theta }\right)=\lambda C_{k-1}^{\left(i\right)}+{\Vert {\mathbf{z}}_k-f_y\left(f_x\left(x_{k-1}^{\left(i\right)}\right)\right)\Vert }^q$$, q=1,2$$; i=1,\dots ,N_s$$, {\mathbf{\pi }}_k^{\left(i\right)}\propto \mu \left({\mathcal{R}}_k^{\left(i\right)}\right)=\frac{1}{{\left({\mathcal{R}}_k^{\left(i\right)}-\min {\left\{{\mathcal{R}}_k^{\left(i\right)}\right\}}_{i=1}^{N_s}+\delta \right)}^{\beta }}$

5.      $\mathrm{Resampling} {\widehat{x}}_{k-1}^{}={\left\{{\widehat{x}}_{k-1}^{\left(i\right)},{\widehat{C}}_{k-1}^{\left(i\right)}\right\}}_{i=1}^{N_s}$$\mathrm{according} \mathrm{to} {\mathbf{\pi }}_k^{\left(i\right)}$

Particle Propagation and Waveform Selection

6.       $x_k^{\left(i\right)}\sim p_k\left(x_k^{}|{\widehat{x}}_{k-1}^{\left(i\right)}\right)$$, \mathrm{compute} \mathrm{the} \mathrm{cost} C_k^{\left(i\right)}\left(\mathbf{\theta }\right)=\lambda C_{k-1}^{\left(i\right)}+{\Vert {\mathbf{z}}_k-f_y\left(x_k^{\left(i\right)}\right)\Vert }^q$

7.      $\mathrm{Compute} \mathrm{the} \cos \mathrm{t} \mathrm{function} {\mathbf{\theta }}_k^{*}=\underset{\mathbf{\theta }}{\arg \min }{C}_{\mathbf{\Theta }}^{\left(i\right)}\left({\mathbf{\theta }}_k|{\mathbf{z}}_{1:k-1}^{};{\mathsf{\Theta }}_{k-1}\right)+C_{\mathbf{\Theta }}\left({\mathbf{\theta }}_k\right)$

8.    $\mathrm{Select} \mathrm{the} \mathrm{optimal} \mathrm{waveform} \mathrm{parameter} {\mathbf{\theta }}_{\mathrm{k}}={\mathbf{\theta }}_k^{*}$

Particle and Measurement Recursive Update

9.      ${\widehat{x}}_k^{\left(i\right)}={\widehat{x}}_k^{\left(i\right)}$$, {\mathbf{z}}_k=f_y\left(x_k^{},w_k;{\mathbf{\theta }}_{\mathrm{k}}\right)$$, C_k^{\left(i\right)}=C_k^{\left(i\right)}\left({\mathbf{\theta }}_{\mathrm{k}}\right)$

10.      ${\mathbf{\sigma }}_k^{2,\left(i\right)}\left({\mathbf{\theta }}_{\mathrm{k}}\right)=\{\begin{array}{ll}{\sigma }_{k-1}^{2,\left(i\right)}& \mathrm{t}\le 10\\ \frac{k-1}{k}{\sigma }_{k-1}^{2,\left(i\right)}+\frac{{\Vert x_k^{\left(i\right)}-f_x\left({\widehat{x}}_{k-1}^{\left(i\right)}\right)\Vert }^2}{k\times \dim \left[x\right]}& \mathrm{t}>10\end{array}$, i = 1, …, Ns

Information Update and State Estimation

11.    ${\widehat{\mathbf{\pi }}}_k^{\left(i\right)}\propto {\mu }_2\left(C_k^{\left(i\right)};{\mathbf{\theta }}_{\mathrm{k}}\right)=\frac{1}{{\left(C_k^{\left(i\right)}-\min {\left\{C_k^{\left(i\right)}\right\}}_{i=1}^{N_s}+\delta \right)}^{\beta }},$$\mathrm{where} \alpha ,\beta >0. \mathrm{Normalize} \mathrm{the} \mathrm{PMF}.$

12.       ${\widehat{x}}_k={\widehat{x}}_k^{mean}={\displaystyle {\sum }_{i=1}^{N_s}{\widehat{\mathbf{\pi }}}_k^{\left(i\right)}{\widehat{x}}_k^{\left(i\right)}}$$. \mathrm{Save} \mathrm{the} {\widehat{\mathbf{x}}}_k$$\mathrm{and} {\widehat{\mathbf{P}}}_k^{}\left({\mathbf{\theta }}_{\mathrm{k}}\right).$

3.3. Convergence of CCRPF Algorithm

Due to the change of waveform parameter, resulting in the change of the measurement noise, convergence results regarding CRPFs may be not valid for the proposed algorithm. We assess the convergence. More details about the preliminary definitions and the derivation of the conditions are referred to in [26].

Lemma 2.

Let $\Delta C\left(x_k^{}|y_k^{}\right)$ be a cost function. With a fixed waveform parameter $\mathbf{\theta }$, if the three following conditions are met:

(1) The set function ${\mu }_k\left(A\subseteq {\left\{x_k^{\left(i\right)}\right\}}_{i=1}^M\right)={\displaystyle \sum _{x\in A}^{}\mu \left(\Delta C\left(x_k|y_k^{}\right)\right)}$ satisfies

$$\underset{M\to \infty }{\lim }\mathrm{Pr}\left[1-\frac{{\mu }_k\left(S^M\left(x_k^{opt}\left(\mathbf{\theta }\right),\epsilon \right)\right)}{{\mu }_k\left({\left\{x_k^{\left(i\right)}\left(\mathbf{\theta }\right)\right\}}_{i=1}^M\right)}\ge \delta \right]=0\hspace{1em}\forall \delta >0,$$

(44)

where $\mathrm{Pr}[\cdot ]$ denotes probability, ${\left\{x_k^{\left(i\right)}\right\}}_{i=1}^M$ is a set of particles drawn at time step $k$, and $S^M\left\{x_k^{opt}\left(\mathbf{\theta }\right),\epsilon \right\}=\left\{x\in {\left\{x_k^{\left(i\right)}\right\}}_{i=1}^M:\Vert x-x_k^{opt}\Vert <\epsilon ;\mathbf{\theta }\right\}$.

(2) The mean incremental cost $\overline{\Delta C_k}={\displaystyle \sum _{i=1}^M{\omega }_k^{\left(i\right)}\Delta C\left(x_k^{\left(i\right)}|y_k^{}\right)}$ converges to the minimal incremental cost

$$\underset{M\to \infty }{\lim }\left|\Delta C\left(x_k^{opt}|y_k^{};\mathbf{\theta }\right)-\overline{\Delta C_k}\right|=0\hspace{1em}(i.p.).$$

(45)

(3) The mean cost estimate is asymptotically optimal,

$$\underset{M\to \infty }{\lim }\left|\Delta C\left({\tilde{x}}_k^{mean}|y_k^{};\mathbf{\theta }\right)-\Delta C\left(x_k^{opt}|y_k^{};\mathbf{\theta }\right)\right|=0\hspace{1em}(i.p.),$$

(46)

where${\tilde{x}}_k^{mean}={\displaystyle \sum _{i=1}^M{\mathbf{\pi }}_k^{\left(i\right)}{x}_k^{\left(i\right)}}$, then with the adaptive parameter ${\mathbf{\theta }}_k$, the set function satisfies:

$$\underset{M\to \infty }{\lim }\left[{\mu }_k\left(S^M\left(x_k^{opt}\left({\mathbf{\theta }}_k\right),\epsilon \right)\right)/{\mu }_k\left({\left\{x_k^{\left(i\right)}\left({\mathbf{\theta }}_k\right)\right\}}_{i=1}^M\right)\right]=1\hspace{1em}(i.p.),$$

(47)

and the mean cost estimate is asymptotically optimal, that is:

$$\underset{M\to \infty }{\lim }\left|\Delta C\left({\tilde{x}}_k^{mean}|y_k^{};{\mathbf{\theta }}_k\right)-\Delta C\left(x_k^{opt}|y_k^{};{\mathbf{\theta }}_k\right)\right|=0\hspace{1em}(i.p.).$$

(48)

See Appendix B for a proof.

4. Numerical Results and Discussion

4.1. Dynamic Model

The evolution of the target state and the corresponding measurements are described by a known discrete-time stochastic model separately [47]. Consider the target moving in the $\mathrm{x}-\mathrm{y}$ plane according to the model:

$$x_{k+1}=F_k{x}_k+v_k,$$

(49)

where ${\left[x_{1,k}^{},x_{2,k}^{},x_{3,k}^{},x_{4,k}^{}\right]}^T$ denotes the state vector X_k of the target, $\left(x_{1,k}^{}{x}_{2,k}^{}\right)$, and $\left(x_{3,k}^{},x_{4,k}^{}\right)$ denotes the position and the velocity, respectively. $v_k={\left[0 0 v_x v_y\right]}^T$ is the system noise. We assume that the sampling period is short enough for the velocity to be a constant during the period, so the state transition matrix is defined as:

$$F_k=\left[\begin{array}{ll}{F}_{11}& F_{12}\\ 0& F_{22}\end{array}\right],F_{11}=F_{22}=\left[\begin{array}{l}1 0\\ 0 1\end{array}\right],F_{12}=\left[\begin{array}{l}T 0\\ 0 T\end{array}\right].$$

(50)

The measurement equation is considered to be highly nonlinear, and as described as follows, ${\mathbf{z}}_k={\left[r_k,{\dot{r}}_k,{\beta }_k\right]}^T$ is selected as the observable vector, and the observation matrix is denoted as:

$${\mathbf{z}}_k={\left[r_k, {\dot{r}}_k, {\beta }_k\right]}^T+w_k,$$

(51)

where the range $r_k$, range-rate ${\dot{r}}_k$, and bearing ${\beta }_k$ compose the observable vector, that is:

$$r_k={\left(x_{1,k}^2+x_{2,k}^2\right)}^{1/2},\hspace{1em}{\dot{r}}_k={\left(x_{3,k}^2+x_{4,k}^2\right)}^{1/2},\hspace{1em}{\beta }_k=\arctan \left(x_{2,k}^{}/x_{1,k}^{}\right).$$

(52)

Select Equation (15) to be the model of $R_k$. The fixed-waveform parameter is set as ${\lambda }_0=50\times {10}^{-6}\mathrm{s}$, $b_0=60\times {10}^9\mathrm{yad}/{\mathrm{s}}^2$, and the initial SNR is set as ${\eta }_0=16$, so $R_0=R\left({\lambda }_0,b_0,{\eta }_0\right)$. The initial true state of the target is ${\left[10, -5, -0.2, 0.2\right]}^T$, and the given initial state is set as ${\left[6 -2 0 0\right]}^T$, ${\mathbf{\sigma }}_0^{2,\left(i\right)}={\left[1.5 10 0.1 1.5\right]}^T$. The sampling period is $T=1 \mathrm{s}$. The generally used heavy-tailed distributions include Laplace distribution, t distribution, uniform distribution, and Gaussian distribution with large variance. Considering that the mixture Gaussian is applied, the sum of the different weighted Gaussian noises with different parameters can be used as the modeling of flicker noise $w_k\left({\mathbf{\theta }}_{k-1}\right)$. The pdf of flicker noise can be denoted as:

$$v_{x,y}\sim 0.1\mathcal{N}\left(0,1\right)+4\times {10}^{-4}\mathcal{N}\left(0,1\right)+{10}^{-6}\mathcal{N}\left(0,1\right),$$

(53)

$$w_r\sim \sqrt{0.1}\mathcal{N}\left(0,{\sigma }_r\right)+0.04\times {10}^{-4}\mathcal{N}\left(0,{\sigma }_r\right)+2.5\times {10}^{-5}\mathcal{N}\left(0,{\sigma }_r\right),$$

(54)

$$w_{\beta }\sim 0.1\mathcal{N}\left(0,{\sigma }_{\beta }\right)+4\times {10}^{-4}\mathcal{N}\left(0,{\sigma }_{\beta }\right)+1.6\times {10}^{-7}\mathcal{N}\left(0,{\sigma }_{\beta }\right),$$

(55)

where $\mathcal{N}\left({\mu }_i,{\Sigma }_i\right)$ denotes the Gaussian distribution with the mean of ${\mu }_i$ and the variance of ${\Sigma }_i$. The particle number is 400.

For a particular scenario and parameters, the overall performance of a filter is evaluated using the metric that is the ARMSE. We use it to estimate the tracking accuracy of the algorithms, and it is defined as:

$$\mathrm{ARMSE}=\sqrt{\frac{1}{N_{m}K}{\displaystyle \sum _{i=1}^{N_m}{\displaystyle \sum _{k=1}^K\left[{\left({\widehat{x}}_{1,k}^i-x_{1,k}\right)}^2+{\left({\widehat{x}}_{3,k}^i-x_{3,k}\right)}^2\right]}}},$$

(56)

where $K$ is the total number of time steps and $N_m$ is the total number of independent MC runs. Let $\left(x_{1,k},x_{3,k}\right)$ and $\left({\widehat{x}}_{1,k}^i,{\widehat{x}}_{3,k}^i\right)$ denote the true and the estimated target position at time $k$ at the $i\text{-}th$ MC run, respectively. Then, the RMSE of the estimated states at $k$ can be computed as [48]:

$$\mathrm{RMSE}\left(k\right)=\sqrt{\frac{1}{N_m}{\displaystyle \sum _{i=1}^{N_m}\left[{\left({\widehat{x}}_{1,k}^i-x_{1,k}\right)}^2+{\left({\widehat{x}}_{2,k}^i-x_{2,k}\right)}^2\right]}}.$$

(57)

The upper and lower limits of the waveform parameters are determined according to the transmitter specifications, and the waveform library can be obtained as:

$$\mathcal{P}=\left\{\lambda \in \left[{\lambda }_{\min }:\Delta \lambda :{\lambda }_{\max }\right],b\in \left[b_{\min }:\Delta b:b_{\max }\right]\right\},$$

(58)

where Δλ denotes the step-size of the envelope duration, and Δb denotes the step-size of the chirp rate.

4.2. Selection of Function μ and Parameters δ

Before implementing the simulation experiments, we check the behavior of the proposed method with different function $\mu$ and different values of tuning parameter $\delta$ in Equations (33), (34) and (41), to select suitable $\mu$ and $\delta$.

The selection of the tuning parameters is investigated in Figure 3. Figure 3a shows four cases of mean absolute deviation after the tracking, with varying values set of $\mu$ and $q$. The performance in the case ($\mu 1, q=1$) seems well at beginning but shows a rapid degradation soon. The performances of the other cases are similar, wherein case ($\mu 2, q=2$) shows a slightly better overall performance. Three cases with varying values of $\delta$ are considered in Figure 3b. It is clear that a small $\delta$ may result in a worse improvement. Therefore, it is reasonable to select $\mu 2, q=2, \delta =1$ for the rest of the simulations according to the analysis and conclusion.

4.3. Simulation Results

4.3.1. Scenario 1: Unknown Statistics in General Environment

We assume that the process noise and measurement noise is independent and temporally white in this case, with the mixture Gaussian pdf as Equations (53)–(55) show, but the given methods are mismatched with the dynamic systems and models.

We apply SPF, CRPF, and CCRPF for target tracking in a two-dimensional space for performance assessment. Figure 4 displays the true path against the tracks of these filters. The other two tracks the target with an error smaller than the SPF, wherein the CCRPF outperforms CRPF slightly. Note how the trajectory of CCRPF deviates from the true path slightly due to the cost function and cognitive structure. EKF does not show the performance of robustness, and we do not show the result here.

All error curves corresponding to the above three filters were obtained by simulation runs. The results for the estimation error of the range and range-rate are shown in Figure 5a,b respectively. Observe from Figure 5 that the SPF introduces large bias in the estimation. CRPF and CCRPF are unbiased throughout the observation period. As for the comparison of CRPF and CCRPF, the overall performance of CCRPF is considered better than that of CRPF. Particularly, the CCRPF estimation is approaching zero bias except for the very beginning stage, indicating a fairly robust performance.

In term of the convergence speed, the convergence of CRPF and CCRPF shows more rapidly than SPF and keeps stable; thus, they are more robust than SPF in tracking performance. The CRPF converges the error asymptotically, while the proposed method achieves convergence in less than 10 time steps, due to the real-time perception of environment in the cognitive algorithm.

CRPF does not have a bad overall performance, but it provides higher peak errors than CCRPF since the cost reference mechanism can only use the fixed waveform, that is to say, it can only be adaptive in the constant environment conceived by the fixed waveform. However, cognitive radar based on CCRPF can not only perceive the dynamic environment by the dynamic waveform, but can also stay adaptive during the recursion after the waveform is chosen. Thus, CCRPF can suppresses the peak error on the fly and limits the overall magnitude of the error, showing the superiority.

A summary of these results and the detailed comparison are also tabulated in

Table 1. Performance comparison of algorithms.

ARMSE	SPF-1 1	SPF-2 2	CRPF-1 3	CRPF-2 4	CCRPF
Range	1.54	1.09	1.26	0.96	0.93
Range-rate	0.37	0.24	0.17	0.16	0.14

¹ SPF-1: SPF with fixed waveform 1, ² SPF-2: SPF with fixed waveform 2, ³ CRPF-1: CRPF with fixed waveform 1, ⁴ CRPF-2: CRPF with fixed waveform 2.

. The parameter of fixed waveform 1 is $\lambda =10\times {10}^{-6}$, $b=10\times {10}^9$, and the parameter of fixed waveform 2 is $\lambda =50\times {10}^{-6}$, $b=60\times {10}^9$. The first row of the table lists the metrics corresponding to the ranging performance by MC simulations with 40 independent runs, which is computed using Equation (56). The second row gives the performance on estimation of range-rate.

From this table, more concretely, we see that the SPF shows a degradation in performance, it has the lowest tracking accuracy of the compared methods. From the comparison of PFs using two different fixed waveforms, we find that PF with waveform 2 is better than that with waveform 1, due to the effects of the measurement noise obtained by different waveforms, and the same conclusion can be found in CRPFs. It is clear that the best filters for this case was CCRPF, which achieves an excellent performance of a 15% improvement over the SPF-2 and a 3% improvement over the CRPF-2.

The unknown environment has a huge influence on SPF for tracking. The mismatch of the noise is the main cause result in the performance degradation. CRPF is less affected since it takes advantage of the cost reference mechanism, so it has better tracking accuracy than SPF. The mastery of environmental information in time has an extra contribution in the proposed algorithm, which has more obvious performance improvement compared with others because of the robust module working cooperatively with the cognitive structure. This also indicates that it is a good supplement to the cognitive radar tracking approach.

To study how the cognition process evolves across time, we have plotted the waveform selection for both the chirp rate and the duration of the pulse envelope in Figure 6a,b. We observe that the transition of the chirp rate is switched from maximum up-sweep to maximum down-sweep.

4.3.2. Scenario 2: Unknown Statistics in Dynamic Noise Environment

Unknown and abrupt change of noise would be considered in this case. The main parameters are the same with Scenario 1. Measurement noise covariance $R$ is known a priori, and the actual process noise covariance $Q$ is set as:

$$\{\begin{array}{ll}{Q}_0=Q^{\circ },& \mathrm{k}=1,\dots ,60\\ Q_0=80\times Q^{\circ },& \mathrm{k}=61,\dots ,80\end{array}.$$

(59)

Figure 7 shows the RMSE curves corresponding to the four filters. When $\mathrm{k}=1,\dots ,60$ with a priori known noise, all algorithms can track the target successfully. In the subsequent phase when $Q$ is increasing from the time step 60 abruptly, although SPF could not be seen as a failure to track the target, its convergence effect and speed are all the worst. CRPF shows a similar convergence speed to CCRPF, because of the similar robust module, but its overall performance during the dynamic noise stage is obviously inferior to the latter. The proposed approach shows the superiority again, with the lowest error and fastest convergence speed.

A summary of the results and the detailed comparison are also tabulated in

Table 2. Performance comparison of algorithms.

ARMSE	SPF	CRPF	CCRPF
$k\le 30$	0.51	0.58	0.23
$k>30$	1.49	0.20	0.19

. We do not use a different waveform in SPF and CRPF anymore because similar work has already been done in Scenario 1.

5. Conclusions

We have concentrated on the tracking method of the cognitive radar based on particle filter. The CCRPF for cognitive radar is proposed as a significant step toward random dynamic systems with unknown statistics. In this work, the mathematical model of cognitive radar has been derived, CRLB and the corresponding cost function have been designed, the cognitive PF algorithm has been presented by completely reconstructing the propagation and update process, and the proposed algorithm has been developed by updating the cost function and noise variance. Moreover, the convergence of the approach has been proofed. The simulation results illustrate that the state error prediction is more adjacent to the CRLB, and the proposed method showed a good outperforming result over the existing methods in accuracy and robustness on highly nonlinear dynamic systems with unknown statistics or a complicated environment. The application of cognitive radar and the adaptation of a traditional filter have been expanded to a wider scope.

Future work: (1) PF is approximately globally optimal, while cognitive estimation is biased, so the closed-form of the optimal solution is not easy to obtain in the PF-based cognitive tracking method, and we can only use MC to get the solution or proof the convergence. However, the cost function in CRPF might have the optimal solution; thus, the cost function design is deserves to have further study. (2) Using novel SMC methods as the cognitive framework may have better prospects in dealing with dynamic system problems with unknown or uncertain statistics. (3) The lower limit and dynamic change of particle number could also be determined.