1 Introduction

As a result of the closed-loop feedback from the receiver to the transmitter, the cognition radar [10] adapts its waveform to the dynamic environment, to enhance the tracking performance at the system level. Therefore, this issue has received significant research in recent years [11, 12].

The waveform agility in tracking originated from [15], where the Kalman filter was used for single target tracking. In addition, the Cramér–Rao lower bound (CRLB) on the measurement error covariance was derived from the Fisher information of the transmitted waveform. Therefore, the transmitted waveform and the tracking filter were connected. References [16, 17] extended this work to the clutter environment; however, the linear observation model was still utilized. The impacts of different waveforms on tracking performance were studied in [19, 21, 22, 32], but the dynamic waveform configuration was not considered. The problem of the joint beam and waveform scheduling was investigated in [5, 23]. One-step-ahead and two-step-ahead algorithms were proposed to select the waveform and sample interval. To optimize the detection threshold and the transmitted waveform jointly, reference [13] choses the cumulative probability of track loss and the state covariance as the cost function. In [14, 24], the fractional Fourier transform was exploited to rotate the ambiguity function (AF) of the transmitted waveform. Therefore, the waveform library was established. Additionally, the interacting multiple model (IMM) algorithm was used for the maneuvering target tracking. References [25, 26, 29] formed the frequency-modulated (FM) waveform library and used the particle filter (PF) to deal with the nonlinear measurement. References [6, 7] added the LFM library utilizing the IMM probabilistic data association filter (IMM-PDAF) and achieved urban terrain tracking in high clutter. Reference [9] presented a novel Kalman filter, which was embedded into the extended kernel recursive least squares Kalman filter (Ex-KRLS-KF) algorithm to further improve the tracking performance. The proposed algorithm improved the tracking performance effectively compared to the state-of-the-art algorithms. Reference [31] proposed a new diffusion sign subband adaptive filtering algorithm with an individual weighting factor (IWF-DSSAF) for distributed estimation in the impulsive noise environment, which achieved better convergence performance than their counterparts. Reference [35] concerned the application of Huber-based robust unscented Kalman filter (HRUKF) in a nonlinear system. An adaptive strategy was proposed to improve filtering performance with non-Gaussian measurement noise, which had a better performance than the traditional ones. In reference [27], an adaptive kernel Kalman filter-based belief propagation algorithm is presented to tracking targets in the case of clutter and false alarms. The proposed algorithm has a better tracking performance and lower computation cost compared with other algorithms. Reference [28] enhances the performance of the interactive multiple model-integrated probabilistic data association algorithm (IMM-IPDA) by the fixed lag smoothing algorithm. Compared with the other recent algorithms in the literature, the algorithm is better in the root-mean-square error (RMSE), true track rate (TTR), and mode probabilities. Reference [8] introduces a novel tracking algorithm, which can significantly reduce the estimation error during non-maneuvering periods. Therefore, it is very suitable for tracking low maneuvering targets.

While the previous research has made seminal contributions to the tracking problem and waveform selection, there are still some issues to be addressed.

  1. (1)

    Most research [4,5,6,7, 13, 15,16,17, 23, 24] has focused on the linear system. As a result, the nonlinear observation model, in particular the maneuvering target tracking in the clutter environment, has received minimal attention.

  2. (2)

    Though [25, 26, 29] discussed the tracking problem in the nonlinear system and the clutter environment, the PF was used for dealing with the nonlinear transform. This approach resulted in high computational complexity and an impractical application.

  3. (3)

    More importantly, to select the optimal waveform, each parameter needed to be traversed in [5,6,7, 11, 15,16,17, 23, 25, 26, 29]. These traversals inevitably resulted in a heavier computational burden.

Therefore, this paper considers solving the maneuvering target tracking problem in clutter by utilizing an algorithm with a simpler structure and improved accuracy. Based on the CS model, which was proposed in [33], the modified PDAF [30] and the SCKF [1] are utilized to handle the clutter and the nonlinear transform. Additionally, a novel waveform agile scheduling method is proposed to structure the adaptive MCS-MPDA-SCKF. The main aspects of this paper are summarized as follows.

  1. (1)

    The MPDAF is utilized to handle the clutter. It corrects the state prediction covariance in the event of missing detection in the PDAF [3, 4], to achieve a much lower estimation error. Simultaneously, and different from the PDAF, the Riccati equation is the accurate solution and thus is directly utilized for waveform optimization.

  2. (2)

    The square-root cubature filter (SCKF), which is based on the third-degree spherical–radial cubature rule, is utilized for high estimation precision as well as low computational complexity. As a result, it is integrated with the MPDAF as the tracker.

  3. (3)

    A waveform scheduling algorithm, which is based on the fractional Fourier transform (FRFT), is designed. By exploiting the dynamic waveform selection, not only the posterior state error is minimized, but also traversal in the waveform library is avoided. As a result, both the performance and the efficiency are improved.

  4. (4)

    The simulation results showed that, compared with the algorithms with fixed waveform, the proposed algorithm achieved a lower track loss percentage as well as higher tracking precision. At the same time, the proposed algorithm possesses the advantages of simpler structure and improved accuracy compared with the two state-of-the-art algorithms.

