Fixed lag smoothing target tracking in clutter for a high pulse repetition frequency radar

A new method to smooth the target hybrid state with Gaussian mixture measurement likelihood-integrated track splitting (GMM-ITS) in the presence of clutter for a high pulse repetition frequency (HPRF) radar is proposed. This method smooths the target state at fixed lag N and considers all feasible multi-scan target existence sequences in the temporal window of scans in order to smooth the target hybrid state. The smoothing window can be of any length N. The proposed method to smooth the target hybrid state at fixed lag is also applied to the enhanced multiple model (EMM) tracking algorithm. Simulation results indicate that the performance of fixed lag smoothing GMM-ITS significantly improves false track discrimination and root mean square errors (RMSEs).


Introduction
It is well known that, when a Pulse-Doppler radar operates in high pulse repetition frequency (HPRF) mode, the target range information is ambiguous due to the aliasing range [1]. In the absence of measurement noise, multiple possible target ranges project onto the same range measurement.
Nagel and Hommel [2] propose the range-gated HPRF method to resolve range ambiguity; however, because of the limited number of bursts, the number of range gates is often not sufficient, especially in a cluttered environment with uncertain target detection. The authors in [3] present a multiple model algorithm to eliminate the range ambiguity problem in an HPRF radar. These references estimate the trajectory state without providing a track quality measure for false track discrimination (FTD). The Gaussian mixture measurement likelihoodintegrated track splitting algorithm (GMM-ITS) [4] and the enhanced multiple model algorithm (EMM) (it incorporates the track quality measure in a multiple model algorithm (MM) [3]) are investigated in [5] for singletarget tracking in clutter using an HPRF radar; both algorithms are capable of trajectory estimation and false track discrimination.
The application of smoothing is quite effective in the situation awareness and threat assessment applications. The past state of the target is updated using all measurement information received till the current scan. Smoothing produces a reduction in estimation error and improves the FTD. The fundamental techniques of smoothing in estimation are provided in [6].
The idea of fixed lag smoothing is proposed in [7]. Later, its application in target tracking is considered in different ways. [8] consider the multi-scan data association for tracking a target and applies fixed lag smoothing; however, it does not consider the track quality measure. [9] applies the RTS to the MHT algorithm, but it also ignores the track quality measure. Chakravorty and Challa [10] propose augmented state integrated probabilistic data association (ASIPDA). In ASIPDA, however, the smoothing probability of target existence uses the smoothing innovation obtained only in the current scan. In sIPDA [11], the authors use the Fraser-Potter [12] approach for smoothing with IPDA. Another extension of smoothing considers track splitting for single-target smoothing [13].
This paper presents a new method to smooth the target hybrid state at fixed lag N. It applies the proposed smoothing algorithm on both GMM-ITS and EMM algorithms to obtain the smoothing benefits for an HPRF radar. In this technique, a relatively small window of scans is defined in each smoothing interval. The target trajectory state and target existence state are smoothed at fixed lag N in the smoothing interval. This technique considers all feasible multi-scan target existence events and smoothed state estimates (using the augmented state GMM-ITS update) calculated at all intermediate scans in the smoothing interval in order to smooth the target hybrid state at fixed lag N.
The fixed lag smoothing-based tracking algorithms proposed in this work also provide a significant improvement over existing online algorithms in terms of false track discrimination and root mean square errors (RMSEs). This paper is organized as follows. The basic models used are discussed in Section 2, while in Section 3 an overview of GMM-ITS and EMM algorithms is provided. The GMM-ITS algorithm for the augmented target state is extended in Section 4. The fixed lag smoothing GMM-ITS (FLs GMM-ITS) algorithm is proposed next in Section 5. The simulation results are presented in Section 7, followed by concluding remarks in Section 8.

Mathematical model
This section introduces the models used for target motion and sensor. The fundamental nomenclature used in this paper is provided in Table 1.

