V-Spline: An Adaptive Smoothing Spline for Trajectory Reconstruction

Trajectory reconstruction is the process of inferring the path of a moving object between successive observations. In this paper, we propose a smoothing spline—which we name the V-spline—that incorporates position and velocity information and a penalty term that controls acceleration. We introduce an adaptive V-spline designed to control the impact of irregularly sampled observations and noisy velocity measurements. A cross-validation scheme for estimating the V-spline parameters is proposed, and, in simulation studies, the V-spline shows superior performance to existing methods. Finally, an application of the V-spline to vehicle trajectory reconstruction in two dimensions is given, in which the penalty term is allowed to further depend on known operational characteristics of the vehicle.


Introduction
Global Positioning System (GPS) technology has become an essential tool in a wide range of applications involving moving vehicles, from transport management [1] and traffic safety studies [2] to modern precision farming [3]. Nevertheless, the accuracy of GPS tracking seems to be neglected in many applications [4,5]. Even if an accurate GPS device is utilized, GPS remains subject to various systematic errors due to the number of satellites in view, uncertainty in satellite orbits, clock and receiver issues, etc. [6,7]. These measurements are usually irregularly recorded, leading to what is known as irregularly spaced or intermittent data. Reconstruction or forecasting based on irregularly spaced data is usually more complicated and less accurate than that based on regularly spaced data [8].
To be fit for purpose, trajectory reconstruction must be accurate and robust. Two key issues for reconstruction are (i) how to handle observations that are inherently noisy measurements of the truth, and (ii) how to interpolate appropriately between observations, also known as path interpolation. In this context, statistical smoothing techniques can be useful processing tools because they are designed to minimize the impact of random error, and still typically require less time to detect random errors than visual inspection [9].
It has been shown previously that, if kinematic information such as velocity and acceleration can be included, interpolation and hence trajectory reconstruction can be greatly improved. The authors in [10,11] used B-splines to give a closed-form expression for a trajectory with continuous second derivatives that passes through the position points smoothly while ignoring outliers. The authors in [12] presented a quintic spline trajectory reconstruction algorithm connecting a series of reference knots that produces continuous position, velocity, and acceleration profiles in the context of computer (or computerized) numerical control (CNC). The authors in [13] gave a piecewise cubic reconstruction found by matching the observed position and velocity at the endpoints of each interval; this is essentially a Hermite spline. The authors in [14] also used Hermite interpolation to fit position, velocity and acceleration with given kinematic constraints. The authors in [15] implemented spline-based trajectories in order to overcome parametric singularities that occur in some reconstruction methods. The author in [16] proposed the kinematic interpolation approach that uses a set of kinematic equations to describe the motion of an object in terms of polynomial splines. Based on an adaptive cubic spline interpolation, the authors in [17] proposed an approach that, in the context of the Aircraft Communication Addressing and Reporting System, improves the smoothness and precision of trajectory reconstruction.
These approaches focus on optimal paths that are typically the shortest in either distance or time between starting and end points. Additionally, in some approaches, the moving object is assumed to be a point-like object. In this case, the object can rotate about itself to orient along the path in the direction of the goal point [15]. This assumption is unlikely to be appropriate for a real vehicle or vessel, particularly a tractor, which is the motivating example in our study.
Modern farming relies on the precise application of fertilizers, pesticides and irrigation. Large commercial farms typically operate a fleet of farm vehicles for these tasks, and it is of crucial importance for economic, environmental and regulatory reasons that the location and operational characteristics of these vehicles are recorded systematically and accurately. In order to do this, it is becoming standard to equip farm vehicles with GPS units to record the location of the vehicle on the farm. It is the goal of this study to develop an appropriate tool to reconstruct vehicle trajectories from such data, particularly when it is intermittent and noisy.
In this study, we assume that we have independent records of the position and velocity of a moving object at a sequence of observation times. Traditional methods often assume motion with constant speed between two observations times, but this will not work well in our case. Additionally, motivated by the fact that tractors often work in open fields, we assume no further information is available to constrain the position of the object. Initially, we constructed trajectories in terms of a Hermite cubic spline basis [18,19]. In each interval, the reconstruction is clearly continuous, as are its first and second derivatives. The goal is then to connect the piecewise splines keeping the trajectory and its first derivative continuous at the interior knots. In this approach, the trajectory is not required to pass through each knot and the main objective is the smoothness of the path, not a shortest or minimum-time path. To formalize this procedure, we propose a new objective function that incorporates velocity information and includes an adaptive penalty term. The penalty term utilises information about the distance and travel time on each interval. We dub the proposed smoothing spline the V-spline because it incorporates velocity information and can be applied to vehicle and vessel tracking. We show that the V-spline works better than other methods in simulation studies and that it produces satisfactory outcomes in a real-world application.
The structure of this paper is as follows: in Section 2, we introduce the basis functions and the V-spline objective function that depends both on position residuals y i − f (t i ) and velocity residuals v i − f (t i ). A new parameter γ in the objective function controls the degree to which the velocity information is used in the reconstruction. We show that the V-spline can be written in terms of modified Hermite spline basis functions. We also introduce a particular adaptive V-spline that seeks to control the impact of irregularly sampled observations and noisy velocity measurements. In Section 3, a cross-validation scheme for estimating the V-spline parameters is given. Section 4 details the performance of the V-spline on simulated data based on the Blocks, Bumps, HeaviSine and Doppler test signals [20]. Finally, an application of the V-spline to a two-dimensional data set is presented in Section 5. R code for implementing V-spline and reproducing our outcomes is provided as Appendices A-C at the end of the manuscript.