The rest of this paper is organized as follows. The state and measurement model, the waveform model, and the clutter model are described in Sect. 2. In Sect. 3, the MPDAF is integrated with the SCKF to structure the MPDA-SCKF. A waveform library based on the FRFT is established, and a direct waveform selection method is presented in Sect. 4. Extensive simulation results to verify the effectiveness of the proposed algorithm are shown in Sect. 5. Section 6 contains the conclusions.

2 Modeling

2.1 Waveform Model

The typical transmitted waveform in a narrowband environment is as follows:

$$ s_{T} (t) = \sqrt 2 {\text{Re}} \left[ {\sqrt {E_{T} } \tilde{s}(t)\exp (j2\pi f_{c} t)} \right] $$
(1)

where ET is the waveform energy, and fc is the carrier frequency. \(\tilde{s}(t)\) is the complex envelope of the transmitted pulse.

The received waveform model is as follows:

$$ s_{R} (t) = \sqrt 2 {\text{Re}} \left[ {\sqrt {E_{R} } \tilde{s}(t - \tau - \frac{{2\dot{r}}}{c}t)\exp (j2\pi f_{c} (t - \frac{{2\dot{r}}}{c}t))} \right] + n(t) $$
(2)

where ER is the received energy, and n(t) is additive white Gaussian noise. \(\tau = 2r/c\), the time delay. \(r\) is the range between radar and target. \(\dot{r}\) is the range rate, and \(\dot{r} < < c\) . When the time-bandwidth product of the waveform satisfies the narrowband condition, sR(t) can be approximated by:

$$ s_{R} (t) \approx \sqrt 2 {\text{Re}} \left[ {\sqrt {E_{R} } \tilde{s}(t - \tau )\exp (j2\pi f_{c} (t - \frac{{2\dot{r}}}{c}t))} \right] + n(t) $$
(3)

2.2 State and Measurement Model

The MCS model [33] and the nonlinear measurement originated from the target are as follows:

$$ \left[ \begin{gathered} {\mathbf{X}}_{k + 1} \hfill \\ {\overline{\mathbf{a}}}_{k + 1}^{(1)} \hfill \\ \end{gathered} \right] = \left[ {\begin{array}{*{20}c} {{\mathbf{F}}_{{{\text{ACS}}}} } & {{\mathbf{U}}_{{{\text{ACS}}}} } \\ {\mathbf{0}} & {\mathbf{I}} \\ \end{array} } \right]\left[ \begin{gathered} {\mathbf{X}}_{k} \hfill \\ {\overline{\mathbf{a}}}_{k}^{(1)} \hfill \\ \end{gathered} \right] + \left[ {\begin{array}{*{20}c} {{\mathbf{W}}_{k} } \\ {\mathbf{0}} \\ \end{array} } \right] $$
(4)
$$ {\mathbf{z}}_{k} = h({\mathbf{X}}_{k} ) + {\mathbf{V}}_{k} $$
(5)

where \({\overline{\mathbf{a}}}_{k}^{(1)}\) is the first derivative acceleration mean. \(h(.)\) is the nonlinear transformation function. zk is the measurement matrix from the target. \({\mathbf{W}}_{k} \sim ({\mathbf{0}},{\mathbf{Q}}_{k} )\) is the process noise. \({\mathbf{V}}_{k} \sim ({\mathbf{0}},{\mathbf{R}}_{k} )\) is the measurement noise. They are mutually independent. Qk = 2ασa2qcs. For more details, please refer to [33]. Assuming that the target moves in a two-dimensional plane, and the range, range rate, as well as orientation are measured simultaneously. As a result, \({\mathbf{X}}_{k} = \left[ {x_{k} ,\dot{x}_{k} ,\ddot{x}_{k} ,y_{k} ,\dot{y}_{k} ,\ddot{y}_{k} } \right]^{{\text{T}}}\),\({\mathbf{z}}_{k} = \left[ {r_{k} ,\dot{r}_{k} ,\theta_{k} } \right]^{{\text{T}}}\) where \(\theta_{k} = \arctan (y/x)\) and \(\dot{r}_{k} = (\dot{x}_{k} x_{k} + \dot{y}_{k} y_{k} )/r_{k}\). FACS and UACS are as follows:

$$ {\mathbf{F}}_{{{\text{ACS}}}} = \left[ {\begin{array}{*{20}c} 1 & T & {T^{2} /2} \\ 0 & 1 & T \\ 0 & 0 & 1 \\ \end{array} } \right] $$
(6)
$$ {\mathbf{U}}_{{{\text{ACS}}}} = \left[ {\begin{array}{*{20}l} {T^{3} /6 - (2 - 2\alpha T + \alpha^{2} T^{2} - 2e^{ - \alpha T} )/2\alpha^{3} } \hfill \\ {T^{2} /2 - (e^{ - \alpha T} - 1 + \alpha T)/\alpha^{2} } \hfill \\ {T - (1 - e^{ - \alpha T} )/\alpha } \hfill \\ \end{array} } \right] $$
(7)

