Ground Moving Target Tracking Filter Considering Terrain and Kinematics

This paper addresses ground target tracking (GTT) for airborne radar. Digital terrain elevation data (DTED) are widely used for GTT as prior information under the premise that ground targets are constrained on terrain. Existing works fuse DTED to a tracking filter in a way that adopts only the assumption that the position of the target is constrained on the terrain. However, by kinematics, it is natural that the velocity of the moving ground target is constrained as well. Furthermore, DTED provides neither continuous nor accurate measurement of terrain elevation. To overcome such limitations, we propose a novel soft terrain constraint and a constraint-aided particle filter. To resolve the difficulties in applying the DTED to the GTT, first, we reconstruct the ground-truth terrain elevation using a Gaussian process and treat DTED as a noisy observation of it. Then, terrain constraint is formulated as joint soft constraints of position and velocity. Finally, we derive a Soft Terrain Constrained Particle Filter (STC-PF) that propagates particles while approximately satisfying the terrain constraint in the prediction step. In the numerical simulations, STC-PF outperforms the Smoothly Constrained Kalman Filter (SCKF) in terms of tracking performance because SCKF can only incorporate hard constraints.


Introduction
Ground tracking radars mounted on airborne platforms play a key role in many applications, especially those for military purposes; surveillance, airstrike, and escort missions done by aircraft commonly require precise tracking of ground targets. In several modern military campaigns, ground moving target indicator (GMTI) radar on-board the Joint Surveillance Target Attack Radar System (STARS) has been proven strategically and tactically significant [1]. Accordingly, algorithms that track ground targets running on radars are becoming more important. Although there have been great advances in target tracking, tracking ground targets is still a challenging problem. The reason is that the characteristics of ground target tracking are different from those of tracking other types of targets. (e.g., high clutter, terrain obscuration, etc.) [2].
Because exploiting appropriate assumptions other than the state-space model can help to improve the statistical inferences of the system [3], many studies have tried to introduce useful assumptions to ground target tracking. They can be classified based on two criteria of what or how assumptions were applied.
Based on the first criteria, existing studies can be classified into two further categories. The first category considers the behavior characteristics of the ground target that are distinguished from those of airborne targets. Fosbury [4] and Kastella [5] each created the terms 'trafficability' and 'hospitability of maneuver', which represent how easily a

•
We propose a sample-efficient particle filter to which the terrain constraint can be applied. The proposed algorithm is named Soft Terrain Constrained Particle Filter (STC-PF). Given the assumption of target motion, STC-PF performs sampling in a direction for which the state satisfies the constraint during the propagation step. As a result, STC-PF is more sample-efficient than scPF. Furthermore, in the numerical simulations, STC-PF using soft terrain constraint outperforms Smoothly Constrained Kalman Filter (SCKF) [30] using hard constraint in terms of tracking performance. • Using a Gaussian process, terrain constraint is formulated as a soft position constraint along with a soft velocity constraint. Because kinematics states that position and velocity is not independent, a constraint on the position of a target implies that the velocity of the target will be constrained as well. Therefore, terrain constraint includes both position constraint and velocity constraint. Furthermore, terrain constraint requires exact terrain elevation and its gradient at an arbitrary position, but DTED (Digital Terrain Elevation Data) [36] cannot provide them. To overcome this issue, we model the ground-truth terrain elevation with a Gaussian process (GP) and treat DTED as a noisy observation [37] of it.
Technically, we used SRTM (Shuttle Radar Topography Mission). However, we will use the term DTED and SRTM interchangeably as they both are data that map terrain elevation of the entire globe.
The structure of this paper is as follows: In Section 2, tracking of a ground target with a terrain constraint is formulated. Section 3 presents the proposed algorithm, STC-PF. Section 4 provides detailed explanations, the results, and a discussion of the numerical simulation. Finally, in Section 5, we conclude.