Target model
The existence of the target in any scan k is a random event and is denoted by χ k . The target propagates following the Markov Chain One model [14]. The probability that the target exists between scan k − 1 and scan k is given by where T k−1,k denotes the time interval between scan k − 1 and scan k and T ave T k−1,k is the average duration of target existence [15]. The target existence χ k and the target non-existenceχ k are mutually exclusive and exhaustive events. If the target does not exist in any scan, then the probability that it will continue its state of non-existence in the subsequent scans is The trajectory of the target is modeled as where v k−1 is the plant noise with zero mean and known covariance Q k−1 . F k−1 denotes the state propagation matrix and is assumed to be known. These matrices are as follows: and where the scalar q denotes the root mean square of the acceleration plant noise and T is the sampling time.

Measurement model
At each scan k, the sensor returns a set of measurements Z k . Let Z k,i denote the i-th measurement in Z k . The origin of each measurement is unknown (target detection or clutter). Each measurement Z k,i at time k consists of azimuth θ k,i , range r k,i , and Doppler velocity d k,i . After decorrelation between the range and Doppler of each measurement [16], every measurement becomes Z k,i

Target measurement
Due to the ambiguous range, GMM-ITS regenerates a set of Gaussian measurement components to update the track. Measurement components G k,i with respect to each  [5]. The measurement is the function of trajectory state and sensor noise and is equal to where ε θ k , ε r k , and ε w k , respectively, denote the zeromean white Gaussian measurement noise with respect to azimuth, range, and decorrelated Doppler measurement. h θ k , h r k , and h w k are the target azimuth, range, and decorrelated Doppler measurement functions, respectively, along with the range and decorrelated Doppler Jacobians H r (x k ) and H w (x k ) [16] where x denotes the norm of x. x p k denotes the position vector relative to sensor. with where x v k denotes the velocity vector relative to the sensor,

Clutter measurement
At each scan k, the number of clutter measurements follows the Poisson distribution. The intensity of the Poisson process at each surveillance space point Z k,i is termed as clutter measurement density and is denoted as where ρ p denotes the clutter measurement position density and ρ d is the Doppler component density of clutter measurement [5]. The polar measurement is transformed into Cartesian position measurement as [17] where y g,p k,i is the position component of measurement y

An overview on GMM-ITS algorithm and EMM algorithm
This section presents a brief description about the GMM-ITS algorithm and the EMM tracking algorithm.

Gaussian mixture measurement ITS
The GMM-ITS algorithm is first introduced in [4], later it is applied to the problem of single-target tracking in clutter using an HPRF radar [5,18]. In GMM-ITS algorithm, the non-linear (non-Gaussian) measurement likelihood is approximated by a Gaussian mixture of measurement components, which corresponds to the ambiguous measurement components due to the aliasing target range in the application of an HPRF radar.

GMM model for an HPRF radar
As presented in Section 2.2.1, each measurement Z k,i received from the HPRF radar regenerates G k,i Gaussian measurement components defined by mean y g k,i , covariance R g k,i , and relative weight γ g k,i . Thus, the conditional likelihood of the measurement Z k,i is approximated by a Gaussian mixture of G k,i measurement components. where

GMM model integration with ITS algorithm
Once the GMM model for HPRF radar is formed, the GMM-ITS algorithm absorbs the integrated track splitting (ITS) to proceed subsequent tracking procedure. The ITS [15] algorithm is a multi-scan tracking algorithm, which updates the target trajectory state and the target existence state at each scan k using multi-scan data association. The GMM-ITS algorithm uses ITS for data association to obtain the improvement in the performance of tracker in clutter using an HPRF radar [5].

Extended multiple model (EMM)
The multiple model (MM) tracking algorithm for an HPRF radar is proposed in [3], where each model proceeds in the probabilistic data association (PDA) [19] sense independently. The MM tracking algorithm is modified by incorporating the probability of target existence update to develop enhanced multiple model (EMM) tracking algorithm [5].