where T is the sample interval. Note that there is a close relationship between Rk and εk, which denotes the selected waveform (or waveform parameters) in time k. (It is shown in Sect. 4). Therefore, for a specific waveform, a corresponding Rk can be denoted as Rk (εk).

2.3 Clutter Model

Owing to the clutter, the measurement in time k not only consists of the target echo but also includes false alarms. As a result, the measurement is given by:

$$ {\mathbf{Z}}_{k} = \left[ {z_{k}^{1} ,z_{k}^{2} ,...,z_{k}^{i} ,...z_{k}^{{m_{k} }} } \right] $$
(8)

where mk is the total number of measurements obtained by the radar in time k. Each zik consists of the range, the range rate, and the angle information. Assuming that the number of false alarms follows the Poisson distribution with the average ρVk, where ρ is the density of false measurement, and Vk is the validation gate volume. The probability mass function of false measurement is as follows:

$$ \mu (m_{k} ) = \frac{{\exp ( - \rho V_{k} )(\rho V_{k} )^{{m_{k} }} }}{{m_{k} !}} $$
(9)

Assuming that the clutter is uniformly distributed in the validation gate. Since the noise is additive, the test statistics according to the noise-only and target-present hypotheses follow an exponential distribution [25]. Simultaneously, the peak value of the ambiguity function is used to estimate the time delay and Doppler shift. As a result, the probability of detection in time k is as follows:

$$ P_{dk} = P_{f}^{{1/(1 + \eta_{k} )}} $$
(10)

where Pf is the desired probability of false alarms, and ηk is the SNR in time k.

3 MPDA-SCKF

3.1 MPDAF

The conventional PDAF [3, 4] calculates the probability of each measurement being from the target to determine its weight. Then, the sum of weighted measurements is regarded as the equivalent measurement and used for the state update. However, owing to the influence of clutter, the number of valid measurements is uncertain. In the PDAF, the current measurement is a priority for updating the state error covariance. As a result, the prediction error covariance cannot be directly utilized for the waveform optimization in the next time interval. To solve this problem, reference [17] utilized the modified Riccati equation as the state error covariance in the next time interval. For the PDAF, the modified Riccati equation was utilized as the approximated solution in the event of a large association area. By contrast, for the MPDAF, the modified Riccati equation is the accurate solution. Additionally, there is no need for a large association area. As a result, the MPDAF is chosen as the update tracker.

Assuming that βik is the associated probability that the measurement zik in time k is from the target. β0k is the associated probability that no measurement is from the target. The posterior state error covariance in MPDAF is as follows:

$$ {\mathbf{P}}_{k + 1|k + 1} = {\mathbf{P}}_{k + 1|k} - \left[ {1 - (\alpha + 1)\beta_{k}^{0} } \right]{\mathbf{K}}_{k + 1} {\mathbf{S}}_{k + 1} {\mathbf{K}}_{k + 1}^{{\text{T}}} + {\tilde{\mathbf{P}}}_{k + 1} $$
(11)

where \({\mathbf{K}}_{k + 1} = {\mathbf{P}}_{k + 1|k} {\mathbf{H}}^{{\text{T}}} {\mathbf{S}}_{k + 1}^{ - 1}\) is the Kalman gain. \({\mathbf{S}}_{k + 1} = {\mathbf{HP}}_{k + 1|k} {\mathbf{H}}^{{\text{T}}} + {\mathbf{R}}_{k + 1}\) is the measurement innovation covariance. \({\mathbf{P}}_{k + 1|k}\) is the prediction covariance. \({\tilde{\mathbf{P}}}_{k + 1}\) is the mean dispersion. α is an influencing factor of the filtering error covariance:

$$ \alpha = \frac{{P_{d} P_{g} (1 - c_{\Gamma } )}}{{1 - P_{d} P_{g} }} $$
(12)

where Pd is the detection probability, and Pg is the gate probability. cΓ implies the influence of the association gate on the innovation covariance Sk+1. When the measurement dimension nz = 3:

$$ c_{\Gamma } = 1 - \frac{{2(\gamma /2)^{3/2} e^{ - \gamma /2} }}{{3\left[ {(\sqrt \pi /2){\text{erf}}\left[ {(\gamma /2)^{1/2} } \right] - (\gamma /2)^{1/2} e^{ - \gamma /2} } \right]}} $$
(13)

where γ is the association threshold. γ = g2. g is the association area. \({\text{erf}}(x) = (2/\sqrt \pi )\int_{ - \infty }^{x} {e^{{ - y^{2} }} dy}\), which is defined as an error function. Compared with the traditional PDAF, α is additional in the MPDAF, which is the result of correctly considering the state prediction error covariance.

It is noticeable that since Rk+1 corresponds to a specific waveform εk+1 (or waveform parameter), Pk+1|k+1 is also related to εk+1, namely, it can be denoted as Pk+1|k+1(εk+1). Aiming to improve the performance, how to dynamically select the optimal waveform is paramount. For the Kalman filter, the posterior state error covariance is obtained by the transmitted waveform [15]. However, Eq. (11) shows that Pk+1|k+1 cannot be directly obtained by the primary transmitted waveform, since β0k and \({\tilde{\mathbf{P}}}_{k + 1}\) are dependent on the measurement set. Owing to the influence of clutter, the measurement set in time k + 1 is uncertain and random. However, reference [30] analyzed the stable performance of MPDAF by the modified Riccati equation. The stochastic items β0k and \({\tilde{\mathbf{P}}}_{k + 1}\) are replaced by their expectations so that the influence of measurement is eliminated. When the modified Riccati equation is adopted to estimate Pk+1|k+1