Problem Formulation
In this section, tracking of a ground target with terrain constraint is formulated as a constrained state estimation problem.
Consider a system described by the following state-space model: where x k is the system state vector at time k, y k the measurement vector, f the system function, g the observation function, w k the process noise vector, and n k the measurement noise vector. The system state vector x k ∈ R 6 consists of the position (x k , y k , z k ) and the velocity (v x,k , v y,k , v z,k ) in local Cartesian coordinates at time k. The system function is a possibly nonlinear function but is assumed to be a constant velocity model in this paper. y k ∈ R 3 is the measurement, which consists of range, azimuth angle, and elevation angle measured from the radar. w k ∼ N(0, Q) is white Gaussian process noise, and n k ∼ N(0, R) is white Gaussian measurement noise. Subsequently, Equations (1) and (2) are realized as follows: The final goal of the state estimation problem is to infer the state sequence of the dynamical system x 0:k from the series of observations y 1:k . Now, the terrain constraint can come into play to incorporate the additional information that the state-space model cannot reflect. The terrain constraint not only represents the assumption that the position of a ground target should be located on the terrain surface but also that the velocity vector of the target should be tangent to the terrain surface. Both assumptions can be transformed into state constraints as follows: where λ k , ϕ k , and h k are the latitude, longitude, and altitude (LLA) of the target at time k. h(λ, ϕ) is ground-truth terrain elevation at latitude λ and longitude ϕ. Note that we do not have direct access toh, but only noisy observations, such that DTED(λ, ϕ) =h(λ, ϕ) + (λ, ϕ).

Soft Terrain Constrained Particle Filter
In this section, the newly proposed algorithm, Soft Terrain Constrained Particle Filter (STC-PF) is derived. In Section 3.1, mathematical modeling of ground-truth terrain elevation is presented. Then, we propose a technique for the transformation of velocity between the LLA coordinates and the local Cartesian coordinates in Section 3.2. Necessary assumptions required for algorithm derivation are described in Sections 3.3. After the algorithm derivation in Section 3.4, we show the similarity between STC-PF and scPF [35] in Section 3.5.

Modeling of Ground-Truth Terrain Elevation
Although the terrain constraint (Equation (5)) requires the ground-truth elevation, it is almost impossible in practice to retrieve it at an arbitrary position. The reason is that DTED provides neither accurate ground-truth terrain elevation (Equation (7)) nor terrain elevation at arbitrary positions. (Equation (6)) This challenge can be met by reconstructing the ground-truth terrain elevation with a Gaussian process (GP) and treating the DTED as independent observations:h ∼ GP(m(λ, ϕ), k((λ, ϕ), (λ , ϕ ))) DTED(λ, ϕ) =h(λ, ϕ) + ∼ N(0, σ DTED ) (8) where the observation noise σ DTED can be estimated from the work of Rodriguez [37]. (see Appendix B) Because GP assigns a probability for each possible terrain, the terrain constraint becomes stochastic. An advantage of this approach is that it enables us to compute the gradient ofh analytically, which is required to apply the velocity constraint. (Equation (5)) More strictly, joint predictive distribution for ground-truth terrain elevation and its gradient can be expressed in a closed-form, (detailed description is in Appendix A) h, ∇h | DTED ∼ N(μ,Σ) (9) provided that the kernel function is differentiable. Figure 1 shows an example of prediction results when zero mean function and squared exponential kernel are utilized.

Velocity Transformation
Another major challenge when applying the terrain constraint to the filter is that the conversion of velocity between the local Cartesian coordinates and the LLA coordinates is not straightforward. More specifically, the terrain constraint (Equation (5)) requires the velocity in LLA coordinates.
This challenge can be met by multiplying the Jacobian, which is obtained by numerical differentiation. Additionally, because velocity in local Cartesian coordinates is relative while that in LLA coordinates is absolute, the velocity of the radar V lla,ownship should be added (or subtracted) after (or before) multiplying by the Jacobian.
Conversion from LLA to local Cartesian can be done in a converse way.

Assumptions
Regarding the motion of the target, we assume the followings: 1.
The vertical position (h) can be determined provided that the horizontal position (λ, ϕ) is fixed.