Augmented state GMM-ITS
This section presents one complete recursion cycle for augmented state GMM-ITS (AS GMM-ITS). At each scan, the augmented state is smoothed, and the results are later used in Section 5.1 to determine the smoothed target state at a fixed lag of N using the proposed method.
The GMM-ITS algorithm [4] details are omitted, and only its application to the augmented state is discussed in order to minimize the complexity at this stage. Let k τ = k − N + r, (0 ≤ r ≤ N) be the variable to address each scan in the smoothing interval. Each measurement in each scan not selected by the already established track initializes a new track [5].
At scan k − N, the track state is initialized as

State augmentation
At scan k τ = k − N (r = 0) in the smoothing interval, the track trajectory state and its associated error covariance matrixX AS k−N|k−N and P AS k−N|k−N are augmented, respectively, aŝ and The superscript AS stands for augmented state, where the (−) sign on the right hand side of (20) corresponds to cross covariance terms of state elements of the augmented state, not detailed here for reasons of clarity. The augmented state propagation matrix is The augmented plant noise covariance matrix is The augmented measurement matrix with respect to the position component of measurement becomes where The linearized augmented measurement coefficient matrix with respect to the Doppler component of measurement becomes where I n is an n-dimensional identity matrix, and 0 n,n and 0 m,n are matrices of zeros, where n is the order of the target state vector, and m is the order of the position measurement vector. Terms w k−1 , v k , F k−1 , and Q k−1 are defined in Section 2.1. The order of matrices F AS and Q AS is ((N + 1).n × (N + 1).n), while H AS,p has dimensions of (m × (N + 1).n), and H AS,w has dimensions of (1 × (N + 1).n).

State prediction
At any scan k τ = k − N + r, (1 ≤ r ≤ N) inside the smoothing interval, the predicted augmented state conditioned on the component c k τ −1 is obtained by standard Kalman filter where KF P represents the standard Kalman filter prediction. The probability density function (PDF) for the predicted augmented state and error covariance at scan k τ is

Measurement selection and likelihood calculation
To reduce the computational requirements, a subset of validated measurements from all the measurements received by the sensor is selected at each scan k τ in the smoothing interval. A gating procedure [20] is used to select the validated measurements for each measurement component g corresponding to each track component c k τ −1 . A gating test (28) is applied to the position component (y g,p k τ ,i ) of each measurement received at each scan in the smoothing interval.
To simplify the notations in the rest of the paper, denote and where κ is the selection threshold. In simulated twodimensional surveillance, κ is selected as 13.3, which corresponds to the gating probability P g = 0.99. Each i-th validated measurement is represented by z k τ ,i and has G s k τ ,i validated components. Hereafter, the parameters y i are attributed to the selected measurement z k τ ,i . At any scan k τ in the smoothing interval, the measurement likelihood of the selected measurement z k τ ,i becomes where p g k τ ,i is the likelihood of measurement component y where p g,c kτ −1 k τ ,i is the likelihood of measurement component y g k τ ,i with respect to track component c k τ −1 where S s p k τ ,i is as defined in (30). The augmented track component is updated by the measurement position component using standard Kalman filter update represented by KF U The likelihood of the decorrelated Doppler component conditioned on the position component of measurement y g k τ ,i is calculated in (36). A transformation is applied to augmented stateX s p k τ |k τ ,i using transformation matrix T ι to provide only the filtered statex g,c kτ −1 ,p k τ |k τ ,i (the predicted Doppler mean at current scan k requires only the filtered state in (8)(9)(10)(11)(12)(13)) for the calculation of Doppler likelihood of the i-th measurement at current scan k. To simplify the notation in the rest of the paper, denote s w = AS, g, c k τ −1 , w . wherê and where

State update
The decorrelated Doppler measurement component w g k τ ,i is used to update the track using standard extended Kalman filter EKF U where The fixed lag smoothed trajectory state in Section 5 requires the Gaussian sum of all the track components to obtain the augmented track state at scan k τ in the smoothing window using (44) and (45).
The component probabilityc k τ is updated as The likelihood ratio at scan k τ is defined as The track augmented trajectory stateX AS k τ and its augmented error covariance matrix P AS k τ provide filtered and smoothed trajectory states ˆx k τ |k τ ,ˆx k τ −1|k τ , . . . ,ˆx k τ − N|k τ } and their corresponding error covariance matrices P k τ |k τ , P k τ −1|k τ , . . . , P k τ −N|k τ at scan k τ in the smoothing interval.
The target existence state at scan k τ updates as [15] P

Fixed lag smoothing GMM-ITS
The fixed lag smoothing GMM-ITS (FLs GMM-ITS) algorithm for any size of smoothing lag N is proposed.
The fixed lag smoothing IPDA (FLs-IPDA) [21] uses a single-scan data association in the smoothing interval. It considers only two scans in the smoothing interval and does not provide any generalization to smooth the target state at any size of fixed lag.
In the FLs GMM-ITS algorithm, the concept of fixed lag smoothing is extended for a more general case of smoothing for any size of lag N, and it utilizes the benefits of multi-scan data association. The results of the augmented state GMM-ITS algorithm derived in Section 4 are used in the smoothing interval. The smoothed state with respect to scan k − N, obtained in (44) and (45), at each scan is weighted by the probabilistic weights calculated using multi-scan target existence events. In Sections 5.1 and 5.2, the FLs GMM-ITS algorithm updates the target trajectory state p x k−N |χ k−N , Z k and target existence state P χ k−N |Z k in the smoothing interval.
In the next smoothing interval, the target state with respect to scan k−N is ignored (as it was updated in the last smoothing interval), and a future scan k + 1 is added in the smoothing interval to smooth the track state at scan k − N + 1.
In the smoothing interval, there are N +1 feasible multi-scan target existence events.
The conditions for feasible multi-scan target existence events are: (1) Target exists at all scans k τ in the smoothing interval, where 1 ≤ r ≤ N.
(2) Non-existence of target at any scan k τ implies target non-existence at all following scans.
(3) Target does not exist at any scan k τ in the smoothing interval.
These conditions considers N + 1 feasible multi-scan target existence events in the smoothing interval. The number of feasible multi-scan target existence events increases in a linear manner with the increase in smoothing lag. The next sections present the formulas for smoothed hybrid target state for any fixed lag N.

Smoothed target trajectory state at fixed lag N
The target trajectory state is composed of position and velocity components. The information of N future scans are used to smooth the target trajectory state at scan k − N. The feasible multi-scan target existence events are used to weight the smoothed target state estimates received at each scan in the interval. The smoothed state estimates are obtained by augmented trajectory state update (44) and augmented state error covariance matrix update (45) at each scan k τ in the interval.
At the last scan k (r = N) in any smoothing interval, the fixed lag smoothed trajectory state estimate is weighted Gaussian sum of smoothed trajectory state estimates (from Section 4) obtained at each scan k τ . These weights are calculated using feasible multi-scan target existence events calculated based on the above conditions. The smoothed target trajectory state estimate at fixed lag k − N is In general, (49) becomes (50) The first term in the summation on the right-hand side is the smoothed trajectory state estimateˆx k−N|k τ at each scan in the smoothing interval k τ =[ k−N, k−N +1, . . . , k] in Section 4 using (44) and (45). ϒ(s) are the probabilistic weights calculated using the multi-scan target existence events in the smoothing interval.
To maintain the clarity, the index variable s is used to address each scan in the interval, 1 ≤ s ≤ N + 1. The factor ϒ(s) is The weighting factor ϒ(s) satisfies N+1 s=1 ϒ (s) = 1.
The term in (51) for s = 1 is The term in the numerator becomes The first joint likelihood term on the right hand side implies that all measurements belong to the clutter; thus, it does not correct the track estimate. The second joint likelihood term contributes in the track update and becomes and for s = N + 1, Each likelihood term in (54), (57), and (58) is calculated as where p ρ,k−s+1 is the likelihood that a measurement belongs to clutter, and μ F m k−s+1 is the Poisson distribution function for the clutter measurements. k−s+1 is calculated using (47) by replacing k τ with k − s + 1.
The denominator term p Z k , Z k−1 . . . Z k−N+1 |χ k−N , Z k−N is the normalizing factor. The weighting factor ϒ(s) in (50) for any s in the smoothing interval is defined as The factors Num( ) and Den( ) are defined as where and p 11 is the state transition probability.