$$ {\mathbf{P}}_{k + 1|k + 1} = {\mathbf{P}}_{k + 1|k} - q_{2} (\rho V_{k + 1} ,P_{{{\text{d}}k + 1}} ){\mathbf{K}}_{k + 1} {\mathbf{S}}_{k + 1} {\mathbf{K}}_{k + 1}^{{\text{T}}} $$
(14)

where q2 is a scalar between 0 and 1, which depends on the clutter density ρ, the association gate Vk+1, and the probability of detection Pdk+1 in time k + 1. q2 is approximately fitted by [18]:

$$ q_{2} (\rho V_{k + 1} ,P_{dk + 1} ) = \frac{{0.9932P_{d} }}{{1 + (0.44P_{d}^{ - 0.9868} - 0.1696P_{d} )\lambda V_{k + 1} }} $$
(15)

when nz = 3 and Vk+1 = (4π/3)g3|Sk+1|1/2.

3.2 MPDAF-SCKF

The SCKF [30], which transformed the nonlinear filtering problem with a Gaussian distribution into an integration calculation by utilizing a three-degree spherical–radical rule, was adopted for various estimation problems. Moreover, it avoided the square-root operation to the state error covariance matrix and achieved improved estimation accuracy. As a result, the SCKF was used to handle the nonlinear transformation problem and was integrated with the MPDAF to create the MPDA-SCKF.

Before the procedure of MPDA-SCKF, the following definition was needed: R is the upper triangular matrix obtained from the QR decomposition on matrix AT; therefore, the QR decomposition operation to A can be expressed as S = Tria(A) = RT, where S is the lower triangular matrix. The procedure of MPDA-SCKF can be described as follows [1, 33]:

Step 1 Time update

Assume that state error covariance is Pk|k in time kT, then.

Step 1.1 Factorize

$$ {\mathbf{P}}_{k|k} = {\mathbf{S}}_{k|k} ({\mathbf{S}}_{k|k} )^{{\text{T}}} $$
(16)

Step 1.2 Evaluate the cubature points and the propagated cubature points

$$ {\mathbf{X}}_{k|k}^{i} = {\hat{\mathbf{X}}}_{k|k} + {\mathbf{S}}_{k|k} {{\varvec{\upxi}}}_{i} ,i = 1,2...m $$
(17)
$$ {\mathbf{X}}_{k + 1|k}^{i*} = f({\mathbf{X}}_{k|k}^{i} ) $$
(18)

where Xik|k and Xi*k+1|k are the cubature point and the propagated cubature point, respectively. m is the total number of cubature points, and n is the dimensional number of the state vector X, which satisfies m = 2n. m = 2n. ξi = (m/2)1/2[1, 1] is the full permutation and reverse of the n-dimensional unit vector:

$$ [{\mathbf{1}}] = \left\{ {\left( \begin{gathered} 1 \hfill \\ 0 \hfill \\ ... \hfill \\ 0 \hfill \\ \end{gathered} \right),\left( \begin{gathered} 0 \hfill \\ 1 \hfill \\ ... \hfill \\ 0 \hfill \\ \end{gathered} \right),...,\left( \begin{gathered} 0 \hfill \\ 0 \hfill \\ ... \hfill \\ 1 \hfill \\ \end{gathered} \right),\left( \begin{gathered} - 1 \hfill \\ 0 \hfill \\ ... \hfill \\ 0 \hfill \\ \end{gathered} \right),\left( \begin{gathered} 0 \hfill \\ - 1 \hfill \\ ... \hfill \\ 0 \hfill \\ \end{gathered} \right),...,\left( \begin{gathered} 0 \hfill \\ 0 \hfill \\ ... \hfill \\ - 1 \hfill \\ \end{gathered} \right)} \right\} $$
(19)

Step 1.3 Estimate the predicted state and the square-root factor of the predicted error covariance.

$$ {\hat{\mathbf{X}}}_{k + 1|k} = \frac{1}{m}\sum\limits_{i = 1}^{m} {{\mathbf{X}}_{k + 1|k}^{i*} } $$
(20)
$$ {\mathbf{S}}_{k + 1|k} = {\text{Tria}}([{\mathbf{X}}_{k + 1|k}^{*} ,{\text{Chol}}({\mathbf{Q}}_{k} )]) $$
(21)

where the weighted, centered matrix

$$ {\mathbf{X}}_{k + 1|k}^{*} = \frac{1}{\sqrt m }\left[ {{\mathbf{X}}_{k + 1|k}^{1*} - {\hat{\mathbf{X}}}_{k + 1|k} ,{\mathbf{X}}_{k + 1|k}^{2*} - {\hat{\mathbf{X}}}_{k + 1|k} ,...,{\mathbf{X}}_{k + 1|k}^{m*} - {\hat{\mathbf{X}}}_{k + 1|k} } \right] $$
(22)

Step 2 Measurement update

Step 2.1 Evaluate the cubature points and the propagated cubature points

$$ {\mathbf{X}}_{k + 1|k}^{i} = {\hat{\mathbf{X}}}_{k + 1|k} + {\mathbf{S}}_{k + 1|k} {{\varvec{\upxi}}}_{i} $$
(23)
$$ {\mathbf{Z}}_{k + 1|k}^{i} = h({\mathbf{X}}_{k + 1|k}^{i} ) $$
(24)

Step 2.2 Estimate the predicted measurement

$$ {\hat{\mathbf{Z}}}_{k + 1|k} = \frac{1}{m}\sum\limits_{i = 1}^{m} {{\mathbf{Z}}_{k + 1|k}^{i} } $$
(25)

Step 2.3 Estimate the square-root of the residual error (innovation) covariance matrix

$$ {\mathbf{S}}_{k + 1|k}^{ZZ} = {\text{Tria}}([{\mathbf{Z}}_{k + 1|k} ,{\text{Chol}}({\mathbf{R}}_{k + 1} )]) $$
(26)

where the weighted, centered matrix

$$ {\mathbf{Z}}_{k + 1|k} = \frac{1}{\sqrt m }\left[ {{\mathbf{Z}}_{k + 1|k}^{1} - {\hat{\mathbf{Z}}}_{k + 1|k} ,{\mathbf{Z}}_{k + 1|k}^{2} - {\hat{\mathbf{Z}}}_{k + 1|k} ,...,{\mathbf{Z}}_{k + 1|k}^{m} - {\hat{\mathbf{Z}}}_{k + 1|k} } \right] $$
(27)

Step 2.4 Calculate the residual error covariance matrix

$$ {\mathbf{P}}_{k + 1|k}^{ZZ} = {\mathbf{S}}_{k + 1|k}^{ZZ} ({\mathbf{S}}_{k + 1|k}^{ZZ} )^{{\text{T}}} $$
(28)

Step 2.5 Estimate the cross-covariance matrix between the state prediction error vector and the residual error vector

$$ {\mathbf{P}}_{k + 1|k}^{XZ} = {\mathbf{X}}_{k + 1|k} ({\mathbf{Z}}_{k + 1|k} )^{{\text{T}}} $$
(29)

where the weighted, centered matrix

$$ {\mathbf{X}}_{k + 1|k} = \frac{1}{\sqrt m }\left[ {{\mathbf{X}}_{k + 1|k}^{1} - {\hat{\mathbf{X}}}_{k + 1|k} ,{\mathbf{X}}_{k + 1|k}^{2} - {\hat{\mathbf{X}}}_{k + 1|k} ,...,{\mathbf{X}}_{k + 1|k}^{m} - {\hat{\mathbf{X}}}_{k + 1|k} } \right] $$
(30)

Step 2.6 Estimate the Kalman gain

$$ {\mathbf{K}}_{k + 1} = {\mathbf{P}}_{k + 1|k}^{XZ} ({\mathbf{P}}_{k + 1|k}^{ZZ} )^{ - 1} $$
(31)

Step 2.7 Estimate the square-root factor of the corresponding error covariance matrix and the posterior estimation error covariance matrix

$$ {\mathbf{S}}_{k + 1|k + 1}^{c} = {\text{Tria}}([{\mathbf{X}}_{k + 1|k} - {\mathbf{K}}_{k + 1} {\mathbf{Z}}_{k + 1|k} ,{\mathbf{K}}_{k + 1} {\text{Chol}}({\mathbf{R}}_{k + 1} )]) $$
(32)
$$ {\mathbf{P}}_{k + 1|k + 1} = {\mathbf{S}}_{k + 1|k + 1}^{c} ({\mathbf{S}}_{k + 1|k + 1}^{c} )^{{\text{T}}} - q_{2} {\mathbf{K}}_{k + 1} {\mathbf{P}}_{k + 1|k}^{ZZ} ({\mathbf{K}}_{k + 1} )^{{\text{T}}} $$
(33)

Step 2.8 Update the state matrix and the corresponding error covariance matrix

$$ {\hat{\mathbf{X}}}_{k + 1|k + 1} = {\hat{\mathbf{X}}}_{k + 1|k} $$
(34)
$$ {\hat{\mathbf{X}}}_{k + 1|k + 1} = {\hat{\mathbf{X}}}_{k + 1|k} + {\mathbf{K}}_{k + 1} {{\varvec{\upnu}}}_{k + 1} $$
(35)
$$ {{\varvec{\upnu}}}_{k + 1} = \sum\limits_{i = 1}^{{m_{k} }} {\beta_{k + 1}^{i} {{\varvec{\upnu}}}_{k + 1}^{i} } $$
(36)
$$ {\mathbf{S}}_{k + 1|k + 1} = {\text{Chol}}({\mathbf{P}}_{k + 1|k + 1} ) $$
(37)

If no measurement is correct, Eq. (34) and Eq. (37) are used for the update; otherwise, Eq. (35), Eq. (36), and Eq. (37) are used for the update. For more information on βik+1 and νk+1, please refer to [3, 4].

4 Waveform Scheduling Through the FRFT

4.1 CRLB of Measurement Error

As mentioned before, the measurement noise covariance matrix Rk+1 is related to the transmitted waveform εk+1 in time k + 1 (the selected waveform in time k). Assuming that A(τ,fd) is AF of transmitted waveform s(t), then,

$$ A(\tau ,f_{{\text{d}}} ) = \int_{ - \infty }^{ + \infty } {s(t - \frac{\tau }{2})} s^{*} (t + \frac{\tau }{2})e^{{j2\pi f_{d} t}} dt $$
(38)

The Fisher information matrix for time delay τ and Doppler shift fd is as follows [15]:

$$ J = \left. { - \eta \left[ {\begin{array}{*{20}c} {\frac{{\partial^{2} A}}{{\partial \tau^{2} }}} & {\frac{{\partial^{2} A}}{{\partial \tau \partial f_{d} }}} \\ {\frac{{\partial^{2} A}}{{\partial f_{d} \partial \tau }}} & {\frac{{\partial^{2} A}}{{\partial f_{d}^{2} }}} \\ \end{array} } \right]} \right|_{{\tau = 0,f_{d} = 0}} $$
(39)

where η is the SNR. The relationship between Rk+1 and A(τ,fd) is as follows:

$$ {\mathbf{R}}_{k + 1} ({{\varvec{\upvarepsilon}}}_{k + 1} ) = {\mathbf{TJ}}^{ - 1} ({{\varvec{\upvarepsilon}}}_{k + 1} ){\mathbf{T}} $$
(40)

where T = diag(c/2,c/(2fc)), the transformation matrix of τ and fd to the range and the range rate. J−1(εk+1) is the CRLB on the selected waveform for unbiased estimation. Equation (40) shows that the selected waveform in time k determines Rk+1. In addition, Rk+1 is related to the posterior estimation error. Therefore, the optimal waveform should be selected to minimize the estimation error in time k + 1 [11]:

$$ {{\varvec{\upvarepsilon}}}_{k + 1}^{*} = \arg \min \det ({\mathbf{P}}_{k + 1|k + 1} ) $$
(41)

Equation (38) ~ Eq. (41) show that the selected waveform in time k affects the measurement noise covariance matrix, as well as the state estimation error in time k + 1. As a result, by utilizing the waveform scheduling, the optimal waveform was obtained. Simultaneously, the performance of target tracking in clutter was further improved.

4.2 Optimal Waveform Selection Based on FRFT

Assuming that the base transmitted waveform (e.g., a rectangular pulse) is s0(t), the AF is A0(τ,fd), the Fisher information matrix is J0, and the corresponding measurement noise covariance matrix is R0. The FRFT [2, 20], where the fractional factor is φk+1, was applied to the base transmitted waveform to achieve the orthogonality between the measurement error ellipse and the state prediction error ellipse and to satisfy Eq. (41).

The FRFT was regarded as a rotation operation to the coordinate system. When the FRFT was applied to a transmitted waveform with φk+1, the AF of the waveform was rotated by φk+1, and a novel waveform was obtained. Samples of waveforms obtained by φk+1 are shown in Fig. 1, where the base waveform was a rectangular pulse.

Fig. 1
figure 1

a Waveform when φk + 1 = 0.1π. b Waveform when φk + 1 = 0.6π. c Waveform when φk + 1 = 0.9π

Full size image

The obtained waveform is characterized by

$$ {\mathbf{J}}_{k + 1} = [{\mathbf{L}}(\varphi_{k + 1} )]^{T} {\mathbf{J}}_{0} {\mathbf{L}}(\varphi_{k + 1} ) $$
(42)
$$ {\mathbf{R}}_{k + 1} = [{\mathbf{L}}(\varphi_{k + 1} )]^{T} {\mathbf{R}}_{0} {\mathbf{L}}(\varphi_{k + 1} ) $$
(43)

where Jk+1 and Rk+1 are the Fisher information matrix and the covariance matrix of measurement error after the rotation, respectively. L(φk+1) is the rotation matrix, satisfying: \({\mathbf{L}}(\varphi_{k + 1} ) = \left[ {\begin{array}{*{20}c} {\cos \varphi_{k + 1} } & { - \sin \varphi_{k + 1} } & 0 \\ {\sin \varphi_{k + 1} } & {\cos \varphi_{k + 1} } & 0 \\ 0 & 0 & 1 \\ \end{array} } \right]\). Since the AF has no relationship with the angle dimension, there is no rotation in the angle. By rotation, orthogonality is achieved. As a result, this method is also called the error ellipse orthogonality method. Compared with the traversal of parameters in the waveform library [5,6,7, 11, 15,16,17, 23, 25, 26, 29], it had the advantage of low computational burden and intuitively physical meaning. Figure 2a and b shows two situations where the prediction error ellipse and the measurement error ellipse are not orthogonal, whereas Fig. 2c shows that the two error ellipses are orthogonal. The minimal overlapped area is only achieved when the two error ellipses are orthogonal. Meanwhile, the minimum tracking error is achieved [14]. The rotation angle φk+1 can be obtained by Eq. (44).