Objective Function
Conventional smoothing spline estimates of f (t) appear as a solution to the following minimization problem: findf ∈ C (2) [a, b] that minimizes the penalized residual sum of squares, for a pre-specified value λ > 0 [21][22][23]. The objective function combines goodness-of-fit to the data with a measure of roughness [24]. For V-splines, we consider the situation of paired position data y = {y 1 , . . . , y n } and velocity data v = {v 1 , . . . , v n } at a sequence of times satisfying where the second derivative of f (t) is piecewise continuous, we define the objective function where γ > 0, and we have chosen the penalty function λ(t) to be a piecewise constant function on interior intervals, i.e for t ∈ [t i , t i+1 ), i = 1, . . . , n − 1, In fact, each f i ∈ C (2) [t i ,t i+1 ] is a Hermite spline which satisfies the properties of a cubic spline. The complete spline function f , which connects all f i s, has piecewise continuous second derivative, and will be continuous if a particular condition is met. The second derivative f is zero on the exterior intervals [a, t 1 ] and [t n , b]. From now on, we will understand λ(t) to be piecewise constant (3), and we will often use λ to refer to the set of λ i . The proof of Theorem 1 is in Appendix B. Remark: In the language of splines, the points t 1 , . . . , t n are the interior knots of the V-spline, and a = t 0 , b = t n+1 are the exterior or boundary knots.

Basis Functions
The cubic Hermite spline where the basis functions are h (i) For V-splines, a slightly more convenient basis is given by Therefore, any f ∈ C (2) p.w. [a, b] can then be represented in the form where {θ k } 2n k=1 are parameters corresponding with the "true" position f (t i ) and velocity f (t i ) at the observation points.