Smoothed target existence state at fixed lag N
All feasible N +1 multi-scan target existence events (Section 5) in the smoothing interval are calculated to determine the smoothed target existence state at fixed lag N. The target existence at fixed lag N is the result of total probability theorem and assumes all above events in the multi-scan event space (64) Compared to (64), in (51), the multi-scan target existence probabilities are conditioned on the target existence event χ k−N , as target trajectory state for non-existence event does not mean anything.
At the start of recursion at scan k − N, the known target existence probability is P{χ k−N |Z k−N }. Equation (64) is expanded following the similar procedure that is used to obtain the ϒ(s) in (51), to propose a procedure to obtain the smoothed target existence state at any scan k τ = k − N + r in the smoothing interval for 0 ≤ r ≤ N: where Num(s) = and where k−n is calculated using (47) by replacing k τ with k − n. In the next smoothing interval, the recursion starts at scan k − N + 1 with target existence probability P{χ k−N+1|k−N+1 } calculated by (48) by replacing k τ with k − N + 1. The target existence state is then smoothed using (65) and (66).
The smoothed trajectory state and error covariance matrix [ˆx k−N|k−s+1 , P k−N|k−s+1 ] for 1 ≤ s ≤ N + 1 in the smoothing interval are obtained in Section 4 using augmented state smoothing algorithm. At the final scan k in the smoothing interval, all these smoothed trajectory state estimates are weighted by the smoothing weights ϒ(s) to obtain the smoothed target trajectory state at fixed lag N using (50). The smoothed target existence state at fixed lag N is also calculated (64), which is an outcome of total probability theorem.

Fixed lag smoothing enhanced MM HPRF tracker
The fixed lag smoothing approach presented above is also applied to extended multiple models (EMM) HPRF tracker. The EMM tracker incorporated the track quality measure in multiple model HPRF tracker [3]. The augmented state is also calculated for the EMM algorithm, and then fixed lag smoothing is used to obtain the smoothed target hybrid state at any scan. The fixed lag smoothed enhanced multiple model HPRF tracker (FLs EMM) is also simulated in the proposed work to compare its performance.