Fig. 2
figure 2

Different intersection situations of two error ellipses. a The two error ellipses are not orthogonal. b The two error ellipses are not orthogonal. c The two error ellipse are orthogonal

Full size image

Assuming that the covariance matrix of state prediction error in time k is Pk+1|k, the fractional factor φk+1 is:

$$ \varphi_{k + 1} = \arctan [\frac{{\nu_{{\text{P}}} (2)}}{{\nu_{{\text{P}}} (1)}}] + \frac{\pi }{2} - \arctan [\frac{{\nu_{{\text{R}}} (2)}}{{\nu_{{\text{R}}} (1)}}] $$
(44)

where νP(i) and νR(i) are the ith element in the eigenvector corresponding to the maximum eigenvalue of P and Rk, respectively. \({\mathbf{P}} = \left[ {\begin{array}{*{20}c} {P_{11} } & {P_{12} } \\ {P_{12} } & {P_{22} } \\ \end{array} } \right]\), \(P_{11} = \det \left( {\left[ {\begin{array}{*{20}c} {P_{k + 1|k} (x,x)} & {P_{k + 1|k} (x,y)} \\ {P_{k + 1|k} (y,x)} & {P_{k + 1|k} (y,y)} \\ \end{array} } \right]} \right)\), \(P_{12} = \det \left( {\left[ {\begin{array}{*{20}c} {P_{k + 1|k} (x,\dot{x})} & {P_{k + 1|k} (x,\dot{y})} \\ {P_{k + 1|k} (y,\dot{x})} & {P_{k + 1|k} (y,\dot{y})} \\ \end{array} } \right]} \right)\), and \(P_{22} = \det \left( {\left[ {\begin{array}{*{20}c} {P_{k + 1|k} (\dot{x},\dot{x})} & {P_{k + 1|k} (\dot{x},\dot{y})} \\ {P_{k + 1|k} (\dot{y},\dot{x})} & {P_{k + 1|k} (\dot{y},\dot{y})} \\ \end{array} } \right]} \right)\). Substitute φk+1 into Eq. (43), Rk+1 is obtainable. Then, the posterior state error covariance matrix Pk+1|k+1 is also obtainable. Such iteration can be applied to the optimal waveform selection for the next time. The procedure of the MCS-MPDA-SCKF based on the waveform scheduling is shown in Fig. 3.

Fig. 3
figure 3

Procedure for the adaptive MPDA-SCKF

Full size image

5 Simulations

The simulations were run on a single Intel (R) Core (TM) i7-4790CPU (3.6 GHz) processor with 4 GB memory, Windows 7 OS, and MATLAB 2014a. The IMM-PDA-PF [29], the MCS-PDA-SCKF with waveform scheduling (MCS-PDA-SCKF-WWS), and two MCS-MPDA-SCKFs with fixed waveform (MCS-MPDA-SCKF with φ = 0 and MCS-MPDA-SCKF with φ = π/3) were used as the baseline to compare with the proposed algorithm. Two MCS-MPDA-SCKFs with fixed waveform were used to demonstrate the effectiveness of the proposed waveform scheduling, the MCS-PDA-SCKF-WWS was used to show the effectiveness of MPDAF, and the IMM-PDA-PF focused on the entire performance of the proposed algorithm. In the IMM, the target motion models were two second-order white noise acceleration models with different noise levels, whose standard deviations were 0.1 and 10, respectively. In this scenario, the clutter parameters were ρ = 0.004, γ = 16, Pf = 0.01, Pg = 0.997, ηk = (r0/rk)4η0, and η0 = 40 dB. The aircraft maneuvered in a two-dimensional plane, where X0 = [4.32 × 105, -2.2552 × 103, 0.0017, 8.8 × 104, -3.9765 × 102, -9.8097]T and R0 = [10000, 100, 0.01]T. The acceleration changes are shown in Fig. 4, when the aircraft was in a weak maneuver during 0 ~ 42 s, 52 ~ 102 s, and 112 ~ 150 s, and in a strong maneuver during 42 ~ 52 s and 102 ~ 112 s. The base waveform was a rectangular pulse. To assess the tracking quality in clutter, the track loss percentage (TLP), the RMSE, and the mean error were utilized as evaluation metrics. A track loss was considered in case of the estimation falling outside the ten-sigma region centered around the true position for ten continuous tracks [34].

$$ [h(x_{k|k} ) - h(x_{k} )]^{{\text{T}}} R_{k}^{ - 1} [h(x_{k|k} ) - h(x_{k} )] > 100 $$
(45)
$$ E_{{{\text{RMSE}},k}} = \left\{ {\frac{1}{M}\sum\limits_{j = 1}^{M} {\left\| {{\mathbf{X}}_{i,k} - {\hat{\mathbf{X}}}_{i,k}^{j} } \right\|_{2}^{2} } } \right\}^{1/2} $$
(46)
$$ E_{{{\text{ME}}}} = \frac{1}{N}\sum\limits_{k = 1}^{N} {E_{{{\text{RMSE}},k}} } $$
(47)