Computing the V-Spline
In terms of the basis functions in the previous section, the objective function (2) is given by where B and C are n × 2n matrices with components and Ω λ is a 2n × 2n matrix with components [ In the following, we reserve the use of boldface for n × 1 vectors and n × n matrices.
The detailed structure of Ω λ is presented in Appendix A. It is convenient to write It is then evident that Ω λ is a bandwidth four matrix.
Since Equation (10) is a quadratic form in terms of θ, it is straightforward to establish that the objective function is minimized at which can be identified as a generalized ridge regression. The fitted V-spline is then given byf (t) = ∑ 2n k=1 N k (t)θ k . The V-spline is an example of a linear smoother [25]. This is because the estimated parameters in Equation (13) are a linear combination of y and v. Denoting byf andf the vector of fitted valuesf (t i ) andf (t i ) at the training points t i , we havê where S λ,γ , T λ,γ , U λ,γ and V λ,γ are smoother matrices that depend only on t i , λ(t) and γ.
It is not hard to show that S λ,γ and V λ,γ are symmetric, positive semi-definite matrices. Note that T λ,γ = U λ,γ .

Adaptive V-Spline
Until now, we have not explicitly considered the impact of irregularly sampled observations of noisy measurements of velocity on trajectory reconstruction. In order to do this, it is instructive to evaluate the contribution to the penalty term from the interval [t i , t i+1 ). Using (4), it is relatively straightforward to show that where is the average velocity over the interval. The ε ± i can be interpreted as the difference at time t i and t − i+1 respectively between the velocity implied by an interpolating Hermite spline and the velocity implied by a straight line reconstruction.
The contribution to the penalty term is then where As a consequence of (17), larger time intervals will tend to contribute less to the penalty term (other things being equal). However, this is exactly when we would expect the velocity at the endpoints of the interval to provide less useful information about the trajectory over the interval. In the case when the observed change in position is small, i.e., when y i+1 − y i =v i ∆T i ≈ 0, over-reliance on noisy measurements of velocity will result in "wiggly" reconstructions. In these two instances-graphically depicted in Figure 1a-we would like the V-spline to adapt and to favor straighter reconstructions; this is a deliberate design choice. We can achieve this by choosing where η is a parameter to be estimated. The penalty term then takes a particularly compelling form: the contribution from the interval [t i , t i+1 ) (17) is proportional to discrepancy in velocity average velocity 2 (19) for all i. We call the resulting spline the adaptive V-spline. The spline when λ i = λ 0 or, more accurately, when λ i is independent of ∆T i andv i , we call it the non-adaptive V-spline. Figure 1. Comparing cubic Hermite spline reconstruction and straight line reconstruction. When ∆T i = t i+1 − t i is large orv i ∆T i = y i+1 − y i is small, the adaptive V-spline favours straighter reconstructions.

Parameter Selection and Cross-Validation
The issue of choosing the smoothing parameter is ubiquitous in curve estimation and there are two different philosophical approaches to the problem. The first is to regard the free choice of smoothing parameter as an advantageous feature of the procedure. The second is to let the data determine the parameter [22,26], using a procedure such as crossvalidation (CV) or generalized cross-validation (GCV) [21]. We prefer the latter and use the data with GCV to train our model and find the best parameters.
In standard regression, which assumes the mean of the observation errors is zero, the true regression curve f (t) has the property that, if an observation y i is omitted at time point t i , the value f (t i ) is the best predictor of y i in terms of mean squared error [22]. We use this observation to motivate a leave-one-out cross-validation scheme to estimate λ and γ for both the non-adaptive and the adaptive V-splines.
The following theorem establishes that we can compute the cross-validation score without knowing thef (−i) (t, λ, γ): Theorem 2. The cross-validation score of a V-spline satisfies arg min λ,γ>0 wheref is the V-spline smoother calculated from the full data set with smoothing parameter λ and γ, and S ii = [S λ,γ ] ii , etc.
The proof of Theorem 2 is in Appendix C.

Simulation Study
In this section, we give an extensive comparison of methods for equal-spaced data. The comparison is based on the ability to reconstruct trajectories derived from Blocks, Bumps, HeaviSine and Doppler, which were used in [20,27,28] to mimic problematic features in imaging, spectroscopy and other types of signal processing.
Letting g(t) denote any one of Blocks, Bumps, HeaviSine or Doppler, we treat g(t) as the instantaneous velocity of the trajectory f (t) at time t, i.e., f (t) = g(t). Setting f (t 1 ) = 0, the position is then updated in terms of the average velocity over each interval: which is accurate to the second order in t i+1 − t i . Finally, the observed position and velocity are found by adding i.i.d. zero-mean Gaussian noise: where ε , σ f is the standard deviation of the positions f (t i ), σ g is the standard deviation of the velocities g(t i ), and SNR is the signal-tonoise ratio, which we take to be 3 or 7.
We compare the performance of the adaptive V-spline with a spatially adaptive penalized spline known as the P-spline with the function asp2 from the package AdaptFi-tOS [29][30][31], a generalized additive model gam from the package mgcv [32,33], the kinematic interpolation approach (KI) by [16], as well as the adaptive V-spline with γ = 0, which becomes a conventional spline with Hermite basis functions, and the non-adaptive V-spline where λ 0 is a constant. It is important to note that only the KI approach, the non-adaptive and adaptive V-splines incorporate velocity information. The V-spline parameters are obtained by minimizing the cross-validation score (22). In the gam model, we use tp basis functions with 1024 knots. For the KI approach, the position at time t i is interpolated from the two neighbouring points at t i−1 and t i+1 . (The positions at t 1 and t n are interpolated from points at (t 1 , t 2 ) and (t n−1 , t n ), respectively.) Following [34], we fix n = 1024 in the simulations.
To examine the performance of the adaptive V-spline, we compute the true mean squared error for each of the reconstructions via: and the Modified Nash-Sutcliffe efficiency (mNSE) [35] via: The results are shown in Tables 1 and 2. The V-spline, either adaptive or non-adaptive, returns the best solution in all cases. The reason for the poor performance of kinematic interpolation is two-fold: first, KI assumes v i is a good approximation to the velocity over the entire interval [t i−1 , t i+1 ). Second, KI is not a true smoother so it is prone to errors in the observations. In contrast, the V-spline successfully smooths and interpolates in the presence of noise. Table 3 shows the ability of the adaptive V-spline to retrieve the true SNR: for reconstructionf , it is estimated by σf /σ (f −y) . Table 3 shows that the estimates from the V-splinê f are very close to the true values. In summary, the simulation study has shown the ability of V-splines to accurately reconstruct trajectories from noisy and potentially problematic velocity profiles. The V-spline outperforms methods that do not use velocity information, and its smoothing strategy appears to be vastly superior to that of kinematic interpolation.

Inference of Tractor Trajectory
In this section, we apply the V-spline to a data set obtained from a GPS unit mounted on a tractor working in a horticultural setting. The motivating problem in this context is to accurately record where pesticide has been applied to ensure that neither over-spraying or under-spraying has occurred.
GPS units in vehicles provide y t , noisy measurements of the actual position x t , and v t , noisy measurements of the actual velocity u t , for a sequence of times t ∈ T, which is irregularly recorded with highly variable time differences ∆T i . These data may also be augmented with information on operating characteristics of the vehicle, b t , in this case data on whether the tractor boom was in a raised or lowered position. The trajectory reconstruction problem is the problem of estimating x s , for an arbitrary time s, given a subset of the observations {y t , v t , b t | t ∈ T}. Note that, in this definition of trajectory reconstruction, we are not explicitly interested in estimating u s .
The original data set consists of n = 928 records of longitude, latitude, speed, bearing and the status of the tractor's boom sprayer. The boom status, "up" and "down", denotes the operational state of the tractor, and indicates different types of trajectories. For example, if boom status is "down", the tractor is probably sowing, watering or harvesting on the farm. In this scenario, the speed is stable and its variance is low. On the contrary, when it is "up", the speed could be high because the driver is travelling between jobs, it could be zero because the driver is having a break, or it might indicate the tractor is turning. In this last situation, however, the acceleration could be high. For this reason, we add further complexity to the model by allowing the penalty parameter to depend on boom status.
For trajectory reconstruction, this data set was converted from longitude and latitude in degrees ( • ) into easting and northing in meters (m) by the Universal Transverse Mercator (UTM) coordinate system. The speed and bearing were converted into velocities (m/s) in those directions as well. See Figure 2.

The V-Spline in d-Dimensions
To generalize the V-spline to d-dimensions, we consider the situation preceding Equation (2) but where now y i , v i ∈ R d . Then, the function f : [a, b] → R d is a ddimensional V-spline if it minimizes: where · 2 is the Euclidean norm in d-dimensions. For each direction α = 1, . . . , d, the fitted V-spline has the formf α (t) = ∑ 2n k=1 N k (t)θ α k , wherê The parameters λ and γ are estimated by minimizing the cross-validation score: arg min λ,γ>0 In what follows, we allow the non-adaptive and adaptive V-splines to depend on the boom status. This is to demonstrate that our method can simply and usefully also incorporate known covariates. In this application, letting b i = 0 denote boom "up", b i = 1 denote boom "down", andv i = y i+1 − y i 2 /∆T i be the average velocity on the interval [t i , t i+1 ), the penalty term for the non-adaptive V-spline is and, for the adaptive V-spline, it is Optimization in (29) is now simply with respect to positive λ d , λ u and γ.

Two-Dimensional Trajectory Reconstruction
The V-spline reconstruction from the tractor data is shown in Figure 3. The parameters λ d , λ u and γ are found by our proposed cross-validation scheme using the stats ::optim function in R [36]. It is immediately evident from the trajectory that the tractor has been moving up and down rows of an orchard or travelling between parts of the orchard. It is instructive to compare the performance of the adaptive V-spline to a line-based approach that simply and unrealistically connects observations by a straight line, kinematic interpolation which also utilizes velocity information, and the non-adaptive V-spline. Figure 4 shows finer detail of the tractor trajectory given by these reconstructions. A feature of the KI method is the hugely unrealistic excursions near the turn-around points at the end of each row as shown in Figure 4b. On the contrary, the adaptive V-spline, see Figure 4d, adapts to the information based on observed velocity discrepancy to avoid such excursions. Without the adaptive term (18), the non-adaptive V-spline performs in a similar way to KI, which can be seen from Figure 4c; this proves the power of the adaptive penalty.  KI relies on misleading velocity information that generates unrealistic trajectories. The non-adaptive V-spline also generates unrealistic trajectories at sharp-turning and braking points. However, the adaptive V-spline adapts to the information based on observed velocity discrepancy and generates plausible trajectories.

Discussion
In this paper, a smoothing spline called the V-spline is proposed that minimizes an objective function which incorporates both position and velocity information. Given n knots, the V-spline has 2n effective degrees of freedom corresponding to n − 1 cubic polynomials with their value and first derivative matched at the n − 2 interior knots. The effective degrees of freedom are then fixed by n position observations and n velocity observations. Note that, in the limit γ → 0, the V-spline reduces to having n effective degrees of freedom. An adaptive version of the V-spline is also introduced that seeks to control the impact of irregularly sampled observations and noisy velocity measurements.
The computational complexity of the V-spline method is equivalent to any smoothing spline that uses a cross-validation procedure to estimate the tuning parameters. The essential difference is that the V-spline incorporates 2n data points (in each dimension), as opposed to n. The impact of this shows up in the time to solve forθ in (13). Thus, the computation time of the V-spline is the same as a standard smoothing spline with 2n observations. Modest computational gains can possibly be made by improving the CV parameter estimation step, but Theorem 2 already assures us that this step is highly efficient. Future research directions for the V-spline include application to ship tracking [18] and development of a fast filtering algorithm.

Conflicts of Interest:
The authors declare no conflict of interest.