2.
Then, the vertical velocity (v h ) can be also determined when the horizontal velocity (v λ , v ϕ ) is fixed.
In Figure 2, assumptions 1 and 2 correspond to the red arrows that inbound to h and v h , respectively. They comprise the 'elevation model'. Due to the recursive Markovian structure, it is possible to infer the current latent state from the previously inferred latent state and the current measurement. Mathematically, by Bayes' rule, the joint distribution of x 0:k given y 1:k can be expressed as The dynamic model P(λ k , ϕ k , v λ,k , v ϕ,k |x k−1 ) and the likelihood model P(y k |λ k , ϕ k , h k ) are found in the above equation. Respectively, they correspond to the blue arrows and the green arrows in Figure 2. Note that the measurement y k is only affected by the position of the target (λ k , ϕ k , h k ), as stated in Equation (4).

Algorithm
The proposed algorithm is based on the SIR(Sequential Importance Resampling) particle filter. In the SIR algorithm, which forms the basis of most sequential Monte Carlo (MC) filters [38], the posterior probability density function P(x 0:k |y 1:k ) is characterized by the set of support points {x i 0:k } N p i=1 (or particles) and the corresponding weights , where N p is the number of particles [39]. The posterior density at time k is approximated as such that where δ(·) represents the Dirac delta function. We assume that the particles are sampled from a well-known proposal distribution, Then, by the principle of importance sampling, the corresponding weight is calculated as Because we have freedom to choose the proposal distribution, we consider a proposal distribution that has a form of q(x 0:k |y 1:k ) = q(x k |x 0:k−1 , y 1:k )q(x 0:k−1 |y 1:k−1 ). (18) In other words, one can draw new support points x i 0:k ∼ q(x 0:k |y 1:k ) by augmenting each of the previous support points Starting from Equation (17), Together with Equation (18) and the recursive relation (Equation (13)), the weight update equation can be simplified.
A further assumption regarding the proposal distribution, This means that the weight of each particle is updated proportionally to its corresponding likelihood. Note that the above weight update equation implicitly includes the normalization given by Equation (15).
The proposed algorithm, STC-PF, is summarized in Algorithm 1. In a vanilla SIR PF, the next state is propagated through the dynamic model only. In contrast, in STC-PF, the next state is propagated through the dynamic model first and then propagated through the elevation model (line numbers 5 and 6 in Algorithm 1).
In Figure 3, a detailed implementation of the elevation model propagation is shown. It is worth mentioning that elevation model propagation can be accelerated by two techniques: parallelization and use of local data during the GP inference. Because the propagation process for each particle does not require information on other particles, it can be parallelized. Furthermore, during the GP inference, only the neighborhood data of DTED are utilized. The range of the neighborhood is defined by the spatial window size L. propagate through the dynamic model propagate through the elevation model   h

Remark on an Existing Work
As mentioned in Section 1, from a mathematical perspective, the proposed algorithm (STC-PF) is similar to scPF (soft-constrained Particle Filter) [35]. Similar to STC-PF, scPF is based on the SIR particle filter; however, the two differ in the sense that scPF utilizes generalized likelihood.ŵ where P(C k |x i k ) is a pseudo-measurement that represents how much the given state x i k satisfies the constraint. If Equation (21) is replaced by then the weight update rule is also changed.
Thus, the generalized likelihood function can be identified by equating the elevation model with the pseudo-measurement. As a result, scPF can be reduced to STC-PF as long as the assumption for target motion holds.

Scenario and Parameter Settings
To evaluate STC-PF, numerical experiments are performed with the following scenario: The radar is mounted on an aircraft that flies at a speed of 70 m/s at a height of 2500 m. The radar tracks a single target that moves along the surface at a speed of 25 m/s. (see Figure 4) The simulation runs for 100 s. Furthermore, to reflect the uncertainty in DTED, a noisy version of DTED is created. More specifically, iid zero-mean Gaussian noise with variance σ DTED is sampled and added for each data entry in DTED. Because it is reasonable to bound the uncertainty of DTED, sampled noise is clipped to 50 m if its absolute value exceeds 50 m. Values of parameters used in the simulation are listed in Table 1. Detailed explanation about the choice of GP hyper-parameters is in the Appendix B. The simulations are performed with two settings that differ in the value of σ DTED . The reasonable value for σ DTED is 3.77 m, which is inferred from [37]. However, another setting whose σ DTED is 1.89 m is also used to observe the sensitivity of the key parameter.

Name
Value