Simulation study
This simulation sections compares the tracking performance of smoothing algorithms (FLs GMM-ITS and FLs EMM) with online tracking algorithms (GMM-ITS and EMM). The results provide the information on the relative benefits of the proposed fixed lag smoothing algorithm. A two-dimensional single-target tracking scenario is presented in Fig. 1. The target can only be detected in the area bounded by the minimum and maximum range and an azimuth between 0 and π/2. The target is moving with a uniform velocity of v T = 40 m/s and with a heading angle of 0 • . The initial position of the target is described by the polar coordinates (15, 000 , 1.0472 rad); the maximum possible velocity for the target is assumed to be v max = 400 m/s, which is used for the one-point track initialization. The HPRF radar is stationery in the origin, which returns a set of measurements (which consists of azimuth, range, and Doppler) at each scan; the minimum and maximum of the detection target range are R max = 20 km and R min = 0 km, respectively. The radar measurement is corrupted by white Gaussian noise with zero mean and standard deviation σ θ = 0.1 • , σ r = 20 m, and σ d = 1 m/s, respectively. The range measurement is correlated with Doppler measurement by coefficient η = −0.8. The radar provides the target position measurement with detection probability P d = 0.8. The uniform clutter position measurement density is ρ p = 10 −7 m −2 (31 clutter measurements per scan on average). The clutter Doppler measurement density is assumed to have a Gaussian probability density function with zero mean and standard deviation of σ c d = 10 m/s. The staggered PRI [22] is applied in the simulations, and the maximum unambiguous range of the HPRF radar is R u (k) = 4005 m, which results in five measurement components for each radar returned measurement on average. Each simulation experiment consists of 1000 runs, and each run simulates 40 scans with sampling time T = 2 m/s. In Fig. 2, a second target is added in the surveillance area crossing the first target, which depicts the multi-target tracking scenario. The second target is moving with a uniform velocity of v T = 46 m/s and with a heading angle of 330 • . The initial position of the second target is (15, 800 , 1.0764 rad). The overlapping region of both the targets is between scans 17 and 28, where both targets share the measurements dominantly in this region.
The Markov Chain One model [14] of target existence is assumed with transition probabilities where p ij denotes the transition probability from event i to event j. Event 1 and event 2 represents χ k andχ k , respectively. The confirmation thresholds for the algorithms are tuned to deliver the same number of confirmed false tracks (CFT) (≈ 7) for all 1000 Monte Carlo simulation runs. In both simulation scenarios, the CFT are set to be same for a fair comparison. Both smoothing algorithms consider the lag size of 2. The confirmed true track rate for a single target in Fig. 3 presents a comparison between FLs GMM-ITS and other algorithms simulated in this work for single-target tracking. FLs GMM-ITS performs better than both online (GMM-ITS, EMM) algorithms and offline (FLs EMM) tracking algorithm under the given surveillance environment.
The RMSEs of confirmed true tracks for single-target tracking scenario are presented in Fig. 4. The smoothing algorithms perform better than online algorithms. At the last scan k = 40, no future information exists, so all algorithms have an identical response in terms of RMSEs.
FLs GMM-ITS and FLs EMM algorithms looks to have similar RMSEs, although one expects better performance for the FLs GMM-ITS algorithm. The reason is that the true tracks which are confirmed earlier in the FLs GMM-ITS compared to FLs EMM generate less favorable performance in terms of estimation errors.
In Fig. 5, the cumulative confirmed true track rate for the multi-target simulation scenario is presented. From scan 1 to scan 20, the confirmed true track rate reaches to almost 100 %, and the FLs GMM-ITS algorithm shows the fastest response. From scan 20, the online algorithms loose almost 50 % of the targets; in this period, the FLs GMM-ITS algorithm performs better compared to the other algorithms and maintains a high confirmed true track rate. The drop in confirmed true track is expected; when the targets share measurements, the singletarget tracking algorithms are expected to loose tracks at a higher rate. The FLs GMM-ITS algorithm recover much faster as compared to other algorithms from scan 28 till the last scan, which indicates the high efficiency of the proposed FLs GMM-ITS algorithm. Figures 6 and 7 present the RMSEs for target 1 and target 2, respectively, for the multi-target simulation scenario depicted in Fig. 2. It is evident that the smoothing algorithms perform significantly better as compared to the online tracking algorithms. Between scan 17 and scan 28, where the true targets share measurements more prominently, the estimation errors are increased because of association of confirmed true track with measurements from other targets. The smoothing algorithms perform much better as compared to the online algorithms and have reduced RMSEs in this interval. The smoothing and filtering algorithms have identical results at the final scan as no future information is available to smooth the target state.
The computational time for the EMM algorithm is set as a reference to determine the percentage of extra computational time needed by other algorithms. The total sampling time is 80, 000 s, which is greater than the computation time of the aforementioned tracking algorithms; therefore, all algorithms are capable of working in real time.
It is evident from the simulation results that use of the smoothing algorithm improves the performance of the tracker in terms of both RMSEs and false track discrimination, at the cost of some delay in both single-target and multi- target tracking simulation environments. Table 2 provides the comparison of computational time for different tracking algorithms.

Conclusions
This paper provides a new procedure to calculate the smoothed target hybrid state at fixed lag N. The benefits of the smoothing algorithm are compared for target tracking algorithms in the clutter using the HPRF radar. The smoothed augmented states (obtained at each scan in the smoothing interval) and their respective weights (calculated using multi-scan target existence events) are used to obtain the fixed lag smoothed trajectory state estimate for the target. The target existence state is also smoothed at fixed lag N using all feasible multi-scan target existence events in the smoothing interval.  The FLs GMM-ITS performs much better than FLs EMM, GMM-ITS, and EMM algorithms to track the target using an HPRF radar in terms of RMSEs and false track discrimination. The smoothed target hybrid state provides small estimation errors and produces excellent false track discrimination in the simulation conditions used in the proposed work.