where Xi,k and \({\hat{\mathbf{X}}}_{i,k}^{j}\) are the true values and the estimation values of ith component of the aircraft state vector in time kT in the jth simulation, respectively. M is the total number of simulations, and N is the total number of steps in a simulation. The statistical results were obtained through 100 Monte Carlo simulations, and they are shown in Figs. 5, 6 and Table 1.

Fig. 4
figure 4

Aircraft acceleration

Full size image
Fig. 5
figure 5

a RMSE on position. b RMSE on velocity. c RMSE on acceleration

Full size image
Fig. 6
figure 6

Rotation angle along with the simulation

Full size image
Table 1 Mean error, track loss percentage, and runtime
Full size table

5.1 Simulation Results and Analysis

From Fig. 5a–c, the proposed algorithm achieved the best performance among all the algorithms whenever the aircraft was in the strong maneuver and weak maneuver. In particular, when faced with abrupt state change, the proposed algorithm achieved the lowest RMSEs and the fastest convergence rate. Figure 6 shows the rotation angle of the waveform via the FRFT in the proposed algorithm. It can be seen that the rotation angle was dynamically configured to select the optimal waveform. As a result, the minimum tracking error was achieved from time to time, and the performance was improved. By contrast, the two filters with fixed waveform failed to exploit the waveform scheduling, leading to a less satisfactory performance than the proposed algorithm. Additionally, it should be noted that the simulation was in a situation of high SNR. Thus, the increase in distance had little effect on the SNR and thus on RMSEs. The increase in RMSE along with the distance, shown in [17] [29], was not apparent.

Table 1 shows the mean error track loss percentage and runtime of all the algorithms. It can be seen again that the proposed algorithm achieved the lowest mean error as well as the lowest track loss percentage. Additionally, the track loss percentage of IMM-PDA-PF, MCS-PDA-SCKF-WWS, and the proposed algorithm was close. By contrast, the track loss percentage of two filters with fixed waveform was twice as large as the proposed algorithm. It appears that the waveform scheduling contributed to the tracking performance improvement in clutter. However, it should be noted again that the scenario was in a high SNR environment. As a result, the track loss percentage of all the algorithms was low. Additionally, it can be seen that the IMM-PDA-PF provided not only poor mean error but also the longest runtime, which is more than three times as large as the proposed algorithm. It was, therefore, impractical in this application. By contrast, the proposed algorithm achieved the lowest mean error while maintaining an efficient operation. By utilizing two filters with fixed waveform, time was marginally extended owing to the additional waveform scheduling module. However, the performance was improved significantly. In particular, compared with the MCS-MPDA-SCKF with φ = 0, the proposed algorithm increased the estimation precision on position by 32.46%, the estimation precision on velocity by 25.62%, and the estimation precision on acceleration by 10.37%. The runtime was increased by 13.23%. Compared with the MCS-MPDA-SCKF with φ = π/3, the proposed algorithm increased the estimation precision on position by 34.83%, the estimation precision on velocity by 35.73%, and the estimation precision on acceleration by 24.17%. At the same time, the runtime was increased by 12.98%. As a result, it achieved higher tracking precision with only a small additional increase in computational complexity.

5.2 Computational Complexity Analysis

The runtime comparison is shown in Sect. 5.1. The qualitative computational complexity is analyzed here. However, the five filters were with different motion models as well as different nonlinear filters. As a result, the MCS-MPDA-SCKF with φ = 0 was chosen as the baseline. Since the two filters were both with fixed waveform, the MCS-MPDA-SCKF with φ = π/3 had approximate computational complexity compared with the baseline. The waveform selection module was added to the proposed algorithm; as a result, its computational complexity is marginally higher than the baseline. The main difference between the MCS-PDA-SCKF-WWS and the proposed algorithm was the MPDAF. However, the computational burden of MPDAF and PDAF was similar. As a result, the MCS-PDA-SCKF-WWS provided a similar runtime compared with the proposed algorithm. As for the IMM-PDA-PF, two filters were needed to work in parallel, and the interacting input and output were additional, and the PF provided more than three times the computational burden compared with the SCKF. Additionally, the waveform parameters required traversing in the library, which resulted in an increased burden compared with the proposed waveform scheduling method. Thus, the computational complexity of the IMM-PDA-PF was about three times as large as the proposed algorithm.

6 Conclusion

This paper proposed an efficient algorithm to solve the problem of target tracking in nonlinear measurement and clutter environment. Based on the MCS motion model, the MPDAF and the SCKF were integrated as the tracker. Simultaneously, an efficient waveform selection algorithm was proposed to configure the waveform dynamically. The simulation results showed that the proposed algorithm achieved lower RMSE, mean error, and track loss percentage than the algorithms with fixed waveforms while maintaining a reasonable runtime improvement. Additionally, compared with the two state-of-the-art algorithms, the proposed algorithm possessed a simpler structure and higher estimation precision. However, more research is required to be done such as extending the proposed model to multi-step ahead scheduling or a low SNR environment.