Baseline Methods
To compare STC-PF with other filters that can incorporate nonlinear constraints, the Smoothly Constrained Kalman Filter (SCKF) is implemented as well [30]. Note that 'Smoothly Constrained' in the name of SCKF does not mean soft constraint. Because SCKF can incorporate only deterministic constraints, it requires approximations of ground-truth terrain elevation that require h and ∇h to be fixed to specific values. One approach used for the comparison is to ignore the noise inherent in DTED and use bilinear interpolation to retrieve the terrain elevation at arbitrary positions.
where BL(λ, ϕ | DTED) is bilinear interpolation at (λ, ϕ) given DTED. For the computation of the gradient ∇h, central numerical differentiation is used instead of analytic differentiation to avoid non-differentiable cases.
where ∆ is a small constant. Another method is to use GP mean regression rather than bilinear interpolation. That is, h ∇h T ≈μ (28) whereμ is the GP joint mean ofh and ∇h in Equation (9). This enables us to reconstruct the most probable ground-truth terrain elevation considering the noise of DTED; however, this method still cannot consider the uncertainty of the inferredh and ∇h values, in contrast to STC-PF.
SCKF requires the Jacobian of the constraint functions: However, it is impossible to differentiate E(x k ) and E v (x k ) analytically because they involve coordinate transformation between local Cartesian and WGS84 LLA. Alternatively, the derivative can be obtained using the central numerical difference regardless of the regression method. ∂E where e x is a canonical unit vector whose first component is nonzero. ∂E/∂y k , ∂E/∂z k , and ∂E v /∂ • can be obtained in a similar way. Because E is not a function of v ·,k , corresponding derivatives automatically become zero.

Results
To evaluate STC-PF, SCKF using bilinear regression, and SCKF using GP mean regression, 100 Monte-Carlo simulations were carried out for each σ DTED value. Tracking performance is assessed based on timewise RMS (Root Mean Squared) error. For example, timewise RMS for local Cartesian x position error at time k is where N MC is the number of repetitions (i.e., 100), x n k the filter mean value for x position at time k in the n th trial, andx k the ground-truth x position at time k. The time average (10 ≤ k ≤ 90) for timewise RMS is also computed for evaluation. Figure 5 shows the timewise RMS for local Cartesian position error and velocity error. In the figures, SCKF using bilinear regression shows the worst tracking performance. In terms of time average of RMS position error, as shown in Table 2, the superiority of STC-PF over SCKF using GP mean regression is clear, although it cannot be identified in Figure 5. In terms of RMS velocity error, STC-PF distinctly outperforms the other two methods. This trend also holds for the different parameter setting, namely σ DTED = 1.89 m, as shown in Figure 6 and Table 3.    On the other hand, the speed of the algorithms is assessed based on the average processing time for a single timestep. STC-PF and SCKF both were implemented in MATLAB and run on an Intel Core i7-8565U with 16.0GB RAM. Table 4 shows that the baseline method runs nearly in real time, while STC-PF does not. Nevertheless, with parallel computing, STC-PF gets much faster and shows the possibility of real time applications.

Discussion
According to the simulation results, the tracking performance of SCKF depends on the regression method. More specifically, SCKF using GP mean regression has smaller RMS errors than SCKF using bilinear regression. This result could be due to the regression method used for target trajectory generation: GP regression, which is suitable for fitting of a smooth curve, might result in a trajectory closer to the ground-truth trajectory generated by bicubic spline interpolation.
Meanwhile, the superiority of STC-PF over SCKF in terms of tracking performance could stem from two factors. The first is that a particle filter is more expressive than a simple Gaussian filter. Because SCKF assumes that the posterior distribution is a simple Gaussian, SCKF adjusts its state to meet the constraint at one point. However, unlike SCKF, particle filters do not assume the form of a posterior distribution. Therefore, STC-PF can estimate the distribution of the velocity hypothesis for each position hypothesis, so that the combined hypothesis independently meets the terrain constraint. As a result, the filter mean value, a weighted sum of each hypothesis, is less biased.
The second reason is that the state estimation with soft constraint is less sensitive to the uncertainty of the constraint. Conversely, state estimation with hard constraint is very sensitive to the uncertainty of the constraint. Figure 7 shows that the terrain constraint (Equation (5)) holds almost perfectly assuming perfect knowledge of the position and DTED. However, the hard constraint, the GP mean regression, deviates from the groundtruth value if small amounts of position uncertainty (approximately 30 m in the longitudinal direction) and DTED noise (σ DTED = 3.77 m) are introduced (see Figure 8). In other words, uncertainty of the horizontal position and DTED can result in a catastrophic state error if a hard constraint is applied. On the other hand, STC-PF can absorb the error to some degree as most of the ground-truth value resides inside the two-sigma bound.  Regarding the speed of the algorithm, it is natural that STC-PF consumes more computational resources than SCKF. This is because STC-PF requires expensive computation for every particle; GP inference whose time complexity is O(n 3 ), where n is the number of observation points involved in the GP inference. Fortunately, time consumption is drastically reduced by the virtue of parallel computing, and it can be reduced further by two measures. One is reducing the number of GP observation points (namely, the spatial window L in Figure 3) as far as the performance degradation is negligible, and the other is increasing the number of cores dedicated to the filtering.

Conclusions
To sum up, we have proposed a particle filter to improve the performance of ground target tracking. To estimate the velocity more accurately, not only a position constraint but also a velocity constraint has been introduced in the terrain constraint. Although DTED provides terrain elevation of the entire globe, it provides inaccurate values at discrete positions. Thus, the ground-truth terrain elevation included in the terrain constraint has been modeled with a Gaussian process, and DTED has been regarded as noisy observations of it. As a result, terrain constraint has become a soft constraint that can reflect the uncertainty of DTED. Finally, we have proposed a particle filter, STC-PF, given the assumption of the motion of the target. STC-PF is based on SIR PF, but the major difference is that STC-PF uses the elevation model. Due to the elevation model, knowledge of the horizontal position and velocity of a target enables us to infer the vertical position and velocity more precisely. In the numerical simulation, STC-PF has been compared with SCKF which can incorporate hard constraints only. Furthermore, to reflect the uncertainty in DTED, filters have made use of DTED contaminated by noise, whereas the ground-truth trajectory of the target is generated by the original DTED. The simulation results showed that STC-PF outperforms SCKF in terms of RMS error, for two possible reasons. The first is that particle filters are more expressive than simple Gaussian filters. The second is that the state estimation with soft constraint is less sensitive to uncertainty of the constraint than that with hard constraint.

Informed Consent Statement: Not applicable.
Data Availability Statement: SRTM data was obtained from USGS and are available at https:// earthexplorer.usgs.gov/ with the permission of USGS. Detection data is not publicly available because this study was carried out as a defense-related project.

Conflicts of Interest:
The authors have generic conflicts with individuals in the same affiliation. The fourth and the fifth authors have general conflicts with individuals in the defense contractors in Korea.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Derivative of Gaussian Process
A Gaussian process is a collection of random variables, any finite number of which have a joint Gaussian distribution. A Gaussian process is completely specified by its mean function and covariance function. The mean function and covariance function are defined as respectively [40]. Following the most common choice, zero mean function and squared exponential kernel are used throughout this paper. Namely, where x Γ = √ x T Γx. In most practical applications, we do not have direct access to function values themselves, only noisy versions thereof: where ∼ N(0, σ). Let X and Y be the concatenation of all observation points and corresponding measurements, respectively. Given the observation set (X, Y), the predictive distribution of the function value f * at arbitrary test points X * can be derived. Starting from the joint distribution of the observation set (X, Y) and the test set (X * , f * ), k(X, X) + σ 2 I k(X, X * ) k(X * , X) k(X * , X * ) .
Furthermore, consider the derivative of the given GP. As differentiation is a linear operator, the derivative of a Gaussian process remains a Gaussian process as long as the kernel function is differentiable [41]. To find the joint probability of the observation set (X, Y) and the derivative observation set (X d , Y d ), covariance between the function value and the derivative value as well as covariance among the derivative values should be described. First, let the derivative of the underlying function be where D is the dimension of x. Then, the explicit expressions of the new covariance functions are Using basic arithmetic operations, the above expression is reduced to vector form: